modulecomputationcore.F90

Go to the documentation of this file.

! vim:fdm=marker

!---------------------------------------------------------------------------------------
MODULE modulecomputationcore
IMPLICIT NONE

CONTAINS

!----------------------------------------------------------------------------------------
SUBROUTINE evaluate_error(Delta_CVec,Delta_Psi,rho_jk_zero,CVecError,PsiError,nConf)
USE moduleinputvariables
!---{{{
IMPLICIT NONE
INTEGER*4,INTENT(IN)            :: nConf
COMPLEX*16,INTENT(IN)           :: Delta_CVec(nConf),Delta_Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(IN)           :: rho_jk_zero(MOrb,MOrb)
REAL*8,INTENT(INOUT)            :: CVecError,PsiError
COMPLEX*16                      :: complexPsiError,O_jk

INTEGER                         :: j,k

CVecError               =  0.d0
PsiError                =  0.d0
complexPsiError         = (0.d0,0.d0)

!compute the error coming from the coefficient vector
CALL normvxz(Delta_CVec,CVecError,nConf)
CVecError = 1.d0 / 4.d0 * CVecError**2

!compute the error coming from the orbitals
DO j=1,MOrb
   
   DO k=1,MOrb
      
      CALL vvaxzz (Delta_Psi(:,j),Delta_Psi(:,k),O_jk,NDX*NDY*NDZ) ! O_jk = <Delta phi_j|Delta phi_k>
      complexPsiError =  complexPsiError + O_jk * rho_jk_zero(j,k)

   END DO

END DO

PsiError = DBLE(complexPsiError)


IF (ABS(IMAG(complexPsiError))>1.0D-14) THEN
   
   write(*,*)"complexPsiError",complexPsiError

END IF
!---}}}
END SUBROUTINE evaluate_error


!----------------------------------------------------------------------------------
SUBROUTINE FuncOrbs(AbsTime,Psi_WORK,DtPsi_WORK,CData,RData,IData,LData)
USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb,Dims,lambda0
!---{{{  evaluates DtPsi of the orbitals, i.e. the right hand side of the
!        orbital EOM of MCTDHB
!
!  |phi_j^dot> = prefactor P [ h |phi_j>}  + Sum_ksql (rho^-1)_jk rho_ksql W_sl |phi_k> ]
!  for j=1,...,MOrb
!  and where P = [1  -  Sum_{n=1}^MOrb |phi_{n}><phi_{n}| ]

!  if LData = .true.  THEN prefactor = -1
!  if LData = .false. THEN prefactor = -i


IMPLICIT NONE
REAL*8,INTENT(IN)          :: AbsTime
COMPLEX*16, DIMENSION(NDX*NDY*NDZ*MOrb), INTENT(INOUT):: Psi_WORK
COMPLEX*16, DIMENSION(NDX*NDY*NDZ*MOrb), INTENT(OUT)  :: DtPsi_WORK
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb)               :: Psi
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb)               :: DtPsi
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb)               :: prefactor_hPsi,Chi
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb)               :: prefactor_W_jjj
COMPLEX*16, DIMENSION(NDX*NDY*NDZ)                    :: Wsl,Wsl_phi_q
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb)               :: projected_Chi
COMPLEX*16, INTENT(IN)     :: CData(MOrb,MOrb,MOrb,MOrb,2)
REAL*8                     :: RData(10)
INTEGER,INTENT(IN)         :: IData(30)
logical,INTENT(IN)         :: LData
logical                    :: Relax

INTEGER                    :: s,l,q,j,k,m,error,i,propDirection
COMPLEX*16                 :: rho_jk(MOrb,MOrb),rho_ijkl(MOrb,MOrb,MOrb,MOrb)
COMPLEX*16                 :: Inv_rho_jk(MOrb,MOrb), WorkMat(MOrb,MOrb)
COMPLEX*16                 :: dicht4(MOrb,MOrb),workr(3*MOrb)
REAL*8                     :: dicht3(MOrb,2)

COMPLEX*16                 :: prefactor
COMPLEX*16                 :: defect,check(MOrb),Ovlp
REAL*8                     :: norm


propDirection   = IData(1)
Relax           = LData
error           = 0
Chi             = (0.d0,0.d0)
projected_Chi   = (0.d0,0.d0)


Psi             = RESHAPE(Psi_WORK,(/NDX*NDY*NDZ,MOrb/))
DtPsi           = RESHAPE(DtPsi_WORK,(/NDX*NDY*NDZ,MOrb/))

!maybe remove this later, may not be neccessary (but check!)
CALL schmidtortho (Psi,NDX*NDY*NDZ,MOrb,error)
CALL schmidtortho (Psi,NDX*NDY*NDZ,MOrb,error)

