Allow scaling of dispersion energy in GFN-FF (#359)

- add $gfn/dispscale=<real> keyword
- store dispersion scaling in calc%param%dispscale
- pass dispersion scaling to d3_gradient routine of GFN-FF
- also remove troublesome PGI preprocessor from hessian (again)

man/xcontrol.7.adoc

@@ -224,6 +224,9 @@ $gfn
 *method*='int'::
     version of the GFN Hamiltonian
 
+*dispscale*='real'::
+    Scale dispersion energy of GFN-FF
+
 $hess
 ~~~~~
 *sccacc*='real'::

src/gfnff/data.f90

@@ -87,6 +87,9 @@ module xtb_gfnff_data
       !> max CN cut-off
       real(wp) :: cnmax
 
+      !> D3 scaling
+      real(wp) :: dispscale
+
       !> Constant data
       real(wp), allocatable :: en(:)
       real(wp), allocatable :: rad(:)

src/gfnff/gdisp0.f90

@@ -237,7 +237,7 @@ subroutine get_atomic_c6_d4(dispm, nat, atoms, gwvec, gwdcn, c6, dc6dcn)
 end subroutine get_atomic_c6_d4
 
 subroutine d3_gradient(dispm, nat, at, xyz, npair, pairlist, zeta_scale, radii, &
-      & r4r2, weighting_factor, cn, dcndr, energy, gradient)
+      & r4r2, weighting_factor, dispscale, cn, dcndr, energy, gradient)
    use xtb_disp_dftd3param
    use xtb_type_dispersionmodel, only : TDispersionModel
 
