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

[Paolo Giannozzi]

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?

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.

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