From 8354d4060c3f60ad9ead03003cab1bc2a9648569 Mon Sep 17 00:00:00 2001
From: Thomas Rose <39367840+Thomas3R@users.noreply.github.com>
Date: Fri, 24 Mar 2023 16:13:00 +0100
Subject: [PATCH 1/3] Add condition for new torsion potential

---
 src/gfnff/gfnff_ini.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/gfnff/gfnff_ini.f90 b/src/gfnff/gfnff_ini.f90
index 6649f46bb..41f605f67 100644
--- a/src/gfnff/gfnff_ini.f90
+++ b/src/gfnff/gfnff_ini.f90
@@ -1912,12 +1912,12 @@ subroutine specialTorsList(nst, mol, topo, sTorsList)
             ! at this point we know that i and nbi are carbons bonded through triple bond
             ! check C2 and C3
             do k=1, 2  ! C2 is other nb of Ci
-              if (topo%nb(k,i).ne.nbi) then
+              if (topo%nb(k,i).ne.nbi.and.mol%at(topo%nb(k,i)).eq.6) then
                 jj=topo%nb(k,i)
               endif
             enddo
             do k=1, 2  ! C3 is other nb of Cnbi
-              if (topo%nb(k,nbi).ne.i) then
+              if (topo%nb(k,nbi).ne.i.and.mol%at(topo%nb(k,nbi)).eq.6) then
                 kk=topo%nb(k,nbi)
               endif
             enddo

From 1ea8e52fbb16a557153b5f99c169976590bfbbd8 Mon Sep 17 00:00:00 2001
From: Thomas Rose <39367840+Thomas3R@users.noreply.github.com>
Date: Thu, 24 Aug 2023 14:07:57 +0200
Subject: [PATCH 2/3] Fix #849: Initialize sTorsl indices

---
 src/gfnff/gfnff_ini.f90 | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/gfnff/gfnff_ini.f90 b/src/gfnff/gfnff_ini.f90
index d88b3e6b0..35df19e35 100644
--- a/src/gfnff/gfnff_ini.f90
+++ b/src/gfnff/gfnff_ini.f90
@@ -1910,7 +1910,13 @@ subroutine specialTorsList(nst, mol, topo, sTorsList)
   integer, intent(inout) :: sTorsList(6, nst)
   integer :: i,j,k,ii,jj,kk,ll,idx
   logical :: iiok, llok
+  ! initialize variables
   idx=0
+  ii=-1
+  jj=-1
+  kk=-1
+  ll=-1
+
   do i=1, mol%n
     ! carbon with two neighbors bonded to other carbon* with two neighbors
     if (mol%at(i).eq.6.and.topo%nb(20,i).eq.2) then
@@ -1931,6 +1937,9 @@ subroutine specialTorsList(nst, mol, topo, sTorsList)
                 kk=topo%nb(k,nbi)
               endif
             enddo
+            if (jj.eq.-1.or.kk.eq.-1) then
+              exit ! next atom i
+            endif
             ! check C1 through C4 are sp2 carbon
             if (topo%hyb(jj).eq.2.and.topo%hyb(kk).eq.2 &
             &   .and.mol%at(jj).eq.6.and.mol%at(kk).eq.6) then

From b6833016538e3c94f8bf0161c4fb7c2029054ab5 Mon Sep 17 00:00:00 2001
From: Thomas Rose <39367840+Thomas3R@users.noreply.github.com>
Date: Tue, 19 Sep 2023 18:18:05 +0200
Subject: [PATCH 3/3] Add calculation of dipole moment for GFN-FF

---
 src/aespot.f90        | 28 ++++++++++++++++++++++++++++
 src/main/property.F90 | 19 +++++++++++++++++++
 src/prog/main.F90     |  2 ++
 3 files changed, 49 insertions(+)

diff --git a/src/aespot.f90 b/src/aespot.f90
index 86ff62c03..c10d90b67 100644
--- a/src/aespot.f90
+++ b/src/aespot.f90
@@ -741,6 +741,34 @@ subroutine molmom(iunit,n,xyz,q,dipm,qp,dip,d3)
 
 end subroutine molmom
 
+! molqdip: computes molecular dipole moments from charge only
+! n           : # of atoms
+! xyz(3,n)    : cartesian coordinates
+! q(n)        : atomic partial charges
+subroutine molqdip(iunit,n,xyz,q)
+   use xtb_mctc_convert
+   implicit none
+   integer, intent(in) :: iunit
+   integer, intent(in) :: n
+   real(wp), intent(in)  :: xyz(:,:),q(:)
+   real(wp) rr1(3), dip
+   integer i,j
+   rr1 = 0.0_wp
+   write(iunit,'(a)')
+   do i = 1,n
+      do j = 1,3
+         rr1(j) = rr1(j)+q(i)*xyz(j,i)
+      enddo
+   enddo
+   dip = sqrt(rr1(1)**2+rr1(2)**2+rr1(3)**2)
+   write(iunit,'(a)',advance='yes')'molecular dipole:'
+   write(iunit,'(a)',advance='no')'                 '
+   write(iunit,'(a)',advance='yes') &
+      & 'x           y           z       tot (Debye)'
+   write(iunit,'(a,4f12.3)') ' q only: ',rr1(1:3),dip*autod
+
+end subroutine molqdip
+
 ! gradient evaluation from
 ! cumulative atomic multipole moment interactions: all interactions up to r**-3
 ! nat         : # of atoms
diff --git a/src/main/property.F90 b/src/main/property.F90
index dcf9f265f..bbeaca3be 100644
--- a/src/main/property.F90
+++ b/src/main/property.F90
@@ -250,6 +250,25 @@ subroutine main_property &
 
 end subroutine main_property
 
+subroutine gfnff_property(iunit, n, xyz, topo, nlist)
+  use xtb_gfnff_topology, only: TGFFTopology 
+  use xtb_gfnff_neighbourlist, only: TGFFNeighbourList
+  use xtb_aespot, only: molqdip
+  !! ========================================================================
+  !  global storage of options, parameters and basis set
+  use xtb_setparam
+  integer, intent(in) :: iunit, n
+  real(wp), intent(in) :: xyz(3,n)
+  type(TGFFTopology), intent(in) :: topo
+  type(TGFFNeighbourList), intent(in) :: nlist
+
+  ! dipole moment from charge
+   if (set%pr_dipole) then
+     call molqdip(iunit, n, xyz, nlist%q)
+   endif
+
+end subroutine gfnff_property
+
 subroutine main_cube &
       (lverbose,mol,wfx,basis,res)
 
diff --git a/src/prog/main.F90 b/src/prog/main.F90
index c319615be..447f93797 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -892,6 +892,8 @@ subroutine xtbMain(env, argParser)
          call main_property(iprop,env,mol,chk%wfn,calc%basis,calc%xtbData,res, &
             & calc%solvation,set%acc)
          call main_cube(set%verbose,mol,chk%wfn,calc%basis,res)
+      type is(TGFFCalculator)
+         call gfnff_property(iprop,mol%n,mol%xyz,calc%topo,chk%nlist)
       end select
    endif