Final adjustments of GFN-FF for Ln/An

src/approxrab.f90

@@ -46,7 +46,7 @@ module xtb_approxrab
       2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898,& ! 80
       1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095,&
       2.34224386,&
-      2.52387890,& !@thomas TODO
+      2.52387890,&
       2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000,& ! 92
       2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,&
       2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,& ! 102
@@ -71,7 +71,7 @@ module xtb_approxrab
       2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774,&
       2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552,&
       2.63264944,&
-      4.24537037,& !@thomas TODO
+      4.24537037,&
       3.66542289, 4.20000000, 4.20000000, 4.20000000, 4.20000000,& ! 92
       4.20000000, 4.20000000, 4.20000000, 4.20000000, 4.20000000,&
       4.20000000, 4.20000000, 4.20000000, 4.20000000, 4.20000000,& ! 102
@@ -96,7 +96,7 @@ module xtb_approxrab
       0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502,& ! 90
       0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675,&
       0.04800662,&
-      0.04582912,& !@thomas TODO
+      0.04582912,&
       0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441,& ! 92
       0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969,&
      -0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000,& ! 102

src/geoopt_driver.f90

@@ -162,7 +162,7 @@ subroutine geometry_optimization &
      ! create new Limited-memory BFGS optimizer 
      call new_lbfgs_optimizer(lbfgs_opt, env, opt_input, filter)
      ! run optimization
