Merge pull request #60 from efposadac/NOCI+DFT

Changes to perform SCF core ionized states with the goal of doing XES:

bin/lowdin

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

examples/FCI_Charry2018/H-.FCI-DZ.lowdin

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

examples/FCI_Charry2018/PsH.FCI-DZ.lowdin

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

examples/FCI_Charry2018/e+H-H-.FCI-DZ.lowdin

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

lib/potentials/GRIBAKIN-POTENTIAL/generatePotential.py

@@ -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)) :
@@ -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()
@@ -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]

src/CI/ConfigurationInteraction.f90

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