|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 !vim:fdm=marker: 00002 !############################################################################################ 00003 !# # 00004 !# ___ __ __ ____ _____ ____ _ _ ____ # 00005 !# / _ ` _ __ ___ _ __ | `/ |/ ___|_ _| _ `| | | | __ ) # 00006 !# | | | | '_ ` / _ ` '_ `| |`/| | | | | | | | | |_| | _ ` # 00007 !# | |_| | |_) | __/ | | | | | | |___ | | | |_| | _ | |_) | # 00008 !# `___/| .__/ `___|_| |_|_| |_|`____| |_| |____/|_| |_|____/ # 00009 !# |_| # 00010 !# Version 2.3 # 00011 !# # 00012 !#Kaspar Sakmann (head of project), Axel Lode, Alexej Streltsov, Ofir Alon, Lorenz Cederbaum# 00013 !# # 00014 !############################################################################################ 00120 00121 PROGRAM OpenMCTDHB 00122 00123 00124 USE moduleparameters 00125 USE moduleinputvariables 00126 USE moduledefaults 00127 USE modulederivesystemvariables 00128 USE moduleallocatables 00129 USE moduleinputoutput 00130 USE modulecomputationcore 00131 USE modulecorrelationfunctions 00132 USE modulegrid 00133 USE moduleextensions 00134 USE moduleintegration 00135 00136 IMPLICIT NONE 00137 external AbsABMError 00138 00139 ! READ INPUT PARAMETERS --------------------------------- 00140 ! the file runpars.inp will be read in. the file runpars.out 00141 ! will be (over)written after every printStep time units 00142 ! and can then be used for restarts 00143 00144 open(110,FILE='runpars.inp',DELIM='APOSTROPHE') 00145 read(110,NML=SystemParameters) 00146 read(110,NML=IMESTParameters) 00147 read(110,NML=GridParameters) 00148 read(110,NML=PrintParameters) 00149 read(110,NML=RunParameters) 00150 close(110) 00151 00152 ! CONSISTENCY CHECKS OF THE INPUT PARAMETERS ----------------- 00153 00154 CALL check_grid(Dims,NDX,NDY,NDZ) ! CHECK THE GRID AND DIMENSIONS 00155 00156 ! INITIALIZE VARIABLES THAT DEPEND ON THE INPUT VARIABLES 00157 00158 tag = TRIM(ADJUSTL(tag)) 00159 CALL get_nConf(NPar,MOrb,nConf) 00160 CALL get_nSuperDiags(NPar,MOrb,nSuperDiags) ! for M=1 there are no superdiagonals 00161 CALL get_totalNoGridPts(NDX,NDY,NDZ,Dims,totalNoGridPts) ! compute the number of grid points 00162 CALL get_identifier(MOrb,NPar,tag,identifier) ! determine an identifier string for the output files. 00163 CALL get_weight(xi,xf,yi,yf,zi,zf,NDX,NDY,NDZ,Dims,weight) 00164 00165 ! the input variables tinyStep, printStep and printAllStep control the intervalls 00166 ! output files are written: 00167 ! tinyStep controls scalar output variables that can be monitored continuously, 00168 ! e.g. the natural orbial occupation numbers. 00169 ! printNext controls when the one-body density/momentum distribution and a 00170 ! restart file is written. 00171 ! printAllNext controls when the first and second order correlation functions 00172 ! are written. 00173 printNextTiny = FLOOR(AbsTime/tinyStep +1.2d0) * tinyStep 00174 printNext = FLOOR(AbsTime/printStep + 1.2d0) * printStep 00175 printAllNext = FLOOR(AbsTime/printAllStep + 1.2d0) * printAllStep 00176 00177 00178 ! PRINT INFO ABOUT THIS RUN TO THE SCREEN 00179 00180 CALL print_header 00181 CALL print_runInfo(lambda0,NPar,MOrb,Dims,NDX,NDY,NDZ,xi,xf,yi,yf,zi,zf,tinyStep,printStep,printAllStep, & 00182 printCorrelationFuncs,AbsTime,maxTime,TolError,restart,Relax,readPotential,propDirection,tag) 00183 00184 00185 ! ALLOCATE AND (INITIALIZE TO ZERO) ALL SYSTEM SIZE DEPENDENT ARRAYS -------------------- 00186 ! ADD MORE ALLOCATABLE ARRAYS TO MODULEALLOCATABLES AND MODULEALLOCATE 00187 00188 CALL allocate_systemArrays(NDX,NDY,NDZ,MOrb,nConf,nSuperDiags) 00189 00190 CALL get_p2over2m(NDX,NDY,NDZ,xi,xf,yi,yf,zi,zf,Dims,p2over2m) !--the array p^2/2m for all p values 00191 CALL get_positionGrid(xi,xf,NDX,gridx) ! the x-grid 00192 CALL get_positionGrid(yi,yf,NDY,gridy) ! the y-grid 00193 CALL get_positionGrid(zi,zf,NDZ,gridz) ! the z-grid 00194 00195 ! compute the fourier transform of the interaction potential * (2pi)^0.5 / dV in FFT order 00196 IF (useIMEST) THEN 00197 CALL get_vtilde(NDX,NDY,NDZ,Dims,xi,xf,yi,yf,zi,zf,vtilde) 00198 END IF 00199 00200 00201 IData(1) = propDirection !propagation direction 1 forward, -1 backward 00202 RData(1) = lambda0 00203 RData(2) = TolError 00204 LData = Relax 00205 iFail = 0 00206 00207 !--------------------------------------------------------------------------- 00208 ! Read initial CVec_zero and Psi, or start from a default 00209 IF (restart) THEN 00210 00211 CALL read_CVecFile(CVec_zero,nConf) 00212 CALL read_PsiFromFile(Psi_zero,NDX,NDY,NDZ,MOrb) 00213 00214 ELSE 00215 00216 ! This usually works. All particles start in the first orbital 00217 CVec_zero = ZERO 00218 CVec_zero(1) = (1.d0,0.d0) 00219 ! Please set the default initial orbitals in get_defaultOrbitals yourself 00220 CALL get_defaultOrbitals(NDX,NDY,NDZ,Dims,gridx,gridy,gridz,MOrb,Psi_zero) 00221 00222 END IF 00223 00224 DO k=1,MOrb 00225 CALL normvxz(Psi_zero(:,k),norm,NDX*NDY*NDZ) 00226 Psi_zero(:,k) = Psi_zero(:,k) / norm 00227 END DO 00228 00229 CALL write_PsiToFile(AbsTime,Psi_zero,NDX,NDY,NDZ,MOrb,0,identifier) 00230 CALL write_PsiToFile(AbsTime,Psi_zero,NDX,NDY,NDZ,MOrb,1,identifier) 00231 00232 !-------------------------------------------------------------------------------- 00233 ! This controls what external potential is used. There are two options: the 00234 ! default is to read in a potential from the file ./restart/V_ext.dat. This 00235 ! allows to compile the code once and work with different potentials. 00236 ! However, if a time-dependent external potential is to be used, it is 00237 ! possible to specify it in the subroutine get_V_ext. The code will have 00238 ! to be recompiled in this case. 00239 IF (readPotential) THEN 00240 00241 CALL read_potentialFromFile(V_ext,NDX,NDY,NDZ) 00242 00243 ELSE 00244 00245 WRITE(*,*) "USE INTERNALLY DEFINED POTENTIAL get_V_ext" 00246 CALL get_V_ext(V_ext,NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime) 00247 00248 END IF 00249 00250 !-------------------------------------------------------------------------------- 00251 ! Before the real computation starts, we would like to write out data about 00252 ! the initial state. To do that a few quantities need to be computed. 00253 00254 ! Compute the Hamiltonian matrix. 00255 CALL get_FullHamiltonianBandMatrix(NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,& 00256 nSuperDiags,lambda0,Psi_zero,BandMat) 00257 00258 ! Compute the reduced first and second order density matrix elements 00259 ! and save them in CData, they will be needed in the computation. 00260 CALL get_rho_jk(nConf,MOrb,NPar,CVec_zero,rho_jk_zero) 00261 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_zero,rho_ijkl_zero) 00262 CData(:,:,1,1,1) = rho_jk_zero 00263 CData(:,:,:,:,2) = rho_ijkl_zero 00264 00265 ! compute the chemical potential of the initial guess 00266 CALL get_mu_kj(AbsTime,mu_kj,Psi_zero,rho_jk_zero,rho_ijkl_zero) 00267 CALL write_mu_kj(AbsTime,mu_kj,identifier) 00268 00269 ! Compute the natural orbitals and their occupations, also initialize 00270 ! the file that lists the times when a certain fragmentation is 00271 ! reached for the first time. 00272 CALL get_naturalOrbitals(MOrb,NDX,NDY,NDZ,rho_jk_zero,Psi_zero,noccs,norbs) 00273 CALL write_noccsToFile(AbsTime,MOrb,noccs,identifier) 00274 CALL write_norbsToFile(AbsTime,norbs,NDX,NDY,NDZ,MOrb,identifier) 00275 CALL write_firstFragmentationTimes(AbsTime,MOrb,NPar,noccs,lambda0,identifier) 00276 00277 00278 ! Compute the cumulative density from -infty up to x=0.d0, coded for 1D only so far 00279 ! and write it to a file. Also the mean density and the Lieb Liniger 00280 ! parameter are computed 00281 IF (Dims == 1) THEN 00282 CALL get_cumulativeProbability(Psi_zero,rho_jk_zero,0.d0,cumulativeProbability) 00283 CALL write_probabilityLeft(AbsTime,cumulativeProbability,MOrb,NPar,identifier) 00284 00285 CALL get_meanDensity(Psi_zero,rho_jk_zero,meanDensity) 00286 CALL write_liebLinigerParameter(AbsTime,lambda0,NPar,meanDensity,identifier) 00287 END IF 00288 00289 ! Compute the initial one-body density and write it to a file 00290 ! todo:2d 3d 00291 CALL get_oneBodyRDMDiagonal(AbsTime,Psi_zero,rho_jk_zero,1,1) 00292 00294 CALL write_correlations(Dims,printCorrelationFuncs,MOMSPACE2D,REALSPACE2D,& 00295 x1const,x1slice,y1const,y1slice,& 00296 x2const,x2slice,y2const,y2slice,& 00297 AbsTime,Psi_zero,rho_jk_zero,rho_ijkl_zero,& 00298 Pnot,xstart,xstop) 00299 00300 00301 ! Write the initial guess to a file 00302 CALL write_PsiToFile(AbsTime,Psi_zero,NDX,NDY,NDZ,MOrb,0,identifier) 00303 CALL write_PsiToFile(AbsTime,Psi_zero,NDX,NDY,NDZ,MOrb,1,identifier) 00304 CALL write_CVecFile(AbsTime,CVec_zero,0,identifier) 00305 CALL write_CVecFile(AbsTime,CVec_zero,1,identifier) 00306 00307 ! Compute the energy of the initial guess and write it to a file 00308 CALL get_energy(CVec_zero,BandMat,energy,nConf) 00309 CALL write_energy(AbsTime,lambda0,NPar,energy,identifier) 00310 00311 ! integrate the density over some coordinates (sometimes useful for visualization in 2d/3d) 00312 ! remove this if you don't want it. 00313 CALL get_integratedDensity(rho_jk_zero,AbsTime,Psi_zero,writemode=3) !yz integrated out 00314 CALL get_integratedDensity(rho_jk_zero,AbsTime,Psi_zero,writemode=4) !x integrated out 00315 CALL get_integratedDensity(rho_jk_zero,AbsTime,Psi_zero,writemode=1) !z integrated out 00316 CALL get_integratedDensity(rho_jk_zero,AbsTime,Psi_zero,writemode=2) !y integrated out 00317 00318 !-------------------------------------------------------------------------------- 00319 00320 tau = 0.01d0 00321 IntPeriod = tau 00322 InitStep = tau 00323 Stepsize = tau ! output of silstep: integration step performed, can be < IntPeriod 00324 Steps = 0 00325 RepeatedSteps = 0 00326 IntegratorTolError = 5.0D-1 * TolError 00327 00328 00329 do k=1,MOrb 00330 CALL normvxz(Psi_zero(:,k),norm,NDX*NDY*NDZ) 00331 Psi_zero(:,k) = Psi_zero(:,k) / norm 00332 !---{{{error output 00333 if (ABS(norm-1.d0)>TolError) then 00334 write(*,*)k,"norm of initial Psi_zero(:,k) ",norm 00335 end if 00336 !---}}} 00337 end do 00338 00339 Psi_tau = Psi_zero ! initialize Psi_tau for the propagation loop 00340 CALL normvxz(CVec_zero,norm,nConf) 00341 CVec_zero = CVec_zero / norm 00342 !---{{{error output 00343 if (ABS(norm-1.d0)>TolError) then 00344 write(*,*)"norm of initial CVec ",norm 00345 end if 00346 !---}}} 00347 00348 ! initialize CVec_tau for the propagation loop 00349 CVec_tau = CVec_zero 00350 00351 !----------------------------- PROPAGATION LOOP ------------------------------------ 00352 propagation_loop: do 00353 00354 CVec_zero = CVec_tau 00355 Psi_zero = Psi_tau 00356 00357 tau1 = tau 00358 tau2 = tau 00359 tau3 = tau 00360 tausave = tau 00361 iprintNextTiny = 0 00362 iprintNext = 0 00363 iprintAllNext = 0 00364 00365 if (AbsTime+tau > printNextTiny) then 00366 iprintNextTiny = 1 00367 tau1 = printNextTiny-AbsTime 00368 end if 00369 00370 if (AbsTime+tau > printNext) then 00371 iprintNextTiny = 1 00372 iprintNext = 1 00373 tau2 = printNext-AbsTime 00374 end if 00375 00376 if (AbsTime+tau > printAllNext) then 00377 iprintNextTiny = 1 00378 iprintNext = 1 00379 iprintAllNext = 1 00380 tau3 = printAllNext-AbsTime 00381 end if 00382 00383 tau = MIN(ABS(tausave),ABS(tau1),ABS(tau2),ABS(tau3)) 00384 if (ABS(tau) <= 1E-6) tau = 1E-6 00385 00386 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00387 !+++++++++++++++++++++++++++++++++++ STEP 1 +++++++++++++++++++++++++++++++++++++++++ 00388 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00389 CALL do_stepOne(nConf,MOrb,NPar,CVec_zero,CData,readPotential,V_ext,& 00390 NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime,nSuperDiags,lambda0,& 00391 Psi_zero,BandMat,rho_jk_zero,rho_ijkl_zero) 00392 00393 00394 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00395 !+++++++++++++++++++++++++++++++++++ STEP 2 +++++++++++++++++++++++++++++++++++++++++ 00396 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00397 00398 CALL do_stepTwo(AbsTime,CVec_zero,DtCVec,BandMat,RData,IData,LData,& 00399 CVec_tauhalf,NPar,nConf,MOrb,IntPeriod,SILmaxIntOrder,& 00400 IntegratorTolError,Relax,CData,tau, & 00401 SILIntPeriod1,nSuperDiags,& 00402 rho_ijkl_tauhalf,rho_jk_tauhalf,propDirection) 00403 00404 00405 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00406 !+++++++++++++++++++++++++++++++++++ STEP 3 +++++++++++++++++++++++++++++++++++++++++ 00407 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00408 00409 00410 CALL do_stepThree(Psi_tauhalf,Psi_zero,Psi_WORK,AbsTime,DtPsi_WORK,& 00411 CData,RData,IData,LData,NDX,NDY,NDZ,MOrb,Relax,tau,& 00412 IntegratorTolError,restartABM,SILIntPeriod1,& 00413 ABMIntOrder) 00414 00415 00416 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00417 00418 propagateif: if (.NOT. Relax) then ! the following Steps are not needed for Relaxation 00419 00420 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00421 !+++++++++++++++++++++++++++++++++++ STEP 4 +++++++++++++++++++++++++++++++++++++++++ 00422 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00423 CALL do_stepFour(rho_jk_zero,rho_ijkl_zero,CData,Rdata,IData,LData,& 00424 Psi_tildetauhalf,Psi_zero,restartABM,SILIntPeriod1,& 00425 AbsTime,IntegratorTolError,NDX,NDY,NDZ,MOrb,ABMIntOrder) 00426 00427 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00428 !+++++++++++++++++++++++++++++++++++ STEP 5 +++++++++++++++++++++++++++++++++++++++++ 00429 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00430 CALL do_stepFive(CData,RData,IData,LData,rho_jk_tauhalf,rho_ijkl_tauhalf,Psi_tau,Psi_tauhalf,& 00431 readPotential,V_ext,NDX,NDY,NDZ,MOrb,nConf,nPar,Dims,gridx,gridy,gridz,AbsTime,SILIntPeriod1,& 00432 ABMIntOrder,IntegratorTolError,restartABM,nSuperDiags,BandMat,lambda0) 00433 00434 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00435 !+++++++++++++++++++++++++++++++++++ STEP 6 +++++++++++++++++++++++++++++++++++++++++ 00436 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00437 CALL do_stepSix(CVec_tau,CVec_tauhalf,SILIntPeriod1,BandMat,CData,RData,& 00438 IData,LData,nConf,SILmaxIntOrder,IntegratorTolError,Relax,RestartSIL,& 00439 StdForm,SILIntPeriod2,AbsTime,iFail,MOrb,nSuperDiags,tau) 00440 00441 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00442 !+++++++++++++++++++++++++++++++++++ STEP 7 +++++++++++++++++++++++++++++++++++++++++ 00443 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00444 CALL do_stepSeven(CVec_tildezero,CVec_tauhalf,AbsTime,SILIntPeriod1,CVec_tau,BandMat,& 00445 CData,RData,IData,LData,nConf,MOrb,propDirection,nSuperDiags,SILmaxIntOrder,& 00446 IntegratorTolError,Stepsize,iFail,Relax) 00447 00448 00449 00450 end if propagateif ! Relaxation will continue hereafter 00451 00452 ! make sure CVec_tildezero is normalized 00453 CALL normvxz(CVec_tildezero,norm,nConf) 00454 CVec_tildezero = CVec_tildezero / norm 00455 00456 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00457 00458 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00459 !+++++++++++++++++++++++++++++++++++ ERROR EVALUATION +++++++++++++++++++++++++++++++ 00460 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 00461 00462 if (.NOT. Relax) then 00463 00464 ! The error of the propagation is composed of two parts: one from the 00465 ! propagation of the coefficient vector CVec and one from the propagation of the 00466 ! orbitals 00467 Delta_CVec = CVec_zero - CVec_tildezero 00468 Delta_Psi = Psi_tauhalf-Psi_tildetauhalf 00469 00470 else 00471 00472 CVec_tau = CVec_tauhalf 00473 Psi_tau = Psi_tauhalf 00474 rho_jk_tau = rho_jk_tauhalf 00475 rho_ijkl_tau = rho_ijkl_tauhalf 00476 00477 IF (.NOT. readPotential) THEN 00478 CALL get_V_ext(V_ext,NDX,NDY,NDZ,Dims,gridx,gridy,gridz,AbsTime+SILIntPeriod1) 00479 END IF 00480 00481 CALL get_FullHamiltonianBandMatrix(NPar,MOrb,NDX,NDY,NDZ,Dims,nConf,& 00482 nSuperDiags,lambda0,Psi_tau,BandMat) 00483 00484 Delta_CVec = CVec_tau - CVec_zero 00485 Delta_Psi = Psi_tau - Psi_zero 00486 00487 end if 00488 00489 CALL evaluate_error(Delta_CVec,Delta_Psi,rho_jk_zero,CVecError,PsiError,nConf) 00490 00491 delta = 2*MAX(CVecError,PsiError) 00492 if (delta > TolError) iFail = 3 00493 00494 REJECTSTEPIF: if (iFail.ne.0) then 00495 00496 CVec_tau = CVec_zero 00497 Psi_tau = Psi_zero 00498 tau = 0.5d0 * tau * (TolError/delta)**0.25 00499 iFail = 0 00500 00501 else 00502 00503 AbsTime = AbsTime + tau 00504 00505 ! normalize CVec_tau if the step is accepted 00506 CALL normvxz(CVec_tau,norm,nConf) 00507 CVec_tau = CVec_tau / norm 00508 00509 ! normalize Psi_tau if the step is accepted 00510 do k=1,MOrb 00511 CALL normvxz(Psi_tau(:,k),norm,NDX*NDY*NDZ) 00512 Psi_tau(:,k) = Psi_tau(:,k) / norm 00513 end do 00514 00515 ! determine new step size 00516 tau = tau * MIN((TolError/delta)**0.25,1.5d0) 00517 00518 00519 if (iprintNextTiny==1) then 00520 00521 CALL get_energy(CVec_tau,BandMat,energy,nConf) 00522 CALL write_energy(AbsTime,lambda0,NPar,energy,identifier) 00523 00524 if (Relax) then 00525 CALL get_defect(nConf,nSuperDiags,BandMat,CVec_tau,energy,defect) 00526 WRITE(*,1111)"AbsTime:",AbsTime," tau:",tau," PsiError:",PsiError,& 00527 " CVecError:",CVecError," defect:",defect 00528 1111 FORMAT(A,F10.4,A,F10.5,A,ES14.6,A,ES14.6,A,ES14.6) 00529 end if 00530 00531 CALL get_rho_jk(nConf,MOrb,NPar,CVec_tau,rho_jk_tau) 00532 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_tau,rho_ijkl_tau) 00533 00535 IF (Dims==1) THEN 00536 CALL get_cumulativeProbability(Psi_tau,rho_jk_tau,0.d0,cumulativeProbability) 00537 CALL write_probabilityLeft(AbsTime,cumulativeProbability,MOrb,NPar,identifier) 00538 00539 CALL get_meanDensity(Psi_tau,rho_jk_tau,meanDensity) 00540 CALL write_liebLinigerParameter(AbsTime,lambda0,NPar,meanDensity,identifier) 00541 END IF 00542 00543 00544 CALL get_mu_kj(AbsTime,mu_kj,Psi_tau,rho_jk_tau,rho_ijkl_tau) 00545 CALL write_mu_kj(AbsTime,mu_kj,identifier) 00546 00547 CALL get_rho_jk(nConf,MOrb,NPar,CVec_tau,rho_jk_tau) 00548 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_tau,rho_ijkl_tau) 00549 CALL get_naturalOrbitals(MOrb,NDX,NDY,NDZ,rho_jk_tau,Psi_tau,noccs,norbs) 00550 CALL write_noccsToFile(AbsTime,MOrb,noccs,identifier) 00551 00552 CALL write_firstFragmentationTimes(AbsTime,MOrb,NPar,noccs,lambda0,identifier) 00553 00554 printNextTiny = FLOOR(AbsTime/tinyStep + 1.2d0) * tinyStep 00555 iprintNextTiny = 0 00556 00557 if (iprintNext==1) then 00558 write(*,*)"AbsTime = ",AbsTime 00559 open(111,FILE='runpars.out',DELIM='APOSTROPHE') 00560 write(UNIT=111,NML=SystemParameters) 00561 write(UNIT=111,NML=IMESTParameters) 00562 write(UNIT=111,NML=GridParameters) 00563 write(UNIT=111,NML=PrintParameters) 00564 write(UNIT=111,NML=RunParameters) 00565 close(111) 00566 00567 CALL write_PsiToFile(AbsTime,Psi_tau,NDX,NDY,NDZ,MOrb,1,identifier) 00568 CALL write_CVecFile(AbsTime,CVec_tau,1,identifier) 00569 CALL get_oneBodyRDMDiagonal(AbsTime,Psi_tau,rho_jk_tau,1,1) !-write 1B density 00570 00572 IF(Dims==1) THEN 00573 CALL pedestrian_FT(Psi_tau,FTPsi,1) 00574 CALL get_dilation(dilation,FTPsi,rho_jk_tau) 00575 CALL pedestrian_FT(Psi_tau,FTPsi,dilation) 00576 CALL get_oneBodyRDMDiagonal(AbsTime,FTPsi,rho_jk_tau,2,dilation) !-write 1B momentum distribution 00577 END IF 00578 00579 CALL get_rho_jk(nConf,MOrb,NPar,CVec_tau,rho_jk_tau) 00580 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_tau,rho_ijkl_tau) 00581 CALL get_naturalOrbitals(MOrb,NDX,NDY,NDZ,rho_jk_tau,Psi_tau,noccs,norbs) 00582 CALL write_noccsToFile(AbsTime,MOrb,noccs,identifier) 00583 00584 printNext = FLOOR(AbsTime/printStep + 1.2d0) * printStep 00585 iprintNext = 0 00586 00587 if (iprintAllNext==1) then 00588 00589 CALL write_PsiToFile(AbsTime,Psi_tau,NDX,NDY,NDZ,MOrb,0,identifier) 00590 CALL write_CVecFile(AbsTime,CVec_tau,0,identifier) 00591 CALL write_norbsToFile(AbsTime,norbs,NDX,NDY,NDZ,MOrb,identifier) 00592 00593 CALL write_correlations(Dims,printCorrelationFuncs,MOMSPACE2D,REALSPACE2D,& 00594 x1const,x1slice,y1const,y1slice,& 00595 x2const,x2slice,y2const,y2slice,& 00596 AbsTime,Psi_tau,rho_jk_tau,rho_ijkl_tau,& 00597 Pnot,xstart,xstop) 00598 00599 iprintAllNext = 0 00600 printAllNext = FLOOR(AbsTime/printAllStep + 1.2d0) * printAllStep 00601 00602 end if 00603 00604 end if 00605 00606 tau = MAX(tausave,tau) 00607 00608 end if 00609 00610 end if REJECTSTEPIF 00611 00612 if (AbsTime.ge.maxTime) then 00613 write(*,*)" maxtime reached" 00614 exit propagation_loop 00615 end if 00616 00617 00618 end do propagation_loop 00619 00620 CALL get_rho_jk(nConf,MOrb,NPar,CVec_tau,rho_jk_tau) 00621 CALL get_rho_ijkl(nConf,MOrb,NPar,CVec_tau,rho_ijkl_tau) 00622 CALL get_naturalOrbitals(MOrb,NDX,NDY,NDZ,rho_jk_tau,Psi_tau,noccs,norbs) 00623 CALL write_norbsToFile(AbsTime,norbs,NDX,NDY,NDZ,MOrb,identifier) 00624 CALL write_noccsToFile(AbsTime,MOrb,noccs,identifier) 00625 CALL write_firstFragmentationTimes(AbsTime,MOrb,NPar,noccs,lambda0,identifier) 00626 00627 CALL write_PsiToFile(AbsTime,Psi_tau,NDX,NDY,NDZ,MOrb,0,identifier) 00628 CALL write_PsiToFile(AbsTime,Psi_tau,NDX,NDY,NDZ,MOrb,1,identifier) 00629 CALL write_CVecFile(AbsTime,CVec_tau,0,identifier) 00630 CALL write_CVecFile(AbsTime,CVec_tau,1,identifier) 00631 00632 ! integrate the density over some coordinates (sometimes useful for visualization in 2d/3d) 00633 ! remove this if you don't want it. 00634 CALL get_integratedDensity(rho_jk_tau,AbsTime,Psi_tau,writemode=3) !yz integrated out 00635 CALL get_integratedDensity(rho_jk_tau,AbsTime,Psi_tau,writemode=4) !x integrated out 00636 CALL get_integratedDensity(rho_jk_tau,AbsTime,Psi_tau,writemode=1) !z integrated out 00637 CALL get_integratedDensity(rho_jk_tau,AbsTime,Psi_tau,writemode=2) !y integrated out 00638 00639 00641 IF (Dims == 1) THEN 00642 CALL get_oneBodyRDMDiagonal(AbsTime,Psi_tau,rho_jk_tau,1,1) !-write 1B density 00643 END IF 00644 00645 CALL deallocate_systemArrays 00646 00647 write(*,*)"END OF PROGRAM REACHED, ADIEU." 00648 00649 end program OpenMCTDHB 00650 00651
1.7.3
