Added C API option to call gradient calculations on point charges

[pultar]

As of now, the C API of xtb does not grant access to functions regarding gradient calculations on point charges. We implemented this feature for our purposes and would like to share this little addition with you.

• the code compiles with GCC 8.2.0 and openblas as backend on a local workstation
• using the same compiler settings, the code can also be compiled on the cluster of ETH
• all tests pass with either set-up
• the new API function can successfully be called from a C++ program through via an include of xtb.h and linking to libxtb.so.6

Signed-off-by: Felix Pultar [email protected]

[awvwgk]

Thanks for sharing.

I think this makes a nice addition, however we need proper error checking in the API. Also, I have a hard time believing that the proposed implementation is actually working, since the function prototype is not matching the implementation (a unit test would be really nice).

Further, I'm not completely happy with the design, the point charge gradient should probably be queried from the results object rather than the calculator. Instead of passing the calculator it would be prudent to extend the single point procedure to save the point charge gradient in the results object and write a simple query to return the point charge gradient from there.

[pultar]

I understand your concern and agree we should add also a unit test, I'll set up something. Speaking of which, is there a reason that the API is only being tested with Intel's compiler?

xtb/test/api/meson.build

Line 17 in 9908849

if cc.get_id() == 'intel'

[awvwgk]

Speaking of which, is there a reason that the API is only being tested with Intel's compiler?

That's an artifact to avoid testing GCC / Apple-Clang mixtures on MacOS because Apple-Clang would run into problems with missing the OpenMP runtime, however I don't recall the reason why I checked for Intel rather than using the former condition.

[pultar]

Thank you, I will extend the test to GCC then once it is implemented correctly from the Fortran side. Another question about your second point. I agree that it is better to store gradients on point charges in the result object. I actually found the point charge gradients might otherwise be 0 when called via the calculator, even if I initialize them with random values like that:

double pcgrad[npc*7] =
      {1.0, 2.0, 3.0,
       1.0, 2.0, 3.0,
       1.0, 2.0, 3.0,
       1.0, 2.0, 3.0,
       1.0, 2.0, 3.0,
       1.0, 2.0, 3.0,
       1.0, 2.0, 3.0};

Maybe the calculator object actually does not store the calculation results?

Just to summarize the necessary steps to do it in line with other API functions (assuming I started fresh from branch main):

I would need to

(1) adapt the function prototype in xtb.h:

/// Query singlepoint results object for pc gradient in Hartree / Bohr
extern XTB_API_ENTRY void XTB_API_CALL
xtb_getPCGradient(xtb_TEnvironment /* env */,
                  xtb_TResults /* res */,
                  double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_1_0_0;

(2) add getPCGradient_api as a public member in src/api/results.f90

(3) implement the subroutine in the same file in an identical fashion to getGradient_api. Do some checks for correct allocations, and then get the point charge gradients like this: dptr(1:3, 1:size(res%pcem%grd, 2)) = res%pcem%grd

(4) in src/xtb/calculator.f90 :: subroutine singlepoint: make sure results of the point charge gradient calculator are stored in the field pcem of results

Would I also need to add pcem as a new field here, similar to gradient / gnorm?

xtb/src/type/data.f90

Line 28 in 7dd7710

type :: scc_results

This is the code area where I would have gotten the calculation results from. However, in my minimal example, the if statement is evaluated as false even though point charges have been provided.

xtb/src/xtb/calculator.f90

Line 190 in 7dd7710

!> check for external point charge field

I am not sure, however, if I need to change the signature of the singlepoint subroutine. During this process, a found a minor mistake in the C API example in the docs. I'll open another PR to address that one.

[pultar]

I added your suggested changes for now while I am working on the cleaner solution that goes via storing the results in a results object and querying them from there. In the meantime, I added a minimal working example as proof of concept:

https://github.com/pultar/xtb-example

I will also implement the unit test in a similar way.

[pultar]

I think I refactored everything the way you described it. The minimal example prints the expected gradients:

https://github.com/pultar/xtb-example/tree/point_charges_refactoring

I included gcc and clang for the C API test and also added my minimal example to it. All tests pass with cmake / make (Mac OS X with gcc-9 and CentOS 7 with gcc-8.2.0). However, if the tests are built with meson / ninja the assert in my test fails because the whole point charge array is filled with zeros. I made xtb more verbose (1e7d11a) and found that when I run the created test binary directly, everything is as expected.

Do you have any idea what could cause this behavior?

[awvwgk]

The chemical hardness for element 99 (Es) is not defined with neither GFN1-xTB nor GFN2-xTB and will therefore be zero leading to an infinity and therefore removing the interaction with the external charges effectively from the calculation.

Usually we use 7 (nitrogen) as dummy element for external charges.

[pultar]

The chemical hardness for element 99 (Es) is not defined with neither GFN1-xTB nor GFN2-xTB and will therefore be zero leading to an infinity and therefore removing the interaction with the external charges effectively from the calculation.

Usually we use 7 (nitrogen) as dummy element for external charges.

I was taking the value from this page of the docs: https://github.com/grimme-lab/xtb_docs/blob/main/source/pcem.rst

For future reference, is nitrogen your general recommendation if one does not know where the point charges come from?

[awvwgk]

They are no real point charges, rather Slater like densities parametrized with a hardness. Setting the hardness to a large value makes them point charge like.

There is a subtle issue there, we specify element symbols to get the hardness of the element, however if a number is given it is taken as hardness in Eh/e^2. In the API we don't bother dealing with strings and therefore use atomic numbers to identify elements, this leaves no room for overwriting the hardness for each element.

Recommendations are always difficult, we try to choose something that doesn't break most of the time, whatever this meant when the original implementation was done. For an actual calculation you should carefully check whether the hardnesses provided by the method work, you can use a stub hardness or better go for a point charge like behavior.

[awvwgk]

I noticed that the gradient norm is huge with a value of 14.0 Eh/a0, this is usually not wanted for tests because systems far away from the equilibrium will introduce some instability in the SCF solution. Since the repulsion energy with 26.7 Eh is also huge I would presume the cartesian coordinates are provided in Ångström rather than in Bohr, right?

[pultar]

I noticed that the gradient norm is huge with a value of 14.0 Eh/a0, this is usually not wanted for tests because systems far away from the equilibrium will introduce some instability in the SCF solution. Since the repulsion energy with 26.7 Eh is also huge I would presume the cartesian coordinates are provided in Ångström rather than in Bohr, right?

Nice catch, I copied the results from an orca calculation and forgot to change the units. The tests are running through now. I don't understand though how this can be dependent on the build system used.

Is the rest of my code correct and in line with the rest of the code base and design you had in mind? I hope I didn't oversee a check for allocation etc.

[awvwgk]

Nice catch, I copied the results from an orca calculation and forgot to change the units. The tests are running through now. I don't understand though how this can be dependent on the build system used.

I have no idea why it only shows up with meson, but that is the reason why we do testing on multiple platforms and build systems to also catch corner cases. One possible cause for the CMake build to pass would be the -DNDEBUG flag which removes the assert macro from the test code and renders it useless.

[pultar]

Nice catch, I copied the results from an orca calculation and forgot to change the units. The tests are running through now. I don't understand though how this can be dependent on the build system used.

I have no idea why it only shows up with meson, but that is the reason why we do testing on multiple platforms and build systems to also catch corner cases. One possible cause for the CMake build to pass would be the -DNDEBUG flag which removes the assert macro from the test code and renders it useless.

I see - I was worried I had introduced a memory leak somewhere. Good to see everything passes now.

[awvwgk]

Thanks for sharing. The implementation looks good to me.