PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 04:13:52 +01:00
Code Releases Activity

MCQC under serious work

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-03-13 16:15:53 +01:00
parent 355e5026d9
commit 3059085d57
16 changed files with 673 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
input/basis
View File

@ -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

1
input/geminal Normal file
View File

@ -0,0 +1 @@
1.0

6
input/int
View File

@ -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

4
input/methods
View File

@ -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

4
input/molecule
View File

@ -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

4
input/options
View File

@ -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

29
input/weight
View File

@ -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

8
src/IntPak/Compute3eInt.f90
View File

@ -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

13
src/IntPak/IntPak.f90
View File

@ -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))

170
src/IntPak/elements.f90 Normal file
View File

@ -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

176
src/IntPak/read_auxiliary_basis.f90 Normal file
View File

@ -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

118
src/IntPak/read_basis.f90 Normal file
View File

@ -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

39
src/IntPak/read_molecule.f90 Normal file
View File

@ -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

122
src/IntPak/read_options.f90 Normal file
View File

@ -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

7
src/MCQC/MCQC.f90
View File

@ -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

11
src/MCQC/read_F12_integrals.f90
View File

@ -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