LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-08-31 00:23:45 +02:00
Code Releases Activity

remove old methods and options

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-11-29 09:54:18 +01:00
parent 16fe319ef3
commit 2a46cc4617
10 changed files with 175 additions and 293 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -16,6 +16,4 @@
T F F F F
# G0T0 evGT qsGT
F F F
# MCMP2
F
# * unrestricted version available

2
input/options
View File

@ -16,5 +16,3 @@
F T T
# BSE: BSE dBSE dTDA evDyn ppBSE BSE2
T F T F F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

9
src/AOtoMO/AOtoMO_transform.f90
View File

@ -9,19 +9,10 @@ subroutine AOtoMO_transform(nBas,c,A)
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
integer :: i,j
! Output variables
double precision,intent(inout):: A(nBas,nBas)
A = matmul(transpose(c),matmul(A,c))
! do j=1,nBas
! do i=1,nBas
! write(10,'(I5,I5,F16.10)') i,j,A(i,j)
! enddo
! enddo
end subroutine AOtoMO_transform

57
src/AOtoMO/natural_orbital.f90
View File

@ -1,57 +0,0 @@
subroutine natural_orbital(nBas,nO,cHF,c)
! Compute natural orbitals and natural occupancies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cHF(nBas,nBas),c(nBas,nBas)
! Local variables
integer :: i,j,k
double precision,allocatable :: eNO(:),cNO(:,:),P(:,:)
! Allocate
allocate(eNO(nBas),cNO(nBas,nBas),P(nBas,nBas))
! Compute density matrix
P = matmul(transpose(cHF),cHF)
call matout(nBas,nBas,P)
cNO = 0d0
do i=1,nBas
do j=1,nBas
do k=1,1
cNO(i,j) = cNO(i,j) + 2d0*P(i,k)*P(j,k)
enddo
enddo
enddo
! cNO(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
! cNO = matmul(transpose(cHF),matmul(cNO,cHF))
call diagonalize_matrix(nBas,cNO,eNO)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Natural orbitals '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cNO)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Natural occupancies'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eNO)
write(*,*)
end subroutine natural_orbital

2
src/GF/BSE2_dynamic_perturbation_iterative.f90
View File

@ -120,7 +120,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) &
- dot_product(Y,matmul(Am_dyn - A_sta,Y)) &
+ dot_product(X,matmul(B_dyn - B_sta,Y)) &
+ dot_product(X,matmul(B_dyn - B_sta,Y)) &
- dot_product(Y,matmul(B_dyn - B_sta,X))
end if

8
src/GW/qsGW.f90
View File

@ -301,14 +301,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
! Compute the overlap between HF and GW orbitals
! call overlap(nBas,cHF,c)
! Compute natural orbitals and occupancies
! call natural_orbital(nBas,nO,cHF,c)
! Did it actually converge?
if(nSCF == maxSCF+1) then

180
src/QuAcK/QuAcK.f90
View File