IF (error.ne.0) THEN

  write(*,*)"schmidtortho failed",error
  STOP

END IF 

IF (Relax) THEN
   prefactor = propDirection * (-1.d0,0.d0)
ELSE 
   prefactor = propDirection * (0.d0,-1.d0)
END IF

!---------regularize the 1-particle density matrix if neccessary---------
!-- compute the inverse of the first order RDM 

rho_jk   = CData(:,:,1,1,1)
rho_ijkl = CData(:,:,:,:,2)

    error = 0
    WorkMat  = rho_jk
    CALL density(Inv_rho_jk,WorkMat,dicht3,dicht4,workr)

    IF (error .ne. 0) THEN
       write(*,*) 'rho_jk not positive definite error',error
       STOP
    END IF

!--- make sure that every orbital is normalized. It should be anyway!----
   DO j=1,MOrb

     CALL normvxz(Psi(:,j),norm,NDX*NDY*NDZ) 
!     write(*,*)j,"norm",norm
     IF (ABS(norm-1.d0)<1.d-14 ) THEN
        Psi(:,j) = Psi(:,j) / norm
     ELSE
        Psi(:,j) = Psi(:,j) / norm
        CALL normvxz(Psi(:,j),norm,NDX*NDY*NDZ) 
     END IF 

   END DO

CALL get_hPsi(MOrb,NDX,NDY,NDZ,Dims,Psi,prefactor_hPsi)
prefactor_hPsi=prefactor*prefactor_hPsi

!--------- |chi_j> = prefactor * h |phi_j> ------------------------
Chi   = prefactor_hPsi


     DO s=1,MOrb
       DO l=1,MOrb
         CALL get_Wsl(Wsl,s,l,Psi)
         DO q=1,MOrb
   
           CALL get_Wsl_phi_q(Wsl,Psi(:,q),Wsl_phi_q) 

           DO j=1,MOrb
             DO k=1,MOrb
   
                Chi(:,j) = Chi(:,j) + prefactor* Inv_rho_jk(j,k)*rho_ijkl(k,s,q,l)*Wsl_phi_q

             END DO  
           END DO

         END DO
       END DO
     END DO

   CALL projector_MCTDHB(Chi,DtPsi,Psi)
   Chi   = DtPsi                               !--DO another projection 
   CALL projector_MCTDHB(Chi,DtPsi,Psi)
   Chi   = DtPsi                               !--DO another projection 
   CALL projector_MCTDHB(Chi,DtPsi,Psi)
 
!---- test the orthogonality of Psi and DtPsi
DO j=1,MOrb

!  CALL normvxz(DtPsi(:,j),norm,NDX*NDY*NDZ) 
  DO k=1,MOrb 

     CALL vvaxzz (Psi(:,j),DtPsi(:,k),defect,NDX*NDY*NDZ)
     IF (ABS(defect)>1.d-10) THEN
     write(*,*)j,k,"<Psi_j|DtPsi_k>",defect
!     read(*,*)
     END IF 

   END DO

END DO

Psi_WORK        = RESHAPE(Psi,(/NDX*NDY*NDZ*MOrb/))
DtPsi_WORK      = RESHAPE(DtPsi,(/NDX*NDY*NDZ*MOrb/))


!write(*,*)"-----leave FuncOrbs-------------"
!---}}}
END SUBROUTINE FuncOrbs

!----------------------------------------------------------------------------------
SUBROUTINE FuncCVecBandMat(AbsTime,CVec,DtCVec,BandMat,RData,IData,LData)
USE   moduleparameters
USE   moduledefaults,ONLY:nConf,nSuperDiags
!---{{{
IMPLICIT NONE

COMPLEX*16,parameter    :: alpha=(1.d0,0.d0)
COMPLEX*16,parameter    :: beta=(0.d0,0.d0)
INTEGER,parameter       :: incx=1
INTEGER,parameter       :: incy=1

INTEGER,   INTENT(IN)   :: IData(30)
REAL*8,    INTENT(IN)   :: AbsTime,RData(10)
COMPLEX*16,INTENT(IN)   :: CVec(nConf),BandMat(nSuperDiags+1,nConf)
logical,   INTENT(IN)   :: LData
logical                 :: Relax
 
COMPLEX*16,INTENT(OUT)  :: DtCVec(nConf)
REAL*8                  :: PropDirection

IF (Abs(IData(1)).ne.1) THEN 
  write(*,*)"illegal propagation direction in FuncCVecBandMat",IData(1)
END IF

PropDirection   = IData(1) * 1.d0
Relax           = LData

CALL ZHBMV ( 'L', nConf, nSuperDiags, alpha, BandMat, nSuperDiags+1, CVec, incx, beta, DtCVec, incy )

IF(Relax) THEN
   DtCVec = DtCVec * (-1.d0,0.d0)  
ELSE
   DtCVec = DtCVec * (0.d0,-1.d0) 
END IF

DtCVec = PropDirection * DtCVec 

!---}}}
END SUBROUTINE FuncCVecBandMat

