[Pw_forum] charge density with denser grid and wrong cube file

[Thomas Gruber]

> Not sure what the problem with the "cube" file is.
> FFT arrays obey periodic boundary conditions: the
> index runs from 1 to nr, and f(nr+1) =>  f(1) . Is
> this the problem?
Yes this is the problem. In the cube file the last data point should be 
the same like the first one (cell boundry), as it is implemented in for 
the xsf file. Just adding the first data point along each direction an 
the end to get a grid that is +1 in each direction.
>
> About the denser grid: I don't think it is that simple.
> Either you perform a Fourier interpolation (have a look
> at code "voronoy.f90" in v.4.1.3) or you perform a local
> interpolation in real space.
I can note find the  "voronoy.f90". When using nx,ny,nz for a plane oder 
orthorhombic box, is it also an interpolation in real space? 
("Otherwise, the required 3D grid is generated from the Fourier 
components (may be VERY slow)" from the pp.x documentation) Can these 
parameters be used to calculate a data grid that is denser than the FFT 
grid? It is ok if it is a multiple of the FFT grid.
> P.
>
> On Mon, 2013-02-18 at 16:56 +0100, Thomas Gruber wrote:
>> >  Hallo all,
>> >  
>> >  I have two problems:
>> >  1. I like to have the charge_density.xsf/*.cube with output_format=5/6,
>> >  but with a denser grid than the FFT grid. nx,ny,nz are ignored in this
>> >  case. They work for output_format=3 with an orthorhombic box, but I want
>> >  my unit cell.
>> >  
>> >  2. The output format for *.xsf works fine, but not for *.cube. For xsf
>> >  the data_grid is x+1,y+1,z+1, but not for cube (x,y,z). The *.cube also
>> >  needs x+1,y+1,z+1 to get the write volume, atom positions and no offset
>> >  along the unit cell surface for the charge density. I tried to correct
>> >  the cube.f90 in PP but the charge density data is still not correct. I
>> >  increased nr# by 1 in each case:
>> >  
>> >      WRITE(ounit,*) 'Cubfile created from PWScf calculation'
>> >      WRITE(ounit,*) ' Total SCF Density'
>> >  !                        origin is forced to (0.0,0.0,0.0)
>> >      WRITE(ounit,'(I5,3F12.6)') nat, 0.0d0, 0.0d0, 0.0d0
>> >      WRITE(ounit,'(I5,3F12.6)') nr1+1, (alat*at(i,1)/dble(nr1),i=1,3)
>> >      WRITE(ounit,'(I5,3F12.6)') nr2+1, (alat*at(i,2)/dble(nr2),i=1,3)
>> >      WRITE(ounit,'(I5,3F12.6)') nr3+1, (alat*at(i,3)/dble(nr3),i=1,3)
>> >  
>> >      DO i=1,nat
>> >         nt = ityp(i)
>> >         ! find atomic number for this atom.
>> >         at_num = atomic_number(trim(atm(nt)))
>> >         at_chrg= dble(at_num)
>> >         ! at_chrg could be alternatively set to valence charge
>> >         ! positions are in cartesian coordinates and a.u.
>> >         !
>> >         ! wrap coordinates back into cell.
>> >         tpos = matmul( transpose(bg), tau(:,i) )
>> >         tpos = tpos - nint(tpos - 0.5d0)
>> >         inpos = alat * matmul( at, tpos )
>> >         WRITE(ounit,'(I5,5F12.6)') at_num, at_chrg, inpos
>> >      ENDDO
>> >  
>> >      DO i1=1,nr1+1
>> >         DO i2=1,nr2+1
>> >            WRITE(ounit,'(6E13.5)') (rho(i1,i2,i3),i3=1,nr3+1)
>> >         ENDDO
>> >      ENDDO
>> >  
>> >  Can someone give me the correction for cube.f90 and a way to increase
>> >  the data_grid density for the output_format=5/6?
>> >  
>> >  Sincerely,
>> >  
>> >  Thomas Gruber
>> >  
> -- Paolo Giannozzi, IOM-Democritos and University of Udine, Italy

Sincerely,

Thomas Gruber

-- 
Thomas Gruber
Institut für Theoretische Physik
TU Bergakademie Freiberg
Leipziger Str. 23, 09599 Freiberg, Germany
Tel: +49 (0)3731 392006
Email: thomas.gruber at physik.tu-freiberg.de