diff --git a/src/AOtoMO/AOtoMO_ERI_GHF.f90 b/src/AOtoMO/AOtoMO_ERI_GHF.f90 index ecc4b06..8a12916 100644 --- a/src/AOtoMO/AOtoMO_ERI_GHF.f90 +++ b/src/AOtoMO/AOtoMO_ERI_GHF.f90 @@ -36,63 +36,4 @@ subroutine AOtoMO_ERI_GHF(nBas,nBas2,c1,c2,ERI_AO,ERI_MO) call dgemm('T','N',nBas2**3,nBas2,nBas,1d0,scr,nBas,c2(1,1),nBas,0d0,ERI_MO,nBas2**3) -! Four-index transform via semi-direct O(N^5) algorithm - -! scr(:,:,:,:) = 0d0 - -! do l=1,nBas2 -! do si=1,nBas -! do la=1,nBas -! do nu=1,nBas -! do mu=1,nBas -! scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO(mu,nu,la,si)*c2(si,l) -! enddo -! enddo -! enddo -! enddo -! enddo - -! ERI_MO(:,:,:,:) = 0d0 - -! do l=1,nBas2 -! do la=1,nBas -! do nu=1,nBas -! do i=1,nBas2 -! do mu=1,nBas -! ERI_MO(i,nu,la,l) = ERI_MO(i,nu,la,l) + c1(mu,i)*scr(mu,nu,la,l) -! enddo -! enddo -! enddo -! enddo -! enddo - - -! scr(:,:,:,:) = 0d0 - -! do l=1,nBas2 -! do k=1,nBas2 -! do la=1,nBas -! do nu=1,nBas -! do i=1,nBas2 -! scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO(i,nu,la,l)*c1(la,k) -! enddo -! enddo -! enddo -! enddo -! enddo - -! ERI_MO(:,:,:,:) = 0d0 - -! do l=1,nBas2 -! do k=1,nBas2 -! do j=1,nBas2 -! do i=1,nBas2 -! do nu=1,nBas -! ERI_MO(i,j,k,l) = ERI_MO(i,j,k,l) + c2(nu,j)*scr(i,nu,k,l) -! enddo -! enddo -! enddo -! enddo -! enddo - end subroutine diff --git a/src/AOtoMO/AOtoMO_oooa.f90 b/src/AOtoMO/AOtoMO_oooa.f90 index ce1a559..c8b60ea 100644 --- a/src/AOtoMO/AOtoMO_oooa.f90 +++ b/src/AOtoMO/AOtoMO_oooa.f90 @@ -35,11 +35,11 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) do si=1,nBas do x=1,nA scr1(mu,nu,la,x) = scr1(mu,nu,la,x) + O(mu,nu,la,si)*cA(si,x) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr2 = 0d0 do mu=1,nBas @@ -48,11 +48,11 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) do i=1,nO do x=1,nA scr2(i,nu,la,x) = scr2(i,nu,la,x) + cO(mu,i)*scr1(mu,nu,la,x) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr1 = 0d0 do nu=1,nBas @@ -61,11 +61,11 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) do k=1,nO do x=1,nA scr1(i,nu,k,x) = scr1(i,nu,k,x) + scr2(i,nu,la,x)*cO(la,k) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ooOoa = 0d0 do nu=1,nBas @@ -74,11 +74,11 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) do k=1,nO do x=1,nA ooOoa(i,j,k,x) = ooOoa(i,j,k,x) + cO(nu,j)*scr1(i,nu,k,x) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do deallocate(scr1,scr2) diff --git a/src/AOtoMO/AOtoMO_oooo.f90 b/src/AOtoMO/AOtoMO_oooo.f90 index e177290..75286bd 100644 --- a/src/AOtoMO/AOtoMO_oooo.f90 +++ b/src/AOtoMO/AOtoMO_oooo.f90 @@ -35,11 +35,11 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) do si=1,nBas do l=1,nO scr1(mu,nu,la,l) = scr1(mu,nu,la,l) + O(mu,nu,la,si)*cO(si,l) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr2 = 0d0 do mu=1,nBas @@ -48,11 +48,11 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) do i=1,nO do l=1,nO scr2(i,nu,la,l) = scr2(i,nu,la,l) + cO(mu,i)*scr1(mu,nu,la,l) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr1 = 0d0 do nu=1,nBas @@ -61,11 +61,11 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) do k=1,nO do l=1,nO scr1(i,nu,k,l) = scr1(i,nu,k,l) + scr2(i,nu,la,l)*cO(la,k) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ooOoo = 0d0 do nu=1,nBas @@ -74,11 +74,11 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) do k=1,nO do l=1,nO ooOoo(i,j,k,l) = ooOoo(i,j,k,l) + cO(nu,j)*scr1(i,nu,k,l) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do deallocate(scr1,scr2) diff --git a/src/AOtoMO/AOtoMO_oovv.f90 b/src/AOtoMO/AOtoMO_oovv.f90 index e82a19f..ad5e5d3 100644 --- a/src/AOtoMO/AOtoMO_oovv.f90 +++ b/src/AOtoMO/AOtoMO_oovv.f90 @@ -29,11 +29,11 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) do si=1,nBas do b=1,nV scr1(mu,nu,la,b) = scr1(mu,nu,la,b) + O(mu,nu,la,si)*cV(si,b) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr2 = 0d0 do mu=1,nBas @@ -42,11 +42,11 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) do i=1,nO do b=1,nV scr2(i,nu,la,b) = scr2(i,nu,la,b) + cO(mu,i)*scr1(mu,nu,la,b) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do scr1 = 0d0 do nu=1,nBas @@ -55,11 +55,11 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) do a=1,nV do b=1,nV scr1(i,nu,a,b) = scr1(i,nu,a,b) + scr2(i,nu,la,b)*cV(la,a) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ooOvv = 0d0 do nu=1,nBas @@ -68,10 +68,10 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) do a=1,nV do b=1,nV ooOvv(i,j,a,b) = ooOvv(i,j,a,b) + cO(nu,j)*scr1(i,nu,a,b) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do end subroutine diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 87ef787..e53efe4 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -26,9 +26,9 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H) do la=1,nBas do si=1,nBas H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 6568030..1a29038 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -25,9 +25,9 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K) do la=1,nBas do mu=1,nBas K(mu,nu) = K(mu,nu) - P(la,si)*ERI(mu,la,si,nu) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/CC/CCD.f90 b/src/CC/CCD.f90 index c44e00e..120d369 100644 --- a/src/CC/CCD.f90 +++ b/src/CC/CCD.f90 @@ -213,7 +213,7 @@ subroutine CCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -231,7 +231,7 @@ subroutine CCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/CCD_correlation_energy.f90 b/src/CC/CCD_correlation_energy.f90 index deb272c..6aad3fa 100644 --- a/src/CC/CCD_correlation_energy.f90 +++ b/src/CC/CCD_correlation_energy.f90 @@ -26,10 +26,10 @@ subroutine CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD) EcCCD = EcCCD + OOVV(i,j,a,b)*t2(i,j,a,b) - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcCCD = 0.25d0*EcCCD diff --git a/src/CC/CCGW.f90 b/src/CC/CCGW.f90 index 80fb8fd..b99644d 100644 --- a/src/CC/CCGW.f90 +++ b/src/CC/CCGW.f90 @@ -316,7 +316,7 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e) write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',eGW(nO)*HaToeV,'|',eGW(nO+1)*HaToeV,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -334,7 +334,7 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e) stop - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,*)' CCGW calculation ' @@ -346,7 +346,7 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',e(p)*HaToeV,'|',(eGW(p)-e(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' end subroutine diff --git a/src/CC/CCSD_correlation_energy.f90 b/src/CC/CCSD_correlation_energy.f90 index ace59a4..62c4f15 100644 --- a/src/CC/CCSD_correlation_energy.f90 +++ b/src/CC/CCSD_correlation_energy.f90 @@ -26,10 +26,10 @@ subroutine CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCCSD) EcCCSD = EcCCSD + OOVV(i,j,a,b)*tau(i,j,a,b) - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcCCSD = 0.5d0*EcCCSD diff --git a/src/CC/DCD.f90 b/src/CC/DCD.f90 index 91344e8..8fc9272 100644 --- a/src/CC/DCD.f90 +++ b/src/CC/DCD.f90 @@ -133,10 +133,10 @@ subroutine DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) EcCC = EcCC + (2d0*OOVV(i,j,a,b) - OOVV(i,j,b,a))*t(i,j,a,b) - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Dump results @@ -156,10 +156,10 @@ subroutine DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) do a=1,nV-nR do b=1,nV-nR tt(i,j,a,b) = 2d0*t(i,j,a,b) - t(i,j,b,a) - enddo - enddo - enddo - enddo + end do + end do + end do + end do xV(:,:) = 0d0 do a=1,nV-nR @@ -262,7 +262,7 @@ subroutine DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) if(abs(rcond) < 1d-15) n_diis = 0 - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -280,7 +280,7 @@ subroutine DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) stop - endif + end if ! Testing zone diff --git a/src/CC/GCCD.f90 b/src/CC/GCCD.f90 index 2f5869a..9d43f70 100644 --- a/src/CC/GCCD.f90 +++ b/src/CC/GCCD.f90 @@ -190,7 +190,7 @@ subroutine GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -208,7 +208,7 @@ subroutine GCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/MP3_correlation_energy.f90 b/src/CC/MP3_correlation_energy.f90 index 693013b..5bdc6ea 100644 --- a/src/CC/MP3_correlation_energy.f90 +++ b/src/CC/MP3_correlation_energy.f90 @@ -28,10 +28,10 @@ subroutine MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3) EcMP3 = EcMP3 + OOVV(i,j,a,b)*(t2(i,j,a,b) + v(i,j,a,b)/delta_OOVV(i,j,a,b)) - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcMP3 = 0.25d0*EcMP3 diff --git a/src/CC/antisymmetrize_ERI.f90 b/src/CC/antisymmetrize_ERI.f90 index f8e74cc..76616b1 100644 --- a/src/CC/antisymmetrize_ERI.f90 +++ b/src/CC/antisymmetrize_ERI.f90 @@ -24,10 +24,10 @@ subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI) do k=1,nBas do l=1,nBas db_ERI(i,j,k,l) = 2d0*ERI(i,j,k,l) - ERI(i,j,l,k) - enddo - enddo - enddo - enddo + end do + end do + end do + end do elseif(ispin == 2) then @@ -36,11 +36,11 @@ subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI) do k=1,nBas do l=1,nBas db_ERI(i,j,k,l) = ERI(i,j,k,l) - ERI(i,j,l,k) - enddo - enddo - enddo - enddo + end do + end do + end do + end do - endif + end if end subroutine diff --git a/src/CC/crCCD.f90 b/src/CC/crCCD.f90 index ffc1881..2f42f46 100644 --- a/src/CC/crCCD.f90 +++ b/src/CC/crCCD.f90 @@ -175,7 +175,7 @@ subroutine crCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,EN write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -193,7 +193,7 @@ subroutine crCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,EN stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/crGCCD.f90 b/src/CC/crGCCD.f90 index b1058aa..66bc58a 100644 --- a/src/CC/crGCCD.f90 +++ b/src/CC/crGCCD.f90 @@ -151,7 +151,7 @@ subroutine crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -169,7 +169,7 @@ subroutine crGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/drCCD.f90 b/src/CC/drCCD.f90 index 0a46948..6374c6f 100644 --- a/src/CC/drCCD.f90 +++ b/src/CC/drCCD.f90 @@ -168,7 +168,7 @@ subroutine drCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,EN write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -186,7 +186,7 @@ subroutine drCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,EN stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/drGCCD.f90 b/src/CC/drGCCD.f90 index 0692d63..39313fd 100644 --- a/src/CC/drGCCD.f90 +++ b/src/CC/drGCCD.f90 @@ -143,7 +143,7 @@ subroutine drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -161,7 +161,7 @@ subroutine drGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/form_EOM_one_body.f90 b/src/CC/form_EOM_one_body.f90 index 77973d6..99a579c 100644 --- a/src/CC/form_EOM_one_body.f90 +++ b/src/CC/form_EOM_one_body.f90 @@ -72,8 +72,8 @@ subroutine form_EOM_one_body(nO,nV,foo,fov,fvv,t1,t2,OOOV,OOVV,OVVV,cFoo,cFov,cF end do end do - enddo - enddo + end do + end do ! OO block diff --git a/src/CC/form_FB.f90 b/src/CC/form_FB.f90 index 99cec79..9dcc922 100644 --- a/src/CC/form_FB.f90 +++ b/src/CC/form_FB.f90 @@ -34,11 +34,11 @@ subroutine form_FB(nC,nO,nV,nR,foo,fvv,fov,OOOV,OOVV,OVVV,t,FooB,FvvB,FovB) do a=1,nV-nR do b=1,nV-nR FooB(i,k) = FooB(i,k) + 0.5d0*OOVV(k,j,a,b)*t(i,j,a,b) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ! Virtual-virtual block @@ -49,11 +49,11 @@ subroutine form_FB(nC,nO,nV,nR,foo,fvv,fov,OOOV,OOVV,OVVV,t,FooB,FvvB,FovB) do j=1,nO-nC do c=1,nV-nR FvvB(a,c) = FvvB(a,c) - 0.5d0*OOVV(i,j,c,b)*t(i,j,a,b) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ! Occupied-virtual block @@ -64,26 +64,26 @@ subroutine form_FB(nC,nO,nV,nR,foo,fvv,fov,OOOV,OOVV,OVVV,t,FooB,FvvB,FovB) do j=1,nO-nC do b=1,nV-nR FovB(i,a) = FovB(i,a) - fov(j,b)*t(i,j,a,b) - enddo - enddo + end do + end do do j=1,nO-nC do b=1,nV-nR do c=1,nV-nR FovB(i,a) = FovB(i,a) + 0.5d0*OVVV(j,a,b,c)*t(i,j,b,c) - enddo - enddo - enddo + end do + end do + end do do j=1,nO-nC do k=1,nO-nC do b=1,nV-nR FovB(i,a) = FovB(i,a) - 0.5d0*OOOV(j,k,i,b)*t(j,k,a,b) - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do end subroutine diff --git a/src/CC/form_X.f90 b/src/CC/form_X.f90 index 8aa82c2..af7c70a 100644 --- a/src/CC/form_X.f90 +++ b/src/CC/form_X.f90 @@ -38,12 +38,12 @@ subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) do c=1,nV-nR do d=1,nV-nR X1(k,l,i,j) = X1(k,l,i,j) + OOVV(k,l,c,d)*t2(i,j,c,d) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do ! Build X2 @@ -53,11 +53,11 @@ subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) do l=1,nO-nC do d=1,nV-nR X2(b,c) = X2(b,c) + OOVV(k,l,c,d)*t2(k,l,b,d) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ! Build X3 @@ -67,11 +67,11 @@ subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) do c=1,nV-nR do d=1,nV-nR X3(k,j) = X3(k,j) + OOVV(k,l,c,d)*t2(j,l,c,d) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do ! Build X4 @@ -82,11 +82,11 @@ subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) do k=1,nO-nC do c=1,nV-nR X4(i,l,a,d) = X4(i,l,a,d) + OOVV(k,l,c,d)*t2(i,k,a,c) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_delta_OOOVVV.f90 b/src/CC/form_delta_OOOVVV.f90 index 7069ab3..e5acc75 100644 --- a/src/CC/form_delta_OOOVVV.f90 +++ b/src/CC/form_delta_OOOVVV.f90 @@ -27,11 +27,11 @@ subroutine form_delta_OOOVVV(nC,nO,nV,nR,eO,eV,delta) delta(i,j,k,a,b,c) = eV(a) + eV(b) + eV(c) - eO(i) - eO(j) - eO(k) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_delta_OOVV.f90 b/src/CC/form_delta_OOVV.f90 index 8b41e89..00647ea 100644 --- a/src/CC/form_delta_OOVV.f90 +++ b/src/CC/form_delta_OOVV.f90 @@ -25,9 +25,9 @@ subroutine form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta) delta(i,j,a,b) = eV(a) + eV(b) - eO(i) - eO(j) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_delta_OV.f90 b/src/CC/form_delta_OV.f90 index e819511..012eb63 100644 --- a/src/CC/form_delta_OV.f90 +++ b/src/CC/form_delta_OV.f90 @@ -24,7 +24,7 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta) do i=1,nO-nC do a=1,nV-nR delta(i,a) = eV(a) - eO(i) - enddo - enddo + end do + end do end subroutine diff --git a/src/CC/form_tau.f90 b/src/CC/form_tau.f90 index 9c6207b..bc1714d 100644 --- a/src/CC/form_tau.f90 +++ b/src/CC/form_tau.f90 @@ -26,9 +26,9 @@ subroutine form_tau(nC,nO,nV,nR,t1,t2,tau) tau(i,j,a,b) = 0.5d0*t2(i,j,a,b) + t1(i,a)*t1(j,b) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_tau_nc.f90 b/src/CC/form_tau_nc.f90 index 4e910c1..c57e42b 100644 --- a/src/CC/form_tau_nc.f90 +++ b/src/CC/form_tau_nc.f90 @@ -26,9 +26,9 @@ subroutine form_tau_nc(nO,nV,t1,t2,tau) tau(i,j,a,b) = t2(i,j,a,b) + t1(i,a)*t1(j,b) - t1(i,b)*t1(j,a) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_taus_nc.f90 b/src/CC/form_taus_nc.f90 index c49b2d5..d590431 100644 --- a/src/CC/form_taus_nc.f90 +++ b/src/CC/form_taus_nc.f90 @@ -26,9 +26,9 @@ subroutine form_taus_nc(nO,nV,t1,t2,taus) taus(i,j,a,b) = t2(i,j,a,b) + 0.5d0*(t1(i,a)*t1(j,b) - t1(i,b)*t1(j,a)) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_u.f90 b/src/CC/form_u.f90 index ef612c2..d481d76 100644 --- a/src/CC/form_u.f90 +++ b/src/CC/form_u.f90 @@ -30,12 +30,12 @@ subroutine form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) do c=1,nV-nR do d=1,nV-nR u(i,j,a,b) = u(i,j,a,b) + 0.5d0*VVVV(a,b,c,d)*t2(i,j,c,d) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do i=1,nO-nC do j=1,nO-nC @@ -44,12 +44,12 @@ subroutine form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) do a=1,nV-nR do b=1,nV-nR u(i,j,a,b) = u(i,j,a,b) + 0.5d0*OOOO(k,l,i,j)*t2(k,l,a,b) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do i=1,nO-nC do j=1,nO-nC @@ -61,11 +61,11 @@ subroutine form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) + OVOV(k,a,j,c)*t2(i,k,b,c) & - OVOV(k,a,i,c)*t2(j,k,b,c) & + OVOV(k,b,i,c)*t2(j,k,a,c) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do end subroutine diff --git a/src/CC/form_v.f90 b/src/CC/form_v.f90 index 52a26b8..9379ae1 100644 --- a/src/CC/form_v.f90 +++ b/src/CC/form_v.f90 @@ -31,12 +31,12 @@ subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) do k=1,nO-nC do l=1,nO-nC v(i,j,a,b) = v(i,j,a,b) + 0.25d0*X1(k,l,i,j)*t2(k,l,a,b) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do i=1,nO-nC do j=1,nO-nC @@ -44,11 +44,11 @@ subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) do b=1,nV-nR do c=1,nV-nR v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X2(b,c)*t2(i,j,a,c) + X2(a,c)*t2(i,j,c,b)) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do do i=1,nO-nC do j=1,nO-nC @@ -56,11 +56,11 @@ subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) do b=1,nV-nR do k=1,nO-nC v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X3(k,j)*t2(i,k,a,b) + X3(k,i)*t2(k,j,a,b)) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do do i=1,nO-nC do j=1,nO-nC @@ -69,11 +69,11 @@ subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) do k=1,nO-nC do c=1,nV-nR v(i,j,a,b) = v(i,j,a,b) + (X4(i,k,a,c)*t2(j,k,b,c) + X4(i,k,b,c)*t2(k,j,a,c)) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do end subroutine diff --git a/src/CC/lCCD.f90 b/src/CC/lCCD.f90 index 440933f..7886c5e 100644 --- a/src/CC/lCCD.f90 +++ b/src/CC/lCCD.f90 @@ -188,7 +188,7 @@ subroutine lCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -206,7 +206,7 @@ subroutine lCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/lGCCD.f90 b/src/CC/lGCCD.f90 index 3c4a3d3..ea70196 100644 --- a/src/CC/lGCCD.f90 +++ b/src/CC/lGCCD.f90 @@ -164,7 +164,7 @@ subroutine lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eH write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -182,7 +182,7 @@ subroutine lGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eH stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 916d9af..5306944 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -195,7 +195,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -212,7 +212,7 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF stop - endif + end if if(dotest) then diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index bbd6444..59d3bc1 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -178,7 +178,7 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -196,7 +196,7 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CC/rGCCD.f90 b/src/CC/rGCCD.f90 index 6aea413..826a9ed 100644 --- a/src/CC/rGCCD.f90 +++ b/src/CC/rGCCD.f90 @@ -154,7 +154,7 @@ subroutine rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eH write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -172,7 +172,7 @@ subroutine rGCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eH stop - endif + end if write(*,*) write(*,*)'----------------------------------------------------' diff --git a/src/CI/CID.f90 b/src/CI/CID.f90 index ef8d65e..22399f0 100644 --- a/src/CI/CID.f90 +++ b/src/CI/CID.f90 @@ -215,6 +215,6 @@ subroutine CID(dotest,singlet,triplet,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0) print*,'Singlet CID transition vectors' call matout(nH,nH,H) write(*,*) - endif + end if end subroutine diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 38bce1d..35401fd 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -303,6 +303,6 @@ subroutine CISD(dotest,singlet,triplet,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0 print*,'Singlet CISD transition vectors' call matout(nH,nH,H) write(*,*) - endif + end if end subroutine diff --git a/src/CI/CIS_D.f90 b/src/CI/CIS_D.f90 index 629d0e1..bc7c4e2 100644 --- a/src/CI/CIS_D.f90 +++ b/src/CI/CIS_D.f90 @@ -257,11 +257,11 @@ subroutine CIS_D(ispin,nBasin,nCin,nOin,nVin,nRin,nSin,maxS,eHF,ERI,w,X) wD = wD - 0.25d0*u(i,j,a,b)**2/(delta(i,j,a,b) - w(m)) - enddo - enddo + end do + end do wD = wD + r(i,a)*v(i,a) - enddo - enddo + end do + end do wD = 0.5d0*wD ! Flush results diff --git a/src/CI/RCIS.f90 b/src/CI/RCIS.f90 index 80a7614..84489e3 100644 --- a/src/CI/RCIS.f90 +++ b/src/CI/RCIS.f90 @@ -54,7 +54,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in print*,'CIS matrix (singlet state)' call matout(nS,nS,A) write(*,*) - endif + end if call diagonalize_matrix(nS,A,Om) call print_excitation_energies('CIS@RHF','singlet',nS,Om) @@ -64,7 +64,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in print*,'Singlet CIS transition vectors' call matout(nS,nS,A) write(*,*) - endif + end if ! Compute CIS(D) correction @@ -79,7 +79,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in end if - endif + end if if(triplet) then @@ -90,7 +90,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in print*,'CIS matrix (triplet state)' call matout(nS,nS,A) write(*,*) - endif + end if call diagonalize_matrix(nS,A,Om) call print_excitation_energies('CIS@RHF','triplet',nS,Om) @@ -100,7 +100,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in print*,'Triplet CIS transition vectors' call matout(nS,nS,A) write(*,*) - endif + end if ! Compute CIS(D) correction @@ -115,6 +115,6 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in end if - endif + end if end subroutine diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index 735bb82..0b0234d 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -76,7 +76,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI print*,'CIS matrix (spin-conserved transitions)' call matout(nS_sc,nS_sc,A_sc) write(*,*) - endif + end if call diagonalize_matrix(nS_sc,A_sc,Om_sc) A_sc(:,:) = transpose(A_sc) @@ -88,7 +88,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI print*,'Spin-conserved CIS transition vectors' call matout(nS_sc,nS_sc,A_sc) write(*,*) - endif + end if ! Testing zone @@ -100,7 +100,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI deallocate(A_sc,Om_sc) - endif + end if !-----------------------! ! Spin-flip transitions ! @@ -124,7 +124,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI print*,'CIS matrix (spin-conserved transitions)' call matout(nS_sf,nS_sf,A_sf) write(*,*) - endif + end if call diagonalize_matrix(nS_sf,A_sf,Om_sf) A_sf(:,:) = transpose(A_sf) @@ -136,7 +136,7 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI print*,'Spin-flip CIS transition vectors' call matout(nS_sf,nS_sf,A_sf) write(*,*) - endif + end if ! Testing zone @@ -148,6 +148,6 @@ subroutine UCIS(dotest,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI deallocate(A_sf,Om_sf) - endif + end if end subroutine diff --git a/src/GF/RG0F3.f90 b/src/GF/RG0F3.f90 index 4995415..6bc6c7b 100644 --- a/src/GF/RG0F3.f90 +++ b/src/GF/RG0F3.f90 @@ -54,12 +54,12 @@ App(p,1) = App(p,1) & - (2d0*V(p,k,p,j) - V(p,k,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,k,i)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -74,12 +74,12 @@ App(p,2) = App(p,2) & + (2d0*V(p,c,p,b) - V(p,c,b,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(j,i,c,a)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -94,12 +94,12 @@ App(p,3) = App(p,3) & + (2d0*V(p,c,p,j) - V(p,c,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,c,i)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do App(:,4) = App(:,3) @@ -116,12 +116,12 @@ App(p,5) = App(p,5) & - (2d0*V(p,b,p,k) - V(p,b,k,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(i,j,k,a)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do App(:,6) = App(:,5) @@ -143,10 +143,10 @@ Bpp(p,1) = Bpp(p,1) & + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -158,10 +158,10 @@ Bpp(p,2) = Bpp(p,2) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Total second-order Green function @@ -184,12 +184,12 @@ Cpp(p,1) = Cpp(p,1) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,c,d)*V(p,i,c,d)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -204,12 +204,12 @@ Cpp(p,2) = Cpp(p,2) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,j,k)*V(p,i,j,k)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Cpp(:,3) = Cpp(:,2) @@ -225,12 +225,12 @@ Cpp(p,4) = Cpp(p,4) & + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(i,j,b,c)*V(p,a,b,c)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Cpp(:,5) = Cpp(:,4) @@ -246,12 +246,12 @@ Cpp(p,6) = Cpp(p,6) & - (2d0*V(p,a,k,l) - V(p,a,l,k))*V(k,l,i,j)*V(p,a,i,j)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do ! Frequency-dependent third-order contribution: "D" terms @@ -275,12 +275,12 @@ + V(p,i,b,a)*(V(a,j,i,c)*(4d0*V(p,j,b,c) - 2d0*V(p,j,c,b)) & + V(a,j,c,i)*( V(p,j,c,b) - 2d0*V(p,j,b,c)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -300,12 +300,12 @@ + V(p,i,a,c)*(V(a,b,i,j)*( V(p,b,j,c) - 2d0*V(p,b,c,j)) & + V(a,b,j,i)*( V(p,b,c,j) - 2d0*V(p,b,j,c)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Dpp(:,3) = Dpp(:,2) @@ -327,12 +327,12 @@ + V(p,a,j,k)*(V(j,i,a,b)*( V(p,i,b,k) - 2d0*V(p,i,k,b)) & + V(j,i,b,a)*( V(p,i,k,b) - 2d0*V(p,i,b,k)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Dpp(:,5) = Dpp(:,4) @@ -354,12 +354,12 @@ - V(p,a,i,k)*(V(i,b,a,j)*( V(p,b,j,k) - 2d0*V(p,b,k,j)) & + V(i,b,j,a)*( V(p,b,k,j) - 2d0*V(p,b,j,k)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do ! Compute renormalization factor (if required) @@ -421,7 +421,7 @@ Sig3(:) = Z(:)*Sig3(:) - endif + end if ! Total third-order Green function diff --git a/src/GF/UGF2_reg_self_energy.f90 b/src/GF/UGF2_reg_self_energy.f90 index d9374f9..c92939f 100644 --- a/src/GF/UGF2_reg_self_energy.f90 +++ b/src/GF/UGF2_reg_self_energy.f90 @@ -69,9 +69,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -87,9 +87,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: aa @@ -106,9 +106,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -124,12 +124,12 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do !------------------! ! Spin-down sector ! @@ -153,9 +153,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -171,9 +171,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: bb @@ -190,9 +190,9 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -208,12 +208,12 @@ subroutine UGF2_reg_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eG SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do Z(:,:) = 1d0/(1d0 - Z(:,:)) diff --git a/src/GF/UGF2_reg_self_energy_diag.f90 b/src/GF/UGF2_reg_self_energy_diag.f90 index 0cc0497..7072b06 100644 --- a/src/GF/UGF2_reg_self_energy_diag.f90 +++ b/src/GF/UGF2_reg_self_energy_diag.f90 @@ -67,9 +67,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -85,9 +85,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: aa @@ -104,9 +104,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -122,11 +122,11 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo + end do !------------------! ! Spin-down sector ! @@ -149,9 +149,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -167,9 +167,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: bb @@ -186,9 +186,9 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -204,11 +204,11 @@ subroutine UGF2_reg_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo + end do Z(:,:) = 1d0/(1d0 - Z(:,:)) diff --git a/src/GF/UGF2_self_energy.f90 b/src/GF/UGF2_self_energy.f90 index 8c99cd2..d846b9a 100644 --- a/src/GF/UGF2_self_energy.f90 +++ b/src/GF/UGF2_self_energy.f90 @@ -58,9 +58,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -74,9 +74,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: aa @@ -91,9 +91,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -107,12 +107,12 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,1) = SigC(p,q,1) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do !------------------! ! Spin-down sector ! @@ -134,9 +134,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -150,9 +150,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: bb @@ -167,9 +167,9 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -183,12 +183,12 @@ subroutine UGF2_self_energy(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,eGF2,S SigC(p,q,2) = SigC(p,q,2) + num*eps/(eps**2 + eta**2) if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do Z(:,:) = 1d0/(1d0 - Z(:,:)) diff --git a/src/GF/UGF2_self_energy_diag.f90 b/src/GF/UGF2_self_energy_diag.f90 index 7d5972a..1a4cf13 100644 --- a/src/GF/UGF2_self_energy_diag.f90 +++ b/src/GF/UGF2_self_energy_diag.f90 @@ -56,9 +56,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -72,9 +72,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: aa @@ -89,9 +89,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -105,11 +105,11 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,1) = SigC(p,1) + num*eps/(eps**2 + eta**2) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo + end do !------------------! ! Spin-down sector ! @@ -130,9 +130,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Addition part: ab @@ -146,9 +146,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: bb @@ -163,9 +163,9 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do ! Removal part: ab @@ -179,11 +179,11 @@ subroutine UGF2_self_energy_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_bb,eHF,e SigC(p,2) = SigC(p,2) + num*eps/(eps**2 + eta**2) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + end do + end do + end do - enddo + end do Z(:,:) = 1d0/(1d0 - Z(:,:)) diff --git a/src/GF/evRGF3.f90 b/src/GF/evRGF3.f90 index 852cbc1..c056076 100644 --- a/src/GF/evRGF3.f90 +++ b/src/GF/evRGF3.f90 @@ -62,12 +62,12 @@ App(p,1) = App(p,1) & - (2d0*V(p,k,p,j) - V(p,k,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,k,i)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -82,12 +82,12 @@ App(p,2) = App(p,2) & + (2d0*V(p,c,p,b) - V(p,c,b,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(j,i,c,a)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -102,12 +102,12 @@ App(p,3) = App(p,3) & + (2d0*V(p,c,p,j) - V(p,c,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,c,i)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do App(:,4) = App(:,3) @@ -124,12 +124,12 @@ App(p,5) = App(p,5) & - (2d0*V(p,b,p,k) - V(p,b,k,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(i,j,k,a)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do App(:,6) = App(:,5) @@ -167,10 +167,10 @@ Bpp(p,1) = Bpp(p,1) & + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -182,10 +182,10 @@ Bpp(p,2) = Bpp(p,2) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Total second-order Green function @@ -208,12 +208,12 @@ Cpp(p,1) = Cpp(p,1) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,c,d)*V(p,i,c,d)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -228,12 +228,12 @@ Cpp(p,2) = Cpp(p,2) & + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,j,k)*V(p,i,j,k)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Cpp(:,3) = Cpp(:,2) @@ -249,12 +249,12 @@ Cpp(p,4) = Cpp(p,4) & + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(i,j,b,c)*V(p,a,b,c)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Cpp(:,5) = Cpp(:,4) @@ -270,12 +270,12 @@ Cpp(p,6) = Cpp(p,6) & - (2d0*V(p,a,k,l) - V(p,a,l,k))*V(k,l,i,j)*V(p,a,i,j)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do ! Frequency-dependent third-order contribution: "D" terms @@ -299,12 +299,12 @@ + V(p,i,b,a)*(V(a,j,i,c)*(4d0*V(p,j,b,c) - 2d0*V(p,j,c,b)) & + V(a,j,c,i)*( V(p,j,c,b) - 2d0*V(p,j,b,c)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do do p=nC+1,nBas-nR do i=nC+1,nO @@ -324,12 +324,12 @@ + V(p,i,a,c)*(V(a,b,i,j)*( V(p,b,j,c) - 2d0*V(p,b,c,j)) & + V(a,b,j,i)*( V(p,b,c,j) - 2d0*V(p,b,j,c)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Dpp(:,3) = Dpp(:,2) @@ -351,12 +351,12 @@ + V(p,a,j,k)*(V(j,i,a,b)*( V(p,i,b,k) - 2d0*V(p,i,k,b)) & + V(j,i,b,a)*( V(p,i,k,b) - 2d0*V(p,i,b,k)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do Dpp(:,5) = Dpp(:,4) @@ -378,12 +378,12 @@ - V(p,a,i,k)*(V(i,b,a,j)*( V(p,b,j,k) - 2d0*V(p,b,k,j)) & + V(i,b,j,a)*( V(p,b,k,j) - 2d0*V(p,b,j,k)))/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do ! Compute renormalization factor (if required) @@ -445,7 +445,7 @@ Sig3(:) = Z(:)*Sig3(:) - endif + end if ! Total third-order Green function @@ -474,7 +474,7 @@ nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main SCF loop !------------------------------------------------------------------------ @@ -491,6 +491,6 @@ stop - endif + end if end subroutine diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index d19c5d0..1735ffa 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -160,7 +160,7 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved if(minval(rcond(:)) < 1d-15) n_diis = 0 - endif + end if ! Save quasiparticles energy for next cycle @@ -170,7 +170,7 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -187,7 +187,7 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved stop - endif + end if ! Deallocate memory @@ -199,6 +199,6 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved print*,'!!! BSE2 NYI for evUGF2 !!!' - endif + end if end subroutine diff --git a/src/GF/print_G0F3.f90 b/src/GF/print_G0F3.f90 index 12cbba8..6096dbf 100644 --- a/src/GF/print_G0F3.f90 +++ b/src/GF/print_G0F3.f90 @@ -29,7 +29,7 @@ subroutine print_G0F3(nBas,nO,eHF,Z,eGF3) do x=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' write(*,'(2X,A27,F15.6)') 'G0F3 HOMO energy (eV):',eGF3(HOMO)*HaToeV diff --git a/src/GF/print_RG0F2.f90 b/src/GF/print_RG0F2.f90 index 428220c..2371147 100644 --- a/src/GF/print_RG0F2.f90 +++ b/src/GF/print_RG0F2.f90 @@ -38,7 +38,7 @@ subroutine print_RG0F2(nBas,nO,eHF,Sig,eGF,Z,ENuc,ERHF,Ec) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',eGF(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'G0F2 HOMO energy =',eGF(HOMO)*HaToeV,' eV' diff --git a/src/GF/print_UG0F2.f90 b/src/GF/print_UG0F2.f90 index c727f03..694283b 100644 --- a/src/GF/print_UG0F2.f90 +++ b/src/GF/print_UG0F2.f90 @@ -53,7 +53,7 @@ subroutine print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' diff --git a/src/GF/print_evGF3.f90 b/src/GF/print_evGF3.f90 index c8cb2ef..61ed1e2 100644 --- a/src/GF/print_evGF3.f90 +++ b/src/GF/print_evGF3.f90 @@ -29,7 +29,7 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',Z(p),'|',eGF3(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GF/print_evGGF2.f90 b/src/GF/print_evGGF2.f90 index a93be0e..c1312c6 100644 --- a/src/GF/print_evGGF2.f90 +++ b/src/GF/print_evGGF2.f90 @@ -40,7 +40,7 @@ subroutine print_evGGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF,ENuc,ERHF,Ec) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',eGF(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GF/print_evRGF2.f90 b/src/GF/print_evRGF2.f90 index 5d723dc..25c5b75 100644 --- a/src/GF/print_evRGF2.f90 +++ b/src/GF/print_evRGF2.f90 @@ -40,7 +40,7 @@ subroutine print_evRGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF,ENuc,ERHF,Ec) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',eGF(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GF/print_evUGF2.f90 b/src/GF/print_evUGF2.f90 index eefe183..625c2a7 100644 --- a/src/GF/print_evUGF2.f90 +++ b/src/GF/print_evUGF2.f90 @@ -45,7 +45,7 @@ subroutine print_evUGF2(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec) write(*,'(1X,A21,I1,A2,A12)')' Self-consistent evG',nSCF,'F2',' calculation' else write(*,'(1X,A21,I2,A2,A12)')' Self-consistent evG',nSCF,'F2',' calculation' - endif + end if write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & @@ -59,7 +59,7 @@ subroutine print_evUGF2(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' diff --git a/src/GF/print_qsRGF2.f90 b/src/GF/print_qsRGF2.f90 index 3d98502..845c5a0 100644 --- a/src/GF/print_qsRGF2.f90 +++ b/src/GF/print_qsRGF2.f90 @@ -47,7 +47,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, write(*,'(1X,A21,I1,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' else write(*,'(1X,A21,I2,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' @@ -56,7 +56,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, do q=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF @@ -110,7 +110,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, call matout(nBas,1,eGF) write(*,*) - endif + end if end subroutine diff --git a/src/GF/print_qsUGF2.f90 b/src/GF/print_qsUGF2.f90 index 165e7da..130e269 100644 --- a/src/GF/print_qsUGF2.f90 +++ b/src/GF/print_qsUGF2.f90 @@ -72,7 +72,7 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K, write(*,'(1X,A21,I1,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' else write(*,'(1X,A21,I2,A2,A12)')' Self-consistent qsG',nSCF,'F2',' calculation' - endif + end if write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & @@ -86,7 +86,7 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K, write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' @@ -173,6 +173,6 @@ subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K, call matout(nBas,nBas,cGF2(:,:,2)) write(*,*) - endif + end if end subroutine diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 545ba10..8222e56 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -233,7 +233,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -250,7 +250,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si stop - endif + end if ! Deallocate memory diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 7c8a066..6fa0bdb 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -306,7 +306,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,SigCp,Z,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -323,7 +323,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved stop - endif + end if ! Deallocate memory diff --git a/src/GT/GTeh_SigC.f90 b/src/GT/GTeh_SigC.f90 index 0cac37a..e616939 100644 --- a/src/GT/GTeh_SigC.f90 +++ b/src/GT/GTeh_SigC.f90 @@ -37,8 +37,8 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR) eps = w - e(i) + Om(m) num = rhoL(i,p,m)*rhoR(i,p,m) GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2) - enddo - enddo + end do + end do ! Virtual part of the correlation self-energy @@ -47,7 +47,7 @@ double precision function GTeh_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR) eps = w - e(a) - Om(m) num = rhoL(p,a,m)*rhoR(p,a,m) GTeh_SigC = GTeh_SigC + num*eps/(eps**2 + eta**2) - enddo - enddo + end do + end do end function diff --git a/src/GT/GTeh_dSigC.f90 b/src/GT/GTeh_dSigC.f90 index 513c23e..1bc8afc 100644 --- a/src/GT/GTeh_dSigC.f90 +++ b/src/GT/GTeh_dSigC.f90 @@ -37,8 +37,8 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR) eps = w - e(i) + Om(m) num = rhoL(i,p,m)*rhoR(i,p,m) GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo + end do + end do ! Virtual part of the correlation self-energy @@ -47,7 +47,7 @@ double precision function GTeh_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR) eps = w - e(a) - Om(m) num = rhoL(p,a,m)*rhoR(p,a,m) GTeh_dSigC = GTeh_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo + end do + end do end function diff --git a/src/GT/GTeh_excitation_density.f90 b/src/GT/GTeh_excitation_density.f90 index f2ab3d5..8c0916c 100644 --- a/src/GT/GTeh_excitation_density.f90 +++ b/src/GT/GTeh_excitation_density.f90 @@ -46,11 +46,11 @@ subroutine GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,XmY,rhoL,rhoR) rhoR(p,q,m) = rhoR(p,q,m) + (2d0*ERI(p,j,b,q) - ERI(p,j,q,b))*X + (2d0*ERI(p,b,j,q) - ERI(p,b,q,j))*Y - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do !$OMP END DO !$OMP END PARALLEL diff --git a/src/GT/GTeh_plot_self_energy.f90 b/src/GT/GTeh_plot_self_energy.f90 index 3245a06..25b5c9c 100644 --- a/src/GT/GTeh_plot_self_energy.f90 +++ b/src/GT/GTeh_plot_self_energy.f90 @@ -53,7 +53,7 @@ subroutine GTeh_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR) do g=1,nGrid w(g) = wmin + dble(g)*dw - enddo + end do ! Occupied part of the self-energy and renormalization factor @@ -73,8 +73,8 @@ subroutine GTeh_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR) do g=1,nGrid do p=nC+1,nBas-nR S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2) - enddo - enddo + end do + end do S(:,:) = S(:,:)/pi @@ -90,7 +90,7 @@ subroutine GTeh_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR) write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR) write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR) write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR) - enddo + end do ! Closing files diff --git a/src/GT/GTeh_regularization.f90 b/src/GT/GTeh_regularization.f90 index 5e3f7e3..54f5b6d 100644 --- a/src/GT/GTeh_regularization.f90 +++ b/src/GT/GTeh_regularization.f90 @@ -41,16 +41,16 @@ subroutine GTeh_regularization(nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR) kappa = 1d0 - exp(-Dpim*Dpim*s) rhoL(i,p,m) = kappa*rhoL(i,p,m) rhoR(i,p,m) = kappa*rhoR(i,p,m) - enddo + end do do a=nO+1,nBas-nR Dpam = e(p) - e(a) - Om(m) kappa = 1d0 - exp(-Dpam*Dpam*s) rhoL(p,a,m) = kappa*rhoL(p,a,m) rhoR(p,a,m) = kappa*rhoR(p,a,m) - enddo + end do - enddo - enddo + end do + end do end subroutine diff --git a/src/GT/GTpp_SigC.f90 b/src/GT/GTpp_SigC.f90 index c511524..3616111 100644 --- a/src/GT/GTpp_SigC.f90 +++ b/src/GT/GTpp_SigC.f90 @@ -43,14 +43,14 @@ double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt do cd=1,nVVs eps = w + e(i) - Om1s(cd) GTpp_SigC = GTpp_SigC + rho1s(p,i,cd)**2*eps/(eps**2 + eta**2) - enddo + end do do cd=1,nVVt eps = w + e(i) - Om1t(cd) GTpp_SigC = GTpp_SigC + rho1t(p,i,cd)**2*eps/(eps**2 + eta**2) - enddo + end do - enddo + end do !---------------------------------------------- ! Virtual part of the T-matrix self-energy @@ -61,13 +61,13 @@ double precision function GTpp_SigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt do kl=1,nOOs eps = w + e(a) - Om2s(kl) GTpp_SigC = GTpp_SigC + rho2s(p,a,kl)**2*eps/(eps**2 + eta**2) - enddo + end do do kl=1,nOOt eps = w + e(a) - Om2t(kl) GTpp_SigC = GTpp_SigC + rho2t(p,a,kl)**2*eps/(eps**2 + eta**2) - enddo + end do - enddo + end do end function diff --git a/src/GT/GTpp_dSigC.f90 b/src/GT/GTpp_dSigC.f90 index 1251ab4..a804a3f 100644 --- a/src/GT/GTpp_dSigC.f90 +++ b/src/GT/GTpp_dSigC.f90 @@ -43,14 +43,14 @@ double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV do cd=1,nVVs eps = w + e(i) - Om1s(cd) GTpp_dSigC = GTpp_dSigC - rho1s(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do cd=1,nVVt eps = w + e(i) - Om1t(cd) GTpp_dSigC = GTpp_dSigC - rho1t(p,i,cd)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo + end do !---------------------------------------------- ! Virtual part of the T-matrix self-energy @@ -62,13 +62,13 @@ double precision function GTpp_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVV do kl=1,nOOs eps = w + e(a) - Om2s(kl) GTpp_dSigC = GTpp_dSigC - rho2s(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do kl=1,nOOt eps = w + e(a) - Om2t(kl) GTpp_dSigC = GTpp_dSigC - rho2t(p,a,kl)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo + end do end function diff --git a/src/GT/GTpp_phBSE_static_kernel_A.f90 b/src/GT/GTpp_phBSE_static_kernel_A.f90 index fc51198..f3721a3 100644 --- a/src/GT/GTpp_phBSE_static_kernel_A.f90 +++ b/src/GT/GTpp_phBSE_static_kernel_A.f90 @@ -51,19 +51,19 @@ subroutine GTpp_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Ome eps = + Omega1(cd) chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) - enddo + end do do kl=1,nOO eps = - Omega2(kl) chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) - enddo + end do KA(ia,jb) = lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do !$omp end parallel do diff --git a/src/GT/GTpp_phBSE_static_kernel_B.f90 b/src/GT/GTpp_phBSE_static_kernel_B.f90 index 2d82a95..dbde257 100644 --- a/src/GT/GTpp_phBSE_static_kernel_B.f90 +++ b/src/GT/GTpp_phBSE_static_kernel_B.f90 @@ -48,18 +48,18 @@ subroutine GTpp_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Ome do cd=1,nVV eps = + Omega1(cd) chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) - enddo + end do do kl=1,nOO eps = - Omega2(kl) chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) - enddo + end do KB(ia,jb) = lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GT/GTpp_plot_self_energy.f90 b/src/GT/GTpp_plot_self_energy.f90 index 3d3b8d5..5672e3b 100644 --- a/src/GT/GTpp_plot_self_energy.f90 +++ b/src/GT/GTpp_plot_self_energy.f90 @@ -57,7 +57,7 @@ subroutine GTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om do g=1,nGrid w(g) = wmin + dble(g)*dw - enddo + end do ! Occupied part of the self-energy and renormalization factor @@ -77,8 +77,8 @@ subroutine GTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om do g=1,nGrid do p=nC+1,nBas-nR S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2) - enddo - enddo + end do + end do S(:,:) = S(:,:)/pi @@ -94,7 +94,7 @@ subroutine GTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR) write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR) write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR) - enddo + end do ! Closing files diff --git a/src/GT/GTpp_regularization.f90 b/src/GT/GTpp_regularization.f90 index 15fefd2..dd8dbfb 100644 --- a/src/GT/GTpp_regularization.f90 +++ b/src/GT/GTpp_regularization.f90 @@ -54,6 +54,6 @@ subroutine GTpp_regularization(nBas,nC,nO,nV,nR,nOO,nVV,e,Om1,rho1,Om2,rho2) end do end do - enddo + end do end subroutine diff --git a/src/GT/GTpp_self_energy.f90 b/src/GT/GTpp_self_energy.f90 index 49cdefb..45e940c 100644 --- a/src/GT/GTpp_self_energy.f90 +++ b/src/GT/GTpp_self_energy.f90 @@ -51,18 +51,18 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1 num = rho1s(p,i,cd)*rho1s(q,i,cd) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do cd=1,nVVt eps = e(p) + e(i) - Om1t(cd) num = rho1t(p,i,cd)*rho1t(q,i,cd) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo - enddo - enddo + end do + end do + end do !---------------------------------------------- ! Virtual part of the T-matrix self-energy @@ -77,18 +77,18 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1 num = rho2s(p,a,kl)*rho2s(q,a,kl) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do kl=1,nOOt eps = e(p) + e(a) - Om2t(kl) num = rho2t(p,a,kl)*rho2t(q,a,kl) Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo - enddo - enddo + end do + end do + end do !---------------------------------------------- ! Galitskii-Migdal correlation energy @@ -101,16 +101,16 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1 eps = e(i) + e(j) - Om1s(cd) num = rho1s(i,j,cd)*rho1s(i,j,cd) EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + end do do cd=1,nVVt eps = e(i) + e(j) - Om1t(cd) num = rho1t(i,j,cd)*rho1t(i,j,cd) EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + end do - enddo - enddo + end do + end do do a=nO+1,nBas-nR do b=nO+1,nBas-nR @@ -119,16 +119,16 @@ subroutine GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s,rho1 eps = e(a) + e(b) - Om2s(kl) num = rho2s(a,b,kl)*rho2s(a,b,kl) EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + end do do kl=1,nOOt eps = e(a) + e(b) - Om2t(kl) num = rho2t(a,b,kl)*rho2t(a,b,kl) EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + end do - enddo - enddo + end do + end do Z(:) = 1d0/(1d0 - Z(:)) diff --git a/src/GT/GTpp_self_energy_diag.f90 b/src/GT/GTpp_self_energy_diag.f90 index 9cd2a63..8675f52 100644 --- a/src/GT/GTpp_self_energy_diag.f90 +++ b/src/GT/GTpp_self_energy_diag.f90 @@ -51,17 +51,17 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s num = rho1s(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do cd=1,nVVt eps = e(p) + e(i) - Om1t(cd) num = rho1t(p,i,cd)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo - enddo + end do + end do !---------------------------------------------- ! Virtual part of the T-matrix self-energy @@ -75,17 +75,17 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s num = rho2s(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do do kl=1,nOOt eps = e(p) + e(a) - Om2t(kl) num = rho2t(p,a,kl)**2 Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do - enddo - enddo + end do + end do !---------------------------------------------- ! Galitskii-Migdal correlation energy @@ -98,16 +98,16 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s eps = e(i) + e(j) - Om1s(cd) num = rho1s(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + end do do cd=1,nVVt eps = e(i) + e(j) - Om1t(cd) num = rho1t(i,j,cd)**2 EcGM = EcGM + num*eps/(eps**2 + eta**2) - enddo + end do - enddo - enddo + end do + end do do a=nO+1,nBas-nR do b=nO+1,nBas-nR @@ -116,16 +116,16 @@ subroutine GTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,Om1s eps = e(a) + e(b) - Om2s(kl) num = rho2s(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + end do do kl=1,nOOt eps = e(a) + e(b) - Om2t(kl) num = rho2t(a,b,kl)**2 EcGM = EcGM - num*eps/(eps**2 + eta**2) - enddo + end do - enddo - enddo + end do + end do Z(:) = 1d0/(1d0 - Z(:)) diff --git a/src/GT/evRGTpp.f90 b/src/GT/evRGTpp.f90 index 0eee598..4df94b4 100644 --- a/src/GT/evRGTpp.f90 +++ b/src/GT/evRGTpp.f90 @@ -243,7 +243,7 @@ subroutine evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ diff --git a/src/GT/evUGTpp.f90 b/src/GT/evUGTpp.f90 index 0a0061d..2f8cff4 100644 --- a/src/GT/evUGTpp.f90 +++ b/src/GT/evUGTpp.f90 @@ -275,7 +275,7 @@ subroutine evUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main loop diff --git a/src/GT/print_RG0T0eh.f90 b/src/GT/print_RG0T0eh.f90 index c57051c..49a45b5 100644 --- a/src/GT/print_RG0T0eh.f90 +++ b/src/GT/print_RG0T0eh.f90 @@ -36,7 +36,7 @@ subroutine print_RG0T0eh(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'RG0T0eh HOMO energy =',eGT(HOMO)*HaToeV,' eV' diff --git a/src/GT/print_RG0T0pp.f90 b/src/GT/print_RG0T0pp.f90 index de5e911..776a434 100644 --- a/src/GT/print_RG0T0pp.f90 +++ b/src/GT/print_RG0T0pp.f90 @@ -45,7 +45,7 @@ subroutine print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigT(p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'G0T0pp@RHF HOMO energy = ',eGT(HOMO)*HaToeV,' eV' diff --git a/src/GT/print_UG0T0.f90 b/src/GT/print_UG0T0.f90 index 19df30c..942b1dc 100644 --- a/src/GT/print_UG0T0.f90 +++ b/src/GT/print_UG0T0.f90 @@ -48,7 +48,7 @@ subroutine print_UG0T0(nBas,nO,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigT(p,1)*HaToeV,SigT(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGT(p,1)*HaToeV,eGT(p,2)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F15.6,A3)') 'UG0T0 HOMO energy (eV) =',maxval(HOMO(:))*HaToeV,' eV' diff --git a/src/GT/print_evRGTeh.f90 b/src/GT/print_evRGTeh.f90 index 51a0780..88552c5 100644 --- a/src/GT/print_evRGTeh.f90 +++ b/src/GT/print_evRGTeh.f90 @@ -32,7 +32,7 @@ subroutine print_evRGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent evG',nSCF,'Teh',nSCF,' calculation' else write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent evG',nSCF,'Teh',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' @@ -41,7 +41,7 @@ subroutine print_evRGTeh(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGT,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GT/print_evRGTpp.f90 b/src/GT/print_evRGTpp.f90 index 4453cc3..f1d40cf 100644 --- a/src/GT/print_evRGTpp.f90 +++ b/src/GT/print_evRGTpp.f90 @@ -34,7 +34,7 @@ subroutine print_evRGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent evG',nSCF,'Tpp',nSCF,' calculation' else write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent evG',nSCF,'Tpp',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' @@ -43,7 +43,7 @@ subroutine print_evRGTpp(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigT(p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GT/print_evUGT.f90 b/src/GT/print_evUGT.f90 index ba7c4c6..7003c19 100644 --- a/src/GT/print_evUGT.f90 +++ b/src/GT/print_evUGT.f90 @@ -44,7 +44,7 @@ subroutine print_evUGT(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation' else write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|' @@ -54,7 +54,7 @@ subroutine print_evUGT(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigT,Z,eGT,EcGM,EcRPA) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigT(p,1)*HaToeV,SigT(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGT(p,1)*HaToeV,eGT(p,2)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GT/print_qsRGTeh.f90 b/src/GT/print_qsRGTeh.f90 index 75aefd3..9f8f266 100644 --- a/src/GT/print_qsRGTeh.f90 +++ b/src/GT/print_qsRGTeh.f90 @@ -52,7 +52,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent qsG',nSCF,'Teh',nSCF,' calculation' else write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent qsG',nSCF,'Teh',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' @@ -61,7 +61,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF @@ -118,6 +118,6 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ call vecout(nBas,eGT) write(*,*) - endif + end if end subroutine diff --git a/src/GT/print_qsRGTpp.f90 b/src/GT/print_qsRGTpp.f90 index c0a1ad6..c3bebfa 100644 --- a/src/GT/print_qsRGTpp.f90 +++ b/src/GT/print_qsRGTpp.f90 @@ -52,7 +52,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent qsG',nSCF,'Tpp',nSCF,' calculation' else write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent qsG',nSCF,'Tpp',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' @@ -61,7 +61,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF @@ -118,6 +118,6 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ call vecout(nBas,eGT) write(*,*) - endif + end if end subroutine diff --git a/src/GT/print_qsUGT.f90 b/src/GT/print_qsUGT.f90 index 919ecca..af4280c 100644 --- a/src/GT/print_qsUGT.f90 +++ b/src/GT/print_qsUGT.f90 @@ -51,7 +51,7 @@ subroutine print_qsUGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigT,Z,ENuc,ET,EV,EJ,E write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' else write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' - endif + end if write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|' @@ -61,7 +61,7 @@ subroutine print_qsUGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigT,Z,ENuc,ET,EV,EJ,E write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigT(p,p,1)*HaToeV,SigT(p,p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGT(p,1)*HaToeV,eGT(p,2)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90 index dce484f..b814001 100644 --- a/src/GT/qsRGTeh.f90 +++ b/src/GT/qsRGTeh.f90 @@ -257,7 +257,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -274,7 +274,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d stop - endif + end if ! Deallocate memory diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 8211468..2ff7165 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -300,7 +300,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -317,7 +317,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d stop - endif + end if ! Deallocate memory diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index 752e755..5f539bb 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -385,7 +385,7 @@ write(*,*) 'EcGM', EcGM(1) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsUGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigTp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -402,7 +402,7 @@ write(*,*) 'EcGM', EcGM(1) stop - endif + end if ! Free memory diff --git a/src/GW/GGW_phBSE_static_kernel_A.f90 b/src/GW/GGW_phBSE_static_kernel_A.f90 index 6d38245..9ebf100 100644 --- a/src/GW/GGW_phBSE_static_kernel_A.f90 +++ b/src/GW/GGW_phBSE_static_kernel_A.f90 @@ -44,13 +44,13 @@ subroutine GGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,K do kc=1,nS eps = Om(kc)**2 + eta**2 chi = chi + rho(i,j,kc)*rho(a,b,kc)*Om(kc)/eps - enddo + end do KA(ia,jb) = 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GGW_phBSE_static_kernel_B.f90 b/src/GW/GGW_phBSE_static_kernel_B.f90 index 5afc33f..c5081df 100644 --- a/src/GW/GGW_phBSE_static_kernel_B.f90 +++ b/src/GW/GGW_phBSE_static_kernel_B.f90 @@ -44,13 +44,13 @@ subroutine GGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,K do kc=1,nS eps = Om(kc)**2 + eta**2 chi = chi + rho(i,b,kc)*rho(a,j,kc)*Om(kc)/eps - enddo + end do KB(ia,jb) = 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GW_excitation_density.f90 b/src/GW/GW_excitation_density.f90 index a6ea714..1b26f66 100644 --- a/src/GW/GW_excitation_density.f90 +++ b/src/GW/GW_excitation_density.f90 @@ -36,11 +36,11 @@ subroutine GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) jb = jb + 1 do ia=1,nS rho(p,q,ia) = rho(p,q,ia) + ERI(p,j,q,b)*XpY(ia,jb) - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do !$OMP END DO !$OMP END PARALLEL diff --git a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 index 64cadc6..868a4a7 100644 --- a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 @@ -58,8 +58,8 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,KA_dyn KA_dyn(ia,jb) = KA_dyn(ia,jb) - num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - enddo - enddo + end do + end do do c=nO+1,nBas-nR do d=nO+1,nBas-nR @@ -70,8 +70,8 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,KA_dyn KA_dyn(ia,jb) = KA_dyn(ia,jb) + num*dem/(dem**2 + eta**2) ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 - enddo - enddo + end do + end do do k=nC+1,nO do l=nC+1,nO @@ -85,10 +85,10 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,KA_dyn end do end do - enddo - enddo - enddo - enddo + end do + end do + end do + end do !$omp end parallel do end subroutine diff --git a/src/GW/GW_phBSE2_dynamic_kernel_B.f90 b/src/GW/GW_phBSE2_dynamic_kernel_B.f90 index bc7d5b5..3710480 100644 --- a/src/GW/GW_phBSE2_dynamic_kernel_B.f90 +++ b/src/GW/GW_phBSE2_dynamic_kernel_B.f90 @@ -54,8 +54,8 @@ subroutine GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn) KB_dyn(ia,jb) = KB_dyn(ia,jb) - num*dem/(dem**2 + eta**2) - enddo - enddo + end do + end do do c=nO+1,nBas-nR do d=nO+1,nBas-nR @@ -65,8 +65,8 @@ subroutine GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn) KB_dyn(ia,jb) = KB_dyn(ia,jb) + num*dem/(dem**2 + eta**2) - enddo - enddo + end do + end do do k=nC+1,nO do l=nC+1,nO @@ -79,10 +79,10 @@ subroutine GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn) end do end do - enddo - enddo - enddo - enddo + end do + end do + end do + end do !$omp end parallel do end subroutine diff --git a/src/GW/GW_phBSE_dynamic_kernel_A.f90 b/src/GW/GW_phBSE_dynamic_kernel_A.f90 index aaa528e..f56dbdf 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90 @@ -58,7 +58,7 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) - enddo + end do KA_dyn(ia,jb) = KA_dyn(ia,jb) - 2d0*lambda*chi @@ -71,14 +71,14 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do !$omp end parallel do diff --git a/src/GW/GW_phBSE_dynamic_kernel_B.f90 b/src/GW/GW_phBSE_dynamic_kernel_B.f90 index 7f82509..34801f2 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_B.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_B.f90 @@ -52,13 +52,13 @@ subroutine GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh eps = - OmRPA(kc) - (eGW(a) - eGW(j)) chi = chi + rho(i,b,kc)*rho(j,a,kc)*eps/(eps**2 + eta**2) - enddo + end do KB(ia,jb) = KB(ia,jb) - 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GW_phBSE_static_kernel.f90 b/src/GW/GW_phBSE_static_kernel.f90 index 1154f0a..2f00373 100644 --- a/src/GW/GW_phBSE_static_kernel.f90 +++ b/src/GW/GW_phBSE_static_kernel.f90 @@ -46,13 +46,13 @@ subroutine GW_phBSE_static_kernel(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,W) do m=1,nS dem = Om(m)**2 + eta**2 chi = chi + rho(p,q,m)*rho(r,s,m)*Om(m)/dem - enddo + end do W(p,s,q,r) = - lambda*ERI(p,s,q,r) + 4d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GW_phBSE_static_kernel_A.f90 b/src/GW/GW_phBSE_static_kernel_A.f90 index fb3ec9d..c2de92c 100644 --- a/src/GW/GW_phBSE_static_kernel_A.f90 +++ b/src/GW/GW_phBSE_static_kernel_A.f90 @@ -44,13 +44,13 @@ subroutine GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KA do kc=1,nS eps = Om(kc)**2 + eta**2 chi = chi + rho(i,j,kc)*rho(a,b,kc)*Om(kc)/eps - enddo + end do KA(ia,jb) = 4d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GW_phBSE_static_kernel_B.f90 b/src/GW/GW_phBSE_static_kernel_B.f90 index 3554de1..3258493 100644 --- a/src/GW/GW_phBSE_static_kernel_B.f90 +++ b/src/GW/GW_phBSE_static_kernel_B.f90 @@ -44,13 +44,13 @@ subroutine GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KB do kc=1,nS eps = Om(kc)**2 + eta**2 chi = chi + rho(i,b,kc)*rho(a,j,kc)*Om(kc)/eps - enddo + end do KB(ia,jb) = 4d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine diff --git a/src/GW/GW_plot_self_energy.f90 b/src/GW/GW_plot_self_energy.f90 index 44b0dea..5c60440 100644 --- a/src/GW/GW_plot_self_energy.f90 +++ b/src/GW/GW_plot_self_energy.f90 @@ -52,7 +52,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) do g=1,nGrid w(g) = wmin + dble(g)*dw - enddo + end do ! Occupied part of the self-energy and renormalization factor @@ -72,8 +72,8 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) do g=1,nGrid do p=nC+1,nBas-nR S(p,g) = eta/((w(g) - eHF(p) - SigC(p,g))**2 + eta**2) - enddo - enddo + end do + end do S(:,:) = S(:,:)/pi @@ -89,7 +89,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho) write(9 ,*) w(g)*HaToeV,((w(g)-eHF(p))*HaToeV,p=nC+1,nBas-nR) write(10,*) w(g)*HaToeV,(Z(p,g),p=nC+1,nBas-nR) write(11,*) w(g)*HaToeV,(S(p,g),p=nC+1,nBas-nR) - enddo + end do ! Closing files diff --git a/src/GW/GW_ppBSE_static_kernel_B.f90 b/src/GW/GW_ppBSE_static_kernel_B.f90 index 2752782..7da42cb 100644 --- a/src/GW/GW_ppBSE_static_kernel_B.f90 +++ b/src/GW/GW_ppBSE_static_kernel_B.f90 @@ -57,7 +57,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda eps = Om(m)**2 + eta**2 chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps & + rho(i,b,m)*rho(a,j,m)*Om(m)/eps - enddo + end do KB(ab,ij) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) @@ -88,7 +88,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda eps = Om(m)**2 + eta**2 chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps & + rho(i,b,m)*rho(a,j,m)*Om(m)/eps - enddo + end do KB(ab,ij) = 2d0*lambda*chi @@ -119,7 +119,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda eps = Om(m)**2 + eta**2 chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps & + rho(i,b,m)*rho(a,j,m)*Om(m)/eps - enddo + end do KB(ab,ij) = lambda*chi diff --git a/src/GW/GW_ppBSE_static_kernel_C.f90 b/src/GW/GW_ppBSE_static_kernel_C.f90 index bd015e3..82c6e58 100644 --- a/src/GW/GW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GW_ppBSE_static_kernel_C.f90 @@ -56,7 +56,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + rho(a,d,m)*rho(b,c,m)*Om(m)/eps - enddo + end do KC(ab,cd) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) @@ -87,7 +87,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + rho(a,d,m)*rho(b,c,m)*Om(m)/eps - enddo + end do KC(ab,cd) = 2d0*lambda*chi @@ -118,7 +118,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & + rho(a,d,m)*rho(b,c,m)*Om(m)/eps - enddo + end do KC(ab,cd) = lambda*chi diff --git a/src/GW/GW_ppBSE_static_kernel_D.f90 b/src/GW/GW_ppBSE_static_kernel_D.f90 index f15db35..e0838c7 100644 --- a/src/GW/GW_ppBSE_static_kernel_D.f90 +++ b/src/GW/GW_ppBSE_static_kernel_D.f90 @@ -56,7 +56,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + rho(i,l,m)*rho(j,k,m)*Om(m)/eps - enddo + end do KD(ij,kl) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) @@ -87,7 +87,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + rho(i,l,m)*rho(j,k,m)*Om(m)/eps - enddo + end do KD(ij,kl) = 2d0*lambda*chi @@ -118,7 +118,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI eps = Om(m)**2 + eta**2 chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + rho(i,l,m)*rho(j,k,m)*Om(m)/eps - enddo + end do KD(ij,kl) = lambda*chi diff --git a/src/GW/GW_regularization.f90 b/src/GW/GW_regularization.f90 index 6e33c60..ed77ae4 100644 --- a/src/GW/GW_regularization.f90 +++ b/src/GW/GW_regularization.f90 @@ -39,15 +39,15 @@ subroutine GW_regularization(nBas,nC,nO,nV,nR,nS,e,Om,rho) Dpim = e(p) - e(i) + Om(m) kappa = 1d0 - exp(-Dpim*Dpim*s) rho(p,i,m) = kappa*rho(p,i,m) - enddo + end do do a=nO+1,nBas-nR Dpam = e(p) - e(a) - Om(m) kappa = 1d0 - exp(-Dpam*Dpam*s) rho(p,a,m) = kappa*rho(p,a,m) - enddo + end do - enddo - enddo + end do + end do end subroutine diff --git a/src/GW/RG0W0.f90 b/src/GW/RG0W0.f90 index 051e3ce..cba3ac9 100644 --- a/src/GW/RG0W0.f90 +++ b/src/GW/RG0W0.f90 @@ -56,6 +56,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA double precision,allocatable :: XmY(:,:) double precision,allocatable :: rho(:,:,:) + double precision,allocatable :: W(:,:,:,:) + double precision,allocatable :: eGWlin(:) double precision,allocatable :: eGW(:) @@ -161,10 +163,6 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA call print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) -! Deallocate memory - - deallocate(SigC,Z,Om,XpY,XmY,rho) - ! Perform BSE calculation if(dophBSE) then @@ -235,6 +233,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA end if +! if(.true.) then + +! allocate(W(nBas,nBas,nBas,nBas)) +! call GW_phBSE_static_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,Om,rho,W) +! call pCCD(dotest,264,1d-7,5,nBas,nC,nO,nV,nR,ERI,W,ERHF,eGW) +! deallocate(W) + +! end if + ! Testing zone if(dotest) then diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 274f928..b0965fe 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -283,7 +283,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -300,7 +300,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, stop - endif + end if print *, "Wall time for Fock and exchange build", tt print *, "Wall Time for AO to MO", tao diff --git a/src/GW/UGW_excitation_density.f90 b/src/GW/UGW_excitation_density.f90 index 571051b..c6d2461 100644 --- a/src/GW/UGW_excitation_density.f90 +++ b/src/GW/UGW_excitation_density.f90 @@ -47,9 +47,9 @@ subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ER rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aaaa(p,j,q,b)*XpY(ia,jb) - enddo - enddo - enddo + end do + end do + end do ! Opposite-spin contribution do ia=1,nSt @@ -60,12 +60,12 @@ subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ER rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aabb(p,j,q,b)*XpY(ia,jb) - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do !------------! ! Beta block ! @@ -83,9 +83,9 @@ subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ER rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb) - enddo - enddo - enddo + end do + end do + end do ! Same-spin contribution do ia=1,nSt @@ -96,11 +96,11 @@ subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ER rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bbbb(p,j,q,b)*XpY(ia,jb) - enddo - enddo - enddo + end do + end do + end do - enddo - enddo + end do + end do end subroutine diff --git a/src/GW/UGW_phBSE_dynamic_kernel_A.f90 b/src/GW/UGW_phBSE_dynamic_kernel_A.f90 index 06bc3de..5231a4d 100644 --- a/src/GW/UGW_phBSE_dynamic_kernel_A.f90 +++ b/src/GW/UGW_phBSE_dynamic_kernel_A.f90 @@ -61,7 +61,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ chi = 0d0 do kc=1,nS_sc chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - enddo + end do A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi @@ -74,7 +74,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,1)) chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*eps/(eps**2 + eta**2) - enddo + end do A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi @@ -87,14 +87,14 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,1)) chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! bbbb block @@ -110,7 +110,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ chi = 0d0 do kc=1,nS_sc chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - enddo + end do A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi @@ -123,7 +123,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,2)) chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*eps/(eps**2 + eta**2) - enddo + end do A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi @@ -136,14 +136,14 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,2)) chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end if @@ -167,7 +167,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ chi = 0d0 do kc=1,nS_sc chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - enddo + end do A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi @@ -180,7 +180,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*eps/(eps**2 + eta**2) - enddo + end do A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi @@ -193,7 +193,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi @@ -216,7 +216,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ chi = 0d0 do kc=1,nS_sc chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - enddo + end do A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi @@ -229,7 +229,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*eps/(eps**2 + eta**2) - enddo + end do A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi @@ -242,7 +242,7 @@ subroutine UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_ eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo + end do ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi diff --git a/src/GW/UGW_phBSE_static_kernel_A.f90 b/src/GW/UGW_phBSE_static_kernel_A.f90 index b37bb4b..fbd04b1 100644 --- a/src/GW/UGW_phBSE_static_kernel_A.f90 +++ b/src/GW/UGW_phBSE_static_kernel_A.f90 @@ -58,14 +58,14 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps - enddo + end do A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! bbbb block @@ -82,14 +82,14 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps - enddo + end do A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end if @@ -114,7 +114,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps - enddo + end do A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi @@ -138,7 +138,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps - enddo + end do A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi diff --git a/src/GW/UGW_phBSE_static_kernel_B.f90 b/src/GW/UGW_phBSE_static_kernel_B.f90 index 81cb378..c357ecc 100644 --- a/src/GW/UGW_phBSE_static_kernel_B.f90 +++ b/src/GW/UGW_phBSE_static_kernel_B.f90 @@ -57,14 +57,14 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps - enddo + end do B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! bbbb block @@ -82,14 +82,14 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps - enddo + end do B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi - enddo - enddo - enddo - enddo + end do + end do + end do + end do end if @@ -115,7 +115,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps - enddo + end do B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi @@ -139,7 +139,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s do kc=1,nS_sc eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps - enddo + end do B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index c0cbde1..63e323c 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -225,7 +225,7 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ diff --git a/src/GW/print_GG0W0.f90 b/src/GW/print_GG0W0.f90 index 0bf83ef..acf4142 100644 --- a/src/GW/print_GG0W0.f90 +++ b/src/GW/print_GG0W0.f90 @@ -36,7 +36,7 @@ subroutine print_GG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'G0W0@GHF HOMO energy = ',eGW(HOMO)*HaToeV,' eV' diff --git a/src/GW/print_RG0W0.f90 b/src/GW/print_RG0W0.f90 index bc9ea28..a786d26 100644 --- a/src/GW/print_RG0W0.f90 +++ b/src/GW/print_RG0W0.f90 @@ -36,7 +36,7 @@ subroutine print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') 'G0W0@RHF HOMO energy = ',eGW(HOMO)*HaToeV,' eV' diff --git a/src/GW/print_UG0W0.f90 b/src/GW/print_UG0W0.f90 index caafa6e..9ec1c25 100644 --- a/src/GW/print_UG0W0.f90 +++ b/src/GW/print_UG0W0.f90 @@ -54,7 +54,7 @@ subroutine print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' diff --git a/src/GW/print_evGGW.f90 b/src/GW/print_evGGW.f90 index 2fd8550..9ff3b39 100644 --- a/src/GW/print_evGGW.f90 +++ b/src/GW/print_evGGW.f90 @@ -43,7 +43,7 @@ subroutine print_evGGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GW/print_evRGW.f90 b/src/GW/print_evRGW.f90 index 21e9edc..072b755 100644 --- a/src/GW/print_evRGW.f90 +++ b/src/GW/print_evRGW.f90 @@ -43,7 +43,7 @@ subroutine print_evRGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF diff --git a/src/GW/print_evUGW.f90 b/src/GW/print_evUGW.f90 index e182663..55e24b9 100644 --- a/src/GW/print_evUGW.f90 +++ b/src/GW/print_evUGW.f90 @@ -62,7 +62,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA,EcGM) write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' diff --git a/src/GW/print_qsGGW.f90 b/src/GW/print_qsGGW.f90 index 19f4721..c4136d8 100644 --- a/src/GW/print_qsGGW.f90 +++ b/src/GW/print_qsGGW.f90 @@ -128,7 +128,7 @@ subroutine print_qsGGW(nBas,nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,S,SigC,Z,ENuc,ET do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF @@ -185,7 +185,7 @@ subroutine print_qsGGW(nBas,nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,S,SigC,Z,ENuc,ET call vecout(nBas,eGW) write(*,*) - endif + end if end subroutine diff --git a/src/GW/print_qsRGW.f90 b/src/GW/print_qsRGW.f90 index 533d9c6..d09a670 100644 --- a/src/GW/print_qsRGW.f90 +++ b/src/GW/print_qsRGW.f90 @@ -62,7 +62,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF @@ -119,6 +119,6 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E call vecout(nBas,eGW) write(*,*) - endif + end if end subroutine diff --git a/src/GW/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 index 2078dde..3638c59 100644 --- a/src/GW/print_qsUGW.f90 +++ b/src/GW/print_qsUGW.f90 @@ -83,7 +83,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,ENuc,ET,EV,EJ,Ex,Ec write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' - enddo + end do write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90 index b12f1a4..230bd73 100644 --- a/src/GW/qsGGW.f90 +++ b/src/GW/qsGGW.f90 @@ -361,7 +361,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsGGW(nBas,nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -378,7 +378,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop stop - endif + end if ! Perform BSE calculation diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index ed9d762..6e10f75 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -249,7 +249,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -266,7 +266,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop stop - endif + end if ! Deallocate memory diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index e176cc7..ea57d9e 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -332,7 +332,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,EK,EcGM,EcRPA(ispin),EqsGW,SigCp,Z,dipole) - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -349,7 +349,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE stop - endif + end if ! Deallocate memory diff --git a/src/GW/ufBSE.f90 b/src/GW/ufBSE.f90 index 7034624..24568f3 100644 --- a/src/GW/ufBSE.f90 +++ b/src/GW/ufBSE.f90 @@ -207,7 +207,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) if(Z(s) > 1d-7) & write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|' - enddo + end do write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 index aa7362f..5a3cdea 100644 --- a/src/HF/GHF_search.f90 +++ b/src/HF/GHF_search.f90 @@ -170,7 +170,7 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then diff --git a/src/HF/GHF_stability.f90 b/src/HF/GHF_stability.f90 index f993147..a69cf2f 100644 --- a/src/HF/GHF_stability.f90 +++ b/src/HF/GHF_stability.f90 @@ -55,7 +55,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then @@ -90,7 +90,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then diff --git a/src/HF/MOM_overlap.f90 b/src/HF/MOM_overlap.f90 index 0e97fb4..4eb66ec 100644 --- a/src/HF/MOM_overlap.f90 +++ b/src/HF/MOM_overlap.f90 @@ -27,8 +27,8 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON) do i=1,nBas do j=1,nBas pOv(j) = pOv(j) + Ov(i,j)**2 - enddo - enddo + end do + end do pOv(:) = sqrt(pOV(:)) @@ -41,7 +41,7 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON) ploc = maxloc(pOv,nBas) ON(ploc) = 1d0 pOv(ploc) = 0d0 - enddo + end do ! print*,'--- Occupation numbers ---' diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index d44e555..27393a4 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -141,7 +141,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index 7a72560..d7b5385 100644 --- a/src/HF/RHF_stability.f90 +++ b/src/HF/RHF_stability.f90 @@ -55,7 +55,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then @@ -92,7 +92,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then @@ -132,7 +132,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om(:)) < 0d0) then diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 index 9d48dd6..a49746d 100644 --- a/src/HF/RMOM.f90 +++ b/src/HF/RMOM.f90 @@ -149,14 +149,14 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) Gap = 0d0 - endif + end if ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -174,7 +174,7 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) stop - endif + end if write(*,*) write(*,*) ' --- Final MO occupations --- ' diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index 6b16f62..0f597ee 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -166,7 +166,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu do ia=1,min(nS_sc,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om_sc(:)) < 0d0) then diff --git a/src/HF/UHF_stability.f90 b/src/HF/UHF_stability.f90 index 0d18e14..22fa9af 100644 --- a/src/HF/UHF_stability.f90 +++ b/src/HF/UHF_stability.f90 @@ -68,7 +68,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) do ia=1,min(nS_sc,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om_sc(:)) < 0d0) then @@ -103,7 +103,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) do ia=1,min(nS_sc,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om_sc(:)) < 0d0) then @@ -153,7 +153,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) do ia=1,min(nS_sf,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Om_sf(ia),'|',Om_sf(ia)*HaToeV,'|' - enddo + end do write(*,*)'-------------------------------------------------------------' if(minval(Om_sf(:)) < 0d0) then diff --git a/src/HF/density.f90 b/src/HF/density.f90 index a365405..aacc6c5 100644 --- a/src/HF/density.f90 +++ b/src/HF/density.f90 @@ -28,8 +28,8 @@ subroutine density(nGrid,nBas,P,AO,rho) rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG) - enddo - enddo - enddo + end do + end do + end do end subroutine diff --git a/src/HF/density_matrix.f90 b/src/HF/density_matrix.f90 index 5d40950..b7902f7 100644 --- a/src/HF/density_matrix.f90 +++ b/src/HF/density_matrix.f90 @@ -23,8 +23,8 @@ subroutine density_matrix(nBas,ON,c,P) do nu=1,nBas do mu=1,nBas P(mu,nu) = P(mu,nu) + 2d0*ON(i)*c(mu,i)*c(nu,i) - enddo - enddo - enddo + end do + end do + end do end subroutine diff --git a/src/HF/dipole_moment.f90 b/src/HF/dipole_moment.f90 index 656a39b..8805851 100644 --- a/src/HF/dipole_moment.f90 +++ b/src/HF/dipole_moment.f90 @@ -45,9 +45,9 @@ subroutine dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) do nu=1,nBas do mu=1,nBas dipole(ixyz) = dipole(ixyz) - P(mu,nu)*dipole_int(mu,nu,ixyz) - enddo - enddo + end do + end do - enddo + end do end subroutine diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index 80efaa8..658c853 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -39,8 +39,8 @@ subroutine huckel_guess(nBas,S,Hc,X,c) F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) F(nu,mu) = F(mu,nu) - enddo - enddo + end do + end do call core_guess(nBas,F,X,c) diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index 6637651..fef8d87 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -43,6 +43,6 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) print*,'Wrong guess option' stop - endif + end if end subroutine diff --git a/src/HF/print_GHF_spin.f90 b/src/HF/print_GHF_spin.f90 index ea436c4..88f525c 100644 --- a/src/HF/print_GHF_spin.f90 +++ b/src/HF/print_GHF_spin.f90 @@ -31,8 +31,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) do j = 1, nBas Ca(j,i) = C( j,i) Cb(j,i) = C(nBas+j,i) - enddo - enddo + end do + end do ! TODO DGEMM allocate(Paa(nO,nO), Pab(nO,nO), Pba(nO,nO), Pbb(nO,nO)) @@ -48,22 +48,22 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) do i = 1, nO Na = Na + Paa(i,i) Nb = Nb + Pbb(i,i) - enddo + end do nonco_z = dble(nO) do j = 1, nO do i = 1, nO nonco_z = nonco_z - (Paa(i,j) - Pbb(i,j))**2 - enddo - enddo + end do + end do nonco_z = 0.25d0 * nonco_z contam_ghf = 0.d0 do j = 1, nO do i = 1, nO contam_ghf = contam_ghf + (Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)) - enddo - enddo + end do + end do Sz = 0.5d0 * (Na - Nb) Sz2 = Sz*Sz + nonco_z @@ -85,7 +85,7 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) Sc_x = Sc_x + (+0.5d0,0.d0) * (Pab(i,i) + Pba(i,i)) Sc_y = Sc_y + (0.d0,-0.5d0) * (Pab(i,i) - Pba(i,i)) Sc_z = Sc_z + (+0.5d0,0.d0) * (Paa(i,i) - Pbb(i,i)) - enddo + end do write(*,'(A15,2F10.6)') ' < Sx > = ',Sc_x write(*,'(A15,2F10.6)') ' < Sy > = ',Sc_y write(*,'(A15,2F10.6)') ' < Sz > = ',Sc_z @@ -99,8 +99,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) Sc_xx = Sc_xx - zabs((+0.5d0,0.d0) * (Pab(i,j) + Pba(i,j)))**2 Sc_yy = Sc_yy - zabs((0.d0,-0.5d0) * (Pab(i,j) - Pba(i,j)))**2 Sc_zz = Sc_zz - zabs((+0.5d0,0.d0) * (Paa(i,j) - Pbb(i,j)))**2 - enddo - enddo + end do + end do write(*,'(A15,2F10.6)') ' < Sx^2 > = ',Sc_xx write(*,'(A15,2F10.6)') ' < Sy^2 > = ',Sc_yy write(*,'(A15,2F10.6)') ' < Sz^2 > = ',Sc_zz @@ -114,8 +114,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) do j = 1, nO Sc_xy = Sc_xy - (+0.5d0,0.d0) * (Pab(i,j) + Pba(i,j)) * (0.d0,-0.5d0) * (Pab(j,i) - Pba(j,i)) Sc_yx = Sc_yx - (+0.5d0,0.d0) * (Pab(j,i) + Pba(j,i)) * (0.d0,-0.5d0) * (Pab(i,j) - Pba(i,j)) - enddo - enddo + end do + end do write(*,'(A15,2F10.6)') ' < Sx.Sy > = ',Sc_xy write(*,'(A15,2F10.6)') ' < Sy.Sx > = ',Sc_yx @@ -127,8 +127,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) do j = 1, nO Sc_xz = Sc_xz - (+0.5d0,0.d0) * (Pab(i,j) + Pba(i,j)) * (+0.5d0,0.d0) * (Paa(j,i) - Pbb(j,i)) Sc_zx = Sc_zx - (+0.5d0,0.d0) * (Pab(j,i) + Pba(j,i)) * (+0.5d0,0.d0) * (Paa(i,j) - Pbb(i,j)) - enddo - enddo + end do + end do write(*,'(A15,2F10.6)') ' < Sx.Sz > = ',Sc_xz write(*,'(A15,2F10.6)') ' < Sz.Sx > = ',Sc_zx @@ -140,8 +140,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) do j = 1, nO Sc_yz = Sc_yz - (0.d0,-0.5d0) * (Pab(i,j) - Pba(i,j)) * (+0.5d0,0.d0) * (Paa(j,i) - Pbb(j,i)) Sc_zy = Sc_zy - (0.d0,-0.5d0) * (Pab(j,i) - Pba(j,i)) * (+0.5d0,0.d0) * (Paa(i,j) - Pbb(i,j)) - enddo - enddo + end do + end do write(*,'(A15,2F10.6)') ' < Sy.Sz > = ',Sc_yz write(*,'(A15,2F10.6)') ' < Sz.Sy > = ', Sc_zy write(*,*) @@ -163,8 +163,8 @@ subroutine print_GHF_spin(nBas,nBas2,nO,C,S) Mc(2,2) = Mc(2,2) - 0.25d0 * (Pba(i,j) - Pab(i,j))**2 Mc(3,3) = Mc(3,3) - 0.25d0 * (Paa(i,j) - Pbb(i,j))**2 Mc(1,3) = Mc(1,3) - 0.25d0 * (Pab(i,j) + Pba(i,j))*(Paa(j,i) - Pbb(j,i)) - enddo - enddo + end do + end do Mc(3,1) = Mc(1,3) write(*,*) 'The collinearity matrix is' diff --git a/src/MP/GMP2.f90 b/src/MP/GMP2.f90 index 7c95a02..72f4628 100644 --- a/src/MP/GMP2.f90 +++ b/src/MP/GMP2.f90 @@ -94,10 +94,10 @@ subroutine GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF,EcMP2) E2xs2 = E2xs2 + num*fs2 E2xk = E2xk + num*fk - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcMP2 = E2d + E2x EcsMP2 = E2ds + E2xs diff --git a/src/MP/GMP3.f90 b/src/MP/GMP3.f90 index 70ced8b..5691371 100644 --- a/src/MP/GMP3.f90 +++ b/src/MP/GMP3.f90 @@ -81,10 +81,10 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcMP2 = 0.25d0*E2 @@ -104,12 +104,12 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) E3a = E3a + OOVV(i,j,a,b)*OOOO(k,l,i,j)*VVOO(a,b,k,l)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do E3b = 0d0 @@ -125,12 +125,12 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) E3b = E3b + OOVV(i,j,a,b)*VVVV(a,b,c,d)*VVOO(c,d,i,j)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do E3c = 0d0 @@ -146,12 +146,12 @@ subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) E3c = E3c + OOVV(i,j,a,b)*OVVO(k,b,c,j)*VVOO(a,c,i,k)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c diff --git a/src/MP/RMP2.f90 b/src/MP/RMP2.f90 index e16fb10..1d5f4a6 100644 --- a/src/MP/RMP2.f90 +++ b/src/MP/RMP2.f90 @@ -91,10 +91,10 @@ subroutine RMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF,EcMP2) E2xs2 = E2xs2 - ERI(i,j,a,b)*ERI(i,j,b,a)*fs2 E2xk = E2xk - ERI(i,j,a,b)*ERI(i,j,b,a)*fk - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcMP2 = 2d0*E2d - E2x EcsMP2 = 2d0*E2ds - E2xs diff --git a/src/MP/RMP3.f90 b/src/MP/RMP3.f90 index fc7b817..788737b 100644 --- a/src/MP/RMP3.f90 +++ b/src/MP/RMP3.f90 @@ -105,10 +105,10 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e) E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do EcMP2 = 0.25d0*E2 @@ -128,12 +128,12 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e) E3a = E3a + OOVV(i,j,a,b)*OOOO(k,l,i,j)*VVOO(a,b,k,l)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do E3b = 0d0 @@ -149,12 +149,12 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e) E3b = E3b + OOVV(i,j,a,b)*VVVV(a,b,c,d)*VVOO(c,d,i,j)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do E3c = 0d0 @@ -170,12 +170,12 @@ subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e) E3c = E3c + OOVV(i,j,a,b)*OVVO(k,b,c,j)*VVOO(a,c,i,k)/(eps1*eps2) - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c diff --git a/src/MP/UMP2.f90 b/src/MP/UMP2.f90 index fee5982..98f4084 100644 --- a/src/MP/UMP2.f90 +++ b/src/MP/UMP2.f90 @@ -67,10 +67,10 @@ subroutine UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EUHF,eHF,Ec) Exaa = Exaa - 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,b,a)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do Ecaa = Edaa + Exaa Ec(1) = Ecaa @@ -93,10 +93,10 @@ subroutine UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EUHF,eHF,Ec) Edab = Edab + ERI_ab(i,j,a,b)*ERI_ab(i,j,a,b)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do Ecab = Edab + Exab Ec(2) = Ecab @@ -121,10 +121,10 @@ subroutine UMP2(dotest,nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EUHF,eHF,Ec) Exbb = Exbb - 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,b,a)/eps - enddo - enddo - enddo - enddo + end do + end do + end do + end do Ecbb = Edbb + Exbb Ec(3) = Ecbb diff --git a/src/RPA/crRRPA.f90 b/src/RPA/crRRPA.f90 index 2dc09aa..01febce 100644 --- a/src/RPA/crRRPA.f90 +++ b/src/RPA/crRRPA.f90 @@ -76,7 +76,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, call print_excitation_energies('crRPA@RHF','singlet',nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if ! Triplet manifold @@ -91,7 +91,7 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, call print_excitation_energies('crRPA@RHF','triplet',nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if if(exchange_kernel) then diff --git a/src/RPA/phACFDT_correlation_energy.f90 b/src/RPA/phACFDT_correlation_energy.f90 index 48cfaa5..16955eb 100644 --- a/src/RPA/phACFDT_correlation_energy.f90 +++ b/src/RPA/phACFDT_correlation_energy.f90 @@ -65,10 +65,10 @@ subroutine phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) - delta_Kx*ERI(i,j,b,a) - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Compute Tr(K x P_lambda) diff --git a/src/RPA/phRRPA.f90 b/src/RPA/phRRPA.f90 index c39eead..27a12ab 100644 --- a/src/RPA/phRRPA.f90 +++ b/src/RPA/phRRPA.f90 @@ -78,7 +78,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, call print_excitation_energies('phRPA@RHF','singlet',nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if ! Triplet manifold @@ -93,7 +93,7 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, call print_excitation_energies('phRPA@RHF','triplet',nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if if(exchange_kernel) then diff --git a/src/RPA/phRRPAx.f90 b/src/RPA/phRRPAx.f90 index f4414f7..2f27eb5 100644 --- a/src/RPA/phRRPAx.f90 +++ b/src/RPA/phRRPAx.f90 @@ -77,7 +77,7 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO call print_excitation_energies('phRPAx@RHF','singlet',nS,Om) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if ! Triplet manifold @@ -92,7 +92,7 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO call print_excitation_energies('phRPAx@RHF','triplet',nS,Om) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) - endif + end if if(exchange_kernel) then diff --git a/src/RPA/phURPA.f90 b/src/RPA/phURPA.f90 index ca5063e..9c9e991 100644 --- a/src/RPA/phURPA.f90 +++ b/src/RPA/phURPA.f90 @@ -94,7 +94,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB deallocate(Aph,Bph,Om,XpY,XmY) - endif + end if ! Spin-flip transitions @@ -119,7 +119,7 @@ subroutine phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nB deallocate(Aph,Bph,Om,XpY,XmY) - endif + end if if(exchange_kernel) then diff --git a/src/RPA/phURPAx.f90 b/src/RPA/phURPAx.f90 index 8dffc71..df18238 100644 --- a/src/RPA/phURPAx.f90 +++ b/src/RPA/phURPAx.f90 @@ -92,7 +92,7 @@ subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,n deallocate(Aph,Bph,Om,XpY,XmY) - endif + end if ! Spin-flip transitions @@ -117,7 +117,7 @@ subroutine phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,n deallocate(Aph,Bph,Om,XpY,XmY) - endif + end if if(exchange_kernel) then diff --git a/src/RPA/ppRRPA.f90 b/src/RPA/ppRRPA.f90 index 9c6a330..fb5c3c4 100644 --- a/src/RPA/ppRRPA.f90 +++ b/src/RPA/ppRRPA.f90 @@ -83,7 +83,7 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF, deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp) - endif + end if ! Triplet manifold @@ -115,7 +115,7 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF, deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp) - endif + end if EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/utils/elements.f90 b/src/utils/elements.f90 index 5ecf30f..a8a35a2 100644 --- a/src/utils/elements.f90 +++ b/src/utils/elements.f90 @@ -72,9 +72,9 @@ function element_core(zval,zatom) else write(*,*) '!!! not imlemented in element_core !!!' stop - endif + end if - endif + end if end function element_core diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index 42ecb67..37d2b21 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -48,11 +48,11 @@ subroutine orthogonalization_matrix(nBas,S,X) write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i) - endif + end if Uval(i) = 1d0/sqrt(Uval(i)) - enddo + end do call ADAt(nBas,Uvec,Uval,X) @@ -75,9 +75,9 @@ subroutine orthogonalization_matrix(nBas,S,X) write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' - endif + end if - enddo + end do call AD(nBas,Uvec,Uval) X = Uvec @@ -96,13 +96,13 @@ subroutine orthogonalization_matrix(nBas,S,X) ! Uval(i) = 1d0/sqrt(Uval(i)) ! else ! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization' -! endif -! enddo +! end if +! end do ! ! call AD(nBas,Uvec,Uval) ! X = Uvec - endif + end if ! Print results @@ -114,6 +114,6 @@ subroutine orthogonalization_matrix(nBas,S,X) call matout(nBas,nBas,X) write(*,*) - endif + end if end subroutine orthogonalization_matrix diff --git a/src/utils/read_auxiliary_basis.f90 b/src/utils/read_auxiliary_basis.f90 index cf4f6ed..2819348 100644 --- a/src/utils/read_auxiliary_basis.f90 +++ b/src/utils/read_auxiliary_basis.f90 @@ -50,7 +50,7 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & nShell = nShell + 1 do k=1,3 CenterShell(nShell,k) = XYZAtoms(iAt,k) - enddo + end do ! Shell type and contraction degree @@ -76,7 +76,7 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & elseif(shelltype == "I") then TotAngMomShell(nShell) = 6 write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif + end if ! Read exponents and contraction coefficients @@ -84,10 +84,10 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & do k=1,Kshell(nShell) read(2,*) ExpShell(nShell,k),DShell(nShell,k) write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo + end do + end do write(*,'(A28)') '------------------' - enddo + end do ! Total number of shells @@ -124,7 +124,7 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & nShell = nShell + 1 do k=1,3 CenterShell(nShell,k) = XYZAtoms(iAt,k) - enddo + end do ! Shell type and contraction degree @@ -150,7 +150,7 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & elseif(shelltype == "I") then TotAngMomShell(nShell) = 6 write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif + end if ! Read exponents and contraction coefficients @@ -158,10 +158,10 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & do k=1,Kshell(nShell) read(3,*) ExpShell(nShell,k),DShell(nShell,k) write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo + end do + end do write(*,'(A28)') '------------------' - enddo + end do ! Total number of shells diff --git a/src/utils/read_basis.f90 b/src/utils/read_basis.f90 index 75a9f49..c2efad3 100644 --- a/src/utils/read_basis.f90 +++ b/src/utils/read_basis.f90 @@ -70,7 +70,7 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh do k=1,ncart CenterShell(nShell,k) = rNuc(iNuc,k) - enddo + end do ! Shell type and contraction degree @@ -130,14 +130,14 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh do k=1,Kshell(nShell) read(2,*) kk,ExpShell(nShell,k),DShell(nShell,k) write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo + end do min_exponent(iNuc,TotAngMomShell(nShell)+1) & = min(min_exponent(iNuc,TotAngMomShell(nShell)+1),minval(ExpShell(nShell,1:KShell(nShell)))) max_exponent(iNuc) = max(max_exponent(iNuc),maxval(ExpShell(nShell,:))) max_ang_mom(iNuc) = max(max_ang_mom(iNuc),TotAngMomShell(nShell)) - enddo + end do !------------------------------------------------------------------------ ! End loop over shells !------------------------------------------------------------------------ @@ -151,7 +151,7 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh ! print*,'minimum exponent for atom n. ',iNuc,' = ' ! print*,min_exponent(iNuc,1:max_ang_mom(iNuc)+1) - enddo + end do !------------------------------------------------------------------------ ! End loop over atoms !------------------------------------------------------------------------ @@ -171,7 +171,7 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh nBas = 0 do iShell=1,nShell nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 - enddo + end do write(*,'(A28)') '------------------' write(*,'(A28,1X,I16)') 'Number of basis functions',NBas diff --git a/src/utils/read_dipole_integrals.f90 b/src/utils/read_dipole_integrals.f90 index 03e5284..553b46c 100644 --- a/src/utils/read_dipole_integrals.f90 +++ b/src/utils/read_dipole_integrals.f90 @@ -33,21 +33,21 @@ subroutine read_dipole_integrals(nBas,R) read(21,*,end=21) mu,nu,Dip R(mu,nu,1) = Dip R(nu,mu,1) = Dip - enddo + end do 21 close(unit=21) do read(22,*,end=22) mu,nu,Dip R(mu,nu,2) = Dip R(nu,mu,2) = Dip - enddo + end do 22 close(unit=22) do read(23,*,end=23) mu,nu,Dip R(mu,nu,3) = Dip R(nu,mu,3) = Dip - enddo + end do 23 close(unit=23) ! Print results @@ -66,6 +66,6 @@ subroutine read_dipole_integrals(nBas,R) write(*,'(A28)') '(mu|z|nu) integrals' write(*,'(A28)') '----------------------' call matout(nBas,nBas,R(:,:,3)) - endif + end if end subroutine read_dipole_integrals diff --git a/src/utils/read_geometry.f90 b/src/utils/read_geometry.f90 index d2f27bd..5c2da13 100644 --- a/src/utils/read_geometry.f90 +++ b/src/utils/read_geometry.f90 @@ -38,7 +38,7 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) read(10,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3) write(11,'(A3,1X,3F16.10)') El,rNuc(i,1)*BoToAn,rNuc(i,2)*BoToAn,rNuc(i,3)*BoToAn ZNuc(i) = dble(element_number(El)) - enddo + end do ! Compute nuclear repulsion energy @@ -48,8 +48,8 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) do j=i+1,nNuc RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2 ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB)) - enddo - enddo + end do + end do ! Close file with geometry specification close(unit=10) @@ -63,7 +63,7 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) write(*,'(A28,1X,I16)') 'Atom n. ',i write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i) write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rNuc(i,j),j=1,ncart) - enddo + end do write(*,*) write(*,'(A28)') '------------------' write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index 02c91a7..3b2baeb 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -48,7 +48,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) read(8,*,end=8) mu,nu,Ov S(mu,nu) = Ov S(nu,mu) = Ov - enddo + end do 8 close(unit=8) ! Read kinetic integrals @@ -58,7 +58,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) read(9,*,end=9) mu,nu,Kin T(mu,nu) = Kin T(nu,mu) = Kin - enddo + end do 9 close(unit=9) ! Read nuclear integrals @@ -68,7 +68,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) read(10,*,end=10) mu,nu,Nuc V(mu,nu) = Nuc V(nu,mu) = Nuc - enddo + end do 10 close(unit=10) ! Define core Hamiltonian @@ -98,7 +98,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) G(nu,mu,si,la) = ERI ! <43|21> G(si,la,nu,mu) = ERI - enddo + end do 11 close(unit=11) @@ -125,9 +125,9 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) do la=1,nBas do si=1,nBas call matout(nBas,nBas,G(1,1,la,si)) - enddo - enddo + end do + end do write(*,*) - endif + end if end subroutine read_integrals diff --git a/src/utils/sort.f90 b/src/utils/sort.f90 index 87d3239..36793ae 100644 --- a/src/utils/sort.f90 +++ b/src/utils/sort.f90 @@ -33,7 +33,7 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) end do if (i >= j) then exit - endif + end if tmp = x(i) x(i) = x(j) x(j) = tmp @@ -42,22 +42,22 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) iorder(j) = itmp i=i+1 j=j-1 - enddo + end do if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then if (first < i-1) then call rec_quicksort(x, iorder, isize, first, i-1,level/2) - endif + end if if (j+1 < last) then call rec_quicksort(x, iorder, isize, j+1, last,level/2) - endif + end if else if (first < i-1) then call rec_quicksort(x, iorder, isize, first, i-1,level/2) - endif + end if if (j+1 < last) then call rec_quicksort(x, iorder, isize, j+1, last,level/2) - endif - endif + end if + end if end subroutine subroutine set_order(x,iorder,isize,jsize) diff --git a/src/utils/spatial_to_spin_ERI.f90 b/src/utils/spatial_to_spin_ERI.f90 index bdb2919..8b82d85 100644 --- a/src/utils/spatial_to_spin_ERI.f90 +++ b/src/utils/spatial_to_spin_ERI.f90 @@ -27,9 +27,9 @@ subroutine spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) * Kronecker_delta(mod(q,2),mod(s,2)) & * ERI((p+1)/2,(q+1)/2,(r+1)/2,(s+1)/2) - enddo - enddo - enddo - enddo + end do + end do + end do + end do end subroutine spatial_to_spin_ERI diff --git a/src/utils/spatial_to_spin_MO_energy.f90 b/src/utils/spatial_to_spin_MO_energy.f90 index 0a35a8c..158d6ad 100644 --- a/src/utils/spatial_to_spin_MO_energy.f90 +++ b/src/utils/spatial_to_spin_MO_energy.f90 @@ -21,7 +21,7 @@ subroutine spatial_to_spin_MO_energy(nBas,e,nBas2,se) se(p) = e((p+1)/2) - enddo + end do ! print*,'MO energies in spinorbital basis' ! call matout(nBas2,1,se) diff --git a/src/utils/spatial_to_spin_fock.f90 b/src/utils/spatial_to_spin_fock.f90 index 614c693..db694ab 100644 --- a/src/utils/spatial_to_spin_fock.f90 +++ b/src/utils/spatial_to_spin_fock.f90 @@ -23,7 +23,7 @@ subroutine spatial_to_spin_fock(nBas,F,nBas2,sF) sF(p,q) = Kronecker_delta(mod(p,2),mod(q,2))*F((p+1)/2,(q+1)/2) - enddo - enddo + end do + end do end subroutine spatial_to_spin_fock diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index cf341be..efd3087 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -17,7 +17,7 @@ function Kronecker_delta(i,j) result(delta) delta = 1d0 else delta = 0d0 - endif + end if end function @@ -40,7 +40,7 @@ function KroneckerDelta(i,j) result(delta) delta = 1 else delta = 0 - endif + end if end function @@ -133,13 +133,13 @@ subroutine matout(m,n,A) do i=1,m do j=ilower,iupper B(j-ilower+1) = A(i,j) - enddo + end do do j=1,num if(abs(B(j)) < small) B(j) = 0d0 - enddo + end do write(*,'(I7,10F15.8)') i,(B(j),j=1,num) - enddo - enddo + end do + end do end subroutine @@ -181,7 +181,7 @@ subroutine trace_vector(n,v,Tr) Tr = 0d0 do i=1,n Tr = Tr + v(i) - enddo + end do end subroutine @@ -208,7 +208,7 @@ function trace_matrix(n,A) result(Tr) Tr = 0d0 do i=1,n Tr = Tr + A(i,i) - enddo + end do end function @@ -255,7 +255,7 @@ subroutine identity_matrix(N,A) do i=1,N A(i,i) = 1d0 - enddo + end do end subroutine @@ -286,9 +286,9 @@ subroutine prepend(N,M,A,b) do i=1,N do j=M-1,1,-1 A(i,j+1) = A(i,j) - enddo + end do A(i,1) = b(i) - enddo + end do end subroutine @@ -315,9 +315,9 @@ subroutine append(N,M,A,b) do i=1,N do j=2,M A(i,j-1) = A(i,j) - enddo + end do A(i,M) = b(i) - enddo + end do end subroutine @@ -348,9 +348,9 @@ subroutine AtDA(N,A,D,B) do j=1,N do k=1,N B(i,k) = B(i,k) + A(j,i)*D(j)*A(j,k) - enddo - enddo - enddo + end do + end do + end do end subroutine @@ -381,9 +381,9 @@ subroutine ADAt(N,A,D,B) do j=1,N do k=1,N B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j) - enddo - enddo - enddo + end do + end do + end do end subroutine !------------------------------------------------------------------------ @@ -402,8 +402,8 @@ subroutine DA(N,D,A) do i=1,N do j=1,N A(i,j) = D(i)*A(i,j) - enddo - enddo + end do + end do end subroutine @@ -423,8 +423,8 @@ subroutine AD(N,A,D) do i=1,N do j=1,N A(i,j) = A(i,j)*D(j) - enddo - enddo + end do + end do end subroutine @@ -466,8 +466,8 @@ subroutine CalcTrAB(n,A,B,Tr) do i=1,n do j=1,n Tr = Tr + A(i,j)*B(j,i) - enddo - enddo + end do + end do end subroutine @@ -488,7 +488,7 @@ function EpsilonSwitch(i,j) result(delta) delta = 1 else delta = -1 - endif + end if end function @@ -534,7 +534,7 @@ subroutine CalcInv3(a,det) b(i,1) = a(i,1) b(i,2) = a(i,2) b(i,3) = a(i,3) - enddo + end do a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2) a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3) @@ -551,8 +551,8 @@ subroutine CalcInv3(a,det) do i=1,3 do j=1,3 a(i,j) = a(i,j)/det - enddo - enddo + end do + end do end subroutine @@ -584,7 +584,7 @@ subroutine CalcInv4(a,det) b(2,i) = a(2,i) b(3,i) = a(3,i) b(4,i) = a(4,i) - enddo + end do a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)) a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)) @@ -609,8 +609,8 @@ subroutine CalcInv4(a,det) do i=1,4 do j=1,4 a(i,j) = a(i,j)/det - enddo - enddo + end do + end do end subroutine @@ -621,7 +621,7 @@ subroutine wall_time(t) integer*8, save :: rate = 0 if (rate == 0) then CALL SYSTEM_CLOCK(count_rate=rate) - endif + end if CALL SYSTEM_CLOCK(count=c) t = dble(c)/dble(rate) end subroutine diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 3801928..d63b722 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -33,7 +33,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR) if(info /= 0) then print*,'Problem in diagonalize_general_matrix (dgeev)!!' - endif + end if end subroutine @@ -70,7 +70,7 @@ subroutine diagonalize_matrix(N,A,e) if(info /= 0) then print*,'Problem in diagonalize_matrix (dsyev)!!' - endif + end if end subroutine @@ -114,7 +114,7 @@ subroutine svd(N,A,U,D,Vt) if (info /= 0) then print *, info, ': SVD failed' stop - endif + end if end @@ -144,7 +144,7 @@ subroutine inverse_matrix(N,A,B) print*,info stop 'error in inverse (dgetrf)!!' - endif + end if call dgetri(N,B,N,ipiv,work,lwork,info) @@ -153,7 +153,7 @@ subroutine inverse_matrix(N,A,B) print *, info stop 'error in inverse (dgetri)!!' - endif + end if deallocate(ipiv,work) @@ -195,7 +195,7 @@ subroutine linear_solve(N,A,b,x,rcond) ! print *, info ! stop 'error in linear_solve (dsysvx)!!' -! endif +! end if end subroutine @@ -226,7 +226,7 @@ subroutine easy_linear_solve(N,A,b,x) print *, info stop 'error in linear_solve (dsysv)!!' - endif + end if end subroutine