!----------------------------------------------------------------------------------
SUBROUTINE get_FullHamiltonianBandMatrix(NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,&
                                         nSuperDiags,lambda0,Psi,BandMat)
!---{{{
IMPLICIT NONE
INTEGER,INTENT(IN)      :: NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,nSuperDiags
REAL*8,INTENT(IN)       :: lambda0
COMPLEX*16,INTENT(IN)   :: Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16              :: h_jk(MOrb,MOrb)
COMPLEX*16,INTENT(OUT)  :: BandMat(nSuperDiags+1,nConf)
INTEGER                 :: m,n,n1,n2,L
REAL*8                  :: dn1,dn2
COMPLEX*16              :: W1112,W1222,W1122,W1212,W1221
COMPLEX*16              :: W1111,W2222,H_mn

BandMat = (0.d0,0.d0)
H_mn=(0.d0,0.d0)

IF((MOrb>2).or.(MOrb<1)) THEN
  WRITE(*,*) "SORRY, YOU WILL HAVE TO IMPLEMENT THIS IN get_FullHamiltonianBandMatrix, MOrb=",MOrb
  STOP
END IF

CALL get_h_jk(Psi,h_jk,NDX,NDY,NDZ,Dims,MOrb)
CALL get_Wijkl(W1111,1,1,1,1,Psi)

IF (MOrb==1) THEN

   BandMat(1,1) = DBLE(NPar*h_jk(1,1)  + 0.5d0*W1111*NPar*(NPar-1.d0) )

ELSE IF (MOrb==2) THEN

  CALL get_Wijkl(W2222,2,2,2,2,Psi)
  CALL get_Wijkl(W1112,1,1,1,2,Psi)
  CALL get_Wijkl(W1212,1,2,1,2,Psi)
  CALL get_Wijkl(W1221,1,2,2,1,Psi)
  CALL get_Wijkl(W1222,1,2,2,2,Psi)
  CALL get_Wijkl(W1122,1,1,2,2,Psi)

  DO m=1,nConf

     L = 1 - m

     DO n=m,MIN(nConf,m+nSuperDiags)
        !n1 and n2 are the occupation numbers of  |R> in <L|H|R>  
        n2  = n - 1 
        n1  = NPar - n2 
        dn1 = DBLE(n1) 
        dn2 = DBLE(n2) 
     
        IF(n==m) THEN

            BandMat(L+n,m) = DBLE(dn1*h_jk(1,1) + dn2*h_jk(2,2) &
            + 0.5d0*W1111*dn1*(dn1-1.d0)     &
            + 0.5d0*W2222*dn2*(dn2-1.d0)     &
            + dn1*dn2 * (W1212+W1221))

        ELSE IF (n==(m+1)) THEN
                     
            H_mn  = h_jk(1,2)*DSQRT(1.d0 * (dn1+1.d0) * dn2)   &
                   + dn1 * DSQRT( dn2*(dn1+1.d0) ) * W1112            &
                   + (dn2-1.d0)* DSQRT( (dn1+1.d0)*dn2 )  * W1222
            BandMat(L+n,m) =  DCONJG(H_mn)

        ELSE IF (n==(m+2)) THEN
        
            H_mn   =  0.5d0*DSQRT(dn2*(dn2-1.d0)*(dn1+1.d0)*(dn1+2.d0) ) * W1122
            BandMat(L+n,m) =  DCONJG(H_mn)

        ELSE 

          write(*,*)"I am out of bounds m,n:",m,n
          STOP

        END IF

     END DO 
  END DO 

END IF


!---}}}
END SUBROUTINE get_FullHamiltonianBandMatrix

