Skip to content

Commit

Permalink
Changes to perform SCF core ionized states with the goal of doing XES:
Browse files Browse the repository at this point in the history 
Ionize selected orbitals with fractional occupation zero - to perform
Corrected "ionizeSpecies" spelling in input
Separated single SCF iterate into smaller subroutines
Updated exchange orbitals convergence method - required to ensure ionization of the desied state
Added core ionized H2O molecule example

On going: NOCI Franck Condon factors and transition dipole integrals
Matrix power and other functions computed from SVD
Separated NOCI generate densities into smaller subroutines
Small changes to print/write/read different molecular systems
Integer Vector read/write to file

Other:
One body energy components for CI density
Fixed output printing densities below 1E-100
  • Loading branch information
@fsmoncadaa
fsmoncadaa committed Feb 9, 2024
1 parent 83217e3 commit d3e8e76
Show file tree
Hide file tree
Showing 54 changed files with 2,915 additions and 1,051 deletions.
6 changes: 4 additions & 2 deletions bin/lowdin
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ if [ $extFile="lowdin" ]; then
'{
printf("\n&Output\n")
for(i=1;i<=NF;i++){
if($i=="specie"){
if($i=="species"){
printf("\tOutput_species = %s\n",$toupper((i+2)) )
}
else if($i=="orbital"){
Expand Down Expand Up @@ -361,7 +361,8 @@ if [ $extFile="lowdin" ]; then
cp $nameFile*.dens $LOWDIN_SCRATCH/$nameFile &> /dev/null
cp $nameFile*.sup $LOWDIN_SCRATCH/$nameFile &> /dev/null
cp $nameFile*.der $LOWDIN_SCRATCH/$nameFile &> /dev/null
cp $nameFile*.NOCI.coords $LOWDIN_SCRATCH/$nameFile &> /dev/null
cp $nameFile*.NOCI* $LOWDIN_SCRATCH/$nameFile &> /dev/null
cp $nameFile*.refNOCI* $LOWDIN_SCRATCH/$nameFile &> /dev/null
mv $nameFile*.dens $LOWDIN_SCRATCH/$nameFile &> /dev/null
mv $nameFile*.coup $LOWDIN_SCRATCH/$nameFile &> /dev/null
mv $nameFile*.over $LOWDIN_SCRATCH/$nameFile &> /dev/null
Expand Down Expand Up @@ -444,6 +445,7 @@ if [ $extFile="lowdin" ]; then
mv $LOWDIN_SCRATCH/$nameFile/*.plainvec $currentPath &> 2
# mv $LOWDIN_SCRATCH/$nameFile/*.val $currentPath &> 2
mv $LOWDIN_SCRATCH/$nameFile/*.NOCI.coords $currentPath &> 2
mv $LOWDIN_SCRATCH/$nameFile/*.NOCI.s* $currentPath &> 2
mv $LOWDIN_SCRATCH/$nameFile/*.cub $currentPath &> 2
mv $LOWDIN_SCRATCH/$nameFile/*.dens $currentPath &> 2
mv $LOWDIN_SCRATCH/$nameFile/*.eps $currentPath &> 2
Expand Down
56 changes: 56 additions & 0 deletions examples/FCI_Charry2018/H-.FCI-DZ.lowdin
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
!The goal of this calculation is to compute the binding energy of a positron bound complex
!Reported:
!E(H-): -0.524029

SYSTEM_DESCRIPTION='H- from Charry 2018 (10.1002/anie.201800914)'

GEOMETRY
e-(H) AUG-CC-PVDZ 0.00 0.00 0.00 addParticles=1
H dirac 0.00 0.00 0.00
END GEOMETRY

!method to solve the SCF - CI only works for unrestricted reference
!CI level strings chooses the desired excitations to be included. FCI is all possible excitations

TASKS
method = "UHF"
configurationInteractionLevel ="FCI"
!configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ"
END TASKS

!Compute only the "numberOfCIstates" states. Here we only need the ground state
!Compute the density matrix for "CIstatesToPrint" states, for density outputs
!Generate the natural orbitals, for visualization in molden files
!The Davidson diagonalization implemented in JADAMILU is the recomended method.
!For small systems, full matrix diagonalization with DSYEVX is possible
!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed
!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing
!Strict SCF convergence improves the quality of the CI results (not required for the FCI)

CONTROL
numberOfCIstates=1
CIStatesToPrint=1
CINaturalOrbitals=T
CIdiagonalizationMethod = "JADAMILU"
!CIdiagonalizationMethod = "DSYEVX"
!CIdiagonalizationMethod = "ARPACK"
CIPrintEigenVectorsFormat = "OCCUPIED" !"NONE","ORBITALS"
CIPrintThreshold = 5e-2
!totalEnergyTolerance=1E-12
END CONTROL

!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space
INPUT_CI
species="E-ALPHA" core=0 active=0
species="E-BETA" core=0 active=0
END INPUT_CI

!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results
OUTPUTS
moldenFile state=1
densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1
END OUTPUTS




58 changes: 58 additions & 0 deletions examples/FCI_Charry2018/PsH.FCI-DZ.lowdin
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
!The goal of this calculation is to compute the binding energy of a positron bound complex
!Reported:
!E(PsH): -0.734559

SYSTEM_DESCRIPTION='PsH from Charry 2018 (10.1002/anie.201800914)'

GEOMETRY
e-(H) AUG-CC-PVDZ 0.00 0.00 0.00 addParticles=1
e+ E+-H-AUG-CC-PVDZ 0.00 0.00 0.00
H dirac 0.00 0.00 0.00
END GEOMETRY

!method to solve the SCF - CI only works for unrestricted reference
!CI level strings chooses the desired excitations to be included. FCI is all possible excitations

TASKS
method = "UHF"
configurationInteractionLevel ="FCI"
!configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ"
END TASKS

!Compute only the "numberOfCIstates" states. Here we only need the ground state
!Compute the density matrix for "CIstatesToPrint" states, for density outputs
!Generate the natural orbitals, for visualization in molden files
!The Davidson diagonalization implemented in JADAMILU is the recomended method.
!For small systems, full matrix diagonalization with DSYEVX is possible
!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed
!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing
!Strict SCF convergence improves the quality of the CI results (not required for the FCI)

CONTROL
numberOfCIstates=1
CIStatesToPrint=1
CINaturalOrbitals=T
CIdiagonalizationMethod = "JADAMILU"
!CIdiagonalizationMethod = "DSYEVX"
!CIdiagonalizationMethod = "ARPACK"
CIPrintEigenVectorsFormat = "OCCUPIED" !"NONE","ORBITALS"
CIPrintThreshold = 5e-2
!totalEnergyTolerance=1E-12
END CONTROL

!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space
INPUT_CI
species="E-ALPHA" core=0 active=0
species="E-BETA" core=0 active=0
species="POSITRON" core=0 active=0
END INPUT_CI

!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results
OUTPUTS
moldenFile state=1
densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1
END OUTPUTS




64 changes: 64 additions & 0 deletions examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
!The goal of this calculation is to compute the binding energy of a positron bound complex
!Reported:
!E(e+H2^2-): -1.279680 a.u.

SYSTEM_DESCRIPTION='PsH from Charry 2018 (10.1002/anie.201800914)'

!add two electrons (one for each hydrogen anion)
!remove one positron
GEOMETRY
e-(H) AUG-CC-PVDZ 0.00 0.00 -1.6 addParticles=1
e-(H) AUG-CC-PVDZ 0.00 0.00 1.6 addParticles=1
e+ E+-H-AUG-CC-PVDZ 0.00 0.00 -1.6
e+ E+-H-AUG-CC-PVDZ 0.00 0.00 1.6 addParticles=-1
H dirac 0.00 0.00 -1.6
H dirac 0.00 0.00 1.6
END GEOMETRY

!method to solve the SCF - CI only works for unrestricted reference
!CI level strings chooses the desired excitations to be included. FCI is all possible excitations

TASKS
method = "UHF"
configurationInteractionLevel ="FCI"
!configurationInteractionLevel ="CIS","CISD","CISDT","CISDTQ"
END TASKS

!Compute only the "numberOfCIstates" states. Here we only need the ground state
!Compute the density matrix for "CIstatesToPrint" states, for density outputs
!Generate the natural orbitals, for visualization in molden files
!The Davidson diagonalization implemented in JADAMILU is the recomended method.
!For small systems, full matrix diagonalization with DSYEVX is possible
!CI EigenVectors with coefficient higher than "CIPrintThreshold" are printed
!Printing format "OCCUPIED" shows the coefficients, "ORBITALS" shows the strings, "NONE" skips printing
!Strict SCF convergence improves the quality of the CI results (not required for the FCI)

CONTROL
numberOfCIstates=1
CIStatesToPrint=1
CINaturalOrbitals=T
CIdiagonalizationMethod = "JADAMILU"
!CIdiagonalizationMethod = "DSYEVX"
CIPrintEigenVectorsFormat = "OCCUPIED"
!CIPrintEigenVectorsFormat = "NONE","ORBITALS"
CIPrintThreshold = 5e-2
!totalEnergyTolerance=1E-12
END CONTROL

!INPUT_CI block help us define the frozen core and active virtuals orbitals. Here we are not restricting the excitation space
INPUT_CI
species="E-ALPHA" core=0 active=0
species="E-BETA" core=0 active=0
species="POSITRON" core=0 active=0
END INPUT_CI

!With CI, moldenFiles, 1D and 2D density slices and density cubes are good ways to visualize the density results
OUTPUTS
moldenFile state=1
densityPlot dimensions=2 point1=0.0 0.0 -6.0 point2=0.0 0.0 6.0 state=1
densityPlot dimensions=3 point1=0.0 -6.0 -6.0 point2=0.0 -6.0 6.0 point3=0.0 6.0 -6.0 state=1
END OUTPUTS




8 changes: 4 additions & 4 deletions lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@

polarizabilities = [(angstromToBohr**3)*x for x in polarizabilities] # a.u.
maxR = int(maxR/step)
print maxR
print(maxR)

angularMoment = list()
for i in range(0,len(exponents)) :
Expand All @@ -45,14 +45,14 @@

ii = 0
for atom in range(0,len(origin)):
print origin[atom]
print(origin[atom])
realPotentialFileName = potentialName + "."+ str(origin[atom][2]) +".data"

alpha = polarizabilities[atom]
r = origin[atom][2]
realPotentialFile = open (realPotentialFileName, "w")

for i in xrange(1,maxR+1,1):
for i in range(1,maxR+1,1):
i = i*step
realPotentialFile.write(str(i) + " "+ str( -1.0*alpha/(2.0*(i-r)**4) * (1 - math.exp(-((i-r)**6)/(cutoff**6))) ) + "\n" )
realPotentialFile.close()
Expand Down Expand Up @@ -108,7 +108,7 @@
gnuplotOutput.close()

if len(coeff) == 0:
print "Fitting error, results not found"
print("Fitting error, results not found")
continue

coefficients = [float(k) for k in coeff]
Expand Down
80 changes: 58 additions & 22 deletions src/CI/ConfigurationInteraction.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1142,6 +1142,7 @@ subroutine ConfigurationInteraction_densityMatrices()
character(50) :: file, wfnfile, speciesName, auxstring
character(50) :: arguments(2)
type(matrix), allocatable :: coefficients(:), atomicDensityMatrix(:,:), ciDensityMatrix(:,:), auxDensMatrix(:,:)
type(matrix), allocatable :: kineticMatrix(:), attractionMatrix(:), externalPotMatrix(:)
integer numberOfSpecies

type(matrix) :: auxdensityEigenVectors
Expand Down Expand Up @@ -1213,8 +1214,13 @@ subroutine ConfigurationInteraction_densityMatrices()
write (*,*) "=============================="
write (*,*) ""

allocate( coefficients(numberOfSpecies), atomicDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), &
ciDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), auxDensMatrix(numberOfSpecies,ConfigurationInteraction_instance%nproc) )
allocate( coefficients(numberOfSpecies), &
kineticMatrix(numberOfSpecies), &
attractionMatrix(numberOfSpecies), &
externalPotMatrix(numberOfSpecies), &
atomicDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), &
ciDensityMatrix(numberOfSpecies,CONTROL_instance%CI_STATES_TO_PRINT), &
auxDensMatrix(numberOfSpecies,ConfigurationInteraction_instance%nproc) )

wfnFile = "lowdin.wfn"
wfnUnit = 20
Expand All @@ -1228,33 +1234,45 @@ subroutine ConfigurationInteraction_densityMatrices()
! numberOfOrbitals = ConfigurationInteraction_instance%numberOfOrbitals%values(species)
numberOfOccupiedOrbitals = ConfigurationInteraction_instance%numberOfOccupiedOrbitals%values(species)

arguments(2) = speciesName
arguments(1) = "COEFFICIENTS"
! print *, "trolo", numberOfOrbitals, numberOfContractions, numberOfOccupiedOrbitals
arguments(2) = speciesName
! print *, "trolo", numberOfOrbitals, numberOfContractions, numberOfOccupiedOrbitals

coefficients(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))
arguments(1) = "COEFFICIENTS"
coefficients(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))

! print *, "trololo"
arguments(1) = "KINETIC"
kineticMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))

arguments(1) = "ATTRACTION"
attractionMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))

arguments(1) = "EXTERNAL_POTENTIAL"
if( CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) &
externalPotMatrix(species) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))
! print *, "trololo"

do state=1, CONTROL_instance%CI_STATES_TO_PRINT
do state=1, CONTROL_instance%CI_STATES_TO_PRINT

call Matrix_constructor ( ciDensityMatrix(species,state) , &
int(numberOfContractions,8), &
int(numberOfContractions,8), 0.0_8 )
call Matrix_constructor ( ciDensityMatrix(species,state) , &
int(numberOfContractions,8), &
int(numberOfContractions,8), 0.0_8 )

do k=1, numberOfOccupiedOrbitals
ciDensityMatrix(species,state)%values( k, k)=1.0_8
end do
do k=1, numberOfOccupiedOrbitals
ciDensityMatrix(species,state)%values( k, k)=1.0_8
end do

end do
end do

do n=1, ConfigurationInteraction_instance%nproc
do n=1, ConfigurationInteraction_instance%nproc

call Matrix_constructor ( auxDensMatrix(species,n) , &
int(numberOfContractions,8), &
int(numberOfContractions,8), 0.0_8 )
end do
call Matrix_constructor ( auxDensMatrix(species,n) , &
int(numberOfContractions,8), &
int(numberOfContractions,8), 0.0_8 )
end do
end do

close(wfnUnit)
Expand Down Expand Up @@ -1566,7 +1584,25 @@ subroutine ConfigurationInteraction_densityMatrices()
end do
end do


write(*,*) ""
write(*,*) "==============================="
write(*,*) " ONE BODY ENERGY CONTRIBUTIONS:"
write(*,*) ""
do state=1, CONTROL_instance%CI_STATES_TO_PRINT
write(*,*) " STATE: ", state
do species=1, molecularSystem_instance%numberOfQuantumSpecies
write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name ) // &
" Kinetic energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*kineticMatrix(species)%values)
write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name ) // &
"/Fixed interact. energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*attractionMatrix(species)%values)
if( CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL) &
write(*,"(A38,F25.12)") trim( MolecularSystem_instance%species(species)%name) // &
" Ext Pot energy = ", sum(transpose(atomicDensityMatrix(species,state)%values)*externalPotMatrix(species)%values)
print *, ""
end do
print *, ""
end do

!! Natural orbitals

if (CONTROL_instance%CI_NATURAL_ORBITALS) then
Expand Down
Loading

0 comments on commit d3e8e76

Please sign in to comment.