diff --git a/input/basis b/input/basis index c12e4d6..b9ca7b5 100644 --- a/input/basis +++ b/input/basis @@ -1,20 +1,9 @@ -1 5 -S 6 1.00 - 1264.5857000 0.0019448 - 189.9368100 0.0148351 - 43.1590890 0.0720906 - 12.0986630 0.2371542 - 3.8063232 0.4691987 - 1.2728903 0.3565202 -S 3 1.00 - 3.1964631 -0.1126487 - 0.7478133 -0.2295064 - 0.2199663 1.1869167 -S 1 1.00 - 0.0823099 1.0000000 -P 3 1.00 - 3.1964631 0.0559802 - 0.7478133 0.2615506 - 0.2199663 0.7939723 -P 1 1.00 - 0.0823099 1.0000000 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 +S 1 1.00 + 0.2976000 1.0000000 +P 1 1.00 + 1.2750000 1.0000000 diff --git a/input/geminal b/input/geminal new file mode 100644 index 0000000..d3827e7 --- /dev/null +++ b/input/geminal @@ -0,0 +1 @@ +1.0 diff --git a/input/int b/input/int index 6f9394e..5fe40bc 100644 --- a/input/int +++ b/input/int @@ -1,13 +1,13 @@ # Debuggin mode? F # Chemist notation for two-electron integral? - T + F # Exposant of the Slater geminal 1.0 # One-electron integrals: Ov Kin Nuc - F F F + T T T # Two-electron integrals: ERI F12 Yuk Erf - T F F F + T T T T # Three-electron integrals: Type1 Type2 Type3 T F F # Four-electron integrals: Type1 Type2 Type3 diff --git a/input/methods b/input/methods index efe6a18..9eb4c24 100644 --- a/input/methods +++ b/input/methods @@ -1,13 +1,13 @@ # HF MOM T F # MP2 MP3 MP2-F12 - F F T + T F T # CIS TDHF ADC F F F # GF2 GF3 F F # G0W0 evGW qsGW - T F F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index 67ee9ac..7d017f4 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEl nCore nRyd - 1 4 0 0 + 1 2 0 0 # Znuc x y z - Be 0.0 0.0 0.0 + He 0.0 0.0 0.0 diff --git a/input/options b/input/options index 6d94128..2c2d513 100644 --- a/input/options +++ b/input/options @@ -1,7 +1,7 @@ # RHF: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.0000001 T 5 1 1 -# MP: MP2-F12 - T +# MP: + # CIS/TDHF: singlet triplet T F # GF: maxSCF thresh DIIS n_diis renormalization diff --git a/input/weight b/input/weight index c12e4d6..b9ca7b5 100644 --- a/input/weight +++ b/input/weight @@ -1,20 +1,9 @@ -1 5 -S 6 1.00 - 1264.5857000 0.0019448 - 189.9368100 0.0148351 - 43.1590890 0.0720906 - 12.0986630 0.2371542 - 3.8063232 0.4691987 - 1.2728903 0.3565202 -S 3 1.00 - 3.1964631 -0.1126487 - 0.7478133 -0.2295064 - 0.2199663 1.1869167 -S 1 1.00 - 0.0823099 1.0000000 -P 3 1.00 - 3.1964631 0.0559802 - 0.7478133 0.2615506 - 0.2199663 0.7939723 -P 1 1.00 - 0.0823099 1.0000000 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 +S 1 1.00 + 0.2976000 1.0000000 +P 1 1.00 + 1.2750000 1.0000000 diff --git a/src/IntPak/Compute3eInt.f90 b/src/IntPak/Compute3eInt.f90 index c787602..008fdca 100644 --- a/src/IntPak/Compute3eInt.f90 +++ b/src/IntPak/Compute3eInt.f90 @@ -266,11 +266,11 @@ subroutine Compute3eInt(debug,iType,nShell, & if(abs(c3eInt) > 1d-15) then nSigc3eInt = nSigc3eInt + 1 t_c3eInt = end_c3eInt - start_c3eInt - write(iFile,'(F20.15,I6,I6,I6,I6,I6,I6)') & - c3eInt,iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3 + write(iFile,'(I9,I9,I9,I9,I9,I9,F25.15)') & + iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3,c3eInt if(.true.) then - write(*,'(A15,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6)') & - '(a1a2a3|b1b2b3) = ',c3eInt,iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3 + 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 endif endif diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90 index fb23ec2..8bab7eb 100644 --- a/src/IntPak/IntPak.f90 +++ b/src/IntPak/IntPak.f90 @@ -19,7 +19,7 @@ program IntPak logical :: do4eInt(n4eInt) integer :: nNuc,nBas,iType - integer :: nEl,nO,nV,nCore,nRyd + integer :: nEl,nO,nV,nC,nR double precision :: ExpS double precision :: ENuc integer :: KG @@ -56,9 +56,6 @@ program IntPak call read_options(debug,chemist_notation,ExpS,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt) -! Which integrals do you want? - - !------------------------------------------------------------------------ ! Read input information !------------------------------------------------------------------------ @@ -69,7 +66,7 @@ program IntPak ! nBas = number of basis functions (see below) ! = nO + nV - call read_molecule(nNuc,nEl,nO,nCore,nRyd) + call read_molecule(nNuc,nEl,nO,nC,nR) allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3)) @@ -246,7 +243,7 @@ program IntPak call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,iType,nShell, & + call Compute2eInt(debug,chemist_notation,iType,nShell, & ExpS,KG,DG,ExpG, & CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) @@ -281,7 +278,7 @@ program IntPak ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,iType,nShell, & + call Compute2eInt(debug,chemist_notation,iType,nShell, & ExpS,KG,DG,ExpG, & CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) @@ -317,7 +314,7 @@ program IntPak ExpS = ExpS*ExpS call cpu_time(start_2eInt(iType)) - call Compute2eInt(debug,iType,nShell, & + call Compute2eInt(debug,chemist_notation,iType,nShell, & ExpS,KG,DG,ExpG, & CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & np2eInt(iType),nSigp2eInt(iType),nc2eInt(iType),nSigc2eInt(iType)) diff --git a/src/IntPak/elements.f90 b/src/IntPak/elements.f90 new file mode 100644 index 0000000..22953dc --- /dev/null +++ b/src/IntPak/elements.f90 @@ -0,0 +1,170 @@ +function element_number(element_name) + + implicit none + + integer,parameter :: nelement_max = 103 + character(len=2),intent(in) :: element_name + integer :: element_number + character(len=2),parameter :: element_list(nelement_max) = & + (/' H', 'He', & ! 2 + 'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10 + 'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18 + ' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36 + 'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54 + 'Cs','Ba', & ! 56 + 'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70 + 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86 + 'Fr','Ra', & ! 88 + 'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102 + 'Lr' & ! 103 + /) + +!===== + integer :: ielement +!===== + + ielement=1 + do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) ) + if( ielement == nelement_max ) then + write(*,'(a,a)') ' Input symbol ',element_name + write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements' + write(*,*) '!!! element symbol not understood !!!' + stop + endif + ielement = ielement + 1 + enddo + + element_number = ielement + +end function element_number + +function element_core(zval,zatom) + implicit none + double precision,intent(in) :: zval + double precision,intent(in) :: zatom + integer :: element_core +!===== + + ! + ! If zval /= zatom, this is certainly an effective core potential + ! and no core states should be frozen. + if( ABS(zval - zatom) > 1d0-3 ) then + element_core = 0 + else + + if( zval <= 4.00001d0 ) then ! up to Be + element_core = 0 + else if( zval <= 12.00001d0 ) then ! up to Mg + element_core = 1 + else if( zval <= 30.00001d0 ) then ! up to Ca + element_core = 5 + else if( zval <= 48.00001d0 ) then ! up to Sr + element_core = 9 + else + write(*,*) '!!! not imlemented in element_core !!!' + stop + endif + + endif + + +end function element_core + +function element_covalent_radius(zatom) + +! Return covalent radius of an atom + + implicit none + include 'parameters.h' + + integer,intent(in) :: zatom + double precision :: element_covalent_radius + + ! + ! Data from Cambridge Structural Database + ! http://en.wikipedia.org/wiki/Covalent_radius + ! + ! Values are first given in picometer + ! They will be converted in bohr later on + select case(zatom) + case( 1) + element_covalent_radius = 31. + case( 2) + element_covalent_radius = 28. + case( 3) + element_covalent_radius = 128. + case( 4) + element_covalent_radius = 96. + case( 5) + element_covalent_radius = 84. + case( 6) + element_covalent_radius = 73. + case( 7) + element_covalent_radius = 71. + case( 8) + element_covalent_radius = 66. + case( 9) + element_covalent_radius = 57. + case(10) ! Ne. + element_covalent_radius = 58. + case(11) + element_covalent_radius = 166. + case(12) + element_covalent_radius = 141. + case(13) + element_covalent_radius = 121. + case(14) + element_covalent_radius = 111. + case(15) + element_covalent_radius = 107. + case(16) + element_covalent_radius = 105. + case(17) + element_covalent_radius = 102. + case(18) ! Ar. + element_covalent_radius = 106. + case(19) + element_covalent_radius = 203. + case(20) + element_covalent_radius = 176. + case(21) + element_covalent_radius = 170. + case(22) + element_covalent_radius = 160. + case(23) + element_covalent_radius = 153. + case(24) + element_covalent_radius = 139. + case(25) + element_covalent_radius = 145. + case(26) + element_covalent_radius = 145. + case(27) + element_covalent_radius = 140. + case(28) + element_covalent_radius = 124. + case(29) + element_covalent_radius = 132. + case(30) + element_covalent_radius = 122. + case(31) + element_covalent_radius = 120. + case(32) + element_covalent_radius = 119. + case(34) + element_covalent_radius = 120. + case(35) + element_covalent_radius = 120. + case(36) ! Kr. + element_covalent_radius = 116. + case default + write(*,*) '!!! covalent radius not available !!!' + stop + end select + + ! pm to bohr conversion + element_covalent_radius = element_covalent_radius*pmtoau + + +end function element_covalent_radius + diff --git a/src/IntPak/read_auxiliary_basis.f90 b/src/IntPak/read_auxiliary_basis.f90 new file mode 100644 index 0000000..cf4f6ed --- /dev/null +++ b/src/IntPak/read_auxiliary_basis.f90 @@ -0,0 +1,176 @@ +subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & + TotAngMomShell,KShell,DShell,ExpShell) + +! Read auxiliary basis set information + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: NAtoms + double precision,intent(in) :: XYZAtoms(NAtoms,3) + +! Local variables + + integer :: nShAt,iAt + integer :: i,j,k + character :: shelltype + +! Output variables + + integer,intent(out) :: nShell + double precision,intent(out) :: CenterShell(maxShell,3) + integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell) + double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) + +!------------------------------------------------------------------------ +! Primary basis set information +!------------------------------------------------------------------------ + +! Open file with basis set specification + + open(unit=2,file='input/basis') + +! Read basis information + + write(*,'(A28)') 'Gaussian basis set' + write(*,'(A28)') '------------------' + + nShell = 0 + do i=1,NAtoms + read(2,*) iAt,nShAt + write(*,'(A28,1X,I16)') 'Atom n. ',iAt + write(*,'(A28,1X,I16)') 'number of shells ',nShAt + write(*,'(A28)') '------------------' + +! Basis function centers + + do j=1,nShAt + nShell = nShell + 1 + do k=1,3 + CenterShell(nShell,k) = XYZAtoms(iAt,k) + enddo + +! Shell type and contraction degree + + read(2,*) shelltype,KShell(nShell) + if(shelltype == "S") then + TotAngMomShell(nShell) = 0 + write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) + elseif(shelltype == "P") then + TotAngMomShell(nShell) = 1 + write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) + elseif(shelltype == "D") then + TotAngMomShell(nShell) = 2 + write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) + elseif(shelltype == "F") then + TotAngMomShell(nShell) = 3 + write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) + elseif(shelltype == "G") then + TotAngMomShell(nShell) = 4 + write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) + elseif(shelltype == "H") then + TotAngMomShell(nShell) = 5 + write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) + elseif(shelltype == "I") then + TotAngMomShell(nShell) = 6 + write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) + endif + +! Read exponents and contraction coefficients + + write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' + do k=1,Kshell(nShell) + read(2,*) ExpShell(nShell,k),DShell(nShell,k) + write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) + enddo + enddo + write(*,'(A28)') '------------------' + enddo + +! Total number of shells + + write(*,'(A28,1X,I16)') 'Number of shells in OBS',nShell + write(*,'(A28)') '------------------' + write(*,*) + +! Close file with basis set specification + + close(unit=2) + +!------------------------------------------------------------------------ +! Auxiliary basis set information +!------------------------------------------------------------------------ + +! Open file with auxilairy basis specification + + open(unit=3,file='input/auxbasis') + +! Read basis information + + write(*,'(A28)') 'Auxiliary basis set' + write(*,'(A28)') '-------------------' + + do i=1,NAtoms + read(3,*) iAt,nShAt + write(*,'(A28,1X,I16)') 'Atom n. ',iAt + write(*,'(A28,1X,I16)') 'number of shells ',nShAt + write(*,'(A28)') '------------------' + +! Basis function centers + + do j=1,nShAt + nShell = nShell + 1 + do k=1,3 + CenterShell(nShell,k) = XYZAtoms(iAt,k) + enddo + +! Shell type and contraction degree + + read(3,*) shelltype,KShell(nShell) + if(shelltype == "S") then + TotAngMomShell(nShell) = 0 + write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) + elseif(shelltype == "P") then + TotAngMomShell(nShell) = 1 + write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) + elseif(shelltype == "D") then + TotAngMomShell(nShell) = 2 + write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) + elseif(shelltype == "F") then + TotAngMomShell(nShell) = 3 + write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) + elseif(shelltype == "G") then + TotAngMomShell(nShell) = 4 + write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) + elseif(shelltype == "H") then + TotAngMomShell(nShell) = 5 + write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) + elseif(shelltype == "I") then + TotAngMomShell(nShell) = 6 + write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) + endif + +! Read exponents and contraction coefficients + + write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' + do k=1,Kshell(nShell) + read(3,*) ExpShell(nShell,k),DShell(nShell,k) + write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) + enddo + enddo + write(*,'(A28)') '------------------' + enddo + +! Total number of shells + + write(*,'(A28,1X,I16)') 'Number of shells in ABS',nShell + write(*,'(A28)') '------------------' + write(*,*) + +! Close file with basis set specification + + close(unit=3) + +end subroutine read_auxiliary_basis diff --git a/src/IntPak/read_basis.f90 b/src/IntPak/read_basis.f90 new file mode 100644 index 0000000..075e1d6 --- /dev/null +++ b/src/IntPak/read_basis.f90 @@ -0,0 +1,118 @@ +subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) + +! Read basis set information + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nNuc,nO + double precision,intent(in) :: rNuc(nNuc,ncart) + +! Local variables + + integer :: nShAt,iNuc,iShell + integer :: i,j,k + character :: shelltype + +! Output variables + + integer,intent(out) :: nShell,nBas + double precision,intent(out) :: CenterShell(maxShell,ncart) + integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell) + double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) + integer,intent(out) :: nV + +!------------------------------------------------------------------------ +! Primary basis set information +!------------------------------------------------------------------------ + +! Open file with basis set specification + + open(unit=2,file='input/basis') + +! Read basis information + + write(*,'(A28)') 'Gaussian basis set' + write(*,'(A28)') '------------------' + + nShell = 0 + do i=1,nNuc + read(2,*) iNuc,nShAt + write(*,'(A28,1X,I16)') 'Atom n. ',iNuc + write(*,'(A28,1X,I16)') 'number of shells ',nShAt + write(*,'(A28)') '------------------' + +! Basis function centers + + do j=1,nShAt + nShell = nShell + 1 + do k=1,ncart + CenterShell(nShell,k) = rNuc(iNuc,k) + enddo + +! Shell type and contraction degree + + read(2,*) shelltype,KShell(nShell) + if(shelltype == "S") then + TotAngMomShell(nShell) = 0 + write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) + elseif(shelltype == "P") then + TotAngMomShell(nShell) = 1 + write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) + elseif(shelltype == "D") then + TotAngMomShell(nShell) = 2 + write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) + elseif(shelltype == "F") then + TotAngMomShell(nShell) = 3 + write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) + elseif(shelltype == "G") then + TotAngMomShell(nShell) = 4 + write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) + elseif(shelltype == "H") then + TotAngMomShell(nShell) = 5 + write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) + elseif(shelltype == "I") then + TotAngMomShell(nShell) = 6 + write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) + endif + +! Read exponents and contraction coefficients + + write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' + do k=1,Kshell(nShell) + read(2,*) ExpShell(nShell,k),DShell(nShell,k) + write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) + enddo + enddo + write(*,'(A28)') '------------------' + enddo + +! Total number of shells + + write(*,'(A28,1X,I16)') 'Number of shells',nShell + write(*,'(A28)') '------------------' + write(*,*) + +! Close file with basis set specification + + close(unit=2) + +! Calculate number of basis functions + + nBas = 0 + do iShell=1,nShell + nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 + enddo + + write(*,'(A28)') '------------------' + write(*,'(A28,1X,I16)') 'Number of basis functions',NBas + write(*,'(A28)') '------------------' + write(*,*) + +! Number of virtual orbitals + + nV = nBas - nO + +end subroutine read_basis diff --git a/src/IntPak/read_molecule.f90 b/src/IntPak/read_molecule.f90 new file mode 100644 index 0000000..4d0ef7c --- /dev/null +++ b/src/IntPak/read_molecule.f90 @@ -0,0 +1,39 @@ +subroutine read_molecule(nAt,nEl,nO,nC,nR) + +! Read number of atoms nAt and number of electrons nEl + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(out) :: nAt,nEl,nO,nC,nR + +! Open file with geometry specification + + open(unit=1,file='input/molecule') + +! Read number of atoms and number of electrons + + read(1,*) + read(1,*) nAt,nEl,nC,nR + + nO = nEl/2 + +! Print results + + write(*,'(A28)') '----------------------' + write(*,'(A28,1X,I16)') 'Number of atoms',nAt + write(*,'(A28)') '----------------------' + write(*,*) + write(*,'(A28)') '----------------------' + write(*,'(A28,1X,I16)') 'Number of electrons',nEl + write(*,'(A28)') '----------------------' + write(*,*) + +! Close file with geometry specification + + close(unit=1) + +end subroutine read_molecule diff --git a/src/IntPak/read_options.f90 b/src/IntPak/read_options.f90 new file mode 100644 index 0000000..d4bbfec --- /dev/null +++ b/src/IntPak/read_options.f90 @@ -0,0 +1,122 @@ +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/MCQC/MCQC.f90 b/src/MCQC/MCQC.f90 index fb9315f..96c00d5 100644 --- a/src/MCQC/MCQC.f90 +++ b/src/MCQC/MCQC.f90 @@ -108,7 +108,7 @@ program MCQC ! nBas = number of basis functions (see below) ! = nO + nV - call read_molecule(nNuc,nEl,nC,nO,nR) + call read_molecule(nNuc,nEl,nO,nC,nR) allocate(ZNuc(nNuc),rNuc(nNuc,3)) ! Read geometry @@ -122,8 +122,7 @@ program MCQC ! Read basis set information !------------------------------------------------------------------------ - call read_basis(nNuc,rNuc,nBas,nC,nO,nV,nR,nS, & - nShell,TotAngMomShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell) + call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) !------------------------------------------------------------------------ ! Read auxiliary basis set information @@ -243,7 +242,7 @@ program MCQC allocate(F12(nBas,nBas,nBas,nBas),Yuk(nBas,nBas,nBas,nBas),FC(nBas,nBas,nBas,nBas,nBas,nBas)) ! Read integrals call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC) - call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,EcMP2(1),cHF,cHF,EcMP2F12) + call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,EcMP2(1),cHF,EcMP2F12) call cpu_time(end_MP2F12) t_MP2F12 = end_MP2F12 - start_MP2F12 diff --git a/src/MCQC/read_F12_integrals.f90 b/src/MCQC/read_F12_integrals.f90 index 0846be2..93fe2e5 100644 --- a/src/MCQC/read_F12_integrals.f90 +++ b/src/MCQC/read_F12_integrals.f90 @@ -13,7 +13,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y,FC) logical :: debug integer :: mu,nu,la,si,ka,ta - double precision :: ERI,F12,Yuk,ExpS + double precision :: ERI,F12,Yuk,F13C12,ExpS ! Output variables @@ -107,13 +107,15 @@ subroutine read_F12_integrals(nBas,S,C,F,Y,FC) FC = 0d0 do - read(31,*,end=31) mu,nu,la,si,ka,ta,FC + read(31,*,end=31) mu,nu,la,si,ka,ta,F13C12 + FC(mu,nu,la,si,ka,ta) = F13C12 enddo 31 close(unit=31) - ! Print results + if(debug) then + write(*,'(A28)') '----------------------' write(*,'(A28)') 'Electron repulsion integrals' write(*,'(A28)') '----------------------' @@ -123,6 +125,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y,FC) enddo enddo write(*,*) + write(*,'(A28)') '----------------------' write(*,'(A28)') 'F12 integrals' write(*,'(A28)') '----------------------' @@ -132,6 +135,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y,FC) enddo enddo write(*,*) + write(*,'(A28)') '----------------------' write(*,'(A28)') 'Yukawa integrals' write(*,'(A28)') '----------------------' @@ -141,6 +145,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y,FC) enddo enddo write(*,*) + endif ! Read exponent of Slater geminal