Merge pull request #85 from efposadac/tcp

Multipole moment changes and QDO addition

bin/lowdin

@@ -203,6 +203,12 @@ if [ $extFile="lowdin" ]; then
 			}
                         else if($i=="m"){
 				printf("\tInputParticle_mass = %s\n",toupper($(i+2)) )
+			}
+                        else if($i=="omega"){
+				printf("\tInputParticle_omega = %s\n",toupper($(i+2)) )
+			}
+	                else if($i=="qdoCenterOf"){
+				printf("\tInputParticle_qdoCenterOf = \"%s\"\n",toupper($(i+2)) )
 			};
 		};
 		printf("\tInputParticle_origin = %15.12E %15.12E %15.12E \n/\n",$3,$4,$5)

lib/basis/AUG-CC-PVTZ

@@ -65,7 +65,29 @@ O-HYDROGEN EA- (AUG-CC-PVTZ) BASIS TYPE: 1
 9 2 1
 0.24700000 1.00000000
 
-
+O-HYDROGEN EB- (AUG-CC-PVTZ) BASIS TYPE: 1
+# 
+9
+1 0 3
+33.87000000 0.00606800
+5.09500000 0.04530800
+1.15900000 0.20282200
+2 0 1
+0.32580000 1.00000000
+3 0 1
+0.10270000 1.00000000
+4 0 1
+0.02526000 1.00000000
+5 1 1
+1.40700000 1.00000000
+6 1 1
+0.38800000 1.00000000
+7 1 1
+0.10200000 1.00000000
+8 2 1
+1.05700000 1.00000000
+9 2 1
+0.24700000 1.00000000
 
 O-HELIUM HE (AUG-CC-PVTZ) BASIS TYPE: 1
 # 

lib/basis/SHARON-E+6S2P

@@ -17,3 +17,43 @@ O-POSITRON E+ (sharon) BASIS TYPE: 2
 0.27492  1.0
 8 1 1
 0.084534  1.0
+
+O-POSITRON E+A (sharon) BASIS TYPE: 2
+#
+8
+1 0 1
+4.9231  1.0
+2 0 1
+0.92126 1.0
+3 0 1
+0.32058 1.0
+4 0 1
+0.14169  1.0
+5 0 1
+0.061803  1.0
+6 0 1
+0.0267577 1.0
+7 1 1
+0.27492  1.0
+8 1 1
+0.084534  1.0
+
+O-POSITRON E+B (sharon) BASIS TYPE: 2
+#
+8
+1 0 1
+4.9231  1.0
+2 0 1
+0.92126 1.0
+3 0 1
+0.32058 1.0
+4 0 1
+0.14169  1.0
+5 0 1
+0.061803  1.0
+6 0 1
+0.0267577 1.0
+7 1 1
+0.27492  1.0
+8 1 1
+0.084534  1.0

src/CalcProp/CalculateProperties.f90

@@ -91,7 +91,7 @@ module CalculateProperties_
        CalculateProperties_getPopulation, &
        CalculateProperties_showContributionsToElectrostaticMoment, &
        CalculateProperties_getDipoleOfPuntualCharges, &
-       CalculateProperties_getDipoleOfQuantumSpecie
+       CalculateProperties_getDipoleOfQuantumSpecies
   !     CalculateProperties_expectedR2, &
   !    CalculateProperties_polarizability, &
   !     CalculateProperties_showExpectedR2, &
@@ -538,7 +538,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
     write (6,"(T19,4A13)") "<Dx>","<Dy>", "<Dz>"," |D|"
 
     do i=1, numberOfSpecies
-       dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecie(this, i)
+       dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecies(this, i)
        totalDipole(:)=totalDipole(:)+dipole(i,:)
        write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecie( i )), dipole(i,:)
     end do
@@ -556,7 +556,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
 
     totalDipole=0.0_8
     do i=1, numberOfSpecies
