|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 ! vim:fdm=marker: 00002 00009 MODULE moduleintegration 00010 IMPLICIT NONE 00011 00012 CONTAINS 00019 SUBROUTINE do_stepOne(nConf,MOrb,NPar,CVec_zero,CData,readPotential,V_ext,& 00020 NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime,nSuperDiags,lambda0,& 00021 Psi_zero,BandMat,rho_jk_zero,rho_ijkl_zero) 00022 00023 00024 USE modulecomputationcore 00025 USE moduleextensions 00026 00027 INTEGER,INTENT(IN) :: NDX,NDY,NDZ,MOrb,nConf,nSuperDiags,NPar,Dims 00028 LOGICAL,INTENT(IN) :: ReadPotential 00029 REAL*8,INTENT(IN) :: lambda0,gridx(NDX),gridy(NDY),gridz(NDZ),AbsTime 00030 REAL*8,INTENT(INOUT) :: V_ext(NDX*NDY*NDZ) 00031 COMPLEX*16,INTENT(INOUT) :: Psi_Zero(NDX*NDY*NDZ,MOrb),CData(MOrb,MOrb,MOrb,MOrb,2) 00032 COMPLEX*16,INTENT(INOUT) :: rho_ijkl_zero(MOrb,MOrb,MOrb,MOrb),rho_jk_zero(MOrb,MOrb) 00033 00034 00035 COMPLEX*16,INTENT(INOUT) :: BandMat(nSuperDiags+1,nConf),CVec_zero(nConf) 00036 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00037 !+++++++++++++++++++++++++++++++++++ STEP 1 +++++++++++++++++++++++++++++++++++++++++ 00038 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00039 ! Given the orbitals Psi at t=AbsTime and the coefficient vector 00040 ! CVec(AbsTime) compute rho_ijkl and build the Hamiltonian matrix 00041 ! write(*,*)"AbsTime",AbsTime 00042 ! write(*,*)"++++++++ STEP 1 +++++++++++" 00043 00044 CALL get_rho_jk(nConf,MOrb,NPar,CVec_zero,rho_jk_zero) 00045 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_zero,rho_ijkl_zero) 00046 CData(:,:,1,1,1) = rho_jk_zero 00047 CData(:,:,:,:,2) = rho_ijkl_zero 00048 00049 ! If a time-dependent external potential is used, it needs to be updated here 00050 IF (.NOT. readPotential) THEN 00051 CALL get_V_ext(V_ext,NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime) 00052 END IF 00053 00054 ! Build the CI matrix H(AbsTime) 00055 CALL get_FullHamiltonianBandMatrix(NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,& 00056 nSuperDiags,lambda0,Psi_zero,BandMat) 00057 END SUBROUTINE do_stepOne 00058 00065 SUBROUTINE do_stepTwo(AbsTime,CVec_zero,DtCVec,BandMat,RData,IData,LData,& 00066 CVec_tauhalf,NPar,nConf,MOrb,IntPeriod,SILmaxIntOrder,& 00067 IntegratorTolError,Relax,CData,tau, & 00068 SILIntPeriod1,nSuperDiags,& 00069 rho_ijkl_tauhalf,rho_jk_tauhalf,propDirection) 00070 00071 00072 USE modulecomputationcore 00073 00074 ! EXTERNAL FuncCVecBandMat 00075 00076 INTEGER, INTENT(IN) :: nSuperDiags,SILmaxIntOrder,MOrb,nConf,NPar,propDirection 00077 COMPLEX*16 :: Krylov(nConf,2:SILmaxIntOrder) 00078 COMPLEX*16,INTENT(INOUT) :: CVec_tauhalf(nConf),DtCVec(nConf) 00079 COMPLEX*16,INTENT(INOUT) :: CData(MOrb,MOrb,MOrb,MOrb,2) 00080 COMPLEX*16,INTENT(INOUT) :: BandMat(nSuperDiags+1,nConf),CVec_zero(nConf) 00081 COMPLEX*16,INTENT(OUT) :: rho_ijkl_tauhalf(MOrb,MOrb,MOrb,MOrb),rho_jk_tauhalf(MOrb,MOrb) 00082 00083 00084 REAL*8, INTENT(INOUT) :: AbsTime,SILIntPeriod1,IntPeriod 00085 REAL*8, INTENT(INOUT) :: RData(10),tau 00086 INTEGER*4, INTENT(INOUT) :: IData(30) 00087 REAL*8, INTENT(IN) :: IntegratorTolError 00088 logical, INTENT(INOUT) :: LData 00089 logical, INTENT(IN) :: Relax 00090 logical :: RestartSIL,StdForm 00091 INTEGER :: TrueOrder,ErrorCode,Steps 00092 REAL*8 :: EigenVector(0:SILmaxIntOrder,0:SILmaxIntOrder+2), 00093 Diagonal(0:SILmaxIntOrder),OffDiag(SILmaxIntOrder), 00094 EigenVal(0:SILmaxIntOrder),OffDg2(SILmaxIntOrder+1), 00095 norm 00096 00097 !------------- INITIALIZE VARIABLES-------------------- 00098 EigenVector = 0.d0 00099 Diagonal = 0.d0 00100 OffDiag = 0.d0 00101 EigenVal = 0.d0 00102 OffDg2 = 0.d0 00103 ErrorCode=0 00104 TrueOrder=0 00105 Steps=0 00106 RestartSIL=.TRUE. 00107 StdForm = .TRUE. 00108 00109 00110 ! get DtCVec for SILStep 00111 IData(1) = propDirection 00112 CALL FuncCVecBandMat(AbsTime,CVec_zero,DtCVec,BandMat,RData,IData,LData) 00113 00114 ! integrate t = 0 -> tau/2 00115 CVec_tauhalf = CVec_zero 00116 00117 IntPeriod = tau / 2.d0 00118 CALL SILStep(CVec_tauhalf,DtCVec,nConf,IntPeriod,SILmaxIntOrder, & 00119 IntegratorTolError,Relax,RestartSIL,StdForm,Steps, & 00120 Krylov,SILIntPeriod1, & 00121 TrueOrder,ErrorCode,AbsTime,FuncCVecBandMat, & 00122 EigenVector,EigenVal,Diagonal,OffDg2,OffDiag, & 00123 BandMat,RData,IData,LData) 00124 00125 00126 if (ErrorCode.ne.0) then 00127 write(*,*)"++++++++ STEP 2 +++++++++++" 00128 write(*,*)"SILStep ErrorCode =",ErrorCode 00129 STOP 00130 end if 00131 00132 ! the integrated time interval may be shorter than IntPeriod if so, adjust tau 00133 IntPeriod = SILIntPeriod1 00134 tau = 2.d0 * SILIntPeriod1 00135 00136 ! normalize CVec to get rid of numerical artefacts 00137 CALL normvxz(CVec_tauhalf,norm,nConf) 00138 CVec_tauhalf=CVec_tauhalf/norm 00139 00140 00141 ! The vector CVec_tauhalf contains CVec(tau/2) where tau/2 = Stepsize 00142 ! Note that Stepsize can be less than IntPeriod! 00143 ! now get rho_jk(tau/2) and rho_jjjj(tau/2) using CVec(tau/2) 00144 ! and update CData 00145 00146 CALL get_rho_jk(nConf,MOrb,NPar,CVec_tauhalf,rho_jk_tauhalf) 00147 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_tauhalf,rho_ijkl_tauhalf) 00148 CData(:,:,1,1,1) = rho_jk_tauhalf 00149 CData(:,:,:,:,2) = rho_ijkl_tauhalf 00150 00151 00152 00153 END SUBROUTINE do_stepTwo 00154 00161 SUBROUTINE do_stepThree(Psi_tauhalf,Psi_zero,Psi_WORK,AbsTime,DtPsi_WORK,& 00162 CData,RData,IData,LData,NDX,NDY,NDZ,MOrb,Relax,tau,& 00163 IntegratorTolError,restartABM,SILIntPeriod1,& 00164 ABMIntOrder) 00165 00166 USE modulecomputationcore 00167 00168 EXTERNAL AbsABMError 00169 00170 COMPLEX*16,INTENT(INOUT) :: Psi_Zero(NDX*NDY*NDZ,MOrb),CData(MOrb,MOrb,MOrb,MOrb,2) 00171 COMPLEX*16,INTENT(INOUT) :: Psi_tauhalf(NDX*NDY*NDZ,MOrb) 00172 00173 00174 REAL*8, INTENT(INOUT) :: RData(10),SILIntPeriod1,tau,AbsTime,IntegratorTolError 00175 INTEGER*4, INTENT(INOUT) :: IData(30) 00176 LOGICAL, INTENT(INOUT) :: LData 00177 LOGICAL, INTENT(IN) :: Relax,restartABM 00178 00179 INTEGER, INTENT(IN) :: NDX,NDY,NDZ,ABMIntOrder,MOrb 00180 00181 COMPLEX*16 :: Psi_WORK(NDX*NDY*NDZ*MOrb), AuxPsi(NDX*NDY*NDZ*MOrb,ABMIntOrder+1) 00182 COMPLEX*16 :: DtPsi_WORK(NDX*NDY*NDZ*MOrb) 00183 INTEGER :: Steps,totalNoGridPts,RepeatedSteps,ErrorCode,k 00184 REAL*8 :: InitStep,ABMIntPeriod,norm 00185 00186 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00187 !+++++++++++++++++++++++++++++++++++ STEP 3 +++++++++++++++++++++++++++++++++++++++++ 00188 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00189 ! Propagate Psi(0)-> Psi(tau/2) 00190 ! using rho_jk(tau/2) and rho_ijkl(tau/2) with the ABM integrator 00191 ! write(*,*)"++++++++ STEP 3 +++++++++++" 00192 totalNoGridPts=NDX*NDY*NDZ 00193 Psi_tauhalf = Psi_zero 00194 Psi_WORK = RESHAPE(Psi_tauhalf,(/NDX*NDY*NDZ*MOrb/)) 00195 00196 CALL FuncOrbs(AbsTime,Psi_WORK,DtPsi_WORK,CData,RData,IData,LData) ! get dPsi/dt 00197 00198 ! do the Relaxation/propagation of the orbitals 00199 if (Relax) then 00200 InitStep = tau 00201 ABMIntPeriod = tau 00202 else 00203 InitStep = SILIntPeriod1 00204 ABMIntPeriod = SILIntPeriod1 00205 end if 00206 00207 CALL ABM(Psi_WORK,DtPsi_WORK,totalNoGridPts*MOrb,ABMIntPeriod,AbsTime,ABMIntOrder, & 00208 InitStep,IntegratorTolError,restartABM,Steps,RepeatedSteps, & 00209 ErrorCode,AuxPsi,FuncOrbs,AbsABMError,CData, & 00210 RData,IData,LData) 00211 if(ErrorCode.ne.0) then 00212 !---{{{ error output 00213 write(*,*)" ERROR OUTPUT" 00214 write(6,*)" totalNoGridPts = ",totalNoGridPts 00215 write(6,*)" ABMIntPeriod = ",ABMIntPeriod 00216 write(6,*)" AbsTime = ",AbsTime 00217 write(6,*)" ABMIntOrder = ",ABMIntOrder 00218 write(6,*)" TolError = ", 2.d0*IntegratorTolError 00219 write(6,*)" restartABM = ",restartABM 00220 write(6,*)" ABM Steps = ",Steps 00221 write(6,*)" ABM RepeatedSteps = ",RepeatedSteps 00222 write(6,*)" ABM ErrorCode = ",ErrorCode 00223 read(*,*) 00224 !----}}} 00225 end if 00226 Psi_tauhalf = RESHAPE(Psi_WORK,(/NDX*NDY*NDZ,MOrb/)) ! Psi_WORK -> Psi_tauhalf 00227 00228 do k=1,MOrb 00229 00230 CALL normvxz(Psi_tauhalf(:,k),norm,NDX*NDY*NDZ) 00231 Psi_tauhalf(:,k) = Psi_tauhalf(:,k) / norm 00232 00233 end do 00234 00235 00236 END SUBROUTINE do_stepThree 00243 SUBROUTINE do_stepFour(rho_jk_zero,rho_ijkl_zero,CData,Rdata,IData,LData,& 00244 Psi_tildetauhalf,Psi_zero,restartABM,SILIntPeriod1,& 00245 AbsTime,IntegratorTolError,NDX,NDY,NDZ,MOrb,ABMIntOrder) 00246 00247 USE modulecomputationcore 00248 00249 COMPLEX*16,INTENT(INOUT) :: Psi_tildetauhalf(NDX*NDY*NDZ,MOrb) 00250 COMPLEX*16,INTENT(INOUT) :: Psi_zero(NDX*NDY*NDZ,MOrb) 00251 COMPLEX*16,INTENT(INOUT) :: rho_ijkl_zero(MOrb,MOrb,MOrb,MOrb),rho_jk_zero(MOrb,MOrb) 00252 COMPLEX*16,INTENT(INOUT) :: CData(MOrb,MOrb,MOrb,MOrb,2) 00253 00254 REAL*8, INTENT(INOUT) :: RData(10),SILIntPeriod1,IntegratorTolError 00255 INTEGER*4, INTENT(INOUT) :: IData(30) 00256 LOGICAL, INTENT(INOUT) :: LData 00257 LOGICAL, INTENT(IN) :: restartABM 00258 00259 REAL*8, INTENT(INOUT) :: AbsTime 00260 INTEGER, INTENT(IN) :: NDX,NDY,NDZ,MOrb,ABMIntOrder 00261 00262 00263 COMPLEX*16 :: Psi_WORK(NDX*NDY*NDZ*MOrb), AuxPsi(NDX*NDY*NDZ*MOrb,ABMIntOrder+1) 00264 COMPLEX*16 :: DtPsi_WORK(NDX*NDY*NDZ*MOrb) 00265 INTEGER :: Steps,totalNoGridPts,RepeatedSteps,ErrorCode,k 00266 REAL*8 :: InitStep,ABMIntPeriod,norm 00267 00268 EXTERNAL AbsABMError 00269 00270 totalNoGridPts=NDX*NDY*NDZ 00271 InitStep=SILIntPeriod1 00272 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00273 !+++++++++++++++++++++++++++++++++++ STEP 4 +++++++++++++++++++++++++++++++++++++++++ 00274 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00275 ! propagate Psi_zero->Psi_tildetauhalf using rho_ij_zero and rho_ijkl_zero 00276 ! write(*,*)"++++++++ STEP 4 +++++++++++" 00277 ! change CData accordingly: 00278 00279 CData(:,:,1,1,1) = rho_jk_zero 00280 CData(:,:,:,:,2) = rho_ijkl_zero 00281 00282 ! the ABM integrator requires the time derivative of Psi_zero at the beginning 00283 ! FuncOrbs normalizes Psi_zero and returns DtPsi 00284 Psi_tildetauhalf = Psi_zero 00285 00286 00287 ! the ABM integrator takes and returns one-dimensional arrays only 00288 ! therefore it is necessary to copy the 2d array into a 1d array 00289 Psi_WORK = RESHAPE(Psi_tildetauhalf,(/NDX*NDY*NDZ*MOrb/)) 00290 00291 CALL FuncOrbs(AbsTime,Psi_WORK,DtPsi_WORK,CData,RData,IData,LData) ! get dPsi/dt 00292 ! write(*,*) SILIntPeriod1 00293 CALL ABM(Psi_WORK,DtPsi_WORK,totalNoGridPts*MOrb,SILIntPeriod1,AbsTime,ABMIntOrder,& 00294 InitStep,IntegratorTolError,restartABM,Steps,RepeatedSteps,& 00295 ErrorCode,AuxPsi,FuncOrbs,AbsABMError,CData, & 00296 RData,IData,LData) 00297 if(ErrorCode.ne.0) then 00298 !---{{{ error output 00299 write(*,*)"++++++++ STEP 4 +++++++++++" 00300 write(*,*)" ERROR OUTPUT" 00301 write(6,*)" totalNoGridPts = ",totalNoGridPts 00302 write(6,*)" SILIntPeriod1 = ",SILIntPeriod1 00303 write(6,*)" AbsTime = ",AbsTime 00304 write(6,*)" ABMIntOrder = ",ABMIntOrder 00305 write(6,*)" TolError = ",2.d0*IntegratorTolError 00306 write(6,*)" restartABM = ",restartABM 00307 write(6,*)" ABM Steps = ",Steps 00308 write(6,*)" ABM RepeatedSteps = ",RepeatedSteps 00309 write(6,*)" ABM ErrorCode = ",ErrorCode 00310 read(*,*) 00311 !----}}} 00312 end if 00313 Psi_tildetauhalf = RESHAPE(Psi_WORK,(/NDX*NDY*NDZ,MOrb/)) ! Psi_WORK -> Psi_tildetauhalf 00314 00315 ! normalize the orbitals 00316 00317 do k=1,MOrb 00318 00319 CALL normvxz(Psi_tildetauhalf(:,k),norm,NDX*NDY*NDZ) 00320 Psi_tildetauhalf(:,k) = Psi_tildetauhalf(:,k) / norm 00321 !---{{{error output 00322 if (ABS(norm-1.d0)>2.d0*IntegratorTolError) then 00323 write(*,*)"++++++++ STEP 4 +++++++++++" 00324 write(*,*)k,"after ABM in STEP 4 norm of Psi_tildetauhalf(:,k)",norm 00325 end if 00326 !---}}} 00327 00328 end do 00329 00330 END SUBROUTINE do_stepFour 00338 SUBROUTINE do_stepFive(CData,RData,IData,LData,rho_jk_tauhalf,rho_ijkl_tauhalf,Psi_tau,Psi_tauhalf,& 00339 readPotential,V_ext,NDX,NDY,NDZ,MOrb,nConf,nPar,Dims,gridx,gridy,gridz,AbsTime,SILIntPeriod1,& 00340 ABMIntOrder,IntegratorTolError,restartABM,nSuperDiags,BandMat,lambda0) 00341 00342 USE modulecomputationcore 00343 USE moduleextensions 00344 00345 LOGICAL,INTENT(IN) :: ReadPotential 00346 LOGICAL,INTENT(IN) :: restartABM 00347 00348 COMPLEX*16,INTENT(INOUT) :: Psi_tauhalf(NDX*NDY*NDZ,MOrb) 00349 COMPLEX*16,INTENT(INOUT) :: Psi_tau(NDX*NDY*NDZ,MOrb) 00350 COMPLEX*16,INTENT(INOUT) :: rho_ijkl_tauhalf(MOrb,MOrb,MOrb,MOrb),rho_jk_tauhalf(MOrb,MOrb) 00351 COMPLEX*16,INTENT(INOUT) :: CData(MOrb,MOrb,MOrb,MOrb,2) 00352 00353 COMPLEX*16,INTENT(INOUT) :: BandMat(nSuperDiags+1,nConf) 00354 REAL*8, INTENT(INOUT) :: RData(10),SILIntPeriod1,IntegratorTolError 00355 INTEGER*4, INTENT(INOUT) :: IData(30) 00356 LOGICAL, INTENT(INOUT) :: LData 00357 00358 REAL*8, INTENT(INOUT) :: AbsTime 00359 INTEGER, INTENT(IN) :: NDX,NDY,NDZ,MOrb,nPar,nConf,nSuperDiags,Dims 00360 INTEGER, INTENT(IN) :: ABMIntOrder 00361 REAL*8, INTENT(INOUT) :: gridx(NDX),gridy(NDY),gridz(NDZ),V_ext(NDX*NDY*NDZ),lambda0 00362 00363 COMPLEX*16 :: Psi_WORK(NDX*NDY*NDZ*MOrb), AuxPsi(NDX*NDY*NDZ*MOrb,ABMIntOrder+1) 00364 COMPLEX*16 :: DtPsi_WORK(NDX*NDY*NDZ*MOrb) 00365 INTEGER :: Steps,totalNoGridPts,RepeatedSteps,ErrorCode,k 00366 REAL*8 :: InitStep,ABMIntPeriod,norm 00367 00368 00369 00370 EXTERNAL AbsABMError 00371 00372 totalNoGridPts=NDX*NDY*NDZ 00373 InitStep=SILIntPeriod1 00374 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00375 !+++++++++++++++++++++++++++++++++++ STEP 5 +++++++++++++++++++++++++++++++++++++++++ 00376 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00377 ! Propagate Psi(tau/2) -> Psi(tau) using rho_jk(tau/2),rho_ijkl(tau/2) with the ABM integrator 00378 ! write(*,*)"++++++++ STEP 5 +++++++++++" 00379 00380 ! first adjust CData accordingly 00381 CData(:,:,1,1,1) = rho_jk_tauhalf 00382 CData(:,:,:,:,2) = rho_ijkl_tauhalf 00383 00384 ! the ABM integrator requires the time derivative of Psi_zero at the beginning 00385 ! FuncOrbs normalizes Psi_zero and returns DtPsi 00386 Psi_tau = Psi_tauhalf 00387 00388 ! the ABM integrator takes and returns one-dimensional arrays only 00389 ! therefore it is necessary to copy the 2d array into a 1d array 00390 Psi_WORK = RESHAPE(Psi_tau,(/NDX*NDY*NDZ*MOrb/)) 00391 00392 ! If a time-dependent external potential is used, it needs to be updated here 00393 IF (.NOT. readPotential) THEN 00394 CALL get_V_ext(V_ext,NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime+SILIntPeriod1) 00395 END IF 00396 CALL FuncOrbs(AbsTime+SILIntPeriod1,Psi_WORK,DtPsi_WORK,CData,RData,IData,LData) ! get dPsi/dt 00397 00398 CALL ABM(Psi_WORK,DtPsi_WORK,totalNoGridPts*MOrb,SILIntPeriod1,AbsTime+SILIntPeriod1,& 00399 ABMIntOrder,InitStep,IntegratorTolError,restartABM, & 00400 Steps,RepeatedSteps,ErrorCode,AuxPsi, & 00401 FuncOrbs,AbsABMError,CData,RData,IData,LData) 00402 00403 if(ErrorCode.ne.0) then 00404 !---{{{ error output 00405 write(*,*)"++++++++ STEP 5 +++++++++++" 00406 write(*,*)" ERROR OUTPUT" 00407 write(6,*)" totalNoGridPts = ",totalNoGridPts 00408 write(6,*)" SILIntPeriod1 = ",SILIntPeriod1 00409 write(6,*)" AbsTime = ",AbsTime 00410 write(6,*)" ABMIntOrder = ",ABMIntOrder 00411 write(6,*)" TolError = ",2.d0*IntegratorTolError 00412 write(6,*)" restartABM = ",restartABM 00413 write(6,*)" ABM Steps = ",Steps 00414 write(6,*)" ABM RepeatedSteps = ",RepeatedSteps 00415 write(6,*)" ABM ErrorCode = ",ErrorCode 00416 read(*,*) 00417 !----}}} 00418 end if 00419 00420 Psi_tau = RESHAPE(Psi_WORK,(/NDX*NDY*NDZ,MOrb/)) ! Psi_WORK -> Psi_tau 00421 00422 ! normalize the orbitals 00423 do k=1,MOrb 00424 CALL normvxz(Psi_tau(:,k),norm,NDX*NDY*NDZ) 00425 Psi_tau(:,k) = Psi_tau(:,k) / norm 00426 !---{{{error output 00427 if (ABS(norm-1.d0)>2.d0*IntegratorTolError) then 00428 write(*,*)"++++++++ STEP 5 +++++++++++" 00429 write(*,*)k,"after ABM in STEP 5 norm of Psi_tau(:,k)",norm 00430 end if 00431 !---}}} 00432 end do 00433 00434 ! If a time-dependent external potential is used, it needs to be updated here 00435 IF (.NOT. readPotential) THEN 00436 CALL get_V_ext(V_ext,NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime+2*SILIntPeriod1) 00437 END IF 00438 CALL get_FullHamiltonianBandMatrix(NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,& 00439 nSuperDiags,lambda0,Psi_tau,BandMat) 00440 00441 00442 00443 END SUBROUTINE do_stepFive 00444 00451 SUBROUTINE do_stepSix(CVec_tau,CVec_tauhalf,SILIntPeriod1,BandMat,CData,RData,& 00452 IData,LData,nConf,SILmaxIntOrder,IntegratorTolError,Relax,RestartSIL,& 00453 StdForm,SILIntPeriod2,AbsTime,iFail,MOrb,nSuperDiags,tau) 00454 00455 00456 USE modulecomputationcore 00457 00458 00459 INTEGER, INTENT(IN) :: nSuperDiags,SILmaxIntOrder,MOrb,nConf 00460 INTEGER, INTENT(INOUT) :: iFail 00461 COMPLEX*16 :: Krylov(nConf,2:SILmaxIntOrder),DtCVec(nConf) 00462 COMPLEX*16,INTENT(INOUT) :: CVec_tauhalf(nConf) 00463 COMPLEX*16,INTENT(INOUT) :: CVec_tau(nConf) 00464 COMPLEX*16,INTENT(INOUT) :: CData(MOrb,MOrb,MOrb,MOrb,2) 00465 COMPLEX*16,INTENT(INOUT) :: BandMat(nSuperDiags+1,nConf) 00466 00467 REAL*8, INTENT(INOUT) :: AbsTime,SILIntPeriod1,SILIntPeriod2 00468 REAL*8, INTENT(INOUT) :: RData(10),tau 00469 INTEGER*4, INTENT(INOUT) :: IData(30) 00470 REAL*8, INTENT(IN) :: IntegratorTolError 00471 logical, INTENT(INOUT) :: LData 00472 logical, INTENT(IN) :: Relax 00473 00474 logical :: RestartSIL,StdForm 00475 INTEGER :: TrueOrder,ErrorCode,Steps 00476 REAL*8 :: EigenVector(0:SILmaxIntOrder,0:SILmaxIntOrder+2), 00477 Diagonal(0:SILmaxIntOrder),OffDiag(SILmaxIntOrder), 00478 EigenVal(0:SILmaxIntOrder),OffDg2(SILmaxIntOrder+1), 00479 norm,IntPeriod 00480 00481 !------------- INITIALIZE VARIABLES-------------------- 00482 EigenVector = 0.d0 00483 Diagonal = 0.d0 00484 OffDiag = 0.d0 00485 EigenVal = 0.d0 00486 OffDg2 = 0.d0 00487 ErrorCode=0 00488 TrueOrder=0 00489 Steps=0 00490 RestartSIL=.TRUE. 00491 StdForm = .TRUE. 00492 00493 00494 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00495 !+++++++++++++++++++++++++++++++++++ STEP 6 +++++++++++++++++++++++++++++++++++++++++ 00496 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00497 ! Propagate CVec(tau/2)->CVec(tau) using H(tau) 00498 ! write(*,*)"++++++++ STEP 6 +++++++++++" 00499 CVec_tau = CVec_tauhalf 00500 CALL FuncCVecBandMat(AbsTime+2*SILIntPeriod1,CVec_tau,DtCVec,BandMat,RData,IData,LData) 00501 00502 ! integrate t = tau/2 -> tau 00503 IntPeriod = SILIntPeriod1 00504 CALL SILStep(CVec_tau,DtCVec,nConf,IntPeriod,SILmaxIntOrder,IntegratorTolError, & 00505 Relax,RestartSIL,StdForm,Steps,Krylov,SILIntPeriod2, & 00506 TrueOrder,ErrorCode,AbsTime,FuncCVecBandMat, & 00507 EigenVector,EigenVal,Diagonal,OffDg2,OffDiag, & 00508 BandMat,RData,IData,LData) 00509 00510 if (ErrorCode.ne.0) then 00511 write(*,*)"++++++++ STEP 6 +++++++++++" 00512 write(*,*)"SILStep ErrorCode =",ErrorCode 00513 STOP 00514 end if 00515 00516 ! Note that SILIntPeriod2 can be less than IntPeriod! 00517 ! compare the integration periods of the SIL integrations 00518 if(Abs(SILIntPeriod1-SILIntPeriod2)>1.d-10) then 00519 00520 write(*,*)"The SIL integration Periods of Step 2 and 6 are different" 00521 write(*,*)"SILIntPeriod1",SILIntPeriod1 00522 write(*,*)"SILIntPeriod2",SILIntPeriod2 00523 write(*,*)"repeat the step with a smaller step size" 00524 iFail = 1 00525 read(*,*) 00526 00527 end if 00528 00529 ! make sure CVec_tau is normalized 00530 CALL normvxz(CVec_tau,norm,nConf) 00531 CVec_tau = CVec_tau / norm 00532 00533 00534 END SUBROUTINE do_stepSix 00535 00545 SUBROUTINE do_stepSeven(CVec_tildezero,CVec_tauhalf,AbsTime,SILIntPeriod1,CVec_tau,BandMat,& 00546 CData,RData,IData,LData,nConf,MOrb,propDirection,nSuperDiags,SILmaxIntOrder,& 00547 IntegratorTolError,Stepsize,iFail,Relax) 00548 00549 use modulecomputationcore 00550 00551 INTEGER, INTENT(IN) :: nSuperDiags,SILmaxIntOrder,MOrb,nConf 00552 COMPLEX*16 :: Krylov(nConf,2:SILmaxIntOrder),DtCVec(nConf) 00553 COMPLEX*16,INTENT(INOUT) :: CVec_tauhalf(nConf),CVec_tau(nConf) 00554 COMPLEX*16,INTENT(INOUT) :: CVec_tildezero(nConf) 00555 COMPLEX*16,INTENT(INOUT) :: CData(MOrb,MOrb,MOrb,MOrb,2) 00556 COMPLEX*16,INTENT(INOUT) :: BandMat(nSuperDiags+1,nConf) 00557 REAL*8, INTENT(INOUT) :: AbsTime,SILIntPeriod1 00558 REAL*8, INTENT(INOUT) :: RData(10),Stepsize 00559 INTEGER*4, INTENT(INOUT) :: IData(30) 00560 REAL*8, INTENT(IN) :: IntegratorTolError 00561 logical, INTENT(INOUT) :: LData,Relax 00562 INTEGER, INTENT(INOUT) :: iFail,propDirection 00563 00564 logical :: RestartSIL,StdForm 00565 INTEGER :: TrueOrder,ErrorCode,Steps 00566 REAL*8 :: EigenVector(0:SILmaxIntOrder,0:SILmaxIntOrder+2), 00567 Diagonal(0:SILmaxIntOrder),OffDiag(SILmaxIntOrder), 00568 EigenVal(0:SILmaxIntOrder),OffDg2(SILmaxIntOrder+1), 00569 norm,IntPeriod 00570 !------------- INITIALIZE VARIABLES-------------------- 00571 EigenVector = 0.d0 00572 Diagonal = 0.d0 00573 OffDiag = 0.d0 00574 EigenVal = 0.d0 00575 OffDg2 = 0.d0 00576 ErrorCode=0 00577 TrueOrder=0 00578 Steps=0 00579 RestartSIL=.TRUE. 00580 StdForm = .TRUE. 00581 00582 00583 00584 00585 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00586 !+++++++++++++++++++++++++++++++++++ STEP 7 +++++++++++++++++++++++++++++++++++++++++ 00587 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00588 ! Propagate backwards CVec(tau/2)->CVe_tilde(0) using H(tau) for error 00589 ! estimation 00590 ! write(*,*)"++++++++ STEP 7 +++++++++++" 00591 ! get DtCVec for SILStep 00592 CVec_tildezero = CVec_tauhalf 00593 00594 IData(1) = (-1)*propDirection !backward propagation 00595 CALL FuncCVecBandMat(AbsTime+2.d0*SILIntPeriod1,CVec_tau,DtCVec,BandMat,RData,IData,LData) 00596 00597 ! integrate t = tau/2 -> 0 00598 IntPeriod = SILIntPeriod1 00599 CALL SILStep(CVec_tildezero,DtCVec,nConf,IntPeriod, & 00600 SILmaxIntOrder,IntegratorTolError, & 00601 Relax,RestartSIL,StdForm,Steps,Krylov,Stepsize, & 00602 TrueOrder,ErrorCode,AbsTime,FuncCVecBandMat, & 00603 EigenVector,EigenVal,Diagonal,OffDg2,OffDiag, & 00604 BandMat,RData,IData,LData) 00605 !---{{{ error output 00606 if (ErrorCode.ne.0) then 00607 write(*,*)"++++++++ STEP 7 +++++++++++" 00608 write(*,*)"SILStep ErrorCode =",ErrorCode 00609 STOP 00610 end if 00611 00612 00613 if (Abs(Stepsize-SILIntPeriod1)>1.d-10) then 00614 write(*,*)"++++++++ STEP 7 +++++++++++" 00615 write(*,*)"backward integration had a different length than the forward integration in STEP 2" 00616 write(*,*)"repeat the step with a smaller step size" 00617 write(*,*)"SILIntPeriod1",SILIntPeriod1 00618 write(*,*)"Stepsize",Stepsize 00619 iFail = 2 00620 end if 00621 !---}}} 00622 00623 00624 END SUBROUTINE do_stepSeven 00625 00626 END MODULE moduleintegration
1.7.3
