expose eeqbc standalone via calculator

Author: pprcht

src/calculator/api_engrad.f90

@@ -35,6 +35,7 @@ module api_engrad
   use gfn0_api
   use gfnff_api
   use crest_solvation,only:solvation_setup,solvation_core
+  use crest_electrostatic,only:electrostatic_core
   use libpvol_api
   use lj
   use approxg_module
@@ -54,6 +55,7 @@ module api_engrad
   public :: mlip_engrad
   public :: preinit_mlip_parallel
   public :: solvation_engrad
+  public :: electrostatic_engrad
 
 !=========================================================================================!
 !=========================================================================================!
@@ -157,6 +159,37 @@ subroutine solvation_engrad(mol,calc,energy,grad,iostatus)
     return
   end subroutine solvation_engrad
 
+!========================================================================================!
+
+  subroutine electrostatic_engrad(mol,calc,energy,grad,iostatus)
+!*********************************************************************
+!* Interface singlepoint call for the charge-equilibration calculator
+!*********************************************************************
+    implicit none
+    type(coord) :: mol
+    type(calculation_settings) :: calc
+    real(wp),intent(inout) :: energy
+    real(wp),intent(inout) :: grad(3,mol%nat)
+    integer,intent(out) :: iostatus
+    character(len=:),allocatable :: model
+
+    iostatus = 0
+    if (.not.allocated(calc%eeq)) then
+      iostatus = 1
+      return
+    end if
+    model = 'eeqbc'
+    if (allocated(calc%eeq%charge_model)) model = calc%eeq%charge_model
+
+!>--- energy, gradient and atomic charges (kept on the settings)
+    if (.not.allocated(calc%qat)) allocate (calc%qat(mol%nat))
+    call electrostatic_core(mol,calc%chrg,model,energy,grad,calc%qat,iostatus)
+    if (iostatus /= 0) return
+    if (.not.calc%prstdout) call api_print_e_grd(.false.,calc%prch,mol,energy,grad)
+
+    return
+  end subroutine electrostatic_engrad
+
 !========================================================================================!
 
   subroutine gfn0_engrad(mol,calc,g0calc,energy,grad,iostatus)

src/calculator/calc_type.f90

@@ -25,6 +25,7 @@ module calc_type
   use tblite_api
   use gfn0_api
   use crest_solvation,only:solvation_data
+  use crest_electrostatic,only:electrostatic_data
   use gfnff_api,only:gfnff_data
   use libpvol_api,only:libpvol_calculator
 !>--- other types
@@ -58,10 +59,11 @@ module calc_type
     integer :: penalty   = 13
     integer :: mlip      = 14
     integer :: solvation = 15
+    integer :: electrostatic = 16
   end type enum_jobtype
   type(enum_jobtype), parameter,public :: jobtype = enum_jobtype()
 
-  character(len=45),parameter,private :: jobdescription(16) = [ &
+  character(len=45),parameter,private :: jobdescription(17) = [ &
      & 'Unknown calculation type                    ', &
      & 'xTB calculation via external binary         ', &
      & 'Generic script execution                    ', &
@@ -77,7 +79,8 @@ module calc_type
      & 'Approximate free energy computation         ', &
      & 'Empirical penalty function                  ', &
      & 'MLIP via persistent python socket           ', &
-     & 'Standalone implicit solvation contribution  ']
+     & 'Standalone implicit solvation contribution  ', &
+     & 'Charge-equilibration electrostatics         ']
 !&>
 
 !=========================================================================================!
@@ -188,8 +191,11 @@ module calc_type
     real(wp) :: pvradscal = 1.0_wp     !> Scaling factor for SAS radii
     type(libpvol_calculator),allocatable :: libpvol
 
-!>--- implicit solvation composite data      
-    type(solvation_data),allocatable :: solv 
+!>--- implicit solvation composite data
+    type(solvation_data),allocatable :: solv
+
+!>--- charge-equilibration electrostatics data
+    type(electrostatic_data),allocatable :: eeq
 
 !>--- approxg data
     type(approxg_params) :: ag
@@ -1164,6 +1170,7 @@ subroutine calculation_settings_deallocate(self)
     if (allocated(self%tblite)) deallocate (self%tblite)
     if (allocated(self%g0calc)) deallocate (self%g0calc)
     if (allocated(self%solv)) deallocate (self%solv)
+    if (allocated(self%eeq)) deallocate (self%eeq)
     if (allocated(self%ff_dat)) deallocate (self%ff_dat)
     if (allocated(self%libpvol)) deallocate (self%libpvol)
 
