New orbital removal routine, to remove the lowest eigenvalues (previously was comparing to OVERLAP_EIGEN_THRESHOLD)

Author: fsmoncadaa

src/scf/WaveFunction.f90

@@ -2739,36 +2739,48 @@ subroutine Wavefunction_removeOrbitalsBelowEigenThreshold(this)
     type(WaveFunction) :: this
 
     integer(8) :: numberOfContractions
-    real(8) :: normCheck    
+    type(Vector) :: normCheck, normCheckSorted
+    real(8) :: threshold
     integer :: i, j, mu, nu, index
 
+    if(this%removedOrbitals .eq. 0) return
+    
     numberOfContractions = MolecularSystem_getTotalnumberOfContractions(this%species,this%molSys)
+    call Vector_constructor(normCheck, int(numberOfContractions,4), 0.0_8)
 
-    i=0
-    do index = 1 , numberOfContractions
-       i=i+1
-       normCheck=0.0
+    do i=1, numberOfContractions
        do mu = 1 , numberOfContractions
           do nu = 1 , numberOfContractions
-             normCheck=normCheck+this%waveFunctionCoefficients%values(mu,i)*&
+             normCheck%values(i)=normCheck%values(i)+this%waveFunctionCoefficients%values(mu,i)*&
                   this%waveFunctionCoefficients%values(nu,i)*&
                   this%overlapMatrix%values(mu,nu)
           end do
        end do
-       if ( normCheck .lt. CONTROL_instance%OVERLAP_EIGEN_THRESHOLD) then
-          if (  CONTROL_instance%DEBUG_SCFS) &
-               print *, "shifting eigenvector no.", i, "with normCheck", normCheck, "to the end of the coefficients matrix"
+    end do
+
+    normCheckSorted=normCheck
+    call Vector_reverseSortElements(normCheckSorted)
+    threshold=normCheckSorted%values(this%removedOrbitals)
 
+    i=0
+    do index = 1 , numberOfContractions
+       i=i+1
+       if ( normCheck%values(i) .le. threshold) then
+          if (  CONTROL_instance%DEBUG_SCFS) &
+               print *, "shifting eigenvector no.", i, "with normCheck", normCheck%values(i), "to the end of the coefficients matrix"
           do j = i , numberOfContractions-1
              this%molecularOrbitalsEnergy%values(j)=this%molecularOrbitalsEnergy%values(j+1)
-             this%waveFunctionCoefficients%values(:,j) = this%waveFunctionCoefficients%values(:,j+1)
+             this%waveFunctionCoefficients%values(1:numberOfContractions,j) = this%waveFunctionCoefficients%values(1:numberOfContractions,j+1)
+             normCheck%values(j)=normCheck%values(j+1)
           end do
+          i=i-1
           ! Make eigenenergy a very large number
           this%molecularOrbitalsEnergy%values(numberOfContractions)=1.0E+308_8
-          this%waveFunctionCoefficients%values(:,numberOfContractions)=0.0
-          i=i-1
+          this%waveFunctionCoefficients%values(1:numberOfContractions,numberOfContractions)=0.0_8
+          normCheck%values(numberOfContractions)=1.0E+308_8
        end if
     end do
+    
   end subroutine Wavefunction_removeOrbitalsBelowEigenThreshold
 
 

test/QDO-largeBasis.lowdin

@@ -6,10 +6,12 @@ TASKS
 method = "UHF"
 END TASKS
 CONTROL
+readCoefficients=F
 setQDOenergyZero=T
 totalEnergyTolerance=1E-12
 overlapEigenThreshold=1E-6
 integralStorage="DIRECT"
+debugSCFS=F
 END CONTROL
 BASIS PB6ZQDO
 O-POSITIVE qp+ (PB6ZQDO) BASIS TYPE: 2