Near-linear basis set dependency handling and level shifting for RHF/UHF/DFT/UDFT

src/getEnergy.f90

@@ -95,8 +95,15 @@ subroutine getEnergy(isGuess, ierr)
 #ifdef MPIV
    !-------------- MPI / ALL NODES ----------------------------------
    if (bMPI) then
+      quick_molspec%NBSuse   => NBSuse
+      quick_qm_struct%NBSuse => NBSuse
+
+      call MPI_BCAST(NBSuse,1,mpi_integer,0,MPI_COMM_WORLD,mpierror)
+
+      if(.not. master) call allocate_quick_qm_struct_fullx(quick_qm_struct)
+
       call MPI_BCAST(quick_qm_struct%s,nbasis*nbasis,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
-      call MPI_BCAST(quick_qm_struct%x,nbasis*nbasis,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
+      call MPI_BCAST(quick_qm_struct%x,nbasis*NBSuse,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(quick_qm_struct%Ecore,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
    endif
    !-------------- END MPI / ALL NODES ------------------------------

src/modules/quick_basis_module.f90

@@ -51,6 +51,7 @@ module quick_basis_module
 
         ! total basis number
         integer, pointer :: nbasis
+        integer, pointer :: NBSuse
 
         ! the first and last basis function for an atom
         integer, dimension(:),allocatable :: first_basis_function, last_basis_function
@@ -128,7 +129,7 @@ module quick_basis_module
 
 
    integer,target :: nshell,nprim,jshell,jbasis
-   integer,target :: nbasis
+   integer,target :: nbasis, NBSuse
    integer :: maxcontract
 
    ! used for 2e integral indices

src/modules/quick_calculated_module.f90

@@ -32,9 +32,12 @@ module quick_calculated_module
    ! the elements will be introduced following
    type quick_qm_struct_type
 
-      ! Basis Set Number
+      ! number of basis functions
       integer,pointer :: nbasis
 
+      ! number of basis functions accounting for near-linear dependency
+      integer,pointer :: NBSuse
+
       ! overlap matrix, will be calculated once, independent on
       ! orbital coefficent. Its dimension is nbasis*nbasis
       double precision,dimension(:,:), allocatable :: s
@@ -49,6 +52,12 @@ module quick_calculated_module
       ! operator matrix, the dimension is nbasis*nbasis. For HF, it's Fock Matrix
       double precision,dimension(:,:), allocatable :: o
 
+      ! effective operator matrix if basis functions are eliminated to remove near-linear dependency
+      double precision,dimension(:,:), allocatable :: oeff
+
+      ! effective beta operator matrix if basis functions are eliminated to remove near-linear dependency
+      double precision,dimension(:,:), allocatable :: oeffb
+
       ! matrix for saving XC potential, required for incremental KS build
       double precision,dimension(:,:), allocatable :: oxc 
 
@@ -77,6 +86,14 @@ module quick_calculated_module
       ! the dimension is nbasis*nbasis.
       double precision,dimension(:,:), allocatable :: vec
 
+      ! matrix to hold eigenvectors for level shifting,
+      ! the dimension is NBSuse*NBSuse.
+      double precision,dimension(:,:), allocatable :: oldvec
+
+      ! matrix to hold beta eigenvectors for level shifting,
+      ! the dimension is NBSuse*NBSuse.
+      double precision,dimension(:,:), allocatable :: oldvecb
+
       ! Density matrix, when it's unrestricted system, it presents alpha density
       ! the dimension is nbasis*nbasis.
       double precision,dimension(:,:), allocatable :: dense
@@ -215,6 +232,54 @@ module quick_calculated_module
    !----------------------
 contains
 
+   !--------------------------------------
+   ! subroutine to allocate lists whose
+   ! dimensions depend on near-linear
+   ! dependency
+   !--------------------------------------
+   subroutine allocate_quick_qm_struct_fullx(self)
+       use quick_molspec_module, only: quick_molspec
+       use quick_method_module,  only: quick_method
+
+       implicit none
+
+       type (quick_qm_struct_type) self
+
+       ! alpha fields
+       if(self%NBSuse.ne.self%nbasis)then
+           if(.not. allocated(self%oeff))allocate(self%oeff(self%NBSuse,self%NBSuse))
+           self%oeff = 0.0d0
+       endif
+       if(.not. allocated(self%x))      allocate(self%x(self%nbasis,self%NBSuse))
+       if(.not. allocated(self%vec))    allocate(self%vec(self%NBSuse,self%NBSuse))
+       if(.not. allocated(self%oldvec)) allocate(self%oldvec(self%NBSuse,self%NBSuse))
+       if(.not. allocated(self%co))     allocate(self%co(self%nbasis,self%NBSuse))
+       if(.not. allocated(self%E))      allocate(self%E(self%NBSuse))
+
+       self%x      = 0.0d0
+       self%vec    = 0.0d0
+       self%oldvec = 0.0d0
+       self%co     = 0.0d0
+       self%E      = 0.0d0
+
+       ! beta fields (unrestricted only): cob and Eb are resized to NBSuse here;
+       ! oeffb is only needed when NBSuse < nbasis; oldvecb is always needed.
+       if(quick_method%unrst) then
+          if(self%NBSuse.ne.self%nbasis)then
+             if(.not. allocated(self%oeffb)) allocate(self%oeffb(self%NBSuse,self%NBSuse))
+             self%oeffb = 0.0d0
+          endif
+          if(.not. allocated(self%cob))     allocate(self%cob(self%nbasis,self%NBSuse))
+          if(.not. allocated(self%Eb))      allocate(self%Eb(self%NBSuse))
+          if(.not. allocated(self%oldvecb)) allocate(self%oldvecb(self%NBSuse,self%NBSuse))
+
+          self%cob     = 0.0d0
+          self%Eb      = 0.0d0
+          self%oldvecb = 0.0d0
+       endif
+
+   end subroutine
+
    !--------------------------------------
    ! subroutine to allocate variables
    !---------------------------------------
@@ -238,16 +303,12 @@ subroutine allocate_quick_qm_struct(self)
 
       ! those matrices is necessary for all calculation or the basic of other calculation
       if(.not. allocated(self%s)) allocate(self%s(nbasis,nbasis))
-      if(.not. allocated(self%x)) allocate(self%x(nbasis,nbasis))
       if(.not. allocated(self%oneElecO)) allocate(self%oneElecO(nbasis,nbasis))
       if(.not. allocated(self%o)) allocate(self%o(nbasis,nbasis))
       if(.not. allocated(self%oSave)) allocate(self%oSave(nbasis,nbasis))
-      if(.not. allocated(self%co)) allocate(self%co(nbasis,nbasis))
-      if(.not. allocated(self%vec)) allocate(self%vec(nbasis,nbasis))
       if(.not. allocated(self%dense)) allocate(self%dense(nbasis,nbasis))
       if(.not. allocated(self%denseSave)) allocate(self%denseSave(nbasis,nbasis))
       if(.not. allocated(self%denseOld)) allocate(self%denseOld(nbasis,nbasis))
-      if(.not. allocated(self%E)) allocate(self%E(nbasis))
       if(.not. allocated(self%iDegen)) allocate(self%iDegen(nbasis))
 
       if(.not. allocated(self%Mulliken)) allocate(self%Mulliken(natom))
@@ -274,16 +335,14 @@ subroutine allocate_quick_qm_struct(self)
          if(.not. allocated(self%CPHFB)) allocate(self%CPHFB(idimA,natom*3))
       endif
 
-      ! if unrestricted, some more varibles is required to be allocated
-      if (quick_method%unrst) then
-         if(.not. allocated(self%ob)) allocate(self%ob(nbasis,nbasis))
-         if(.not. allocated(self%obSave)) allocate(self%obSave(nbasis,nbasis))
-         if(.not. allocated(self%denseab)) allocate(self%denseab(nbasis,nbasis))
-         if(.not. allocated(self%densebSave)) allocate(self%densebSave(nbasis,nbasis))
-         if(.not. allocated(self%densebOld)) allocate(self%densebOld(nbasis,nbasis))
-         if(.not. allocated(self%cob)) allocate(self%cob(nbasis,nbasis))
-         if(.not. allocated(self%Eb)) allocate(self%Eb(nbasis))
-      endif
+       ! if unrestricted, some more varibles is required to be allocated
+       if (quick_method%unrst) then
+          if(.not. allocated(self%ob)) allocate(self%ob(nbasis,nbasis))
+          if(.not. allocated(self%obSave)) allocate(self%obSave(nbasis,nbasis))
+          if(.not. allocated(self%denseab)) allocate(self%denseab(nbasis,nbasis))
+          if(.not. allocated(self%densebSave)) allocate(self%densebSave(nbasis,nbasis))
+          if(.not. allocated(self%densebOld)) allocate(self%densebOld(nbasis,nbasis))
+       endif
 
       if (quick_method%unrst .or. quick_method%DFT) then
          if(.not. allocated(self%denseb)) allocate(self%denseb(nbasis,nbasis))
@@ -347,6 +406,7 @@ subroutine dat_quick_qm_struct(self, idatafile)
       integer fail
 
       integer nbasis
+      integer NBSuse
       integer natom
       integer nelec
       integer idimA
@@ -357,30 +417,30 @@ subroutine dat_quick_qm_struct(self, idatafile)
       type (quick_qm_struct_type) self
 
       nbasis=self%nbasis
+      NBSuse=self%NBSuse
       natom=quick_molspec%natom
       nelec=quick_molspec%nelec
       nelecb=quick_molspec%nelecb
 
-
       call write_int_rank0(idatafile, "nbasis", nbasis, fail)
       call write_int_rank0(idatafile, "natom",  natom, fail)
       call write_int_rank0(idatafile, "nelec",  nelec, fail)
       call write_int_rank0(idatafile, "nelecb", nelecb, fail)
       call write_real8_rank3(idatafile, "s", nbasis, nbasis, 1, self%s, fail)
-      call write_real8_rank3(idatafile, "x", nbasis, nbasis, 1, self%x, fail)
+      call write_real8_rank3(idatafile, "x", nbasis, NBSuse, 1, self%x, fail)
       call write_real8_rank3(idatafile, "o", nbasis, nbasis, 1, self%o, fail)
-      call write_real8_rank3(idatafile, "co", nbasis, nbasis, 1, self%co, fail)
-      call write_real8_rank3(idatafile, "vec", nbasis, nbasis, 1, self%vec, fail)
+      call write_real8_rank3(idatafile, "co", nbasis, NBSuse, 1, self%co, fail)
+      call write_real8_rank3(idatafile, "vec", NBSuse, NBSuse, 1, self%vec, fail)
       call write_real8_rank3(idatafile, "dense", nbasis, nbasis, 1, self%dense, fail)
-      call write_real8_rank3(idatafile, "E", nbasis, 1, 1, self%E, fail)
+      call write_real8_rank3(idatafile, "E", NBSuse, 1, 1, self%E, fail)
       call write_int_rank3(idatafile, "iDegen", nbasis, 1, 1, self%iDegen, fail)
       call write_real8_rank3(idatafile, "Mulliken", nbasis, 1, 1, self%Mulliken, fail)
       call write_real8_rank3(idatafile, "Lowdin", nbasis, 1, 1, self%Lowdin, fail)
 
       ! if unrestricted, some more varibles is required to be allocated
       if (quick_method%unrst) then
-         call write_real8_rank3(idatafile, "cob", nbasis, nbasis, 1, self%cob, fail)
-         call write_real8_rank3(idatafile, "Eb", nbasis, 1, 1, self%Eb, fail)
+         call write_real8_rank3(idatafile, "cob", nbasis, NBSuse, 1, self%cob, fail)
+         call write_real8_rank3(idatafile, "Eb", NBSuse, 1, 1, self%Eb, fail)
       endif
 
       if (quick_method%unrst .or. quick_method%DFT) then
@@ -398,7 +458,7 @@ subroutine deallocate_quick_qm_struct(self)
       implicit none
       integer io
 
-      integer nbasis
+      integer nbasis, NBSuse
       integer natom
       integer nelec
       integer idimA
@@ -407,14 +467,19 @@ subroutine deallocate_quick_qm_struct(self)
       type (quick_qm_struct_type) self
 
       nullify(self%nbasis)
+      nullify(self%NBSuse)
       ! those matrices is necessary for all calculation or the basic of other calculation
       if (allocated(self%s)) deallocate(self%s)
       if (allocated(self%x)) deallocate(self%x)
       if (allocated(self%oneElecO)) deallocate(self%oneElecO)
       if (allocated(self%o)) deallocate(self%o)
+      if (allocated(self%oeff)) deallocate(self%oeff)
+      if (allocated(self%oeffb)) deallocate(self%oeffb)
       if (allocated(self%oSave)) deallocate(self%oSave)
       if (allocated(self%co)) deallocate(self%co)
       if (allocated(self%vec)) deallocate(self%vec)
+      if (allocated(self%oldvec)) deallocate(self%oldvec)
+      if (allocated(self%oldvecb)) deallocate(self%oldvecb)
       if (allocated(self%dense)) deallocate(self%dense)
       if (allocated(self%denseSave)) deallocate(self%denseSave)
       if (allocated(self%denseOld)) deallocate(self%denseOld)
@@ -496,6 +561,7 @@ subroutine broadcast_quick_qm_struct(self)
       call MPI_BCAST(self%oSave,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(self%co,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(self%vec,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
+      call MPI_BCAST(self%oldvec,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(self%dense,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(self%denseSave,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
       call MPI_BCAST(self%denseOld,nbasis2,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
@@ -568,16 +634,18 @@ subroutine init_quick_qm_struct(self)
       nelecb=quick_molspec%nelecb
 
       call zeroMatrix(self%s,nbasis)
-      call zeroMatrix(self%x,nbasis)
+!      call zeroMatrix(self%x,nbasis)
       call zeroMatrix(self%oneElecO,nbasis)
       call zeroMatrix(self%o,nbasis)
+!      call zeroMatrix(self%oeff,nbasis)
       call zeroMatrix(self%oSave,nbasis)
-      call zeroMatrix(self%co,nbasis)
-      call zeroMatrix(self%vec,nbasis)
+!      call zeroMatrix(self%co,nbasis)
+!      call zeroMatrix(self%vec,nbasis)
+!      call zeroMatrix(self%oldvec,nbasis)
       call zeroMatrix(self%dense,nbasis)
       call zeroMatrix(self%denseSave,nbasis)
       call zeroMatrix(self%denseOld,nbasis)
-      call zeroVec(self%E,nbasis)
+!      call zeroVec(self%E,nbasis)
       call zeroiVec(self%iDegen,nbasis)
       call zeroVec(self%Mulliken,natom)
       call zeroVec(self%Lowdin,natom)
@@ -603,10 +671,10 @@ subroutine init_quick_qm_struct(self)
 
       ! if unrestricted, some more varibles is required to be allocated
       if (quick_method%unrst) then
-         call zeroMatrix(self%cob,nbasis)
+         if (allocated(self%cob))  call zeroMatrix(self%cob,nbasis)
          call zeroMatrix(self%denseab,nbasis)
          call zeroMatrix(self%denseb,nbasis)
-         call zeroVec(self%Eb,nbasis)
+         if (allocated(self%Eb))   call zeroVec(self%Eb,self%NBSuse)
       endif
 
    end subroutine

src/modules/quick_method_module.f90

@@ -94,6 +94,11 @@ module quick_method_module
         ! this is DFT grid
         integer :: iSG = 1             ! =0. SG0, =1. SG1(DEFAULT)
 
+        ! Level shift
+        integer :: LShift_cycle = 3              ! After what cycle allow Level shifting
+        double precision :: LShift_err = 0.1d0   ! Minimum error for allowing Level shifting
+        double precision :: LShift_gap = 0.2d0   ! HOMO-LUMO gap after Level shifting
+
         ! Initial guess part
         logical :: SAD = .true.        ! SAD initial guess(default)
         logical :: MFCC = .false.      ! MFCC
@@ -126,6 +131,8 @@ module quick_method_module
         !tol
         double precision :: pmaxrms        = 1.0d-6   ! density matrix convergence criteria
         double precision :: basisCutoff    = 1.0d-6  ! basis set cutoff
+        double precision :: overlapCutoff  = 1.0d-5  ! cutoff for near-linear dependency
+        double precision :: ovmatelems     = 1.0d-6  ! cutoff to consider overlap matrix elements
         !signif
 
         ! following are some gradient cutoff criteria
@@ -455,6 +462,10 @@ subroutine print_quick_method(self,io,ierr)
                 if (self%iSG .eq. 1) write(io,'(" STANDARD GRID = SG1")')
             endif
 
+           write(io,'(" Level shifting allowed after cycle ", I4)') self%LShift_cycle
+           write(io,'(" DIIS error must exceed ", F5.3," for level shifting")') self%LShift_err
+           write(io,'(" HOMO-LUMO gap after level shifting = ", F5.3)') self%LShift_gap
+
             if (self%opt) then
                 write(io,'(" GEOMETRY OPTIMIZATION")',advance="no")
                 if (self%diisOpt)   write(io,'(" USE DIIS FOR GEOMETRY OPTIMIZATION")')
@@ -524,10 +535,11 @@ subroutine print_quick_method(self,io,ierr)
 
             ! cutoff size
             write (io,'(" COMPUTATIONAL CUTOFF: ")')
-            write (io,'("      TWO-e INTEGRAL   = ",E10.3)') self%integralcutoff
-            write (io,'("      BASIS SET PRIME  = ",E10.3)') self%primLimit
-            write (io,'("      MATRIX ELEMENTS  = ",E10.3)') self%DMCutoff
-            write (io,'("      BASIS FUNCTION   = ",E10.3)') self%basisCutoff
+            write (io,'("      TWO-e INTEGRAL         = ",E10.3)') self%integralcutoff
+            write (io,'("      BASIS SET PRIME        = ",E10.3)') self%primLimit
+            write (io,'("      MATRIX ELEMENTS        = ",E10.3)') self%DMCutoff
+            write (io,'("      BASIS FUNCTION         = ",E10.3)') self%basisCutoff
+            write (io,'("      NEAR-LINEAR DEPENDENCY = ",E10.3)') self%overlapCutoff
             if (self%grad) then
                 write (io,'("      GRADIENT CUTOFF  = ",E10.3)') self%gradCutoff
             endif
@@ -772,6 +784,11 @@ subroutine read_quick_method(self,keywd,ierr)
                 self%primLimit=self%integralCutoff*1.0d-1
             endif
 
+            ! Overlap-cutoff
+            if (index(keywd,'OVCUT') /= 0) then
+                call read(keywd, 'OVCUT', self%overlapCutoff)
+            endif
+
             ! Grad cutoff
             if (index(keywd,'GRADCUTOFF') /= 0) then
                 call read(keywd, 'GRADCUTOFF', self%gradCutoff)
@@ -875,6 +892,15 @@ subroutine read_quick_method(self,keywd,ierr)
                self%extgrid_angstrom=.true.
                self%ext_grid=.true.
            endif
+           if (index(keyWD,'LSHIFT_CYCLE').ne.0) then
+               call read(keywd,'LSHIFT_CYCLE', self%LShift_cycle)
+           endif
+           if (index(keyWD,'LSHIFT_ERR').ne.0) then
+               call read(keywd,'LSHIFT_ERR', self%LShift_err)
+           endif
+           if (index(keyWD,'LSHIFT_GAP').ne.0) then
+               call read(keywd,'LSHIFT_GAP', self%LShift_gap)
+           endif
         end subroutine read_quick_method
 
 
@@ -925,6 +951,10 @@ subroutine init_quick_method(self,ierr)
             self%efield_grid = .false.     ! Electric field (EFIELD) evaluated on grid
             self%efg_grid = .false.        ! Electric field gradient (EFG)
 
+            self%LShift_cycle = 3     ! After what cycle allow Level shifting
+            self%LShift_err = 0.1d0   ! Minimum error for allowing Level shifting
+            self%LShift_gap = 0.2d0   ! HOMO-LUMO gap after Level shifting
+
             self%diisOpt =  .false.  ! DIIS Optimization
             self%core =  .false.     !
             self%annil =  .false.    !
@@ -972,6 +1002,8 @@ subroutine init_quick_method(self,ierr)
 
             self%pmaxrms        = 1.0d-6   ! density matrix convergence criteria
             self%basisCutoff    = 1.0d-6  ! basis set cutoff
+            self%overlapCutoff  = 1.0d-5  ! Near-linear dependency check cutoff
+            self%ovmatelems     = 1.0d-6  ! cutoff to consider overlap matrix elements
 
             self%stepMax        = .1d0/0.529177249d0
                                            ! max change of one step

src/modules/quick_molspec_module.f90

@@ -85,6 +85,7 @@ module quick_molspec_module
 
       ! basis set number
       integer, pointer:: nbasis
+      integer, pointer:: NBSuse
 
    end type quick_molspec_type