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

[Gerrit Groenhof]

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

I agree, but I am not singling out a contribution here. I'll try to explain what I want to do. There may well be a flaw in our thinking:

For this specific QM methods self-consistent charge DFTB one varies the charges on the atoms variationally. Thus in the SCF cycle, one starts with a set of guess charges on the atoms, and constructs the fock matrix with all one and two electron integrals. Most of these are parameters. 

However, the interaction with the environment, in a periodic system, should include the interaction of the QM atoms with their own periodic images.

(This is similar to bonded interactions. The atoms see their bonding partners inside a box via the bonded potential, but interaction with the periodic images of themselves and their bonding partners via the coulomb interactions.)

Then, if we have the potential at the sites of the QM atoms, we diagonilze the fock matrix to get the molecular orbitals. However, with these new orbitals we get changes in the QM charges (and electronic integrals), so that we need to construct the fock matrix again. The electrostatic term should be updated as well, as the charges on the QM atoms in the periodic images have changed.

After reaching the variational minimum, we still would need to compute the QM/MM forces, and these I thought we get from a normal pme call with the final charges.

best,

Gerrit




> 
> 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
>>>> 
>>> 
>> 
> 
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.