!----------------------------------------------------------------------------------
SUBROUTINE projector_MCTDHB(Chi,projectedChi,Psi)
USE   moduleinputvariables,ONLY:MOrb,NDX,NDY,NDZ
!---{{{

IMPLICIT NONE
COMPLEX*16,INTENT(IN)     :: Psi(NDX*NDY*NDZ,MOrb),Chi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(OUT)    :: projectedChi(NDX*NDY*NDZ,MOrb)
COMPLEX*16                :: Eta(NDX*NDY*NDZ)
COMPLEX*16                :: Chi_new(NDX*NDY*NDZ,MOrb)
COMPLEX*16                :: check(MOrb),O_jj(MOrb)
REAL*8                    :: norm(MOrb)
INTEGER                   :: j,n,k

norm            =       0.d0
check           =       (0.d0,0.d0)
projectedChi    =       (0.d0,0.d0)
O_jj            =       (0.d0,0.d0)

!----check whether the Psi_k functions are normalized
DO k=1,MOrb

   CALL normvxz(Psi(:,k),norm(k),NDX*NDY*NDZ) 
   IF (ABS(norm(k)-1.d0)>1.0d-14) THEN 
      write(*,*)k,"-th phi norm:",norm(k)
   END IF 

END DO

projectedChi    =       Chi 

DO k=1,MOrb

   Eta             =       (0.d0,0.d0)
   DO j=1,MOrb
     
      CALL vvaxzz (Psi(:,j),Chi(:,k),O_jj(j),NDX*NDY*NDZ)
      Eta = Eta + Psi(:,j)*O_jj(j)

   END DO 
   
   projectedChi(:,k) = projectedChi(:,k) - Eta 

END DO


!---}}}
END SUBROUTINE projector_MCTDHB


!----------------------------------------------------------------------------------
SUBROUTINE get_hPsi(MOrb,NDX,NDY,NDZ,Dims,Psi,hPsi)
USE moduleallocatables,ONLY:p2over2m,V_ext
USE modulefft
!---{{{
IMPLICIT NONE

intEGER,INTENT(IN)              :: MOrb,NDX,NDY,NDZ,Dims
COMPLEX*16,INTENT(IN)           :: Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(OUT)          :: hPsi(NDX*NDY*NDZ,MOrb)

INTEGER                         :: m,n,INFO,noGridPoints
COMPLEX*16,DIMENSION(NDX*NDY*NDZ,MOrb) :: chi,phi
REAL*8                          :: norm

noGridPoints = NDX*NDY*NDZ
INFO = 0
!--------------------------------------------------

!---compute |chi> = p^2/(2m)|Psi>--------{{{

chi     = Psi 

if (Dims==1) then 
!--- {{{ 1d 
  if(noGridPoints.ne.NDX) then
    write(*,*)"this is not good 1d" 
    STOP
  endif
  
  DO m=1,MOrb

    CALL get_FFT1D(-1,NDX,chi(:,m),INFO)

    DO n=1,noGridPoints
      chi(n,m)= chi(n,m) * p2over2m(n) ! multiply by (p**2)/2m (where m=1)
    END DO    
    chi(:,m)=chi(:,m)/(1.d0*noGridPoints)  

    CALL get_FFT1D(1,NDX,chi(:,m),INFO)

  END DO 
!---}}}

ELSE IF (Dims==2) THEN

!--- {{{ 2d 
  IF(noGridPoints.ne.NDX*NDY) THEN
    write(6,*)"this is not good 2d" 
    STOP
  END IF

  DO m=1,MOrb

    CALL get_FFT2D(-1,NDX,NDY,chi(:,m),INFO)

    DO n=1,noGridPoints
      chi(n,m)= chi(n,m) * p2over2m(n) ! multiply by (p**2)/2m (where m=1)
    END DO    
    chi(:,m)=chi(:,m)/(1.d0*noGridPoints)  

    CALL get_FFT2D(1,NDX,NDY,chi(:,m),INFO)

  END DO 
!---}}}

else if (Dims==3) then

!--- {{{ 3d 
  if(noGridPoints.ne.NDX*NDY*NDZ) then
    write(*,*)"this is not good 3d" 
    STOP
  endif

  DO m=1,MOrb

    CALL get_FFT3D(-1,NDX,NDY,NDZ,chi(:,m),INFO)

    DO n=1,noGridPoints
      chi(n,m)= chi(n,m) * p2over2m(n) ! multiply by (p**2)/2m (where m=1)
    END DO    
    chi(:,m)=chi(:,m)/(noGridPoints*1.d0)  

    CALL get_FFT3D(1,NDX,NDY,NDZ,chi(:,m),INFO)


  END DO 

!---}}}

else
  write(6,*)"more than 3d?"
  STOP
endif

!}}}

!---compute |phi> = V|Psi> 
do m=1,MOrb
  do n=1,noGridPoints
    phi(n,m)=Psi(n,m)*V_ext(n)
  end do 
end do 

hPsi = chi + phi


!---}}}
END SUBROUTINE get_hPsi


