Force on external point charge or external electric field by electron density?

[chen850512337]

Is your feature request related to a problem? Please describe.
I am trying to write a QM/MM interface of xtb and charmm MD package, but I cannot find a proper way to calculated the force on MM point charge.

Describe the solution you'd like
If the force on MM point charge is given directly is better.
Another way is to give the electric field on each point charge site, so the force can be calculated easily.
Hope this can be in the C API.

Describe alternatives you've considered
Using the cube file to calculate the force on each point can also work, I just want to know if there is a better way to go.

[awvwgk]

The API already supports setting point charges via

xtb/include/xtb.h

Lines 187 to 199 in 0245411

	/// Add a external charge potential to calculator (only supported in GFN1/2-xTB)
	extern XTB_API_ENTRY void XTB_API_CALL
	xtb_setExternalCharges(xtb_TEnvironment /* env */,
	xtb_TCalculator /* calc */,
	int* /* n */,
	int* /* numbers [n] */,
	double* /* charges [n] */,
	double* /* positions [n][3] */) XTB_API_SUFFIX__VERSION_1_0_0;
	/// Unset the external charge potential
	extern XTB_API_ENTRY void XTB_API_CALL
	xtb_releaseExternalCharges(xtb_TEnvironment /* env */,
	xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_1_0_0;

The gradient will also be calculated but cannot be retrieved from the API since there is no call for this yet.

---

Some notes on the current state of the external point charges, in case anyone want to give it a try. The external point charges are wrapped in a derived type here inside the calculator object

xtb/src/xtb/calculator.f90

Lines 73 to 74 in 0245411

	!> External potential
	type(tb_pcem) :: pcem

which has a field for the gradients which will be populated after the singlepoint calculation

xtb/src/type/pcem.f90

Line 30 in 0245411

real(wp),allocatable :: grd(:,:)

The only access via the API to the external point charge field is its deconstructor here

xtb/src/api/calculator.f90

Lines 502 to 532 in 0245411

	!> Unset the external charge potential
	subroutine releaseExternalCharges_api(venv, vcalc) &
	& bind(C, name="xtb_releaseExternalCharges")
	character(len=*), parameter :: source = 'xtb_api_releaseExternalCharges'
	type(c_ptr), value :: venv
	type(VEnvironment), pointer :: env
	type(c_ptr), value :: vcalc
	type(VCalculator), pointer :: calc
	if (c_associated(venv)) then
	call c_f_pointer(venv, env)
	call checkGlobalEnv
	if (.not.c_associated(vcalc)) then
	call env%ptr%error("Singlepoint calculator is not allocated", source)
	return
	end if
	call c_f_pointer(vcalc, calc)
	! There are only limited cases when we need to perform any action.
	! Deconstructors should behave nicely, therefore no errors on wrong input.
	if (allocated(calc%ptr)) then
	select type(xtb => calc%ptr)
	type is(TxTBCalculator)
	call xtb%pcem%deallocate
	end select
	end if
	end if
	end subroutine releaseExternalCharges_api

But the logic to retrieve the gradient in a getter function would be quite similar.

[Lennard94]

Hi,
a quick-hack for the those running into the same issue:
In xtb/src/api/calculator.f90, you can add the function

subroutine getPCGradient_api(vcalc, dptr) &
      & bind(C, name="xtb_getPCGradient")
   character(len=*), parameter :: source = "xtb_api_getPCGradient"
   type(c_ptr), value :: vcalc
   type(VCalculator), pointer :: calc
   real(c_double), intent(inout) :: dptr(3, *)
   call c_f_pointer(vcalc, calc)
   select type(xtb => calc%ptr)

   type is(TxTBCalculator)
   dptr(1:3, 1:size(xtb%pcem%grd, 2)) = xtb%pcem%grd

   end select
end subroutine getPCGradient_api

and in the xtb.h file

xtb_getPCGradient(xtb_TCalculator,
                double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_1_0_0;

Best
Lennard

Edit: Obvious mistake in the header file, corrected

[awvwgk]

Close, but the C-API export doesn't match the provided header. Be careful to check for the existence of the elements you are accessing, otherwise the application will segfault.

Here is something that might work (I didn't compile nor test it with xtb, but if it works, feel free to submit a patch).

 !> Retrieve point charge gradient
 subroutine getPCGradient_api(venv, vcalc, gradient) & 
       & bind(C, name="xtb_getPCGradient") 
    character(len=*), parameter :: source = 'xtb_getPCGradient'
    type(c_ptr), value :: venv 
    type(VEnvironment), pointer :: env 
    type(c_ptr), value :: vcalc 
    type(VCalculator), pointer :: calc 
    real(c_double), intent(inout) :: gradient(*)
  
    if (c_associated(venv)) then 
       call c_f_pointer(venv, env) 
       call checkGlobalEnv 
  
       if (.not.c_associated(vcalc)) then 
          call env%ptr%error("Singlepoint calculator is not allocated", source) 
          return 
       end if 
       call c_f_pointer(vcalc, calc) 
  
       if (.not.allocated(calc%ptr)) then
          call env%ptr%error("Singlepoint calculator is not initialized", source) 
          return
       end if

       select type(xtb => calc%ptr) 
       type is(TxTBCalculator)
          if (.not.allocated(xtb%pcem%grd)) then
             call env%ptr%error("Point charge gradient not available due to missing point charge data", source)
             return
          end if
          gradient(1:size(xtb%pcem%grd)) = reshape(xtb%pcem%grd, [size(xtb%pcem%grd)])
       class default
          call env%ptr%error("Calculator does not provide point charge gradient", source) 
       end select 
 
    end if  
 end subroutine releaseExternalCharges_api 

Which provides the C-API

extern XTB_API_ENTRY void XTB_API_CALL 
xtb_getPCGradient(xtb_TEnvironment /* env */,
                  xtb_TCalculator /* calc */,
                  double* /* gradient [n][3] */);

[Lennard94]

The code snippets are actually tested (without segfaults), however, I did not copy pasted the complete changes of the files to this thread. Anyway, I guess the pointers you gave in your first response (in addition with your last comment) are sufficient to make the required changes.

Edit: You are obviously right, I pasted the wrong version
For the header file

xtb_getPCGradient(xtb_TCalculator,
                double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_1_0_0;

Best