From a60aa8d8ee795f74b13e0a2aab6b5086a375f0d9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 27 Oct 2023 13:35:10 +0200 Subject: [PATCH] GMP3 and GG0F2 --- input/methods | 4 +- src/GF/GG0F2.f90 | 120 ++++++++++++++++++++++ src/GF/GGF.f90 | 109 ++++++++++++++++++++ src/GF/GGF2_QP_graph.f90 | 78 ++++++++++++++ src/GF/GGF2_SigC.f90 | 50 +++++++++ src/GF/GGF2_dSigC.f90 | 52 ++++++++++ src/GF/GGF2_self_energy_diag.f90 | 72 +++++++++++++ src/HF/print_GHF.f90 | 41 ++++---- src/MP/GMP.f90 | 22 +++- src/MP/GMP2.f90 | 5 +- src/MP/GMP3.f90 | 170 +++++++++++++++++++++++++++++++ src/MP/MP2.f90 | 5 +- src/MP/RMP.f90 | 3 +- src/MP/UMP.f90 | 3 +- src/MP/UMP2.f90 | 3 +- src/QuAcK/GQuAcK.f90 | 22 ++-- src/QuAcK/QuAcK.f90 | 2 +- src/RPA/GRPA.f90 | 2 +- src/RPA/ppGRPA.f90 | 100 ++++++++++++++++++ 19 files changed, 815 insertions(+), 48 deletions(-) create mode 100644 src/GF/GG0F2.f90 create mode 100644 src/GF/GGF.f90 create mode 100644 src/GF/GGF2_QP_graph.f90 create mode 100644 src/GF/GGF2_SigC.f90 create mode 100644 src/GF/GGF2_dSigC.f90 create mode 100644 src/GF/GGF2_self_energy_diag.f90 create mode 100644 src/MP/GMP3.f90 create mode 100644 src/RPA/ppGRPA.f90 diff --git a/input/methods b/input/methods index a90066f..7f0c8cf 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF GHF ROHF F F T F # MP2 MP3 - F F + T T # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD @@ -13,6 +13,6 @@ # G0F2 evGF2 qsGF2 G0F3 evGF3 F F F F F # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW - T F F F F F + F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh F F F F F F diff --git a/src/GF/GG0F2.f90 b/src/GF/GG0F2.f90 new file mode 100644 index 0000000..eceb702 --- /dev/null +++ b/src/GF/GG0F2.f90 @@ -0,0 +1,120 @@ +subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + +! Perform a one-shot second-order Green function calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nC + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + double precision :: Ec + double precision :: EcBSE(nspin) + double precision,allocatable :: eGFlin(:) + double precision,allocatable :: eGF(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot second-order Green function |' + write(*,*)'************************************************' + write(*,*) + +! Memory allocation + + allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas)) + +! Frequency-dependent second-order contribution + + if(regularize) then + +! call GF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z) + + else + + call GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eHF,ERI,SigC,Z) + + end if + + eGFlin(:) = eHF(:) + Z(:)*SigC(:) + + if(linearize) then + + write(*,*) '*** Quasiparticle energies obtained by linearization ***' + + eGF(:) = eGFlin(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z) + + end if + + ! Print results + + call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) + call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) + +! Perform BSE2 calculation + +! if(dophBSE) then +! +! call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:)) +! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:)) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! Perform ppBSE2 calculation + +! if(doppBSE) then +! +! call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine diff --git a/src/GF/GGF.f90 b/src/GF/GGF.f90 new file mode 100644 index 0000000..9c63e33 --- /dev/null +++ b/src/GF/GGF.f90 @@ -0,0 +1,109 @@ +subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & + nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) + +! Green's function module + + implicit none + include 'parameters.h' + +! Input variables + + logical :: doG0F2 + logical :: doevGF2 + logical :: doqsGF2 + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + + double precision,intent(in) :: EHF + double precision,intent(in) :: epsHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + double precision :: start_GF ,end_GF ,t_GF + +!------------------------------------------------------------------------ +! Compute G0F2 electronic binding energies +!------------------------------------------------------------------------ + + if(doG0F2) then + + call wall_time(start_GF) + call GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + call wall_time(end_GF) + + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute evGF2 electronic binding energies +!------------------------------------------------------------------------ + + if(doevGF2) then + + call wall_time(start_GF) +! call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & +! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & +! ERI,dipole_int,epsHF) + call wall_time(end_GF) + + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform qsGF2 calculation +!------------------------------------------------------------------------ + + if(doqsGF2) then + + call wall_time(start_GF) +! call qsGGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & +! nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) + call wall_time(end_GF) + + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds' + write(*,*) + + end if + +end subroutine diff --git a/src/GF/GGF2_QP_graph.f90 b/src/GF/GGF2_QP_graph.f90 new file mode 100644 index 0000000..0a4fb8e --- /dev/null +++ b/src/GF/GGF2_QP_graph.f90 @@ -0,0 +1,78 @@ +subroutine GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z) + +! Compute the graphical solution of the GF2 QP equation + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGFlin(nBas) + double precision,intent(in) :: eOld(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision,external :: GGF2_SigC,GGF2_dSigC + double precision :: SigC,dSigC + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGF(nBas) + double precision,intent(out) :: Z(nBas) + +! Run Newton's algorithm to find the root + + write(*,*)'-----------------------------------------------------' + write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GFlin (eV)','e_GF (eV)','Z' + write(*,*)'-----------------------------------------------------' + + do p=nC+1,nBas-nR + + w = eGFlin(p) + nIt = 0 + f = 1d0 + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + SigC = GGF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI) + dSigC = GGF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI) + f = w - eHF(p) - SigC + df = 1d0/(1d0 - dSigC) + + w = w - df*f + + end do + + if(nIt == maxIt) then + + eGF(p) = eGFlin(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGFlin(p)*HaToeV,eGF(p)*HaToeV,Z(p),'Cvg Failed!' + + else + + eGF(p) = w + Z(p) = df + + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGFlin(p)*HaToeV,eGF(p)*HaToeV,Z(p) + + end if + + end do + +end subroutine diff --git a/src/GF/GGF2_SigC.f90 b/src/GF/GGF2_SigC.f90 new file mode 100644 index 0000000..82bdd2a --- /dev/null +++ b/src/GF/GGF2_SigC.f90 @@ -0,0 +1,50 @@ +double precision function GGF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + + GGF2_SigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = w + eHF(a) - eHF(i) - eHF(j) + GGF2_SigC = GGF2_SigC + 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2*eps/(eps**2 + eta**2) + + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = w + eHF(i) - eHF(a) - eHF(b) + GGF2_SigC = GGF2_SigC + 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2*eps/(eps**2 + eta**2) + + end do + end do + end do + +end function diff --git a/src/GF/GGF2_dSigC.f90 b/src/GF/GGF2_dSigC.f90 new file mode 100644 index 0000000..085842c --- /dev/null +++ b/src/GF/GGF2_dSigC.f90 @@ -0,0 +1,52 @@ +double precision function GGF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps + +! Initialize + + GGF2_dSigC = 0d0 + + ! Occupied part of the correlation self-energy + + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = w + eHF(a) - eHF(i) - eHF(j) + GGF2_dSigC = GGF2_dSigC - 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = w + eHF(i) - eHF(a) - eHF(b) + GGF2_dSigC = GGF2_dSigC - 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + +end function diff --git a/src/GF/GGF2_self_energy_diag.f90 b/src/GF/GGF2_self_energy_diag.f90 new file mode 100644 index 0000000..0fb610c --- /dev/null +++ b/src/GF/GGF2_self_energy_diag.f90 @@ -0,0 +1,72 @@ +subroutine GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z) + +! Compute diagonal part of the GF2 self-energy and its renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + integer :: p + double precision :: eps + double precision :: num + +! Output variables + + double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: Z(nBas) + +! Initialize + + SigC(:) = 0d0 + Z(:) = 0d0 + +! Compute GF2 self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = e(p) + e(a) - e(i) - e(j) + num = 0.5d0*(ERI(p,a,i,j) - ERI(p,a,j,i))**2 + + SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + + do p=nC+1,nBas-nR + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = e(p) + e(i) - e(a) - e(b) + num = 0.5d0*(ERI(p,i,a,b) - ERI(p,i,b,a))**2 + + SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index fa4b92f..9490602 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -24,7 +24,7 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) ! Local variables integer :: ixyz - integer :: i,j + integer :: mu,nu integer :: HOMO integer :: LUMO double precision :: Gap @@ -49,42 +49,37 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) allocate(Paa(nBas2,nBas2),Pab(nBas2,nBas2),Pba(nBas2,nBas2),Pbb(nBas2,nBas2)) -! Paa(:,:) = P( 1:nBas , 1:nBas ) -! Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) -! Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) -! Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) + Paa(:,:) = P( 1:nBas , 1:nBas ) + Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) + Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) + Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) allocate(Ca(nBas,nBas2),Cb(nBas,nBas2)) Ca(:,:) = C( 1:nBas ,1:nBas2) Cb(:,:) = C(nBas+1:nBas2,1:nBas2) - Paa = matmul(transpose(Ca),matmul(P( 1:nBas , 1:nBas ),Ca)) - Pab = matmul(transpose(Ca),matmul(P( 1:nBas ,nBas+1:nBas2),Cb)) - Pba = matmul(transpose(Cb),matmul(P(nBas+1:nBas2, 1:nBas ),Ca)) - Pbb = matmul(transpose(Cb),matmul(P(nBas+1:nBas2,nBas+1:nBas2),Cb)) - ! Compute expectation values of S^2 - Sx2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab+Pba)**2 - do i=1,nBas2 - do j=1,nBas2 - Sx2 = Sx2 - 0.5d0*(Paa(i,j)*Pbb(j,i) + Pab(i,j)*Pab(j,i)) + Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 + do mu=1,nBas + do nu=1,nBas + Sx2 = Sx2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) + Pab(mu,nu)*Pab(nu,mu)) end do end do - Sy2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) - 0.25d0*trace_matrix(nBas2,Pab+Pba)**2 - do i=1,nBas2 - do j=1,nBas2 - Sy2 = Sy2 - 0.5d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i)) + Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2 + do mu=1,nBas + do nu=1,nBas + Sy2 = Sy2 - 0.5d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) end do end do - Sz2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab-Pba)**2 - do i=1,nBas2 - do j=1,nBas2 - Sz2 = Sz2 - 0.25d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i)) - Sz2 = Sz2 + 0.25d0*(Pab(i,j)*Pba(j,i) - Pba(i,j)*Pab(j,i)) + Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2 + do mu=1,nBas + do nu=1,nBas + Sz2 = Sz2 - 0.25d0*(Paa(mu,nu)*Pbb(nu,mu) - Pab(mu,nu)*Pab(nu,mu)) + Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu)) end do end do diff --git a/src/MP/GMP.f90 b/src/MP/GMP.f90 index 2fe338d..eeb239e 100644 --- a/src/MP/GMP.f90 +++ b/src/MP/GMP.f90 @@ -1,4 +1,4 @@ -subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) +subroutine GMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) ! Moller-Plesset module @@ -8,6 +8,7 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) ! Input variables logical,intent(in) :: doMP2 + logical,intent(in) :: doMP3 logical,intent(in) :: regularize integer,intent(in) :: nBas @@ -23,6 +24,7 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) ! Local variables double precision :: start_MP ,end_MP ,t_MP + double precision :: Ec ! Output variables @@ -33,7 +35,23 @@ subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) if(doMP2) then call wall_time(start_MP) - call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec) + call wall_time(end_MP) + + t_MP = end_MP - start_MP + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute MP3 energy +!------------------------------------------------------------------------ + + if(doMP3) then + + call wall_time(start_MP) + call GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) call wall_time(end_MP) t_MP = end_MP - start_MP diff --git a/src/MP/GMP2.f90 b/src/MP/GMP2.f90 index ee00ed8..48be3b7 100644 --- a/src/MP/GMP2.f90 +++ b/src/MP/GMP2.f90 @@ -1,4 +1,4 @@ -subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) +subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Perform second-order Moller-Plesset calculation with and without regularizers @@ -30,10 +30,11 @@ subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) double precision :: E2x,E2xs,E2xs2,E2xk double precision :: EcsMP2,Ecs2MP2,EckMP2 - double precision :: EcMP2 ! Output variables + double precision,intent(out) :: EcMP2 + ! Hello world write(*,*) diff --git a/src/MP/GMP3.f90 b/src/MP/GMP3.f90 new file mode 100644 index 0000000..5f2e494 --- /dev/null +++ b/src/MP/GMP3.f90 @@ -0,0 +1,170 @@ +subroutine GMP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) + +! Perform third-order Moller-Plesset calculation + + implicit none + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc,EHF + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eps,E2,EcMP2 + double precision :: eps1,eps2,E3a,E3b,E3c + double precision :: EcMP3 + + integer :: i,j,k,l,a,b,c,d + + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VVOO(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Moller-Plesset third-order calculation |' + write(*,*)'************************************************' + write(*,*) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,ERI,dbERI) + +! Form energy denominator + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVOO(nV,nV,nO,nO),VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) + VVOO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas, 1:nO , 1:nO ) + VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas) + + deallocate(dbERI) + + allocate(eO(nO),eV(nV)) + + eO(:) = e(1:nO) + eV(:) = e(nO+1:nBas) + +! Compute MP2 energy + + E2 = 0d0 + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + eps = eO(i) + eO(j) - eV(a) - eV(b) + + E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps + + enddo + enddo + enddo + enddo + + EcMP2 = 0.25d0*E2 + +! Compute MP3 energy + + E3a = 0d0 + + do i=nC+1,nO + do j=nC+1,nO + do k=nC+1,nO + do l=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + eps1 = eO(i) + eO(j) - eV(a) - eV(b) + eps2 = eO(k) + eO(l) - eV(a) - eV(b) + + 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 + + E3b = 0d0 + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + do c=1,nV-nR + do d=1,nV-nR + + eps1 = eO(i) + eO(j) - eV(a) - eV(b) + eps2 = eO(i) + eO(j) - eV(c) - eV(d) + + 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 + + E3c = 0d0 + + do i=nC+1,nO + do j=nC+1,nO + do k=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + do c=1,nV-nR + + eps1 = eO(i) + eO(j) - eV(a) - eV(b) + eps2 = eO(i) + eO(k) - eV(a) - eV(c) + + 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 + + EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c + + write(*,*) + write(*,'(A32)') '-----------------------' + write(*,'(A32)') ' MP3 calculation ' + write(*,'(A32)') '-----------------------' + write(*,'(A32,1X,F16.10)') ' MP2 contribution ',EcMP2 + write(*,'(A32,1X,F16.10)') ' MP3 contribution ',EcMP3 + write(*,'(A32)') '-----------------------' + write(*,'(A32,1X,F16.10)') ' MP3 correlation energy', EcMP2 + EcMP3 + write(*,'(A32,1X,F16.10)') ' MP3 total energy',ENuc + EHF + EcMP2 + EcMP3 + write(*,'(A32)') '-----------------------' + write(*,*) + +end subroutine diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 1c431db..baf714c 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -1,4 +1,4 @@ -subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) +subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Perform second-order Moller-Plesset calculation with and without regularizers @@ -29,10 +29,11 @@ subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e) double precision :: E2x,E2xs,E2xs2,E2xk double precision :: EcsMP2,Ecs2MP2,EckMP2 - double precision :: EcMP2 ! Output variables + double precision,intent(out) :: EcMP2 + ! Hello world write(*,*) diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90 index 032ba9f..e048efd 100644 --- a/src/MP/RMP.f90 +++ b/src/MP/RMP.f90 @@ -24,6 +24,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) ! Local variables double precision :: start_MP ,end_MP ,t_MP + double precision :: Ec ! Output variables @@ -34,7 +35,7 @@ subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) if(doMP2) then call wall_time(start_MP) - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,Ec) call wall_time(end_MP) t_MP = end_MP - start_MP diff --git a/src/MP/UMP.f90 b/src/MP/UMP.f90 index 14edae4..6d1dc46 100644 --- a/src/MP/UMP.f90 +++ b/src/MP/UMP.f90 @@ -26,6 +26,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb ! Local variables double precision :: start_MP ,end_MP ,t_MP + double precision :: Ec(nsp) ! Output variables @@ -36,7 +37,7 @@ subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbb if(doMP2) then call wall_time(start_MP) - call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF) + call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF,Ec) call wall_time(end_MP) t_MP = end_MP - start_MP diff --git a/src/MP/UMP2.f90 b/src/MP/UMP2.f90 index 2eb2fd1..d7a361c 100644 --- a/src/MP/UMP2.f90 +++ b/src/MP/UMP2.f90 @@ -29,10 +29,11 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) double precision :: Edab,Exab,Ecab double precision :: Edbb,Exbb,Ecbb double precision :: Ed,Ex - double precision :: Ec(nsp) ! Output variables + double precision,intent(out) :: Ec(nsp) + ! Hello world write(*,*) diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 1e1e2ac..844f593 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -1,4 +1,4 @@ -subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & +subroutine GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & @@ -11,6 +11,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs logical,intent(in) :: doGHF logical,intent(in) :: dostab logical,intent(in) :: doMP2 + logical,intent(in) :: doMP3 logical,intent(in) :: dophRPA,dophRPAx,doppRPA logical,intent(in) :: doG0F2,doevGF2,doqsGF2 logical,intent(in) :: doG0W0,doevGW,doqsGW @@ -180,7 +181,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doMP) then call wall_time(start_MP) - call GMP(doMP2,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call GMP(doMP2,doMP3,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -198,8 +199,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doRPA) then call wall_time(start_RPA) - call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, & - nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S) + call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas2,nC,nO,nV,nR,nS,ENuc,EHF, & + ERI_MO,dipole_int_MO,epsHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -217,10 +218,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doGF) then call wall_time(start_GF) -! call GGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & -! dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, & -! nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & -! dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + call GGF(doG0F2,doevGF2,doqsGF2,maxSCF_GF,thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,lin_GF,eta_GF,reg_GF, & + nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -238,10 +237,9 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doGW) then call wall_time(start_GW) - call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS, & - EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) + call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,exchange_kernel,doXBS, & + dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc, & + nBas2,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 18b7005..acc97e9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -215,7 +215,7 @@ program QuAcK !--------------------------! if(doGQuAcK) & - call GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + call GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & diff --git a/src/RPA/GRPA.f90 b/src/RPA/GRPA.f90 index 88fda48..ce90a38 100644 --- a/src/RPA/GRPA.f90 +++ b/src/RPA/GRPA.f90 @@ -70,7 +70,7 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO, if(doppRPA) then call wall_time(start_RPA) -! call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,epsHF) + call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA diff --git a/src/RPA/ppGRPA.f90 b/src/RPA/ppGRPA.f90 new file mode 100644 index 0000000..45c2ad8 --- /dev/null +++ b/src/RPA/ppGRPA.f90 @@ -0,0 +1,100 @@ +subroutine ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,e) + +! Perform ppGRPA calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + integer :: nOO + integer :: nVV + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + double precision,allocatable :: Om1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: Om2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + + double precision :: EcRPA + +! Hello world + + write(*,*) + write(*,*)'****************************************' + write(*,*)'| particle-particle GRPA calculation |' + write(*,*)'****************************************' + write(*,*) + +! Initialization + + EcRPA = 0d0 + + ispin = 4 + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) + + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,e,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,e,ERI,Dpp) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) + +! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + call print_excitation_energies('ppRPA@GHF (N+2)',ispin,nVV,Om1) + call print_excitation_energies('ppRPA@GHF (N-2)',ispin,nOO,Om2) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',EcRPA + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EHF + EcRPA + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '--------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of ppRPA correlation energy' +! write(*,*) '--------------------------------------------------------' +! write(*,*) + +! call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcRPA) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcRPA(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcRPA(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcRPA(1) + EcRPA(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + EHF + EcRPA(1) + EcRPA(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine