[gmx-developers] how to retrieve the electrostatic potential from pme?

[Gerrit Groenhof]

  Dear All,

To include also the effect of coulomb interactions between MM and QM 
atoms, I need the electrostatic potential V(r). For the moment I need it 
only at the positions of the QM nuclei.

But, I do not understand what is going on in pme.c.

What I learned from browsing pme.c is that after computing the potential 
on the grid in k space, it is transformed back on the real space grid 
and then forces are interpolated  by gather_f_bsplines.

However, I do not understand if there are already forces on the grid 
points that are interpolated, or that the potential is interpolated at 
the site of the atom, after which forces on that atom are computed, and 
the interpolated potential is not further used. The comment after the 
routine suggests the latter is done, but that's a comment of 13 years 
ago.  Furhtermore, in the spline interpolation, the output looks like a 
force already: it is different in x,y and z, and is called f?

Or should I be using the gather_energy_bslpine instead? Is it possible 
then to vary the position within the gridcell? Or shoud I stick to the 
position of atom n?

Then, if we manage to get v(r), I thoughT of testing it the following way:

1) We compute the total (MM) electrostatic energy, with the QM atoms 
exluded among themselved. In that case, I think that the QM atoms see 
their periodic counterparts (charges are assiggend to the atoms at every 
step of the SCF), but not the QM atoms in the central box.

2) we add a probe charge of +1 at one of the QM nuclei (and exclude its 
interaction from the QM atoms) and recompute the total MM electrostatic 
energy.

3) we remove all atoms (MM + QM) from the box, and only keep the probe 
charge, and comute the total electrostatic energy of the probe charge in 
the periodic system

4) The potential at the position of the probe charge is then: E(2) 
-E(1)-E(3)?

Are there any objections for this?


Best wishes,

Gerrit