diff --git a/input/methods b/input/methods index 3560ea4..36fcf41 100644 --- a/input/methods +++ b/input/methods @@ -16,6 +16,4 @@ T F F F F # G0T0 evGT qsGT F F F -# MCMP2 - F # * unrestricted version available diff --git a/input/options b/input/options index f9037df..9683845 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/AOtoMO/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 index a4cd5a3..7919084 100644 --- a/src/AOtoMO/AOtoMO_transform.f90 +++ b/src/AOtoMO/AOtoMO_transform.f90 @@ -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 diff --git a/src/AOtoMO/natural_orbital.f90 b/src/AOtoMO/natural_orbital.f90 deleted file mode 100644 index 14717a9..0000000 --- a/src/AOtoMO/natural_orbital.f90 +++ /dev/null @@ -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 diff --git a/src/GF/BSE2_dynamic_perturbation_iterative.f90 b/src/GF/BSE2_dynamic_perturbation_iterative.f90 index 9147bcd..5f72a5f 100644 --- a/src/GF/BSE2_dynamic_perturbation_iterative.f90 +++ b/src/GF/BSE2_dynamic_perturbation_iterative.f90 @@ -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 diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index ec35dc0..7040ee9 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d7ce9aa..9ddfc32 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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 !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 9b909e9..d025073 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -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) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 2478088..a67f579 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -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) diff --git a/src/utils/read_F12_integrals.f90 b/src/utils/read_F12_integrals.f90 new file mode 100644 index 0000000..93fe2e5 --- /dev/null +++ b/src/utils/read_F12_integrals.f90 @@ -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