-       dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecie(this, i)*2.54174619
+       dipole(i,:)=CalculateProperties_getDipoleOfQuantumSpecies(this, i)*2.54174619
        totalDipole(:)=totalDipole(:)+dipole(i,:)
        write (6,"(T5,A15,3F13.8)") trim(MolecularSystem_getNameOfSpecie( i )), dipole(i,:)
     end do
@@ -577,7 +577,7 @@ subroutine CalculateProperties_showContributionsToElectrostaticMoment(this)
     write (6,"(T19,6A13)") "<xx>","<yy>", "<zz>", "<xy>","<xz>","<yz>"
 
     do i=1, numberOfSpecies
-       quadrupole(i,:)=CalculateProperties_getQuadrupoleOfQuantumSpecie(this, i)*2.54174619*0.52917720859
+       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,:)
     end do
@@ -648,7 +648,7 @@ end function CalculateProperties_getQuadrupoleOfPuntualCharges
   !<
   !! @brief calcula el aporte al dipolo debido a particulas no fijas
   !>
-  function calculateproperties_getdipoleofquantumspecie( this, i ) result( output )
+  function CalculateProperties_getDipoleOfQuantumSpecies( this, i ) result( output )
     implicit none
     type(calculateproperties) :: this
     integer :: i !specieid
@@ -660,13 +660,13 @@ function calculateproperties_getdipoleofquantumspecie( this, i ) result( output
 
     output = output * molecularsystem_getcharge( i )
 
-  end function calculateproperties_getdipoleofquantumspecie
+  end function CalculateProperties_getDipoleOfQuantumSpecies
 
 
   !<
   !! @brief calcula el aporte al dipolo debido a particulas no fijas
   !>
-  function calculateproperties_getquadrupoleofquantumspecie( this, i ) result( output )
+  function CalculateProperties_getQuadrupoleOfQuantumSpecies( this, i ) result( output )
     implicit none
     type(calculateproperties) :: this
     integer :: i !specieid
@@ -681,7 +681,7 @@ function calculateproperties_getquadrupoleofquantumspecie( this, i ) result( out
 
     output = output * molecularsystem_getcharge( i )
 
-  end function calculateproperties_getquadrupoleofquantumspecie
+  end function CalculateProperties_getQuadrupoleOfQuantumSpecies
 
   subroutine CalculateProperties_exception( typeMessage, description, debugDescription)
     implicit none

src/core/CONTROL.f90

@@ -101,7 +101,7 @@ module CONTROL_
      logical :: HF_PRINT_EIGENVALUES
      character(20) :: HF_PRINT_EIGENVECTORS
      real(8) :: OVERLAP_EIGEN_THRESHOLD
-     real(8) :: ELECTRIC_FIELD(6)
+     real(8) :: ELECTRIC_FIELD(3)
      integer :: MULTIPOLE_ORDER
 
      !!***************************************************************************
@@ -442,7 +442,7 @@ module CONTROL_
   logical :: LowdinParameters_HFprintEigenvalues
   character(20) :: LowdinParameters_HFprintEigenvectors
   real(8) :: LowdinParameters_overlapEigenThreshold
-  real(8) :: LowdinParameters_electricField(6)
+  real(8) :: LowdinParameters_electricField(3)
   integer :: LowdinParameters_multipoleOrder
 
   !!***************************************************************************

src/core/InputManager.f90

@@ -344,6 +344,8 @@ subroutine InputManager_loadGeometry()
     real(8):: InputParticle_origin(3)
     real(8) :: InputParticle_charge
     real(8) :: InputParticle_mass
+    real(8) :: InputParticle_omega
+    character(15):: InputParticle_qdoCenterOf
     character(3):: InputParticle_fixedCoordinates
     integer:: InputParticle_addParticles
     real(8):: InputParticle_multiplicity    
@@ -357,6 +359,8 @@ subroutine InputManager_loadGeometry()
          InputParticle_basisSetName, &
          InputParticle_charge, &
          InputParticle_mass, &
+         InputParticle_omega, &
+         InputParticle_qdoCenterOf, &
          InputParticle_origin, &
          InputParticle_fixedCoordinates, &
          InputParticle_multiplicity, &
@@ -534,6 +538,8 @@ subroutine InputManager_loadGeometry()
        InputParticle_basisSetName = "NONE"
        InputParticle_charge=0.0_8
        InputParticle_mass=0.0_8
+       InputParticle_omega=0.0_8
+       InputParticle_qdoCenterOf = "NONE"
        InputParticle_origin=0.0_8
        InputParticle_fixedCoordinates = "NON"
        InputParticle_multiplicity = 1.0_8
@@ -614,7 +620,8 @@ subroutine InputManager_loadGeometry()
                   spin="ALPHA", &
                   id = particlesID(speciesID), &
                   charge = InputParticle_charge, &
-                  mass = InputParticle_mass )
+                  mass = InputParticle_mass, & 
+                  omega = InputParticle_omega )
 
              !!BETA SET
              speciesID = speciesID + 1