!--------------------------------------------------------------------
SUBROUTINE get_Wsl_IMEST(Wsl,psi_s,psi_l)
USE moduleinputvariables,ONLY: NDX,NDY,NDZ,Dims
USE moduleallocatables, ONLY: vtilde
USE modulefft
!---{{{
IMPLICIT NONE
COMPLEX*16,INTENT(out)  :: Wsl(NDX*NDY*NDZ)
COMPLEX*16,INTENT(in)   :: psi_s(NDX*NDY*NDZ),psi_l(NDX*NDY*NDZ)

INTEGER                 :: i,j,k,m,n
COMPLEX*16              :: fsltilde(NDX*NDY*NDZ)

!---------------FFT-------------------------------
INTEGER                 :: info,noGridPoints

noGridPoints = NDX*NDY*NDZ

!compute f_sl = dx phi_s^* phi_l
DO m=1,noGridPoints
  fsltilde(m)=DCONJG(psi_s(m))*psi_l(m)
END DO

!get f_sl-tilde dx from f_sl dx
IF (Dims==1) THEN 
!---{{{ 1d 
  IF(noGridPoints.ne.NDX) THEN
    write(6,*)"this is not good wsl" 
    stop
  END IF

  CALL get_FFT1D(-1,NDX,fsltilde,INFO)
  fsltilde = 1.d0/DSQRT(1.d0*noGridPoints) * fsltilde

!---}}}

ELSE IF (Dims==2) THEN
!--- {{{ 2d 
  IF(noGridPoints.ne.NDX*NDY) THEN
    write(6,*)"this is not good 2d" 
    stop
  END IF

  CALL get_FFT2D(-1,NDX,NDY,fsltilde,INFO)
  fsltilde = 1.d0/DSQRT(1.d0*noGridPoints) * fsltilde

!---}}}

ELSE IF (Dims==3) THEN
!--- {{{ 3d 
  IF(noGridPoints.ne.NDX*NDY*NDZ) THEN
    write(6,*)"this is not good 3d" 
    stop
  END IF

  CALL get_FFT3D(-1,NDX,NDY,NDZ,fsltilde,INFO)
  fsltilde = 1.d0/DSQRT(1.d0*noGridPoints) * fsltilde
      
!---}}}
END IF

!compute W_sl-tilde(m) = f_sl-tilde(m)*v-tilde(m)*dx
DO m=1,noGridPoints
  Wsl(m)=fsltilde(m)*vtilde(m)
END DO

!W_sl(i)=F^+(Wsl-tilde)
IF (Dims==1) THEN 
!--- {{{ 1d 
  IF(noGridPoints.ne.NDX) THEN
    write(6,*)"this is not good wsl" 
    stop
  END IF

  CALL get_FFT1D(1,NDX,Wsl,INFO)
  Wsl = 1.d0/DSQRT(1.d0*noGridPoints) * Wsl

!---}}}

ELSE IF (Dims==2) THEN
!--- {{{ 2d 
  IF(noGridPoints.ne.NDX*NDY) THEN
    write(6,*)"this is not good 2d" 
    stop
  END IF

  CALL get_FFT2D(1,NDX,NDY,Wsl,INFO)
  Wsl = 1.d0/DSQRT(1.d0*noGridPoints) * Wsl

!---}}}

ELSE IF (Dims==3) THEN
!--- {{{ 3d 
  IF(noGridPoints.ne.NDX*NDY*NDZ) THEN
    write(6,*)"this is not good 3d" 
    stop
  END IF

  CALL get_FFT3D(1,NDX,NDY,NDZ,Wsl,INFO)
  Wsl = 1.d0/DSQRT(1.d0*noGridPoints) * Wsl

!---}}}
END IF

!---}}}
END SUBROUTINE get_Wsl_IMEST

