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

Gerrit Groenhof ggroenh at gwdg.de
Wed May 30 10:36:26 CEST 2012

  Hi Berk,

> Hi,
>
> From what you write I seems to be much simpler.
>
> If you have charges for the qm atoms you can simply do PME as usual,
> you would only need to subtract the interactions between the qm atoms,
> which would be done in the same way as normal exclusions.
> Or am I missing something?

Maybe, but I am nor sure myself.

If I simply use gmx_pme_do, I won't get the potential at the atoms this
way, only the forces I think.

We need to add to the QM hamiltonian terms like
\frac{1}{2}S_{\mu\nu}^{AB}*(V(R_A)+V(R_B)), where \mu and \nu are
basisfunctions located on atom A and B. V(R) is the electrostatic
potential due to the MM atoms, as well as the QM atoms's own periodic images
S is an overlap integral, and R_A/B the position of the atoms. A and B
can be the same.

Gerrit

>
> Cheers,
>
> Berk
>
> On 05/30/2012 10:11 AM, Gerrit Groenhof wrote:
>>  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
>>
>