Implement dipole moment in GFN-FF (#869)

* uses EEQ charges and atom positions

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

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)
 

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