modulegrid.F90

Go to the documentation of this file.

! vim:fdm=marker:
!------------------------------------------------------------------------
MODULE modulegrid
IMPLICIT NONE

CONTAINS

!------------------------------------------------------------------------
SUBROUTINE get_positionGrid(pos_i,pos_f,gdim,grid)
!---{{{ 
IMPLICIT NONE

INTEGER,INTENT(IN)    :: gdim   !< number of grid points in this direction
REAL*8,INTENT(IN)     :: pos_i  !< left border of box
REAL*8,INTENT(IN)     :: pos_f  !< right border of box( =last grid point + dx)
REAL*8,INTENT(OUT)    :: grid(gdim) !< the position grid   

REAL*8                :: dx
INTEGER               :: i

dx             =     (pos_f-pos_i)/(1.d0*gdim)
DO i=1,gdim
  grid(i) = pos_i+(i-1)*dx
END DO
!---}}}
END SUBROUTINE get_positionGrid

!------------------------------------------------------------------------
SUBROUTINE get_momentumGrid(xl,xu,gdim,momentumGrid)
USE moduleparameters, ONLY: PI
!---{{{ returns momentum grid in ascending order
IMPLICIT NONE

INTEGER,INTENT(IN)              :: gdim !< number of grid points in this direction
REAL*8,INTENT(IN)     :: xl   !< left border of box
REAL*8,INTENT(IN)     :: xu   !< right border of box( =last grid point + dx)
REAL*8,INTENT(OUT)    :: momentumGrid(gdim) !< array of momentum space grid points in ascending order

REAL*8                :: dp
INTEGER                         :: i

IF (gdim>1) then 
   dp             =     2.d0*PI/(xu-xl)

   IF (MOD(gdim,2)==0) then !gdim even

      DO i=1,gdim
         momentumGrid(i)=(i-gdim/2)*dp  ! e.g.: (-dp,0,dp,2dp)
      END DO

   ELSE ! gdim uneven

      DO i=1,gdim
         momentumGrid(i)=(i-gdim/2-1)*dp ! e.g.: (-2dp,-dp,0,dp,2dp)
      END DO

   END IF 

ELSE IF (gdim==1) then

   momentumGrid(1) = 0.d0

ELSE

   write(*,*)"ERROR in dimension"
   STOP

END IF

!---}}}
END SUBROUTINE get_momentumGrid


!------------------------------------------------------------------------
SUBROUTINE get_FFTmomentumGrid(pos_i,pos_f,gdim,FFTmomentumGrid)
USE moduleparameters,ONLY:PI
!---{{{ returns momentum grid in fft order
IMPLICIT NONE

INTEGER,INTENT(IN)    :: gdim  !< the number of grid points in this direction
REAL*8,INTENT(IN)     :: pos_i !< left border of the box
REAL*8,INTENT(IN)     :: pos_f !< right border of the box ( =last grid point + dx)
REAL*8,INTENT(OUT)    :: FFTmomentumGrid(gdim) !< array of momentum space grid points in FFT order

REAL*8                :: dp
INTEGER               :: i

dp             =     2.d0*pi/(pos_f-pos_i)
DO i=1,gdim
  FFTmomentumGrid(i) = (i-1)*dp
END DO

DO i=gdim/2+2,gdim !put momenta in FFT order (0,dp,...,N/2*dp,(-N/2+1)*dp,...,-dp)
  FFTmomentumGrid(i)=1.d0*(i-gdim-1)*dp 
END DO

!---}}}
END SUBROUTINE  get_FFTmomentumGrid

