Update handling of charges in GFN-FF (#558)

- only use PDB charges from residues if they match system charge
  (previous behaviour: PDB charges always overwrite total system charge)

src/gfnff/gfnff_ini.f90

@@ -17,7 +17,7 @@
 module xtb_gfnff_ini
 contains
 
-subroutine gfnff_ini(env,pr,makeneighbor,mol,ichrg,gen,param,topo,accuracy)
+subroutine gfnff_ini(env,pr,makeneighbor,mol,gen,param,topo,accuracy)
       use xtb_mctc_accuracy, only : wp, sp
       use xtb_type_molecule
       use xtb_type_environment, only : TEnvironment
@@ -41,7 +41,6 @@ subroutine gfnff_ini(env,pr,makeneighbor,mol,ichrg,gen,param,topo,accuracy)
       type(TGFFData), intent(in) :: param
       real(wp), intent(in) :: accuracy
 
-      integer, intent(in) :: ichrg         ! mol. charge
       logical, intent(in) :: pr            ! print flag
       logical, intent(in) :: makeneighbor  ! make a neigbor list or use existing one?
 !--------------------------------------------------------------------------------------------------
@@ -445,7 +444,7 @@ integer function itabrow6(i)
          if (is_iostat_end(err) .and. i == mol%n) err = 0
          call close_file(ich)
          if (err == 0) then
-            if (i < mol%n .or. abs(sum(qtmp) - ichrg) > 1.0e-3_wp) then
+            if (i < mol%n .or. abs(sum(qtmp) - mol%chrg) > 1.0e-3_wp) then
                call env%warning("Rejecting external charges input due to missmatch", source)
             else
                topo%qfrag=dnint(qtmp)
@@ -460,7 +459,7 @@ integer function itabrow6(i)
             return
          end if
       endif
-      if(mol%n.lt.100.and.topo%nfrag.gt.2.and.ichrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
+      if(mol%n.lt.100.and.topo%nfrag.gt.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
          itmp=0
          do i=1,mol%n
             itmp(topo%fraglist(i))=itmp(topo%fraglist(i))+1
@@ -471,22 +470,22 @@ integer function itabrow6(i)
          call env%error('fragment charge input required', source)
          return
       endif
-      if(mol%n.ge.100.and.topo%nfrag.gt.2.and.ichrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
-         topo%qfrag(1)=ichrg
+      if(mol%n.ge.100.and.topo%nfrag.gt.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
+         topo%qfrag(1)=mol%chrg
          topo%qfrag(2:topo%nfrag)=0
       endif
-      if(topo%nfrag.eq.2.and.ichrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
+      if(topo%nfrag.eq.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
          write(env%unit,*) 'trying auto detection of charge on 2 fragments:'
          topo%qfrag(1)=0
-         topo%qfrag(2)=dble(ichrg)
+         topo%qfrag(2)=mol%chrg
          call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,dum1,topo)
          call env%check(exitRun)
          if (exitRun) then
             call env%error("Failed to generate charges", source)
             return
          end if
          topo%qfrag(2)=0
-         topo%qfrag(1)=dble(ichrg)
+         topo%qfrag(1)=mol%chrg
          call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,dum2,topo)
          call env%check(exitRun)
          if (exitRun) then
@@ -495,7 +494,7 @@ integer function itabrow6(i)
          end if
          if(dum1.lt.dum2) then
             topo%qfrag(1)=0
-            topo%qfrag(2)=dble(ichrg)
+            topo%qfrag(2)=mol%chrg
          endif
          write(env%unit,*) 'dEes      :',dum1-dum2
          write(env%unit,*) 'charge 1/2:',topo%qfrag(1:2)

src/gfnff/gfnff_setup.f90

@@ -31,7 +31,7 @@ subroutine gfnff_setup(env,verbose,restart,mol,gen,param,topo,accuracy,version)
   use xtb_type_environment, only : TEnvironment
   use xtb_type_molecule, only : TMolecule
   use xtb_gfnff_param, only : ini, gfnff_set_param
-  use xtb_setparam, only : ichrg, dispscale
+  use xtb_setparam, only : dispscale
   implicit none
   character(len=*), parameter :: source = 'gfnff_setup'
 ! Dummy
@@ -77,7 +77,7 @@ subroutine gfnff_setup(env,verbose,restart,mol,gen,param,topo,accuracy,version)
      end if
   end if
 
-  call gfnff_ini(env,verbose,ini,mol,ichrg,gen,param,topo,accuracy)
+  call gfnff_ini(env,verbose,ini,mol,gen,param,topo,accuracy)
 
   call env%check(exitRun)
   if (exitRun) then
@@ -98,7 +98,6 @@ subroutine gfnff_input(env, mol, topo)
   use xtb_type_molecule
   use xtb_mctc_filetypes, only : fileType
   use xtb_gfnff_param
-  use xtb_setparam, only : ichrg
   implicit none
   ! Dummy
   type(TMolecule),intent(in) :: mol
@@ -149,8 +148,14 @@ subroutine gfnff_input(env, mol, topo)
         & + dble(mol%pdb(iatom)%charge)
       topo%qpdb(iatom) = mol%pdb(iatom)%charge
     end do
-    ichrg=idint(sum(topo%qfrag(1:topo%nfrag)))
-    write(env%unit,'(10x,"charge from pdb residues: ",i0)') ichrg
+    if (abs(mol%chrg - sum(topo%qfrag(1:topo%nfrag))) < sqrt(epsilon(1.0_wp))) then
+       write(env%unit,'(10x,"charge from pdb residues: ",i0)') &
+          & nint(sum(topo%qfrag(1:topo%nfrag)))
+    else
+       if (allocated(topo%qpdb)) deallocate(topo%qpdb)
+       topo%qfrag(1:topo%nfrag) = 0.0_wp
+       topo%nfrag = 0
+    end if
   !--------------------------------------------------------------------
   ! SDF case
   case(fileType%sdf,fileType%molfile)
@@ -203,12 +208,18 @@ subroutine gfnff_input(env, mol, topo)
     ini = .true.
     call open_file(ich,'.CHRG','r')
     if (ich.ne.-1) then
-      read(ich,'(a)')atmp
+      read(ich,'(a)')  ! first line contains total charge
+      read(ich,'(a)')atmp  ! second line contains fragment charges
       call close_file(ich)
       call readline(atmp,floats,s,ns,nf)
       topo%qfrag(1:nf)=floats(1:nf)
-      ichrg=int(sum(topo%qfrag(1:nf)))
-      topo%qfrag(nf+1:mol%n)=9999
+      if (abs(mol%chrg - sum(topo%qfrag(1:nf))) < sqrt(epsilon(1.0_wp))) then
+         write(env%unit,'(10x,"charge from .CHRG file: ",i0)') &
+            & nint(sum(topo%qfrag(1:nf)))
+      else
+         ! ignore fragment charges if they are not consistent with the total charge
+         topo%qfrag(1:nf) = 0.0_wp
+      end if
     else
       topo%qfrag(1)=mol%chrg
       topo%qfrag(2:mol%n)=0