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

[Gerrit Groenhof]

  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
>>
>