diff --git a/input/basis b/input/basis deleted file mode 100644 index 79a747a..0000000 --- a/input/basis +++ /dev/null @@ -1,9 +0,0 @@ -1 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 diff --git a/input/molecule b/input/molecule deleted file mode 100644 index e2a4fd3..0000000 --- a/input/molecule +++ /dev/null @@ -1,5 +0,0 @@ -# nAt nEla nElb nCore nRyd - 2 4 3 0 0 -# Znuc x y z - C 0. 0. -0.16245872 - H 0. 0. 1.93436816 diff --git a/input/molecule.xyz b/input/molecule.xyz deleted file mode 100644 index 7a4f218..0000000 --- a/input/molecule.xyz +++ /dev/null @@ -1,4 +0,0 @@ - 2 - - C 0.0000000000 0.0000000000 -0.0859694585 - H 0.0000000000 0.0000000000 1.0236236215 diff --git a/input/options b/input/options index d9b23c7..de8da3b 100644 --- a/input/options +++ b/input/options @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/input/weight b/input/weight index fb05e68..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,18 +1,9 @@ 1 3 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 0.1220000 1.0000000 + 1 0.2976000 1.0000000 P 1 - 1 0.7270000 1.0000000 -2 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 + 1 1.2750000 1.0000000 diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..9681dd8 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,2 @@ +*. +*/Makefile diff --git a/src/QuAcK/ADC.f90 b/src/ADC/ADC.f90 similarity index 100% rename from src/QuAcK/ADC.f90 rename to src/ADC/ADC.f90 diff --git a/src/QuAcK/ADC2.f90 b/src/ADC/ADC2.f90 similarity index 100% rename from src/QuAcK/ADC2.f90 rename to src/ADC/ADC2.f90 diff --git a/src/ADC/Makefile b/src/ADC/Makefile new file mode 100644 index 0000000..291c256 --- /dev/null +++ b/src/ADC/Makefile @@ -0,0 +1,2 @@ +TARGET=ADC.a +include ../Makefile.include diff --git a/src/QuAcK/AOtoMO_integral_transform.f90 b/src/AOtoMO/AOtoMO_integral_transform.f90 similarity index 100% rename from src/QuAcK/AOtoMO_integral_transform.f90 rename to src/AOtoMO/AOtoMO_integral_transform.f90 diff --git a/src/QuAcK/AOtoMO_oooa.f90 b/src/AOtoMO/AOtoMO_oooa.f90 similarity index 100% rename from src/QuAcK/AOtoMO_oooa.f90 rename to src/AOtoMO/AOtoMO_oooa.f90 diff --git a/src/QuAcK/AOtoMO_oooo.f90 b/src/AOtoMO/AOtoMO_oooo.f90 similarity index 100% rename from src/QuAcK/AOtoMO_oooo.f90 rename to src/AOtoMO/AOtoMO_oooo.f90 diff --git a/src/QuAcK/AOtoMO_oooooo.f90 b/src/AOtoMO/AOtoMO_oooooo.f90 similarity index 100% rename from src/QuAcK/AOtoMO_oooooo.f90 rename to src/AOtoMO/AOtoMO_oooooo.f90 diff --git a/src/QuAcK/AOtoMO_oovv.f90 b/src/AOtoMO/AOtoMO_oovv.f90 similarity index 100% rename from src/QuAcK/AOtoMO_oovv.f90 rename to src/AOtoMO/AOtoMO_oovv.f90 diff --git a/src/QuAcK/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 similarity index 100% rename from src/QuAcK/AOtoMO_transform.f90 rename to src/AOtoMO/AOtoMO_transform.f90 diff --git a/src/QuAcK/Coulomb_matrix_AO_basis.f90 b/src/AOtoMO/Coulomb_matrix_AO_basis.f90 similarity index 100% rename from src/QuAcK/Coulomb_matrix_AO_basis.f90 rename to src/AOtoMO/Coulomb_matrix_AO_basis.f90 diff --git a/src/QuAcK/Coulomb_matrix_MO_basis.f90 b/src/AOtoMO/Coulomb_matrix_MO_basis.f90 similarity index 100% rename from src/QuAcK/Coulomb_matrix_MO_basis.f90 rename to src/AOtoMO/Coulomb_matrix_MO_basis.f90 diff --git a/src/QuAcK/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 similarity index 100% rename from src/QuAcK/Hartree_matrix_AO_basis.f90 rename to src/AOtoMO/Hartree_matrix_AO_basis.f90 diff --git a/src/QuAcK/Hartree_matrix_MO_basis.f90 b/src/AOtoMO/Hartree_matrix_MO_basis.f90 similarity index 100% rename from src/QuAcK/Hartree_matrix_MO_basis.f90 rename to src/AOtoMO/Hartree_matrix_MO_basis.f90 diff --git a/src/QuAcK/MOtoAO_transform.f90 b/src/AOtoMO/MOtoAO_transform.f90 similarity index 100% rename from src/QuAcK/MOtoAO_transform.f90 rename to src/AOtoMO/MOtoAO_transform.f90 diff --git a/src/AOtoMO/Makefile b/src/AOtoMO/Makefile new file mode 100644 index 0000000..030aef4 --- /dev/null +++ b/src/AOtoMO/Makefile @@ -0,0 +1,2 @@ +TARGET=AOtoMO.a +include ../Makefile.include diff --git a/src/QuAcK/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 similarity index 100% rename from src/QuAcK/exchange_matrix_AO_basis.f90 rename to src/AOtoMO/exchange_matrix_AO_basis.f90 diff --git a/src/QuAcK/natural_orbital.f90 b/src/AOtoMO/natural_orbital.f90 similarity index 100% rename from src/QuAcK/natural_orbital.f90 rename to src/AOtoMO/natural_orbital.f90 diff --git a/src/QuAcK/LSD_SR.f b/src/BasCor/LSD_SR.f similarity index 100% rename from src/QuAcK/LSD_SR.f rename to src/BasCor/LSD_SR.f diff --git a/src/BasCor/Makefile b/src/BasCor/Makefile new file mode 100644 index 0000000..cb32779 --- /dev/null +++ b/src/BasCor/Makefile @@ -0,0 +1,2 @@ +TARGET=BasCor.a +include ../Makefile.include diff --git a/src/QuAcK/basis_correction.f90 b/src/BasCor/basis_correction.f90 similarity index 100% rename from src/QuAcK/basis_correction.f90 rename to src/BasCor/basis_correction.f90 diff --git a/src/QuAcK/dft_grid.f b/src/BasCor/dft_grid.f similarity index 100% rename from src/QuAcK/dft_grid.f rename to src/BasCor/dft_grid.f diff --git a/src/QuAcK/ec_srlda.f90 b/src/BasCor/ec_srlda.f90 similarity index 100% rename from src/QuAcK/ec_srlda.f90 rename to src/BasCor/ec_srlda.f90 diff --git a/src/QuAcK/ecmd_lda.f90 b/src/BasCor/ecmd_lda.f90 similarity index 100% rename from src/QuAcK/ecmd_lda.f90 rename to src/BasCor/ecmd_lda.f90 diff --git a/src/QuAcK/f_grid.f90 b/src/BasCor/f_grid.f90 similarity index 100% rename from src/QuAcK/f_grid.f90 rename to src/BasCor/f_grid.f90 diff --git a/src/QuAcK/fc_srlda.f90 b/src/BasCor/fc_srlda.f90 similarity index 100% rename from src/QuAcK/fc_srlda.f90 rename to src/BasCor/fc_srlda.f90 diff --git a/src/QuAcK/mu_grid.f90 b/src/BasCor/mu_grid.f90 similarity index 100% rename from src/QuAcK/mu_grid.f90 rename to src/BasCor/mu_grid.f90 diff --git a/src/QuAcK/quadrature_grid.f90 b/src/BasCor/quadrature_grid.f90 similarity index 100% rename from src/QuAcK/quadrature_grid.f90 rename to src/BasCor/quadrature_grid.f90 diff --git a/src/QuAcK/read_F12_integrals.f90 b/src/BasCor/read_F12_integrals.f90 similarity index 100% rename from src/QuAcK/read_F12_integrals.f90 rename to src/BasCor/read_F12_integrals.f90 diff --git a/src/QuAcK/read_grid.f90 b/src/BasCor/read_grid.f90 similarity index 100% rename from src/QuAcK/read_grid.f90 rename to src/BasCor/read_grid.f90 diff --git a/src/QuAcK/CCD.f90 b/src/CC/CCD.f90 similarity index 100% rename from src/QuAcK/CCD.f90 rename to src/CC/CCD.f90 diff --git a/src/QuAcK/CCD_correlation_energy.f90 b/src/CC/CCD_correlation_energy.f90 similarity index 100% rename from src/QuAcK/CCD_correlation_energy.f90 rename to src/CC/CCD_correlation_energy.f90 diff --git a/src/QuAcK/CCSD.f90 b/src/CC/CCSD.f90 similarity index 100% rename from src/QuAcK/CCSD.f90 rename to src/CC/CCSD.f90 diff --git a/src/QuAcK/CCSDT.f90 b/src/CC/CCSDT.f90 similarity index 100% rename from src/QuAcK/CCSDT.f90 rename to src/CC/CCSDT.f90 diff --git a/src/QuAcK/CCSD_Ec_nc.f90 b/src/CC/CCSD_Ec_nc.f90 similarity index 100% rename from src/QuAcK/CCSD_Ec_nc.f90 rename to src/CC/CCSD_Ec_nc.f90 diff --git a/src/QuAcK/CCSD_correlation_energy.f90 b/src/CC/CCSD_correlation_energy.f90 similarity index 100% rename from src/QuAcK/CCSD_correlation_energy.f90 rename to src/CC/CCSD_correlation_energy.f90 diff --git a/src/QuAcK/MP3_correlation_energy.f90 b/src/CC/MP3_correlation_energy.f90 similarity index 100% rename from src/QuAcK/MP3_correlation_energy.f90 rename to src/CC/MP3_correlation_energy.f90 diff --git a/src/CC/Makefile b/src/CC/Makefile new file mode 100644 index 0000000..df8eebf --- /dev/null +++ b/src/CC/Makefile @@ -0,0 +1,2 @@ +TARGET=CC.a +include ../Makefile.include diff --git a/src/QuAcK/antisymmetrize_ERI.f90 b/src/CC/antisymmetrize_ERI.f90 similarity index 100% rename from src/QuAcK/antisymmetrize_ERI.f90 rename to src/CC/antisymmetrize_ERI.f90 diff --git a/src/QuAcK/drCCD.f90 b/src/CC/drCCD.f90 similarity index 100% rename from src/QuAcK/drCCD.f90 rename to src/CC/drCCD.f90 diff --git a/src/QuAcK/form_EOM_one_body.f90 b/src/CC/form_EOM_one_body.f90 similarity index 100% rename from src/QuAcK/form_EOM_one_body.f90 rename to src/CC/form_EOM_one_body.f90 diff --git a/src/QuAcK/form_EOM_two_body.f90 b/src/CC/form_EOM_two_body.f90 similarity index 100% rename from src/QuAcK/form_EOM_two_body.f90 rename to src/CC/form_EOM_two_body.f90 diff --git a/src/QuAcK/form_T.f90 b/src/CC/form_T.f90 similarity index 100% rename from src/QuAcK/form_T.f90 rename to src/CC/form_T.f90 diff --git a/src/QuAcK/form_X.f90 b/src/CC/form_X.f90 similarity index 100% rename from src/QuAcK/form_X.f90 rename to src/CC/form_X.f90 diff --git a/src/QuAcK/form_abh.f90 b/src/CC/form_abh.f90 similarity index 100% rename from src/QuAcK/form_abh.f90 rename to src/CC/form_abh.f90 diff --git a/src/QuAcK/form_cF_nc.f90 b/src/CC/form_cF_nc.f90 similarity index 100% rename from src/QuAcK/form_cF_nc.f90 rename to src/CC/form_cF_nc.f90 diff --git a/src/QuAcK/form_cW_nc.f90 b/src/CC/form_cW_nc.f90 similarity index 100% rename from src/QuAcK/form_cW_nc.f90 rename to src/CC/form_cW_nc.f90 diff --git a/src/QuAcK/form_delta_OOOVVV.f90 b/src/CC/form_delta_OOOVVV.f90 similarity index 100% rename from src/QuAcK/form_delta_OOOVVV.f90 rename to src/CC/form_delta_OOOVVV.f90 diff --git a/src/QuAcK/form_delta_OOVV.f90 b/src/CC/form_delta_OOVV.f90 similarity index 100% rename from src/QuAcK/form_delta_OOVV.f90 rename to src/CC/form_delta_OOVV.f90 diff --git a/src/QuAcK/form_delta_OV.f90 b/src/CC/form_delta_OV.f90 similarity index 100% rename from src/QuAcK/form_delta_OV.f90 rename to src/CC/form_delta_OV.f90 diff --git a/src/QuAcK/form_g.f90 b/src/CC/form_g.f90 similarity index 100% rename from src/QuAcK/form_g.f90 rename to src/CC/form_g.f90 diff --git a/src/QuAcK/form_h.f90 b/src/CC/form_h.f90 similarity index 100% rename from src/QuAcK/form_h.f90 rename to src/CC/form_h.f90 diff --git a/src/QuAcK/form_ladder_r.f90 b/src/CC/form_ladder_r.f90 similarity index 100% rename from src/QuAcK/form_ladder_r.f90 rename to src/CC/form_ladder_r.f90 diff --git a/src/QuAcK/form_r1.f90 b/src/CC/form_r1.f90 similarity index 100% rename from src/QuAcK/form_r1.f90 rename to src/CC/form_r1.f90 diff --git a/src/QuAcK/form_r1_nc.f90 b/src/CC/form_r1_nc.f90 similarity index 100% rename from src/QuAcK/form_r1_nc.f90 rename to src/CC/form_r1_nc.f90 diff --git a/src/QuAcK/form_r2.f90 b/src/CC/form_r2.f90 similarity index 100% rename from src/QuAcK/form_r2.f90 rename to src/CC/form_r2.f90 diff --git a/src/QuAcK/form_r2_nc.f90 b/src/CC/form_r2_nc.f90 similarity index 100% rename from src/QuAcK/form_r2_nc.f90 rename to src/CC/form_r2_nc.f90 diff --git a/src/QuAcK/form_ring_r.f90 b/src/CC/form_ring_r.f90 similarity index 100% rename from src/QuAcK/form_ring_r.f90 rename to src/CC/form_ring_r.f90 diff --git a/src/QuAcK/form_tau.f90 b/src/CC/form_tau.f90 similarity index 100% rename from src/QuAcK/form_tau.f90 rename to src/CC/form_tau.f90 diff --git a/src/QuAcK/form_tau_nc.f90 b/src/CC/form_tau_nc.f90 similarity index 100% rename from src/QuAcK/form_tau_nc.f90 rename to src/CC/form_tau_nc.f90 diff --git a/src/QuAcK/form_taus_nc.f90 b/src/CC/form_taus_nc.f90 similarity index 100% rename from src/QuAcK/form_taus_nc.f90 rename to src/CC/form_taus_nc.f90 diff --git a/src/QuAcK/form_u.f90 b/src/CC/form_u.f90 similarity index 100% rename from src/QuAcK/form_u.f90 rename to src/CC/form_u.f90 diff --git a/src/QuAcK/form_ub.f90 b/src/CC/form_ub.f90 similarity index 100% rename from src/QuAcK/form_ub.f90 rename to src/CC/form_ub.f90 diff --git a/src/QuAcK/form_ubb.f90 b/src/CC/form_ubb.f90 similarity index 100% rename from src/QuAcK/form_ubb.f90 rename to src/CC/form_ubb.f90 diff --git a/src/QuAcK/form_v.f90 b/src/CC/form_v.f90 similarity index 100% rename from src/QuAcK/form_v.f90 rename to src/CC/form_v.f90 diff --git a/src/QuAcK/lCCD.f90 b/src/CC/lCCD.f90 similarity index 100% rename from src/QuAcK/lCCD.f90 rename to src/CC/lCCD.f90 diff --git a/src/QuAcK/pCCD.f90 b/src/CC/pCCD.f90 similarity index 100% rename from src/QuAcK/pCCD.f90 rename to src/CC/pCCD.f90 diff --git a/src/QuAcK/rCCD.f90 b/src/CC/rCCD.f90 similarity index 100% rename from src/QuAcK/rCCD.f90 rename to src/CC/rCCD.f90 diff --git a/src/QuAcK/CIS.f90 b/src/CI/CIS.f90 similarity index 100% rename from src/QuAcK/CIS.f90 rename to src/CI/CIS.f90 diff --git a/src/QuAcK/CISD.f90 b/src/CI/CISD.f90 similarity index 100% rename from src/QuAcK/CISD.f90 rename to src/CI/CISD.f90 diff --git a/src/CI/Makefile b/src/CI/Makefile new file mode 100644 index 0000000..22a571c --- /dev/null +++ b/src/CI/Makefile @@ -0,0 +1,2 @@ +TARGET=CI.a +include ../Makefile.include diff --git a/src/QuAcK/UCIS.f90 b/src/CI/UCIS.f90 similarity index 100% rename from src/QuAcK/UCIS.f90 rename to src/CI/UCIS.f90 diff --git a/src/QuAcK/MOM.f90 b/src/HF/MOM.f90 similarity index 100% rename from src/QuAcK/MOM.f90 rename to src/HF/MOM.f90 diff --git a/src/QuAcK/MOM_overlap.f90 b/src/HF/MOM_overlap.f90 similarity index 100% rename from src/QuAcK/MOM_overlap.f90 rename to src/HF/MOM_overlap.f90 diff --git a/src/HF/Makefile b/src/HF/Makefile new file mode 100644 index 0000000..e33698c --- /dev/null +++ b/src/HF/Makefile @@ -0,0 +1,2 @@ +TARGET=HF.a +include ../Makefile.include diff --git a/src/QuAcK/RHF.f90 b/src/HF/RHF.f90 similarity index 100% rename from src/QuAcK/RHF.f90 rename to src/HF/RHF.f90 diff --git a/src/QuAcK/UHF.f90 b/src/HF/UHF.f90 similarity index 100% rename from src/QuAcK/UHF.f90 rename to src/HF/UHF.f90 diff --git a/src/QuAcK/density.f90 b/src/HF/density.f90 similarity index 100% rename from src/QuAcK/density.f90 rename to src/HF/density.f90 diff --git a/src/QuAcK/density_matrix.f90 b/src/HF/density_matrix.f90 similarity index 100% rename from src/QuAcK/density_matrix.f90 rename to src/HF/density_matrix.f90 diff --git a/src/QuAcK/huckel_guess.f90 b/src/HF/huckel_guess.f90 similarity index 100% rename from src/QuAcK/huckel_guess.f90 rename to src/HF/huckel_guess.f90 diff --git a/src/QuAcK/mo_guess.f90 b/src/HF/mo_guess.f90 similarity index 100% rename from src/QuAcK/mo_guess.f90 rename to src/HF/mo_guess.f90 diff --git a/src/HF/obj/.gitignore b/src/HF/obj/.gitignore new file mode 100644 index 0000000..e69de29 diff --git a/src/QuAcK/print_RHF.f90 b/src/HF/print_RHF.f90 similarity index 100% rename from src/QuAcK/print_RHF.f90 rename to src/HF/print_RHF.f90 diff --git a/src/QuAcK/print_UHF.f90 b/src/HF/print_UHF.f90 similarity index 100% rename from src/QuAcK/print_UHF.f90 rename to src/HF/print_UHF.f90 diff --git a/src/IntPak/CalcBoysF.f90 b/src/IntPak/CalcBoysF.f90 deleted file mode 100644 index 68fa43a..0000000 --- a/src/IntPak/CalcBoysF.f90 +++ /dev/null @@ -1,47 +0,0 @@ -!module c_functions -! use iso_c_binding -! interface -! function gsl_sf_gamma_inc_P(a,t) bind(C, name="gsl_sf_gamma_inc_P") -! use iso_c_binding, only: c_double -! real(kind=c_double), value :: a,t -! real(kind=c_double) :: gsl_sf_gamma_inc_P -! end function gsl_sf_gamma_inc_P -! end interface -!end module - -subroutine CalcBoysF(maxm,t,Fm) -! use c_functions -! Comute the generalized Boys function Fm(t) using Slatec library - - implicit none - -! Input variables - - double precision,intent(in) :: t - integer,intent(in) :: maxm - -! Local variables - - integer :: m - double precision :: dm - double precision :: dgami - - -! Output variables - - double precision,intent(inout):: Fm(0:maxm) - - if(t == 0d0) then - do m=0,maxm - dm = dble(m) - Fm(m) = 1d0/(2d0*dm+1d0) - end do - else - do m=0,maxm - dm = dble(m) -! Fm(m) = gamma(dm+0.5d0)*gsl_sf_gamma_inc_P(dm+0.5d0,t)/(2d0*t**(dm+0.5d0)) - Fm(m) = dgami(dm+0.5d0,t)/(2d0*t**(dm+0.5d0)) - end do - end if - -end subroutine CalcBoysF diff --git a/src/IntPak/CalcOm.f90 b/src/IntPak/CalcOm.f90 deleted file mode 100644 index 8b5c8e7..0000000 --- a/src/IntPak/CalcOm.f90 +++ /dev/null @@ -1,40 +0,0 @@ -subroutine CalcOm(maxm,ExpPQi,NormPQSq,Om) - -! Comute the 0^m: (00|00)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: ExpPQi,NormPQSq - -! Local variables - - integer :: m - double precision :: pi,dm,t - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om (0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Campute generalized Boys functions - - t = NormPQSq/ExpPQi - call CalcBoysF(maxm,t,Fm) - -! Compute (00|00)^m - - do m=0,maxm - dm =dble(m) - Om(m) = (2d0/sqrt(pi))*(-1d0)**dm*(1d0/ExpPQi)**(dm+0.5d0)*Fm(m) - end do - - deallocate(Fm) - -end subroutine CalcOm diff --git a/src/IntPak/CalcOm3e.f90 b/src/IntPak/CalcOm3e.f90 deleted file mode 100644 index 76a3525..0000000 --- a/src/IntPak/CalcOm3e.f90 +++ /dev/null @@ -1,44 +0,0 @@ -subroutine CalcOm3e(maxm,delta0,delta1,Y1,Y0,Om) - -! Compute the 0^m for ERIs: (00|00)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: delta0,delta1,Y0,Y1 - -! Local variables - - integer :: m - double precision :: pi,t,OG - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om (0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Calculate OG - - OG = (pi**4/delta0)**(3d0/2d0)*exp(-Y0) - -! Campute generalized Boys functions - - t = delta1/(delta1-delta0)*(Y1-Y0) - call CalcBoysF(maxm,t,Fm) - -! Compute (000|000)^m - - do m=0,maxm - Om(m) = (2d0/sqrt(pi))*OG*sqrt(delta0/(delta1-delta0))*(delta1/(delta1-delta0))**m - Om(m) = Om(m)*Fm(m) - end do - - deallocate(Fm) - -end subroutine CalcOm3e diff --git a/src/IntPak/CalcOmERI.f90 b/src/IntPak/CalcOmERI.f90 deleted file mode 100644 index 05c4abc..0000000 --- a/src/IntPak/CalcOmERI.f90 +++ /dev/null @@ -1,39 +0,0 @@ -subroutine CalcOmERI(maxm,ExpY,NormYSq,Om) - -! Compute the 0^m for ERIs: (00|00)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: ExpY,NormYSq - -! Local variables - - integer :: m - double precision :: pi,t - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om (0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Campute generalized Boys functions - - t = ExpY*NormYSq - call CalcBoysF(maxm,t,Fm) - -! Compute (00|00)^m - - do m=0,maxm - Om(m) = (2d0/sqrt(pi))*sqrt(ExpY)*Fm(m) - end do - - deallocate(Fm) - -end subroutine CalcOmERI diff --git a/src/IntPak/CalcOmErf.f90 b/src/IntPak/CalcOmErf.f90 deleted file mode 100644 index ef1535a..0000000 --- a/src/IntPak/CalcOmErf.f90 +++ /dev/null @@ -1,39 +0,0 @@ -subroutine CalcOmErf(maxm,ExpY,fG,NormYSq,Om) - -! Compute the 0^m for the long-range Coulomb operator: (00|erf(mu r)/r|00)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: ExpY,fG,NormYSq - -! Local variables - - integer :: m - double precision :: pi,t - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om (0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Campute generalized Boys functions - - t = fG*NormYSq - call CalcBoysF(maxm,t,Fm) - -! Compute (00|00)^m - - do m=0,maxm - Om(m) = (2d0/sqrt(pi))*sqrt(fG)*(fG/ExpY)**m*Fm(m) - end do - - deallocate(Fm) - -end subroutine CalcOmErf diff --git a/src/IntPak/CalcOmNuc.f90 b/src/IntPak/CalcOmNuc.f90 deleted file mode 100644 index 5f5e730..0000000 --- a/src/IntPak/CalcOmNuc.f90 +++ /dev/null @@ -1,40 +0,0 @@ -subroutine CalcOmNuc(maxm,ExpPQi,NormPQSq,Om) - -! Compute (0|V|0)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: ExpPQi,NormPQSq - -! Local variables - - integer :: m - double precision :: pi,dm,t - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om (0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Campute generalized Boys functions - - t = NormPQSq/ExpPQi - call CalcBoysF(maxm,t,Fm) - -! Compute (00|00)^m - - do m=0,maxm - dm =dble(m) - Om(m) = (2d0/sqrt(pi))*(1d0/ExpPQi)**(dm+0.5d0)*Fm(m) - end do - - deallocate(Fm) - -end subroutine CalcOmNuc diff --git a/src/IntPak/CalcOmYuk.f90 b/src/IntPak/CalcOmYuk.f90 deleted file mode 100644 index 739c10e..0000000 --- a/src/IntPak/CalcOmYuk.f90 +++ /dev/null @@ -1,43 +0,0 @@ -subroutine CalcOmYuk(maxm,ExpG,ExpY,fG,NormYSq,Om) - -! Compute the 0^m for the screened Coulomb operator: (00|f12/r12|00)^m - - implicit none - -! Input variables - - integer,intent(in) :: maxm - double precision,intent(in) :: ExpG,ExpY,fG,NormYSq - -! Local variables - - integer :: m,k - double precision :: pi,t,dbinom - double precision,allocatable :: Fm(:) - -! Output variables - - double precision,intent(inout):: Om(0:maxm) - - allocate(Fm(0:maxm)) - - pi = 4d0*atan(1d0) - -! Campute generalized Boys functions - - t = (ExpY - fG)*NormYSq - call CalcBoysF(maxm,t,Fm) - -! Compute (00|00)^m - - do m=0,maxm - Om(m) = 0d0 - do k=0,m - Om(m) = Om(m) + dbinom(m,k)*(ExpY/ExpG)**k*Fm(k) - end do - Om(m) = (2d0/sqrt(pi))*sqrt(ExpY)*(fG/ExpG)*exp(-fG*NormYSq)*Om(m) - end do - - deallocate(Fm) - -end subroutine CalcOmYuk diff --git a/src/IntPak/Compute2eInt.f90 b/src/IntPak/Compute2eInt.f90 deleted file mode 100644 index 35f684c..0000000 --- a/src/IntPak/Compute2eInt.f90 +++ /dev/null @@ -1,308 +0,0 @@ -subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np2eInt,nSigp2eInt,nc2eInt,nSigc2eInt) - - -! Compute various two-electron integrals - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: debug - logical,intent(in) :: chemist_notation - integer,intent(in) :: iType,nShell - double precision,intent(in) :: ExpS - integer,intent(in) :: KG - double precision,intent(in) :: DG(KG),ExpG(KG) - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -! Local variables - - integer :: KBra(2),KKet(2) - double precision :: CenterBra(2,3),CenterKet(2,3) - integer :: TotAngMomBra(2),TotAngMomKet(2) - integer :: AngMomBra(2,3),AngMomKet(2,3) - integer :: nShellFunctionBra(2),nShellFunctionKet(2) - integer,allocatable :: ShellFunctionA1(:,:),ShellFunctionA2(:,:) - integer,allocatable :: ShellFunctionB1(:,:),ShellFunctionB2(:,:) - double precision :: ExpBra(2),ExpKet(2) - double precision :: DBra(2),DKet(2) - double precision :: norm_coeff - - integer :: iBasA1,iBasA2,iBasB1,iBasB2 - integer :: iShA1,iShA2,iShB1,iShB2 - integer :: iShFA1,iShFA2,iShFB1,iShFB2 - integer :: iKA1,iKA2,iKB1,iKB2 - integer :: iFile - - double precision :: p2eInt,c2eInt - double precision :: start_c2eInt,end_c2eInt,t_c2eInt - -! Output variables - - integer,intent(out) :: np2eInt,nSigp2eInt,nc2eInt,nSigc2eInt - - np2eInt = 0 - nSigp2eInt = 0 - - nc2eInt = 0 - nSigc2eInt = 0 - - iBasA1 = 0 - iBasA2 = 0 - iBasB1 = 0 - iBasB2 = 0 - -! Open file to write down integrals - - iFile = 0 - - if(iType == 1) then - -! Compute two-electron integrals over the Coulomb operator - - write(*,*) '******************************************' - write(*,*) ' Compute two-electron repulsion integrals ' - write(*,*) '******************************************' - write(*,*) - - iFile = 21 - open(unit=iFile,file='int/ERI.dat') - - elseif(iType == 2) then - -! Compute two-electron integrals over Slater geminals - - write(*,*) '****************************************' - write(*,*) ' Compute two-electron geminal integrals ' - write(*,*) '****************************************' - write(*,*) - - iFile = 22 - open(unit=iFile,file='int/F12.dat') - - elseif(iType == 3) then - -! Compute two-electron integrals over the Yukawa operator - - write(*,*) '***************************************' - write(*,*) ' Compute two-electron Yukawa integrals ' - write(*,*) '***************************************' - write(*,*) - - iFile = 23 - open(unit=iFile,file='int/Yuk.dat') - - elseif(iType == 4) then - -! Compute two-electron integrals over the long-range Coulomb operator - - write(*,*) '**************************************' - write(*,*) ' Compute long-range Coulomb integrals ' - write(*,*) '**************************************' - write(*,*) - - iFile = 24 - open(unit=iFile,file='int/Erf.dat') - - end if - -!------------------------------------------------------------------------ -! Loops over shell A1 -!------------------------------------------------------------------------ - do iShA1=1,nShell - - CenterBra(1,1) = CenterShell(iShA1,1) - CenterBra(1,2) = CenterShell(iShA1,2) - CenterBra(1,3) = CenterShell(iShA1,3) - - TotAngMomBra(1) = TotAngMomShell(iShA1) - nShellFunctionBra(1) = (TotAngMomBra(1)*TotAngMomBra(1) + 3*TotAngMomBra(1) + 2)/2 - allocate(ShellFunctionA1(1:nShellFunctionBra(1),1:3)) - call GenerateShell(TotAngMomBra(1),nShellFunctionBra(1),ShellFunctionA1) - - KBra(1) = KShell(iShA1) - - do iShFA1=1,nShellFunctionBra(1) - - iBasA1 = iBasA1 + 1 - AngMomBra(1,1) = ShellFunctionA1(iShFA1,1) - AngMomBra(1,2) = ShellFunctionA1(iShFA1,2) - AngMomBra(1,3) = ShellFunctionA1(iShFA1,3) - -!------------------------------------------------------------------------ -! Loops over shell B1 -!------------------------------------------------------------------------ - do iShB1=1,iShA1 - - CenterKet(1,1) = CenterShell(iShB1,1) - CenterKet(1,2) = CenterShell(iShB1,2) - CenterKet(1,3) = CenterShell(iShB1,3) - - TotAngMomKet(1) = TotAngMomShell(iShB1) - nShellFunctionKet(1) = (TotAngMomKet(1)*TotAngMomKet(1) + 3*TotAngMomKet(1) + 2)/2 - allocate(ShellFunctionB1(1:nShellFunctionKet(1),1:3)) - call GenerateShell(TotAngMomKet(1),nShellFunctionKet(1),ShellFunctionB1) - - KKet(1) = KShell(iShB1) - - do iShFB1=1,nShellFunctionKet(1) - - iBasB1 = iBasB1 + 1 - AngMomKet(1,1) = ShellFunctionB1(iShFB1,1) - AngMomKet(1,2) = ShellFunctionB1(iShFB1,2) - AngMomKet(1,3) = ShellFunctionB1(iShFB1,3) - -!------------------------------------------------------------------------ -! Loops over shell A2 -!------------------------------------------------------------------------ - do iShA2=1,iShA1 - - CenterBra(2,1) = CenterShell(iShA2,1) - CenterBra(2,2) = CenterShell(iShA2,2) - CenterBra(2,3) = CenterShell(iShA2,3) - - TotAngMomBra(2) = TotAngMomShell(iShA2) - nShellFunctionBra(2) = (TotAngMomBra(2)*TotAngMomBra(2) + 3*TotAngMomBra(2) + 2)/2 - allocate(ShellFunctionA2(1:nShellFunctionBra(2),1:3)) - call GenerateShell(TotAngMomBra(2),nShellFunctionBra(2),ShellFunctionA2) - - KBra(2) = KShell(iShA2) - - do iShFA2=1,nShellFunctionBra(2) - - iBasA2 = iBasA2 + 1 - AngMomBra(2,1) = ShellFunctionA2(iShFA2,1) - AngMomBra(2,2) = ShellFunctionA2(iShFA2,2) - AngMomBra(2,3) = ShellFunctionA2(iShFA2,3) - -!------------------------------------------------------------------------ -! Loops over shell B2 -!------------------------------------------------------------------------ - do iShB2=1,iShA2 - - CenterKet(2,1) = CenterShell(iShB2,1) - CenterKet(2,2) = CenterShell(iShB2,2) - CenterKet(2,3) = CenterShell(iShB2,3) - - TotAngMomKet(2) = TotAngMomShell(iShB2) - nShellFunctionKet(2) = (TotAngMomKet(2)*TotAngMomKet(2) + 3*TotAngMomKet(2) + 2)/2 - allocate(ShellFunctionB2(1:nShellFunctionKet(2),1:3)) - call GenerateShell(TotAngMomKet(2),nShellFunctionKet(2),ShellFunctionB2) - - KKet(2) = KShell(iShB2) - - do iShFB2=1,nShellFunctionKet(2) - - iBasB2 = iBasB2 + 1 - AngMomKet(2,1) = ShellFunctionB2(iShFB2,1) - AngMomKet(2,2) = ShellFunctionB2(iShFB2,2) - AngMomKet(2,3) = ShellFunctionB2(iShFB2,3) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - call cpu_time(start_c2eInt) - - c2eInt = 0d0 - - do iKA1=1,KBra(1) - ExpBra(1) = ExpShell(iShA1,iKA1) - DBra(1) = DShell(iShA1,iKA1)*norm_coeff(ExpBra(1),AngMomBra(1,1:3)) - do iKA2=1,KBra(2) - ExpBra(2) = ExpShell(iShA2,iKA2) - DBra(2) = DShell(iShA2,iKA2)*norm_coeff(ExpBra(2),AngMomBra(2,1:3)) - do iKB1=1,KKet(1) - ExpKet(1) = ExpShell(iShB1,iKB1) - DKet(1) = DShell(iShB1,iKB1)*norm_coeff(ExpKet(1),AngMomKet(1,1:3)) - do iKB2=1,KKet(2) - ExpKet(2) = ExpShell(iShB2,iKB2) - DKet(2) = DShell(iShB2,iKB2)*norm_coeff(ExpKet(2),AngMomKet(2,1:3)) - - call S2eInt(debug,iType,np2eInt,nSigp2eInt, & - ExpS,KG,DG,ExpG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet, & - p2eInt) - - c2eInt = c2eInt + DBra(1)*DBra(2)*DKet(1)*DKet(2)*p2eInt - - end do - end do - end do - end do - call cpu_time(end_c2eInt) - - nc2eInt = nc2eInt + 1 - - if(abs(c2eInt) > 1d-15) then - - nSigc2eInt = nSigc2eInt + 1 - t_c2eInt = end_c2eInt - start_c2eInt - - if(chemist_notation) then - - write(iFile,'(I6,I6,I6,I6,F20.15)') iBasA1,iBasB1,iBasA2,iBasB2,c2eInt -! write(iFile,'(F20.15,I6,I6,I6,I6)') c2eInt,iBasA1,iBasB1,iBasA2,iBasB2 - - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & - '(a1b1|a2b2) = ',c2eInt,iBasA1,iBasB1,iBasA2,iBasB2 - end if - - else - - write(iFile,'(I6,I6,I6,I6,F20.15)') iBasA1,iBasA2,iBasB1,iBasB2,c2eInt -! write(iFile,'(F20.15,I6,I6,I6,I6)') c2eInt,iBasA1,iBasA2,iBasB1,iBasB2 - - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & - ' = ',c2eInt,iBasA1,iBasA2,iBasB1,iBasB2 - end if - - end if - end if - -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB2) - end do - iBasB2 = 0 -!------------------------------------------------------------------------ -! End loops over shell B2 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA2) - end do - iBasA2 = 0 -!------------------------------------------------------------------------ -! End loops over shell A2 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB1) - end do - iBasB1 = 0 -!------------------------------------------------------------------------ -! End loops over shell B1 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA1) - end do - iBasA1 = 0 -!------------------------------------------------------------------------ -! End loops over shell A1 -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=iFile) - -end subroutine Compute2eInt diff --git a/src/IntPak/Compute3eInt.f90 b/src/IntPak/Compute3eInt.f90 deleted file mode 100644 index accf6de..0000000 --- a/src/IntPak/Compute3eInt.f90 +++ /dev/null @@ -1,328 +0,0 @@ -subroutine Compute3eInt(debug,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np3eInt,nSigp3eInt,nc3eInt,nSigc3eInt) - - -! Compute long-range Coulomb integrals - - implicit none - include 'parameters.h' - - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: iType,nShell - double precision,intent(in) :: ExpS - integer,intent(in) :: KG - double precision,intent(in) :: DG(KG),ExpG(KG) - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -! Local variables - - integer :: KBra(3),KKet(3) - double precision :: CenterBra(3,3),CenterKet(3,3) - integer :: TotAngMomBra(3),TotAngMomKet(3) - integer :: AngMomBra(3,3),AngMomKet(3,3) - integer :: nShellFunctionBra(3),nShellFunctionKet(3) - integer,allocatable :: ShellFunctionA1(:,:),ShellFunctionA2(:,:),ShellFunctionA3(:,:) - integer,allocatable :: ShellFunctionB1(:,:),ShellFunctionB2(:,:),ShellFunctionB3(:,:) - double precision :: ExpBra(3),ExpKet(3) - double precision :: DBra(3),DKet(3) - double precision :: norm_coeff - - integer :: iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3 - integer :: iShA1,iShA2,iShA3,iShB1,iShB2,iShB3 - integer :: iShFA1,iShFA2,iShFA3,iShFB1,iShFB2,iShFB3 - integer :: iKA1,iKA2,iKA3,iKB1,iKB2,iKB3 - integer :: iFile - - double precision :: p3eInt,c3eInt - double precision :: start_c3eInt,end_c3eInt,t_c3eInt - -! Output variables - - integer,intent(out) :: np3eInt,nSigp3eInt,nc3eInt,nSigc3eInt - -! Compute three-electron integrals - - write(*,*) '**********************************' - write(*,*) ' Compute three-electron integrals ' - write(*,*) '**********************************' - write(*,*) - - np3eInt = 0 - nSigp3eInt = 0 - - nc3eInt = 0 - nSigc3eInt = 0 - - iBasA1 = 0 - iBasA2 = 0 - iBasA3 = 0 - iBasB1 = 0 - iBasB2 = 0 - iBasB3 = 0 - -! Open file to write down integrals - - iFile = 0 - - if(iType == 1) then - iFile = 31 - open(unit=iFile,file='int/3eInt_Type1.dat') - elseif(iType == 2) then - iFile = 32 - open(unit=iFile,file='int/3eInt_Type2.dat') - elseif(iType == 3) then - iFile = 33 - open(unit=iFile,file='int/3eInt_Type3.dat') - end if - -!------------------------------------------------------------------------ -! Loops over shell A1 -!------------------------------------------------------------------------ - do iShA1=1,nShell - - CenterBra(1,1) = CenterShell(iShA1,1) - CenterBra(1,2) = CenterShell(iShA1,2) - CenterBra(1,3) = CenterShell(iShA1,3) - - TotAngMomBra(1) = TotAngMomShell(iShA1) - nShellFunctionBra(1) = (TotAngMomBra(1)*TotAngMomBra(1) + 3*TotAngMomBra(1) + 2)/2 - allocate(ShellFunctionA1(1:nShellFunctionBra(1),1:3)) - call GenerateShell(TotAngMomBra(1),nShellFunctionBra(1),ShellFunctionA1) - - KBra(1) = KShell(iShA1) - - do iShFA1=1,nShellFunctionBra(1) - - iBasA1 = iBasA1 + 1 - AngMomBra(1,1) = ShellFunctionA1(iShFA1,1) - AngMomBra(1,2) = ShellFunctionA1(iShFA1,2) - AngMomBra(1,3) = ShellFunctionA1(iShFA1,3) - -!------------------------------------------------------------------------ -! Loops over shell A2 -!------------------------------------------------------------------------ - do iShA2=1,nShell - - CenterBra(2,1) = CenterShell(iShA2,1) - CenterBra(2,2) = CenterShell(iShA2,2) - CenterBra(2,3) = CenterShell(iShA2,3) - - TotAngMomBra(2) = TotAngMomShell(iShA2) - nShellFunctionBra(2) = (TotAngMomBra(2)*TotAngMomBra(2) + 3*TotAngMomBra(2) + 2)/2 - allocate(ShellFunctionA2(1:nShellFunctionBra(2),1:3)) - call GenerateShell(TotAngMomBra(2),nShellFunctionBra(2),ShellFunctionA2) - - KBra(2) = KShell(iShA2) - - do iShFA2=1,nShellFunctionBra(2) - - iBasA2 = iBasA2 + 1 - AngMomBra(2,1) = ShellFunctionA2(iShFA2,1) - AngMomBra(2,2) = ShellFunctionA2(iShFA2,2) - AngMomBra(2,3) = ShellFunctionA2(iShFA2,3) - -!------------------------------------------------------------------------ -! Loops over shell A3 -!------------------------------------------------------------------------ - do iShA3=1,nShell - - CenterBra(3,1) = CenterShell(iShA3,1) - CenterBra(3,2) = CenterShell(iShA3,2) - CenterBra(3,3) = CenterShell(iShA3,3) - - TotAngMomBra(3) = TotAngMomShell(iShA3) - nShellFunctionBra(3) = (TotAngMomBra(3)*TotAngMomBra(3) + 3*TotAngMomBra(3) + 2)/2 - allocate(ShellFunctionA3(1:nShellFunctionBra(3),1:3)) - call GenerateShell(TotAngMomBra(3),nShellFunctionBra(3),ShellFunctionA3) - - KBra(3) = KShell(iShA3) - - do iShFA3=1,nShellFunctionBra(3) - - iBasA3 = iBasA3 + 1 - AngMomBra(3,1) = ShellFunctionA3(iShFA3,1) - AngMomBra(3,2) = ShellFunctionA3(iShFA3,2) - AngMomBra(3,3) = ShellFunctionA3(iShFA3,3) - -!------------------------------------------------------------------------ -! Loops over shell B1 -!------------------------------------------------------------------------ - do iShB1=1,nShell - - CenterKet(1,1) = CenterShell(iShB1,1) - CenterKet(1,2) = CenterShell(iShB1,2) - CenterKet(1,3) = CenterShell(iShB1,3) - - TotAngMomKet(1) = TotAngMomShell(iShB1) - nShellFunctionKet(1) = (TotAngMomKet(1)*TotAngMomKet(1) + 3*TotAngMomKet(1) + 2)/2 - allocate(ShellFunctionB1(1:nShellFunctionKet(1),1:3)) - call GenerateShell(TotAngMomKet(1),nShellFunctionKet(1),ShellFunctionB1) - - KKet(1) = KShell(iShB1) - - do iShFB1=1,nShellFunctionKet(1) - - iBasB1 = iBasB1 + 1 - AngMomKet(1,1) = ShellFunctionB1(iShFB1,1) - AngMomKet(1,2) = ShellFunctionB1(iShFB1,2) - AngMomKet(1,3) = ShellFunctionB1(iShFB1,3) - -!------------------------------------------------------------------------ -! Loops over shell B2 -!------------------------------------------------------------------------ - do iShB2=1,nShell - - CenterKet(2,1) = CenterShell(iShB2,1) - CenterKet(2,2) = CenterShell(iShB2,2) - CenterKet(2,3) = CenterShell(iShB2,3) - - TotAngMomKet(2) = TotAngMomShell(iShB2) - nShellFunctionKet(2) = (TotAngMomKet(2)*TotAngMomKet(2) + 3*TotAngMomKet(2) + 2)/2 - allocate(ShellFunctionB2(1:nShellFunctionKet(2),1:3)) - call GenerateShell(TotAngMomKet(2),nShellFunctionKet(2),ShellFunctionB2) - - KKet(2) = KShell(iShB2) - - do iShFB2=1,nShellFunctionKet(2) - - iBasB2 = iBasB2 + 1 - AngMomKet(2,1) = ShellFunctionB2(iShFB2,1) - AngMomKet(2,2) = ShellFunctionB2(iShFB2,2) - AngMomKet(2,3) = ShellFunctionB2(iShFB2,3) - -!------------------------------------------------------------------------ -! Loops over shell B3 -!------------------------------------------------------------------------ - do iShB3=1,nShell - - CenterKet(3,1) = CenterShell(iShB3,1) - CenterKet(3,2) = CenterShell(iShB3,2) - CenterKet(3,3) = CenterShell(iShB3,3) - - TotAngMomKet(3) = TotAngMomShell(iShB3) - nShellFunctionKet(3) = (TotAngMomKet(3)*TotAngMomKet(3) + 3*TotAngMomKet(3) + 2)/2 - allocate(ShellFunctionB3(1:nShellFunctionKet(3),1:3)) - call GenerateShell(TotAngMomKet(3),nShellFunctionKet(3),ShellFunctionB3) - - KKet(3) = KShell(iShB3) - - do iShFB3=1,nShellFunctionKet(3) - - iBasB3 = iBasB3 + 1 - AngMomKet(3,1) = ShellFunctionB3(iShFB3,1) - AngMomKet(3,2) = ShellFunctionB3(iShFB3,2) - AngMomKet(3,3) = ShellFunctionB3(iShFB3,3) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - call cpu_time(start_c3eInt) - - c3eInt = 0d0 - - do iKA1=1,KBra(1) - ExpBra(1) = ExpShell(iShA1,iKA1) - DBra(1) = DShell(iShA1,iKA1)*norm_coeff(ExpBra(1),AngMomBra(1,1:3)) - do iKA2=1,KBra(2) - ExpBra(2) = ExpShell(iShA2,iKA2) - DBra(2) = DShell(iShA2,iKA2)*norm_coeff(ExpBra(2),AngMomBra(2,1:3)) - do iKA3=1,KBra(3) - ExpBra(3) = ExpShell(iShA3,iKA3) - DBra(3) = DShell(iShA3,iKA3)*norm_coeff(ExpBra(3),AngMomBra(3,1:3)) - do iKB1=1,KKet(1) - ExpKet(1) = ExpShell(iShB1,iKB1) - DKet(1) = DShell(iShB1,iKB1)*norm_coeff(ExpKet(1),AngMomKet(1,1:3)) - do iKB2=1,KKet(2) - ExpKet(2) = ExpShell(iShB2,iKB2) - DKet(2) = DShell(iShB2,iKB2)*norm_coeff(ExpKet(2),AngMomKet(2,1:3)) - do iKB3=1,KKet(3) - ExpKet(3) = ExpShell(iShB3,iKB3) - DKet(3) = DShell(iShB3,iKB3)*norm_coeff(ExpKet(3),AngMomKet(3,1:3)) - - call S3eInt(debug,iType,np3eInt,nSigp3eInt, & - ExpS,KG,DG,ExpG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet, & - p3eInt) - - c3eInt = c3eInt + DBra(1)*DBra(2)*DBra(3)*DKet(1)*DKet(2)*DKet(3)*p3eInt - - end do - end do - end do - end do - end do - end do - call cpu_time(end_c3eInt) - - nc3eInt = nc3eInt + 1 - if(abs(c3eInt) > 1d-15) then - nSigc3eInt = nSigc3eInt + 1 - t_c3eInt = end_c3eInt - start_c3eInt - write(iFile,'(I9,I9,I9,I9,I9,I9,F25.15)') & - iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3,c3eInt - if(.true.) then - write(*,'(A15,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6,1X,F16.10)') & - '(a1a2a3|b1b2b3) = ',iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3,c3eInt - end if - end if - -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB3) - end do - iBasB3 = 0 -!------------------------------------------------------------------------ -! End loops over shell B3 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB2) - end do - iBasB2 = 0 -!------------------------------------------------------------------------ -! End loops over shell B2 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB1) - end do - iBasB1 = 0 -!------------------------------------------------------------------------ -! End loops over shell B1 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA3) - end do - iBasA3 = 0 -!------------------------------------------------------------------------ -! End loops over shell A3 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA2) - end do - iBasA2 = 0 -!------------------------------------------------------------------------ -! End loops over shell A2 -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA1) - end do - iBasA1 = 0 -!------------------------------------------------------------------------ -! End loops over shell A1 -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=iFile) - -end subroutine Compute3eInt diff --git a/src/IntPak/Compute4eInt.f90 b/src/IntPak/Compute4eInt.f90 deleted file mode 100644 index 3127942..0000000 --- a/src/IntPak/Compute4eInt.f90 +++ /dev/null @@ -1,246 +0,0 @@ -subroutine Compute4eInt(debug,nEl,iType,nShell,ExpS, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - npErf,nSigpErf,ncErf,nSigcErf) - - -! Compute long-range Coulomb integrals - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: nEl,iType,nShell - double precision :: ExpS - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -! Local variables - - integer :: KA,KB,KC,KD - double precision :: CenterA(3),CenterB(3),CenterC(3),CenterD(3) - integer :: TotAngMomA,TotAngMomB,TotAngMomC,TotAngMomD - integer :: AngMomA(3),AngMomB(3),AngMomC(3),AngMomD(3) - integer :: nShellFunctionA,nShellFunctionB, & - nShellFunctionC,nShellFunctionD - integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:), & - ShellFunctionC(:,:),ShellFunctionD(:,:) - double precision :: ExpA,ExpB,ExpC,ExpD - double precision,allocatable :: DA,DB,DC,DD - double precision :: norm_coeff - - integer :: iBasA,iBasB,iBasC,iBasD - integer :: iShA,iShB,iShC,iShD - integer :: iShFA,iShFB,iShFC,iShFD - integer :: iKA,iKB,iKC,iKD - - double precision :: pErf,cErf - double precision :: start_cErf,end_cErf,t_cErf - -! Output variables - - integer,intent(out) :: npErf,nSigpErf,ncErf,nSigcErf - -! Compute two-electron integrals over long-range Coulomb operator - - write(*,*) '**********************************' - write(*,*) ' Compute three-electron integrals ' - write(*,*) '**********************************' - write(*,*) - - npErf = 0 - nSigpErf = 0 - - ncErf = 0 - nSigcErf = 0 - - iBasA = 0 - iBasB = 0 - iBasC = 0 - iBasD = 0 - -! Open file to write down integrals - - open(unit=41,file='int/4eInt_Type1.dat') - -!------------------------------------------------------------------------ -! Loops over shell A -!------------------------------------------------------------------------ - do iShA=1,nShell - - CenterA(1) = CenterShell(iShA,1) - CenterA(2) = CenterShell(iShA,2) - CenterA(3) = CenterShell(iShA,3) - - TotAngMomA = TotAngMomShell(iShA) - nShellFunctionA = (TotAngMomA*TotAngMomA + 3*TotAngMomA + 2)/2 - allocate(ShellFunctionA(1:nShellFunctionA,1:3)) - call GenerateShell(TotAngMomA,nShellFunctionA,ShellFunctionA) - - KA = KShell(iShA) - - do iShFA=1,nShellFunctionA - - iBasA = iBasA + 1 - AngMomA(1) = ShellFunctionA(iShFA,1) - AngMomA(2) = ShellFunctionA(iShFA,2) - AngMomA(3) = ShellFunctionA(iShFA,3) - -!------------------------------------------------------------------------ -! Loops over shell B -!------------------------------------------------------------------------ - do iShB=1,iShA - - CenterB(1) = CenterShell(iShB,1) - CenterB(2) = CenterShell(iShB,2) - CenterB(3) = CenterShell(iShB,3) - - TotAngMomB = TotAngMomShell(iShB) - nShellFunctionB = (TotAngMomB*TotAngMomB + 3*TotAngMomB + 2)/2 - allocate(ShellFunctionB(1:nShellFunctionB,1:3)) - call GenerateShell(TotAngMomB,nShellFunctionB,ShellFunctionB) - - KB = KShell(iShB) - - do iShFB=1,nShellFunctionB - - iBasB = iBasB + 1 - AngMomB(1) = ShellFunctionB(iShFB,1) - AngMomB(2) = ShellFunctionB(iShFB,2) - AngMomB(3) = ShellFunctionB(iShFB,3) - -!------------------------------------------------------------------------ -! Loops over shell C -!------------------------------------------------------------------------ - do iShC=1,iShA - - CenterC(1) = CenterShell(iShC,1) - CenterC(2) = CenterShell(iShC,2) - CenterC(3) = CenterShell(iShC,3) - - TotAngMomC = TotAngMomShell(iShC) - nShellFunctionC = (TotAngMomC*TotAngMomC + 3*TotAngMomC + 2)/2 - allocate(ShellFunctionC(1:nShellFunctionC,1:3)) - call GenerateShell(TotAngMomC,nShellFunctionC,ShellFunctionC) - - KC = KShell(iShC) - - do iShFC=1,nShellFunctionC - - iBasC = iBasC + 1 - AngMomC(1) = ShellFunctionC(iShFC,1) - AngMomC(2) = ShellFunctionC(iShFC,2) - AngMomC(3) = ShellFunctionC(iShFC,3) - -!------------------------------------------------------------------------ -! Loops over shell D -!------------------------------------------------------------------------ - do iShD=1,iShC - - CenterD(1) = CenterShell(iShD,1) - CenterD(2) = CenterShell(iShD,2) - CenterD(3) = CenterShell(iShD,3) - - TotAngMomD = TotAngMomShell(iShD) - nShellFunctionD = (TotAngMomD*TotAngMomD + 3*TotAngMomD + 2)/2 - allocate(ShellFunctionD(1:nShellFunctionD,1:3)) - call GenerateShell(TotAngMomD,nShellFunctionD,ShellFunctionD) - - KD = KShell(iShD) - - do iShFD=1,nShellFunctionD - - iBasD = iBasD + 1 - AngMomD(1) = ShellFunctionD(iShFD,1) - AngMomD(2) = ShellFunctionD(iShFD,2) - AngMomD(3) = ShellFunctionD(iShFD,3) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - call cpu_time(start_cErf) - - cErf = 0d0 - - do iKA=1,KA - ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) - do iKB=1,KB - ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) - do iKC=1,KC - ExpC = ExpShell(iShC,iKC) - DC = DShell(iShC,iKC)*norm_coeff(ExpC,AngMomC) - do iKD=1,KD - ExpD = ExpShell(iShD,iKD) - DD = DShell(iShD,iKD)*norm_coeff(ExpD,AngMomD) - -! Erf module -! call ErfInt(debug,npErf,nSigpErf, & -! ExpS, & -! ExpA,CenterA,AngMomA, & -! ExpB,CenterB,AngMomB, & -! ExpC,CenterC,AngMomC, & -! ExpD,CenterD,AngMomD, & -! pErf) - -! cErf = cErf + DA*DB*DC*DD*pErf - - end do - end do - end do - end do - call cpu_time(end_cErf) - - ncErf = ncErf + 1 - if(abs(cErf) > 1d-15) then - nSigcErf = nSigcErf + 1 - t_cErf = end_cErf - start_cErf - write(41,'(F20.15,I6,I6,I6,I6)') & - cErf,iBasA,iBasB,iBasC,iBasD - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & - '(ab|erf(r)/r|cd) = ',cErf,iBasA,iBasB,iBasC,iBasD - end if - end if - -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionD) - end do - iBasD = 0 -!------------------------------------------------------------------------ -! End loops over shell D -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionC) - end do - iBasC = 0 -!------------------------------------------------------------------------ -! End loops over shell C -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB) - end do - iBasB = 0 -!------------------------------------------------------------------------ -! End loops over shell B -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA) - end do - iBasA = 0 -!------------------------------------------------------------------------ -! End loops over shell A -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=41) - -end subroutine Compute4eInt diff --git a/src/IntPak/ComputeKin.f90 b/src/IntPak/ComputeKin.f90 deleted file mode 100644 index 11e3b34..0000000 --- a/src/IntPak/ComputeKin.f90 +++ /dev/null @@ -1,167 +0,0 @@ -subroutine ComputeKin(debug,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - npKin,nSigpKin,ncKin,nSigcKin) - - -! Compute one-electron kinetic integrals - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: nShell - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -! Local variables - - integer :: KA,KB - double precision :: CenterA(3),CenterB(3) - integer :: TotAngMomA,TotAngMomB - integer :: AngMomA(3),AngMomB(3) - integer :: nShellFunctionA,nShellFunctionB - integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) - double precision :: ExpA,ExpB - double precision,allocatable :: DA,DB - double precision :: norm_coeff - - integer :: iBasA,iBasB - integer :: iShA,iShB - integer :: iShFA,iShFB - integer :: iKA,iKB - - double precision :: pKin,cKin - double precision :: start_cKin,end_cKin,t_cKin - -! Output variables - - integer,intent(out) :: npKin,nSigpKin,ncKin,nSigcKin - -! Compute one-electron integrals - - write(*,*) '****************************************' - write(*,*) ' Compute one-electron kinetic integrals ' - write(*,*) '****************************************' - write(*,*) - - npKin = 0 - nSigpKin = 0 - - ncKin = 0 - nSigcKin = 0 - - iBasA = 0 - iBasB = 0 - -! Open file to write down integrals - - open(unit=9,file='int/Kin.dat') - -!------------------------------------------------------------------------ -! Loops over shell A -!------------------------------------------------------------------------ - do iShA=1,nShell - - CenterA(1) = CenterShell(iShA,1) - CenterA(2) = CenterShell(iShA,2) - CenterA(3) = CenterShell(iShA,3) - - TotAngMomA = TotAngMomShell(iShA) - nShellFunctionA = (TotAngMomA*TotAngMomA + 3*TotAngMomA + 2)/2 - allocate(ShellFunctionA(1:nShellFunctionA,1:3)) - call GenerateShell(TotAngMomA,nShellFunctionA,ShellFunctionA) - - KA = KShell(iShA) - - do iShFA=1,nShellFunctionA - - iBasA = iBasA + 1 - AngMomA(1) = ShellFunctionA(iShFA,1) - AngMomA(2) = ShellFunctionA(iShFA,2) - AngMomA(3) = ShellFunctionA(iShFA,3) - -!------------------------------------------------------------------------ -! Loops over shell B -!------------------------------------------------------------------------ - do iShB=1,nShell - - CenterB(1) = CenterShell(iShB,1) - CenterB(2) = CenterShell(iShB,2) - CenterB(3) = CenterShell(iShB,3) - - TotAngMomB = TotAngMomShell(iShB) - nShellFunctionB = (TotAngMomB*TotAngMomB + 3*TotAngMomB + 2)/2 - allocate(ShellFunctionB(1:nShellFunctionB,1:3)) - call GenerateShell(TotAngMomB,nShellFunctionB,ShellFunctionB) - - KB = KShell(iShB) - - do iShFB=1,nShellFunctionB - - iBasB = iBasB + 1 - AngMomB(1) = ShellFunctionB(iShFB,1) - AngMomB(2) = ShellFunctionB(iShFB,2) - AngMomB(3) = ShellFunctionB(iShFB,3) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - call cpu_time(start_cKin) - - cKin = 0d0 - - do iKA=1,KA - ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) - do iKB=1,KB - ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) - - call KinInt(npKin,nSigpKin, & - ExpA,CenterA,AngMomA, & - ExpB,CenterB,AngMomB, & - pKin) - - cKin = cKin + DA*DB*pKin - - end do - end do - call cpu_time(end_cKin) - - ncKin = ncKin + 1 - if(abs(cKin) > 1d-15) then - nSigcKin = nSigcKin + 1 - t_cKin = end_cKin - start_cKin - write(9,'(I6,I6,F20.15)') iBasA,iBasB,cKin -! write(9,'(F20.15,I6,I6)') cKin,iBasA,iBasB - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|T|b) = ',cKin,iBasA,iBasB - end if - end if -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB) - end do - iBasB = 0 -!------------------------------------------------------------------------ -! End loops over shell B -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA) - end do - iBasA = 0 -!------------------------------------------------------------------------ -! End loops over shell A -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=9) - -end subroutine ComputeKin diff --git a/src/IntPak/ComputeNuc.f90 b/src/IntPak/ComputeNuc.f90 deleted file mode 100644 index fd402b0..0000000 --- a/src/IntPak/ComputeNuc.f90 +++ /dev/null @@ -1,190 +0,0 @@ -subroutine ComputeNuc(debug,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - NAtoms,ZNuc,XYZAtoms, & - npNuc,nSigpNuc,ncNuc,nSigcNuc) - - -! Compute electron repulsion integrals - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: nShell - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - integer :: NAtoms - double precision :: ZNuc(NAtoms),XYZAtoms(NAtoms,3) - -! Local variables - - integer :: KA,KB - double precision :: CenterA(3),CenterB(3),CenterC(3) - integer :: TotAngMomA,TotAngMomB - integer :: AngMomA(3),AngMomB(3) - integer :: nShellFunctionA,nShellFunctionB - integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) - double precision :: ExpA,ExpB,ZC - double precision,allocatable :: DA,DB - double precision :: norm_coeff - - integer :: iBasA,iBasB - integer :: iShA,iShB,iNucC - integer :: iShFA,iShFB - integer :: iKA,iKB - - double precision :: pNuc,cNuc - double precision :: start_cNuc,end_cNuc,t_cNuc - -! Output variables - - integer,intent(out) :: npNuc,nSigpNuc,ncNuc,nSigcNuc - -! Compute one-electron nuclear attraction integrals - - write(*,*) '***************************************************' - write(*,*) ' Compute one-electron nuclear attraction integrals ' - write(*,*) '***************************************************' - write(*,*) - - npNuc = 0 - nSigpNuc = 0 - - ncNuc = 0 - nSigcNuc = 0 - - iBasA = 0 - iBasB = 0 - iNucC = 0 - -! Open file to write down integrals - - open(unit=10,file='int/Nuc.dat') - -!------------------------------------------------------------------------ -! Loops over shell A -!------------------------------------------------------------------------ - do iShA=1,nShell - - CenterA(1) = CenterShell(iShA,1) - CenterA(2) = CenterShell(iShA,2) - CenterA(3) = CenterShell(iShA,3) - - TotAngMomA = TotAngMomShell(iShA) - nShellFunctionA = (TotAngMomA*TotAngMomA + 3*TotAngMomA + 2)/2 - allocate(ShellFunctionA(1:nShellFunctionA,1:3)) - call GenerateShell(TotAngMomA,nShellFunctionA,ShellFunctionA) - - KA = KShell(iShA) - - do iShFA=1,nShellFunctionA - - iBasA = iBasA + 1 - AngMomA(1) = ShellFunctionA(iShFA,1) - AngMomA(2) = ShellFunctionA(iShFA,2) - AngMomA(3) = ShellFunctionA(iShFA,3) - -!------------------------------------------------------------------------ -! Loops over shell B -!------------------------------------------------------------------------ - do iShB=1,nShell - - CenterB(1) = CenterShell(iShB,1) - CenterB(2) = CenterShell(iShB,2) - CenterB(3) = CenterShell(iShB,3) - - TotAngMomB = TotAngMomShell(iShB) - nShellFunctionB = (TotAngMomB*TotAngMomB + 3*TotAngMomB + 2)/2 - allocate(ShellFunctionB(1:nShellFunctionB,1:3)) - call GenerateShell(TotAngMomB,nShellFunctionB,ShellFunctionB) - - KB = KShell(iShB) - - do iShFB=1,nShellFunctionB - - iBasB = iBasB + 1 - AngMomB(1) = ShellFunctionB(iShFB,1) - AngMomB(2) = ShellFunctionB(iShFB,2) - AngMomB(3) = ShellFunctionB(iShFB,3) - -!------------------------------------------------------------------------ -! Loops over nuclear centers -!------------------------------------------------------------------------ - call cpu_time(start_cNuc) - - cNuc = 0d0 - - do iNucC=1,NAtoms - - CenterC(1) = XYZAtoms(iNucC,1) - CenterC(2) = XYZAtoms(iNucC,2) - CenterC(3) = XYZAtoms(iNucC,3) - - ZC = ZNuc(iNucC) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - - do iKA=1,KA - ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) - do iKB=1,KB - ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) - - call NucInt(debug,npNuc,nSigpNuc, & - ExpA,CenterA,AngMomA, & - ExpB,CenterB,AngMomB, & - CenterC, & - pNuc) - - cNuc = cNuc - DA*DB*ZC*pNuc - - end do - end do -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - call cpu_time(end_cNuc) -!------------------------------------------------------------------------ -! End loops over nuclear centers C -!------------------------------------------------------------------------ - - ncNuc = ncNuc + 1 - if(abs(cNuc) > 1d-15) then - nSigcNuc = nSigcNuc + 1 - t_cNuc = end_cNuc - start_cNuc - write(10,'(I6,I6,F20.15)') iBasA,iBasB,cNuc -! write(10,'(F20.15,I6,I6)') cNuc,iBasA,iBasB - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|V|b) = ',cNuc,iBasA,iBasB - write(*,*) - end if - end if - - end do - deallocate(ShellFunctionB) - end do - iBasB = 0 -!------------------------------------------------------------------------ -! End loops over shell B -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA) - end do - iBasA = 0 -!------------------------------------------------------------------------ -! End loops over shell A -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=10) - -end subroutine ComputeNuc diff --git a/src/IntPak/ComputeOv.f90 b/src/IntPak/ComputeOv.f90 deleted file mode 100644 index 5879729..0000000 --- a/src/IntPak/ComputeOv.f90 +++ /dev/null @@ -1,171 +0,0 @@ -subroutine ComputeOv(debug,nBas,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - npOv,nSigpOv,ncOv,nSigcOv,S) - - -! Compute one-electron overlap integrals - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: nBas,nShell - double precision,intent(in) :: CenterShell(maxShell,3) - integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -! Local variables - - integer :: KA,KB - double precision :: CenterA(3),CenterB(3) - integer :: TotAngMomA,TotAngMomB - integer :: AngMomA(3),AngMomB(3) - integer :: nShellFunctionA,nShellFunctionB - integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) - double precision :: ExpA,ExpB - double precision,allocatable :: DA,DB - double precision :: norm_coeff - - integer :: iBasA,iBasB - integer :: iShA,iShB - integer :: iShFA,iShFB - integer :: iKA,iKB - - double precision :: pOv,cOv - double precision :: start_cOv,end_cOv,t_cOv - -! Output variables - - integer,intent(out) :: npOv,nSigpOv,ncOv,nSigcOv - double precision,intent(out) :: S(nBas,nBas) - - -! Compute one-electron integrals - - write(*,*) '****************************************' - write(*,*) ' Compute one-electron overlap integrals ' - write(*,*) '****************************************' - write(*,*) - - npOv = 0 - nSigpOv = 0 - - ncOv = 0 - nSigcOv = 0 - - iBasA = 0 - iBasB = 0 - -! Open file to write down integrals - - open(unit=8,file='int/Ov.dat') - -!------------------------------------------------------------------------ -! Loops over shell A -!------------------------------------------------------------------------ - do iShA=1,nShell - - CenterA(1) = CenterShell(iShA,1) - CenterA(2) = CenterShell(iShA,2) - CenterA(3) = CenterShell(iShA,3) - - TotAngMomA = TotAngMomShell(iShA) - nShellFunctionA = (TotAngMomA*TotAngMomA + 3*TotAngMomA + 2)/2 - allocate(ShellFunctionA(1:nShellFunctionA,1:3)) - call GenerateShell(TotAngMomA,nShellFunctionA,ShellFunctionA) - - KA = KShell(iShA) - - do iShFA=1,nShellFunctionA - - iBasA = iBasA + 1 - AngMomA(1) = ShellFunctionA(iShFA,1) - AngMomA(2) = ShellFunctionA(iShFA,2) - AngMomA(3) = ShellFunctionA(iShFA,3) - -!------------------------------------------------------------------------ -! Loops over shell B -!------------------------------------------------------------------------ - do iShB=1,nShell - - CenterB(1) = CenterShell(iShB,1) - CenterB(2) = CenterShell(iShB,2) - CenterB(3) = CenterShell(iShB,3) - - TotAngMomB = TotAngMomShell(iShB) - nShellFunctionB = (TotAngMomB*TotAngMomB + 3*TotAngMomB + 2)/2 - allocate(ShellFunctionB(1:nShellFunctionB,1:3)) - call GenerateShell(TotAngMomB,nShellFunctionB,ShellFunctionB) - - KB = KShell(iShB) - - do iShFB=1,nShellFunctionB - - iBasB = iBasB + 1 - AngMomB(1) = ShellFunctionB(iShFB,1) - AngMomB(2) = ShellFunctionB(iShFB,2) - AngMomB(3) = ShellFunctionB(iShFB,3) - -!------------------------------------------------------------------------ -! Loops over contraction degrees -!------------------------------------------------------------------------- - call cpu_time(start_cOv) - - cOv = 0d0 - - do iKA=1,KA - ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) - do iKB=1,KB - ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) - - call OvInt(npOv,nSigpOv, & - ExpA,CenterA,AngMomA, & - ExpB,CenterB,AngMomB, & - pOv) - - cOv = cOv + DA*DB*pOv - - end do - end do - call cpu_time(end_cOv) - - ncOv = ncOv + 1 - S(iBasA,iBasB) = cOv - if(abs(cOv) > 1d-15) then - nSigcOv = nSigcOv + 1 - t_cOv = end_cOv - start_cOv - write(8,'(I6,I6,F20.15)') iBasA,iBasB,cOv -! write(8,'(F20.15,I6,I6)') cOv,iBasA,iBasB - if(debug) then - write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|b) = ',cOv,iBasA,iBasB - end if - end if - -!------------------------------------------------------------------------ -! End loops over contraction degrees -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionB) - end do - iBasB = 0 -!------------------------------------------------------------------------ -! End loops over shell B -!------------------------------------------------------------------------ - end do - deallocate(ShellFunctionA) - end do - iBasA = 0 -!------------------------------------------------------------------------ -! End loops over shell A -!------------------------------------------------------------------------ - write(*,*) - -! Close files to write down integrals - - close(unit=8) - -end subroutine ComputeOv diff --git a/src/IntPak/FormVRR3e.f90 b/src/IntPak/FormVRR3e.f90 deleted file mode 100644 index 97904b4..0000000 --- a/src/IntPak/FormVRR3e.f90 +++ /dev/null @@ -1,174 +0,0 @@ -subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) - -! Form stuff we need... - - implicit none - include 'parameters.h' - - -! Input variables - - double precision,intent(in) :: ExpZ(3),ExpG(3,3) - double precision,intent(in) :: CenterZ(3,3) - -! Local variables - - integer :: i,j,k,l - double precision :: ZetaMat(3,3) - double precision :: CMat(3,3),GMat(3,3) - double precision :: Delta0Mat(3,3),Delta1Mat(3,3) - double precision :: InvDelta0Mat(3,3),InvDelta1Mat(3,3) - double precision :: CenterY(3,3,3) - double precision :: YMat(3,3),Y2Mat(3,3) - double precision :: DYMat(3,3,3),D2YMat(3,3,3,3) - double precision :: D0Mat(3,3),D1Mat(3,3) - - double precision :: KappaCross - -! Output variables - - double precision,intent(out) :: DY0(3),DY1(3),D2Y0(3,3),D2Y1(3,3) - double precision,intent(out) :: delta0,delta1,Y0,Y1 - -! Initalize arrays - - ZetaMat = 0d0 - CMat = 0d0 - GMat = 0d0 - YMat = 0d0 - Y2Mat = 0d0 - D0Mat = 0d0 - D1Mat = 0d0 - -! Form the zeta matrix Eq. (15a) - - do i=1,3 - ZetaMat(i,i) = ExpZ(i) - end do - -! print*,'Zeta' -! call matout(3,3,ZetaMat) - -! Form the C matrix Eq. (15a) - - CMat(1,1) = 1d0 - CMat(2,2) = 1d0 - CMat(1,2) = -1d0 - CMat(2,1) = -1d0 - -! print*,'C' -! call matout(3,3,CMat) - -! Form the G matrix Eq. (15b) - - do i=1,3 - do j=1,i-1 - GMat(i,j) = - ExpG(j,i) - end do - do j=i+1,3 - GMat(i,j) = - ExpG(i,j) - end do - end do - - do i=1,3 - do j=1,i-1 - GMat(i,i) = GMat(i,i) + ExpG(j,i) - end do - do j=i+1,3 - GMat(i,i) = GMat(i,i) + ExpG(i,j) - end do - end do - -! print*,'G' -! call matout(3,3,GMat) - -! Form the Y and Y^2 matrices Eq. (16b) - - do i=1,3 - do j=i+1,3 - do k=1,3 - CenterY(i,j,k) = CenterZ(i,k) - CenterZ(j,k) - Y2Mat(i,j) = Y2Mat(i,j) + CenterY(i,j,k)**2 - end do - YMat(i,j) = sqrt(Y2Mat(i,j)) - end do - end do - -! print*,'Y' -! call matout(3,3,YMat) - -! print*,'Y2' -! call matout(3,3,Y2Mat) - -! Form the delta0 and delta1 matrices Eq. (14) - - do i=1,3 - do j=1,3 - Delta0Mat(i,j) = ZetaMat(i,j) + GMat(i,j) - Delta1Mat(i,j) = Delta0Mat(i,j) + CMat(i,j) - end do - end do - -! Form the DY and D2Y matrices - - do i=1,3 - do j=1,3 - do k=1,3 - DYMat(i,j,k) = KappaCross(i,j,k)*YMat(j,k)/ExpZ(i) - do l=1,3 - D2YMat(i,j,k,l) = 0.5d0*KappaCross(i,k,l)*KappaCross(j,k,l)/(ExpZ(i)*ExpZ(j)) - end do - end do - end do - end do - -! Compute the inverse of the Delta0 and Delta1 matrices - -! InvDelta0Mat = Delta0Mat -! InvDelta1Mat = Delta1Mat - do i=1,3 - do j=1,3 - InvDelta0Mat(i,j) = Delta0Mat(i,j) - InvDelta1Mat(i,j) = Delta1Mat(i,j) - end do - end do -! call amove(3,3,Delta0Mat,InvDelta0Mat) -! call amove(3,3,Delta1Mat,InvDelta1Mat) - - call CalcInv3(InvDelta0Mat,delta0) - call CalcInv3(InvDelta1Mat,delta1) - -! call matout(3,3,InvDelta0Mat) -! call matout(3,3,InvDelta1Mat) -! print*, 'delta0,delta1 = ',delta0,delta1 - -! Form the Delta matrix Eq. (16a) - - do i=1,3 - do j=1,3 - do k=1,3 - do l=1,3 - D0Mat(i,j) = D0Mat(i,k) + ZetaMat(i,k)*InvDelta0Mat(k,l)*ZetaMat(l,j) - D1Mat(i,j) = D1Mat(i,k) + ZetaMat(i,k)*InvDelta1Mat(k,l)*ZetaMat(l,j) - end do - end do - end do - end do - -! Form the derivative matrices - - do i=1,3 - call CalcTrAB(3,D0Mat,D2YMat,DY0(i)) - call CalcTrAB(3,D1Mat,D2YMat,DY1(i)) - do j=1,3 - call CalcTrAB(3,D0Mat,D2YMat,D2Y0(i,j)) - call CalcTrAB(3,D1Mat,D2YMat,D2Y1(i,j)) - end do - end do - -! Compute Y0 and Y1 - - call CalcTrAB(3,D0Mat,Y2Mat,Y0) - call CalcTrAB(3,D1Mat,Y2Mat,Y1) - -end subroutine FormVRR3e diff --git a/src/IntPak/G2eInt.f90 b/src/IntPak/G2eInt.f90 deleted file mode 100644 index aa6e90d..0000000 --- a/src/IntPak/G2eInt.f90 +++ /dev/null @@ -1,140 +0,0 @@ -function G2eInt(debug,iType, & - ExpG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet) - -! Compute recursively the primitive two-electron integral [ab|cd] - - implicit none - include 'parameters.h' - - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: iType - double precision,intent(in) :: ExpBra(2),ExpKet(2) - double precision,intent(in) :: ExpG - double precision,intent(in) :: CenterBra(2,3),CenterKet(2,3) - integer,intent(in) :: AngMomBra(2,3),AngMomKet(2,3) - -! Local variables - - integer :: TotAngMomBra(3),TotAngMomKet(3) - double precision :: ExpZi(2),ExpY(2,2) - double precision :: CenterZ(2,3),CenterAB(2,3),CenterZA(2,3),CenterY(2,2,3) - double precision :: NormABSq(2),NormYSq(2,2) - double precision :: GAB(2) - double precision,allocatable :: Om(:) - double precision :: fG - double precision :: HRR2e,VRR2e - double precision :: a1a2b1b2 - - integer :: i,j,k,maxm - double precision :: start_Om,finish_Om,start_RR,finish_RR,t_Om,t_RR - -! Output variables - double precision :: G2eInt - -! Pre-computed shell-pair quantities - - do i=1,2 - ExpZi(i) = 1d0/(ExpBra(i) + ExpKet(i)) - end do - - NormABSq = 0d0 - do j=1,3 - do i=1,2 - CenterZ(i,j) = (ExpBra(i)*CenterBra(i,j) + ExpKet(i)*CenterKet(i,j))*ExpZi(i) - CenterAB(i,j) = CenterBra(i,j) - CenterKet(i,j) - CenterZA(i,j) = CenterZ(i,j) - CenterBra(i,j) - NormABSq(i) = NormABSq(i) + CenterAB(i,j)**2 - end do - end do - - do i=1,2 - GAB(i) = (pi*ExpZi(i))**(1.5d0)*exp(-ExpBra(i)*ExpKet(i)*NormABSq(i)*ExpZi(i)) - end do - -! Pre-computed shell-quartet quantities - - do i=1,2 - do j=1,2 - ExpY(i,j) = 1d0/(ExpZi(i) + ExpZi(j)) - end do - end do - - do i=1,2 - do j=1,2 - NormYSq(i,j) = 0d0 - do k=1,3 - CenterY(i,j,k) = CenterZ(i,k) - CenterZ(j,k) - NormYSq(i,j) = NormYSq(i,j) + CenterY(i,j,k)**2 - end do - end do - end do - -! fG = (ExpZ(1)*ExpZ(2)*ExpG)/(ExpZ(1)*ExpZ(2) + ExpZ(1)*ExpG + ExpZ(2)*ExpG) - fG = 1d0/(ExpZi(1) + 1d0/ExpG + ExpZi(2)) - -! Total angular momemtum - - maxm = 0 - do i=1,2 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - maxm = maxm + TotAngMomBra(i) + TotAngMomKet(i) - end do - -! Pre-compute (00|00)^m - - allocate(Om(0:maxm)) - call cpu_time(start_Om) - - if(iType == 1) then - call CalcOmERI(maxm,ExpY(1,2),NormYSq(1,2),Om) - elseif(iType == 3) then - call CalcOmYuk(maxm,ExpG,ExpY(1,2),fG,NormYSq(1,2),Om) - elseif(iType == 4) then - call CalcOmErf(maxm,ExpY(1,2),fG,NormYSq(1,2),Om) - end if - - call cpu_time(finish_Om) - -! Print (00|00)^m - - if(debug) then - write(*,*) '(00|00)^m' - do i=0,maxm - write(*,*) i,Om(i) - end do - write(*,*) - end if - -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) - - if(TotAngMomKet(1) == 0 .and. TotAngMomKet(2) == 0) then - if(TotAngMomBra(1) == 0 .and. TotAngMomBra(2) == 0) then - a1a2b1b2 = Om(0) - else - a1a2b1b2 = VRR2e(0,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - end if - else - a1a2b1b2 = HRR2e(AngMomBra,AngMomKet,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) - end if - - call cpu_time(finish_RR) - -! Timings - - t_Om = finish_Om - start_Om - t_RR = finish_RR - start_RR - -! Print result - - G2eInt = GAB(1)*GAB(2)*a1a2b1b2 - -end function G2eInt diff --git a/src/IntPak/G3eInt.f90 b/src/IntPak/G3eInt.f90 deleted file mode 100644 index d6926f5..0000000 --- a/src/IntPak/G3eInt.f90 +++ /dev/null @@ -1,124 +0,0 @@ -function G3eInt(debug,iType, & - ExpG13,ExpG23, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet) - -! Compute two-electron integrals over the Yukawa operator - - implicit none - include 'parameters.h' - - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: iType - double precision,intent(in) :: ExpG13,ExpG23 - double precision,intent(in) :: ExpBra(3),ExpKet(3) - double precision,intent(in) :: CenterBra(3,3),CenterKet(3,3) - integer,intent(in) :: AngMomBra(3,3),AngMomKet(3,3) - -! Local variables - - double precision :: ExpG(3,3) - integer :: TotAngMomBra(3),TotAngMomKet(3) - double precision :: ExpZ(3) - double precision :: CenterZ(3,3),CenterAB(3,3),CenterZA(3,3) - double precision :: NormABSq(3) - double precision :: GAB(3) - double precision,allocatable :: Om(:) - double precision :: HRR3e,VRR3e - - double precision :: DY0(3),DY1(3),D2Y0(3,3),D2Y1(3,3) - double precision :: delta0,delta1,Y0,Y1 - - integer :: i,j,maxm - double precision :: start_Om,finish_Om,t_Om,start_RR,finish_RR,t_RR - double precision :: a1a2a3b1b2b3 - -! Output variables - double precision :: G3eInt - -! Gaussian geminal exponents - - ExpG = 0d0 - ExpG(1,3) = ExpG13 - ExpG(2,3) = ExpG23 - -! Pre-computed quantities for shell-pair - - do i=1,3 - ExpZ(i) = ExpBra(i) + ExpKet(i) - end do - - NormABSq = 0d0 - do i=1,3 - do j=1,3 - CenterZ(i,j) = (ExpBra(i)*CenterBra(i,j) + ExpKet(i)*CenterKet(i,j))/ExpZ(i) - CenterAB(i,j) = CenterBra(i,j) - CenterKet(i,j) - CenterZA(i,j) = CenterZ(i,j) - CenterBra(i,j) - NormABSq(i) = NormABSq(i) + CenterAB(i,j)**2 - end do - end do - - do i=1,3 - GAB(i) = (pi/ExpZ(i))**(1.5d0)*exp(-ExpBra(i)*ExpKet(i)*NormABSq(i)/ExpZ(i)) - end do - -! Pre-computed shell-sextet quantities - - call FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) - -! Total angular momemtum - - maxm = 0 - do i=1,3 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - maxm = maxm + TotAngMomBra(i) + TotAngMomKet(i) - end do - -! Pre-compute (000|000)^m - - allocate(Om(0:maxm)) - call cpu_time(start_Om) - call CalcOm3e(maxm,delta0,delta1,Y0,Y1,Om) - call cpu_time(finish_Om) - -! Print (000|000)^m - - if(.false.) then - write(*,*) '(000|000)^m' - do i=0,maxm - write(*,*) i,Om(i) - end do - write(*,*) - end if - -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) - if(TotAngMomKet(1) == 0 .and. TotAngMomKet(2) == 0 .and. TotAngMomKet(3) == 0) then - if(TotAngMomBra(1) == 0 .and. TotAngMomBra(2) == 0 .and. TotAngMomBra(3) == 0) then - a1a2a3b1b2b3 = Om(0) - else - a1a2a3b1b2b3 = VRR3e(0,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - else - a1a2a3b1b2b3 = HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - - - call cpu_time(finish_RR) - -! Timings - - t_Om = finish_Om - start_Om - t_RR = finish_RR - start_RR - -! Print result - - G3eInt = GAB(1)*GAB(2)*GAB(3)*a1a2a3b1b2b3 - -end function G3eInt diff --git a/src/IntPak/GF12Int.f90 b/src/IntPak/GF12Int.f90 deleted file mode 100644 index ac1f555..0000000 --- a/src/IntPak/GF12Int.f90 +++ /dev/null @@ -1,107 +0,0 @@ -function GF12Int(ExpG,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,ExpC,CenterC,AngMomC,ExpD,CenterD,AngMomD) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - -! Input variables - - double precision,intent(in) :: ExpG - double precision,intent(in) :: ExpA,ExpB,ExpC,ExpD - double precision,intent(in) :: CenterA(3),CenterB(3),CenterC(3),CenterD(3) - integer,intent(in) :: AngMomA(3),AngMomB(3),AngMomC(3),AngMomD(3) - - -! Local variables - - double precision :: ExpAi,ExpBi,ExpCi,ExpDi,ExpGi - double precision :: ExpP,ExpQ,ExpPi,ExpQi,ExpPGQi - double precision :: CenterP(3),CenterQ(3),CenterAB(3),CenterCD(3),CenterPQSq(3),CenterRA(3),CenterRC(3) - double precision :: NormABSq,NormCDSq - double precision :: GAB,GCD - double precision :: fP,fG,fQ,gP,gG,gQ - double precision :: HRRF12 - - integer :: i - double precision :: pi - double precision :: start_RR,finish_RR,t_RR - double precision :: Gabcd(3) - -! Output variables - double precision :: GF12Int - - pi = 4d0*atan(1d0) - -! Pre-computed shell quantities - - ExpAi = 1d0/ExpA - ExpBi = 1d0/ExpB - ExpCi = 1d0/ExpC - ExpDi = 1d0/ExpD - ExpGi = 1d0/ExpG - -! Pre-computed quantities for shell-pair AB - - ExpP = ExpA + ExpB - ExpPi = 1d0/ExpP - - NormABSq = 0d0 - Do i=1,3 - CenterP(i) = (ExpA*CenterA(i) + ExpB*CenterB(i))*ExpPi - CenterAB(i) = CenterA(i) - CenterB(i) - NormABSq = NormABSq + CenterAB(i)**2 - Enddo - - GAB = (pi*ExpPi)**(1.5d0)*exp(-NormABSq/(ExpAi+ExpBi)) - -! Pre-computed quantities for shell-pair CD - - ExpQ = ExpC + ExpD - ExpQi = 1d0/ExpQ - - NormCDSq = 0d0 - Do i=1,3 - CenterQ(i) = (ExpC*CenterC(i) + ExpD*CenterD(i))*ExpQi - CenterCD(i) = CenterC(i) - CenterD(i) - NormCDSq = NormCDSq + CenterCD(i)**2 - Enddo - - GCD = (pi*ExpQi)**(1.5d0)*exp(-NormCDSq/(ExpCi+ExpDi)) - -! Pre-computed shell-quartet quantities - - ExpPGQi = ExpPi + ExpGi + ExpQi - - Do i=1,3 - CenterPQSq(i) = (CenterP(i) - CenterQ(i))**2 - Enddo - - fP = ExpPi/ExpPGQi - fG = ExpGi/ExpPGQi - fQ = ExpQi/ExpPGQi - - gP = (1d0 - fP)*0.5d0*ExpPi - gG = fP*0.5d0*expQi - gQ = (1d0 - fQ)*0.5d0*ExpQi - - do i=1,3 - CenterRA(i) = CenterP(i) - CenterA(i) + fP*(CenterQ(i) - CenterP(i)) - CenterRC(i) = CenterQ(i) - CenterC(i) + fQ*(CenterP(i) - CenterQ(i)) - end do -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) -! Loop over cartesian directions - Do i=1,3 - Gabcd(i) = HRRF12(AngMomA(i),AngMomB(i),AngMomC(i),AngMomD(i),fG,gP,gG,gQ,ExpPGQi, & - CenterPQSq(i),CenterRA(i),CenterRC(i),CenterAB(i),CenterCD(i)) - Enddo - call cpu_time(finish_RR) - -! Print result - - GF12Int = GAB*GCD*Gabcd(1)*Gabcd(2)*Gabcd(3) - t_RR = finish_RR - start_RR - -end function GF12Int diff --git a/src/IntPak/GenerateShell.f90 b/src/IntPak/GenerateShell.f90 deleted file mode 100644 index 7be3c00..0000000 --- a/src/IntPak/GenerateShell.f90 +++ /dev/null @@ -1,30 +0,0 @@ -subroutine GenerateShell(atot,nShellFunction,ShellFunction) - - implicit none - -! Input variables - - integer,intent(in) :: atot,nShellFunction - -! Local variables - - integer :: ax,ay,az,ia - -! Output variables - - integer,intent(out) :: ShellFunction(nShellFunction,3) - - ia = 0 - do ax=atot,0,-1 - do az=0,atot - ay = atot - ax - az - if(ay >= 0) then - ia = ia + 1 - ShellFunction(ia,1) = ax - ShellFunction(ia,2) = ay - ShellFunction(ia,3) = az - end if - end do - end do - -end subroutine GenerateShell diff --git a/src/IntPak/HRR2e.f90 b/src/IntPak/HRR2e.f90 deleted file mode 100644 index 1ab9126..0000000 --- a/src/IntPak/HRR2e.f90 +++ /dev/null @@ -1,101 +0,0 @@ -recursive function HRR2e(AngMomBra,AngMomKet, & - maxm,Om,ExpZi,ExpY, & - CenterAB,CenterZA,CenterY) & - result(a1a2b1b2) - -! Horintal recurrence relations for two-electron integrals - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: AngMomBra(2,3),AngMomKet(2,3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm),ExpZi(2),ExpY(2,2) - double precision,intent(in) :: CenterAB(2,3),CenterZA(2,3),CenterY(2,2,3) - -! Local variables - - logical :: NegAngMomKet(2) - integer :: TotAngMomBra(2),TotAngMomKet(2) - integer :: a1p(2,3),b1m(2,3),a2p(2,3),b2m(2,3) - integer :: i,j,xyz - double precision :: VRR2e - -! Output variables - - double precision :: a1a2b1b2 - - do i=1,2 - NegAngMomKet(i) = AngMomKet(i,1) < 0 .or. AngMomKet(i,2) < 0 .or. AngMomKet(i,3) < 0 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - end do - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ -! if(NegAngMomKet(1) .or. NegAngMomKet(2)) then -! a1a2b1b2 = 0d0 -!------------------------------------------------------------------------ -! 1st and 2nd vertical recurrence relations: -!------------------------------------------------------------------------ -! elseif(TotAngMomKet(1) == 0 .and. TotAngMomKet(2) == 0) then - if(TotAngMomKet(1) == 0 .and. TotAngMomKet(2) == 0) then - a1a2b1b2 = VRR2e(0,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) -!------------------------------------------------------------------------ -! 1st horizontal recurrence relation (2 terms): -!------------------------------------------------------------------------ - elseif(TotAngMomKet(2) == 0) then - do i=1,2 - do j=1,3 - a1p(i,j) = AngMomBra(i,j) - b1m(i,j) = AngMomKet(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomKet(1,1) > 0) then - xyz = 1 - elseif(AngMomKet(1,2) > 0) then - xyz = 2 - elseif(AngMomKet(1,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRR2e!' - end if -! End loop over cartesian directions - a1p(1,xyz) = a1p(1,xyz) + 1 - b1m(1,xyz) = b1m(1,xyz) - 1 - a1a2b1b2 = HRR2e(a1p,b1m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) & - + CenterAB(1,xyz)*HRR2e(AngMomBra,b1m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) -!------------------------------------------------------------------------ -! 2nd horizontal recurrence relation (2 terms): -!------------------------------------------------------------------------ - else - do i=1,2 - do j=1,3 - a2p(i,j) = AngMomBra(i,j) - b2m(i,j) = AngMomKet(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomKet(2,1) > 0) then - xyz = 1 - elseif(AngMomKet(2,2) > 0) then - xyz = 2 - elseif(AngMomKet(2,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRR2e!' - end if -! End loop over cartesian directions - a2p(2,xyz) = a2p(2,xyz) + 1 - b2m(2,xyz) = b2m(2,xyz) - 1 - a1a2b1b2 = HRR2e(a2p,b2m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) & - + CenterAB(2,xyz)*HRR2e(AngMomBra,b2m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) - end if - -end function HRR2e diff --git a/src/IntPak/HRR3e.f90 b/src/IntPak/HRR3e.f90 deleted file mode 100644 index 6d597ca..0000000 --- a/src/IntPak/HRR3e.f90 +++ /dev/null @@ -1,128 +0,0 @@ -recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) & - result(a1a2a3b1b2b3) - -! Horizontal recurrence relations for three-electron integrals - - implicit none - include 'parameters.h' - - -! Input variables - - integer,intent(in) :: AngMomBra(3,3),AngMomKet(3,3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm),ExpZ(3),CenterAB(3,3),CenterZA(3,3) - double precision,intent(in) :: DY0(3),DY1(3),D2Y0(3,3),D2Y1(3,3) - -! Local variables - - logical :: NegAngMomKet(3) - integer :: TotAngMomBra(3),TotAngMomKet(3) - integer :: a1p(3,3),b1m(3,3),a2p(3,3),b2m(3,3),a3p(3,3),b3m(3,3) - integer :: i,j,xyz - double precision :: VRR3e - -! Output variables - - double precision :: a1a2a3b1b2b3 - - do i=1,3 - NegAngMomKet(i) = AngMomKet(i,1) < 0 .or. AngMomKet(i,2) < 0 .or. AngMomKet(i,3) < 0 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - end do - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ - if(NegAngMomKet(1) .or. NegAngMomKet(2) .or. NegAngMomKet(3)) then - a1a2a3b1b2b3 = 0d0 -!------------------------------------------------------------------------ -! 1st and 2nd vertical recurrence relations: -!------------------------------------------------------------------------ - elseif(TotAngMomKet(1) == 0 .and. TotAngMomKet(2) == 0 .and. TotAngMomKet(3) == 0) then - a1a2a3b1b2b3 = VRR3e(0,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) -!------------------------------------------------------------------------ -! 1st horizontal recurrence relation (2 terms): -!------------------------------------------------------------------------ - elseif(TotAngMomKet(2) == 0 .and. TotAngMomKet(3) == 0) then - do i=1,3 - do j=1,3 - a1p(i,j) = AngMomBra(i,j) - b1m(i,j) = AngMomKet(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomKet(1,1) > 0) then - xyz = 1 - elseif(AngMomKet(1,2) > 0) then - xyz = 2 - elseif(AngMomKet(1,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRR3e!' - end if -! End loop over cartesian directions - a1p(1,xyz) = a1p(1,xyz) + 1 - b1m(1,xyz) = b1m(1,xyz) - 1 - a1a2a3b1b2b3 = HRR3e(a1p,b1m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + CenterAB(1,xyz)* & - HRR3e(AngMomBra,b1m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) -!------------------------------------------------------------------------ -! 2nd horizontal recurrence relation (2 terms): -!------------------------------------------------------------------------ - elseif(TotAngMomKet(3) == 0) then - do i=1,3 - do j=1,3 - a2p(i,j) = AngMomBra(i,j) - b2m(i,j) = AngMomKet(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomKet(2,1) > 0) then - xyz = 1 - elseif(AngMomKet(2,2) > 0) then - xyz = 2 - elseif(AngMomKet(2,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRR3e!' - end if -! End loop over cartesian directions - a2p(2,xyz) = a2p(2,xyz) + 1 - b2m(2,xyz) = b2m(2,xyz) - 1 - a1a2a3b1b2b3 = HRR3e(a2p,b2m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + CenterAB(2,xyz)* & - HRR3e(AngMomBra,b2m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) -!------------------------------------------------------------------------ -! 3rd horizontal recurrence relation (2 terms): -!------------------------------------------------------------------------ - else - do i=1,3 - do j=1,3 - a3p(i,j) = AngMomBra(i,j) - b3m(i,j) = AngMomKet(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomKet(3,1) > 0) then - xyz = 1 - elseif(AngMomKet(3,2) > 0) then - xyz = 2 - elseif(AngMomKet(3,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRR3e!' - end if -! End loop over cartesian directions - a3p(3,xyz) = a3p(3,xyz) + 1 - b3m(3,xyz) = b3m(3,xyz) - 1 - a1a2a3b1b2b3 = HRR3e(a3p,b3m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + CenterAB(3,xyz)* & - HRR3e(AngMomBra,b3m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - -end function HRR3e diff --git a/src/IntPak/HRRF12.f90 b/src/IntPak/HRRF12.f90 deleted file mode 100644 index 6ad2efe..0000000 --- a/src/IntPak/HRRF12.f90 +++ /dev/null @@ -1,40 +0,0 @@ -recursive function HRRF12(AngMomA,AngMomB,AngMomC,AngMomD,fG,gP,gG,gQ,ExpPGQi, & - CenterPQSq,CenterRA,CenterRC,CenterAB,CenterCD) & - result(Gabcd) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - -! Input variables - integer,intent(in) :: AngMomA,AngMomB,AngMomC,AngMomD - double precision,intent(in) :: ExpPGQi - double precision,intent(in) :: fG,gP,gG,gQ - double precision,intent(in) :: CenterPQSq,CenterRA,CenterRC - double precision,intent(in) :: CenterAB,CenterCD - -! Local variables - double precision :: VRRF12 - double precision :: Gabcd - - If(AngMomB < 0 .or. AngMomD < 0) then - Gabcd = 0d0 - Else - If(AngMomB == 0 .and. AngMomD == 0) then - Gabcd = VRRF12(AngMomA,AngMomC,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) - Else - If(AngMomD == 0) then - Gabcd = HRRF12(AngMomA+1,AngMomB-1,AngMomC,AngMomD,fG,gP,gG,gQ,ExpPGQi, & - CenterPQSq,CenterRA,CenterRC,CenterAB,CenterCD) & - + CenterAB*HRRF12(AngMomA,AngMomB-1,AngMomC,AngMomD,fG,gP,gG,gQ, & - ExpPGQi,CenterPQSq,CenterRA,CenterRC,CenterAB,CenterCD) - Else - Gabcd = HRRF12(AngMomA,AngMomB,AngMomC+1,AngMomD-1,fG,gP,gG,gQ,ExpPGQi, & - CenterPQSq,CenterRA,CenterRC,CenterAB,CenterCD) & - + CenterCD*HRRF12(AngMomA,AngMomB,AngMomC,AngMomD-1,fG,gP,gG,gQ, & - ExpPGQi,CenterPQSq,CenterRA,CenterRC,CenterAB,CenterCD) - EndIf - EndIf - EndIf - -end function HRRF12 diff --git a/src/IntPak/HRRNuc.f90 b/src/IntPak/HRRNuc.f90 deleted file mode 100644 index 50039ee..0000000 --- a/src/IntPak/HRRNuc.f90 +++ /dev/null @@ -1,71 +0,0 @@ -recursive function HRRNuc(AngMomA,AngMomB,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - result(Gab) - -! Horizontal recurrence relation for one-electron nuclear attraction integrals - - implicit none - -! Input variables - - integer,intent(in) :: AngMomA(3),AngMomB(3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm) - double precision,intent(in) :: ExpPi - double precision,intent(in) :: CenterAB(3),CenterPA(3),CenterPC(3) - -! Local variables - - logical :: NegAngMomB - integer :: TotAngMomA,TotAngMomB - integer :: xyz,ap(3),bm(3) - integer :: i - double precision :: VRRNuc - -! Output variables - - double precision :: Gab - - NegAngMomB = AngMomB(1) < 0 .or. AngMomB(2) < 0 .or. AngMomB(3) < 0 - - TotAngMomA = AngMomA(1) + AngMomA(2) + AngMomA(3) - TotAngMomB = AngMomB(1) + AngMomB(2) + AngMomB(3) - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ - if(NegAngMomB) then - Gab = 0d0 - else -!------------------------------------------------------------------------ -! Vertical recurrence relations: (a|0) -!------------------------------------------------------------------------ - if(TotAngMomB == 0) then - Gab = VRRNuc(0,AngMomA,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - else -!------------------------------------------------------------------------ -! 1st horizontal recurrence relation (2 terms): (a|b+) -!------------------------------------------------------------------------ - do i=1,3 - ap(i) = AngMomA(i) - bm(i) = AngMomB(i) - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomB(1) > 0) then - xyz = 1 - elseif(AngMomB(2) > 0) then - xyz = 2 - elseif(AngMomB(3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in HRRNuc!' - end if -! End loop over cartesian directions - ap(xyz) = ap(xyz) + 1 - bm(xyz) = bm(xyz) - 1 - Gab = HRRNuc(ap,bm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - + CenterAB(xyz)*HRRNuc(AngMomA,bm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - end if - end if - -end function HRRNuc diff --git a/src/IntPak/HRROv.f90 b/src/IntPak/HRROv.f90 deleted file mode 100644 index cc7311f..0000000 --- a/src/IntPak/HRROv.f90 +++ /dev/null @@ -1,28 +0,0 @@ -recursive function HRROv(AngMomA,AngMomB,ExpPi,CenterAB,CenterPA) & - result(Gab) - -! Horizontal recurrence relations for one-electron overlap integrals - - implicit none - -! Input variables - integer,intent(in) :: AngMomA,AngMomB - double precision,intent(in) :: ExpPi - double precision,intent(in) :: CenterAB,CenterPA - -! Local variables - double precision :: VRROv - double precision :: Gab - - if(AngMomB < 0) then - Gab = 0d0 - else - if(AngMomB == 0) then - Gab = VRROv(AngMomA,ExpPi,CenterPA) - else - Gab = HRROv(AngMomA+1,AngMomB-1,ExpPi,CenterAB,CenterPA) & - + CenterAB*HRROv(AngMomA,AngMomB-1,ExpPi,CenterAB,CenterPA) - end if - end if - -end function HRROv diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90 deleted file mode 100644 index 9310168..0000000 --- a/src/IntPak/IntPak.f90 +++ /dev/null @@ -1,553 +0,0 @@ -program IntPak - - implicit none - include 'parameters.h' - - logical :: debug - logical :: chemist_notation - - logical :: doOv - logical :: doKin - logical :: doNuc - - logical :: doERI - logical :: doF12 - logical :: doYuk - logical :: doErf - - logical :: do3eInt(n3eInt) - logical :: do4eInt(n4eInt) - - integer :: nNuc,nBas,iType - integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin),nS(nspin) - double precision :: ExpS - double precision :: ENuc - integer :: KG - double precision,allocatable :: DG(:),ExpG(:) - double precision,allocatable :: ZNuc(:),rNuc(:,:) - - integer :: nShell - integer,allocatable :: TotAngMomShell(:),KShell(:) - double precision,allocatable :: CenterShell(:,:),DShell(:,:),ExpShell(:,:) - - double precision :: start_1eInt(n1eInt),end_1eInt(n1eInt),t_1eInt(n1eInt) - double precision :: start_2eInt(n2eInt),end_2eInt(n2eInt),t_2eInt(n2eInt) - double precision :: start_3eInt(n3eInt),end_3eInt(n3eInt),t_3eInt(n3eInt) - double precision :: start_4eInt(n4eInt),end_4eInt(n4eInt),t_4eInt(n4eInt) - - integer :: np1eInt(n1eInt),nSigp1eInt(n1eInt),nc1eInt(n1eInt),nSigc1eInt(n1eInt) - integer :: np2eInt(n2eInt),nSigp2eInt(n2eInt),nc2eInt(n2eInt),nSigc2eInt(n2eInt) - integer :: np3eInt(n3eInt),nSigp3eInt(n3eInt),nc3eInt(n3eInt),nSigc3eInt(n3eInt) - integer :: np4eInt(n4eInt),nSigp4eInt(n4eInt),nc4eInt(n4eInt),nSigc4eInt(n4eInt) - - double precision,allocatable :: S(:,:) - - -! Hello World - - write(*,*) - write(*,*) '********************************' - write(*,*) '* IntPak *' - write(*,*) '* Integral Package for dummies *' - write(*,*) '********************************' - write(*,*) - -! Read options for integral calculations - - call read_options(debug,chemist_notation,ExpS,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt) - -!------------------------------------------------------------------------ -! Read input information -!------------------------------------------------------------------------ - -! Read number of atoms, number of electrons of the system -! nO = number of occupied orbitals -! nV = number of virtual orbitals (see below) -! nBas = number of basis functions (see below) -! = nO + nV - - call read_molecule(nNuc,nEl,nO,nC,nR) - - allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3)) - -! Read geometry - - call read_geometry(nNuc,ZNuc,rNuc,ENuc) - - allocate(CenterShell(1:maxShell,1:3),TotAngMomShell(1:maxShell),KShell(1:maxShell), & - DShell(1:maxShell,1:maxK),ExpShell(1:maxShell,1:maxK)) - - call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) - -!------------------------------------------------------------------------ -! Memory allocation -!------------------------------------------------------------------------ - allocate(S(1:nBas,1:nBas)) - -!------------------------------------------------------------------------ -! Compute one-electron overlap integrals -!------------------------------------------------------------------------ - if(doOv) then - - iType = 1 - - call cpu_time(start_1eInt(iType)) - call ComputeOv(debug,nBas,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType),S) - call cpu_time(end_1eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive overlap integrals = ',np1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive overlap integrals = ',nSigp1eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted overlap integrals = ',nc1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted overlap integrals = ',nSigc1eInt(iType) - - write(*,*) - - t_1eInt(iType) = end_1eInt(iType) - start_1eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute one-electron kinetic integrals -!------------------------------------------------------------------------ - - if(doKin) then - - iType = 2 - - call cpu_time(start_1eInt(iType)) - call ComputeKin(debug,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType)) - call cpu_time(end_1eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive kinetic integrals = ',np1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive kinetic integrals = ',nSigp1eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted kinetic integrals = ',nc1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted kinetic integrals = ',nSigc1eInt(iType) - - write(*,*) - - t_1eInt(iType) = end_1eInt(iType) - start_1eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute one-electron nuclear attraction integrals -!------------------------------------------------------------------------ - - if(doNuc) then - - iType = 3 - - call cpu_time(start_1eInt(iType)) - call ComputeNuc(debug,nShell, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - nNuc,ZNuc,rNuc, & - np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType)) - call cpu_time(end_1eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive nuclear integrals = ',np1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive nuclear integrals = ',nSigp1eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted nuclear integrals = ',nc1eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted nuclear integrals = ',nSigc1eInt(iType) - - write(*,*) - - t_1eInt(iType) = end_1eInt(iType) - start_1eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' - write(*,*) - - end if - -!------------------------------------------------------------------------ -! Compute ERIs -!------------------------------------------------------------------------ - - if(doERI) then - - iType = 1 - KG = 1 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 1d0 /) - ExpG = (/ 0d0 /) - - call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,chemist_notation,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) - call cpu_time(end_2eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive ERIs = ',np2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive ERIs = ',nSigp2eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted ERIs = ',nc2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted ERIs = ',nSigc2eInt(iType) - - write(*,*) - - t_2eInt(iType) = end_2eInt(iType) - start_2eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_2eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute F12 two-electron integrals -!------------------------------------------------------------------------ - - if(doF12) then - - iType = 2 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - -! KG = 10 -! allocate(DG(1:KG),ExpG(1:KG)) - -! DG = (/ 220.983854141, 18.52358977132, 4.81060044582, 1.892812227999, & -! 0.920641976732, 0.505281134191, 0.295757471525, 0.1753021140139, & -! 0.0969611396173, 0.0386163391551 /) -! ExpG = (/ 5722.54799330, 191.0413784782, 27.4417708701, 6.39987966572, & -! 1.82203908762, 0.548835646170, 0.156252937904, 0.036440796942, & -! 0.0052344680925, 0.00017474733304 /) - -! KG = 20 -! allocate(DG(1:KG),ExpG(1:KG)) - -! DG = (/ 841.88478132, 70.590185207, 18.3616020768, 7.2608642093, & -!3.57483416444, 2.01376031082, 1.24216542801, 0.81754348620, & -!0.564546514023, 0.404228610699, 0.297458536575, 0.223321219537, & -!0.169933732064, 0.130190978230, 0.099652303426, 0.075428246546, & -!0.0555635614051, 0.0386791283055, 0.0237550435652, 0.0100062783874 /) - -! ExpG = (/84135.654509, 2971.58727634, 474.716025959, 130.676724560, & -!47.3938388887, 20.2078651631, 9.5411021938, 4.8109546955, & -!2.52795733067, 1.35894103210, 0.73586710268, 0.39557629706, & -!0.20785895177, 0.104809693858, 0.049485682527, 0.021099788990, & -!0.007652472186, 0.0021065225215, 0.0003365204879, 0.00001188556749 /) - - - - call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,chemist_notation,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) - call cpu_time(end_2eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive geminal integrals = ',np2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive geminal integrals = ',nSigp2eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted geminal integrals = ',nc2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted geminal integrals = ',nSigc2eInt(iType) - - write(*,*) - - t_2eInt(iType) = end_2eInt(iType) - start_2eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_2eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute Yukawa two-electron integrals -!------------------------------------------------------------------------ - - if(doYuk) then - - iType = 3 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,chemist_notation,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) - call cpu_time(end_2eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive Yukawa integrals = ',np2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive Yukawa integrals = ',nSigp2eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted Yukawa integrals = ',nc2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted Yukawa integrals = ',nSigc2eInt(iType) - - write(*,*) - - t_2eInt(iType) = end_2eInt(iType) - start_2eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_2eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute long-range Coulomb two-electron integrals -!------------------------------------------------------------------------ - - if(doErf) then - - iType = 4 - KG = 1 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 1d0 /) - ExpG = (/ 1d0 /) -! ExpS = ExpS*ExpS - - call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,chemist_notation,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) - call cpu_time(end_2eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive long-range Coulomb integrals = ',np2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive long-range Coulomb integrals = ',nSigp2eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted long-range Coulomb integrals = ',nc2eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted long-range Coulomb integrals = ',nSigc2eInt(iType) - - write(*,*) - - t_2eInt(iType) = end_2eInt(iType) - start_2eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_2eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute three-electron integrals: Type 1 => chain C12 S23 -!------------------------------------------------------------------------ - - if(do3eInt(1)) then - - iType = 1 -! KG = 1 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) -! DG = (/ 1d0 /) -! ExpG = (/ 1d0 /) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_3eInt(iType)) - call Compute3eInt(debug,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np3eInt(iType),nSigp3eInt(iType),nc3eInt(iType),nSigc3eInt(iType)) - call cpu_time(end_3eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f23/r12 integrals = ',np3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f23/r12 integrals = ',nSigp3eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f23/r12 integrals = ',nc3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f23/r12 integrals = ',nSigc3eInt(iType) - - write(*,*) - - t_3eInt(iType) = end_3eInt(iType) - start_3eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_3eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute three-electron integrals: Type 2 => cyclic C12 S13 S23 -!------------------------------------------------------------------------ - - if(do3eInt(2)) then - - iType = 2 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_3eInt(iType)) - call Compute3eInt(debug,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np3eInt(iType),nSigp3eInt(iType),nc3eInt(iType),nSigc3eInt(iType)) - call cpu_time(end_3eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f13.f23/r12 integrals = ',np3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f13.f23/r12 integrals = ',nSigp3eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f13.f23/r12 integrals = ',nc3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f13.f23/r12 integrals = ',nSigc3eInt(iType) - - write(*,*) - - t_3eInt(iType) = end_3eInt(iType) - start_3eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_3eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute three-electron integrals: Type 3 => chain S13 S23 -!------------------------------------------------------------------------ - - if(do3eInt(3)) then - - iType = 3 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_3eInt(iType)) - call Compute3eInt(debug,iType,nShell, & - ExpS,KG,DG,ExpG, & - CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - np3eInt(iType),nSigp3eInt(iType),nc3eInt(iType),nSigc3eInt(iType)) - call cpu_time(end_3eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f13.f23 integrals = ',np3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f13.f23 integrals = ',nSigp3eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f13.f23 integrals = ',nc3eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f13.f23 integrals = ',nSigc3eInt(iType) - - write(*,*) - - t_3eInt(iType) = end_3eInt(iType) - start_3eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_3eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute four-electron integrals: Type 1 => chain C12 S14 S23 -!------------------------------------------------------------------------ - - if(do4eInt(1)) then - - iType = 1 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_4eInt(iType)) -! call Compute4eInt(debug,iType,nShell,ExpS, & -! ExpS,KG,DG,ExpG, & -! CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & -! np4eInt(iType),nSigp4eInt(iType),nc4eInt(iType),nSigc4eInt(iType)) - call cpu_time(end_4eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f14.f23/r12 integrals = ',np4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f14.f23/r12 integrals = ',nSigp4eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f14.f23/r12 integrals = ',nc4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f14.f23/r12 integrals = ',nSigc4eInt(iType) - - write(*,*) - - t_4eInt(iType) = end_4eInt(iType) - start_4eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_4eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute four-electron integrals: Type 2 => trident C12 S13 S14 -!------------------------------------------------------------------------ - - if(do4eInt(2)) then - - iType = 2 - KG = 6 - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_4eInt(iType)) -! call Compute4eInt(debug,iType,nShell,ExpS, & -! ExpS,KG,DG,ExpG, & -! CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & -! np4eInt(iType),nSigp4eInt(iType),nc4eInt(iType),nSigc4eInt(iType)) - call cpu_time(end_4eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f13.f14/r12 integrals = ',np4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f13.f14/r12 integrals = ',nSigp4eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f13.f14/r12 integrals = ',nc4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f13.f14/r12 integrals = ',nSigc4eInt(iType) - - write(*,*) - - t_4eInt(iType) = end_4eInt(iType) - start_4eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_4eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if - -!------------------------------------------------------------------------ -! Compute four-electron integrals: Type 3 => chain C12 S13 S34 -!------------------------------------------------------------------------ - - if(do4eInt(3)) then - - iType = 3 - KG = 6 - allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) - ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) - - call cpu_time(start_4eInt(iType)) -! call Compute4eInt(debug,iType,nShell, & -! ExpS,KG,DG,ExpG, & -! CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & -! np4eInt(iType),nSigp4eInt(iType),nc4eInt(iType),nSigc4eInt(iType)) - call cpu_time(end_4eInt(iType)) - - write(*,'(A65,1X,I9)') 'Total number of primitive f13.f34/r12 integrals = ',np4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant primitive f13.f34/r12 integrals = ',nSigp4eInt(iType) - - write(*,'(A65,1X,I9)') 'Total number of contracted f13.f34/r12 integrals = ',nc4eInt(iType) - write(*,'(A65,1X,I9)') 'Number of significant contracted f13.f34/r12 integrals = ',nSigc4eInt(iType) - - write(*,*) - - t_4eInt(iType) = end_4eInt(iType) - start_4eInt(iType) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_4eInt(iType),' seconds' - write(*,*) - - deallocate(DG,ExpG) - - end if -!------------------------------------------------------------------------ -! End of IntPak -!------------------------------------------------------------------------ -end program IntPak diff --git a/src/IntPak/KinInt.f90 b/src/IntPak/KinInt.f90 deleted file mode 100644 index c7f892d..0000000 --- a/src/IntPak/KinInt.f90 +++ /dev/null @@ -1,76 +0,0 @@ -subroutine KinInt(npKin,nSigpKin,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,pKin) - -! Compute one-electron kinetic integrals - - implicit none - -! Input variables - - double precision,intent(in) :: ExpA,ExpB - double precision,intent(in) :: CenterA(3),CenterB(3) - integer,intent(in) :: AngMomA(3),AngMomB(3) - - -! Local variables - - double precision :: ExpAi,ExpBi - double precision :: ExpP,ExpPi - double precision :: CenterP(3),CenterAB(3),CenterPA(3) - double precision :: NormABSq - double precision :: GAB - double precision :: HRROv,RRKin - - integer :: i - double precision :: pi - double precision :: start_RR,finish_RR,t_RR - double precision :: s(3),k(3) - -! Output variables - - integer,intent(inout) :: npKin,nSigpKin - double precision,intent(out) :: pKin - - pi = 4d0*atan(1d0) - -! Pre-computed shell quantities - - ExpAi = 1d0/ExpA - ExpBi = 1d0/ExpB - -! Pre-computed quantities for shell-pair AB - - ExpP = ExpA + ExpB - ExpPi = 1d0/ExpP - - NormABSq = 0d0 - Do i=1,3 - CenterP(i) = (ExpA*CenterA(i) + ExpB*CenterB(i))*ExpPi - CenterPA(i) = CenterP(i) - CenterA(i) - CenterAB(i) = CenterA(i) - CenterB(i) - NormABSq = NormABSq + CenterAB(i)**2 - Enddo - - GAB = (pi*ExpPi)**(1.5d0)*exp(-NormABSq/(ExpAi+ExpBi)) - -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) -! Loop over cartesian directions - Do i=1,3 - s(i) = HRROv(AngMomA(i),AngMomB(i),ExpPi,CenterAB(i),CenterPA(i)) - k(i) = RRKin(AngMomA(i),AngMomB(i),ExpA,ExpB,ExpPi,CenterAB(i),CenterPA(i)) - Enddo - call cpu_time(finish_RR) - - pKin = k(1)*s(2)*s(3) + s(1)*k(2)*s(3) + s(1)*s(2)*k(3) - pKin = GAB*pKin - t_RR = finish_RR - start_RR - -! Print result - npKin = npKin + 1 - if(abs(pKin) > 1d-15) then - nSigpKin = nSigpKin + 1 - end if - -end subroutine KinInt diff --git a/src/IntPak/Makefile b/src/IntPak/Makefile deleted file mode 100644 index ee03122..0000000 --- a/src/IntPak/Makefile +++ /dev/null @@ -1,35 +0,0 @@ -IDIR =../../include -BDIR =../../bin -ODIR = obj -SDIR =. -FC = gfortran -I$(IDIR) -ifeq ($(DEBUG),1) -FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -else -FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3 -endif - -LIBS = ../../lib/*.a -LIBS += -lblas -llapack - -SRCF90 = $(wildcard *.f90) - -SRC = $(wildcard *.f) - -OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC)) - -$(ODIR)/%.o: %.f90 - $(FC) -c -o $@ $< $(FFLAGS) - -$(ODIR)/%.o: %.f - $(FC) -c -o $@ $< $(FFLAGS) - -$(BDIR)/IntPak: $(OBJ) - $(FC) -o $@ $^ $(FFLAGS) $(LIBS) - -debug: - DEBUG=1 make $(BDIR)/IntPak -#DEBUG=1 make clean $(BDIR)/IntPak - -clean: - rm -f $(ODIR)/*.o $(BDIR)/IntPak $(BDIR)/debug diff --git a/src/IntPak/NucInt.f90 b/src/IntPak/NucInt.f90 deleted file mode 100644 index 6e18276..0000000 --- a/src/IntPak/NucInt.f90 +++ /dev/null @@ -1,114 +0,0 @@ -subroutine NucInt(debug,npNuc,nSigpNuc, & - ExpA,CenterA,AngMomA, & - ExpB,CenterB,AngMomB, & - CenterC, & - pNuc) - -! Compute recursively the primitive one-electron nuclear attraction integrals - - implicit none - -! Input variables - - logical,intent(in) :: debug - double precision,intent(in) :: ExpA,ExpB - double precision,intent(in) :: CenterA(3),CenterB(3),CenterC(3) - integer,intent(in) :: AngMomA(3),AngMomB(3) - -! Local variables - - double precision :: ExpAi,ExpBi - integer :: TotAngMomA,TotAngMomB - double precision :: ExpP,ExpPi - double precision :: CenterP(3),CenterAB(3),CenterPA(3),CenterPC(3) - double precision :: NormABSq,NormPCSq - double precision :: G - double precision,allocatable :: Om(:) - double precision :: HRRNuc - double precision :: Gab - - double precision :: pi - integer :: i,maxm - double precision :: start_Om,finish_Om,start_RR,finish_RR,t_Om,t_RR - -! Output variables - - integer,intent(inout) :: npNuc,nSigpNuc - double precision,intent(out) :: pNuc - - pi = 4d0*atan(1d0) - -! Pre-computed shell quantities - - ExpAi = 1d0/ExpA - ExpBi = 1d0/ExpB - -! Pre-computed quantities for shell-pair AB - - ExpP = ExpA + ExpB - ExpPi = 1d0/ExpP - - NormABSq = 0d0 - NormPCSq = 0d0 - do i=1,3 - CenterP(i) = (ExpA*CenterA(i) + ExpB*CenterB(i))*ExpPi - CenterAB(i) = CenterA(i) - CenterB(i) - CenterPA(i) = CenterP(i) - CenterA(i) - CenterPC(i) = CenterP(i) - CenterC(i) - NormABSq = NormABSq + CenterAB(i)**2 - NormPCSq = NormPCSq + CenterPC(i)**2 - end do - - G = (pi*ExpPi)**(1.5d0)*exp(-NormABSq/(ExpAi+ExpBi)) - -! Total angular momemtum - - TotAngMomA = AngMomA(1) + AngMomA(2) + AngMomA(3) - TotAngMomB = AngMomB(1) + AngMomB(2) + AngMomB(3) - - maxm = TotAngMomA + TotAngMomB - -! Pre-compute (0|V|0)^m - - allocate(Om(0:maxm)) - call cpu_time(start_Om) - call CalcOmNuc(maxm,ExpPi,NormPCSq,Om) - call cpu_time(finish_Om) - -! Print (0|V|0)^m - - if(debug) then - write(*,*) '(0|V|0)^m' - do i=0,maxm - write(*,*) i,Om(i) - end do - write(*,*) - end if - -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) - Gab = HRRNuc(AngMomA,AngMomB,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - call cpu_time(finish_RR) - -! Timings - - t_Om = finish_Om - start_Om - t_RR = finish_RR - start_RR - -! Print result - - pNuc = G*Gab - - npNuc = npNuc + 1 - if(abs(pNuc) > 1d-15) then - nSigpNuc = nSigpNuc + 1 -! write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '[a|V|b] = ',pNuc - end if - -! Deallocate arrays - - deallocate(Om) - -end subroutine NucInt diff --git a/src/IntPak/OvInt.f90 b/src/IntPak/OvInt.f90 deleted file mode 100644 index c8f930c..0000000 --- a/src/IntPak/OvInt.f90 +++ /dev/null @@ -1,74 +0,0 @@ -subroutine OvInt(npOv,nSigpOv,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,pOv) - -! Compute one-electron overlap integrals - - implicit none - -! Input variables - - double precision,intent(in) :: ExpA,ExpB - double precision,intent(in) :: CenterA(3),CenterB(3) - integer,intent(in) :: AngMomA(3),AngMomB(3) - - -! Local variables - - double precision :: ExpAi,ExpBi - double precision :: ExpP,ExpPi - double precision :: CenterP(3),CenterAB(3),CenterPA(3) - double precision :: NormABSq - double precision :: G - double precision :: HRROv - - integer :: i - double precision :: pi - double precision :: start_RR,finish_RR,t_RR - double precision :: Gab(3) - -! Output variables - - integer,intent(inout) :: npOv,nSigpOv - double precision,intent(out) :: pOv - - pi = 4d0*atan(1d0) - -! Pre-computed shell quantities - - ExpAi = 1d0/ExpA - ExpBi = 1d0/ExpB - -! Pre-computed quantities for shell-pair AB - - ExpP = ExpA + ExpB - ExpPi = 1d0/ExpP - - NormABSq = 0d0 - Do i=1,3 - CenterP(i) = (ExpA*CenterA(i) + ExpB*CenterB(i))*ExpPi - CenterPA(i) = CenterP(i) - CenterA(i) - CenterAB(i) = CenterA(i) - CenterB(i) - NormABSq = NormABSq + CenterAB(i)**2 - Enddo - - G = (pi*ExpPi)**(1.5d0)*exp(-NormABSq/(ExpAi+ExpBi)) - -!------------------------------------------------------------------------ -! Launch reccurence relations! -!------------------------------------------------------------------------ - call cpu_time(start_RR) -! Loop over cartesian directions - Do i=1,3 - Gab(i) = HRROv(AngMomA(i),AngMomB(i),ExpPi,CenterAB(i),CenterPA(i)) - Enddo - call cpu_time(finish_RR) - - pOv = G*Gab(1)*Gab(2)*Gab(3) - t_RR = finish_RR - start_RR - -! Print result - npOv = npOv + 1 - if(abs(pOv) > 1d-15) then - nSigpOv = nSigpOv + 1 - end if - -end subroutine OvInt diff --git a/src/IntPak/RRKin.f90 b/src/IntPak/RRKin.f90 deleted file mode 100644 index 30993f2..0000000 --- a/src/IntPak/RRKin.f90 +++ /dev/null @@ -1,29 +0,0 @@ -function RRKin(AngMomA,AngMomB,ExpA,ExpB,ExpPi,CenterAB,CenterPA) & - result(Gab) - -! Recurrence relation for one-electron kinetic integrals - - implicit none - -! Input variables - integer,intent(in) :: AngMomA,AngMomB - double precision,intent(in) :: ExpA,ExpB,ExpPi - double precision,intent(in) :: CenterAB,CenterPA - -! Local variables - double precision :: HRROv - double precision :: a,b,s1,s2,s3,s4 - double precision :: Gab - - a = dble(AngMomA) - b = dble(AngMomB) - - s1 = HRROv(AngMomA-1,AngMomB-1,ExpPi,CenterAB,CenterPA) - s2 = HRROv(AngMomA+1,AngMomB-1,ExpPi,CenterAB,CenterPA) - s3 = HRROv(AngMomA-1,AngMomB+1,ExpPi,CenterAB,CenterPA) - s4 = HRROv(AngMomA+1,AngMomB+1,ExpPi,CenterAB,CenterPA) - - Gab = 0.5d0*a*b*s1 - ExpA*b*s2 - a*ExpB*s3 + 2d0*ExpA*ExpB*s4 - - -end function RRKin diff --git a/src/IntPak/S2eInt.f90 b/src/IntPak/S2eInt.f90 deleted file mode 100644 index e3475d0..0000000 --- a/src/IntPak/S2eInt.f90 +++ /dev/null @@ -1,70 +0,0 @@ -subroutine S2eInt(debug,iType,np2eInt,nSigp2eInt, & - ExpS,KG,DG,ExpG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet, & - p2eInt) - -! Perform contraction over the operator for two-electron integrals - - implicit none - include 'parameters.h' - - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: iType - double precision,intent(in) :: ExpS - integer,intent(in) :: KG - double precision,intent(in) :: DG(KG),ExpG(KG) - double precision,intent(in) :: ExpBra(2),ExpKet(2) - double precision,intent(in) :: CenterBra(2,3),CenterKet(2,3) - integer,intent(in) :: AngMomBra(2,3),AngMomKet(2,3) - -! Local variables - - double precision :: ExpSG - double precision :: G2eInt,GF12Int - - integer :: k - -! Output variables - - integer,intent(out) :: np2eInt,nSigp2eInt - double precision :: p2eInt - - p2eInt = 0d0 - -! Gaussian geminal - - if(iType == 2) then - do k=1,KG - ExpSG = ExpG(k)*ExpS**2 - p2eInt = p2eInt & - + DG(k)*GF12Int(ExpSG, & - ExpBra(1),CenterBra(1,1:3),AngMomBra(1,1:3), & - ExpKet(1),CenterKet(1,1:3),AngMomKet(1,1:3), & - ExpBra(2),CenterBra(2,1:3),AngMomBra(2,1:3), & - ExpKet(2),CenterKet(2,1:3),AngMomKet(2,1:3)) - end do - else - do k=1,KG - ExpSG = ExpG(k)*ExpS**2 - p2eInt = p2eInt & - + DG(k)*G2eInt(debug,iType, & - ExpSG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet) - end do - end if - -! Print result - - np2eInt = np2eInt + 1 - - if(abs(p2eInt) > 1d-15) then - nSigp2eInt = nSigp2eInt + 1 - if(.false.) write(*,'(A15,1X,F16.10)') '[a1a2|b1b2] = ',p2eInt - end if - -end subroutine S2eInt diff --git a/src/IntPak/S3eInt.f90 b/src/IntPak/S3eInt.f90 deleted file mode 100644 index 997862c..0000000 --- a/src/IntPak/S3eInt.f90 +++ /dev/null @@ -1,76 +0,0 @@ -subroutine S3eInt(debug,iType,np3eInt,nSigp3eInt, & - ExpS,KG,DG,ExpG, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet, & - p3eInt) - -! Perform contraction over the operators for three-electron integrals - - implicit none - include 'parameters.h' - - -! Input variables - - logical,intent(in) :: debug - integer,intent(in) :: iType - double precision,intent(in) :: ExpS - integer,intent(in) :: KG - double precision,intent(in) :: DG(KG),ExpG(KG) - double precision,intent(in) :: ExpBra(3),ExpKet(3) - double precision,intent(in) :: CenterBra(3,3),CenterKet(3,3) - integer,intent(in) :: AngMomBra(3,3),AngMomKet(3,3) - -! Local variables - - double precision :: ExpSG13,ExpSG23 - double precision :: G3eInt - - integer :: k,l - -! Output variables - - integer,intent(out) :: np3eInt,nSigp3eInt - double precision :: p3eInt - - p3eInt = 0d0 - - if(iType == 1) then - - do k=1,KG - ExpSG13 = ExpG(k)*ExpS**2 - p3eInt = p3eInt & - + DG(k)*G3eInt(debug,iType, & - ExpSG13,ExpSG23, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet) - end do - - end if - - if(iType == 2 .or. iType == 3) then - - do k=1,KG - do l=1,KG - ExpSG13 = ExpG(k)*ExpS**2 - ExpSG23 = ExpG(l)*ExpS**2 - p3eInt = p3eInt & - + DG(k)*DG(l)*G3eInt(debug,iType, & - ExpSG13,ExpSG23, & - ExpBra,CenterBra,AngMomBra, & - ExpKet,CenterKet,AngMomKet) - end do - end do - - end if - -! Print result - - np3eInt = np3eInt + 1 - - if(abs(p3eInt) > 1d-15) then - nSigp3eInt = nSigp3eInt + 1 - if(.false.) write(*,'(A15,1X,F16.10)') '[a1a2a3|b1b2b3] = ',p3eInt - end if - -end subroutine S3eInt diff --git a/src/IntPak/VRR2e.f90 b/src/IntPak/VRR2e.f90 deleted file mode 100644 index 1b2189d..0000000 --- a/src/IntPak/VRR2e.f90 +++ /dev/null @@ -1,130 +0,0 @@ -recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - result(a1a2) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: m - integer,intent(in) :: AngMomBra(2,3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm),ExpZi(2),ExpY(2,2) - double precision,intent(in) :: CenterZA(2,3),CenterY(2,2,3) - -! Local variables - - logical :: NegAngMomBra(2) - integer :: TotAngMomBra(2) - integer :: a1m(2,3),a2m(2,3) - integer :: a1mm(2,3),a2mm(2,3) - integer :: a1m2m(2,3) - double precision :: fZ(2) - integer :: i,j,xyz - -! Output variables - - double precision :: a1a2 - - do i=1,2 - NegAngMomBra(i) = AngMomBra(i,1) < 0 .or. AngMomBra(i,2) < 0 .or. AngMomBra(i,3) < 0 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - end do - - fZ(1) = ExpY(1,2)*ExpZi(1) - fZ(2) = ExpY(1,2)*ExpZi(2) - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ -! if(NegAngMomBra(1) .or. NegAngMomBra(2)) then -! a1a2 = 0d0 -!------------------------------------------------------------------------ -! Fundamental integral: (00|00)^m -!------------------------------------------------------------------------ -! elseif(TotAngMomBra(1) == 0 .and. TotAngMomBra(2) == 0) then - if(TotAngMomBra(1) == 0 .and. TotAngMomBra(2) == 0) then - a1a2 = Om(m) -!------------------------------------------------------------------------ -! 1st vertical recurrence relation (4 terms): (a+0|00)^m -!------------------------------------------------------------------------ - elseif(TotAngMomBra(2) == 0) then - do i=1,2 - do j=1,3 - a1m(i,j) = AngMomBra(i,j) - a1mm(i,j) = AngMomBra(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomBra(1,1) > 0) then - xyz = 1 - elseif(AngMomBra(1,2) > 0) then - xyz = 2 - elseif(AngMomBra(1,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRR2e!' - end if -! End loop over cartesian directions - a1m(1,xyz) = a1m(1,xyz) - 1 - a1mm(1,xyz) = a1mm(1,xyz) - 2 - if(AngMomBra(1,xyz) <= 0) then - a1a2 = 0d0 - elseif(AngMomBra(1,xyz) == 1) then - a1a2 = CenterZA(1,xyz)*VRR2e(m,a1m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - - fZ(1)*CenterY(1,2,xyz)*VRR2e(m+1,a1m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - else - a1a2 = CenterZA(1,xyz)*VRR2e(m,a1m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - - fZ(1)*CenterY(1,2,xyz)*VRR2e(m+1,a1m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - + 0.5d0*dble(AngMomBra(1,xyz)-1)*ExpZi(1)*( & - VRR2e(m,a1mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - - fZ(1)*VRR2e(m+1,a1mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY)) - end if -!------------------------------------------------------------------------ -! 2nd vertical recurrence relation (5 terms): (a0|c+0)^m -!------------------------------------------------------------------------ - else - do i=1,2 - do j=1,3 - a2m(i,j) = AngMomBra(i,j) - a2mm(i,j) = AngMomBra(i,j) - a1m2m(i,j) = AngMomBra(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomBra(2,1) > 0) then - xyz = 1 - elseif(AngMomBra(2,2) > 0) then - xyz = 2 - elseif(AngMomBra(2,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRR2e!' - end if -! End loop over cartesian directions - a2m(2,xyz) = a2m(2,xyz) - 1 - a2mm(2,xyz) = a2mm(2,xyz) - 2 - a1m2m(1,xyz) = a1m2m(1,xyz) - 1 - a1m2m(2,xyz) = a1m2m(2,xyz) - 1 - if(AngMomBra(2,xyz) <= 0) then - a1a2 = 0d0 - elseif(AngMomBra(2,xyz) == 1) then - a1a2 = CenterZA(2,xyz)*VRR2e(m,a2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - + fZ(2)*CenterY(1,2,xyz)*VRR2e(m+1,a2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - else - a1a2 = CenterZA(2,xyz)*VRR2e(m,a2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - + fZ(2)*CenterY(1,2,xyz)*VRR2e(m+1,a2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - + 0.5d0*dble(AngMomBra(2,xyz)-1)*ExpZi(2)*( & - VRR2e(m,a2mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - - fZ(2)*VRR2e(m+1,a2mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY)) - end if - if(AngMomBra(1,xyz) > 0) & - a1a2 = a1a2 & - + 0.5d0*dble(AngMomBra(1,xyz))*fZ(2)*ExpZi(1)*VRR2e(m+1,a1m2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - end if - -end function VRR2e diff --git a/src/IntPak/VRR3e.f90 b/src/IntPak/VRR3e.f90 deleted file mode 100644 index a866e97..0000000 --- a/src/IntPak/VRR3e.f90 +++ /dev/null @@ -1,174 +0,0 @@ -recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - result(a1a2a3) - -! Vertical recurrence relations for three-electron integrals - - implicit none - include 'parameters.h' - - -! Input variables - - integer,intent(in) :: m - integer,intent(in) :: AngMomBra(3,3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm),ExpZ(3),CenterZA(3,3) - double precision,intent(in) :: DY0(3),DY1(3),D2Y0(3,3),D2Y1(3,3) - -! Local variables - - logical :: NegAngMomBra(3) - integer :: TotAngMomBra(3) - integer :: a1m(3,3),a2m(3,3),a3m(3,3) - integer :: a1mm(3,3),a2mm(3,3),a3mm(3,3) - integer :: a1m2m(3,3),a1m3m(3,3),a2m3m(3,3) - integer :: i,j,xyz - -! Output variables - - double precision :: a1a2a3 - - do i=1,3 - NegAngMomBra(i) = AngMomBra(i,1) < 0 .or. AngMomBra(i,2) < 0 .or. AngMomBra(i,3) < 0 - TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - end do - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ - if(NegAngMomBra(1) .or. NegAngMomBra(2) .or. NegAngMomBra(3)) then - a1a2a3 = 0d0 -!------------------------------------------------------------------------ -! Fundamental integral: (000|000)^m -!------------------------------------------------------------------------ - elseif(TotAngMomBra(1) == 0 .and. TotAngMomBra(2) == 0 .and. TotAngMomBra(3) == 0) then - a1a2a3 = Om(m) -!------------------------------------------------------------------------ -! 1st vertical recurrence relation (4 terms): (a1+00|000)^m -!------------------------------------------------------------------------ - elseif(TotAngMomBra(2) == 0 .and. TotAngMomBra(3) == 0) then - do i=1,3 - do j=1,3 - a1m(i,j) = AngMomBra(i,j) - a1mm(i,j) = AngMomBra(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomBra(1,1) > 0) then - xyz = 1 - elseif(AngMomBra(1,2) > 0) then - xyz = 2 - elseif(AngMomBra(1,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRR3e!' - end if -! End loop over cartesian directions - a1m(1,xyz) = a1m(1,xyz) - 1 - a1mm(1,xyz) = a1mm(1,xyz) - 2 - if(AngMomBra(1,xyz) == 0) then - a1a2a3 = 0d0 - elseif(AngMomBra(1,xyz) == 1) then - a1a2a3 = (CenterZA(1,xyz) - DY0(1))*VRR3e(m, a1m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(1) - DY0(1))*VRR3e(m+1,a1m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - else - a1a2a3 = (CenterZA(1,xyz) - DY0(1))*VRR3e(m, a1m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(1) - DY0(1))*VRR3e(m+1,a1m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + dble(AngMomBra(1,xyz)-1)*(0.5d0/ExpZ(1) - D2Y0(1,1))*VRR3e(m, a1mm,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(1,xyz)-1)*(D2Y1(1,1) - D2Y0(1,1))*VRR3e(m+1,a1mm,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if -!------------------------------------------------------------------------ -! 2nd vertical recurrence relation (6 terms): (a1a2+0|000)^m -!------------------------------------------------------------------------ - elseif(TotAngMomBra(3) == 0) then - do i=1,3 - do j=1,3 - a2m(i,j) = AngMomBra(i,j) - a2mm(i,j) = AngMomBra(i,j) - a1m2m(i,j) = AngMomBra(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomBra(2,1) > 0) then - xyz = 1 - elseif(AngMomBra(2,2) > 0) then - xyz = 2 - elseif(AngMomBra(2,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRR3e!' - end if -! End loop over cartesian directions - a2m(2,xyz) = a2m(2,xyz) - 1 - a2mm(2,xyz) = a2mm(2,xyz) - 2 - a1m2m(1,xyz) = a1m2m(1,xyz) - 1 - a1m2m(2,xyz) = a1m2m(2,xyz) - 1 - if(AngMomBra(2,xyz) == 0) then - a1a2a3 = 0d0 - elseif(AngMomBra(2,xyz) == 1) then - a1a2a3 = (CenterZA(2,xyz) - DY0(2))*VRR3e(m, a2m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(2) - DY0(2))*VRR3e(m+1,a2m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - else - a1a2a3 = (CenterZA(2,xyz) - DY0(2))*VRR3e(m, a2m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(2) - DY0(2))*VRR3e(m+1,a2m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + dble(AngMomBra(2,xyz)-1)*(0.5d0/ExpZ(2) - D2Y0(2,2))*VRR3e(m, a2mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(2,xyz)-1)*(D2Y1(2,2) - D2Y0(2,2))*VRR3e(m+1,a2mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - if(AngMomBra(1,xyz) > 0) & - a1a2a3 = a1a2a3 & - + dble(AngMomBra(1,xyz))*(-D2Y0(2,1))*VRR3e(m, a1m2m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(1,xyz))*(D2Y1(2,1) - D2Y0(2,1))*VRR3e(m+1,a1m2m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) -!------------------------------------------------------------------------ -! 3rd vertical recurrence relation (8 terms): (a1a2a3+|000)^m -!------------------------------------------------------------------------ - else - do i=1,3 - do j=1,3 - a3m(i,j) = AngMomBra(i,j) - a3mm(i,j) = AngMomBra(i,j) - a1m3m(i,j) = AngMomBra(i,j) - a2m3m(i,j) = AngMomBra(i,j) - end do - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomBra(3,1) > 0) then - xyz = 1 - elseif(AngMomBra(3,2) > 0) then - xyz = 2 - elseif(AngMomBra(3,3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRR3e!' - end if -! End loop over cartesian directions - a3m(3,xyz) = a3m(3,xyz) - 1 - a3mm(3,xyz) = a3mm(3,xyz) - 2 - a1m3m(1,xyz) = a1m3m(1,xyz) - 1 - a1m3m(3,xyz) = a1m3m(3,xyz) - 1 - a2m3m(2,xyz) = a2m3m(2,xyz) - 1 - a2m3m(3,xyz) = a2m3m(3,xyz) - 1 - if(AngMomBra(3,xyz) == 0) then - a1a2a3 = 0d0 - elseif(AngMomBra(3,xyz) == 1) then - a1a2a3 = (CenterZA(3,xyz) - DY0(3))*VRR3e(m, a3m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(3) - DY0(3))*VRR3e(m+1,a3m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - else - a1a2a3 = (CenterZA(3,xyz) - DY0(3))*VRR3e(m, a3m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - (DY1(3) - DY0(3))*VRR3e(m+1,a3m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - + dble(AngMomBra(3,xyz)-1)*(0.5d0/ExpZ(3) - D2Y0(3,3))*VRR3e(m, a3mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(3,xyz)-1)*(D2Y1(3,3) - D2Y0(3,3))*VRR3e(m+1,a3mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - if(dble(AngMomBra(1,xyz)) > 0) & - a1a2a3 = a1a2a3 & - + dble(AngMomBra(1,xyz))*(-D2Y0(3,1))*VRR3e(m, a1m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(1,xyz))*(D2Y1(3,1) - D2Y0(3,1))*VRR3e(m+1,a1m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - if(dble(AngMomBra(2,xyz)) > 0) & - a1a2a3 = a1a2a3 & - + dble(AngMomBra(2,xyz))*(-D2Y0(3,2))*VRR3e(m, a2m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - - dble(AngMomBra(2,xyz))*(D2Y1(3,2) - D2Y0(3,2))*VRR3e(m+1,a2m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - end if - -end function VRR3e diff --git a/src/IntPak/VRRF12.f90 b/src/IntPak/VRRF12.f90 deleted file mode 100644 index 0537cf9..0000000 --- a/src/IntPak/VRRF12.f90 +++ /dev/null @@ -1,36 +0,0 @@ -recursive function VRRF12(AngMomA,AngMomC,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) & - result(Gac) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - -! Input variables - - integer,intent(in) :: AngMomA,AngMomC - double precision,intent(in) :: ExpPGQi - double precision,intent(in) :: fG,gP,gG,gQ - double precision,intent(in) :: CenterPQSq,CenterRA,CenterRC - -! Output variables - - double precision :: Gac - - if(AngMomA < 0 .or. AngMomC < 0) then - Gac = 0d0 - else - if(AngMomA == 0 .and. AngMomC == 0) then - Gac = sqrt(fG)*exp(-CenterPQSq/ExpPGQi) - else - If(AngMomC == 0) then - Gac = CenterRA*VRRF12(AngMomA-1,AngMomC,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) & - + dble(AngMomA-1)*gP*VRRF12(AngMomA-2,AngMomC,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) - else - Gac = CenterRC*VRRF12(AngMomA,AngMomC-1,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) & - + dble(AngMomA)*gG*VRRF12(AngMomA-1,AngMomC-1,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) & - + dble(AngMomC-1)*gQ*VRRF12(AngMomA,AngMomC-2,fG,gP,gG,gQ,ExpPGQi,CenterPQSq,CenterRA,CenterRC) - endIf - endIf - endIf - -end function VRRF12 diff --git a/src/IntPak/VRRNuc.f90 b/src/IntPak/VRRNuc.f90 deleted file mode 100644 index 026e7a2..0000000 --- a/src/IntPak/VRRNuc.f90 +++ /dev/null @@ -1,76 +0,0 @@ -recursive function VRRNuc(m,AngMomA,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - result(Ga) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - -! Input variables - - integer,intent(in) :: m - integer,intent(in) :: AngMomA(3) - integer,intent(in) :: maxm - double precision,intent(in) :: Om(0:maxm) - double precision,intent(in) :: ExpPi - double precision,intent(in) :: CenterAB(3),CenterPA(3),CenterPC(3) - -! Local variables - - logical :: NegAngMomA - integer :: TotAngMomA - integer :: xyz,am(3),amm(3) - integer :: i - -! Output variables - - double precision :: Ga - - NegAngMomA = AngMomA(1) < 0 .or. AngMomA(2) < 0 .or. AngMomA(3) < 0 - TotAngMomA = AngMomA(1) + AngMomA(2) + AngMomA(3) - -!------------------------------------------------------------------------ -! Termination condition -!------------------------------------------------------------------------ - - if(NegAngMomA) then - - Ga = 0d0 - - else -!------------------------------------------------------------------------ -! Fundamental integral: (0|0)^m -!------------------------------------------------------------------------ - if(TotAngMomA == 0) then - - Ga = Om(m) - - else -!------------------------------------------------------------------------ -! Vertical recurrence relation (4 terms): (a+|0)^m -!------------------------------------------------------------------------ - do i=1,3 - am(i) = AngMomA(i) - amm(i) = AngMomA(i) - end do -! Loop over cartesian directions - xyz = 0 - if (AngMomA(1) > 0) then - xyz = 1 - elseif(AngMomA(2) > 0) then - xyz = 2 - elseif(AngMomA(3) > 0) then - xyz = 3 - else - write(*,*) 'xyz = 0 in VRRNuc!' - end if -! End loop over cartesian directions - am(xyz) = am(xyz) - 1 - amm(xyz) = amm(xyz) - 2 - Ga = CenterPA(xyz)*VRRNuc(m,am,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - + 0.5d0*dble(am(xyz))*ExpPi*VRRNuc(m,amm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - - CenterPC(xyz)*ExpPi*VRRNuc(m+1,am,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - - 0.5d0*dble(am(xyz))*ExpPi**2*VRRNuc(m+1,amm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - end if - end if - -end function VRRNuc diff --git a/src/IntPak/VRROv.f90 b/src/IntPak/VRROv.f90 deleted file mode 100644 index 10627ca..0000000 --- a/src/IntPak/VRROv.f90 +++ /dev/null @@ -1,28 +0,0 @@ -recursive function VRROv(AngMomA,ExpPi,CenterPA) & - result(Ga) - -! Compute two-electron integrals over Gaussian geminals - - implicit none - -! Input variables - - integer,intent(in) :: AngMomA - double precision,intent(in) :: ExpPi - double precision,intent(in) :: CenterPA - -! Output variables - - double precision :: Ga - - if(AngMomA < 0) then - Ga = 0d0 - else - if(AngMomA == 0) then - Ga = 1d0 - else - Ga = CenterPA*VRROv(AngMomA-1,ExpPi,CenterPA) + 0.5d0*dble(AngMomA-1)*ExpPi*VRROv(AngMomA-2,ExpPi,CenterPA) - end if - end if - -end function VRROv diff --git a/src/IntPak/obj/.gitignore b/src/IntPak/obj/.gitignore deleted file mode 100644 index 5761abc..0000000 --- a/src/IntPak/obj/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.o diff --git a/src/IntPak/read_options.f90 b/src/IntPak/read_options.f90 deleted file mode 100644 index d4bbfec..0000000 --- a/src/IntPak/read_options.f90 +++ /dev/null @@ -1,122 +0,0 @@ -subroutine read_options(debug,chemist_notation,ExpS,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt) - -! Read desired methods - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(out) :: debug - logical,intent(out) :: chemist_notation - - double precision,intent(out) :: ExpS - - logical,intent(out) :: doOv - logical,intent(out) :: doKin - logical,intent(out) :: doNuc - - logical,intent(out) :: doERI - logical,intent(out) :: doF12 - logical,intent(out) :: doYuk - logical,intent(out) :: doErf - - logical,intent(out) :: do3eInt(n3eInt) - logical,intent(out) :: do4eInt(n4eInt) - -! Local variables - - character(len=1) :: answer1,answer2,answer3,answer4 - -! Open file with method specification - - open(unit=1,file='input/int') - -! Read HF options - - debug = .false. - chemist_notation = .false. - ExpS = 1d0 - - doOv = .false. - doKin = .false. - doNuc = .false. - - doERI = .false. - doF12 = .false. - doYuk = .false. - doErf = .false. - - do3eInt(1) = .false. - do3eInt(2) = .false. - do3eInt(3) = .false. - - do4eInt(1) = .false. - do4eInt(2) = .false. - do4eInt(3) = .false. - -! Debugging mode - - read(1,*) - read(1,*) answer1 - - if(answer1 == 'T') debug = .true. - -! Chem or Phys notations? - - read(1,*) - read(1,*) answer1 - - if(answer1 == 'T') chemist_notation = .true. - -! Read exponent of Slater geminal - read(1,*) - read(1,*) ExpS - - write(*,'(A28)') '------------------' - write(*,'(A28,1X,F16.10)') 'Slater geminal exponent',ExpS - write(*,'(A28)') '------------------' - write(*,*) - -! One-electron integrals - - read(1,*) - read(1,*) answer1,answer2,answer3 - - if(answer1 == 'T') doOv = .true. - if(answer2 == 'T') doKin = .true. - if(answer3 == 'T') doNuc = .true. - -! Two-electron integrals - - read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - - if(answer1 == 'T') doERI = .true. - if(answer2 == 'T') doF12 = .true. - if(answer3 == 'T') doYuk = .true. - if(answer4 == 'T') doErf = .true. - -! Three-electron integrals - - read(1,*) - read(1,*) answer1,answer2,answer3 - - if(answer1 == 'T') do3eInt(1) = .true. - if(answer2 == 'T') do3eInt(2) = .true. - if(answer3 == 'T') do3eInt(3) = .true. - -! Four-electron integrals - - read(1,*) - read(1,*) answer1,answer2,answer3 - - if(answer1 == 'T') do4eInt(1) = .true. - if(answer2 == 'T') do4eInt(2) = .true. - if(answer3 == 'T') do4eInt(3) = .true. - -! Close file with options - - close(unit=1) - -end subroutine read_options diff --git a/src/QuAcK/BSE2.f90 b/src/MBPT/BSE2.f90 similarity index 100% rename from src/QuAcK/BSE2.f90 rename to src/MBPT/BSE2.f90 diff --git a/src/QuAcK/BSE2_A_matrix_dynamic.f90 b/src/MBPT/BSE2_A_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/BSE2_A_matrix_dynamic.f90 rename to src/MBPT/BSE2_A_matrix_dynamic.f90 diff --git a/src/QuAcK/BSE2_B_matrix_dynamic.f90 b/src/MBPT/BSE2_B_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/BSE2_B_matrix_dynamic.f90 rename to src/MBPT/BSE2_B_matrix_dynamic.f90 diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/MBPT/BSE2_dynamic_perturbation.f90 similarity index 100% rename from src/QuAcK/BSE2_dynamic_perturbation.f90 rename to src/MBPT/BSE2_dynamic_perturbation.f90 diff --git a/src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 b/src/MBPT/BSE2_dynamic_perturbation_iterative.f90 similarity index 100% rename from src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 rename to src/MBPT/BSE2_dynamic_perturbation_iterative.f90 diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/MBPT/Bethe_Salpeter.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter.f90 rename to src/MBPT/Bethe_Salpeter.f90 diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/MBPT/Bethe_Salpeter_AB_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 rename to src/MBPT/Bethe_Salpeter_AB_matrix_dynamic.f90 diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 b/src/MBPT/Bethe_Salpeter_A_matrix.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_A_matrix.f90 rename to src/MBPT/Bethe_Salpeter_A_matrix.f90 diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 rename to src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 b/src/MBPT/Bethe_Salpeter_B_matrix.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_B_matrix.f90 rename to src/MBPT/Bethe_Salpeter_B_matrix.f90 diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/MBPT/Bethe_Salpeter_ZAB_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 rename to src/MBPT/Bethe_Salpeter_ZAB_matrix_dynamic.f90 diff --git a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/MBPT/Bethe_Salpeter_ZA_matrix_dynamic.f90 similarity index 100% rename from src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 rename to src/MBPT/Bethe_Salpeter_ZA_matrix_dynamic.f90 diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 similarity index 95% rename from src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 rename to src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 index 197f2b2..c960a75 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, & - OmRPA,rho_RPA,OmBSE,XpY,XmY) +subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -55,10 +54,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) -! Print main components of transition vectors - - call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) - if(dTDA) then write(*,*) write(*,*) '*** dynamical TDA activated ***' diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 similarity index 97% rename from src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 rename to src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index c25da2a..bdebd2a 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -54,10 +54,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, if(.not.dTDA) allocate(Am_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) -! Print main components of transition vectors - - call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) - if(dTDA) then write(*,*) write(*,*) '*** dynamical TDA activated ***' diff --git a/src/QuAcK/G0F2.f90 b/src/MBPT/G0F2.f90 similarity index 100% rename from src/QuAcK/G0F2.f90 rename to src/MBPT/G0F2.f90 diff --git a/src/QuAcK/G0F3.f90 b/src/MBPT/G0F3.f90 similarity index 100% rename from src/QuAcK/G0F3.f90 rename to src/MBPT/G0F3.f90 diff --git a/src/QuAcK/G0T0.f90 b/src/MBPT/G0T0.f90 similarity index 100% rename from src/QuAcK/G0T0.f90 rename to src/MBPT/G0T0.f90 diff --git a/src/QuAcK/G0W0.f90 b/src/MBPT/G0W0.f90 similarity index 100% rename from src/QuAcK/G0W0.f90 rename to src/MBPT/G0W0.f90 diff --git a/src/MBPT/Makefile b/src/MBPT/Makefile new file mode 100644 index 0000000..3d135d7 --- /dev/null +++ b/src/MBPT/Makefile @@ -0,0 +1,2 @@ +TARGET=MBPT.a +include ../Makefile.include diff --git a/src/QuAcK/QP_graph.f90 b/src/MBPT/QP_graph.f90 similarity index 100% rename from src/QuAcK/QP_graph.f90 rename to src/MBPT/QP_graph.f90 diff --git a/src/QuAcK/QP_roots.f90 b/src/MBPT/QP_roots.f90 similarity index 100% rename from src/QuAcK/QP_roots.f90 rename to src/MBPT/QP_roots.f90 diff --git a/src/MBPT/Sangalli_dynamic_perturbation.f90 b/src/MBPT/Sangalli_dynamic_perturbation.f90 new file mode 100644 index 0000000..9718415 --- /dev/null +++ b/src/MBPT/Sangalli_dynamic_perturbation.f90 @@ -0,0 +1,124 @@ +subroutine Sangalli_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) + +! Compute dynamical effects via perturbation theory for Sangalli's kernel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dTDA + 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 + integer,intent(in) :: nS + + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: OmBSE(nS) + double precision,intent(in) :: XpY(nS,nS) + double precision,intent(in) :: XmY(nS,nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: ia + + integer,parameter :: maxS = 10 + double precision :: gapGW + + double precision,allocatable :: OmDyn(:) + double precision,allocatable :: ZDyn(:) + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + + double precision,allocatable :: Ap_dyn(:,:) + double precision,allocatable :: ZAp_dyn(:,:) + + double precision,allocatable :: Bp_dyn(:,:) + double precision,allocatable :: ZBp_dyn(:,:) + + double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: ZAm_dyn(:,:) + + double precision,allocatable :: Bm_dyn(:,:) + double precision,allocatable :: ZBm_dyn(:,:) + +! Memory allocation + + allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) + + if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) + +! Print main components of transition vectors + + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + + if(dTDA) then + write(*,*) + write(*,*) '*** dynamical TDA activated ***' + write(*,*) + end if + + gapGW = eGW(nO+1) - eGW(nO) + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static BSE excitation energies with Sangalli kernel ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do ia=1,min(nS,maxS) + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + if(dTDA) then + +! call Sangalli_A_matrix_dynamic (eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) +! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) + + ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) + + else + +! call Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) +! call Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),Am_dyn(:,:)) + +! call Sangalli_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Bp_dyn(:,:)) +! call Sangalli_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),Bm_dyn(:,:)) + +! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) +! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),ZAm_dyn(:,:)) +! call Sangalli_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),ZBp_dyn(:,)) +! call Sangalli_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),ZBm_dyn(:,)) + + ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & + + dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & + + dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) + + OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & + - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & + - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) + + end if + + ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) + OmDyn(ia) = ZDyn(ia)*OmDyn(ia) + + write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & + ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine Sangalli_dynamic_perturbation diff --git a/src/QuAcK/SigmaC.f90 b/src/MBPT/SigmaC.f90 similarity index 100% rename from src/QuAcK/SigmaC.f90 rename to src/MBPT/SigmaC.f90 diff --git a/src/QuAcK/UG0W0.f90 b/src/MBPT/UG0W0.f90 similarity index 100% rename from src/QuAcK/UG0W0.f90 rename to src/MBPT/UG0W0.f90 diff --git a/src/QuAcK/USigmaC.f90 b/src/MBPT/USigmaC.f90 similarity index 100% rename from src/QuAcK/USigmaC.f90 rename to src/MBPT/USigmaC.f90 diff --git a/src/QuAcK/dSigmaC.f90 b/src/MBPT/dSigmaC.f90 similarity index 100% rename from src/QuAcK/dSigmaC.f90 rename to src/MBPT/dSigmaC.f90 diff --git a/src/QuAcK/dUSigmaC.f90 b/src/MBPT/dUSigmaC.f90 similarity index 100% rename from src/QuAcK/dUSigmaC.f90 rename to src/MBPT/dUSigmaC.f90 diff --git a/src/QuAcK/evGF2.f90 b/src/MBPT/evGF2.f90 similarity index 100% rename from src/QuAcK/evGF2.f90 rename to src/MBPT/evGF2.f90 diff --git a/src/QuAcK/evGF3.f90 b/src/MBPT/evGF3.f90 similarity index 100% rename from src/QuAcK/evGF3.f90 rename to src/MBPT/evGF3.f90 diff --git a/src/QuAcK/evGT.f90 b/src/MBPT/evGT.f90 similarity index 100% rename from src/QuAcK/evGT.f90 rename to src/MBPT/evGT.f90 diff --git a/src/QuAcK/evGW.f90 b/src/MBPT/evGW.f90 similarity index 100% rename from src/QuAcK/evGW.f90 rename to src/MBPT/evGW.f90 diff --git a/src/QuAcK/evUGW.f90 b/src/MBPT/evUGW.f90 similarity index 100% rename from src/QuAcK/evUGW.f90 rename to src/MBPT/evUGW.f90 diff --git a/src/QuAcK/exchange_matrix_MO_basis.f90 b/src/MBPT/exchange_matrix_MO_basis.f90 similarity index 100% rename from src/QuAcK/exchange_matrix_MO_basis.f90 rename to src/MBPT/exchange_matrix_MO_basis.f90 diff --git a/src/QuAcK/excitation_density.f90 b/src/MBPT/excitation_density.f90 similarity index 100% rename from src/QuAcK/excitation_density.f90 rename to src/MBPT/excitation_density.f90 diff --git a/src/QuAcK/excitation_density_RI.f90 b/src/MBPT/excitation_density_RI.f90 similarity index 100% rename from src/QuAcK/excitation_density_RI.f90 rename to src/MBPT/excitation_density_RI.f90 diff --git a/src/QuAcK/excitation_density_SOSEX.f90 b/src/MBPT/excitation_density_SOSEX.f90 similarity index 100% rename from src/QuAcK/excitation_density_SOSEX.f90 rename to src/MBPT/excitation_density_SOSEX.f90 diff --git a/src/QuAcK/excitation_density_SOSEX_RI.f90 b/src/MBPT/excitation_density_SOSEX_RI.f90 similarity index 100% rename from src/QuAcK/excitation_density_SOSEX_RI.f90 rename to src/MBPT/excitation_density_SOSEX_RI.f90 diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/MBPT/excitation_density_Tmatrix.f90 similarity index 100% rename from src/QuAcK/excitation_density_Tmatrix.f90 rename to src/MBPT/excitation_density_Tmatrix.f90 diff --git a/src/QuAcK/plot_GW.f90 b/src/MBPT/plot_GW.f90 similarity index 100% rename from src/QuAcK/plot_GW.f90 rename to src/MBPT/plot_GW.f90 diff --git a/src/QuAcK/print_G0F2.f90 b/src/MBPT/print_G0F2.f90 similarity index 100% rename from src/QuAcK/print_G0F2.f90 rename to src/MBPT/print_G0F2.f90 diff --git a/src/QuAcK/print_G0F3.f90 b/src/MBPT/print_G0F3.f90 similarity index 100% rename from src/QuAcK/print_G0F3.f90 rename to src/MBPT/print_G0F3.f90 diff --git a/src/QuAcK/print_G0T0.f90 b/src/MBPT/print_G0T0.f90 similarity index 100% rename from src/QuAcK/print_G0T0.f90 rename to src/MBPT/print_G0T0.f90 diff --git a/src/QuAcK/print_G0W0.f90 b/src/MBPT/print_G0W0.f90 similarity index 100% rename from src/QuAcK/print_G0W0.f90 rename to src/MBPT/print_G0W0.f90 diff --git a/src/QuAcK/print_UG0W0.f90 b/src/MBPT/print_UG0W0.f90 similarity index 100% rename from src/QuAcK/print_UG0W0.f90 rename to src/MBPT/print_UG0W0.f90 diff --git a/src/QuAcK/print_evGF2.f90 b/src/MBPT/print_evGF2.f90 similarity index 100% rename from src/QuAcK/print_evGF2.f90 rename to src/MBPT/print_evGF2.f90 diff --git a/src/QuAcK/print_evGF3.f90 b/src/MBPT/print_evGF3.f90 similarity index 100% rename from src/QuAcK/print_evGF3.f90 rename to src/MBPT/print_evGF3.f90 diff --git a/src/QuAcK/print_evGT.f90 b/src/MBPT/print_evGT.f90 similarity index 100% rename from src/QuAcK/print_evGT.f90 rename to src/MBPT/print_evGT.f90 diff --git a/src/QuAcK/print_evGW.f90 b/src/MBPT/print_evGW.f90 similarity index 100% rename from src/QuAcK/print_evGW.f90 rename to src/MBPT/print_evGW.f90 diff --git a/src/QuAcK/print_evUGW.f90 b/src/MBPT/print_evUGW.f90 similarity index 100% rename from src/QuAcK/print_evUGW.f90 rename to src/MBPT/print_evUGW.f90 diff --git a/src/QuAcK/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 similarity index 100% rename from src/QuAcK/print_qsGW.f90 rename to src/MBPT/print_qsGW.f90 diff --git a/src/QuAcK/qsGW.f90 b/src/MBPT/qsGW.f90 similarity index 100% rename from src/QuAcK/qsGW.f90 rename to src/MBPT/qsGW.f90 diff --git a/src/QuAcK/qsGW_PT.f90 b/src/MBPT/qsGW_PT.f90 similarity index 100% rename from src/QuAcK/qsGW_PT.f90 rename to src/MBPT/qsGW_PT.f90 diff --git a/src/QuAcK/renormalization_factor.f90 b/src/MBPT/renormalization_factor.f90 similarity index 100% rename from src/QuAcK/renormalization_factor.f90 rename to src/MBPT/renormalization_factor.f90 diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/MBPT/renormalization_factor_Tmatrix.f90 similarity index 100% rename from src/QuAcK/renormalization_factor_Tmatrix.f90 rename to src/MBPT/renormalization_factor_Tmatrix.f90 diff --git a/src/QuAcK/self_energy_Tmatrix.f90 b/src/MBPT/self_energy_Tmatrix.f90 similarity index 100% rename from src/QuAcK/self_energy_Tmatrix.f90 rename to src/MBPT/self_energy_Tmatrix.f90 diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/MBPT/self_energy_Tmatrix_diag.f90 similarity index 100% rename from src/QuAcK/self_energy_Tmatrix_diag.f90 rename to src/MBPT/self_energy_Tmatrix_diag.f90 diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/MBPT/self_energy_correlation.f90 similarity index 100% rename from src/QuAcK/self_energy_correlation.f90 rename to src/MBPT/self_energy_correlation.f90 diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/MBPT/self_energy_correlation_diag.f90 similarity index 100% rename from src/QuAcK/self_energy_correlation_diag.f90 rename to src/MBPT/self_energy_correlation_diag.f90 diff --git a/src/QuAcK/self_energy_exchange.f90 b/src/MBPT/self_energy_exchange.f90 similarity index 100% rename from src/QuAcK/self_energy_exchange.f90 rename to src/MBPT/self_energy_exchange.f90 diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/MBPT/unrestricted_Bethe_Salpeter.f90 similarity index 79% rename from src/QuAcK/unrestricted_Bethe_Salpeter.f90 rename to src/MBPT/unrestricted_Bethe_Salpeter.f90 index 1b7566a..0ce2590 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter.f90 @@ -105,22 +105,13 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Compute the dynamical screening at the BSE level !------------------------------------------------- -! if(dBSE) then + if(dBSE) & + call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) -! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) -! -! if(evDyn) then -! -! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & -! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) -! else -! -! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & -! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) -! end if + deallocate(OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) -! end if - end if !-----------------------! @@ -150,22 +141,13 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ! Compute the dynamical screening at the BSE level !------------------------------------------------- -! if(dBSE) then - -! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - -! if(evDyn) then -! -! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & -! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) -! else -! -! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & -! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) -! end if - -! end if + if(dBSE) & + call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sf,nS_sc, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) + deallocate(OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) + end if end subroutine unrestricted_Bethe_Salpeter diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 similarity index 97% rename from src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 rename to src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 index 1b1e415..25e81c5 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & +subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -20,6 +20,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n integer,intent(in) :: nS_sc double precision,intent(in) :: eta double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 new file mode 100644 index 0000000..b985f23 --- /dev/null +++ b/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 @@ -0,0 +1,204 @@ +subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,OmBSE,A_dyn) + +! Compute the extra term for dynamical Bethe-Salpeter equation for linear response in the unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + 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) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + integer,intent(in) :: nS_sc + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: OmRPA(nS_sc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: A_dyn(nSt,nSt) + +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! + + A_dyn(:,:) = 0d0 + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + 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 + + A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*eps/(eps**2 + eta**2) + + 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 + + A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi + + enddo + enddo + enddo + enddo + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + 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 + + A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*eps/(eps**2 + eta**2) + + 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 + + A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi + + enddo + enddo + enddo + enddo + + end if + +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! + + if(ispin == 2) then + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + 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 + + A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*eps/(eps**2 + eta**2) + + eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*eps/(eps**2 + eta**2) + + enddo + + A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + 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 + + A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*eps/(eps**2 + eta**2) + + eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*eps/(eps**2 + eta**2) + + enddo + + A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine unrestricted_Bethe_Salpeter_A_matrix_dynamic diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 similarity index 100% rename from src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 rename to src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 diff --git a/src/QuAcK/unrestricted_QP_graph.f90 b/src/MBPT/unrestricted_QP_graph.f90 similarity index 100% rename from src/QuAcK/unrestricted_QP_graph.f90 rename to src/MBPT/unrestricted_QP_graph.f90 diff --git a/src/QuAcK/unrestricted_excitation_density.f90 b/src/MBPT/unrestricted_excitation_density.f90 similarity index 100% rename from src/QuAcK/unrestricted_excitation_density.f90 rename to src/MBPT/unrestricted_excitation_density.f90 diff --git a/src/QuAcK/unrestricted_renormalization_factor.f90 b/src/MBPT/unrestricted_renormalization_factor.f90 similarity index 100% rename from src/QuAcK/unrestricted_renormalization_factor.f90 rename to src/MBPT/unrestricted_renormalization_factor.f90 diff --git a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 b/src/MBPT/unrestricted_self_energy_correlation_diag.f90 similarity index 100% rename from src/QuAcK/unrestricted_self_energy_correlation_diag.f90 rename to src/MBPT/unrestricted_self_energy_correlation_diag.f90 diff --git a/src/QuAcK/AO_values.f90 b/src/MC/AO_values.f90 similarity index 100% rename from src/QuAcK/AO_values.f90 rename to src/MC/AO_values.f90 diff --git a/src/QuAcK/AO_values_grid.f90 b/src/MC/AO_values_grid.f90 similarity index 100% rename from src/QuAcK/AO_values_grid.f90 rename to src/MC/AO_values_grid.f90 diff --git a/src/QuAcK/Green_function.f90 b/src/MC/Green_function.f90 similarity index 100% rename from src/QuAcK/Green_function.f90 rename to src/MC/Green_function.f90 diff --git a/src/QuAcK/MCMP2.f90 b/src/MC/MCMP2.f90 similarity index 100% rename from src/QuAcK/MCMP2.f90 rename to src/MC/MCMP2.f90 diff --git a/src/QuAcK/MO_values_grid.f90 b/src/MC/MO_values_grid.f90 similarity index 100% rename from src/QuAcK/MO_values_grid.f90 rename to src/MC/MO_values_grid.f90 diff --git a/src/MC/Makefile b/src/MC/Makefile new file mode 100644 index 0000000..656f285 --- /dev/null +++ b/src/MC/Makefile @@ -0,0 +1,2 @@ +TARGET=MC.a +include ../Makefile.include diff --git a/src/QuAcK/MinMCMP2.f90 b/src/MC/MinMCMP2.f90 similarity index 100% rename from src/QuAcK/MinMCMP2.f90 rename to src/MC/MinMCMP2.f90 diff --git a/src/QuAcK/NDrift.f90 b/src/MC/NDrift.f90 similarity index 100% rename from src/QuAcK/NDrift.f90 rename to src/MC/NDrift.f90 diff --git a/src/QuAcK/Newton.f90 b/src/MC/Newton.f90 similarity index 100% rename from src/QuAcK/Newton.f90 rename to src/MC/Newton.f90 diff --git a/src/QuAcK/drift.f90 b/src/MC/drift.f90 similarity index 100% rename from src/QuAcK/drift.f90 rename to src/MC/drift.f90 diff --git a/src/QuAcK/generate_shell.f90 b/src/MC/generate_shell.f90 similarity index 100% rename from src/QuAcK/generate_shell.f90 rename to src/MC/generate_shell.f90 diff --git a/src/QuAcK/initialize_random_generator.f90 b/src/MC/initialize_random_generator.f90 similarity index 100% rename from src/QuAcK/initialize_random_generator.f90 rename to src/MC/initialize_random_generator.f90 diff --git a/src/QuAcK/norm_trial.f90 b/src/MC/norm_trial.f90 similarity index 100% rename from src/QuAcK/norm_trial.f90 rename to src/MC/norm_trial.f90 diff --git a/src/QuAcK/optimize_timestep.f90 b/src/MC/optimize_timestep.f90 similarity index 100% rename from src/QuAcK/optimize_timestep.f90 rename to src/MC/optimize_timestep.f90 diff --git a/src/QuAcK/transition_probability.f90 b/src/MC/transition_probability.f90 similarity index 100% rename from src/QuAcK/transition_probability.f90 rename to src/MC/transition_probability.f90 diff --git a/src/QuAcK/MP2.f90 b/src/MP/MP2.f90 similarity index 100% rename from src/QuAcK/MP2.f90 rename to src/MP/MP2.f90 diff --git a/src/QuAcK/MP2F12.f90 b/src/MP/MP2F12.f90 similarity index 100% rename from src/QuAcK/MP2F12.f90 rename to src/MP/MP2F12.f90 diff --git a/src/QuAcK/MP3.f90 b/src/MP/MP3.f90 similarity index 100% rename from src/QuAcK/MP3.f90 rename to src/MP/MP3.f90 diff --git a/src/MP/Makefile b/src/MP/Makefile new file mode 100644 index 0000000..8e0c07c --- /dev/null +++ b/src/MP/Makefile @@ -0,0 +1,2 @@ +TARGET=MP.a +include ../Makefile.include diff --git a/src/QuAcK/UMP2.f90 b/src/MP/UMP2.f90 similarity index 100% rename from src/QuAcK/UMP2.f90 rename to src/MP/UMP2.f90 diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..abc97eb --- /dev/null +++ b/src/Makefile @@ -0,0 +1,53 @@ +ALL_DIRS=$(filter-out modules,$(patsubst %/,%,$(wildcard */))) + +# Rules for Modules +################### + +MOD_DIRS=numgrid + + +# Rules for Libraries +##################### + +LDIR =../lib +LIB_DIRS=$(filter-out $(MAIN_DIRS), $(ALL_DIRS)) +MAKEFILES=$(patsubst %,%/Makefile, $(ALL_DIRS)) + +FORCE: + +%/Makefile: + @(echo TARGET=$*.a > $*/Makefile && echo include ../Makefile.include >> $@) + +$(LDIR)/%.a: FORCE %/Makefile $(MOD_DIRS) + $(MAKE) -C $* lib + + + +# Rules for executables +####################### + +MAIN_DIRS=QuAcK eDFT + +BDIR=../bin +ALL_EXEC=$(patsubst %, $(BDIR)/%, $(MAIN_DIRS) ) + +$(BDIR)/%: $(patsubst %,$(LDIR)/%.a,$(LIB_DIRS)) %/Makefile + $(MAKE) -C $* ../$@ + + + + +# Rules for both +################ + +default: $(ALL_EXEC) + +.DEFAULT_GOAL := default +.PHONY: default +.PRECIOUS: $(LDIR)/%.a %/Makefile + +clean: + rm -f -- $(patsubst %/,$(LDIR)/%.a,$(wildcard */)) ; \ + rm -f -- $(ALL_EXEC) ; \ + for i in $(MAIN_DIRS) $(LIB_DIRS) ; do rm -f -- $$i/obj/* ; done + diff --git a/src/Makefile.include b/src/Makefile.include new file mode 100644 index 0000000..4e1fb3a --- /dev/null +++ b/src/Makefile.include @@ -0,0 +1,42 @@ +BDIR =../../bin +LDIR =../../lib +ODIR = obj +SDIR =. +include ../Makefile.common + +LIBS = $(filter-out $(LDIR)/$(TARGET), $(wildcard $(LDIR)/*.a)) +LIBS += -lblas -llapack -lc++ + +SRCF90 = $(wildcard *.f90) + +SRC = $(wildcard *.f) + +OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC)) + +FORCE: + +$(ODIR)/%.o: %.f90 + $(FC) -c -o $@ $< $(FFLAGS) + +$(ODIR)/%.o: %.f + $(FC) -c -o $@ $< $(FFLAGS) + +$(LDIR)/$(TARGET): $(patsubst %,$(LDIR)/%,$(DEPEND)) $(OBJ) + $(AR) -static -o $@ $^ + +debug: + DEBUG=1 make $(LDIR)/$(TARGET) + +clean: + rm -f $(ODIR)/*.o $(LDIR)/$(TARGET) + +$(BDIR)/%: $(OBJ) FORCE + $(FC) -o $@ $(FFLAGS) $(LIBS) $(OBJ) + +default: + $(MAKE) -C .. + +lib: $(LDIR)/$(TARGET) + +.DEFAULT_GOAL := default +.PHONY: lib default diff --git a/src/QuAcK/Makefile b/src/QuAcK/Makefile index 8ba0980..235fdec 100644 --- a/src/QuAcK/Makefile +++ b/src/QuAcK/Makefile @@ -1,43 +1,2 @@ -IDIR =../../include -BDIR =../../bin -ODIR = obj -OODIR = ../IntPak/obj -SDIR =. -FC = gfortran -I$(IDIR) -Wall -g -Wno-unused -Wno-unused-dummy-argument -ifeq ($(DEBUG),1) - FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -else - FFLAGS = -O3 -endif - -ifeq ($(PROFIL),1) - FC += -p -fno-inline -endif - -LIBS = ../../lib/*.a -LIBS += -lblas -llapack - -SRCF90 = $(wildcard *.f90) - -SRC = $(wildcard *.f) - -OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC)) - -$(ODIR)/%.o: %.f90 - $(FC) -c -o $@ $< $(FFLAGS) - -$(ODIR)/%.o: %.f - $(FC) -c -o $@ $< $(FFLAGS) - -$(BDIR)/QuAcK: $(OBJ) - $(FC) -o $@ $^ $(FFLAGS) $(LIBS) - -debug: - DEBUG=1 make $(BDIR)/QuAcK -#DEBUG=1 make clean $(BDIR)/QuAcK - -profil: - PROFIL=1 make $(BDIR)/QuAcK - -clean: - rm -f $(ODIR)/*.o $(BDIR)/QuAcK $(BDIR)/debug +TARGET=QuAcK.a +include ../Makefile.include diff --git a/src/QuAcK/dcgw.f90 b/src/QuAcK/dcgw.f90 deleted file mode 100644 index e67ff3b..0000000 --- a/src/QuAcK/dcgw.f90 +++ /dev/null @@ -1,84 +0,0 @@ -function SigC_dcgw(x,y) result(SigC) - -! Degeneracy-corrected GW - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: x,y - -! Local variables - - double precision,parameter :: eta = 0.1d0 - double precision :: r - -! Output variables - - double precision :: SigC - -! Compute the divergence-free term - - r = y/x - -! Bare style - - SigC = y*r - -! DCPT style - -! SigC = -0.5d0*x*(1d0-sqrt(1d0+4d0*r*r)) - -! QDPT style - -! SigC = y*r/sqrt(1d0+4d0*r*r) - -! Infinitesimal - -! SigC = y*y*x/(x*x+eta*eta) - -end function SigC_dcgw - -function Z_dcgw(x,y) result(Z) - - -! Derivative of the degeneracy-corrected GW - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: x,y - -! Local variables - - double precision,parameter :: eta = 1d0 - double precision :: r - -! Output variables - - double precision :: Z - -! Compute the divergence-free term - - r = y/x - -! Bare style - - Z = r*r - -! DCPT style - -! Z = 0.5d0*(1d0-1d0/sqrt(1d0+4d0*r*r)) - -! QDPT style - -! Z = r/sqrt(1d0+4d0*r*r)/(1d0+4d0*r*r) - -! Infinitesimal - -! Z = y*y*(x*x-eta*eta)/(x*x+eta*eta)**2 - -end function Z_dcgw diff --git a/src/QuAcK/ACFDT.f90 b/src/RPA/ACFDT.f90 similarity index 100% rename from src/QuAcK/ACFDT.f90 rename to src/RPA/ACFDT.f90 diff --git a/src/QuAcK/ACFDT_correlation_energy.f90 b/src/RPA/ACFDT_correlation_energy.f90 similarity index 100% rename from src/QuAcK/ACFDT_correlation_energy.f90 rename to src/RPA/ACFDT_correlation_energy.f90 diff --git a/src/RPA/Makefile b/src/RPA/Makefile new file mode 100644 index 0000000..79e1bcd --- /dev/null +++ b/src/RPA/Makefile @@ -0,0 +1,2 @@ +TARGET=RPA.a +include ../Makefile.include diff --git a/src/QuAcK/RPAx.f90 b/src/RPA/RPAx.f90 similarity index 100% rename from src/QuAcK/RPAx.f90 rename to src/RPA/RPAx.f90 diff --git a/src/QuAcK/URPAx.f90 b/src/RPA/URPAx.f90 similarity index 100% rename from src/QuAcK/URPAx.f90 rename to src/RPA/URPAx.f90 diff --git a/src/QuAcK/UdRPA.f90 b/src/RPA/UdRPA.f90 similarity index 100% rename from src/QuAcK/UdRPA.f90 rename to src/RPA/UdRPA.f90 diff --git a/src/QuAcK/dRPA.f90 b/src/RPA/dRPA.f90 similarity index 100% rename from src/QuAcK/dRPA.f90 rename to src/RPA/dRPA.f90 diff --git a/src/QuAcK/linear_response.f90 b/src/RPA/linear_response.f90 similarity index 100% rename from src/QuAcK/linear_response.f90 rename to src/RPA/linear_response.f90 diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/RPA/linear_response_A_matrix.f90 similarity index 100% rename from src/QuAcK/linear_response_A_matrix.f90 rename to src/RPA/linear_response_A_matrix.f90 diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/RPA/linear_response_B_matrix.f90 similarity index 100% rename from src/QuAcK/linear_response_B_matrix.f90 rename to src/RPA/linear_response_B_matrix.f90 diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/RPA/linear_response_B_pp.f90 similarity index 100% rename from src/QuAcK/linear_response_B_pp.f90 rename to src/RPA/linear_response_B_pp.f90 diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/RPA/linear_response_C_pp.f90 similarity index 100% rename from src/QuAcK/linear_response_C_pp.f90 rename to src/RPA/linear_response_C_pp.f90 diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/RPA/linear_response_D_pp.f90 similarity index 100% rename from src/QuAcK/linear_response_D_pp.f90 rename to src/RPA/linear_response_D_pp.f90 diff --git a/src/QuAcK/linear_response_ph.f90 b/src/RPA/linear_response_ph.f90 similarity index 100% rename from src/QuAcK/linear_response_ph.f90 rename to src/RPA/linear_response_ph.f90 diff --git a/src/QuAcK/linear_response_pp.f90 b/src/RPA/linear_response_pp.f90 similarity index 100% rename from src/QuAcK/linear_response_pp.f90 rename to src/RPA/linear_response_pp.f90 diff --git a/src/QuAcK/ppRPA.f90 b/src/RPA/ppRPA.f90 similarity index 100% rename from src/QuAcK/ppRPA.f90 rename to src/RPA/ppRPA.f90 diff --git a/src/QuAcK/print_excitation.f90 b/src/RPA/print_excitation.f90 similarity index 100% rename from src/QuAcK/print_excitation.f90 rename to src/RPA/print_excitation.f90 diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/RPA/print_transition_vectors.f90 similarity index 100% rename from src/QuAcK/print_transition_vectors.f90 rename to src/RPA/print_transition_vectors.f90 diff --git a/src/QuAcK/print_unrestricted_transition_vectors.f90 b/src/RPA/print_unrestricted_transition_vectors.f90 similarity index 100% rename from src/QuAcK/print_unrestricted_transition_vectors.f90 rename to src/RPA/print_unrestricted_transition_vectors.f90 diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 similarity index 100% rename from src/QuAcK/sort_ppRPA.f90 rename to src/RPA/sort_ppRPA.f90 diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/RPA/unrestricted_linear_response.f90 similarity index 99% rename from src/QuAcK/unrestricted_linear_response.f90 rename to src/RPA/unrestricted_linear_response.f90 index e816162..422a863 100644 --- a/src/QuAcK/unrestricted_linear_response.f90 +++ b/src/RPA/unrestricted_linear_response.f90 @@ -61,7 +61,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A) if(BSE) & - call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,A) ! Tamm-Dancoff approximation diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/RPA/unrestricted_linear_response_A_matrix.f90 similarity index 100% rename from src/QuAcK/unrestricted_linear_response_A_matrix.f90 rename to src/RPA/unrestricted_linear_response_A_matrix.f90 diff --git a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 b/src/RPA/unrestricted_linear_response_B_matrix.f90 similarity index 100% rename from src/QuAcK/unrestricted_linear_response_B_matrix.f90 rename to src/RPA/unrestricted_linear_response_B_matrix.f90 diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile index 3efafe3..795c0b4 100644 --- a/src/eDFT/Makefile +++ b/src/eDFT/Makefile @@ -1,45 +1,2 @@ -IDIR =../../include -BDIR =../../bin -ODIR = obj -OODIR = ../IntPak/obj -SDIR =. -FC = gfortran -I$(IDIR) -ifeq ($(DEBUG),1) -FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -else -FFLAGS = -O3 -endif - -##ifeq ($(PROFIL),1) -## FC += -p -fno-inline -##endif - -LIBS = ../../lib/*.a -LIBS += -lblas -llapack -lc++ - -SRCF90 = $(wildcard *.f90) - -SRC = $(wildcard *.F) - -OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.F,$(ODIR)/%.o,$(SRC)) - -$(BDIR)/eDFT: $(OBJ) - $(FC) -o $@ $^ $(FFLAGS) $(LIBS) - -numgrid.mod $(ODIR)/numgrid.o: numgrid.f90 - $(FC) -c -o $(ODIR)/numgrid.o $< $(FFLAGS) - -$(ODIR)/%.o: %.f90 numgrid.mod - $(FC) -c -o $@ $< $(FFLAGS) - -$(ODIR)/%.o: %.F numgrid.mod - $(FC) -c -o $@ $< $(FFLAGS) - -debug: - DEBUG=1 make $(BDIR)/eDFT - -profil: - PROFIL=1 make $(BDIR)/eDFT - -clean: - rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug numgrid.mod +TARGET=eDFT.a +include ../Makefile.include diff --git a/src/eDFT/numgrid.f90 b/src/numgrid/numgrid.f90 similarity index 100% rename from src/eDFT/numgrid.f90 rename to src/numgrid/numgrid.f90 diff --git a/src/numgrid/obj/.gitignore b/src/numgrid/obj/.gitignore new file mode 100644 index 0000000..e69de29 diff --git a/src/utils/Makefile b/src/utils/Makefile index 55227cc..530241b 100644 --- a/src/utils/Makefile +++ b/src/utils/Makefile @@ -1,37 +1,2 @@ -IDIR =../../include -BDIR =../../bin -LDIR =../../lib -ODIR = obj -OODIR = ../IntPak/obj -SDIR =. -FC = gfortran -I$(IDIR) -AR = libtool -ifeq ($(DEBUG),1) -FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -else -FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3 -endif - -LIBS = ../../lib/*.a -LIBS += -lblas -llapack - -SRCF90 = $(wildcard *.f90) - -SRC = $(wildcard *.f) - -OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC)) - -$(ODIR)/%.o: %.f90 - $(FC) -c -o $@ $< $(FFLAGS) - -$(ODIR)/%.o: %.f - $(FC) -c -o $@ $< $(FFLAGS) - -$(LDIR)/utils.a: $(OBJ) - $(AR) -static -o $@ $^ - -debug: - DEBUG=1 make $(LDIR)/utils.a - -clean: - rm -f $(ODIR)/*.o $(LDIR)/utils.a +TARGET=utils.a +include ../Makefile.include diff --git a/src/QuAcK/chem_to_phys_ERI.f90 b/src/utils/chem_to_phys_ERI.f90 similarity index 100% rename from src/QuAcK/chem_to_phys_ERI.f90 rename to src/utils/chem_to_phys_ERI.f90 diff --git a/src/QuAcK/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 similarity index 100% rename from src/QuAcK/orthogonalization_matrix.f90 rename to src/utils/orthogonalization_matrix.f90 diff --git a/src/QuAcK/overlap.f90 b/src/utils/overlap.f90 similarity index 100% rename from src/QuAcK/overlap.f90 rename to src/utils/overlap.f90 diff --git a/src/QuAcK/phys_to_chem_ERI.f90 b/src/utils/phys_to_chem_ERI.f90 similarity index 100% rename from src/QuAcK/phys_to_chem_ERI.f90 rename to src/utils/phys_to_chem_ERI.f90 diff --git a/src/QuAcK/read_MOs.f90 b/src/utils/read_MOs.f90 similarity index 100% rename from src/QuAcK/read_MOs.f90 rename to src/utils/read_MOs.f90 diff --git a/src/utils/read_dipole_integrals.f90 b/src/utils/read_dipole_integrals.f90 new file mode 100644 index 0000000..03e5284 --- /dev/null +++ b/src/utils/read_dipole_integrals.f90 @@ -0,0 +1,71 @@ +subroutine read_dipole_integrals(nBas,R) + +! Read one-electron integrals related to dipole moment from files + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + +! Local variables + + logical :: debug = .false. + integer :: mu,nu + double precision :: Dip + +! Output variables + + double precision,intent(out) :: R(nBas,nBas,ncart) + +! Open file with integrals + + open(unit=21,file='int/x.dat') + open(unit=22,file='int/y.dat') + open(unit=23,file='int/z.dat') + +! Read (x,y,z) integrals + + R(:,:,:) = 0d0 + + do + read(21,*,end=21) mu,nu,Dip + R(mu,nu,1) = Dip + R(nu,mu,1) = Dip + enddo + 21 close(unit=21) + + do + read(22,*,end=22) mu,nu,Dip + R(mu,nu,2) = Dip + R(nu,mu,2) = Dip + enddo + 22 close(unit=22) + + do + read(23,*,end=23) mu,nu,Dip + R(mu,nu,3) = Dip + R(nu,mu,3) = Dip + enddo + 23 close(unit=23) + +! Print results + if(debug) then + write(*,'(A28)') '----------------------' + write(*,'(A28)') '(mu|x|nu) integrals' + write(*,'(A28)') '----------------------' + call matout(nBas,nBas,R(:,:,1)) + write(*,*) + write(*,'(A28)') '----------------------' + write(*,'(A28)') '(mu|y|nu) integrals' + write(*,'(A28)') '----------------------' + call matout(nBas,nBas,R(:,:,2)) + write(*,*) + write(*,'(A28)') '----------------------' + write(*,'(A28)') '(mu|z|nu) integrals' + write(*,'(A28)') '----------------------' + call matout(nBas,nBas,R(:,:,3)) + endif + +end subroutine read_dipole_integrals diff --git a/src/QuAcK/rij.f90 b/src/utils/rij.f90 similarity index 100% rename from src/QuAcK/rij.f90 rename to src/utils/rij.f90 diff --git a/src/QuAcK/spatial_to_spin_ERI.f90 b/src/utils/spatial_to_spin_ERI.f90 similarity index 100% rename from src/QuAcK/spatial_to_spin_ERI.f90 rename to src/utils/spatial_to_spin_ERI.f90 diff --git a/src/QuAcK/spatial_to_spin_MO_energy.f90 b/src/utils/spatial_to_spin_MO_energy.f90 similarity index 100% rename from src/QuAcK/spatial_to_spin_MO_energy.f90 rename to src/utils/spatial_to_spin_MO_energy.f90