@@ -636,7 +643,8 @@ subroutine InputManager_loadGeometry()
                   spin="BETA", &
                   id = particlesID(speciesID), &
                   charge = InputParticle_charge, &
-                  mass = InputParticle_mass )             
+                  mass = InputParticle_mass, &
+                  omega = InputParticle_omega )         
 
           else 
 
@@ -662,7 +670,8 @@ subroutine InputManager_loadGeometry()
                   rotateAround=InputParticle_rotateAround,&                
                   id = particlesID(speciesID), &
                   charge = InputParticle_charge, &
-                  mass = InputParticle_mass )                         
+                  mass = InputParticle_mass, &
+                  omega = InputParticle_omega )                         
 
           end if
 
@@ -684,7 +693,8 @@ subroutine InputManager_loadGeometry()
                   rotationPoint=InputParticle_rotationPoint, &
                   rotateAround=InputParticle_rotateAround,&                
                   id = counter, &
-                  charge = InputParticle_charge )             
+                  charge = InputParticle_charge, &
+                  qdoCenterOf = InputParticle_qdoCenterOf )
 
           else
           !! Loads Particle
@@ -700,7 +710,8 @@ subroutine InputManager_loadGeometry()
                   rotationPoint=InputParticle_rotationPoint, &
                   rotateAround=InputParticle_rotateAround,&                
                   id = counter, &
-                  charge = InputParticle_charge )             
+                  charge = InputParticle_charge, &
+                  qdoCenterOf = InputParticle_qdoCenterOf )
           end if
        end if
     end do

src/core/MolecularSystem.f90

@@ -387,16 +387,17 @@ subroutine MolecularSystem_showParticlesInformation(this)
     !!
     print *,""
     print *,"                INFORMATION OF QUANTUM SPECIES "
-    write (6,"(T5,A70)") "---------------------------------------------------------------------"
-    write (6,"(T10,A2,A4,A8,A12,A4,A5,A6,A5,A4,A5,A12)") "ID", " ","Symbol", " ","mass", " ","charge", " ","spin","","multiplicity"
-    write (6,"(T5,A70)") "---------------------------------------------------------------------"
+    write (6,"(T5,A70)") "---------------------------------------------------------------------------------------------"
+    write (6,"(T10,A2,A4,A8,A12,A4,A5,A6,A5,A6,A5,A4,A5,A12)") "ID", " ","Symbol", " ","mass", " ","charge", " ","omega","","spin","","multiplicity"
+    write (6,"(T5,A70)") "---------------------------------------------------------------------------------------------"
 
     do i = 1, system%numberOfQuantumSpecies
-       write (6,'(T8,I3.0,A5,A10,A5,F10.4,A5,F5.2,A5,F5.2,A5,F5.2)') &
+       write (6,'(T8,I3.0,A5,A10,A5,F10.4,A5,F5.2,A5,F5.2,A5,F5.2,A5,F5.2)') &
             i, " ", &
             trim(system%species(i)%symbol)," ",&
             system%species(i)%mass," ",&
             system%species(i)%charge, " ",&
+            system%species(i)%omega," ",&
             system%species(i)%spin, "",&
             system%species(i)%multiplicity
     end do
