Merge pull request #98 from efposadac/ROCI

Roci
efposadac · Nov 19, 2024 · 7b0b800 · 7b0b800
2 parents f31d392 + 9cd5a9f
commit 7b0b800
Show file tree

Hide file tree

Showing 57 changed files with 6,555 additions and 4,768 deletions.
diff --git a/bin/lowdin b/bin/lowdin
@@ -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
 
@@ -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
@@ -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

diff --git a/lib/basis/NAKAI-CC-PVDZ b/lib/basis/NAKAI-CC-PVDZ
@@ -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
+
diff --git a/src/CI/CIStrings.f90 b/src/CI/CIStrings.f90
@@ -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) 
 

diff --git a/src/CI/CIcore.f90 b/src/CI/CIcore.f90
@@ -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

diff --git a/src/CI/CImod.f90 b/src/CI/CImod.f90
@@ -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: "
@@ -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 )
@@ -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 = &
@@ -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 )
@@ -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)
@@ -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
@@ -1472,7 +1472,7 @@ subroutine CImod_densityMatrices()
                 write(*,*) "-----------------"
 
                 numberOfContractions = MolecularSystem_getTotalNumberOfContractions( species )
-                speciesName = MolecularSystem_getNameOfSpecie(species)
+                speciesName = MolecularSystem_getNameOfSpecies(species)
 
 
                 call Vector_constructor ( auxdensityEigenValues, &
@@ -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 )
 

diff --git a/src/CI/Configuration.f90 b/src/CI/Configuration.f90
@@ -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"

diff --git a/src/CalcProp/CalculateProperties.f90 b/src/CalcProp/CalculateProperties.f90
@@ -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))
@@ -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))
@@ -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))
@@ -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-"
@@ -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
@@ -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
@@ -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"
@@ -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,:)
@@ -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
@@ -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
@@ -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) )
@@ -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>"

diff --git a/src/DFT/DFT.f90 b/src/DFT/DFT.f90
@@ -28,6 +28,7 @@ program DFT
   use MolecularSystem_
   use DensityFunctionalTheory_
   use GridManager_
+  use Functional_
   use String_
   use Matrix_
   use Exception_
@@ -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
@@ -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))
 
@@ -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
@@ -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: "