diff --git a/src/scf_module.F90 b/src/scf_module.F90 index af6c34d8c..d691eec65 100644 --- a/src/scf_module.F90 +++ b/src/scf_module.F90 @@ -173,7 +173,6 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & integer :: ndp,nqp integer,allocatable :: matlist (:,:) - integer,allocatable :: matlist2(:,:) integer,allocatable :: xblist(:,:) real(wp),allocatable :: sqrab(:) real(wp),allocatable :: dcndr(:,:,:) @@ -198,7 +197,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & integer :: ich ! file handle integer :: npr,ii,jj,kk,i,j,k,m,iat,jat,mi,jter,atj,kkk,mj,mm - integer :: ishell,jshell,np,ia,ndimv,l,nmat,nmat2 + integer :: ishell,jshell,np,ia,ndimv,l,nmat integer :: ll,i1,i2,nn,ati,nxb,lin,startpdiag integer :: is,js,gav @@ -326,8 +325,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & & X(basis%nao,basis%nao), & & ves(basis%nshell), & & zsh(basis%nshell),& - & matlist (2,basis%nao*(basis%nao+1)/2), & - & matlist2(2,basis%nao*(basis%nao+1)/2-basis%nao)) + & matlist (2,basis%nao*(basis%nao+1)/2)) allocate(selfEnergy(maxval(xtbData%nshell), mol%n)) allocate(dSEdcn(maxval(xtbData%nshell), mol%n)) @@ -587,24 +585,13 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & ! ------------------------------------------------------------------------ ! prepare matrix indices nmat =0 - nmat2=0 do ii=1,basis%nao iat=basis%aoat2(ii) do jj=1,ii-1 jat=basis%aoat2(jj) - if(abs(S(jj,ii)).lt.neglect) then - S(jj,ii)=0.0_wp - S(ii,jj)=0.0_wp - cycle - endif nmat=nmat+1 matlist(1,nmat)=ii matlist(2,nmat)=jj - if(iat.ne.jat)then - nmat2=nmat2+1 - matlist2(1,nmat2)=ii - matlist2(2,nmat2)=jj - endif enddo ! CB: moved this here so j/i indices from matlist come in a reasonable order ! also setup CN dep. stuff @@ -713,17 +700,21 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & if (mol%npbc == 0) then allocate(H(basis%nao, basis%nao)) H(:, :) = 0.0_wp - do m = 1, nmat2 - i = matlist2(1,m) - j = matlist2(2,m) - k = j+i*(i-1)/2 - !ishell = ao2sh(i) - !jshell = ao2sh(j) - ! SCC terms - !eh1 = autoev*(shellShift(ishell) + shellShift(jshell)) - !H1 = -S(j,i)*eh1*0.5_wp - H(j,i) = H0(k)*evtoau/S(j,i) - H(i,j) = H(j,i) + do i = 1, basis%nao + do j = 1, i-1 + k = j+i*(i-1)/2 + !ishell = ao2sh(i) + !jshell = ao2sh(j) + ! SCC terms + !eh1 = autoev*(shellShift(ishell) + shellShift(jshell)) + !H1 = -S(j,i)*eh1*0.5_wp + if (H0(k) .eq. 0.0_wp) then + H(j,i) = 0.0_wp + else + H(j,i) = H0(k)*evtoau/S(j,i) + end if + H(i,j) = H(j,i) + end do enddo call build_dSDQH0_noreset(xtbData%nShell, xtbData%hamiltonian, selfEnergy, & & dSEdcn, intcut, mol%n, basis%nao, basis%nbf, mol%at, mol%xyz, & diff --git a/test/unit/test_gfn2.f90 b/test/unit/test_gfn2.f90 index 1d18851fa..6860271e1 100644 --- a/test/unit/test_gfn2.f90 +++ b/test/unit/test_gfn2.f90 @@ -299,11 +299,11 @@ subroutine test_gfn2gbsa_api(error) call check_(error, hl_gap, 3.408607724814_wp, thr=1e-5_wp) call check_(error, energy,-22.002501380096_wp, thr=thr) - call check_(error, norm2(gradient),0.19441977481008E-01_wp, thr=thr) + call check_(error, norm2(gradient),0.19366866479022E-01_wp, thr=thr) call check_(error, gradient(1,1), .9038674439127e-02_wp, thr=thr) call check_(error, gradient(3,2),-.1394693523214e-02_wp, thr=thr) call check_(error, gradient(3,11),-gradient(3,10), thr=thr) - call check_(error, gradient(1,8),0.22890674680144E-02_wp, thr=thr) + call check_(error, gradient(1,8),0.17878350605942E-02_wp, thr=thr) end subroutine test_gfn2gbsa_api @@ -1080,7 +1080,7 @@ subroutine test_gfn2_wbo(error) ! check scc ! call check_(error, energy, -33.314958144107_wp, thr=thr) call check_(error, norm2(gradient), 0.018945181484_wp, thr=thr) - call check_(error, hl_gap, 1.805948277212_wp, thr=thr) + call check_(error, hl_gap, 1.805948321662_wp, thr=thr) ! check wbo !