@ -17,9 +17,6 @@ program QuAcK
logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical :: doG0T0,doevGT,doqsGT
logical :: doMCMP2,doMinMCMP2
logical :: doGTGW = .false.
logical :: doBas
integer :: nNuc,nBas,nBasCABS
integer :: nEl(nspin)
@ -29,7 +26,7 @@ program QuAcK
integer :: nR(nspin)
integer :: nS(nspin)
double precision :: ENuc,ERHF,EUHF,Norm
double precision :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
double precision :: EcMP2(3),EcMP3,EcMP2F12(3)
double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
@ -106,9 +103,6 @@ program QuAcK
double precision :: start_MP2 ,end_MP2 ,t_MP2
double precision :: start_MP3 ,end_MP3 ,t_MP3
double precision :: start_MP2F12 ,end_MP2F12 ,t_MP2F12
double precision :: start_MCMP2 ,end_MCMP2 ,t_MCMP2
double precision :: start_MinMCMP2,end_MinMCMP2,t_MinMCMP2
double precision :: start_Bas ,end_Bas ,t_Bas
integer :: maxSCF_HF,n_diis_HF
double precision :: thresh_HF,level_shift
@ -175,8 +169,7 @@ program QuAcK
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
doG0T0,doevGT,doqsGT)
! Read options for methods
@ -188,12 +181,7 @@ program QuAcK
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn,ppBSE,BSE2, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Weird stuff
doMinMCMP2 = .false.
BSE,dBSE,dTDA,evDyn,ppBSE,BSE2)
!------------------------------------------------------------------------
! Read input information
@ -227,20 +215,6 @@ program QuAcK
max_ang_mom,min_exponent,max_exponent)
nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:))
!------------------------------------------------------------------------
! Read auxiliary basis set information
!------------------------------------------------------------------------
! call ReadAuxBasis(nNuc,rNuc,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell)
! Compute the number of basis functions
! call CalcNBasis(nShell,TotAngMomShell,nA)
! Number of virtual orbitals in complete space
! nBasCABS = nA - nBas
!------------------------------------------------------------------------
! Read one- and two-electron integrals
!------------------------------------------------------------------------
@ -1240,154 +1214,6 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Information for Monte Carlo calculations
!------------------------------------------------------------------------
if(doMCMP2 .or. doMinMCMP2) then
! Print simulation details
write(*,'(A32)') '----------------------'
write(*,'(A32,1X,I16)') 'Number of Monte Carlo steps',nMC
write(*,'(A32,1X,I16)') 'Number of equilibration steps',nEq
write(*,'(A32,1X,I16)') 'Number of walkers',nWalk
write(*,'(A32,1X,F16.10)') 'Initial time step',dt
write(*,'(A32,1X,I16)') 'Frequency of ouput',nPrint
write(*,'(A32,1X,I16)') 'Seed for random number generator',iSeed
write(*,'(A32)') '----------------------'
write(*,*)
! Initialize random number generator
call initialize_random_generator(iSeed)
!------------------------------------------------------------------------
! Type of weight function
!------------------------------------------------------------------------
! TrialType = 0 => HF density
! TrialType = 1 => Custom one-electron function
!------------------------------------------------------------------------
TrialType = 0
allocate(cTrial(nBas),gradient(nBas),hessian(nBas,nBas))
end if
!------------------------------------------------------------------------
! Compute MC-MP2 energy
!------------------------------------------------------------------------
if(doMCMP2) then
call cpu_time(start_MCMP2)
! call MCMP2(doDrift,nBas,nC,nO,nV,cHF,eHF,EcMP2, &
! nMC,nEq,nWalk,dt,nPrint, &
! nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
! Norm,EcMCMP2,Err_EcMCMP2,Var_EcMCMP2)
call cpu_time(end_MCMP2)
t_MCMP2 = end_MCMP2 - start_MCMP2
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MC-MP2 = ',t_MCMP2,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Minimize MC-MP2 variance
!------------------------------------------------------------------------
if(doMinMCMP2) then
call cpu_time(start_MinMCMP2)
! call MinMCMP2(nBas,nEl,nC,nO,nV,cHF,eHF,EcMP2, &
! nMC,nEq,nWalk,dt,nPrint, &
! nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
! TrialType,Norm,cTrial,gradient,hessian)
call cpu_time(end_MinMCMP2)
t_MinMCMP2 = end_MinMCMP2 - start_MinMCMP2
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MC-MP2 variance minimization = ',t_MinMCMP2,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Range-separated GT/GW
!------------------------------------------------------------------------
! if(doGTGW) then
! ! Read and transform long-range two-electron integrals
! allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas))
! call read_LR(nBas,ERI_ERF_AO)
! call cpu_time(start_AOtoMO)
! write(*,*)
! write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient'
! write(*,*)
! call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO)
! call cpu_time(end_AOtoMO)
! deallocate(ERI_ERF_AO)
! t_AOtoMO = end_AOtoMO - start_AOtoMO
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
! write(*,*)
! ! Long-range G0W0 calculation
! call cpu_time(start_G0W0)
! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
! dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0)
! call cpu_time(end_G0W0)
!
! t_G0W0 = end_G0W0 - start_G0W0
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds'
! write(*,*)
! ! Short-range G0T0 calculation
! ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:)
! call cpu_time(start_G0T0)
! call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
! singlet,triplet,linGW,eta_GW, &
! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0T0)
! call cpu_time(end_G0T0)
!
! t_G0T0 = end_G0T0 - start_G0T0
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
! write(*,*)
! call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
! end if
!------------------------------------------------------------------------
! Basis set correction
!------------------------------------------------------------------------
doBas = .false.
if(doBas) then
call cpu_time(start_Bas)
call basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
ERI_MO,eHF,cHF,PHF,eG0W0)
call cpu_time(end_Bas)
t_Bas = end_Bas - start_Bas
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for basis set correction = ',t_Bas,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute FCI
!------------------------------------------------------------------------

12
src/QuAcK/read_methods.f90
View File

@ -8,8 +8,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
doG0T0,doevGT,doqsGT)
! Read desired methods
@ -26,7 +25,6 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(out) :: doG0T0,doevGT,doqsGT
logical,intent(out) :: doMCMP2
! Local variables
@ -85,8 +83,6 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doevGT = .false.
doqsGT = .false.
doMCMP2 = .false.
! Read mean-field methods
read(1,*)
@ -168,12 +164,6 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
if(answer2 == 'T') doevGT = .true.
if(answer3 == 'T') doqsGT = .true.
! Read stochastic methods
read(1,*)
read(1,*) answer1
if(answer1 == 'T') doMCMP2 = .true.
! Close file with geometry specification
close(unit=1)

27
src/QuAcK/read_options.f90
View File