!--------------------------------------------------------------------
SUBROUTINE get_vtilde(NDX,NDY,NDZ,Dims,xi,xf,yi,yf,zi,zf,vtilde)
USE moduleparameters
USE moduleextensions
USE modulegrid
USE modulefft
!---{{{
IMPLICIT NONE

INTEGER,INTENT(IN)              :: NDX
INTEGER,INTENT(IN)              :: NDY
INTEGER,INTENT(IN)              :: NDZ
INTEGER,INTENT(IN)              :: Dims
REAL*8,INTENT(IN)               :: xi
REAL*8,INTENT(IN)               :: xf
REAL*8,INTENT(IN)               :: yi
REAL*8,INTENT(IN)               :: yf
REAL*8,INTENT(IN)               :: zi
REAL*8,INTENT(IN)               :: zf
COMPLEX*16, INTENT(OUT)         :: vtilde(NDX*NDY*NDZ)

COMPLEX*16                      :: wtilde(NDX*NDY*NDZ)
REAL*8                          :: x,y,z,xall,yall,zall,dx,dy,dz
INTEGER                         :: m,i,j,k

INTEGER                         :: info,noGridPoints
  
  noGridPoints=NDX*NDY*NDZ
  xall  = xf-xi
  yall  = yf-yi
  zall  = zf-zi
  dx    = xall/(NDX*1.d0)
  dy    = yall/(NDY*1.d0)
  dz    = zall/(NDZ*1.d0)

  DO m=1,noGridPoints

    call get_ijk_from_m(m,NDX,NDY,i,j,k)
    x = (i-1)*dx-(xall/2.d0)
    y = (j-1)*dy-(yall/2.d0)
    z = (k-1)*dz-(zall/2.d0)

    IF (Dims==1) THEN
      wtilde(m) = Wint_IMEST(x,0.d0,0.d0)
    ELSE IF (Dims==2) THEN
      wtilde(m) = Wint_IMEST(x,y,0.d0)
    ELSE IF (Dims==3) THEN
      wtilde(m) = Wint_IMEST(x,y,z)
    ELSE 
      write(6,*)"more than 3d?"
      stop
    END IF

  END DO    

IF (Dims==1) THEN 
!--- {{{ 1d 
  IF(noGridPoints.ne.NDX) THEN
    write(6,*)"this is not good" 
    stop
  END IF

  CALL get_FFT1D(-1,NDX,wtilde,INFO)
  IF (info.ne.0) write(6,*) "info get_FFT1D=",info
 
  DO i=1,NDX
    vtilde(i) = (-1.d0)**(i-1)*wtilde(i) 
  END DO

!---}}}
ELSE IF (Dims==2) THEN
!---{{{ 2d
  IF(noGridPoints.ne.NDX*NDY) THEN
    write(6,*)"this is not good" 
    stop
  END IF

  CALL get_FFT2D(-1,NDX,NDY,wtilde,INFO)
  IF (info.ne.0) write(6,*) "info get_FFT1D=",info

  DO m=1,noGridPoints
    call get_ijk_from_m(m,NDX,NDY,i,j,k)
    vtilde(m) = (-1.d0)**(i-1+j-1)*wtilde(m) 
  END DO


!---}}}
ELSE IF (Dims==3) THEN
!---{{{
  IF(noGridPoints.ne.NDX*NDY*NDZ) THEN
    write(6,*)"this is not good ..." 
    stop
  END IF

  CALL get_FFT3D(-1,NDX,NDY,NDZ,wtilde,INFO)
  IF (info.ne.0) write(6,*) "info get_FFT1D=",info

  DO m=1,noGridPoints
    call get_ijk_from_m(m,NDX,NDY,i,j,k)
    vtilde(m) = (-1.d0)**(i-1 + j-1 + k-1) * wtilde(m) 
!    write(6,*)m,"vtilde(m)",vtilde(m)
  END DO

!---}}}
ELSE
  write(6,*)"more than 3d?"
  stop
END IF

!---}}}
END SUBROUTINE get_vtilde

!----------------------------------------------------------------------------------
SUBROUTINE get_h_jk(Psi,h_jk,NDX,NDY,NDZ,Dims,MOrb)
!---{{{ 
IMPLICIT NONE
COMPLEX*16,INTENT(IN)   :: Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(OUT)  :: h_jk(MOrb,MOrb)
INTEGER,INTENT(IN)      :: MOrb,NDX,NDY,NDZ,Dims
COMPLEX*16              :: hPsi(NDX*NDY*NDZ,MOrb)
INTEGER                 :: i,j


hPsi=(0.0d0,0.d0)
CALL get_hPsi(MOrb,NDX,NDY,NDZ,Dims,Psi,hPsi)

DO i=1,MOrb
  DO j=i,MOrb

   CALL vvaxzz (Psi(:,i),hPsi(:,j),h_jk(i,j),NDX*NDY*NDZ) 
   h_jk(j,i)=DCONJG(h_jk(i,j))

  END DO
  h_jk(i,i)=DBLE(h_jk(i,i))
END DO

!---}}}
END SUBROUTINE get_h_jk

!----------------------------------------------------------------------------------
SUBROUTINE get_Wsl(Wsl,s,l,Psi)
USE moduleinputvariables,ONLY:lambda0,NDX,NDY,NDZ,MOrb,useIMEST
USE modulederivesystemvariables,ONLY:weight
!---{{{ 
IMPLICIT NONE
INTEGER,INTENT(IN)          :: s,l
COMPLEX*16,INTENT(IN)       :: Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(OUT)      :: Wsl(NDX*NDY*NDZ)
COMPLEX*16                  :: PsiS(NDX*NDY*NDZ),PsiL(NDX*NDY*NDZ)
INTEGER                     :: m

PsiS=Psi(:,s)
PsiL=Psi(:,l)
Wsl=(0.d0,0.d0)

IF (useIMEST) THEN

  CALL get_Wsl_IMEST(Wsl,PsiS,PsiL)

ELSE

DO m=1,NDX*NDY*NDZ
   Wsl(m) = lambda0 * DCONJG(PsiS(m)/weight) * PsiL(m)/weight 
END DO

END IF

!---}}}
END SUBROUTINE get_Wsl