@@ -1240,6 +1247,7 @@ subroutine calculation_settings_copy(self,src)
     if (allocated(src%refcharges))     self%refcharges     = src%refcharges
     if (allocated(src%tbliteparam))    self%tbliteparam    = src%tbliteparam
     if (allocated(src%solv))           self%solv           = src%solv
+    if (allocated(src%eeq))            self%eeq            = src%eeq
 
 ! ── gradient settings ────────────────────────────────────────────────────────
     self%numgrad    = src%numgrad
@@ -1452,6 +1460,12 @@ subroutine calculation_settings_shortflag(self)
         if (allocated(self%solv%solvent)) &
         & self%shortflag = self%shortflag//'('//trim(self%solv%solvent)//')'
       end if
+    case (jobtype%electrostatic)
+      self%shortflag = 'EEQ'
+      if (allocated(self%eeq)) then
+        if (allocated(self%eeq%charge_model)) &
+        & self%shortflag = trim(self%eeq%charge_model)
+      end if
     case default
       self%shortflag = 'undefined'
     end select
@@ -1631,6 +1645,15 @@ subroutine calculation_settings_info(self,iunit)
       end if
     end if
 
+    !> charge-equilibration electrostatics details
+    if (self%id == jobtype%electrostatic .and. allocated(self%eeq)) then
+      write (iunit,fmt4) 'Charge-equilibration electrostatics (multicharge)'
+      if (allocated(self%eeq%charge_model)) then
+        write (atmp,*) 'charge model'
+        write (iunit,fmt3) atmp,trim(self%eeq%charge_model)
+      end if
+    end if
+
     if (any((/jobtype%orca,jobtype%xtbsys,jobtype%turbomole, &
     &  jobtype%generic,jobtype%terachem/) == self%id)) then
       if (index(self%binary,'gxtb') .ne. 0) then

src/calculator/calculator.F90

@@ -411,6 +411,9 @@ subroutine potential_core(molptr,calc,id,iostatus)
 
     case (jobtype%solvation)
       call solvation_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)
+
+    case (jobtype%electrostatic)
+      call electrostatic_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)
     case default
       calc%etmp(id) = 0.0_wp
       calc%grdtmp(:,:,id) = 0.0_wp

src/calculator/components/electrostatic.F90

@@ -38,12 +38,18 @@ module crest_electrostatic
   implicit none
   private
 
-  public :: electrostatic_engrad
+  public :: electrostatic_data
+  public :: electrostatic_core
+
+!> Bundled electrostatic settings (stored on calculation_settings)
+  type :: electrostatic_data
+    character(len=:),allocatable :: charge_model  !> 'eeq' | 'eeqbc'
+  end type electrostatic_data
 
 contains
 !========================================================================================!
 
-  subroutine electrostatic_engrad(mol,chrg,model,energy,gradient,qat,iostatus,dqdr)
+  subroutine electrostatic_core(mol,chrg,model,energy,gradient,qat,iostatus,dqdr)
     !***********************************************************************
     !* Charge-equilibration energy, gradient and atomic charges.
     !*
