GFN-FF energy with GBSA (#251)

* fixed charge input

Signed-off-by: Sebastian Spicher <[email protected]>

* using molecule type

Signed-off-by: Sebastian Spicher <[email protected]>

* GBSA total energy fixed + GFN-FF/GBSA Testsuite

Signed-off-by: Sebastian Spicher <[email protected]>

TESTSUITE/gfnff.f90

@@ -156,11 +156,103 @@ subroutine test_gfnff_hb
    call assert_close(res_gff%e_es,   -0.152313816530_wp,thr*10)
    call assert_close(res_gff%e_disp, -0.001251669186_wp,thr)
    call assert_close(res_gff%e_rep,   0.066881023899_wp,thr)
-   call assert_close(res_gff%e_hb,   -0.0068942923371_wp,thr)
-   call assert_close(res_gff%e_xb,   -0.0000000000000_wp,thr)
-   call assert_close(res_gff%e_batm, -0.0000000000000_wp,thr)
+   call assert_close(res_gff%e_hb,   -0.006894292337_wp,thr)
+   call assert_close(res_gff%e_xb,   -0.000000000000_wp,thr)
+   call assert_close(res_gff%e_batm, -0.000000000000_wp,thr)
 
    call mol%deallocate
 
    call terminate(afail)
 end subroutine test_gfnff_hb
+
+subroutine test_gfnff_gbsa
+   use xtb_mctc_accuracy, only : wp
+   use assertion
+   use xtb_mctc_systools
+   use xtb_type_environment
+   use xtb_type_options
+   use xtb_type_molecule
+   use xtb_type_data
+   use xtb_gfnff_param
+   use xtb_gfnff_setup
+   use xtb_gfnff_eg
+   use xtb_gfnff_ini
+   use xtb_setparam
+   use xtb_setmod
+   use xtb_disp_dftd3param
+   use xtb_disp_dftd4
+   use xtb_gfnff_calculator, only : TGFFCalculator
+   use xtb_main_setup, only : newGFFCalculator, addSolvationModel
+   implicit none
+   real(wp),parameter :: thr = 1.0e-10_wp
+   integer, parameter :: nat = 7
+   integer, parameter :: at(nat) = [6,8,1,1,8,1,1]
+   real(wp),parameter :: xyz(3,nat) = reshape(&
+      &[-5.78520874132429_wp,-1.92328475821000_wp,-0.02944611115854_wp, &
+      & -5.57801768832583_wp, 0.17912532844037_wp, 0.72444143178660_wp, &
+      & -4.27822256938673_wp,-2.74845397256109_wp,-1.13038073598642_wp, &
+      & -7.47879539136783_wp,-2.97570121473211_wp, 0.39488815557786_wp, &
+      & -0.83005338399036_wp, 2.43458470560665_wp,-0.78566331969245_wp, &
+      & -0.74201439536855_wp, 4.04199055249898_wp, 0.09144422329636_wp, &
+      & -2.44679415487233_wp, 1.69392751177087_wp,-0.27417668699116_wp],&
+      & shape(xyz))
+   type(scc_options),parameter :: opt = scc_options( &
+      &  prlevel = 2, maxiter = 30, acc = 1.0_wp, etemp = 300.0_wp, grad = .true.,&
+      &  solvent = "h2o")
+   logical, parameter :: restart = .false.
+
+   type(TMolecule)     :: mol
+   type(TEnvironment)  :: env
+   type(scc_results)   :: res_gff
+   type(TGFFCalculator) :: calc
+
+   real(wp) :: etot
+   real(wp), allocatable :: g(:,:)
+   character(len=:),allocatable :: fnv
+   integer  :: ipar
+
+   logical  :: exist
+
+   call init(env)
+   call init(mol,at,xyz)
+
+   call delete_file('charges')
+   call newGFFCalculator(env, mol, calc, '---', .false.)
+   call addSolvationModel(env, calc, opt%solvent)
+
+   call env%checkpoint("GFN-FF parameter setup failed")
+
+   allocate( g(3,mol%n), source = 0.0_wp )
+
+   call assert_eq(calc%topo%nbond,5)
+   call assert_eq(calc%topo%nangl,4)
+   call assert_eq(calc%topo%ntors,1)
+
+   g = 0.0_wp
+   gff_print=.true.
+
+   call gfnff_eg(env,gff_print,mol%n,nint(mol%chrg),mol%at,mol%xyz,make_chrg, &
+      & g,etot,res_gff,calc%param,calc%topo,calc%solv,.true.,calc%version, &
+      & calc%accuracy)
+
+   call assert_close(res_gff%e_total,-0.963965759735_wp,thr)
+   call assert_close(res_gff%gnorm,   0.006620747519_wp,thr)
+   call assert_close(res_gff%e_bond, -0.856707643513_wp,thr)
+   call assert_close(res_gff%e_angl,  0.000579711773_wp,thr)
+   call assert_close(res_gff%e_tors,  0.000000008811_wp,thr)
+   call assert_close(res_gff%e_es,   -0.151138740773_wp,thr*10)
+   call assert_close(res_gff%e_disp, -0.001251669186_wp,thr)
+   call assert_close(res_gff%e_rep,   0.066881023899_wp,thr)
+   call assert_close(res_gff%e_hb,   -0.006894292337_wp,thr)
+   call assert_close(res_gff%e_xb,   -0.000000000000_wp,thr)
+   call assert_close(res_gff%e_batm, -0.000000000000_wp,thr)
+   call assert_close(res_gff%g_solv, -0.015434158338_wp,thr)
+   call assert_close(res_gff%g_sasa,  0.000126368690_wp,thr)
+   call assert_close(res_gff%g_hb,   -0.008444397163_wp,thr)
+   call assert_close(res_gff%g_born, -0.008973572992_wp,thr)
+   call assert_close(res_gff%g_shift, 0.001857443126_wp,thr)
+
+   call mol%deallocate
+
+   call terminate(afail)
+end subroutine test_gfnff_gbsa

TESTSUITE/meson.build

@@ -174,3 +174,4 @@ test('GFN0-xTB: SRB (PBC)', xtb_test, args: ['peeq', 'srb'], env: xtbenv)
 
 test('GFN-FF: SP', xtb_test, args: ['gfnff', 'sp'], env: xtbenv)
 test('GFN-FF: SP (HB)', xtb_test, args: ['gfnff', 'hb'], env: xtbenv)
+test('GFN-FF: GBSA', xtb_test, args: ['gfnff', 'gbsa'], env: xtbenv)

TESTSUITE/tests_peeq.f90

@@ -96,6 +96,7 @@ program peeq_tester
       select case(sec)
       case('sp');  call test_gfnff_sp
       case('hb');  call test_gfnff_hb
+      case('gbsa');call test_gfnff_gbsa
       end select
    case('peeq')
       select case(sec)

src/gfnff/calculator.f90

@@ -166,6 +166,8 @@ subroutine singlepoint(self, env, mol, wfn, printlevel, restart, &
       write(env%unit,'(9x,53(":"))')
       write(env%unit,outfmt) "total energy      ", results%e_total,"Eh   "
       if (.not.silent.and.allocated(self%solv)) then
+         write(env%unit,outfmt) "total w/o Gsolv   ", &
+            &  results%e_total-results%g_solv, "Eh   "
          write(env%unit,outfmt) "total w/o Gsasa/hb", &
             &  results%e_total-results%g_sasa-results%g_hb-results%g_shift, "Eh   "
       endif

src/gfnff/gfnff_eg.f90

@@ -539,7 +539,7 @@ subroutine gfnff_eg(env,pr,n,ichrg,at,xyz,makeq,g,etot,res_gff, &
 !!!!!!!!!!!!!!!!!!
       etot = ees + edisp + erep + ebond &
      &           + eangl + etors + ehb + exb + ebatm + eext &
-     &           + gsolv + gshift + gborn + ghb
+     &           + gsolv
 
 !!!!!!!!!!!!!!!!!!
 ! printout