@@ -659,6 +660,7 @@ subroutine MolecularSystem_saveToFile(targetFilePrefix)
     end do    
     !! Saving Point charges
     write(40,*) MolecularSystem_instance%numberOfPointCharges
+
     !! Saving info of each point charge
     do i = 1, MolecularSystem_instance%numberOfPointCharges
        call Particle_saveToFile(MolecularSystem_instance%pointCharges(i), unit=40)
@@ -1279,6 +1281,17 @@ function MolecularSystem_getCharge( speciesID ) result( output )
 
    end function MolecularSystem_getCharge
 
+   !> @brief Returns the omega frequency of speciesID. Why we have these functions??
+   function MolecularSystem_getOmega( speciesID ) result( output )
+     implicit none
+     integer :: speciesID
+
+     real(8) :: output
+
+     output = MolecularSystem_instance%species(speciesID)%omega
+
+   end function MolecularSystem_getOmega
+
    !> @brief Returns the mass of speciesID
    !! @author E. F. Posada, 2013
    !! @version 1.0   
@@ -1291,6 +1304,30 @@ function MolecularSystem_getMass( speciesID ) result( output )
      output = MolecularSystem_instance%species(speciesID)%mass
 
    end function MolecularSystem_getMass
+
+   !> @brief Returns QDO center of quantum species
+   function MolecularSystem_getQDOcenter( speciesID ) result( origin )
+     implicit none
+     integer :: speciesID
+     integer :: i
+     logical :: centerFound
+     real(8) :: origin(3)
+
+     centerFound = .False.
+     do i = 1 , size( MolecularSystem_instance%pointCharges )
+        if ( trim(MolecularSystem_instance%pointCharges(i)%qdoCenterOf) == trim(MolecularSystem_instance%species(speciesID)%symbol) ) then 
+          origin = MolecularSystem_instance%pointCharges(i)%origin 
+          centerFound = .True.
+          exit
+        end if
+    end do
+    if ( .not. centerFound ) then
+        call MolecularSystem_exception(ERROR, "No QDO center for species: "//MolecularSystem_instance%species(speciesID)%symbol, "MolecularSystem_getQDOcenter"   )
+    end if
+
+
+   end function MolecularSystem_getQDOCenter
+
 
    !> @brief Returns the Factor Of Exchange Integrals
    !! @author E. F. Posada, 2013
@@ -1404,6 +1441,14 @@ function MolecularSystem_getPointChargesEnergy() result( output )
 
         end do
      end do
+
+    !! Point charge potential with the external electric field
+    if ( sum(abs(CONTROL_instance%ELECTRIC_FIELD )) .ne. 0 ) then
+      do i=1, size( MolecularSystem_instance%pointCharges )      
+        output = output + sum(CONTROL_instance%ELECTRIC_FIELD(:) * MolecularSystem_instance%pointCharges(i)%origin(:) )*  MolecularSystem_instance%pointCharges(i)%charge 
+      end do
+    end if
+
 
    end function MolecularSystem_getPointChargesEnergy
 
@@ -1846,6 +1891,7 @@ subroutine MolecularSystem_copyConstructor(this,originalThis)
        this%species(i)%statistics = originalThis%species(i)%statistics
        this%species(i)%charge = originalThis%species(i)%charge
        this%species(i)%mass = originalThis%species(i)%mass
+       this%species(i)%omega = originalThis%species(i)%omega
        this%species(i)%spin = originalThis%species(i)%spin
        this%species(i)%totalCharge = originalThis%species(i)%totalCharge
        this%species(i)%kappa = originalThis%species(i)%kappa
@@ -1888,6 +1934,8 @@ subroutine MolecularSystem_copyConstructor(this,originalThis)
        this%allParticles(i)%particlePtr%origin=originalThis%allParticles(i)%particlePtr%origin
        this%allParticles(i)%particlePtr%charge=originalThis%allParticles(i)%particlePtr%charge
        this%allParticles(i)%particlePtr%mass=originalThis%allParticles(i)%particlePtr%mass
+       this%allParticles(i)%particlePtr%omega=originalThis%allParticles(i)%particlePtr%omega
+       this%allParticles(i)%particlePtr%qdoCenterOf=originalThis%allParticles(i)%particlePtr%qdoCenterOf
        this%allParticles(i)%particlePtr%spin=originalThis%allParticles(i)%particlePtr%spin
        this%allParticles(i)%particlePtr%totalCharge=originalThis%allParticles(i)%particlePtr%totalCharge
        this%allParticles(i)%particlePtr%klamt=originalThis%allParticles(i)%particlePtr%klamt
@@ -2048,6 +2096,7 @@ subroutine MolecularSystem_mergeTwoSystems(mergedThis,thisA,thisB,sysAbasisList,
        mergedThis%species(i)%statistics = thisA%species(i)%statistics
        mergedThis%species(i)%charge = thisA%species(i)%charge
        mergedThis%species(i)%mass = thisA%species(i)%mass
+       mergedThis%species(i)%omega = thisA%species(i)%omega
        mergedThis%species(i)%spin = thisA%species(i)%spin
        mergedThis%species(i)%totalCharge = thisA%species(i)%totalCharge
        mergedThis%species(i)%kappa = thisA%species(i)%kappa
@@ -2090,6 +2139,8 @@ subroutine MolecularSystem_mergeTwoSystems(mergedThis,thisA,thisB,sysAbasisList,
        mergedThis%pointCharges(i)%origin=thisA%pointCharges(i)%origin
        mergedThis%pointCharges(i)%charge=thisA%pointCharges(i)%charge
        mergedThis%pointCharges(i)%mass=thisA%pointCharges(i)%mass
+       mergedThis%pointCharges(i)%omega=thisA%pointCharges(i)%omega
+       mergedThis%pointCharges(i)%qdoCenterOf=thisA%pointCharges(i)%qdoCenterOf
        mergedThis%pointCharges(i)%spin=thisA%pointCharges(i)%spin
        mergedThis%pointCharges(i)%totalCharge=thisA%pointCharges(i)%totalCharge
        mergedThis%pointCharges(i)%klamt=thisA%pointCharges(i)%klamt
@@ -2131,6 +2182,8 @@ subroutine MolecularSystem_mergeTwoSystems(mergedThis,thisA,thisB,sysAbasisList,
                 mergedThis%species(i)%particles(j)%origin=thisA%species(i)%particles(jj)%origin
                 mergedThis%species(i)%particles(j)%charge=thisA%species(i)%particles(jj)%charge
                 mergedThis%species(i)%particles(j)%mass=thisA%species(i)%particles(jj)%mass
+                mergedThis%species(i)%particles(j)%omega=thisA%species(i)%particles(jj)%omega
+                mergedThis%species(i)%particles(j)%qdoCenterOf=thisA%species(i)%particles(jj)%qdoCenterOf
                 mergedThis%species(i)%particles(j)%spin=thisA%species(i)%particles(jj)%spin
                 mergedThis%species(i)%particles(j)%totalCharge=thisA%species(i)%particles(jj)%totalCharge
                 mergedThis%species(i)%particles(j)%klamt=thisA%species(i)%particles(jj)%klamt
@@ -2218,6 +2271,8 @@ subroutine MolecularSystem_mergeTwoSystems(mergedThis,thisA,thisB,sysAbasisList,
              mergedThis%species(i)%particles(j)%origin=thisB%species(i)%particles(jj)%origin
              mergedThis%species(i)%particles(j)%charge=thisB%species(i)%particles(jj)%charge
              mergedThis%species(i)%particles(j)%mass=thisB%species(i)%particles(jj)%mass
+             mergedThis%species(i)%particles(j)%omega=thisB%species(i)%particles(jj)%omega
+             mergedThis%species(i)%particles(j)%qdoCenterOf=thisB%species(i)%particles(jj)%qdoCenterOf
              mergedThis%species(i)%particles(j)%spin=thisB%species(i)%particles(jj)%spin
              mergedThis%species(i)%particles(j)%totalCharge=thisB%species(i)%particles(jj)%totalCharge
              mergedThis%species(i)%particles(j)%klamt=thisB%species(i)%particles(jj)%klamt