Skip to content

Commit

Permalink
Implemented the rotational formula from NOCI elements, as in our rece…
Browse files Browse the repository at this point in the history 
…ntly submitted paper.

Improved NOCI in parallel, by running several SCFs at the same time.
This involved moving away from MolecularSystem_instance and Wavefunction_instance public structures,
and removing public structures from DFT program (Grid, Functionals).
Split NonOrthogonalCI module into smaller modules.
Added two RoCI tests, and one KS test for positron-electron correlation.
Remove duplicate "getNameOfSpecie" routine.
Fixed issues #79 #83
  • Loading branch information
@fsmoncadaa
fsmoncadaa committed Nov 19, 2024
1 parent d13e413 commit 6ca18d3
Show file tree
Hide file tree
Showing 56 changed files with 6,127 additions and 4,850 deletions.
19 changes: 9 additions & 10 deletions bin/lowdin
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,7 +358,7 @@ if [ $extFile="lowdin" ]; then
fi
gawk '($1~/BASIS/ && toupper($2)~/^'$BASIS_NAME'$/){flag=1; next}
($0~/END/){flag=0};
(flag==1){print toupper($0)}' $nameFile > $LOWDIN_DATA/basis/$BASIS_NAME
(flag==1){print toupper($0)}' $nameFile > $BASIS_NAME.$PID
done
fi