@@ -251,6 +251,7 @@ subroutine d3_gradient(dispm, nat, at, xyz, npair, pairlist, zeta_scale, radii,
    real(wp), intent(in) :: r4r2(:)
    real(wp), intent(in) :: radii(:)
 
+   real(wp), intent(in) :: dispscale
    real(wp), intent(in) :: weighting_factor
    real(wp), intent(in) :: cn(:)
    real(wp), intent(in) :: dcndr(:, :, :)
@@ -281,7 +282,8 @@ subroutine d3_gradient(dispm, nat, at, xyz, npair, pairlist, zeta_scale, radii,
 
    !$omp parallel do default(none) schedule(runtime) &
    !$omp reduction(+:energies, gradient, sigma, dEdcn) &
-   !$omp shared(at, xyz, npair, pairlist, zeta_scale, radii, r4r2, c6, dc6dcn) &
+   !$omp shared(at, xyz, npair, pairlist, zeta_scale, dispscale, radii, r4r2, c6, &
+   !$omp&       dc6dcn) &
    !$omp private(ij, img, iat, jat, ati, atj, r2, rij, r4r2ij, r0, t6, t8, t10, &
    !$omp&        d6, d8, d10, disp, ddisp, dE, dG, dS)
    do img = 1, npair
@@ -302,8 +304,8 @@ subroutine d3_gradient(dispm, nat, at, xyz, npair, pairlist, zeta_scale, radii,
       d6 = -6*r2**2*t6**2
       d8 = -8*r2**3*t8**2
 
-      disp = (t6 + 2*r4r2ij*t8) * zeta_scale(ij)
-      ddisp= (d6 + 2*r4r2ij*d8) * zeta_scale(ij)
+      disp = (t6 + 2*r4r2ij*t8) * zeta_scale(ij)*dispscale
+      ddisp= (d6 + 2*r4r2ij*d8) * zeta_scale(ij)*dispscale
 
       dE = -c6(iat, jat)*disp * 0.5_wp
       dG = -c6(iat, jat)*ddisp*rij

src/gfnff/gfnff_eg.f90

@@ -255,7 +255,7 @@ subroutine gfnff_eg(env,pr,n,ichrg,at,xyz,makeq,g,etot,res_gff, &
       if (pr) call timer%measure(5,'D3')
       if(nd3.gt.0) then
          call d3_gradient(topo%dispm, n, at, xyz, nd3, d3list, topo%zetac6, &
-            & param%d3r0, sqrtZr4r2, 4.0d0, cn, dcn, edisp, g)
+            & param%d3r0, sqrtZr4r2, 4.0d0, param%dispscale, cn, dcn, edisp, g)
       endif
       deallocate(d3list)
       if (pr) call timer%measure(5)

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
+  use xtb_setparam, only : ichrg, dispscale
   implicit none
   character(len=*), parameter :: source = 'gfnff_setup'
 ! Dummy
@@ -58,6 +58,8 @@ subroutine gfnff_setup(env,verbose,restart,mol,gen,param,topo,accuracy,version)
   end if
 
   call gfnff_set_param(mol%n, gen, param)
+  param%dispscale = dispscale
+  print*, dispscale
   if (restart) then
      inquire(file='gfnff_topo', exist=ex)
      if (ex) then

src/hessian.F90

@@ -204,7 +204,6 @@ subroutine numhess( &
       ! ascending order
       ! now compute a subblock of the Hessian
 ! PGI 20.7 produces invalid LLVM-IR with the following OpenMP directives
-#ifndef __PGIC__
       !$ nproc = omp_get_num_threads()
       !$omp parallel default(shared) &
       !$omp firstprivate(et,maxiter,acc) &
@@ -215,7 +214,6 @@ subroutine numhess( &
       !$ call mkl_set_num_threads(1)
 #endif
       !$omp do schedule(dynamic)
-#endif
       do a = 1, nonfrozh
          ia=indx(a)
          do ic = 1, 3
@@ -256,21 +254,18 @@ subroutine numhess( &
             !$omp end critical
          endif
       enddo
-#ifndef __PGIC__
       !$omp end do
       !$omp end parallel
       !$ call omp_set_num_threads(nproc)
 #ifdef WITH_MKL
       !$ call mkl_set_num_threads(nproc)
-#endif
 #endif
 
    else
 !! ------------------------------------------------------------------------
 !  normal case
 !! ------------------------------------------------------------------------
 ! PGI 20.7 produces invalid LLVM-IR with the following OpenMP directives
-#ifndef __PGIC__
       !$ nproc = omp_get_num_threads()
       !$omp parallel default(shared) &
       !$omp firstprivate(et,maxiter,acc) &
@@ -281,7 +276,6 @@ subroutine numhess( &
       !$ call mkl_set_num_threads(1)
 #endif
       !$omp do schedule(dynamic)
-#endif
       do ia = 1, mol%n
          do ic = 1, 3
             ii = (ia-1)*3+ic
@@ -331,13 +325,11 @@ subroutine numhess( &
             !$omp end critical
          endif
       enddo
-#ifndef __PGIC__
       !$omp end do
       !$omp end parallel
       !$ call omp_set_num_threads(nproc)
 #ifdef WITH_MKL
       !$ call mkl_set_num_threads(nproc)
-#endif
 #endif
 
    endif

src/set_module.f90

@@ -169,6 +169,7 @@ subroutine write_set_gfn(ictrl)
       write(ictrl,'(3x,"d4=",a)') bool2string(newdisp)
    write(ictrl,'(3x,"scc=",a)') bool2string(solve_scc)
    write(ictrl,'(3x,"periodic=",a)') bool2string(periodic)
+   write(ictrl,'(3x,"dispscale=",g0)') dispscale
 end subroutine write_set_gfn
 
 subroutine write_set_scc(ictrl)
@@ -1322,6 +1323,7 @@ subroutine set_gfn(env,key,val)
    logical,save :: set2 = .true.
    logical,save :: set3 = .true.
    logical,save :: set4 = .true.
+   logical,save :: set5 = .true.
    select case(key)
    case default ! do nothing
       call env%warning("the key '"//key//"' is not recognized by gfn",source)
@@ -1352,6 +1354,9 @@ subroutine set_gfn(env,key,val)
    case('periodic')
       if (getValue(env,val,ldum).and.set4) periodic = ldum
       set4 = .false.
+   case('dispscale')
+      if (getValue(env,val,ddum).and.set5) dispscale = ddum
+      set5 = .false.
    end select
 end subroutine set_gfn
 

src/setparam.f90

@@ -70,6 +70,7 @@ module xtb_setparam
    real(wp) :: etemp = 300.0_wp
 !  damping in Broyden SCC procedure (0.05 for critical cases, autoadjusted)
    real(wp) :: broydamp = 0.40_wp
+   real(wp) :: dispscale = 1.0_wp
 
    integer, parameter :: p_guess_sad = 0
    integer, parameter :: p_guess_gasteiger = 1