@@ -121,7 +127,7 @@ subroutine electrostatic_engrad(mol,chrg,model,energy,gradient,qat,iostatus,dqdr
     qat(:) = 0.0_wp
     if (present(dqdr)) dqdr(:,:,:) = 0.0_wp
 #endif
-  end subroutine electrostatic_engrad
+  end subroutine electrostatic_core
 
 !========================================================================================!
 end module crest_electrostatic

src/calculator/components/solvation.F90

@@ -37,7 +37,7 @@ module crest_solvation
   use crest_parameters
   use strucrd,only:coord
 #ifdef WITH_DDX
-  use crest_electrostatic,only:electrostatic_engrad
+  use crest_electrostatic,only:electrostatic_core
   use crest_ddx_pc,only:ddx_pc_engrad
   use crest_surface,only:surface_engrad
   use mctc_env,only:error_type
@@ -149,7 +149,7 @@ subroutine solvation_core(mol,chrg,solv,energy,gradient,iostatus)
     allocate (qat(nat),dqdr(3,nat,nat),gtmp(3,nat),sasa(nat),dsdr(3,nat,nat))
 
     ! ── charges + dq/dR from the electrostatic component ─────────────────────
-    call electrostatic_engrad(mol,chrg,solv%charge_model,edum,gtmp,qat,iostatus,dqdr=dqdr)
+    call electrostatic_core(mol,chrg,solv%charge_model,edum,gtmp,qat,iostatus,dqdr=dqdr)
     if (iostatus /= 0) return
 
     ! ── polar (ddX) free energy + gradient (incl. dq/dR chain term) ──────────

src/parsing/parse_calcdata.f90

@@ -290,6 +290,10 @@ subroutine parse_setting_auto(env,job,kv,rd)
       case ('solvation','solvator')
         job%id = jobtype%solvation
         if (.not.allocated(job%solv)) allocate (job%solv)
+      case ('eeq','eeqbc','electrostatic')
+        job%id = jobtype%electrostatic
+        if (.not.allocated(job%eeq)) allocate (job%eeq)
+        if (kv%value_c /= 'electrostatic') job%eeq%charge_model = kv%value_c
       case default
         job%id = jobtype%unknown
         !>--- keyword was recognized, but invalid argument supplied
@@ -414,6 +418,11 @@ subroutine parse_setting_auto(env,job,kv,rd)
       if (.not.allocated(job%solv)) allocate (job%solv)
       job%solv%do_hbond = kv%value_b
 
+!>--- charge-equilibration electrostatics calculator settings
+    case ('eeq_model','charge_model')
+      if (.not.allocated(job%eeq)) allocate (job%eeq)
+      job%eeq%charge_model = kv%value_c
+
     case ('refine','refinement')
       select case (kv%value_c)
       case ('sp','singlepoint')

test/test_ddx.F90

@@ -5,7 +5,7 @@ module test_ddx
   use crest_testmol
 #ifdef WITH_DDX
   use crest_ddx_pc,only:ddx_pc_engrad
-  use crest_electrostatic,only:electrostatic_engrad
+  use crest_electrostatic,only:electrostatic_core
   use crest_surface,only:surface_engrad
   use crest_solvation,only:solvation_data,solvation_setup,solvation_core
 #endif
@@ -132,16 +132,16 @@ subroutine test_eeq_fd(error)
     h = 1.0e-4_wp
     call get_testmol('cytosine',mol)
     allocate (grad(3,mol%nat),qat(mol%nat),g(3,mol%nat))
-    call electrostatic_engrad(mol,0,'eeqbc',energy,grad,qat,io)
+    call electrostatic_core(mol,0,'eeqbc',energy,grad,qat,io)
     call check(error,io,0)
     if (allocated(error)) return
     rms = 0.0_wp; n = 0
     do iat = 1,mol%nat
       do k = 1,3
         mol%xyz(k,iat) = mol%xyz(k,iat)+h
-        call electrostatic_engrad(mol,0,'eeqbc',er,g,qat,io)
+        call electrostatic_core(mol,0,'eeqbc',er,g,qat,io)
         mol%xyz(k,iat) = mol%xyz(k,iat)-2.0_wp*h
-        call electrostatic_engrad(mol,0,'eeqbc',el,g,qat,io)
+        call electrostatic_core(mol,0,'eeqbc',el,g,qat,io)
         mol%xyz(k,iat) = mol%xyz(k,iat)+h
         dev = (er-el)/(2.0_wp*h)-grad(k,iat)
         rms = rms+dev*dev; n = n+1
@@ -199,7 +199,7 @@ subroutine test_solv_full_fd(error)
     call get_testmol('cytosine',mol)
     allocate (grad(3,mol%nat),geeq(3,mol%nat),qat(mol%nat),dqdr(3,mol%nat,mol%nat))
     !> analytic full gradient: charges + dq/dR from EEQ-BC, fed into ddX
-    call electrostatic_engrad(mol,0,'eeqbc',eeq,geeq,qat,io,dqdr=dqdr)
+    call electrostatic_core(mol,0,'eeqbc',eeq,geeq,qat,io,dqdr=dqdr)
     call ddx_pc_engrad(mol,qat,'cpcm',eps_water,energy,grad,io,dqdr=dqdr)
     call check(error,io,0)
     if (allocated(error)) return
@@ -229,7 +229,7 @@ subroutine solv_e(mol,energy,io)
     real(wp),allocatable   :: q(:),g(:,:)
     real(wp) :: eeq
     allocate (q(mol%nat),g(3,mol%nat))
-    call electrostatic_engrad(mol,0,'eeqbc',eeq,g,q,io)
+    call electrostatic_core(mol,0,'eeqbc',eeq,g,q,io)
     call ddx_e(mol,q,energy,io)
   end subroutine solv_e