Expand Down Expand Up @@ -393,6 +393,14 @@ if [ $extFile="lowdin" ]; then
mv $nameFile*.over $LOWDIN_SCRATCH/$nameFile &> /dev/null
mv $nameFile*.kin $LOWDIN_SCRATCH/$nameFile &> /dev/null
mv $nameFile*.coeff $LOWDIN_SCRATCH/$nameFile &> /dev/null
#PID to avoid basis duplicates in simultaneous calculations
if [ ${#BASIS_NAMES[@]} -gt "0" ]
then
for BASIS_NAME in ${BASIS_NAMES[@]}
do
mv $BASIS_NAME.$PID $LOWDIN_SCRATCH/$nameFile/$BASIS_NAME &> /dev/null
done
fi

if [ -e $nameFile.gms.bs ]
then
Expand Down Expand Up @@ -492,15 +500,6 @@ if [ $extFile="lowdin" ]; then
rm -rf $LOWDIN_SCRATCH/$nameFile
fi

### Clean custom basis files
if [ ${#BASIS_NAMES[@]} -gt "0" ]
then
for BASIS_NAME in ${BASIS_NAMES[@]}
do
rm $LOWDIN_DATA/basis/$BASIS_NAME
done
fi

else
echo $1 ", this file does not exist. "
exit 1
Expand Down
14 changes: 14 additions & 0 deletions lib/basis/NAKAI-CC-PVDZ
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
O-HYDROGEN H (CC-PVDZ) BASIS TYPE: 1
#
5
1 0 1
13.01000000 1.0
2 0 1
1.96200000 1.0
3 0 1
0.44460000 1.0
4 0 1
0.12200000 1.00000000
5 1 1
0.72700000 1.00000000

2 changes: 1 addition & 1 deletion src/CI/CIStrings.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ subroutine CIStrings_buildStrings()

CIcore_instance%numberOfStrings(i)%values(1) = 1 !! ground

write (*,"(A,A)") " ", MolecularSystem_getNameOfSpecie(i)
write (*,"(A,A)") " ", MolecularSystem_getNameOfSpecies(i)

do cilevel = 1,CIcore_instance%CILevel(i)

Expand Down
2 changes: 1 addition & 1 deletion src/CI/CIcore.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ subroutine CIcore_constructor(level)


do i=1, numberOfSpecies
nameOfSpecie= trim( MolecularSystem_getNameOfSpecie( i ) )
nameOfSpecie= trim( MolecularSystem_getNameOfSpecies( i ) )
numberOfContractions = MolecularSystem_getTotalNumberOfContractions( i )

arguments(2) = nameOfSpecie
Expand Down
18 changes: 9 additions & 9 deletions src/CI/CImod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,9 +111,9 @@ subroutine CImod_run()

write (*,"(A32)",advance="no") "Number of orbitals for species: "
do i = 1, numberOfSpecies-1
write (*,"(A)",advance="no") trim(MolecularSystem_getNameOfSpecie(i))//", "
write (*,"(A)",advance="no") trim(MolecularSystem_getNameOfSpecies(i))//", "
end do
write (*,"(A)",advance="no") trim(MolecularSystem_getNameOfSpecie(numberOfSpecies))
write (*,"(A)",advance="no") trim(MolecularSystem_getNameOfSpecies(numberOfSpecies))
write (*,*) ""

write (*,"(A28)",advance="no") " occupied orbitals: "
Expand Down Expand Up @@ -558,7 +558,7 @@ subroutine CImod_getTransformedIntegrals()
allocate(CIcore_instance%fourIndexArray(numberOfSpecies))

do i=1, numberOfSpecies
nameOfSpecie= trim( MolecularSystem_getNameOfSpecie( i ) )
nameOfSpecie= trim( MolecularSystem_getNameOfSpecies( i ) )
specieID = MolecularSystem_getSpecieID( nameOfSpecie=nameOfSpecie )
ocupationNumber = MolecularSystem_getOcupationNumber( i )
numberOfContractions = MolecularSystem_getTotalNumberOfContractions( i )
Expand All @@ -577,7 +577,7 @@ subroutine CImod_getTransformedIntegrals()

open(unit=wfnUnit, file=trim(wfnFile), status="old", form="unformatted")

arguments(2) = MolecularSystem_getNameOfSpecie(i)
arguments(2) = MolecularSystem_getNameOfSpecies(i)
arguments(1) = "COEFFICIENTS"

coefficients = &
Expand Down Expand Up @@ -651,7 +651,7 @@ subroutine CImod_getTransformedIntegrals()
if ( numberOfSpecies > 1 ) then
do j = 1 , numberOfSpecies
if ( i .ne. j) then
nameOfOtherSpecie = trim( MolecularSystem_getNameOfSpecie( j ) )
nameOfOtherSpecie = trim( MolecularSystem_getNameOfSpecies( j ) )
otherSpecieID = MolecularSystem_getSpecieID( nameOfSpecie=nameOfOtherSpecie )
ocupationNumberOfOtherSpecie = MolecularSystem_getOcupationNumber( j )
numberOfContractionsOfOtherSpecie = MolecularSystem_getTotalNumberOfContractions( j )
Expand Down Expand Up @@ -1077,7 +1077,7 @@ subroutine CImod_densityMatrices()

!Inicializando las matrices
do species=1, numberOfSpecies
speciesName = MolecularSystem_getNameOfSpecie(species)
speciesName = MolecularSystem_getNameOfSpecies(species)

numberOfContractions = MolecularSystem_getTotalNumberOfContractions( species )
! numberOfOrbitals = CIcore_instance%numberOfOrbitals%values(species)
Expand Down Expand Up @@ -1388,7 +1388,7 @@ subroutine CImod_densityMatrices()

!! Building the CI reduced density matrix in the atomic orbital representation
do species=1, numberOfSpecies
speciesName = MolecularSystem_getNameOfSpecie(species)
speciesName = MolecularSystem_getNameOfSpecies(species)
numberOfContractions = MolecularSystem_getTotalNumberOfContractions( species )

do state=1, CONTROL_instance%CI_STATES_TO_PRINT
Expand Down Expand Up @@ -1472,7 +1472,7 @@ subroutine CImod_densityMatrices()
write(*,*) "-----------------"

numberOfContractions = MolecularSystem_getTotalNumberOfContractions( species )
speciesName = MolecularSystem_getNameOfSpecie(species)
speciesName = MolecularSystem_getNameOfSpecies(species)


call Vector_constructor ( auxdensityEigenValues, &
Expand Down Expand Up @@ -1684,7 +1684,7 @@ subroutine CImod_densityMatrices()
! end do

! !Write occupation numbers to file
! write (6,"(T8,A10,A20)") trim(MolecularSystem_getNameOfSpecie(specie)),"OCCUPATIONS:"
! write (6,"(T8,A10,A20)") trim(MolecularSystem_getNameOfSpecies(specie)),"OCCUPATIONS:"

! call Matrix_show ( ciOccupationNumbers )

Expand Down
2 changes: 1 addition & 1 deletion src/CI/Configuration.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1147,7 +1147,7 @@ subroutine Configuration_show(this)


do i=1, numberOfSpecies
print *, "For specie ", MolecularSystem_getNameOfSpecie ( i )
print *, "For specie ", MolecularSystem_getNameOfSpecies( i )
! print *, "Excitations: ", this%order(i)
! print *, "Ndeterminants: ",this%nDeterminants
! print *, "Occupations"
Expand Down
28 changes: 14 additions & 14 deletions src/CalcProp/CalculateProperties.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,10 +146,10 @@ subroutine CalculateProperties_constructor( this,fileName )
open(unit = occupationsUnit, file=trim(occupationsFile), status="old", form="formatted")
do speciesID=1, numberOfSpecies
numberOfContractions = MolecularSystem_getTotalNumberOfContractions (speciesID )
print *, "We are calculating properties for ", trim(MolecularSystem_getNameOfSpecie(speciesID)), &
print *, "We are calculating properties for ", trim(MolecularSystem_getNameOfSpecies(speciesID)), &
" in the CI ground state"
auxstring="1" !ground state
arguments(2) = MolecularSystem_getNameOfSpecie(speciesID)
arguments(2) = MolecularSystem_getNameOfSpecies(speciesID)
arguments(1) = "DENSITYMATRIX"//trim(adjustl(auxstring))
this%densityMatrix(speciesID)= Matrix_getFromFile(unit=occupationsUnit, rows= int(numberOfcontractions,4), &
columns= int(numberOfcontractions,4), binary=.false., arguments=arguments(1:2))
Expand All @@ -159,9 +159,9 @@ subroutine CalculateProperties_constructor( this,fileName )
open(unit=wfnUnit, file=trim(wfnFile), status="old", form="unformatted")
do speciesID=1, numberOfSpecies
numberOfContractions = MolecularSystem_getTotalNumberOfContractions (speciesID )
print *, "We are calculating properties for ", trim(MolecularSystem_getNameOfSpecie(speciesID)), &
print *, "We are calculating properties for ", trim(MolecularSystem_getNameOfSpecies(speciesID)), &
" in the HF/KS ground state"
arguments(2) = MolecularSystem_getNameOfSpecie(speciesID)
arguments(2) = MolecularSystem_getNameOfSpecies(speciesID)
arguments(1) = "DENSITY"
this%densityMatrix(speciesID) = Matrix_getFromFile(unit=wfnUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))
Expand All @@ -180,7 +180,7 @@ subroutine CalculateProperties_constructor( this,fileName )
do speciesID=1, numberOfSpecies
numberOfContractions = MolecularSystem_getTotalNumberOfContractions (speciesID )
! Overlap matrix
arguments(2) = MolecularSystem_getNameOfSpecie(speciesID)
arguments(2) = MolecularSystem_getNameOfSpecies(speciesID)
arguments(1) = "OVERLAP"
this%overlapMatrix(speciesID) = Matrix_getFromFile(unit=integralsUnit, rows= int(numberOfContractions,4), &
columns= int(numberOfContractions,4), binary=.true., arguments=arguments(1:2))
Expand Down Expand Up @@ -294,7 +294,7 @@ subroutine CalculateProperties_showPopulationAnalyses(this)

do type= 1, size(analysis)

speciesName = trim(MolecularSystem_getNameOfSpecie( speciesID ))
speciesName = trim(MolecularSystem_getNameOfSpecies( speciesID ))

if(trim(speciesName) .eq. "E-ALPHA") then
speciesNickname="E-"
Expand Down Expand Up @@ -366,7 +366,7 @@ subroutine CalculateProperties_showPopulationAnalyses(this)

! search_specie: do i = 1, MolecularSystem_getNumberOfQuantumSpecies()
! speciesName=""
! speciesName = trim(MolecularSystem_getNameOfSpecie(i))
! speciesName = trim(MolecularSystem_getNameOfSpecies(i))

! if( scan(trim(speciesName),"E")==1 ) then
! if( scan(trim(speciesName),"-")>1 ) then
Expand Down Expand Up @@ -408,7 +408,7 @@ function CalculateProperties_getPopulation( this, typeOfPopulation, speciesID, t

call Matrix_constructor( auxMatrix, int( numberOfcontractions, 8), int( numberOfcontractions, 8) )

speciesName=trim(MolecularSystem_getNameOfSpecie( speciesID ))
speciesName=trim(MolecularSystem_getNameOfSpecies( speciesID ))
if(trim(speciesName) .eq. "E-ALPHA") then
call Matrix_constructor( output, int( numberOfcontractions, 8), 2_8 )
otherSpeciesID=speciesID+1
Expand Down Expand Up @@ -485,7 +485,7 @@ subroutine CalculateProperties_showExpectedPositions(this)
print *,""
write (6,"(T19,4A9)") "<x>","<y>", "<z>", ""
do i=1, numberOfSpecies
write (6,"(T5,A15,3F9.4)") trim(MolecularSystem_getNameOfSpecie( i )), CalculateProperties_getExpectedPosition(this, i)
write (6,"(T5,A15,3F9.4)") trim(MolecularSystem_getNameOfSpecies( i )), CalculateProperties_getExpectedPosition(this, i)
end do
print *,""
print *,"END EXPECTED POSITIONS"
Expand Down Expand Up @@ -540,7 +540,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
do i=1, numberOfSpecies
dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecies(this, i)
totalDipole(:)=totalDipole(:)+dipole(i,:)
write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecie( i )), dipole(i,:)
write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecies( i )), dipole(i,:)
end do
dipole(numberOfSpecies+1,:)=CalculateProperties_getDipoleOfPuntualCharges()
totalDipole(:)=totalDipole(:)+dipole(numberOfSpecies+1,:)
Expand All @@ -558,7 +558,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
do i=1, numberOfSpecies
dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecies(this, i)*2.54174619
totalDipole(:)=totalDipole(:)+dipole(i,:)
write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecie( i )), dipole(i,:)
write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecies( i )), dipole(i,:)
end do

dipole(numberOfSpecies+1,:)=CalculateProperties_getDipoleOfPuntualCharges()*2.54174619
Expand All @@ -579,7 +579,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
do i=1, numberOfSpecies
quadrupole(i,:)=CalculateProperties_getQuadrupoleOfQuantumSpecies(this, i)*2.54174619*0.52917720859
totalQuadrupole(:)=totalQuadrupole(:)+quadrupole(i,:)
write (6,"(T5,A15,6F14.8)") trim(MolecularSystem_getNameOfSpecie( i )), quadrupole(i,:)
write (6,"(T5,A15,6F14.8)") trim(MolecularSystem_getNameOfSpecies( i )), quadrupole(i,:)
end do

quadrupole(numberOfSpecies+1,:)=CalculateProperties_getQuadrupoleOfPuntualCharges()*2.54174619*0.52917720859
Expand Down Expand Up @@ -755,7 +755,7 @@ end module CalculateProperties_
!
!
! do i=1, numberOfSpecies
! nameOfSpecieSelected = trim( Particle_Manager_getNameOfSpecie( i ) )
! nameOfSpecieSelected = trim( Particle_Manager_getNameOfSpecies( i ) )
! numberOfContractions = Particle_Manager_getTotalNumberOfContractions( i )
! call Matrix_constructor (densityMatrix, int(numberOfContractions,8), int(numberOfContractions,8))
! densityMatrix = MolecularSystem_getDensityMatrix( trim(nameOfSpecieSelected) )
Expand Down Expand Up @@ -796,7 +796,7 @@ end module CalculateProperties_
! print *,""
! write (6,"(T19,A9)") "<R^2>"
! do i=1, numberOfSpecies
! write (6,"(T5,A15,F9.4)") trim(Particle_Manager_getNameOfSpecie( i )), (this%expectedR2%values(i))
! write (6,"(T5,A15,F9.4)") trim(Particle_Manager_getNameOfSpecies( i )), (this%expectedR2%values(i))
! end do
! print *,""
! print *,"END EXPECTED <R^2>"
Expand Down
27 changes: 19 additions & 8 deletions src/DFT/DFT.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ program DFT
use MolecularSystem_
use DensityFunctionalTheory_
use GridManager_
use Functional_
use String_
use Matrix_
use Exception_
Expand All @@ -36,6 +37,7 @@ program DFT

character(50) :: job
character(100) :: densFile
type(Grid), allocatable :: grids(:), gridsCommonPoints(:,:)
type(Matrix), allocatable :: densityMatrix(:)
type(Matrix), allocatable :: exchangeCorrelationMatrix(:)
type(Matrix) :: exchangeCorrelationEnergy
Expand All @@ -62,21 +64,30 @@ program DFT
!!Load the system in lowdin.sys format
call MolecularSystem_loadFromFile( "LOWDIN.SYS" )

call Functional_createFunctionals( )
numberOfSpecies=MolecularSystem_getNumberOfQuantumSpecies()

!! Allocate memory.
if(allocated(grids)) deallocate(grids)
allocate(grids(numberOfSpecies))

if (allocated(gridsCommonPoints)) deallocate(gridsCommonPoints)
allocate(gridsCommonPoints(numberOfSpecies,numberOfSpecies))

do speciesID = 1 , numberOfSpecies
grids(speciesID)%molSys => MolecularSystem_instance
end do

!!!Building grids jobs
select case ( job )
case ("BUILD_SCF_GRID")
call DensityFunctionalTheory_buildSCFGrid()
call DensityFunctionalTheory_buildSCFGrid(grids,gridsCommonPoints)
STOP
case ("BUILD_FINAL_GRID" )
call DensityFunctionalTheory_buildFinalGrid()
call DensityFunctionalTheory_buildFinalGrid(grids,gridsCommonPoints)
STOP
end select

!!!Computing energy and potential jobs
numberOfSpecies=MolecularSystem_getNumberOfQuantumSpecies()

!!!Computing energy and potential jobs
allocate( densityMatrix(numberOfSpecies) , numberOfParticles(numberOfSpecies), &
exchangeCorrelationMatrix(numberOfSpecies))

Expand All @@ -99,7 +110,7 @@ program DFT

select case ( job )
case ("SCF_DFT")
call DensityFunctionalTheory_SCFDFT(densityMatrix, exchangeCorrelationMatrix, exchangeCorrelationEnergy, numberOfParticles)
call DensityFunctionalTheory_SCFDFT(grids,gridsCommonPoints,densityMatrix, exchangeCorrelationMatrix, exchangeCorrelationEnergy, numberOfParticles)
case ("FINAL_DFT")
!read scf information for comparison
do speciesID = 1 , numberOfSpecies
Expand All @@ -124,7 +135,7 @@ program DFT
close(unit=excUnit)
end do

call DensityFunctionalTheory_finalDFT(densityMatrix, exchangeCorrelationMatrix, exchangeCorrelationEnergy, numberOfParticles)
call DensityFunctionalTheory_finalDFT(grids,gridsCommonPoints,densityMatrix, exchangeCorrelationMatrix, exchangeCorrelationEnergy, numberOfParticles)
! case default
! write(*,*) "USAGE: lowdin-DFT.x job "
! write(*,*) "Where job can be: "
Expand Down
Loading

0 comments on commit 6ca18d3

Please sign in to comment.