!------------------------------------------------------------------------
SUBROUTINE get_p2over2m(NDX,NDY,NDZ,xi,xf,yi,yf,zi,zf,Dims,p2over2m) 
! {{{
IMPLICIT NONE

REAL*8,INTENT(IN)     :: xi  !< left border of box
REAL*8,INTENT(IN)     :: xf  !< right border of box( =last grid point + dx)
REAL*8,INTENT(IN)     :: yi  !< left border of box
REAL*8,INTENT(IN)     :: yf  !< right border of box( =last grid point + dy)
REAL*8,INTENT(IN)     :: zi  !< left border of box
REAL*8,INTENT(IN)     :: zf  !< right border of box( =last grid point + dz)
INTEGER,INTENT(IN)    :: NDX !< number of grid points in x
INTEGER,INTENT(IN)    :: NDY !< number of grid points in y
INTEGER,INTENT(IN)    :: NDZ !< number of grid points in z
INTEGER,INTENT(IN)    :: Dims !< number of dimensions
REAL*8,INTENT(OUT)    :: p2over2m(NDX*NDY*NDZ)  !< the array p^2/2m in FFT order

INTEGER               :: i,j,k,m,n
REAL*8                :: FFTmomentumGridx(NDX),FFTmomentumGridy(NDY),FFTmomentumGridz(NDZ),dpx,dpy,dpz


FFTmomentumGridx = 0.d0
FFTmomentumGridy = 0.d0
FFTmomentumGridz = 0.d0

IF(NDX>1) THEN
!---{{{ get momentumGridx in fft order
  CALL  get_FFTmomentumGrid(xi,xf,NDX,FFTmomentumGridx)
!  write(6,*)"FFTmomentumGridx",FFTmomentumGridx
!  write(6,*)"xi=",xi
!  write(6,*)"xf=",xf
!---}}}
END IF

IF (NDY>1) THEN
!---{{{ get momentumGridy in fft order
  CALL  get_FFTmomentumGrid(yi,yf,NDY,FFTmomentumGridy)
!---}}}
END IF

IF (NDZ>1) THEN
!---{{{ get momentumGridz in fft order
  CALL  get_FFTmomentumGrid(zi,zf,NDZ,FFTmomentumGridz)
!---}}}
END IF

DO m=1,NDX*NDY*NDZ

  CALL get_ijk_from_m(m,NDX,NDY,i,j,k)

! multiply by px**2+py**2+pz**2

  IF(NDX>=1) THEN

    p2over2m(m)=FFTmomentumGridx(i)**2 
!    write(6,*)"i",i,"FFTmomentumGridx(i)",FFTmomentumGridx(i)
  END IF

  IF (NDY>1) THEN
    p2over2m(m) = p2over2m(m) + FFTmomentumGridy(j)**2 
  END IF

  IF (NDZ>1) THEN
    p2over2m(m) = p2over2m(m) + FFTmomentumGridz(k)**2 
  END IF

  p2over2m(m) = 0.5d0 * p2over2m(m)

END DO

!---}}}
END SUBROUTINE get_p2over2m

!------------------------------------------------------------------------
SUBROUTINE get_ijk_from_m(m,NDX,NDY,i,j,k)
!---{{{
IMPLICIT NONE
INTEGER,INTENT(IN)      :: m    !< the index that runs over all grid points in 1d,2d and 3d
INTEGER,INTENT(IN)      :: NDX  !< number of grid points in x-direction
INTEGER,INTENT(IN)      :: NDY  !< number of grid points in y-direction
INTEGER,INTENT(OUT)     :: i    !< index of grid point x_i = positionGridx(i)
INTEGER,INTENT(OUT)     :: j    !< index of grid point y_j = positionGridx(j)
INTEGER,INTENT(OUT)     :: k    !< index of grid point z_k = positionGridx(k)


INTEGER                 :: n

  n = m
  k = (n-1)/(NDX*NDY) + 1 
  n =  n-(k-1)*NDX*NDY
  j = (n-1)/NDX + 1
  n =  n-(j-1)*NDX
  i =  n

!---}}}
END SUBROUTINE get_ijk_from_m

