|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 ! vim:fdm=marker: 00002 !----------------------------------------------------------------------------- 00008 MODULE moduleinputoutput 00009 IMPLICIT NONE 00010 00011 00012 CONTAINS 00013 00014 !------------------------------------------------------------------------------ 00019 SUBROUTINE print_header 00020 !---{{{ 00021 IMPLICIT NONE 00022 WRITE(*,*) 00023 WRITE(*,1100)" ##############################################################################################################" 00024 WRITE(*,1100)" # ___ __ __ ____ _____ ____ _ _ ____ #" 00025 WRITE(*,1100)" # / _ ` _ __ ___ _ __ | `/ |/ ___|_ _| _ `| | | | __ ) #" 00026 WRITE(*,1100)" # | | | | '_ ` / _ ` '_ `| |`/| | | | | | | | | |_| | _ ` #" 00027 WRITE(*,1100)" # | |_| | |_) | __/ | | | | | | |___ | | | |_| | _ | |_) | #" 00028 WRITE(*,1100)" # `___/| .__/ `___|_| |_|_| |_|`____| |_| |____/|_| |_|____/ #" 00029 WRITE(*,1100)" # |_| #" 00030 WRITE(*,1100)" # #" 00031 WRITE(*,1100)" # (2007-2012) Version 2.3 #" 00032 WRITE(*,1100)" # #" 00033 WRITE(*,1100)" # Kaspar Sakmann (head of project), Axel Lode, Alexej Streltsov, Ofir Alon, Lorenz Cederbaum #" 00034 WRITE(*,1100)" # #" 00035 WRITE(*,1100)" ##############################################################################################################" 00036 WRITE(*,*) 00037 00038 1100 FORMAT(116A) 00039 !---}}} 00040 END SUBROUTINE print_header 00041 00042 !------------------------------------------------------------------------------- 00048 SUBROUTINE print_runInfo(lambda0,NPar,MOrb,Dims,NDX,NDY,NDZ,xi,xf,yi,yf,zi,zf,tinyStep,printStep,printAllStep, & 00049 printCorrelationFuncs,AbsTime,maxTime,TolError,restart,Relax,readPotential,propDirection,tag) 00050 !---{{{ 00051 IMPLICIT NONE 00052 REAL*8,INTENT(IN) :: lambda0 00053 INTEGER ,INTENT(IN):: NPar 00054 INTEGER,INTENT(IN) :: MOrb 00055 INTEGER,INTENT(IN) :: Dims 00056 00057 INTEGER,INTENT(IN) :: NDX 00058 INTEGER,INTENT(IN) :: NDY 00059 INTEGER,INTENT(IN) :: NDZ 00060 REAL*8,INTENT(IN) :: xi 00061 REAL*8,INTENT(IN) :: xf 00062 REAL*8,INTENT(IN) :: yi 00063 REAL*8,INTENT(IN) :: yf 00064 REAL*8,INTENT(IN) :: zi 00065 REAL*8,INTENT(IN) :: zf 00066 00067 REAL*8,INTENT(IN) :: tinyStep 00068 REAL*8,INTENT(IN) :: printStep 00069 REAL*8,INTENT(IN) :: printAllStep 00070 INTEGER,INTENT(IN) :: printCorrelationFuncs 00071 00072 REAL*8,INTENT(IN) :: AbsTime 00073 REAL*8,INTENT(IN) :: maxTime 00074 REAL*8,INTENT(IN) :: TolError 00075 LOGICAL,INTENT(IN) :: restart 00076 LOGICAL,INTENT(IN) :: Relax 00077 LOGICAL,INTENT(IN) :: readPotential 00078 INTEGER,INTENT(IN) :: propDirection 00079 CHARACTER(55),INTENT(IN) :: tag 00080 00081 WRITE(*,*) 00082 IF (Relax) THEN 00083 00084 IF (propDirection == 1) THEN 00085 WRITE(*,109)"Relaxation run. Propagate in imaginary time to the ground state." 00086 ELSE IF (propDirection == -1) THEN 00087 WRITE(*,112)"Relax = ",Relax," propDirection",propDirection 00088 WRITE(*,109)"RELAXATION IN THE WRONG DIRECTION. THIS WILL BLOW UP. I WILL THEREFORE STOP." 00089 WRITE(*,*) 00090 STOP 00091 ELSE 00092 WRITE(*,110)"ILLEGAL VALUE FOR propagationDirection",propDirection 00093 STOP 00094 END IF 00095 00096 ELSE 00097 IF (propDirection == 1) THEN 00098 WRITE(*,109)"Propagation run. Propagate forward in time." 00099 ELSE IF (propDirection == -1) THEN 00100 WRITE(*,109)"Propagation run. Propagate backwards in time (count forwards)." 00101 ELSE 00102 WRITE(*,110)"ILLEGAL VALUE FOR propagationDirection: ",propDirection 00103 STOP 00104 END IF 00105 00106 END IF 00107 00108 IF (MOrb== 1) THEN 00109 WRITE(*,109)"Gross-Pitaevskii computation" 00110 ELSE IF (MOrb == 2) THEN 00111 WRITE(*,109)"Many-body computation" 00112 ELSE IF (MOrb > 2) THEN 00113 WRITE(*,109)"More than two orbitals will be available in the next release, but not yet." 00114 STOP 00115 ELSE 00116 WRITE(*,110)"ILLEGAL VALUE FOR MOrb",MOrb 00117 STOP 00118 END IF 00119 00120 IF (restart) THEN 00121 WRITE(*,109)"Restart from restart.orbs.inp and restart.coef.inp " 00122 ELSE IF (.NOT. restart) THEN 00123 WRITE(*,109)"Start from default orbitals: harmonic oscillator eigenfunctions" 00124 ELSE 00125 WRITE(*,*)"ILLEGAL VALUE FOR restart",restart 00126 STOP 00127 END IF 00128 00129 IF (readPotential) THEN 00130 WRITE(*,109)"external potential will be read from V_ext.inp " 00131 ELSE IF (.NOT. readPotential) THEN 00132 WRITE(*,109)"external potential defined in get_V_ext" 00133 ELSE 00134 WRITE(*,*)"ILLEGAL VALUE FOR readPotential",readPotential 00135 STOP 00136 END IF 00137 00138 00139 00140 00141 WRITE(*,*) 00142 00143 IF (Dims < 4 .AND. Dims >0) THEN 00144 WRITE(*,113)"Number of dimensions = ",Dims 00145 ELSE 00146 WRITE(*,113)"INVALID NUMBER OF DIMENSIONS ",Dims 00147 STOP 00148 END IF 00149 WRITE(*,113)"Number of bosons = ",NPar 00150 WRITE(*,113)"Number of orbitals = ",MOrb 00151 WRITE(*,111)"Interaction strength lambda0 = ",lambda0 00152 WRITE(*,111)"lambda0(N-1) = ",lambda0*(NPar-1) 00153 00154 WRITE(*,*) 00155 00156 WRITE(*,117)"x box: [",xi,",",xf,") grid points: ",NDX 00157 WRITE(*,117)"y box: [",yi,",",yf,") grid points: ",NDY 00158 WRITE(*,117)"z box: [",zi,",",zf,") grid points: ",NDZ 00159 WRITE(*,*) 00160 WRITE(*,114)"Initial Time = ",AbsTime 00161 WRITE(*,114)"Final Time = ",maxTime 00162 WRITE(*,*) 00163 WRITE(*,115)"Error tolerance = ",TolError 00164 WRITE(*,*) 00165 00166 109 FORMAT(5X,A) 00167 110 FORMAT(5X,A60,I6) 00168 111 FORMAT(5X,A45,F13.9) 00169 112 FORMAT(5X,A15,L,A15,I6) 00170 113 FORMAT(5X,A45,I6) 00171 114 FORMAT(5X,A45,F19.9) 00172 115 FORMAT(5X,A45,ES14.4) 00173 117 FORMAT(13X,A,F8.2,A,F8.2,A,I8) 00174 !---}}} 00175 END SUBROUTINE print_runInfo 00176 00177 !----------------------------------------------------------------------------- 00183 SUBROUTINE get_doubleAsString(x,doubleString) 00184 ! {{{ 00185 IMPLICIT NONE 00186 INTEGER,parameter :: DBL=8 00187 real(kind=DBL),INTENT(IN) :: x 00188 CHARACTER(len=10),INTENT(OUT) :: doubleString 00189 00190 IF (x <= 999999.999d0) WRITE(doubleString,'(F10.3)') x 00191 IF (x <= 99999.9999d0) WRITE(doubleString,'(F10.4)') x 00192 IF (x <= 9999.99999d0) WRITE(doubleString,'(F10.5)') x 00193 IF (x <= 999.999999d0) WRITE(doubleString,'(F10.6)') x 00194 IF (x <= 99.9999999d0) WRITE(doubleString,'(F10.7)') x 00195 IF (x <= 9.99999999d0) WRITE(doubleString,'(F10.8)') x 00196 00197 IF ((x < 0.d0) .OR. (x >= 1000000.0d0)) THEN 00198 WRITE(6,*)"x not implemented" 00199 STOP 00200 END IF 00201 00202 !}}} 00203 END SUBROUTINE get_doubleAsString 00204 00205 !---------------------------------------------------------------------------- 00211 SUBROUTINE get_identifier(MOrb,NPar,tag,identifier) 00212 !---{{{ 00213 IMPLICIT NONE 00214 INTEGER,INTENT(IN) :: MOrb,NPar 00215 CHARACTER(len=255),INTENT(OUT):: identifier 00216 CHARACTER(len=55), INTENT(IN) :: tag 00217 CHARACTER(len=10) :: MOrbAsString 00218 CHARACTER(len=1) :: NparAsString1 00219 CHARACTER(len=2) :: NparAsString2 00220 CHARACTER(len=3) :: NparAsString3 00221 CHARACTER(len=4) :: NparAsString4 00222 CHARACTER(len=5) :: NparAsString5 00223 CHARACTER(len=6) :: NparAsString6 00224 CHARACTER(len=7) :: NparAsString7 00225 00226 IF (MOrb>0 .AND. MOrb<10) THEN 00227 write (MOrbAsString, '(I1)') MOrb 00228 ELSE IF (MOrb>=10 .AND. MOrb<100) THEN 00229 write (MOrbAsString, '(I2)') MOrb 00230 ELSE 00231 write (*,*) "TOO MANY/FEW ORBITALS MOrb: ", MOrb 00232 END IF 00233 00234 00235 IF ((NPar>0) .and. (NPar<10)) THEN 00236 write (NParAsString1, '(I1)') NPar 00237 identifier='N'//trim(NParAsString1)//'M'//trim(MOrbAsString)//trim(tag) 00238 ELSE IF ((NPar>=10) .and. (NPar<100)) THEN 00239 write (NParAsString2, '(I2)') NPar 00240 identifier='N'//trim(NParAsString2)//'M'//trim(MOrbAsString)//trim(tag) 00241 ELSE IF ((NPar>=100) .and. (NPar<1000)) THEN 00242 write (NParAsString3, '(I3)') NPar 00243 identifier='N'//trim(NParAsString3)//'M'//trim(MOrbAsString)//trim(tag) 00244 ELSE IF ((NPar>=1000) .and. (NPar<10000)) THEN 00245 write (NParAsString4, '(I4)') NPar 00246 identifier='N'//trim(NParAsString4)//'M'//trim(MOrbAsString)//trim(tag) 00247 ELSE IF ((NPar>=10000) .and. (NPar<100000)) THEN 00248 write (NParAsString5, '(I5)') NPar 00249 identifier='N'//trim(NParAsString5)//'M'//trim(MOrbAsString)//trim(tag) 00250 ELSE IF ((NPar>=100000) .and. (NPar<1000000)) THEN 00251 write (NParAsString6, '(I6)') NPar 00252 identifier='N'//trim(NParAsString6)//'M'//trim(MOrbAsString)//trim(tag) 00253 ELSE IF ((NPar>=1000000) .and. (NPar<10000000)) THEN 00254 write (NParAsString7, '(I7)') NPar 00255 identifier='N'//trim(NParAsString7)//'M'//trim(MOrbAsString)//trim(tag) 00256 ELSE 00257 write(*,*)"Npar is too big in get_identifier: ",Npar 00258 STOP 00259 END IF 00260 !---}}} 00261 END SUBROUTINE get_identifier 00262 00263 !---------------------------------------------------------------------------------------- 00270 SUBROUTINE get_cumulativeProbability(Psi,rho_jk,x,cumulativeProbability) 00271 USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb,NPar 00272 USE moduleallocatables,ONLY:gridx 00273 ! {{{ 00274 00275 IMPLICIT NONE 00276 REAL*8,INTENT(OUT) :: cumulativeProbability 00277 REAL*8,INTENT(IN) :: x 00278 COMPLEX*16,DIMENSION(NDX*NDY*NDZ,MOrb),INTENT(IN) :: Psi 00279 COMPLEX*16,DIMENSION(MOrb,MOrb),INTENT(IN) :: rho_jk 00280 INTEGER :: n,Iorb,Jorb 00281 COMPLEX*16 :: G1n_n 00282 00283 IF((NDY.ne.1).or.(NDZ.ne.1)) THEN 00284 write(*,*)"ONLY 1d implemented in get_cumulativeProbability" 00285 STOP 00286 END IF 00287 00288 cumulativeProbability=0.d0 00289 00290 DO n=1,NDX*NDY*NDZ 00291 G1n_n = (0.d0,0.d0) 00292 DO Jorb = 1,MOrb 00293 DO Iorb=1,MOrb 00294 00295 G1n_n = G1n_n + rho_jk(Iorb,Jorb) * DCONJG(Psi(n,Iorb))*Psi(n,Jorb) 00296 00297 END DO 00298 END DO 00299 00300 IF (gridx(n)<x) THEN 00301 cumulativeProbability = cumulativeProbability + DBLE(G1n_n) 00302 END IF 00303 00304 END DO 00305 00306 cumulativeProbability = cumulativeProbability / Npar 00307 00308 !---}}} 00309 END SUBROUTINE get_cumulativeProbability 00310 00311 !---------------------------------------------------------------------------- 00318 SUBROUTINE write_probabilityLeft(AbsTime,cumulativeProbability,MOrb,NPar,identifier) 00319 !---{{{ 00320 IMPLICIT NONE 00321 REAL*8,INTENT(IN) :: AbsTime,cumulativeProbability 00322 INTEGER, INTENT(IN) :: MOrb,NPar 00323 CHARACTER(255),INTENT(IN) :: identifier 00324 CHARACTER(255) :: filename 00325 INTEGER,save :: counter=0 00326 00327 filename = 'probabilityLeft'//trim(identifier)//'.dat' 00328 00329 IF (counter == 0 ) THEN 00330 00331 open(unit=1979,file=filename,form='FORMATTED',status='REPLACE') 00332 write(1979,*)"# " 00333 write(1979,*)"# the integral over the probability density in the left half of space " 00334 write(1979,*)"# (accurate up to one grid point)" 00335 write(1979,*)"# " 00336 write(1979,*)"# N = trace(rho_jk) = ",NPar 00337 write(1979,*)"# MOrb = ",MOrb 00338 write(1979,*)"# " 00339 write(1979,*)"# AbsTime int_-infty^x rho^(1)(x,x) dx" 00340 close(1979) 00341 00342 END IF 00343 00344 open(unit=76,file=filename,form='FORMATTED',position = 'APPEND',status = 'OLD') 00345 write(76,'(2(F26.16,a1))')AbsTime," ",cumulativeProbability 00346 close(unit=76) 00347 00348 00349 counter = 1 00350 !---}}} 00351 END SUBROUTINE write_probabilityLeft 00352 00353 !---------------------------------------------------------------------------------------- 00359 SUBROUTINE get_energy(CVec,BandMat,energy,nConf) 00360 USE moduledefaults,ONLY:nSuperDiags 00361 !---{{{ 00362 IMPLICIT NONE 00363 INTEGER,INTENT(IN) :: nConf 00364 COMPLEX*16,INTENT(IN) :: CVec(nConf),BandMat(nSuperDiags+1,nConf) 00365 REAL*8,INTENT(OUT) :: energy 00366 COMPLEX*16 :: HCVec(nConf),cenergy 00367 00368 COMPLEX*16,parameter :: alpha=(1.d0,0.d0) 00369 COMPLEX*16,parameter :: beta=(0.d0,0.d0) 00370 INTEGER,parameter :: incx=1 00371 INTEGER,parameter :: incy=1 00372 00373 HCVec=(0.d0,0.d0) 00374 00375 CALL ZHBMV ( 'L', nConf, nSuperDiags, alpha, BandMat, nSuperDiags+1, CVec, incx, beta, HCVec, incy ) 00376 !CALL vvaxzz (CVec,HCVec,cenergy,nConf) 00377 cenergy = DOT_PRODUCT(Cvec,HCVec) 00378 00379 energy=DBLE(cenergy) 00380 00381 !!---}}} 00382 END SUBROUTINE get_energy 00383 00384 00385 !--------------------------------------------------------------------------- 00391 SUBROUTINE write_energy(AbsTime,lambda0,NPar,energy,identifier) 00392 !---{{{ 00393 IMPLICIT NONE 00394 REAL*8,INTENT(IN) :: AbsTime,lambda0,energy 00395 INTEGER,INTENT(IN) :: NPar 00396 INTEGER,save :: counter=0 00397 CHARACTER(255),INTENT(IN) :: identifier 00398 CHARACTER(255) :: filename 00399 00400 filename = 'energy'//trim(identifier)//'.dat' 00401 00402 IF (counter == 0 ) THEN 00403 00404 open(unit=1979,file=filename,form='FORMATTED',status='REPLACE') 00405 write(1979,*)"#" 00406 write(1979,*)"# the energy of the system as a function of time" 00407 write(1979,*)"#" 00408 write(1979,*)"# N = ",NPar 00409 write(1979,*)"# lambda0 = ",lambda0 00410 write(1979,*)"# lambda = lambda0*(N-1) = ",lambda0*(NPar-1.d0) 00411 write(1979,*)"#" 00412 write(1979,*)"# AbsTime energy = E = <Psi|H|Psi> E/N " 00413 close(1979) 00414 00415 00416 END IF 00417 00418 open(unit=76,file=filename, & 00419 form = 'FORMATTED', & 00420 position = 'APPEND', & 00421 status = 'OLD') 00422 write(76,'(3(F26.16,a1))')AbsTime," ",energy," ",energy/DBLE(NPar) 00423 close(unit=76) 00424 00425 counter = 1 00426 00427 !---}}} 00428 END SUBROUTINE write_energy 00429 00430 00431 00432 !--------------------------------------------------------------------------- 00438 SUBROUTINE write_mu_kj(AbsTime,mu_kj,identifier) 00439 USE moduleinputvariables,ONLY:MOrb 00440 !---{{{ 00441 IMPLICIT NONE 00442 REAL*8,INTENT(IN) :: AbsTime 00443 COMPLEX*16,INTENT(IN) :: mu_kj(MOrb,MOrb) 00444 CHARACTER(255),INTENT(IN) :: identifier 00445 INTEGER,save :: counter=0 00446 CHARACTER(255) :: filename 00447 INTEGER :: m,n 00448 00449 filename = 'mu_kj'//trim(identifier)//'.dat' 00450 00451 IF (counter == 0 ) THEN 00452 00453 open(unit=1979,file=filename,form='FORMATTED',status='REPLACE') 00454 write(1979,*)"#" 00455 write(1979,*)"# the chemical potentials mu_kj as a function of time" 00456 write(1979,*)"#" 00457 write(1979,*)"#" 00458 write(1979,*)"# AbsTime mu_kj " 00459 close(1979) 00460 00461 00462 END IF 00463 00464 open(unit=76,file=filename, & 00465 form = 'FORMATTED', & 00466 position = 'APPEND', & 00467 status = 'OLD') 00468 write(76,'(9(F26.16,a1))')AbsTime," ",( (DBLE(mu_kj(m,n))," ",IMAG(mu_kj(m,n))," ", m=1,MOrb),n=1,MOrb) 00469 close(unit=76) 00470 00471 counter = 1 00472 00473 !---}}} 00474 END SUBROUTINE write_mu_kj 00475 00476 00477 00478 00479 !---------------------------------------------------------------------------------- 00487 SUBROUTINE get_defect(nConf,nSuperDiags,BandMat,CVec,energy,defect) 00488 !---{{{ 00489 IMPLICIT NONE 00490 INTEGER,INTENT(IN) :: nConf,nSuperDiags 00491 COMPLEX*16,INTENT(IN) :: CVec(nConf),BandMat(nSuperDiags+1,nConf) 00492 REAL*8,INTENT(IN) :: energy 00493 00494 REAL*8,INTENT(OUT) :: defect 00495 00496 COMPLEX*16 :: cdefect=(0.d0,0.d0) 00497 COMPLEX*16 :: alpha,beta 00498 COMPLEX*16 :: HCVec(nConf) 00499 INTEGER :: incx,incy 00500 00501 incx=1 00502 incy=1 00503 alpha=(1.d0,0.d0) 00504 beta =(0.d0,0.d0) 00505 00506 HCVec=(0.d0,0.d0) 00507 00508 CALL ZHBMV ( 'L', nConf, nSuperDiags, alpha, BandMat, nSuperDiags+1, CVec, incx, beta,HCVec, incy) 00509 00510 cdefect = DOT_PRODUCT(HCVec,HCVec) 00511 00512 !get d = H|Psi> - E|Psi> 00513 HCVec = HCVec - energy * CVec 00514 00515 cdefect = DOT_PRODUCT(HCVec,HCVec) 00516 defect = SQRT(ABS(cdefect)) 00517 00518 !---}}} 00519 END SUBROUTINE get_defect 00520 00521 00522 !---------------------------------------------------------------------------------- 00530 SUBROUTINE get_meanDensity(Psi,rho_jk,meanDensity) 00531 USE moduleinputvariables 00532 ! computes Int{G1(x',x')^2,(x,-infty,x)} if Psi is provided {{{ 00533 00534 IMPLICIT NONE 00535 REAL*8,INTENT(OUT) :: meanDensity 00536 REAL*8 :: w 00537 COMPLEX*16,DIMENSION(NDX*NDY*NDZ,MOrb),INTENT(IN) :: Psi 00538 COMPLEX*16,DIMENSION(MOrb,MOrb),INTENT(IN) :: rho_jk 00539 INTEGER :: n,Iorb,Jorb 00540 COMPLEX*16 :: cdensxsquared 00541 00542 00543 IF((NDY.ne.1).or.(NDZ.ne.1)) THEN 00544 write(*,*)"ONLY 1d implemented" 00545 STOP 00546 END IF 00547 00548 w=DSQRT((xf-xi)/NDX) 00549 meanDensity=0.d0 00550 00551 DO n=1,NDX*NDY*NDZ 00552 cdensxsquared = (0.d0,0.d0) 00553 DO Jorb = 1,MOrb 00554 DO Iorb=1,MOrb 00555 00556 cdensxsquared = cdensxsquared + rho_jk(Iorb,Jorb) * DCONJG(Psi(n,Iorb))*Psi(n,Jorb) 00557 00558 END DO 00559 END DO 00560 meanDensity = meanDensity + DBLE(cdensxsquared)**2/w**2 00561 00562 END DO 00563 00564 meanDensity = meanDensity / Npar 00565 00566 !---}}} 00567 END SUBROUTINE get_meanDensity 00568 00569 !--------------------------------------------------------------------------- 00575 SUBROUTINE write_liebLinigerParameter(AbsTime,lambda0,NPar,meanDensity,identifier) 00576 !---{{{ 00577 IMPLICIT NONE 00578 REAL*8,INTENT(IN) :: AbsTime,lambda0,meanDensity 00579 INTEGER,INTENT(IN) :: NPar 00580 INTEGER,save :: counter=0 00581 CHARACTER(255),INTENT(IN) :: identifier 00582 CHARACTER(255) :: filename 00583 00584 filename = 'liebLinigerParameter'//trim(identifier)//'.dat' 00585 00586 00587 IF (counter == 0 ) THEN 00588 00589 open(unit=1979,file=filename,form='FORMATTED',status='REPLACE') 00590 write(1979,*)"#" 00591 write(1979,*)"# the Lieb Liniger parameter as a function of time" 00592 write(1979,*)"# cf The mathematics of the Bose-gas, Lieb, Seiringer, Solovej, Yngvason" 00593 write(1979,*)"# we use the definition that reduces to gamma of LL's 1963 paper here " 00594 write(1979,*)"# " 00595 write(1979,*)"# AbsTime gamma=lambda0/(int [rho^(1)(x,x)]^2 / N dx)" 00596 write(1979,*)"# N = ",NPar 00597 write(1979,*)"# lambda0 = ",lambda0 00598 write(1979,*)"# lambda = lambda0*(N-1) = ",lambda0*(NPar-1.d0) 00599 close(1979) 00600 00601 00602 END IF 00603 00604 open(unit=76,file=filename, & 00605 form = 'FORMATTED', & 00606 position = 'APPEND', & 00607 status = 'OLD') 00608 write(76,'(2(F26.16,a1))')AbsTime," ",lambda0/meanDensity 00609 close(unit=76) 00610 00611 counter = 1 00612 00613 !---}}} 00614 END SUBROUTINE write_liebLinigerParameter 00615 00616 !---------------------------------------------------------------------------------- 00622 SUBROUTINE get_naturalOrbitals(MOrb,NDX,NDY,NDZ,rho_jk,Psi,noccs,norbs) 00623 USE moduleparameters 00624 !{{{ 00625 00626 IMPLICIT NONE 00627 INTEGER,INTENT(IN) :: MOrb,NDX,NDY,NDZ 00628 COMPLEX*16,INTENT(IN) :: rho_jk(MOrb,MOrb) 00629 COMPLEX*16,INTENT(IN) :: Psi(NDX*NDY*NDZ,MOrb) 00630 REAL*8 ,INTENT(OUT) :: noccs(MOrb) 00631 COMPLEX*16,INTENT(OUT) :: norbs(NDX*NDY*NDZ,MOrb) 00632 00633 COMPLEX*16 :: A(MOrb,MOrb),WORK(2*MOrb) 00634 REAL*8 :: RWORK(max(1, 3*MOrb-2)) 00635 INTEGER :: INFO 00636 INTEGER :: m,n 00637 00638 RWORK = 0.d0 00639 INFO = 0 00640 00641 ! to get the noccs in descENDing order muliply rho_jk by -1 00642 00643 A = (-1.d0) * rho_jk 00644 CALL ZHEEV('V','U',MOrb,A,MOrb,noccs,WORK,2*MOrb,RWORK,INFO) 00645 IF(INFO.NE.0) THEN 00646 WRITE(*,*)"THE SHIT HAS HIT THE FAN IN get_naturalOrbitals" 00647 WRITE(*,*)"INFO",INFO 00648 STOP 00649 END IF 00650 00651 noccs = noccs * (-1.d0) 00652 ! on output the columns of A are the eigenvecors of rho_jk 00653 ! the eigenvalUSE are stored in noccs 00654 ! with U = A and lambda = noccs (diagonal matrix) 00655 ! rho = U lambda U^dagger 00656 ! the m-th norb is THEN norb(:,m) = Sum_j Psi(:,j) A^*_jm 00657 00658 norbs = ZERO 00659 00660 DO m=1,MOrb 00661 DO n=1,MOrb 00662 00663 norbs(:,m) = norbs(:,m) + Psi(:,n)*DCONJG(A(n,m)) 00664 00665 END DO 00666 END DO 00667 00668 !}}} 00669 END SUBROUTINE get_naturalOrbitals 00670 00671 !--------------------------------------------------------------------------------- 00678 SUBROUTINE write_PsiToFile(AbsTime,Psi,NDX,NDY,NDZ,MOrb,fileID,identifier) 00679 USE moduleallocatables,ONLY:gridx,gridy,gridz,V_ext 00680 ! {{{ 00681 IMPLICIT NONE 00682 REAL*8,INTENT(IN) :: AbsTime 00683 COMPLEX*16,INTENT(IN) :: Psi(NDX*NDY*NDZ,MOrb) 00684 INTEGER,INTENT(IN) :: fileID,NDX,NDY,NDZ,MOrb 00685 00686 CHARACTER (LEN=10) :: doubleAsString 00687 INTEGER :: outputunit 00688 INTEGER :: i,j,k,m,n 00689 00690 CHARACTER(255),INTENT(IN) :: identifier 00691 CHARACTER(255) :: filename 00692 00693 outputunit = 11 00694 CALL get_doubleAsString(AbsTime,doubleAsString) 00695 00696 00697 IF (fileID==0) THEN 00698 00699 filename=trim(doubleAsString//trim(identifier)//'orbs.dat') 00700 00701 ELSE 00702 00703 filename=trim('restart.orbs.out') 00704 00705 END IF 00706 00707 00708 open(unit=outputunit,file=trim(filename),form='formatted') 00709 DO k=1,NDZ ! 00710 DO j=1,NDY ! 00711 DO i=1,NDX ! format: x y z V_ext Re(Psi_1) Im(Psi_1),...,Re(Psi_M) Im(Psi_M) 00712 00713 write(11,*)gridx(i)," ",gridy(j)," ",gridz(k)," ",DBLE(V_ext(i+NDX*(j-1)+NDX*NDY*(k-1)))," ", & 00714 ( DBLE(Psi(i+NDX*(j-1)+NDX*NDY*(k-1),m))," ",& 00715 DIMAG(Psi(i+NDX*(j-1)+NDX*NDY*(k-1),m)), m=1,MOrb) 00716 00717 END DO 00718 END DO 00719 END DO 00720 00721 close(outputunit) 00722 00723 ! }}} 00724 END SUBROUTINE write_PsiToFile 00725 00726 !--------------------------------------------------------------------------- 00732 SUBROUTINE write_norbsToFile(AbsTime,norbs,NDX,NDY,NDZ,MOrb,identifier) 00733 USE moduleallocatables,ONLY:gridx,gridy,gridz,V_ext 00734 ! {{{ write an output file. all vectors contain the weights 00735 IMPLICIT NONE 00736 REAL*8,INTENT(IN) :: AbsTime 00737 COMPLEX*16,INTENT(IN) :: norbs(NDX*NDY*NDZ,MOrb) 00738 INTEGER,INTENT(IN) :: NDX,NDY,NDZ,MOrb 00739 CHARACTER(255),INTENT(IN) :: identifier 00740 CHARACTER (LEN=10) :: doubleAsString 00741 CHARACTER(255) :: filename 00742 00743 INTEGER :: outputunit 00744 INTEGER :: i,j,k,m,n 00745 00746 outputunit = 11 00747 CALL get_doubleAsString(AbsTime,doubleAsString) 00748 00749 filename=trim(doubleAsString//trim(identifier)//'naturalorbitals.dat') 00750 00751 open(unit=outputunit,file=trim(filename),form='formatted') 00752 DO k=1,NDZ ! 00753 DO j=1,NDY ! 00754 DO i=1,NDX ! format: x y z V_ext Re(norbs_1) Im(norbs_1),...,Re(norbs_M) Im(norbs_M) 00755 00756 write(outputunit,*)gridx(i)," ",gridy(j)," ",gridz(k)," ",DBLE(V_ext(i+NDX*(j-1)+NDX*NDY*(k-1)))," ", & 00757 ( DBLE(norbs(i+NDX*(j-1)+NDX*NDY*(k-1),m))," ",& 00758 DIMAG(norbs(i+NDX*(j-1)+NDX*NDY*(k-1),m)), m=1,MOrb) 00759 00760 END DO 00761 END DO 00762 END DO 00763 00764 close(outputunit) 00765 00766 ! }}} 00767 END SUBROUTINE write_norbsToFile 00768 00769 !--------------------------------------------------------------------------- 00777 SUBROUTINE write_CVecFile(AbsTime,CVec,fileID,identifier) 00778 USE moduledefaults,ONLY:nConf 00779 ! {{{ 00780 IMPLICIT NONE 00781 REAL*8,INTENT(IN) :: AbsTime 00782 COMPLEX*16,INTENT(IN) :: CVec(nConf) 00783 INTEGER,INTENT(IN) :: fileID 00784 CHARACTER (LEN=255),INTENT(IN) :: identifier 00785 00786 CHARACTER (LEN=10) :: doubleAsString 00787 CHARACTER (LEN=50) :: filename 00788 00789 INTEGER :: outputunit 00790 INTEGER :: i 00791 00792 outputunit = 11 00793 00794 CALL get_doubleAsString(AbsTime,doubleAsString) 00795 00796 IF (fileID==0) THEN 00797 00798 filename=trim(doubleAsString//trim(identifier)//'coef.dat') 00799 00800 ELSE 00801 00802 filename=trim('restart.coef.out') 00803 00804 END IF 00805 00806 open(unit=outputunit,file=trim(filename),form='formatted') 00807 00808 DO i=1,nConf ! format: i Re(CVec(i)) Im(CVec(i)) 00809 IF (ABS(CVec(i))>0.5d-14) THEN 00810 write(11,912)i," ",dreal(CVec(i))," ",dimag(CVec(i)) 00811 END IF 00812 END DO 00813 00814 close(outputunit) 00815 00816 912 format(I12,a3,f18.15,a3,f18.15) 00817 ! }}} 00818 END SUBROUTINE write_CVecFile 00819 00820 !--------------------------------------------------------------------------- 00829 SUBROUTINE read_potentialFromFile(V_ext,NDX,NDY,NDZ) 00830 ! {{{ 00831 IMPLICIT NONE 00832 INTEGER,INTENT(IN) :: NDX,NDY,NDZ 00833 REAL*8,INTENT(OUT) :: V_ext(NDX*NDY*NDZ) 00834 INTEGER :: inputunit,ierr 00835 INTEGER :: i,j,k,m,nvals 00836 REAL*8 :: x,y,z 00837 00838 inputunit = 11 00839 ierr = 0 00840 x = 0.d0 00841 y = 0.d0 00842 z = 0.d0 00843 V_ext = 0.d0 00844 nvals = 0 00845 00846 WRITE(*,100)"--> Reading external potential from file V_ext.inp" 00847 00848 OPEN(UNIT=inputunit,FILE='V_ext.inp',FORM='formatted',ACTION='READ',IOSTAT=ierr) 00849 00850 openif: IF (ierr == 0) THEN 00851 00852 readloop: DO k=1,NDZ 00853 DO j=1,NDY 00854 DO i=1,NDX ! format: x,y,z, V(x,y,z) 00855 00856 READ(inputunit,*,IOSTAT=ierr)x,y,z,V_ext(i+NDX*(j-1)+NDX*NDY*(k-1)) 00857 IF(ierr /= 0) EXIT 00858 nvals = nvals + 1 00859 00860 END DO 00861 END DO 00862 END DO readloop 00863 00864 readif: IF (ierr > 0) THEN 00865 WRITE(*,1020) nvals+1 00866 1020 FORMAT ('0','AN ERROR OCCURED READING LINE:',I6) 00867 STOP 00868 ELSE IF (ierr < 0) THEN 00869 WRITE(*,1030) nvals,NDX*NDY*NDZ 00870 1030 FORMAT(' ','TRIED READING BEYOND THE LAST LINE:',I6,', EXPECTED',I6,' LINES') 00871 STOP 00872 ELSE 00873 CONTINUE 00874 END IF readif 00875 00876 00877 ELSE openif 00878 WRITE(*,1040) ierr 00879 1040 FORMAT (' ','Error opening file: IOSTAT = ', I6) 00880 END IF openif 00881 00882 CLOSE(inputunit) 00883 00884 !WRITE(*,100)"<-- Reading external potential successful" 00885 00886 100 FORMAT(5X,A) 00887 ! }}} 00888 END SUBROUTINE read_potentialFromFile 00889 00890 !--------------------------------------------------------------------------- 00896 SUBROUTINE read_PsiFromFile(Psi,NDX,NDY,NDZ,MOrb) 00897 USE moduleparameters 00898 ! {{{ read the orbitals from a file and orthonormalize them 00899 00900 IMPLICIT NONE 00901 INTEGER,INTENT(IN) :: NDX,NDY,NDZ,MOrb 00902 COMPLEX*16,INTENT(OUT) :: Psi(NDX*NDY*NDZ,MOrb) 00903 INTEGER :: inputunit,ierr 00904 INTEGER :: i,j,k,m,n,nvals 00905 REAL*8 :: x,y,z,V_xyz,RePsi(MOrb),ImPsi(MOrb),norm 00906 00907 inputunit = 11 00908 ierr = 0 00909 x = 0.d0 00910 y = 0.d0 00911 z = 0.d0 00912 V_xyz = 0.d0 00913 RePsi = 0.d0 00914 ImPsi = 0.d0 00915 nvals = 0 00916 norm = 0.d0 00917 00918 WRITE(*,100)"--> Reading orbitals from restart.orbs.inp" 00919 00920 OPEN(UNIT=inputunit,FILE='restart.orbs.inp',FORM='formatted',ACTION='READ',IOSTAT=ierr) 00921 00922 openif: IF (ierr == 0) THEN 00923 00924 readloop: DO k=1,NDZ ! 00925 DO j=1,NDY ! 00926 DO i=1,NDX ! format: x,y,z, V(x,y,z) Re(Psi_1) Im(Psi_1),...,Re(Psi_MOrb) Im(Psi_MOrb) 00927 00928 READ(inputunit,*,IOSTAT=ierr)x,y,z,V_xyz,(RePsi(m), ImPsi(m),m=1,MOrb) 00929 00930 IF(ierr /= 0) THEN 00931 WRITE(*,*) x,y,z,V_xyz 00932 EXIT 00933 END IF 00934 nvals = nvals + 1 00935 00936 DO m=1,MOrb 00937 Psi(i+NDX*(j-1)+NDX*NDY*(k-1),m) = RePsi(m) + ZONEI * ImPsi(m) 00938 END DO 00939 00940 END DO 00941 END DO 00942 END DO readloop 00943 00944 readif: IF (ierr > 0) THEN 00945 WRITE(*,1020) nvals+1 00946 1020 FORMAT ('0',' IN read_PsiFromFile: AN ERROR OCCURED READING LINE',I6) 00947 STOP 00948 ELSE IF (ierr < 0) THEN 00949 WRITE(*,1030) nvals,NDX*NDY*NDZ 00950 1030 FORMAT(' ',' IN read_PsiFromFile: TRIED READING BEYOND THE LAST LINE OF FILE.',I6, & 00951 ' LINES READ IN, BUT EXPECTED',I6) 00952 STOP 00953 ELSE 00954 CONTINUE 00955 END IF readif 00956 00957 00958 ELSE openif 00959 WRITE(*,1040) ierr 00960 1040 FORMAT (' ','Error opening file: IOSTAT = ', I6) 00961 END IF openif 00962 CLOSE(inputunit) 00963 00964 DO m=1,MOrb 00965 CALL normvxz(Psi(:,m),norm,NDX*NDY*NDZ) 00966 WRITE(*,105)"orbital",m," has norm ",norm 00967 END DO 00968 00969 00970 WRITE(*,110)"Gram-Schmidt orthonormalization of the orbitals" 00971 CALL schmidtortho (Psi,NDX*NDY*NDZ,MOrb,ierr) 00972 IF (Ierr.NE.0) THEN 00973 WRITE(*,*)"orthogonalization failed" 00974 END IF 00975 00976 !WRITE(*,100)"<-- Reading orbitals succesfull" 00977 100 FORMAT (5X,A) 00978 105 FORMAT (9X,A,I3,A,F16.8) 00979 110 FORMAT (9X,A) 00980 00981 ! }}} 00982 END SUBROUTINE read_PsiFromFile 00983 00984 !--------------------------------------------------------------------------- 00990 SUBROUTINE read_CVecFile(CVec,nConf) 00991 ! {{{ 00992 IMPLICIT NONE 00993 INTEGER,INTENT(IN) :: nConf 00994 COMPLEX*16,INTENT(OUT) :: CVec(nConf) 00995 REAL*8 :: RE_CVec(nConf),IM_CVec(nConf),norm 00996 REAL*8 :: RE,IM 00997 00998 INTEGER :: inputunit 00999 INTEGER :: i,j 01000 01001 WRITE(*,90)"--> Reading the coefficient vector" 01002 01003 CVec = (0.d0,0.d0) 01004 RE_CVec = 0.0d0 01005 IM_CVec = 0.0d0 01006 IM = 0.0D0 01007 RE = 0.0D0 01008 inputunit = 11 01009 norm = 0.d0 01010 01011 01012 open(unit=inputunit,file='restart.coef.inp',form='formatted') 01013 01014 DO i=1,nConf ! format: i Re(CVec(i)) Im(CVec(i)) 01015 read(inputunit,*,END=21)j,RE,IM 01016 IF (j<=nConf) THEN 01017 RE_CVec(j)=RE 01018 IM_CVec(j)=IM 01019 ELSE 01020 WRITE(*,*)"THE VECTOR CVec HAS ONLY nConf CONFIGURATIONS, nConf:",nConf 01021 WRITE(*,*)"ERROR, I TRIED TO READ CONFIGURATION NO ",j 01022 STOP 01023 END IF 01024 END DO 01025 21 continue 01026 01027 close(inputunit) 01028 01029 01030 CVec = RE_CVec + (0.d0,1.d0)*IM_CVec 01031 CALL normvxz(CVec,norm,nConf) 01032 WRITE(*,100)" Norm of read in coefficient vector =",norm 01033 WRITE(*,100)" Normalize the coefficient vector" 01034 01035 CVec = CVec/norm 01036 CALL normvxz(CVec,norm,nConf) 01037 01038 90 FORMAT(5X,A) 01039 100 FORMAT(5X,A,F16.12) 01040 01041 !WRITE(*,90)"<-- Reading coefficient vector successful" 01042 ! }}} 01043 END SUBROUTINE read_CVecFile 01044 01045 !--------------------------------------------------------------------------- 01051 SUBROUTINE write_noccsToFile(AbsTime,MOrb,noccs,identifier) 01052 !---{{{ 01053 IMPLICIT NONE 01054 INTEGER, INTENT(IN) :: MOrb 01055 REAL*8, INTENT(IN) :: AbsTime,noccs(MOrb) 01056 CHARACTER(len=255),INTENT(IN) :: identifier 01057 01058 INTEGER,save :: counter=0 01059 CHARACTER (len=255) :: fullName,varfmt,MOrbString 01060 CHARACTER (len=6) :: posString 01061 REAL*8 :: trace 01062 INTEGER :: m 01063 01064 01065 trace = 0.d0 01066 01067 DO m=1,MOrb 01068 trace = trace + noccs(m) 01069 END DO 01070 01071 fullName = 'naturalOccupations'//trim(identifier)//'.dat' 01072 01073 IF (MOrb<9 .AND. MOrb>0) THEN 01074 write(MOrbString,'(I1)') MOrb+1 01075 varfmt = "("//trim(MOrbString)//"(F20.10,a1))" 01076 ELSE IF (MOrb<99 .AND. MOrb > 9) THEN 01077 write(MOrbString,'(I2)') MOrb+1 01078 varfmt = "("//trim(MOrbString)//"(F20.10,a1))" 01079 ELSE 01080 WRITE(*,*)"TOO MANY/FEW ORBITALS: MOrb",MOrb 01081 STOP 01082 END IF 01083 01084 IF(counter==0) THEN 01085 open(unit=12,file=trim(fullName),form='formatted',POSITION='REWIND') 01086 write(12,*)"# (fractional) natural orbital occupations" 01087 write(12,*)"# " 01088 write(12,*)"# N = trace(rho_jk) = ",trace 01089 write(12,*)"# MOrb = ",MOrb 01090 write(12,*)"# " 01091 write(12,*)"# time ",( "n_m/N ", m=1,MOrb) 01092 close(12) 01093 END IF 01094 01095 open(unit=12,file=fullName,form='formatted',POSITION='APPEND') 01096 write(12,FMT=trim(varfmt)) AbsTime," ",(noccs(m)/trace, " ", m=1,MOrb) 01097 close(12) 01098 01099 counter=1 01100 !}}} 01101 END SUBROUTINE write_noccsToFile 01102 01103 !--------------------------------------------------------------------------- 01110 SUBROUTINE write_firstFragmentationTimes(AbsTime,MOrb,NPar,noccs,lambda0,identifier) 01111 !---{{{ 01112 IMPLICIT NONE 01113 INTEGER, INTENT(IN) :: MOrb,NPar 01114 REAL*8, INTENT(IN) :: AbsTime,noccs(MOrb),lambda0 01115 CHARACTER (len=255),INTENT(IN) :: identifier 01116 01117 INTEGER,save :: counter=0,fragcounter=0 01118 CHARACTER (len=255) :: fullName 01119 CHARACTER (len=6) :: posString 01120 REAL*8 :: frag,increment,trace,lambda 01121 INTEGER :: m 01122 01123 trace = 0.d0 01124 DO m =1,MOrb 01125 trace = trace + noccs(m) 01126 END DO 01127 01128 lambda = (NPar-1.d0)*lambda0 01129 01130 increment=0.002d0 01131 01132 fullName = 'firstFragmentationTimes'//trim(identifier)//'.dat' 01133 01134 frag = 1.d0 - noccs(1)/trace ! the fragmentation in percent 1.0 - n^(1)_1 / N 01135 01136 posString='APPEND' 01137 IF(counter==0) THEN 01138 posString='REWIND' 01139 01140 open(unit=12,file=trim(fullName),form='formatted',POSITION=posString) 01141 write(12,*)"# times at which a certain fragmetation is reached for the first time" 01142 write(12,*)"# " 01143 write(12,*)"# trace(rho_jk) = ",trace 01144 write(12,*)"# NPar = ",NPar 01145 write(12,*)"# " 01146 write(12,*)"# fragmentation tbreak N lambda=lambda0*(N-1)" 01147 write(12,222) frag," ",AbsTime," ",DBLE(NPar)," ",lambda 01148 close(12) 01149 posString='APPEND' 01150 END IF 01151 01152 IF (frag>increment*(fragcounter+1)) THEN 01153 01154 while_loop: DO 01155 01156 fragcounter=fragcounter+1 01157 open(unit=12,file=fullName,form='formatted',POSITION=posString) 01158 write(12,222) fragcounter*increment," ",AbsTime," ",DBLE(NPar)," ",lambda 01159 close(12) 01160 01161 IF (frag<(fragcounter +1)*increment) exit while_loop 01162 01163 END DO while_loop 01164 01165 END IF 01166 01167 01168 01169 222 format(4(F22.13,a1)) 01170 counter=1 01171 !}}} 01172 END SUBROUTINE write_firstFragmentationTimes 01173 01174 END MODULE moduleinputoutput
1.7.3
