From 75f483972528dba57355dcd27fd7c69dcbf29c3d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 22 Jul 2023 22:19:46 +0200 Subject: [PATCH] create CC module --- src/CC/CC.f90 | 170 ++++++++++++++++++++++++++++++++++++++++++++ src/HF/HF.f90 | 25 +++++-- src/HF/RHF.f90 | 7 +- src/HF/UHF.f90 | 9 +-- src/MP/MP.f90 | 25 ++++--- src/QuAcK/QuAcK.f90 | 130 +++------------------------------ 6 files changed, 216 insertions(+), 150 deletions(-) create mode 100644 src/CC/CC.f90 diff --git a/src/CC/CC.f90 b/src/CC/CC.f90 new file mode 100644 index 0000000..70460c0 --- /dev/null +++ b/src/CC/CC.f90 @@ -0,0 +1,170 @@ +subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + +! Coupled-cluster module + + implicit none + include 'parameters.h' + +! Input variables + + logical :: doCCD + logical :: dopCCD + logical :: doDCD + logical :: doCCSD + logical :: doCCSDT + logical :: do_drCCD + logical :: do_rCCD + logical :: do_crCCD + logical :: do_lCCD + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + 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) :: epsHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: start_CC ,end_CC ,t_CC + +!------------------------------------------------------------------------ +! Perform CCD calculation +!------------------------------------------------------------------------ + + if(doCCD) then + + call cpu_time(start_CC) + call CCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform DCD calculation +!------------------------------------------------------------------------ + + if(doDCD) then + + call cpu_time(start_CC) + call DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & + ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform CCSD or CCSD(T) calculation +!------------------------------------------------------------------------ + + if(doCCSDT) doCCSD = .true. + + if(doCCSD) then + + call cpu_time(start_CC) + call CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform direct ring CCD calculation +!------------------------------------------------------------------------ + + if(do_drCCD) then + + call cpu_time(start_CC) + call drCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform ring CCD calculation +!------------------------------------------------------------------------ + + if(do_rCCD) then + + call cpu_time(start_CC) + call rCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform crossed-ring CCD calculation +!------------------------------------------------------------------------ + + if(do_crCCD) then + + call cpu_time(start_CC) + call crCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform ladder CCD calculation +!------------------------------------------------------------------------ + + if(do_lCCD) then + + call cpu_time(start_CC) + call lCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform pair CCD calculation +!------------------------------------------------------------------------ + + if(dopCCD) then + + call cpu_time(start_CC) + call pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) + call cpu_time(end_CC) + + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' + write(*,*) + + end if + +end subroutine diff --git a/src/HF/HF.f90 b/src/HF/HF.f90 index baf5a21..4c0408c 100644 --- a/src/HF/HF.f90 +++ b/src/HF/HF.f90 @@ -1,5 +1,5 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, & - nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF) ! Hartree-Fock module @@ -36,6 +36,8 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues ! Local variables + double precision :: start_HF ,end_HF ,t_HF + integer :: nSCF integer :: n_diis double precision :: ET @@ -64,7 +66,6 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues double precision,intent(out) :: epsHF(nBas) double precision,intent(out) :: cHF(nBas,nBas) double precision,intent(out) :: PHF(nBas,nBas) - double precision,intent(out) :: vHF(nBas) double precision,intent(out) :: F(nBas,nBas) !------------------------------------------------------------------------ @@ -81,8 +82,14 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues stop end if + call wall_time(start_HF) call RHF(maxSCF,thresh,n_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds' + write(*,*) end if @@ -95,8 +102,14 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues ! Switch on the unrestricted flag unrestricted = .true. + call wall_time(start_HF) call UHF(maxSCF,thresh,n_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds' + write(*,*) end if @@ -115,7 +128,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues end if ! call RMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, & -! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) +! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF) end if @@ -129,7 +142,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues unrestricted = .true. ! call UMOM(maxSCF,thresh,n_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, & -! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF,vHF) +! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF) end if diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 19e9b60..034a383 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,5 +1,5 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,eps,c,P,Vx) + nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,eps,c,P) ! Perform restricted Hartree-Fock calculation @@ -57,7 +57,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc double precision,intent(out) :: eps(nBas) double precision,intent(out) :: c(nBas,nBas) double precision,intent(out) :: P(nBas,nBas) - double precision,intent(out) :: Vx(nBas) double precision,intent(out) :: F(nBas,nBas) ! Hello world @@ -200,8 +199,4 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call print_RHF(nBas,nO,eps,C,ENuc,ET,EV,EJ,EK,EHF,dipole) -! Compute Vx for post-HF calculations - - call mo_fock_exchange_potential(nBas,c,P,ERI,Vx) - end subroutine diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 2b3d949..dad1b69 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -1,5 +1,5 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx) + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P) ! Perform unrestricted Hartree-Fock calculation @@ -61,7 +61,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision,intent(out) :: e(nBas,nspin) double precision,intent(out) :: c(nBas,nBas,nspin) double precision,intent(out) :: P(nBas,nBas,nspin) - double precision,intent(out) :: Vx(nBas,nspin) ! Hello world @@ -252,10 +251,4 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) -! Compute Vx for post-HF calculations - - do ispin=1,nspin - call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin)) - end do - end subroutine diff --git a/src/MP/MP.f90 b/src/MP/MP.f90 index fabe963..23ffc09 100644 --- a/src/MP/MP.f90 +++ b/src/MP/MP.f90 @@ -27,6 +27,8 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, ! Local variables + double precision :: start_MP ,end_MP ,t_MP + ! Output variables !------------------------------------------------------------------------ @@ -34,16 +36,18 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, !------------------------------------------------------------------------ if(doMP2) then - + + call cpu_time(start_MP) if(unrestricted) then - call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF) - else - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - end if + call cpu_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 @@ -53,16 +57,19 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, if(doMP3) then - if(unrestricted) then + call cpu_time(start_MP) + if(unrestricted) then write(*,*) 'MP3 NYI for UHF reference' stop - else - call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - end if + call cpu_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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index dbd8600..275b87a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -8,7 +8,7 @@ program QuAcK logical :: dostab logical :: doKS logical :: doMP,doMP2,doMP3 - logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical :: doCC,doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD logical :: doCIS,doCIS_D,doCID,doCISD,doFCI logical :: dophRPA,dophRPAx,docrRPA,doppRPA @@ -369,138 +369,26 @@ program QuAcK call cpu_time(end_MP) t_MP = end_MP - start_MP - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds' write(*,*) end if !------------------------------------------------------------------------ -! Perform CCD calculation +! Coupled-cluster module !------------------------------------------------------------------------ - if(doCCD) then + doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT + + if(doCC) then call cpu_time(start_CC) - call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call cpu_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform DCD calculation -!------------------------------------------------------------------------ - - if(doDCD) then - - call cpu_time(start_CC) - call DCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform CCSD or CCSD(T) calculation -!------------------------------------------------------------------------ - - if(doCCSDT) doCCSD = .true. - - if(doCCSD) then - - call cpu_time(start_CC) - call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform direct ring CCD calculation -!------------------------------------------------------------------------ - - if(do_drCCD) then - - call cpu_time(start_CC) - call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform ring CCD calculation -!------------------------------------------------------------------------ - - if(do_rCCD) then - - call cpu_time(start_CC) - call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform crossed-ring CCD calculation -!------------------------------------------------------------------------ - - if(do_crCCD) then - - call cpu_time(start_CC) - call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform ladder CCD calculation -!------------------------------------------------------------------------ - - if(do_lCCD) then - - call cpu_time(start_CC) - call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Perform pair CCD calculation -!------------------------------------------------------------------------ - - if(dopCCD) then - - call cpu_time(start_CC) - call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) - call cpu_time(end_CC) - - t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds' write(*,*) end if