!----------------------------------------------------------------------------------
SUBROUTINE get_Wsl_phi_q(Wsl,phi_q,Wsl_phi_q) 
USE moduleinputvariables,ONLY:NDX,NDY,NDZ
!--- {{{
IMPLICIT NONE 
COMPLEX*16,INTENT(IN)       :: Wsl(NDX*NDY*NDZ),phi_q(NDX*NDY*NDZ)
COMPLEX*16,INTENT(OUT)      :: Wsl_phi_q(NDX*NDY*NDZ)
INTEGER                     :: i

Wsl_phi_q=(0.d0,0.d0)
DO i=1,NDX*NDY*NDZ

   Wsl_phi_q(i) = Wsl(i)*phi_q(i)

END DO 

!}}}
END SUBROUTINE get_Wsl_phi_q

!----------------------------------------------------------------------------------
SUBROUTINE get_Wijkl(Wijkl,i,j,k,l,Psi)
USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb
!---compute Wijkl {{{
IMPLICIT NONE
INTEGER,INTENT(IN)          :: i,j,k,l
COMPLEX*16,INTENT(IN)       :: Psi(NDX*NDY*NDZ,MOrb)
COMPLEX*16,INTENT(OUT)      :: Wijkl

COMPLEX*16                  :: PsiI(NDX*NDY*NDZ),PsiK(NDX*NDY*NDZ)
COMPLEX*16                  :: Wjl(NDX*NDY*NDZ)
INTEGER                     :: m

PsiI = Psi(:,i)
PsiK = Psi(:,k)
Wijkl=(0.d0,0.d0)

CALL get_Wsl(Wjl,j,l,Psi)

DO m=1,NDX*NDY*NDZ
   Wijkl = Wijkl + DCONJG(PsiI(m)) *  Wjl(m) * PsiK(m)
END DO


!---}}}
END SUBROUTINE get_Wijkl

!----------------------------------------------------------------------------------
SUBROUTINE get_rho_jk(nConf,MOrb,NPar,CVec,rho_jk)
!---{{{ 
IMPLICIT NONE
INTEGER,INTENT(IN)       :: nConf,MOrb,NPar
COMPLEX*16,INTENT(IN)    :: CVec(nConf)
COMPLEX*16,INTENT(OUT)   :: rho_jk(MOrb,MOrb)
INTEGER                  :: i,j,n,n1,n2

rho_jk = (0.d0,0.d0)

orbitalif: IF (MOrb == 1) THEN 

  rho_jk(1,1) = 1.0d0*NPar

ELSE IF (MOrb == 2) THEN 

  DO n=1,nConf
      
    n1 = nConf-n             
    n2 = NPar-n1            

    rho_jk(1,1) = rho_jk(1,1) + DCONJG(CVec(n))*CVec(n) * n1
    rho_jk(2,2) = rho_jk(2,2) + DCONJG(CVec(n))*CVec(n) * n2

  END DO

  DO n=1,nConf-1

    n1 = nConf-n             
    n2 = NPar-n1            

    rho_jk(1,2) = rho_jk(1,2) + DCONJG(CVec(n))*CVec(n+1) * DSQRT(1.D0*n1*(n2+1)) 
    rho_jk(2,1) = DCONJG(rho_jk(1,2))

  END DO  

ELSE orbitalif 

  WRITE(*,*)"INVALID NUMBER OF ORBITALS, MOrb = ",MOrb        
  WRITE(*,*)"MORE THAN M=2 ORBITALS ARE NOT IMPLEMENTED IN get_rho_jk"        
  STOP

END IF orbitalif


!---}}}
END SUBROUTINE get_rho_jk

