Skip to content

Commit

Permalink
Solving embedding calculation issue (#105).
Browse files Browse the repository at this point in the history
Added tests for Erkale localization, reading fchk files (for issue #103) and the WF-in-DFT method I implemented in my PhD thesis.
  • Loading branch information
fsmoncadaa committed Dec 13, 2024
1 parent 5568d3a commit 04434ca
Show file tree
Hide file tree
Showing 18 changed files with 4,273 additions and 38 deletions.
6 changes: 4 additions & 2 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,12 @@ install:: bin/lowdin bin/lowdin.x
$(SED) -i 's|COMPILATION_DATE|$(shell date)|g' $(TOPDIR)/lowdinvars.sh
cp -rf $(TOPDIR)/lowdinvars.sh $(PREFIX)/.lowdin2/
cp -rf lib $(PREFIX)/.lowdin2/
cp -rf utilities/erkale/build/erkale/basis $(PREFIX)/.lowdin2/lib/erkaleBasis
if [ -e utilities/erkale/build/erkale/basis ]; then \
cp -rf utilities/erkale/build/erkale/basis $(PREFIX)/.lowdin2/lib/erkaleBasis ; fi
mkdir -p $(PREFIX)/.lowdin2/bin
cp -rf $(TOPDIR)/bin/*.x $(PREFIX)/.lowdin2/bin
cp -rf utilities/erkale/erkale/bin/erkale_fchkpt utilities/erkale/erkale/bin/erkale_loc $(PREFIX)/.lowdin2/bin
if [ -e utilities/erkale/erkale/bin/erkale_loc ]; then \
cp -rf utilities/erkale/erkale/bin/erkale_fchkpt utilities/erkale/erkale/bin/erkale_loc $(PREFIX)/.lowdin2/bin ; fi
cp -rf $(TOPDIR)/bin/lowdin $(TOPDIR)
$(SED) -i 's|PREFIX|$(PREFIX)|g' $(TOPDIR)/lowdin
cp -rf $(TOPDIR)/lowdin $(PREFIX)/lowdin2
Expand Down
5 changes: 3 additions & 2 deletions bin/lowdinvars.sh
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,10 @@ then
export LOWDIN_DATA
fi

if [ -z "$LOWDIN_DATA/erkaleBasis" ]
if [ -z "$ERKALE_LIBRARY" ]
then
export ERKALE_LIBRARY="$LOWDIN_DATA/erkaleBasis"
ERKALE_LIBRARY="$LOWDIN_DATA/erkaleBasis"
export ERKALE_LIBRARY
fi

if [ -z "$PATH" ]
Expand Down
13 changes: 6 additions & 7 deletions src/core/MolecularSystem.f90
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,7 @@ subroutine MolecularSystem_showInformation(this)
print *,""
write (6,"(T5,A16,A)") "DESCRIPTION : ", trim( system%description )
write (6,"(T5,A16,I3)") "CHARGE : ",system%charge
write (6,"(T5,A16,F12.4)") "MASS (m_e) : ", MolecularSystem_getTotalMass(system)
write (6,"(T5,A16,ES10.4)") "MASS (m_e) : ", MolecularSystem_getTotalMass(system)
write (6,"(T5,A16,A4)") "PUNTUAL GROUP : ", "NONE"
print *,""

Expand Down Expand Up @@ -821,9 +821,8 @@ subroutine MolecularSystem_loadFromFile( form, targetPrefix )
do j = 1, size(MolecularSystem_instance%species(i)%particles)

read(40,*) MolecularSystem_instance%species(i)%particles(j)%nickname

MolecularSystem_instance%numberOfQuantumParticles = MolecularSystem_instance%numberOfQuantumParticles + 1
call BasisSet_load(MolecularSystem_instance%species(i)%particles(j)%basis, filename, unit = 40)
call BasisSet_load(MolecularSystem_instance%species(i)%particles(j)%basis, "LOWDIN.BAS", "SUBS. A REDUCED", MolecularSystem_instance%species(i)%name, unit = 40)

end do

Expand Down Expand Up @@ -2018,13 +2017,13 @@ subroutine MolecularSystem_readFchk( fileName, coefficients, densityMatrix, name

open(unit=fchkUnit, file=filename, status="old", form="formatted", access='sequential', action='read')

print *, ""
print *, "reading FCHK orbitals from ", fileName
! print *, ""
! print *, "reading FCHK orbitals from ", fileName
!The first two lines don't matter
read(fchkUnit,"(A100)",iostat=io) info
print *, info
! print *, info
read(fchkUnit,"(A100)",iostat=io) info
print *, info
! print *, info

!Read line by line, the first 40 characters determine what's being read
do
Expand Down
2 changes: 1 addition & 1 deletion src/scf/DensityMatrixSCFGuess.f90
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,6 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio
arguments(2) = nameOfSpecies
arguments(1) = "COEFFICIENTS"

call Matrix_constructor(densityMatrix, int(orderOfMatrix,8), int(orderOfMatrix,8), 0.0_8 )
call Matrix_constructor(orbitals, int(orderOfMatrix,8), int(orderOfMatrix,8), 0.0_8 )

readSuccess=.false.
Expand Down Expand Up @@ -160,6 +159,7 @@ subroutine DensityMatrixSCFGuess_getGuess( speciesID, hcoreMatrix, transformatio
end do
end if

call Matrix_constructor(densityMatrix, int(orderOfMatrix,8), int(orderOfMatrix,8), 0.0_8 )
do i = 1 , orderOfMatrix
do j = 1 , orderOfMatrix
do k = 1 , occupationNumber
Expand Down
14 changes: 4 additions & 10 deletions src/scf/MultiSCF.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1327,7 +1327,6 @@ subroutine MultiSCF_saveWfn(this,wfObjects)
character(30) :: labels(2)

character(50) :: integralsFile
character(50) :: arguments(20)
integer :: integralsUnit

!! Open file for wfn
Expand Down Expand Up @@ -1393,25 +1392,20 @@ subroutine MultiSCF_saveWfn(this,wfObjects)
call Matrix_writeToFile(wfObjects(speciesID)%fockMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "OVERLAP"
call Matrix_writeToFile(wfObjects(speciesID)%overlapMatrix, unit=wfnUnit, binary=.true., arguments = arguments(1:2) )
call Matrix_writeToFile(wfObjects(speciesID)%overlapMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "TRANSFORMATION"
call Matrix_writeToFile(wfObjects(speciesID)%transformationMatrix, unit=wfnUnit, binary=.true., arguments = arguments(1:2) )

if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) then
labels(1) = "EXTERNAL_POTENTIAL"
call Matrix_writeToFile(wfObjects(speciesID)%externalPotentialMatrix, unit=wfnUnit, binary=.true., arguments = labels )
end if
call Matrix_writeToFile(wfObjects(speciesID)%transformationMatrix, unit=wfnUnit, binary=.true., arguments = labels )

if (CONTROL_instance%COSMO) then
labels(1) = "COSMO1"
call Matrix_writeToFile(wfObjects(speciesID)%cosmo1, unit=wfnUnit, binary=.true., arguments = arguments(1:2) )
call Matrix_writeToFile(wfObjects(speciesID)%cosmo1, unit=wfnUnit, binary=.true., arguments = labels )
labels(1) = "COSMO2"
call Matrix_writeToFile(wfObjects(speciesID)%cosmo2, unit=wfnUnit, binary=.true., arguments = labels )
labels(1) = "COSMOCOUPLING"
call Matrix_writeToFile(wfObjects(speciesID)%cosmoCoupling, unit=wfnUnit, binary=.true., arguments = labels )
labels(1) = "COSMO4"
call Matrix_writeToFile(wfObjects(speciesID)%cosmo4, unit=wfnUnit, binary=.true., arguments = arguments(1:2) )
call Matrix_writeToFile(wfObjects(speciesID)%cosmo4, unit=wfnUnit, binary=.true., arguments = labels )
end if

end do
Expand Down
21 changes: 18 additions & 3 deletions src/scf/OrbitalLocalizer.f90
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit

close(30)

call system("erkale_fchkpt erkale.read")
call system("erkale_fchkpt erkale.read >> erkale_fchkpt.log")

!! Localize orbitals
open(unit=30, file="erkale.local", status="replace", form="formatted")
Expand All @@ -155,7 +155,7 @@ subroutine OrbitalLocalizer_erkaleLocal(speciesID,densityMatrix,fockMatrix,orbit

close(30)

call system("erkale_fchkpt erkale.write")
call system("erkale_fchkpt erkale.write >> erkale_fchkpt.log")

!! Read orbital coefficients from fchk files
call MolecularSystem_readFchk(trim(CONTROL_instance%INPUT_FILE)//trim(symbolOfSpecies)//".local.fchk", orbitalCoefficients, densityMatrix, nameOfSpecies )
Expand Down Expand Up @@ -1829,11 +1829,12 @@ subroutine OrbitalLocalizer_levelShiftSubsystemOrbitals()
call MolecularSystem_saveToFile()

!! Recalculate one particle integrals
call system("rm lowdin.opints")
! call system("rm lowdin.opints")
call system("lowdin-ints.x ONE_PARTICLE")

!! Start the wavefunction object
deallocate(WaveFunction_instance)
allocate(WaveFunction_instance(numberOfSpecies))
call WaveFunction_constructor(WaveFunction_instance,numberOfSpecies)

do speciesID=1, numberOfSpecies
Expand Down Expand Up @@ -2127,6 +2128,16 @@ subroutine OrbitalLocalizer_levelShiftSubsystemOrbitals()
labels(1) = "DENSITY"
call Matrix_writeToFile(wavefunction_instance(speciesID)%densityMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "KINETIC"
call Matrix_writeToFile(wavefunction_instance(speciesID)%kineticMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "ATTRACTION"
call Matrix_writeToFile(wavefunction_instance(speciesID)%puntualInteractionMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "EXTERNAL"
if(CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) &
call Matrix_writeToFile(wavefunction_instance(speciesID)%externalPotentialMatrix, unit=wfnUnit, binary=.true., arguments = labels )

labels(1) = "HCORE"
call Matrix_writeToFile(wavefunction_instance(speciesID)%hcoreMatrix, unit=wfnUnit, binary=.true., arguments = labels )

Expand All @@ -2148,10 +2159,14 @@ subroutine OrbitalLocalizer_levelShiftSubsystemOrbitals()
end if

if (CONTROL_instance%COSMO) then
labels(1) = "COSMO1"
call Matrix_writeToFile(WaveFunction_instance(speciesID)%cosmo1, unit=wfnUnit, binary=.true., arguments = labels(1:2) )
labels(1) = "COSMO2"
call Matrix_writeToFile(WaveFunction_instance(speciesID)%cosmo2, unit=wfnUnit, binary=.true., arguments = labels )
labels(1) = "COSMOCOUPLING"
call Matrix_writeToFile(WaveFunction_instance(speciesID)%cosmoCoupling, unit=wfnUnit, binary=.true., arguments = labels )
labels(1) = "COSMO4"
call Matrix_writeToFile(WaveFunction_instance(speciesID)%cosmo4, unit=wfnUnit, binary=.true., arguments = labels(1:2) )
end if

end do
Expand Down
46 changes: 46 additions & 0 deletions test/C2OPs-Embed-in-KS-P2-MP2-CISD.lowdin
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
SYSTEM_DESCRIPTION='C2OPs-TZ.geom'
GEOMETRY
e-[C] STO-3G 1.35434340389923 0.55988211198946 0.05525260812571
e-[H] STO-3G 2.45225753052674 -0.98466341806837 1.07188979363580
e-[C] STO-3G 2.63375141553823 -0.31899480577184 0.15987177742197
e-[H] STO-3G 2.54816069399926 -1.05872534257274 -0.70810872288146
e-[O] AUG-CC-PVDZ 3.78583048843762 0.37815496854981 0.19288336791534 addParticles=1 fragmentNumber=1
e-[H] STO-3G 1.39075957854909 1.17998258883901 -0.86139956877665
e-[H] STO-3G 1.29803871966429 1.25456719019296 0.91570797583479
e-[H] STO-3G 0.41993816938556 -0.04623329315828 0.03157276872451
E+ PSX-DZ 3.78583048843762 0.37815496854981 0.19288336791534 fragmentNumber=1
C dirac 1.35434340389923 0.55988211198946 0.05525260812571
H dirac 2.45225753052674 -0.98466341806837 1.07188979363580
C dirac 2.63375141553823 -0.31899480577184 0.15987177742197
H dirac 2.54816069399926 -1.05872534257274 -0.70810872288146
O dirac 3.78583048843762 0.37815496854981 0.19288336791534
H dirac 1.39075957854909 1.17998258883901 -0.86139956877665
H dirac 1.29803871966429 1.25456719019296 0.91570797583479
H dirac 0.41993816938556 -0.04623329315828 0.03157276872451
END GEOMETRY

TASKS
propagatorTheoryCorrection = 2
mollerPlessetCorrection = 2
configurationInteractionLevel = "CISD"
subsystemEmbedding=.T.
method = "UKS"
END TASKS

CONTROL
electronExchangeCorrelationFunctional="B3LYP"
subsystemBasisThreshold=0.001
subsystemOrbitalThreshold=0.2
erkaleLocalizationMethod="IAO"
forceClosedShell=.T.
readCoefficients=.F.
localizeOrbitals=.T.
totalEnergyTolerance=1E-10
scfGlobalMaxIterations=1000
integralsTransformationMethod = "C"
ionizeMO=1
ionizeSpecies = "POSITRON"
ptJustOneOrbital=T
mpFrozenCoreBoundary=1
CIdiagonalizationMethod = "JADAMILU"
END CONTROL
76 changes: 76 additions & 0 deletions test/C2OPs-Embed-in-KS-P2-MP2-CISD.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
#!/usr/bin/env python
from __future__ import print_function
import os
import sys
from colorstring import *

if len(sys.argv)==2:
lowdinbin = sys.argv[1]
else:
lowdinbin = "lowdin2"

testName = sys.argv[0][:-3]
inputName = testName + ".lowdin"
outputName = testName + ".out"

# Reference values and tolerance

refValues = {
"HF energy" : [-153.536301213424,1E-6],
"Embedded HF energy" : [-153.118946505217,1E-6],
"Embedded CISD" : [-153.310526831872,1E-6],
"Embedded MP2" : [-153.352174869549,1E-6],
"Embedded P2 H_1" : [-5.328000,1E-3]
}

testValues = dict(refValues) #copy
for value in testValues: #reset
testValues[value] = 0 #reset

# Run calculation

status = os.system(lowdinbin + " -i " + inputName)

if status:
print(testName + str_red(" ... NOT OK"))
sys.exit(1)

output = open(outputName, "r")
outputRead = output.readlines()

# Values
flag=1
for i in range(0,len(outputRead)):
line = outputRead[i]
if "TOTAL ENERGY =" in line and flag==1:
testValues["HF energy"] = float(line.split()[3])
flag=2
continue
if "TOTAL ENERGY =" in line and flag==2:
testValues["Embedded HF energy"] = float(line.split()[3])
flag=3
continue
if "STATE: 1 ENERGY =" in line:
testValues["Embedded CISD"] = float(line.split()[4])
if "E(MP2) =" in line:
testValues["Embedded MP2"] = float(line.split()[2])
if "Optimized second order pole:" in line:
testValues["Embedded P2 H_1"] = float(line.split()[4])

passTest = True

for value in refValues:
diffValue = abs(refValues[value][0] - testValues[value])
if ( diffValue <= refValues[value][1] ):
passTest = passTest * True
else :
passTest = passTest * False
print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue))

if passTest :
print(testName + str_green(" ... OK"))
else:
print(testName + str_red(" ... NOT OK"))
sys.exit(1)

output.close()
39 changes: 39 additions & 0 deletions test/Gly.e+-localize-molden.lowdin
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
SYSTEM_DESCRIPTION='glicina con postiron'

GEOMETRY
e-(N) 6-311G 0.92914 0.65159 0.25796
e-(C) 6-311G 0.08973 -0.44202 0.82116
e-(C) 6-311G 0.09416 -0.39610 2.35537
e-(O) 6-311G 0.79215 0.47431 2.87347
e-(O) 6-311G -0.60778 -1.24843 2.89495
e-(H) 6-311G 1.89448 0.53627 0.52036
e-(H) 6-311G 0.48806 -1.38246 0.47040
e-(H) 6-311G -0.91478 -0.32058 0.44545
e+ PSX-TZ -0.60778 -1.24843 2.89495
N dirac 0.92914 0.65159 0.25796
C dirac 0.08973 -0.44202 0.82116
C dirac 0.09416 -0.39610 2.35537
O dirac 0.79215 0.47431 2.87347
O dirac -0.60778 -1.24843 2.89495
H dirac 1.89448 0.53627 0.52036
H dirac 0.48806 -1.38246 0.47040
H dirac -0.91478 -0.32058 0.44545
END GEOMETRY

TASKS
method = "RHF"
END TASKS

CONTROL
readCoefficients=F
localizeOrbitals=.T.
erkaleLocalizationMethod="MU"
END CONTROL

OUTPUTS
moldenFile
END OUTPUTS




Loading

0 comments on commit 04434ca

Please sign in to comment.