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

[Berk Hess]

On 05/30/2012 10:36 AM, Gerrit Groenhof wrote:
>  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
gather_energy_bspline get you the energy at a position of a charge.
If you need it at an arbitrary location in space, that can be 
implemented easily.

But I don't see why you would need the periodic image interactions 
between QM atom alone.
That doesn't make much sense. Periodic interactions are strongly 
screened by all the other atoms.
As all electrostatic interactions are strongly coupled, one can not 
single out some contribution and
apply corrections to that alone.

Cheers,

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