|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 ! vim:fdm=marker: 00002 00003 !----------------------------------------------------------------------- 00010 MODULE modulecorrelationfunctions 00011 00012 IMPLICIT NONE 00013 CONTAINS 00021 subroutine write_correlations(Dims,printCorrelationFuncs,MOMSPACE2D,REALSPACE2D,& 00022 x1const,x1slice,y1const,y1slice,& 00023 x2const,x2slice,y2const,y2slice,& 00024 AbsTime,Psi,rho_jk,rho_ijkl,& 00025 Pnot,xstart,xstop) 00026 USE moduleinputvariables,ONLY:MOrb,NDX,NDY,NDZ 00027 00028 IMPLICIT NONE 00029 00030 INTEGER,intent(in) :: Dims,printCorrelationFuncs 00031 LOGICAL,intent(in) :: MOMSPACE2D, REALSPACE2D 00032 LOGICAL,intent(in) :: x1const, x2const, y1const, y2const 00033 LOGICAL,intent(in) :: Pnot 00034 REAL*8,intent(in) :: xstart,xstop 00035 REAL*8,intent(in) :: x1slice, x2slice, y1slice, y2slice 00036 REAL*8,intent(in) :: AbsTime 00037 COMPLEX*16, DIMENSION(MORB*MORB),intent(in) :: rho_jk 00038 COMPLEX*16, DIMENSION(MORB*MORB*MORB*MORB),intent(in) :: rho_ijkl 00039 COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MORB),intent(in) :: Psi 00040 00041 COMPLEX*16, DIMENSION(NDX*NDY*NDZ,MORB) :: FTPsi 00042 INTEGER :: dilation 00043 00044 dilation=1 00045 00046 IF(Dims==1) THEN 00047 if (printCorrelationFuncs==1) then 00048 00049 CALL get_correlations(AbsTime,Psi,1,1,rho_jk,rho_ijkl) 00050 00051 CALL pedestrian_FT(Psi,FTPsi,dilation=1) 00052 CALL get_dilation(dilation,FTPsi,rho_jk) 00053 CALL pedestrian_FT(Psi,FTPsi,dilation) 00054 CALL get_oneBodyRDMDiagonal(AbsTime,FTPsi,rho_jk,2,dilation) 00055 CALL get_correlations(AbsTime,FTPsi,2,dilation,rho_jk,rho_ijkl) 00056 if (Pnot.eqv..TRUE.) then 00057 CALL get_densityNonescape(xstart,xstop,rho_jk,AbsTime,Psi) 00058 endif 00059 end if 00060 ELSE IF (Dims == 2) THEN 00061 if (MOMSPACE2D.eqv..TRUE.) then 00062 call FT_SORTED(Psi,FTPSI) 00063 call get_corr_slice(Abstime,FTPSI,2,rho_jk,rho_ijkl,& 00064 x1const,x1slice,y1const,y1slice,x2const,x2slice,y2const,y2slice) 00065 endif 00066 if (REALSPACE2D.eqv..TRUE.) then 00067 call get_corr_slice(Abstime,PSI,1,rho_jk,rho_ijkl,& 00068 x1const,x1slice,y1const,y1slice,x2const,x2slice,y2const,y2slice) 00069 endif 00070 END IF 00071 end subroutine write_Correlations 00072 00073 00083 subroutine get_integratedDensity(rho_jk,AbsTime,PSI,writemode) 00084 00085 USE moduleinputoutput,ONLY: get_doubleAsString 00086 USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb 00087 USE moduleparameters,ONLY: PI 00088 USE moduleallocatables,ONLY:gridx,gridy,gridz 00089 00090 IMPLICIT NONE 00091 00092 COMPLEX*16, DIMENSION(NDX*NDY*NDZ,Morb),INTENT(IN) :: PSI 00093 complex*16,DIMENSION(Morb,Morb),intent(in) :: rho_jk 00094 real*8,intent(in) :: AbsTime 00095 INTEGER,INTENT(IN):: writemode 00096 00097 complex*16, DIMENSION(NDX*NDY*NDZ) :: DENSITY 00098 complex*16, DIMENSION(NDX,NDY,NDZ) :: int_density 00099 character(len=10) :: timeAsString 00100 00101 integer :: bit,NOGRIDPOINTS,i,j,k,n 00102 logical :: integratex,integratey,integratez 00103 00104 IF(writemode>7 .OR. writemode<0) THEN 00105 WRITE(*,*)'stupid value for writemode',writemode 00106 STOP 00107 END IF 00108 00109 call get_DoubleAsString(AbsTime,timeAsString) 00110 00111 int_density = (0.d0,0.d0) 00112 integratex=.FALSE. 00113 integratey=.FALSE. 00114 integratez=.FALSE. 00115 00116 bit=writemode 00117 IF(MOD(bit,2)==1) integratez=.TRUE. 00118 bit=bit/2 00119 IF(MOD(bit,2)==1) integratey=.TRUE. 00120 bit=bit/2 00121 IF(MOD(bit,2)==1) integratex=.TRUE. 00122 00123 NOGRIDPOINTS=NDX*NDY*NDZ 00124 ! get the density from the one-body-matrix elements 00125 DENSITY = (0.d0,0.d0) 00126 do n=1,NOGRIDPOINTS 00127 do j = 1,Morb 00128 do k=1,MOrb 00129 DENSITY(n) = DENSITY(n) + rho_jk(j,k) * DCONJG(PSI(n,j))*PSI(n,k) 00130 end do 00131 end do 00132 end do 00133 00134 int_density=(0.d0,0.d0) 00135 00136 IF (integratez) THEN 00137 00138 DO i=1,NDX 00139 DO j=1,NDY 00140 DO k=1,NDZ 00141 int_density(i,j,1)=int_density(i,j,1) + DENSITY(i+NDX*(j-1)+NDX*NDY*(k-1)) 00142 END DO 00143 END DO 00144 END DO 00145 00146 END IF 00147 00148 00149 IF (integratey) THEN 00150 00151 DO i=1,NDX 00152 DO k=1,NDZ 00153 DO j=1,NDY 00154 int_density(i,1,k)=int_density(i,1,k) + DENSITY(i+NDX*(j-1)+NDX*NDY*(k-1)) 00155 END DO 00156 END DO 00157 END DO 00158 00159 END IF 00160 00161 00162 IF (integratex) THEN 00163 00164 DO k=1,NDZ 00165 DO j=1,NDY 00166 DO i=1,NDX 00167 int_density(1,j,k)=int_density(1,j,k) + DENSITY(i+NDX*(j-1)+NDX*NDY*(k-1)) 00168 END DO 00169 END DO 00170 END DO 00171 00172 END IF 00173 00174 00175 00176 00177 00178 00179 6667 format(5(F16.9,a1)) 00180 00181 00182 00183 IF (writemode==0) THEN 00184 00185 OPEN(unit=12,file=timeAsString//'xyzdensity.dat',form='formatted') 00186 00187 do k=1,NDZ 00188 do j=1,NDY 00189 do i=1,NDX 00190 write(12,6667) gridx(i)," ",gridy(j)," ",gridz(k)," ",DBLE(int_density(i,j,k))," ",AbsTime 00191 00192 if(MOD(i,NDX)==0) write(12,*)" " 00193 00194 end do 00195 end do 00196 end do 00197 00198 CLOSE(12) 00199 00200 ELSE IF (writemode==1) THEN 00201 00202 OPEN(unit=12,file=timeAsString//'xydensity.dat',form='formatted') 00203 00204 do j=1,NDY 00205 do i=1,NDX 00206 00207 write(12,6667) gridx(i)," ",gridy(j)," ",gridz(1)," ",DBLE(int_density(i,j,1))," ",AbsTime 00208 00209 if(MOD(i,NDX)==0) write(12,*)" " 00210 00211 end do 00212 end do 00213 00214 CLOSE(12) 00215 00216 00217 00218 ELSE IF (writemode==2) THEN 00219 00220 OPEN(unit=12,file=timeAsString//'xzdensity.dat',form='formatted') 00221 00222 do k=1,NDZ 00223 do i=1,NDX 00224 00225 write(12,6667) gridx(i)," ",gridy(1)," ",gridz(k)," ",DBLE(int_density(i,1,k))," ",AbsTime 00226 00227 if(MOD(i,NDX)==0) write(12,*)" " 00228 00229 end do 00230 end do 00231 00232 CLOSE(12) 00233 00234 00235 ELSE IF (writemode==3) THEN 00236 00237 OPEN(unit=12,file=timeAsString//'xdensity.dat',form='formatted') 00238 00239 do i=1,NDX 00240 00241 write(12,6667) gridx(i)," ",gridy(1)," ",gridz(1)," ",DBLE(int_density(i,1,1))," ",AbsTime 00242 00243 if(MOD(i,NDX)==0) write(12,*)" " 00244 00245 end do 00246 00247 CLOSE(12) 00248 00249 ELSE IF (writemode==4) THEN 00250 00251 OPEN(unit=12,file=timeAsString//'yzdensity.dat',form='formatted') 00252 00253 do k=1,NDZ 00254 do j=1,NDY 00255 00256 write(12,6667) gridx(1)," ",gridy(j)," ",gridz(k)," ",DBLE(int_density(1,j,k))," ",AbsTime 00257 00258 if(MOD(j,NDY)==0) write(12,*)" " 00259 00260 end do 00261 end do 00262 00263 CLOSE(12) 00264 00265 ELSE IF (writemode==5) THEN 00266 00267 OPEN(unit=12,file=timeAsString//'ydensity.dat',form='formatted') 00268 00269 do j=1,NDY 00270 00271 write(12,6667) gridx(1)," ",gridy(j)," ",gridz(1)," ",DBLE(int_density(1,j,1))," ",AbsTime 00272 00273 if(MOD(j,NDY)==0) write(12,*)" " 00274 00275 end do 00276 00277 CLOSE(12) 00278 00279 ELSE IF (writemode==6) THEN 00280 00281 OPEN(unit=12,file=timeAsString//'zdensity.dat',form='formatted') 00282 00283 do k=1,NDZ 00284 00285 write(12,6667) gridx(1)," ",gridy(1)," ",gridz(k)," ",DBLE(int_density(1,1,k))," ",AbsTime 00286 00287 if(MOD(k,NDZ)==0) write(12,*)" " 00288 00289 end do 00290 00291 CLOSE(12) 00292 00293 ELSE IF (writemode==7) THEN 00294 00295 OPEN(unit=12,file=timeAsString//'density.dat',form='formatted') 00296 00297 write(12,6667) gridx(1)," ",gridy(1)," ",gridz(1)," ",DBLE(int_density(1,1,1))," ",AbsTime 00298 00299 CLOSE(12) 00300 00301 END IF 00302 00303 00304 00305 end subroutine get_integratedDensity 00306 00307 00313 00314 subroutine get_densityNonescape(xstart,xstop,rho_jk,time,PSI) 00315 00316 USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb 00317 USE moduleparameters,ONLY: PI 00318 USE moduleallocatables,ONLY:gridx,gridy,gridz 00319 00320 IMPLICIT NONE 00321 00322 COMPLEX*16, DIMENSION(NDX*NDY*NDZ,Morb) :: PSI 00323 complex*16,DIMENSION(Morb,Morb),intent(in) :: rho_jk 00324 real*8, intent(in) :: xstart,xstop 00325 real*8 :: nonescape,time,w 00326 00327 complex*16, DIMENSION(NDX*NDY*NDZ) :: DENSITY 00328 00329 integer :: i,j,k,l,m,n,NOGRIDPOINTS 00330 00331 00332 NOGRIDPOINTS=NDX*NDY*NDZ 00333 ! get the density from the one-body-matrix elements 00334 DENSITY = (0.d0,0.d0) 00335 do n=1,NOGRIDPOINTS 00336 do j = 1,Morb 00337 do k=1,MOrb 00338 DENSITY(n) = DENSITY(n) + rho_jk(j,k) * DCONJG(PSI(n,j))*PSI(n,k) 00339 end do 00340 end do 00341 end do 00342 00343 ! integrate from DENSITY from xi to xf 00344 nonescape=0.d0 00345 do n=1,NOGRIDPOINTS 00346 if ((gridx(n).le.xstop).and.(gridx(n).ge.xstart)) then 00347 nonescape=nonescape +dble(DENSITY(n)) 00348 endif 00349 end do 00350 ! open/write/close output 00351 open(unit=542,file='nonescape',form='formatted',ACCESS='APPEND') 00352 00353 write(542,6545) time,' ',nonescape 00354 00355 6545 format(2(F21.16,a1)) 00356 00357 close(542) 00358 end subroutine get_densityNonescape 00359 00360 00361 00362 00363 !----------------------------------------------------------------------- 00372 subroutine get_maximum(Psi,rho_jk,maximum) 00373 USE moduleinputvariables,ONLY:NDX,NDY,NDZ,MOrb 00374 !---{{{ 00375 00376 implicit none 00377 integer :: noGridPoints 00378 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(in) :: Psi 00379 complex*16,DIMENSION(MOrb,MOrb),intent(in) :: rho_jk 00380 complex*16,intent(out) :: maximum 00381 integer :: m,Jorb,Korb 00382 complex*16 :: G1 00383 noGridPoints=NDX*NDY*NDZ 00384 00385 00386 maximum = (0.d0,0.d0) 00387 do m=1,noGridPoints 00388 00389 !compute G1 00390 G1 = (0.d0,0.d0) 00391 do Korb = 1,MOrb 00392 do Jorb=1,MOrb 00393 00394 G1 = G1 + rho_jk(Jorb,Korb) * DCONJG(Psi(m,Jorb))*Psi(m,Korb) 00395 00396 end do 00397 end do 00398 00399 if (REAL(CDABS(G1)) > REAL(CDABS(maximum))) maximum = G1 00400 00401 end do 00402 !---}}} 00403 end subroutine get_maximum 00404 00405 !---------------------------------------------------------------------------- 00426 subroutine get_oneBodyRDMDiagonal(time,Psi,rho_jk,mode,dilation) 00427 USE modulederivesystemvariables, ONLY: weight 00428 USE moduleinputvariables 00429 USE moduleinputoutput,ONLY: get_doubleAsString 00430 USE moduleparameters,ONLY: PI 00431 USE moduleallocatables,ONLY:gridx,gridy,gridz 00432 USE modulegrid 00433 !---{{{ 00434 00435 implicit none 00436 integer :: noGridPoints 00437 integer,intent(in) :: mode,dilation 00438 real*8 :: time,w,dk,mom_X(NDX),mom_Y(NDY),mom_Z(NDZ) 00439 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(in) :: Psi 00440 complex*16,DIMENSION(MOrb,MOrb),intent(in) :: rho_jk 00441 character(len=2) :: aString 00442 character(len=200) :: fullName 00443 character(len=10) :: timeAsString 00444 character(len=10) :: MOrbAsString 00445 character(len=1) :: NparAsString1 00446 character(len=2) :: NparAsString2 00447 character(len=3) :: NparAsString3 00448 character(len=4) :: NparAsString4 00449 character(len=5) :: NparAsString5 00450 character(len=6) :: NparAsString6 00451 character(len=7) :: NparAsString7 00452 integer :: n,Jorb,Korb,i_n,j_n,k_n 00453 complex*16 :: G1n_n 00454 complex*16 :: maximum 00455 real*8, parameter :: TOL=0.00001d0 00456 real*8 :: threshold 00457 real*8 :: outG1n_n 00458 00459 00460 noGridPoints=NDX*NDY*NDZ 00461 00462 00463 if((NDY>1).or. (NDZ>1)) then 00464 ! write(*,*)"test of 2D and 3D still necessary" 00465 ! STOP 00466 end if 00467 00468 00469 if(mode==1) then 00470 w = weight 00471 else if (mode==2) then 00472 dk = (2.d0*PI)**Dims/ (weight**2 * noGridPoints)/dilation 00473 w = DSQRT( dk ) 00474 00475 CALL get_momentumGrid(xi,xf,NDX,mom_X) !returns momentum grid in strictly ascending order 00476 CALL get_momentumGrid(yi,yf,NDY,mom_Y) 00477 CALL get_momentumGrid(zi,zf,NDZ,mom_Z) 00478 00479 mom_X=mom_X/dilation 00480 mom_Y=mom_Y/dilation 00481 mom_Z=mom_Z/dilation 00482 00483 else 00484 write(*,*) "wrong mode" 00485 end if 00486 00487 write (MOrbAsString, '(I1)') MOrb 00488 call get_doubleAsString(time,timeAsString) 00489 00490 if (mode==1) then 00491 aString = 'x-' 00492 else if (mode==2) then 00493 aString = 'k-' 00494 else 00495 write(*,*)"ERROR in get_correlations" 00496 STOP 00497 endif 00498 if ((NPar>0) .and. (NPar.lt.10)) then 00499 write (NParAsString1, '(I1)') NPar 00500 fullName=timeasString//'N'//trim(NParAsString1)//'M'//trim(MOrbAsString)//aString//'density.dat' 00501 else if ((NPar.ge.10) .and. (NPar.lt.100)) then 00502 write (NParAsString2, '(I2)') NPar 00503 fullName=timeasString//'N'//trim(NParAsString2)//'M'//trim(MOrbAsString)//aString//'density.dat' 00504 else if ((NPar.ge.100) .and. (NPar.lt.1000)) then 00505 write (NParAsString3, '(I3)') NPar 00506 fullName=timeasString//'N'//trim(NParAsString3)//'M'//trim(MOrbAsString)//aString//'density.dat' 00507 else if ((NPar.ge.1000) .and. (NPar.lt.10000)) then 00508 write (NParAsString4, '(I4)') NPar 00509 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MOrbAsString)//aString//'density.dat' 00510 else if ((NPar.ge.10000) .and. (NPar.lt.100000)) then 00511 write (NParAsString5, '(I5)') NPar 00512 fullName=timeasString//'N'//trim(NParAsString5)//'M'//trim(MOrbAsString)//aString//'density.dat' 00513 else if ((NPar.ge.100000) .and. (NPar.lt.1000000)) then 00514 write (NParAsString6, '(I6)') NPar 00515 fullName=timeasString//'N'//trim(NParAsString6)//'M'//trim(MOrbAsString)//aString//'density.dat' 00516 else if ((NPar.ge.1000000) .and. (NPar.lt.10000000)) then 00517 write (NParAsString7, '(I7)') NPar 00518 fullName=timeasString//'N'//trim(NParAsString7)//'M'//trim(MOrbAsString)//aString//'density.dat' 00519 else 00520 NParAsString4='XXXX' 00521 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MOrbAsString)//aString//'density.dat' 00522 end if 00523 ! needs to be edited for other dvrs with nonconstant weights! 00524 call get_maximum(Psi,rho_jk,maximum) 00525 threshold = 0.01d0 * DBLE(maximum)/w**2 ! set the minimal density to be plotted 00526 00527 open(unit=12,file=trim(fullName),form='formatted') 00528 do n=1,noGridPoints 00529 00530 !compute G1n_n 00531 G1n_n = (0.d0,0.d0) 00532 do Jorb = 1,MOrb 00533 do Korb=1,MOrb 00534 00535 G1n_n = G1n_n + rho_jk(Jorb,Korb) * DCONJG(Psi(n,Jorb))*Psi(n,Korb) 00536 00537 end do 00538 end do 00539 G1n_n=G1n_n*1.d0/w**2 00540 outG1n_n=DBLE(G1n_n) 00541 00542 call get_ijk_from_m(n,NDX,NDY,i_n,j_n,k_n) 00543 00544 if (mode==1) then 00545 00546 write(12,6667) gridx(i_n)," ",gridy(j_n)," ",gridz(k_n)," ",outG1n_n," ",time 00547 00548 else if (mode==2) then 00549 00550 write(12,6667) mom_X(i_n)," ",mom_Y(j_n)," ",mom_Z(k_n)," ", outG1n_n," ",time 00551 00552 end if 00553 00554 if(MOD(n,NDX)==0) write(12,*)" " 00555 00556 end do 00557 00558 close(12) 00559 00560 00561 6667 format(3(F8.3,a1),2(F16.9,a1)) 00562 !---}}} 00563 end subroutine get_oneBodyRDMDiagonal 00564 00565 !-------------------------------------------------------------------------------------------- 00588 subroutine get_correlations(time,Psi,mode,dilation,rho_jk,rho_ijkl) 00589 USE moduleinputoutput,ONLY: get_doubleAsString 00590 USE modulegrid 00591 USE moduleinputvariables 00592 USE moduleparameters,ONLY: PI 00593 USE moduleallocatables,ONLY:gridx,gridy,gridz 00594 ! {{{ 00595 00596 implicit none 00597 integer :: noGridPoints 00598 integer,intent(in) :: mode,dilation 00599 real*8,intent(in) :: time 00600 real*8 :: w,dk,mom_X(NDX),mom_Y(NDY),mom_Z(NDZ) 00601 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(in) :: Psi 00602 complex*16,intent(in) :: rho_ijkl(MOrb,MOrb,MOrb,MOrb) 00603 complex*16,intent(in) :: rho_jk(MOrb,MOrb) 00604 character(len=2) :: aString 00605 character(len=200) :: fullName 00606 character(len=10) :: timeAsString 00607 character(len=10) :: MOrbAsString 00608 character(len=1) :: NparAsString1 00609 character(len=2) :: NparAsString2 00610 character(len=3) :: NparAsString3 00611 character(len=4) :: NparAsString4 00612 integer :: m,n,Iorb,Jorb,Korb,Lorb,i_m,j_m,k_m,i_n,j_n,k_n 00613 complex*16 :: G1m_m,G1m_n,G1n_n,G2m_n 00614 00615 complex*16 :: maximum 00616 real*8,parameter :: TOL=0.00001d0 00617 real*8 :: threshold 00618 real*8 :: outG1m_m,outReG1m_n,outImG1m_n,outG1n_n,outG2m_n 00619 00620 noGridPoints=NDX*NDY*NDZ 00621 00622 if((NDY>1).or. (NDZ>1)) then 00623 write(*,*)"only 1d implemented so far" 00624 STOP 00625 end if 00626 00627 00628 if(mode==1) then 00629 w=DSQRT((xf-xi)/NDX) 00630 else if (mode==2) then 00631 dk = 2.d0*PI/(xf-xi)/dilation 00632 w = DSQRT( dk ) 00633 00634 CALL get_momentumGrid(xi,xf,NDX,mom_X) !returns momentum grid in strictly ascending order 00635 CALL get_momentumGrid(yi,yf,NDY,mom_Y) 00636 CALL get_momentumGrid(zi,zf,NDZ,mom_Z) 00637 00638 mom_X=mom_X/dilation 00639 mom_Y=mom_Y/dilation 00640 mom_Z=mom_Z/dilation 00641 00642 else 00643 write(*,*) "wrong mode" 00644 end if 00645 00646 write (MOrbAsString, '(I1)') MOrb 00647 call get_doubleAsString(time,timeAsString) 00648 00649 if (mode==1) then 00650 aString = 'x-' 00651 else if (mode==2) then 00652 aString = 'k-' 00653 else 00654 write(*,*)"ERROR in get_correlations" 00655 STOP 00656 endif 00657 00658 if ((NPar>0) .and. (NPar.lt.10)) then 00659 write (NParAsString1, '(I1)') NPar 00660 fullName=timeasString//'N'//trim(NParAsString1)//'M'//trim(MOrbAsString)//aString//'correlations.dat' 00661 else if ((NPar.ge.10) .and. (NPar.lt.100)) then 00662 write (NParAsString2, '(I2)') NPar 00663 fullName=timeasString//'N'//trim(NParAsString2)//'M'//trim(MOrbAsString)//aString//'correlations.dat' 00664 else if ((NPar.ge.100) .and. (NPar.lt.1000)) then 00665 write (NParAsString3, '(I3)') NPar 00666 fullName=timeasString//'N'//trim(NParAsString3)//'M'//trim(MOrbAsString)//aString//'correlations.dat' 00667 else if ((NPar.ge.1000) .and. (NPar.lt.10000)) then 00668 write (NParAsString4, '(I4)') NPar 00669 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MOrbAsString)//aString//'correlations.dat' 00670 else 00671 NParAsString4='XXXX' 00672 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MOrbAsString)//aString//'correlations.dat' 00673 end if 00674 ! needs to be edited for other dvrs with nonconstant weights! 00675 call get_maximum(Psi,rho_jk,maximum) 00676 threshold = 0.01d0 * DBLE(maximum)/w**2 ! set the minimal density to be plotted 00677 ! threshold = 0.d0 00678 00679 open(unit=12,file=trim(fullName),form='formatted') 00680 do n=1,noGridPoints 00681 00682 !compute G1n_n 00683 G1n_n = (0.d0,0.d0) 00684 do Jorb = 1,MOrb 00685 do Korb=1,MOrb 00686 00687 G1n_n = G1n_n + rho_jk(Jorb,Korb) * DCONJG(Psi(n,Jorb))*Psi(n,Korb) 00688 00689 end do 00690 end do 00691 G1n_n=G1n_n*1.d0/w**2 00692 outG1n_n=DBLE(G1n_n) 00693 !--------------- 00694 do m=1,noGridPoints 00695 00696 !compute G1m_m 00697 G1m_m = (0.d0,0.d0) 00698 do Jorb = 1,MOrb 00699 do Korb=1,MOrb 00700 00701 G1m_m = G1m_m + rho_jk(Jorb,Korb) * DCONJG(Psi(m,Jorb))*Psi(m,Korb) 00702 00703 end do 00704 end do 00705 G1m_m=G1m_m*1.d0/w**2 00706 outG1m_m=DBLE(G1m_m) 00707 !---------------- 00708 00709 !compute G1m_n 00710 G1m_n = (0.d0,0.d0) 00711 do Jorb = 1,MOrb 00712 do Korb=1,MOrb 00713 00714 G1m_n = G1m_n + rho_jk(Jorb,Korb) * DCONJG(Psi(m,Jorb))*Psi(n,Korb) 00715 00716 end do 00717 end do 00718 G1m_n=G1m_n*1.d0/(w**2) 00719 outReG1m_n=DBLE(G1m_n) 00720 outImG1m_n=DIMAG(G1m_n) 00721 !---------------- 00722 00723 00724 !compute G2m_n 00725 G2m_n=(0.d0,0.d0) 00726 do Lorb=1,MOrb 00727 do Korb=1,MOrb 00728 do Jorb=1,MOrb 00729 do Iorb=1,MOrb 00730 ! 00731 G2m_n = G2m_n + rho_ijkl(Iorb,Jorb,Korb,Lorb)* & 00732 DCONJG(Psi(m,Iorb))*DCONJG(Psi(n,Jorb))*Psi(m,Korb)*Psi(n,Lorb) 00733 00734 end do 00735 end do 00736 end do 00737 end do 00738 G2m_n=G2m_n*1.d0/(w**4) 00739 outG2m_n=DBLE(G2m_n) 00740 !---------------- 00741 00742 IF((DABS(outG1m_m).LT. MAX(TOL,threshold)) .OR. & 00743 (DABS(outG1n_n).LT. MAX(TOL,threshold)) ) THEN 00744 outReG1m_n = -600.d0 00745 outImG1m_n = -600.d0 00746 outG2m_n = -1.d0 00747 END IF 00748 00749 call get_ijk_from_m(m,NDX,NDY,i_m,j_m,k_m) 00750 call get_ijk_from_m(n,NDX,NDY,i_n,j_n,k_n) 00751 00752 if (mode==1) then 00753 write(12,6666) gridx(i_m)," ",gridy(j_m)," ",gridz(k_m)," ",gridx(i_n)," ",gridy(j_n)," ",gridz(k_n)," ", & 00754 outG1m_m," ",outReG1m_n," ",outImG1m_n," ",outG1n_n," ",outG2m_n 00755 else if (mode==2) then 00756 00757 write(12,6666) mom_X(i_m)," ",mom_Y(j_m)," ",mom_Z(k_m)," ",mom_X(i_n)," ",mom_Y(j_n)," ",mom_Z(k_n)," ", & 00758 outG1m_m," ",outReG1m_n," ",outImG1m_n," ",outG1n_n," ",outG2m_n 00759 end if 00760 00761 if(MOD(m,NDX)==0) write(12,*)" " 00762 00763 end do 00764 00765 end do 00766 close(12) 00767 00768 00769 6666 format(6(F8.3,a1),6(F16.9,a1)) 00770 !---}}} 00771 end subroutine get_correlations 00772 00773 00774 !-------------------------------------------------------------- 00794 subroutine get_dilation(dilation,Psi,rho_jk) 00795 USE moduleinputvariables,ONLY: NDX,NDY,NDZ,MOrb 00796 ! {{{ 00797 implicit none 00798 integer,intent(inout) :: dilation 00799 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(in) :: Psi 00800 complex*16,DIMENSION(MOrb,MOrb),intent(in) :: rho_jk 00801 00802 integer :: m,iLower,iUpper 00803 integer :: Jorb,Korb 00804 complex*16 :: cmax,G1m_m 00805 real*8 :: G1max,THRESHOLD,G1(NDX*NDY*NDZ) 00806 real*8 :: imin,imax 00807 00808 if((NDY>1) .or. (NDZ>1)) then 00809 write(*,*)"only 1d implemented" 00810 STOP 00811 end if 00812 00813 THRESHOLD = 0.01d0 00814 dilation = 1 00815 G1max = 0.d0 00816 00817 call get_maximum(Psi,rho_jk,cmax) 00818 G1max = DBLE(cmax) 00819 00820 !---compute G1m_m 00821 00822 do m=1,NDX*NDY*NDZ 00823 00824 G1m_m = (0.d0,0.d0) 00825 do Jorb = 1,MOrb 00826 do Korb=1,MOrb 00827 00828 G1m_m = G1m_m + rho_jk(Jorb,Korb) * DCONJG(Psi(m,Jorb))*Psi(m,Korb) 00829 00830 end do 00831 end do 00832 00833 G1(m)=DBLE(G1m_m) 00834 00835 end do 00836 00837 !---------------- 00838 !----find boundaries 00839 00840 do m=1,NDX/2 00841 if(G1(m)/G1max > THRESHOLD) THEN 00842 iLower = m 00843 iUpper = NDX-m+1 00844 exit 00845 end if 00846 if(G1(NDX-m+1)/G1max > THRESHOLD) THEN 00847 iLower = m 00848 iUpper = NDX-m+1 00849 exit 00850 end if 00851 end do 00852 00853 imax = 1.d0*(iUpper-(NDX)/2 - 1) 00854 imin = 1.d0*(iLower-(NDX)/2 - 1) 00855 imax = MAX(DABS(imax),DABS(imin)) 00856 dilation = MAX(INT(NDX/2.d0 /(2.d0*imax)),1) 00857 00858 !!}}} 00859 end subroutine get_dilation 00860 00861 !------------------------------------------------------------------------------ 00872 subroutine pedestrian_FT(Psi,FTPsi,dilation) 00873 USE moduleinputvariables 00874 USE moduleparameters,ONLY: PI 00875 USE moduleallocatables,ONLY:gridx,gridy,gridz 00876 !---{{{ 00877 implicit none 00878 integer :: Iorb,m,n 00879 00880 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(in) :: Psi 00881 complex*16,DIMENSION(NDX*NDY*NDZ,MOrb),intent(out) :: FTPsi 00882 00883 real*8,dimension(NDX*NDY*NDZ) :: kvecs 00884 real*8 :: xall,xfirst,xlast,dk,dx,kmin,k,norm 00885 integer,intent(in) :: dilation 00886 00887 !---{{{ check for consistency 00888 if ((NDY>1) .or. (NDZ>1)) then 00889 write(*,*) "more than 1-D has not been implemented yet" 00890 STOP 00891 end if 00892 !---}}} 00893 00894 xfirst = xi 00895 xall = xf-xi 00896 dx = xall/NDX 00897 dk = 2.d0*PI/(xall*dilation) 00898 kmin = -(NDX/2) * dk 00899 00900 FTPsi=(0.d0,0.d0) 00901 00902 do Iorb=1,MOrb 00903 00904 00905 do m=1,NDX 00906 00907 k = kmin+(m-1)*dk 00908 00909 do n=1,NDX 00910 FTPsi(m,Iorb) = FTPsi(m,Iorb) + dx * Psi(n,Iorb)/DSQRT(dx*2.d0*PI) * CDEXP( (0.d0,-1.d0)*k*gridx(n) ) 00911 end do 00912 00913 end do 00914 00915 call normvxz(FTPsi(:,Iorb),norm,NDX) 00916 call xvixdzo(norm,FTPsi(:,Iorb),NDX) 00917 00918 end do 00919 00920 00921 !---}}} 00922 end subroutine pedestrian_FT 00923 00924 !-------------------------------------------------------------------------------------------- 00942 subroutine get_corr_slice(time,PSI,mode,rho_jk,rho_ijkl,& 00943 x1const,x1slice,y1const,y1slice,x2const,x2slice,y2const,y2slice) 00944 ! {{{ 00945 USE moduleinputoutput,ONLY: get_doubleAsString 00946 USE modulegrid 00947 use modulederivesystemvariables, ONLY: get_weight 00948 USE moduleinputvariables,ONLY: NDX,NDY,NDZ,morb,dims,NPAR,xi,xf,yi,yf,zi,zf 00949 USE moduleparameters,ONLY: PI 00950 USE moduleallocatables,ONLY:gridx,gridy,gridz 00951 implicit none 00952 integer :: NOGRIDPOINTS 00953 integer,intent(in) :: mode 00954 real*8,intent(in) :: time 00955 real*8,intent(in) :: x1slice,x2slice,y1slice,y2slice 00956 logical,intent(in) :: x1const,y1const,x2const,y2const 00957 real*8 :: w,dk,mom_X(NDX),mom_Y(NDY),mom_Z(NDZ) 00958 COMPLEX*16, DIMENSION(NDX*NDY*NDZ,Morb) :: PSI 00959 complex*16,intent(in) :: rho_ijkl(Morb,Morb,Morb,Morb) 00960 complex*16,intent(in) :: rho_jk(Morb,Morb) 00961 character(len=2) :: aString 00962 character(len=200) :: fullName 00963 character(len=10) :: timeAsString 00964 character(len=7) :: slice1AsString 00965 character(len=7) :: slice2AsString 00966 real*8 :: slice1,slice2 00967 integer :: slice1pt,slice2pt 00968 character(len=10) :: MorbAsString 00969 character(len=1) :: NparAsString1 00970 character(len=2) :: NparAsString2 00971 character(len=3) :: NparAsString3 00972 character(len=4) :: NparAsString4 00973 integer :: m,n,Iorb,Jorb,Korb,Lorb,i_m,j_m,k_m,i_n,j_n,k_n,k_k 00974 complex*16 :: G1m_m,G1m_n,G1n_n,G2m_n 00975 Real*8 :: onebdens(NDX*NDY) 00976 00977 complex*16 :: maximum 00978 real*8,parameter :: TOL=0.00001d0 00979 real*8 :: threshold 00980 real*8 :: outG1m_m,outReG1m_n,outImG1m_n,outG1n_n,outG2m_n 00981 00982 00983 if (NDY.eq.1) then 00984 write(*,*) 'The sliced correlations are implemented for >=2D systems only!!' 00985 return 00986 endif 00987 00988 00989 NOGRIDPOINTS=NDX*NDY*NDZ 00990 if(mode.eq.1) then 00991 CALL get_weight(xi,xf,yi,yf,zi,zf,NDX,NDY,NDZ,Dims,w) 00992 else if (mode.eq.2) then 00993 CALL get_momentumGrid(xi,xf,NDX,mom_X) !returns momentum grid in strictly ascending order 00994 CALL get_momentumGrid(yi,yf,NDY,mom_Y) 00995 CALL get_momentumGrid(zi,zf,NDZ,mom_Z) 00996 00997 dk = (mom_Y(2)-mom_y(1))*(mom_X(2)-mom_X(1)) 00998 w = DSQRT( dk ) 00999 else 01000 write(*,*) "wrong mode" 01001 end if 01002 01003 01004 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01005 ! some exceptions 01006 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01007 if (mode.eq.1) then 01008 01009 if ((x1const.eqv..true.).and.((x1slice.lt.xi).or.(x1slice.gt.xf))) then 01010 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01011 return 01012 endif 01013 if ((x2const.eqv..true.).and.((x2slice.lt.xi).or.(x2slice.gt.xf))) then 01014 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01015 return 01016 endif 01017 if ((y1const.eqv..true.).and.((y1slice.lt.yi).or.(y1slice.gt.yf))) then 01018 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01019 return 01020 endif 01021 if ((y2const.eqv..true.).and.((y2slice.lt.xi).or.(y2slice.gt.xf))) then 01022 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01023 return 01024 endif 01025 01026 else if (mode.eq.2) then 01027 01028 if ((x1const.eqv..true.).and.((x1slice.gt.mom_x(NDX)).or.(x1slice.lt.mom_x(1)))) then 01029 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01030 return 01031 endif 01032 if ((x2const.eqv..true.).and.((x2slice.gt.mom_x(ndx)).or.(x2slice.lt.mom_x(1)))) then 01033 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01034 return 01035 endif 01036 if ((y1const.eqv..true.).and.((y1slice.gt.mom_y(ndx)).or.(y1slice.lt.mom_y(1)))) then 01037 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01038 return 01039 endif 01040 if ((y2const.eqv..true.).and.((y2slice.gt.mom_y(ndx)).or.(y2slice.lt.mom_y(1)))) then 01041 write(*,*) 'The wanted x1 slice of the correlations is outside the grid!!!' 01042 return 01043 endif 01044 01045 endif 01046 m=0 01047 01048 01049 if (x1const.eqv..true.) then 01050 m=m+1 01051 call getsliceasstring(1,x1slice,slice1AsString) 01052 slice1=x1slice 01053 endif 01054 if (y1const.eqv..true.) then 01055 m=m+1 01056 if (m.eq.1) then 01057 call getsliceasstring(2,y1slice,slice1AsString) 01058 slice1=y1slice 01059 endif 01060 if (m.eq.2) then 01061 call getsliceasstring(2,y1slice,slice2AsString) 01062 slice2=y1slice 01063 endif 01064 endif 01065 if (x2const.eqv..true.) then 01066 m=m+1 01067 if (m.eq.1) then 01068 call getsliceasstring(3,x2slice,slice1AsString) 01069 slice1=x2slice 01070 endif 01071 if (m.eq.2) then 01072 call getsliceasstring(3,x2slice,slice2AsString) 01073 slice2=x2slice 01074 endif 01075 endif 01076 if (y2const.eqv..true.) then 01077 m=m+1 01078 if (m.eq.1) then 01079 call getsliceasstring(4,y2slice,slice1AsString) 01080 slice1=y2slice 01081 endif 01082 if (m.eq.2) then 01083 call getsliceasstring(4,y2slice,slice2AsString) 01084 slice2=y2slice 01085 endif 01086 endif 01087 01088 if (m.ne.2) then 01089 write(*,*) 'You must keep two of the four coordinates of the correlations constant!!!' 01090 return 01091 endif 01092 01093 01094 if (MORB.lt.10) then 01095 write (MorbAsString, '(I1)') Morb 01096 else if ((MORB.ge.10).and.(MORB.lt.100)) then 01097 write (MorbAsString, '(I2)') Morb 01098 else 01099 write(*,*) 'Bigger orbital number than 100 not implemented!' 01100 endif 01101 01102 call get_DoubleAsString(time,timeAsString) 01103 01104 if (mode.eq.1) then 01105 aString = 'x-' 01106 else if (mode.eq.2) then 01107 aString = 'k-' 01108 else 01109 write(*,*)"ERROR in get_correlations" 01110 stop 01111 endif 01112 ! assemble filename 01113 if ((NPar.gt.0) .and. (NPar.lt.10)) then 01114 write (NParAsString1, '(I1)') NPar 01115 fullName=timeasString//'N'//trim(NParAsString1)//'M'//trim(MorbAsString)//aString//& 01116 slice1Asstring//'-'//slice2Asstring//'-'//'correlations.dat' 01117 else if ((NPar.ge.10) .and. (NPar.lt.100)) then 01118 write (NParAsString2, '(I2)') NPar 01119 fullName=timeasString//'N'//trim(NParAsString2)//'M'//trim(MorbAsString)//aString//& 01120 slice1Asstring//'-'//slice2Asstring//'-'//'correlations.dat' 01121 else if ((NPar.ge.100) .and. (NPar.lt.1000)) then 01122 write (NParAsString3, '(I3)') NPar 01123 fullName=timeasString//'N'//trim(NParAsString3)//'M'//trim(MorbAsString)//aString//& 01124 slice1Asstring//'-'//slice2Asstring//'-'//'correlations.dat' 01125 else if ((NPar.ge.1000) .and. (NPar.lt.10000)) then 01126 write (NParAsString4, '(I4)') NPar 01127 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MorbAsString)//aString//& 01128 slice1Asstring//'-'//slice2Asstring//'-'//'correlations.dat' 01129 else 01130 NParAsString4='XXXX' 01131 fullName=timeasString//'N'//trim(NParAsString4)//'M'//trim(MorbAsString)//aString//& 01132 slice1Asstring//'-'//slice2Asstring//'-'//'correlations.dat' 01133 end if 01134 01135 open(unit=12,file=trim(fullName),form='formatted') 01136 do n=1,NOGRIDPOINTS 01137 !compute one body density onebdens 01138 G1n_n = (0.d0,0.d0) 01139 do Jorb = 1,Morb 01140 do Korb=1,Morb 01141 G1n_n = G1n_n + rho_jk(Jorb,Korb) * DCONJG(PSI(n,Jorb))*PSI(n,Korb) 01142 end do 01143 end do 01144 G1n_n=G1n_n*1.d0/w**2 01145 onebdens(n)=DBLE(G1n_n) 01146 !--------------- 01147 end do 01148 01149 ! determine which points are to be plotted 01150 do n=1,NOGRIDPOINTS 01151 call get_ijk_from_m(n,NDX,NDY,i_n,j_n,k_n) 01152 01153 ! x1slice 01154 if ((mode.eq.1).and.((abs(gridX(i_n)-x1slice).lt.gridX(2)-gridX(1)).and.(x1const.eqv..true.))& 01155 .or.((mode.eq.2).and.(abs(mom_x(i_n)-x1slice).lt.mom_X(2)-mom_X(1)).and.(x1const.eqv..true.))) then 01156 call get_ijk_from_m(n+1,NDX,NDY,i_m,j_m,k_m) 01157 if ((abs(gridX(i_n)-x1slice)).lt.(abs(gridX(i_m)-x1slice))) then 01158 slice1pt=i_n 01159 else 01160 slice1pt=i_m 01161 endif 01162 endif 01163 01164 ! y1slice 01165 if (((mode.eq.1).and.(abs(gridY(j_n)-y1slice).lt.gridY(2)-gridY(1)).and.(y1const.eqv..true.)).or.& 01166 ((mode.eq.2).and.(abs(mom_Y(j_n)-y1slice).lt.mom_Y(2)-mom_Y(1)).and.(y1const.eqv..true.))) then 01167 call get_ijk_from_m(n+1,NDX,NDY,i_m,j_m,k_m) 01168 if ((abs(gridY(j_n)-y1slice)).lt.(abs(gridY(j_m)-y1slice))) then 01169 if (x1const.eqv..false.) slice1pt=j_n 01170 if (x1const.eqv..true.) then 01171 slice2pt=j_n 01172 goto 8855 01173 endif 01174 else 01175 if (x1const.eqv..false.) slice1pt=j_m 01176 if (x1const.eqv..true.) then 01177 slice2pt=j_m 01178 goto 8855 01179 endif 01180 endif 01181 endif 01182 01183 ! x2slice 01184 if (((mode.eq.1).and.(abs(gridX(i_n)-x2slice).lt.gridx(2)-gridx(1)).and.(x2const.eqv..true.)).or.& 01185 ((mode.eq.2).and.(abs(mom_X(i_n)-x2slice).lt.mom_x(2)-mom_x(1)).and.(x2const.eqv..true.))) then 01186 call get_ijk_from_m(n+1,NDX,NDY,i_m,j_m,k_m) 01187 if ((abs(gridx(i_n)-x2slice)).lt.(abs(gridX(i_m)-x2slice))) then 01188 if ((x1const.eqv..false.).and.(y1const.eqv..false.)) slice1pt=i_n 01189 if ((x1const.eqv..true.).or.(y1const.eqv..true.)) then 01190 slice2pt=i_n 01191 goto 8855 01192 endif 01193 else 01194 if ((x1const.eqv..false.).and.(y1const.eqv..false.)) slice1pt=j_m 01195 if ((x1const.eqv..true.).or.(y1const.eqv..true.)) then 01196 slice2pt=j_m 01197 goto 8855 01198 endif 01199 endif 01200 endif 01201 01202 ! y2slice 01203 if (((mode.eq.1).and.(abs(gridY(j_n)-y2slice).lt.gridY(2)-gridy(1)).and.(y2const.eqv..true.)).or.& 01204 ((mode.eq.2).and.(abs(mom_Y(j_n)-y2slice).lt.mom_Y(2)-mom_Y(1)).and.(y2const.eqv..true.))) then 01205 call get_ijk_from_m(n+1,NDX,NDY,i_m,j_m,k_m) 01206 if ((abs(gridY(j_n)-y2slice)).lt.(abs(gridY(j_m)-y2slice))) then 01207 if ((x1const.eqv..true.).or.(y1const.eqv..true.).or.((x2const.eqv..true.))) then 01208 slice2pt=j_n 01209 goto 8855 01210 else 01211 write(*,*) 'Something went wrong with the assignment of the slices!!!' 01212 stop 01213 endif 01214 else 01215 if ((x1const.eqv..true.).or.(y1const.eqv..true.).or.((x2const.eqv..true.))) then 01216 slice2pt=j_m 01217 goto 8855 01218 else 01219 write(*,*) 'Something went wrong with the assignment of the slices!!!' 01220 stop 01221 endif 01222 endif 01223 endif 01224 end do 01225 8855 continue 01226 01227 01228 do n=1,NOGRIDPOINTS 01229 call get_ijk_from_m(n,NDX,NDY,i_n,j_n,k_n) 01230 ! x1const and y1const 01231 if (((x1const.eqv..true.).and.(i_n.eq.slice1pt)).and.& 01232 ((y1const.eqv..true.).and.(j_n.eq.slice2pt))) then 01233 do m=1,NOGRIDPOINTS 01234 call get_ijk_from_m(m,NDX,NDY,i_m,j_m,k_m) 01235 01236 !compute G1m_n 01237 G1m_n = (0.d0,0.d0) 01238 do Jorb = 1,Morb 01239 do Korb=1,Morb 01240 G1m_n = G1m_n + rho_jk(Jorb,Korb) * DCONJG(PSI(m,Jorb))*PSI(n,Korb) 01241 end do 01242 end do 01243 G1m_n=G1m_n*1.d0/(w**2) 01244 outReG1m_n=DBLE(G1m_n) 01245 outImG1m_n=DIMAG(G1m_n) 01246 !---------------- 01247 01248 01249 !compute G2m_n 01250 G2m_n=(0.d0,0.d0) 01251 do Lorb=1,Morb 01252 do Korb=1,Morb 01253 do Jorb=1,Morb 01254 do Iorb=1,Morb 01255 G2m_n = G2m_n + rho_ijkl(Iorb,Jorb,Korb,Lorb)* & 01256 DCONJG(PSI(m,Iorb))*DCONJG(PSI(n,Jorb))*PSI(m,Korb)*PSI(n,Lorb) 01257 01258 end do 01259 end do 01260 end do 01261 end do 01262 G2m_n=G2m_n*1.d0/(w*w*w*w) 01263 outG2m_n=DBLE(G2m_n) 01264 !---------------- 01265 01266 if (mode.eq.1) then 01267 write(12,6666) gridx(i_m),gridy(j_m),gridx(i_n),gridy(j_n), & 01268 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01269 else if (mode.eq.2) then 01270 write(12,6666) mom_X(i_m),mom_Y(j_m),mom_X(i_n),mom_Y(j_n), & 01271 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01272 end if 01273 if (mod(m,ndx).eq.0) write(12,*)" " 01274 01275 end do 01276 ! if x1const and x2const or y2const 01277 else if ((x1const.eqv..true.).and.((x2const.eqv..true.).or.(y2const.eqv..true.))& 01278 .and.(i_n.eq.slice1pt)) then 01279 do m=1,NOGRIDPOINTS 01280 01281 call get_ijk_from_m(m,NDX,NDY,i_m,j_m,k_m) 01282 ! if the point is in the slice 01283 if (((x2const.eqv..true.).and.(i_m.eq.slice2pt))& 01284 .or.((y2const.eqv..true.).and.(j_m.eq.slice2pt))) then 01285 !compute G1m_n 01286 G1m_n = (0.d0,0.d0) 01287 do Jorb = 1,Morb 01288 do Korb=1,Morb 01289 G1m_n = G1m_n + rho_jk(Jorb,Korb) * DCONJG(PSI(m,Jorb))*PSI(n,Korb) 01290 end do 01291 end do 01292 G1m_n=G1m_n*1.d0/(w**2) 01293 outReG1m_n=DBLE(G1m_n) 01294 outImG1m_n=DIMAG(G1m_n) 01295 !---------------- 01296 01297 01298 !compute G2m_n 01299 G2m_n=(0.d0,0.d0) 01300 do Lorb=1,Morb 01301 do Korb=1,Morb 01302 do Jorb=1,Morb 01303 do Iorb=1,Morb 01304 G2m_n = G2m_n + rho_ijkl(Iorb,Jorb,Korb,Lorb)* & 01305 DCONJG(PSI(m,Iorb))*DCONJG(PSI(n,Jorb))*PSI(m,Korb)*PSI(n,Lorb) 01306 01307 end do 01308 end do 01309 end do 01310 end do 01311 G2m_n=G2m_n*1.d0/(w*w*w*w) 01312 outG2m_n=DBLE(G2m_n) 01313 !---------------- 01314 01315 if (mode.eq.1) then 01316 write(12,6666) gridx(i_m),gridy(j_m),gridx(i_n),gridy(j_n), & 01317 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01318 else if (mode.eq.2) then 01319 write(12,6666) mom_X(i_m),mom_Y(j_m),mom_X(i_n),mom_Y(j_n), & 01320 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01321 end if 01322 01323 endif ! loops only executed if point is in the slice 01324 enddo 01325 write(12,*) " " 01326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01327 ! if x1const false, y1const true and x2/y2const true 01328 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01329 elseif ((x1const.eqv..false.).and.(y1const.eqv..true.).and.(j_n.eq.slice1pt)) then 01330 01331 do m=1,NOGRIDPOINTS 01332 call get_ijk_from_m(m,NDX,NDY,i_m,j_m,k_m) 01333 ! if the point is in the slice 01334 if (((x2const.eqv..true.).and.(i_m.eq.slice2pt))& 01335 .or.((y2const.eqv..true.).and.(j_m.eq.slice2pt))) then 01336 !compute G1m_n 01337 G1m_n = (0.d0,0.d0) 01338 do Jorb = 1,Morb 01339 do Korb=1,Morb 01340 G1m_n = G1m_n + rho_jk(Jorb,Korb) * DCONJG(PSI(m,Jorb))*PSI(n,Korb) 01341 end do 01342 end do 01343 G1m_n=G1m_n*1.d0/(w**2) 01344 outReG1m_n=DBLE(G1m_n) 01345 outImG1m_n=DIMAG(G1m_n) 01346 01347 !compute G2m_n 01348 G2m_n=(0.d0,0.d0) 01349 do Lorb=1,Morb 01350 do Korb=1,Morb 01351 do Jorb=1,Morb 01352 do Iorb=1,Morb 01353 G2m_n = G2m_n + rho_ijkl(Iorb,Jorb,Korb,Lorb)* & 01354 DCONJG(PSI(m,Iorb))*DCONJG(PSI(n,Jorb))*PSI(m,Korb)*PSI(n,Lorb) 01355 01356 end do 01357 end do 01358 end do 01359 end do 01360 G2m_n=G2m_n*1.d0/(w*w*w*w) 01361 outG2m_n=DBLE(G2m_n) 01362 !---------------- 01363 01364 if (mode.eq.1) then 01365 write(12,6666) gridx(i_m),gridy(j_m),gridx(i_n),gridy(j_n), & 01366 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01367 else if (mode.eq.2) then 01368 write(12,6666) mom_X(i_m),mom_Y(j_m),mom_X(i_n),mom_Y(j_n), & 01369 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01370 end if 01371 01372 endif ! loops only executed if point is in the slice 01373 enddo 01374 write(12,*)" " 01375 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01376 ! if x1const false, y1const false and x2const and y2const true 01377 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01378 elseif ((x1const.eqv..false.).and.(y1const.eqv..false.)) then 01379 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01380 do m=1,NOGRIDPOINTS 01381 call get_ijk_from_m(m,NDX,NDY,i_m,j_m,k_m) 01382 ! if the point is in the slice 01383 if (((x2const.eqv..true.).and.(i_m.eq.slice1pt))& 01384 .and.((y2const.eqv..true.).and.(j_m.eq.slice2pt))) then 01385 !compute G1m_n 01386 G1m_n = (0.d0,0.d0) 01387 do Jorb = 1,Morb 01388 do Korb=1,Morb 01389 G1m_n = G1m_n + rho_jk(Jorb,Korb) * DCONJG(PSI(m,Jorb))*PSI(n,Korb) 01390 end do 01391 end do 01392 G1m_n=G1m_n*1.d0/(w**2) 01393 outReG1m_n=DBLE(G1m_n) 01394 outImG1m_n=DIMAG(G1m_n) 01395 !---------------- 01396 01397 01398 !compute G2m_n 01399 G2m_n=(0.d0,0.d0) 01400 do Lorb=1,Morb 01401 do Korb=1,Morb 01402 do Jorb=1,Morb 01403 do Iorb=1,Morb 01404 G2m_n = G2m_n + rho_ijkl(Iorb,Jorb,Korb,Lorb)* & 01405 DCONJG(PSI(m,Iorb))*DCONJG(PSI(n,Jorb))*PSI(m,Korb)*PSI(n,Lorb) 01406 01407 end do 01408 end do 01409 end do 01410 end do 01411 G2m_n=G2m_n*1.d0/(w*w*w*w) 01412 outG2m_n=DBLE(G2m_n) 01413 !---------------- 01414 if (mode.eq.1) then 01415 write(12,6666) gridx(i_m),gridy(j_m),gridx(i_n),gridy(j_n), & 01416 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01417 else if (mode.eq.2) then 01418 write(12,6666) mom_X(i_m),mom_Y(j_m),mom_X(i_n),mom_Y(j_n), & 01419 onebdens(m),onebdens(n),outReG1m_n,outImG1m_n,outG2m_n 01420 end if 01421 01422 endif ! loops only executed if point is in the slice 01423 enddo 01424 write(12,*)" " 01425 endif !loop m only if n is in the slice 01426 end do 01427 01428 6666 FORMAT(9(E20.10)) 01429 01430 close(12) 01431 01432 !---}}} 01433 end subroutine get_corr_slice 01440 subroutine FT_SORTED(PSI,FTPSI) 01441 !---{{{ 01442 USE moduleinputvariables 01443 USE moduleparameters,ONLY: PI 01444 USE moduleallocatables,ONLY:gridx,gridy,gridz 01445 USE modulefft 01446 01447 implicit none 01448 01449 integer :: Iorb,m,n,i,j,k,INFO 01450 01451 complex*16,DIMENSION(NDX*NDY*NDZ,Morb),intent(in) :: PSI 01452 complex*16,DIMENSION(NDX*NDY*NDZ,Morb),intent(out) :: FTPSI 01453 complex*16,DIMENSION(NDX*NDY*NDZ) :: ORBTMP 01454 01455 real*8,dimension(NDX*NDY*NDZ) :: kvecs 01456 real*8 :: xall,xfirst,xlast,dk,dx,kmin,norm 01457 real*8 :: nrm 01458 01459 if (NDZ.gt.1) then 01460 write(*,*) "more than 2-D FFT sorting has not been implemented yet" 01461 return 01462 end if 01463 !########################################################## 01464 !########################################################## 01465 ! 1D -CASE 01466 if (dims.eq.1) then 01467 !########################################################## 01468 !########################################################## 01469 FTPSI=PSI 01470 do Iorb=1,Morb 01471 ORBTMP=(FTPSI(:,IORB)) 01472 CALL get_FFT1D(-1,NDX,ORBTMP,INFO) 01473 01474 if (mod(NDX,2).eq.0) then 01475 do m=(NDX/2+1),NDX 01476 FTPSI(m-(NDX/2),IORB)=ORBTMP(m) 01477 FTPSI(m,IORB)=ORBTMP(m-(NDX/2)) 01478 enddo 01479 else 01480 do m=(NDX/2+1),NDX 01481 FTPSI(m-(NDX/2),IORB)=ORBTMP(m) 01482 FTPSI(m,IORB)=ORBTMP(m-(NDX/2)) 01483 enddo 01484 endif 01485 01486 call normvxz(FTPSI(:,Iorb),nrm,NDX) 01487 call xvixdzo(nrm,FTPSI(:,Iorb),NDX) 01488 01489 end do 01490 01491 !########################################################## 01492 !########################################################## 01493 ! 1D -CASE 01494 endif 01495 !########################################################## 01496 !########################################################## 01497 01498 !########################################################## 01499 !########################################################## 01500 ! 2D -CASE 01501 if (DIMs.eq.2) then 01502 !########################################################## 01503 !########################################################## 01504 FTPSI=PSI 01505 do Iorb=1,Morb 01506 ORBTMP=(FTPSI(:,IORB)) 01507 01508 CALL get_FFT2D(-1,NDX,NDY,ORBTMP,INFO) 01509 call normvxz(ORBTMP,nrm,NDX*NDY) 01510 call xvixdzo(nrm,ORBTMP,NDX*NDY) 01511 01512 call sort_FFT_to_ascending_2d(ORBTMP,NDX,NDY) 01513 01514 FTPSI(:,IORB)=ORBTMP 01515 end do 01516 !########################################################## 01517 !########################################################## 01518 ! 2D -CASE 01519 endif 01520 !########################################################## 01521 !########################################################## 01522 01523 01524 !---}}} 01525 end subroutine FT_SORTED 01526 01533 subroutine sort_FFT_to_ascending_2d(FFT,NDX,NDY) 01534 IMPLICIT NONE 01535 integer NDX,NDY,i,j,k 01536 COMPLEX*16 FFT(NDX*NDY),ASC(NDX*NDY) 01537 01538 ASC=FFT 01539 01540 ! copy the first half of each vector in X direction to the corresponding last half 01541 01542 IF (MOD(NDX,2).eq.0) THEN 01543 do k=1,NDY 01544 do i=1,NDX 01545 if (i.le.((NDX)/2)) then 01546 FFT((k-1)*NDX+i)=ASC((k-1)*NDX+i+(NDX)/2) 01547 else if (i.gt.((NDX)/2)) then 01548 FFT((k-1)*NDX+i)=ASC((k-1)*NDX+i-((NDX)/2)) 01549 endif 01550 end do 01551 end do 01552 ENDIF 01553 01554 01555 ! copy the first half of each vector in X direction to the corresponding last half 01556 IF (MOD(NDX,2).ne.0) then 01557 do k=1,NDY 01558 do i=1,NDX 01559 if (i.le.((NDX+1)/2-1)) then 01560 FFT((k-1)*NDX+i)=ASC((k-1)*NDX+i+(NDX+1)/2) 01561 else if (i.gt.((NDX+1)/2)-1) then 01562 FFT((k-1)*NDX+i)=ASC((k-1)*NDX+i-((NDX+1)/2-1)) 01563 endif 01564 end do 01565 end do 01566 ENDIF 01567 01568 IF (MOD(NDY,2).eq.0) then 01569 ! copy y direction from 0 -> kmax to last half 01570 01571 do k=1,NDX*NDY 01572 if (k.le.(NDX*(NDY/2))) then 01573 ASC(k)=FFT(NDX*(NDY/2)+k) 01574 elseif (k.gt.NDX*((NDY)/2)) then 01575 ASC(k)=FFT(k-NDX*(NDY/2)) 01576 endif 01577 end do 01578 ENDIF 01579 IF (MOD(NDY,2).ne.0) then 01580 ! copy y direction from 0 -> kmax to last half 01581 do k=1,NDX*NDY 01582 if (k.le.(NDX*(NDY/2))) then 01583 ASC(k)=FFT(NDX*(NDY/2+1)+k) 01584 elseif (k.gt.NDX*(NDY/2)) then 01585 ASC(k)=FFT(k-NDX*((NDY/2))) 01586 endif 01587 end do 01588 ENDIF 01589 FFT=ASC 01590 end subroutine sort_FFT_to_ascending_2d 01591 01598 subroutine getsliceAsString(d,x,doubleString) 01599 ! {{{ 01600 ! 01601 IMPLICIT NONE 01602 integer :: d 01603 integer,parameter :: DBL=8 01604 real(kind=DBL) :: x 01605 CHARACTER(len=7) :: doubleString 01606 CHARACTER(len=2) :: slc 01607 01608 if (d.eq.1) then 01609 write(slc,'(A2)') 'c1' 01610 elseif (d.eq.2) then 01611 write(slc,'(A2)') 'c2' 01612 elseif (d.eq.3) then 01613 write(slc,'(A2)') 'c3' 01614 elseif (d.eq.4) then 01615 write(slc,'(A2)') 'c4' 01616 else 01617 write(*,*) 'This slice cannot be right!' 01618 stop 01619 endif 01620 01621 if (x .LE. 9999.99999d0) WRITE(doubleString,'(F5.0,A2)') x,slc 01622 if (x .LE. 999.999999d0) WRITE(doubleString,'(F5.1,A2)') x,slc 01623 if (x .LE. 99.9999999d0) WRITE(doubleString,'(F5.2,A2)') x,slc 01624 if (x .LE. 9.99999999d0) WRITE(doubleString,'(F5.3,A2)') x,slc 01625 01626 IF ((x .LT. 0.d0) .OR. (x .GE. 10000.0d0)) THEN 01627 WRITE(6,*)"get slice as string: x not implemented!" 01628 WRITE(doubleString,'(A7)') 'XXXXXXX' 01629 END IF 01630 01631 !}}} 01632 end subroutine getsliceasString 01633 01634 END MODULE modulecorrelationfunctions
1.7.3