-     call relax_pbc(lbfgs_opt, env, mol, wfn, calc, filter, printlevel)
+     call relax_pbc(lbfgs_opt, env, mol, wfn, calc, filter, printlevel, fail)
    case(p_engine_inertial)
       call fire & ! FIRE !
          &(env,ilog,mol,wfn,calc, &

src/gfnff/gdisp0.f90

@@ -130,6 +130,7 @@ subroutine weight_references_d4(dispm, nat, atoms, wf, cn, gwvec, gwdcn)
          norm = norm + gw
          dnorm = dnorm + 2*wf*(dispm%cn(iref, ati) - cn(iat)) * gw
       end do
+      if (norm.eq.0.0_wp) cycle
       norm = 1.0_wp / norm
       do iref = 1, dispm%nref(ati)
          expw = weight_cn(wf, cn(iat), dispm%cn(iref, ati))

src/gfnff/gfnff_ini.f90

@@ -2207,6 +2207,7 @@ subroutine adjust_NB_LnH_AnH(param, mol, topo, neigh)
   type(TNeigh), intent(inout) :: neigh ! main type for introducing PBC
   integer nb_tmp(neigh%numnb)
   integer :: i,iTr,iTrH,idx,inb,l,k,m,nnb,count_idx
+  real(wp), parameter :: qthr = -0.0281_wp ! charge threshold
 
   ! loop over all atoms
   do i=1, mol%n
@@ -2224,7 +2225,7 @@ subroutine adjust_NB_LnH_AnH(param, mol, topo, neigh)
           ! check if neighbor is H
           if (mol%at(inb).eq.1) then
             ! remove hydrogen as neighbor if charge is gt threshold
-            if (topo%qa(inb).gt.-0.0282) then
+            if (topo%qa(inb).gt.qthr) then
               ! setup copy of neighbor list for this An or Ln
               nb_tmp = neigh%nb(:,i,iTr)
               nb_tmp(neigh%numnb) = neigh%nb(neigh%numnb,i,iTr) - 1
@@ -2452,10 +2453,10 @@ subroutine gfnff_topo_changes(env, neigh)
    ! check if hardcoded size of ffnb is still up to date
    if (size(set%ffnb, dim=1).ne.neigh%numnb) call env%error('The array set%ffnb has not been adjusted to changes in neigh%numnb.', source)
    ! only do something if there are changes stored in set%ffnb
-   if(set%ffnb(1,1).ne.0) then
+   if(set%ffnb(1,1).ne.-1) then
       d2=size(set%ffnb, dim=2)
       do i=1, d2
-         if (set%ffnb(1,i).eq.0) exit
+         if (set%ffnb(1,i).eq.-1) exit
          idx=set%ffnb(1,i)
          int_tmp = set%ffnb(2:41,i)
          neigh%nb(1:40,idx,1) = int_tmp

src/gfnff/gfnff_ini2.f90

@@ -43,7 +43,8 @@ subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr, &
       type(TMolecule), intent(in) :: mol
       type(TNeigh), intent(inout) :: neigh ! contains nb, nbf and nbm
       logical, intent(in) :: makeneighbor, nb_call
-      integer at(natoms),natoms
+      integer, intent(in) :: natoms
+      integer at(natoms)
       integer hyb (natoms)
       integer itag(natoms)
       real*8  rab   (natoms*(natoms+1)/2)
@@ -180,7 +181,7 @@ subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr, &
             elseif(nm.eq.1)then  ! distinguish M-CR2-R i.e. not an eta coord.
               ncm=0
               do iTr=1, numctr
-                do k=1,numnbf  !
+                do k=1,neigh%nbf(neigh%numnb,i,iTr)  !
                   if(neigh%nbf(k,i,iTr).ne.im)then ! all neighbors that are not the metal im
                     kk=neigh%nbf(k,i,iTr)
                     do l=1,sum(neigh%nbf(neigh%numnb,kk,:))

src/gfnff/neighbor.f90

@@ -115,10 +115,10 @@ subroutine get_nb(self, mol, rab, rtmp, mchar, icase, f_in, f2_in, param)
       !integer, intent(inout) :: nbf(20,mol%n)
       type(TGFFData), intent(in) :: param
       !real(wp) :: latThresh, maxNBr
-      real(wp) :: dist(mol%n,mol%n,self%numctr)
+      real(wp), allocatable :: dist(:,:,:)
       integer :: i,j,iTr
 
-      dist=0.0_wp
+      allocate(dist(mol%n,mol%n,self%numctr), source=0.0_wp)
       !$omp parallel do collapse(3) default(none) shared(dist,mol,self) &
       !$omp private(iTr,i,j)
       do iTr=1, self%numctr
@@ -404,14 +404,15 @@ end function id2v
     subroutine fillnb(self,n,at,rad,dist,mchar,icase,f,f2,param)
       implicit none
       type(TGFFData), intent(in) :: param
-      class(TNeigh), intent(inout) :: self
+      type(TNeigh), intent(inout) :: self
       integer, intent(in) :: n,at(n)
       real(wp), intent(in) :: rad(n*(n+1)/2)
       real(wp), intent(in) :: dist(n,n,self%numctr)
       real(wp), intent(in) :: mchar(n),f,f2
       integer i,j,k,iTr,nn,icase,hc_crit,nnfi,nnfj,lin
-      integer tag(n,n,self%numctr)
+      integer, allocatable :: tag(:,:,:)
       real(wp) rco,fm
+      allocate(tag(n,n,self%numctr), source=0)
       if(icase.eq.1) then
         if(.not. allocated(self%nbf)) then
           allocate(self%nbf(self%numnb,n,self%numctr),source=0)
@@ -434,7 +435,6 @@ subroutine fillnb(self,n,at,rad,dist,mchar,icase,f,f2,param)
         endif
       endif
 
-      tag = 0
       do iTr=1, self%numctr
         do i=1,n
           do j=1,n

src/lbfgs_anc/driver.f90

@@ -36,7 +36,7 @@ module xtb_pbc_optimizer_driver
 
 !> Driver for performing geometry optimization
 !subroutine    relax(self, ctx, mol, calc, filter, accuracy, verbosity)
-subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel)
+subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel, optfail)
   use xtb_pbc, only: cross_product
   use xtb_gfnff_calculator, only : TGFFCalculator, newGFFCalculator
   !> Instance of the optimization driver