@ -6,8 +6,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn,ppBSE,BSE2, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
BSE,dBSE,dTDA,evDyn,ppBSE,BSE2)
! Read desired methods
@ -78,14 +77,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
logical,intent(out) :: ppBSE
logical,intent(out) :: BSE2
integer,intent(out) :: nMC
integer,intent(out) :: nEq
integer,intent(out) :: nWalk
double precision,intent(out) :: dt
integer,intent(out) :: nPrint
integer,intent(out) :: iSeed
logical,intent(out) :: doDrift
! Local variables
character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8
@ -252,22 +243,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
if(answer5 == 'T') ppBSE = .true.
if(answer6 == 'T') BSE2 = .true.
! Read options for MC-MP2: Monte Carlo steps, number of equilibration steps, number of walkers,
! Monte Carlo time step, frequency of output results, and seed for random number generator
nMC = 100000
nEq = 10000
nWalk = 10
dt = 0.3d0
nPrint = 1000
iSeed = 0
doDrift = .false.
read(1,*)
read(1,*) nMC,nEq,nWalk,dt,nPrint,iSeed,answer1
if(answer1 == 'T') doDrift = .true.
! Close file with options
close(unit=1)

169
src/utils/read_F12_integrals.f90 Normal file
View File

@ -0,0 +1,169 @@
subroutine read_F12_integrals(nBas,S,C,F,Y,FC)
! Read one- and two-electron integrals from files
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
! Local variables
logical :: debug
integer :: mu,nu,la,si,ka,ta
double precision :: ERI,F12,Yuk,F13C12,ExpS
! Output variables
double precision,intent(out) :: C(nBas,nBas,nBas,nBas)
double precision,intent(out) :: F(nBas,nBas,nBas,nBas)
double precision,intent(out) :: Y(nBas,nBas,nBas,nBas)
double precision,intent(out) :: FC(nBas,nBas,nBas,nBas,nBas,nBas)
debug = .false.
! Open file with integrals
open(unit=21,file='int/ERI.dat')
open(unit=22,file='int/F12.dat')
open(unit=23,file='int/Yuk.dat')
open(unit=31,file='int/3eInt_Type1.dat')
! Read 1/r12 integrals
C = 0d0
do
read(21,*,end=21) mu,nu,la,si,ERI
! <12|34>
C(mu,nu,la,si) = ERI
! <32|14>
C(la,nu,mu,si) = ERI
! <14|32>
C(mu,si,la,nu) = ERI
! <34|12>
C(la,si,mu,nu) = ERI
! <41|23>
C(si,mu,nu,la) = ERI
! <23|41>
C(nu,la,si,mu) = ERI
! <21|43>
C(nu,mu,si,la) = ERI
! <43|21>
C(si,la,nu,mu) = ERI
enddo
21 close(unit=21)
! Read f12 integrals
F = 0d0
do
read(22,*,end=22) mu,nu,la,si,F12
! <12|34>
F(mu,nu,la,si) = F12
! <32|14>
F(la,nu,mu,si) = F12
! <14|32>
F(mu,si,la,nu) = F12
! <34|12>
F(la,si,mu,nu) = F12
! <41|23>
F(si,mu,nu,la) = F12
! <23|41>
F(nu,la,si,mu) = F12
! <21|43>
F(nu,mu,si,la) = F12
! <43|21>
F(si,la,nu,mu) = F12
enddo
22 close(unit=22)
! Read f12/r12 integrals
Y = 0d0
do
read(23,*,end=23) mu,nu,la,si,Yuk
! <12|34>
Y(mu,nu,la,si) = Yuk
! <32|14>
Y(la,nu,mu,si) = Yuk
! <14|32>
Y(mu,si,la,nu) = Yuk
! <34|12>
Y(la,si,mu,nu) = Yuk
! <41|23>
Y(si,mu,nu,la) = Yuk
! <23|41>
Y(nu,la,si,mu) = Yuk
! <21|43>
Y(nu,mu,si,la) = Yuk
! <43|21>
Y(si,la,nu,mu) = Yuk
enddo
23 close(unit=23)
! Read f13/r12 integrals
FC = 0d0
do
read(31,*,end=31) mu,nu,la,si,ka,ta,F13C12
FC(mu,nu,la,si,ka,ta) = F13C12
enddo
31 close(unit=31)
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,C(1,1,la,si))
enddo
enddo
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'F12 integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,F(1,1,la,si))
enddo
enddo
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Yukawa integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,Y(1,1,la,si))
enddo
enddo
write(*,*)
endif
! Read exponent of Slater geminal
open(unit=4,file='input/geminal')
read(4,*) ExpS
close(unit=4)
! Transform two-electron integrals
! do mu=1,nBas
! do nu=1,nBas
! do la=1,nBas
! do si=1,nBas
! F(mu,nu,la,si) = (S(mu,la)*S(nu,si) - F(mu,nu,la,si))/ExpS
! Y(mu,nu,la,si) = (C(mu,nu,la,si) - Y(mu,nu,la,si))/ExpS
! enddo
! enddo
! enddo
! enddo
end subroutine read_F12_integrals