!-----------------------------------------------------------------------
SUBROUTINE check_grid(Dims,NDX,NDY,NDZ)
!---{{{
IMPLICIT NONE
INTEGER,INTENT(IN)         :: Dims !< number of dimensions used in the current problem
INTEGER,INTENT(IN)         :: NDX  !< number of dimensions used in direction x
INTEGER,INTENT(IN)         :: NDY  !< number of dimensions used in direction y
INTEGER,INTENT(IN)         :: NDZ  !< number of dimensions used in direction z

IF((NDY>NDX).or. (NDZ>NDY)) THEN
  write(6,*)"bad grid dimensions"
  write(6,*)"correct form : NDX >= NDY >= NDZ"
  write(6,*)"NDX = ",NDX
  write(6,*)"NDY = ",NDY
  write(6,*)"NDZ = ",NDZ
  STOP
END IF

IF ((Dims==1).and.((NDY>1).or.(NDZ>1))) THEN 
  write(6,*)"inconsistency in grid and dimensions"
  write(6,*)"Dims = ",Dims
  write(6,*)"NDX = ",NDX
  write(6,*)"NDY = ",NDY
  write(6,*)"NDZ = ",NDZ
  STOP
ELSE IF ((Dims==2).and.((NDZ>1).or.(NDY<2))) THEN 
  write(6,*)"inconsistency in grid and dimensions"
  write(6,*)"Dims = ",Dims
  write(6,*)"NDX = ",NDX
  write(6,*)"NDY = ",NDY
  write(6,*)"NDZ = ",NDZ
  STOP
ELSE IF ((Dims==3).and.(NDZ<2)) THEN 
  write(6,*)"inconsistency in grid and dimensions"
  write(6,*)"Dims = ",Dims
  write(6,*)"NDX = ",NDX
  write(6,*)"NDY = ",NDY
  write(6,*)"NDZ = ",NDZ
  STOP
END IF

IF (NDX<8) THEN
  write(6,*)"grid too small. NDX must be >= 8"
  STOP
END IF
!---}}}
END SUBROUTINE check_grid

!-----------------------------------------------------------------------
SUBROUTINE  get_defaultOrbitals(NDX,NDY,NDZ,Dims,gridx,gridy,gridz,MOrb,Psi)
!---{{{
IMPLICIT NONE                           
INTEGER,INTENT(IN)              :: NDX  !< number of dimensions used in direction x
INTEGER,INTENT(IN)              :: NDY  !< number of dimensions used in direction y
INTEGER,INTENT(IN)              :: NDZ  !< number of dimensions used in direction z
INTEGER,INTENT(IN)              :: Dims !< number of dimensions used in the current problem
INTEGER,INTENT(IN)              :: MOrb !< number of orbitals used in the current problem
REAL*8,INTENT(IN)               :: gridx(NDX) !< x-position grid
REAL*8,INTENT(IN)               :: gridy(NDY) !< y-position grid
REAL*8,INTENT(IN)               :: gridz(NDZ) !< z-position grid
COMPLEX*16,INTENT(OUT)          :: Psi(NDX*NDY*NDZ,MOrb) !< array of default orbitals (gaussians)

INTEGER                         :: m,n,i,j,k
REAL*8                          :: sigmax,sigmay,sigmaz,x,y,z

sigmax = 1.d0
sigmay = 1.d0
sigmaz = 1.d0

IF (Dims==1) THEN

  DO n=1,MOrb
    DO  m=1,NDX

      x=gridx(m)
      Psi(m,n) =  x**(n-1)* EXP(-x**2/(2.d0*sigmax)**2) 
   
    END DO 
  END DO 

ELSE IF (Dims==2) THEN

  DO n=1,MOrb
    DO  m=1,NDX*NDY
      CALL get_ijk_from_m(m,NDX,NDY,i,j,k)
      x=gridx(i)
      y=gridy(j)
   
      Psi(m,n) =  x**(n-1)* EXP(-x**2/(2.d0*sigmax**2)) * EXP(-y**2/(2.d0*sigmay**2)) 
   
    END DO 
  END DO 


ELSE IF (Dims==3) THEN

  
  DO n=1,MOrb
    DO  m=1,NDX*NDY*NDZ

      CALL get_ijk_from_m(m,NDX,NDY,i,j,k)
      x=gridx(i)
      y=gridy(j)
      z=gridz(k)
  

      Psi(m,n) =  x**(n-1)* EXP(-x**2/(2.d0*sigmax**2)) * EXP(-y**2/(2.d0*sigmay**2)) *  &
                  EXP(-z**2/(2.d0*sigmaz**2)) 


    END DO 
  END DO 

ELSE 

  WRITE(*,*) "ILLEGAL VALUE FOR DIMS IN get_defaulOrbitals:",Dims
  STOP

END IF


!---}}}
END SUBROUTINE  get_defaultOrbitals





END MODULE modulegrid