@@ -64,6 +64,8 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel)
   type(scc_results) :: results
   !> optimization cycle failed
   logical :: fail
+  !> optimization failed
+  logical, intent(out) :: optfail
   !> singlepoint energy
   real(wp) :: energy
   !> gradient and stress tensor
@@ -84,6 +86,7 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel)
   real(wp) :: emin_global, xyz_emin(3,mol%n), latt_emin(3,3)
   character(6), parameter :: fname = "noName"
 
+  optfail = .false.
   nvar = filter%get_dimension()
   allocate(gcurr(nvar), glast(nvar), displ(nvar))
   allocate(gxyz(3, mol%n), sigma(3, 3), source=0.0_wp)
@@ -203,6 +206,7 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel)
 
   ! Output message if not converged after maximum steps
   if (.not.converged) then
+    optfail = .true.
     write(*,*) ''
     write(*,*) 'Could not converge in',step,' steps.'
     call env%error("Could not converge geometry.",source)

src/set_module.f90

@@ -1508,16 +1508,18 @@ subroutine set_ffnb(env,key,val)
    character(len=*),intent(in) :: val
    integer :: i,j,k,l
    integer :: i_start,i_end,i_ffnb
+   logical :: nonb
 
    k=1 ! start at 1 since first entry of ffnb is the atom index that the following NBs belong to
    i_start=1
    i_end=1
    i_ffnb=0
+   nonb = .false. ! expect that the atom should have neighbors
    do i=1, len(val)
       ! get next empty row in ffnb and read atom index into first entry
       if (val(i:i).eq.":") then
          do j=1,size(set%ffnb, dim=2)
-           if (set%ffnb(1,j).eq.0 .and. i_ffnb.eq.0) then
+           if (set%ffnb(1,j).eq.-1 .and. i_ffnb.eq.0) then
               i_ffnb = j ! index of next empty row
               l=i-1
               ! we take care of ":" and "," and trust read() to handle whitespaces
@@ -1531,17 +1533,30 @@ subroutine set_ffnb(env,key,val)
       if (val(i:i).eq.",") then
          k = k + 1
          i_end=i-1
-             read(val(i_start:i_end), *) set%ffnb(k,i_ffnb)
-             i_start=i+1
+         read(val(i_start:i_end), *) set%ffnb(k,i_ffnb)
+         ! if the first neighbor index (k=2) is zero the atom (k=1) has no neighbors
+         if (set%ffnb(k,i_ffnb) == 0 .and. k == 2) then
+           nonb = .true.
+         endif
+         i_start=i+1
       endif
       ! read the last neighbor into ffnb
       if (i.eq.len(val)) then
          k = k + 1
          read(val(i_start:), *) set%ffnb(k,i_ffnb)
+         ! if the first neighbor index (k=2) is zero the atom (k=1) has no neighbors
+         if (set%ffnb(k,i_ffnb) == 0 .and. k == 2) then
+           nonb = .true.
+         endif
+         ! set rest of ffnb to 0
+         set%ffnb(k+1:41, i_ffnb) = 0
       endif
    enddo
-   set%ffnb(42,i_ffnb) = k - 1 ! number of neighbors of atom set%ffnb(1,i_ffnb) 
-
+   if (nonb) then
+     set%ffnb(2:,i_ffnb) = 0
+   else
+     set%ffnb(42,i_ffnb) = k - 1 ! number of neighbors of atom set%ffnb(1,i_ffnb) 
+   endif
 
 end subroutine set_ffnb
 

src/setparam.f90

@@ -519,7 +519,7 @@ module xtb_setparam
    !> GFN-FF manual setup of nb list via xcontrol
    !  allows a maximum of 164 atoms neighbors to be changed
    !  ffnb(42,i) stores the number of neighbors of atom i
-   integer :: ffnb(42,164) = 0
+   integer :: ffnb(42,164) = -1
    end type TSet
 
    type(TSet) :: set