!----------------------------------------------------------------------------------
SUBROUTINE get_rho_ijkl(nConf,MOrb,NPar,CVec,rho_ijkl)
USE moduleparameters      
!---{{{ Evaluate the two-body reduced density matrix elements rho_ijkl for MOrb = 1,2 ONLY
IMPLICIT NONE
INTEGER,INTENT(IN)       :: nConf,MOrb,NPar
COMPLEX*16,INTENT(IN)    :: CVec(nConf)
COMPLEX*16,INTENT(OUT)   :: rho_ijkl(MOrb,MOrb,MOrb,MOrb)
INTEGER                  :: i,j,k,l,n,n1,n2

rho_ijkl = ZERO

orbitalif: IF (MOrb == 1) THEN 

  rho_ijkl(1,1,1,1) = 1.0d0*NPar * (NPar-1.d0)

ELSE IF (MOrb == 2) THEN 

  DO n=1,nConf
      
    n1 = nConf-n             
    n2 = NPar-n1            

! rho_kkkk
    rho_ijkl(1,1,1,1) = rho_ijkl(1,1,1,1) + DCONJG(CVec(n))*CVec(n) * n1 * (n1-1.d0)  
    rho_ijkl(2,2,2,2) = rho_ijkl(2,2,2,2) + DCONJG(CVec(n))*CVec(n) * n2 * (n2-1.d0)  

! rho_ksks
    rho_ijkl(1,2,1,2) = rho_ijkl(1,2,1,2) + DCONJG(CVec(n))*CVec(n) * n1 * n2  

  END DO

  rho_ijkl(1,2,2,1) = rho_ijkl(1,2,1,2)
  rho_ijkl(2,1,2,1) = rho_ijkl(1,2,1,2)
  rho_ijkl(2,1,1,2) = rho_ijkl(1,2,1,2)

  DO n=1,nConf-1

    n1 = nConf-n             
    n2 = NPar-n1            
    rho_ijkl(1,1,1,2) = rho_ijkl(1,1,1,2) + DCONJG(CVec(n))*CVec(n+1) * (n1-1.d0) * DSQRT(1.D0*n1*(n2+1)) 
    rho_ijkl(1,2,2,2) = rho_ijkl(1,2,2,2) + DCONJG(CVec(n))*CVec(n+1) * n2 * DSQRT(1.D0*n1*(n2+1))
    
  END DO  

  rho_ijkl(1,1,2,1) = rho_ijkl(1,1,1,2)
  rho_ijkl(1,2,1,1) = DCONJG(rho_ijkl(1,1,1,2))
  rho_ijkl(2,1,1,1) = DCONJG(rho_ijkl(1,1,1,2))  
  rho_ijkl(2,1,2,2) = rho_ijkl(1,2,2,2)
  rho_ijkl(2,2,2,1) = DCONJG(rho_ijkl(1,2,2,2))
  rho_ijkl(2,2,1,2) = DCONJG(rho_ijkl(1,2,2,2))

  DO n=1,nConf-2

    n1 = nConf-n             
    n2 = NPar-n1            
    rho_ijkl(1,1,2,2) = rho_ijkl(1,1,2,2) + DCONJG(CVec(n))*CVec(n+2) * DSQRT((n1-1.d0) * n1 *(n2+1.d0) * (n2+2.d0)) 
    
  END DO  

  rho_ijkl(2,2,1,1) = DCONJG(rho_ijkl(1,1,2,2))

ELSE orbitalif 

  WRITE(*,*)"INVALID NUMBER OF ORBITALS, MOrb = ",MOrb        
  WRITE(*,*)"MORE THAN M=2 ORBITALS ARE NOT IMPLEMENTED IN get_rho_ijkl"        
  STOP

END IF orbitalif

!       }}}
END SUBROUTINE get_rho_ijkl


!----------------------------------------------------------------------------------
SUBROUTINE get_mu_kj(AbsTime,mu_kj,Psi,rho_jk,rho_ijkl)
USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb,Dims,lambda0,NPar
!---{{{  evaluates 
!
!  mu_kj = Sum_{q=1}^M { rho_kq h_jq + Sum_{s,l=0}^M rho_ksql W_jsql }
!  for k,j=1,...,MOrb

IMPLICIT NONE
REAL*8,INTENT(IN)          :: AbsTime
COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MOrb),INTENT(IN)    :: Psi
COMPLEX*16, INTENT(IN)     :: rho_jk(MOrb,MOrb),rho_ijkl(MOrb,MOrb,MOrb,MOrb)
COMPLEX*16, INTENT(OUT)    :: mu_kj(MOrb,MOrb)
INTEGER                    :: s,l,q,j,k

COMPLEX*16                 :: h_jk(MOrb,MOrb)
COMPLEX*16                 :: Wjsql

mu_kj           = (0.d0,0.d0)
Wjsql           = (0.d0,0.d0)

CALL get_h_jk(Psi,h_jk,NDX,NDY,NDZ,Dims,MOrb)

DO k=1,MOrb
  DO j=1,MOrb
    DO q=1,MOrb
      DO s=1,MOrb
        DO l=1,MOrb

           CALL get_Wijkl(Wjsql,j,s,q,l,Psi)
           mu_kj(k,j)=mu_kj(k,j) + rho_ijkl(k,s,q,l)*Wjsql

        END DO
      END DO
  
      mu_kj(k,j)=mu_kj(k,j) + rho_jk(k,q) * h_jk(j,q)

    END DO
  END DO
END DO

mu_kj = mu_kj/Npar


!---}}}
END SUBROUTINE  get_mu_kj





END MODULE modulecomputationcore