Further search on the CEH gradient

src/tblite/ceh/coupled_perturbed.f90

@@ -116,14 +116,15 @@ subroutine get_density_matrix_gradient(mol,bas,wfn,list,dh0dr,dh0dL,doverlap,dde
                !denom = wfn%emo(jao,1)-wfn%emo(iao,1)
                ! occ-virt terms: u^(1)_qp = (F_mo^(1)_qp - S_mo^(1)_qp * E^(0)_pp) / (E^(0)_pp - E^(0)_qq)
                denom = wfn%emo(iao,1)-wfn%emo(jao,1)
-               u_mo(iao,jao) = u_mo(iao,jao) + (dh0dr_mo(iao,jao) - doverlap_mo(iao,jao) * wfn%emo(iao,1))/denom
-               u_mo(jao,iao) = -u_mo(iao,jao)   
+               u_mo(iao,jao) = u_mo(iao,jao) + (dh0dr_mo(iao,jao) - doverlap_mo(iao,jao) * wfn%emo(iao,1))/denom 
+               u_mo(jao,iao) = u_mo(jao,iao) - (dh0dr_mo(iao,jao) - doverlap_mo(iao,jao) * wfn%emo(iao,1))/denom 
                !denom = wfn%emo(jao,1)-wfn%emo(iao,1)
                !u_mo(jao,iao) = u_mo(jao,iao) + (dh0dr_mo(jao,iao) - doverlap_mo(jao,iao) * wfn%emo(jao,1))/denom
                !u_mo(jao,iao) = (dh0dr_mo(jao,iao) - doverlap_mo(jao,iao) * wfn%emo(iao,1))/denom
 
             end do 
          end do 
+         !call write_2d_matrix(u_mo, "u_mo")
          do iao = 1, bas%nao
             ! Diagonal terms: u^(1)_pp = -0.5 S_mo^(1)_pp
             u_mo(iao,iao) = -0.5_wp * doverlap_mo(iao,iao)

src/tblite/ceh/singlepoint.f90

@@ -20,6 +20,10 @@
 
 !> Implementation of the single point calculation for the CEH model
 module tblite_ceh_singlepoint
+
+
+   use iso_fortran_env, only: output_unit
+
    use mctc_env, only : error_type, wp
    use mctc_io, only: structure_type
    use tblite_adjlist, only : adjacency_list, new_adjacency_list
@@ -101,7 +105,7 @@ subroutine ceh_guess(ctx, calc, mol, error, wfn, accuracy, grad, verbosity)
       real(wp) :: nel, cutoff
       real(wp), allocatable :: tmp(:)
 
-      integer :: i, prlevel
+      integer :: ic, prlevel
 
       ! coordination number related arrays
       real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), cn_en(:), dcn_endr(:, :, :), dcn_endL(:, :, :)
@@ -193,8 +197,7 @@ subroutine ceh_guess(ctx, calc, mol, error, wfn, accuracy, grad, verbosity)
          ! Use the electronegativity-weighted CN as a 0th order guess for the charges
          call get_effective_qat(mol, cn_en, wfn%qat, &
          & dcn_endr, dcn_endL, wfn%dqatdr, wfn%dqatdL)
-         !write(*,*) "in ceh_guess qat", wfn%qat
-         !wfn%qsh = 0.0_wp
+
          call calc%coulomb%update(mol, ccache)
          call calc%coulomb%get_potential(mol, ccache, wfn, pot)
          if(grad) then
@@ -206,10 +209,6 @@ subroutine ceh_guess(ctx, calc, mol, error, wfn, accuracy, grad, verbosity)
       ! Add effective Hamiltonian to wavefunction
       call add_pot_to_h1(calc%bas, ints, pot, wfn%coeff)
 
-      !write(*,*) "pot%vat", pot%vat
-      !write(*,*) "pot%vsh", pot%vsh
-      !write(*,*) "pot%vao", pot%vao
-
       ! Solve the effective Hamiltonian
       call ctx%new_solver(solver, calc%bas%nao)
 
@@ -231,18 +230,21 @@ subroutine ceh_guess(ctx, calc, mol, error, wfn, accuracy, grad, verbosity)
 
       call timer%pop
       ttime = timer%get("wall time CEH")
-      
+
       call timer%push("wall time CEH gradient")
       if (grad) then
-         allocate(dh0dr(3,calc%bas%nao,calc%bas%nao),&
-         & dh0dL(3,calc%bas%nao,calc%bas%nao),doverlap(3,calc%bas%nao,calc%bas%nao),doverlap_diat(3,calc%bas%nao,calc%bas%nao), &
-         & source = 0.0_wp) !,&
-         !& ddensitydr(3,calc%bas%nao,calc%bas%nao,1),ddensitydL(3,calc%bas%nao,calc%bas%nao,1))
-
+         allocate(dh0dr(3,calc%bas%nao,calc%bas%nao), dh0dL(3,calc%bas%nao,calc%bas%nao), &
+         & doverlap(3,calc%bas%nao,calc%bas%nao),doverlap_diat(3,calc%bas%nao,calc%bas%nao), &
+         & source = 0.0_wp)
+
          ! Get the derivative of the Fock matrix elements
          call get_hamiltonian_gradient(mol, lattr, list, calc%bas, calc%h0, selfenergy, &
             & dsedr, dsedL, pot, doverlap, doverlap_diat, dh0dr, dh0dL)
-
+         !do ic = 1, 3
+         !   write(*,*) "ic", ic
+         !   call write_2d_matrix(doverlap(ic, :, :), "doverlap")
+         !   call write_2d_matrix(doverlap(ic, :, :), "dh0dr")
+         !end do
          ! Use the matrix element derivatives (F + S) to get the density matrix graidient
          ! based on the coupled-perturbed formalism
          call get_density_matrix_gradient(mol,calc%bas,wfn,list,dh0dr,dh0dL,doverlap,&
@@ -262,5 +264,51 @@ subroutine ceh_guess(ctx, calc, mol, error, wfn, accuracy, grad, verbosity)
       call timer%pop()
 
    end subroutine ceh_guess
+
+
+   subroutine write_2d_matrix(matrix, name, unit, step)
+      implicit none
+      real(wp), intent(in) :: matrix(:, :)
+      character(len=*), intent(in), optional :: name
+      integer, intent(in), optional :: unit
+      integer, intent(in), optional :: step
+      integer :: d1, d2
+      integer :: i, j, k, l, istep, iunit
+
+      d1 = size(matrix, dim=1)
+      d2 = size(matrix, dim=2)
+
+      if (present(unit)) then
+        iunit = unit
+      else
+        iunit = output_unit
+      end if
+
+      if (present(step)) then
+        istep = step
+      else
+        istep = 6
+      end if
+
+      if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name
+
+      do i = 1, d2, istep
+        l = min(i + istep - 1, d2)
+        write (iunit, '(/,6x)', advance='no')
+        do k = i, l
+          write (iunit, '(6x,i7,3x)', advance='no') k
+        end do
+        write (iunit, '(a)')
+        do j = 1, d1
+          write (iunit, '(i6)', advance='no') j
+          do k = i, l
+            write (iunit, '(1x,f15.12)', advance='no') matrix(j, k)
+          end do
+          write (iunit, '(a)')
+        end do
+      end do
+
+    end subroutine write_2d_matrix
+
 
 end module tblite_ceh_singlepoint