[gmx-developers] Implementation of fluctuating charge into GROMACS: Pressure calculation

[Lee-Ping Wang]

Hi there,

 

I've gotten the QTPIE interaction to work with PBC using a custom-designed
shift function; the matrix elements are taken smoothly to zero at the cutoff
radius.  This ensures that there is no force discontinuity, which I've found
causes some artifacts in the dynamics and energy minimization. 

 

Implementing the Ewald summation has been a bit problematic, because as you
said, there needs to be a pairwise decomposition of the electrostatic energy
in computing the matrix elements.  I think I have something that works, but
the forces aren't implemented yet, and it's expensive because it requires a
double loop in the k-space summation (the normal Ewald summation only needs
a single loop).  Hopefully, the shift function will perform well enough and
I won't worry so much about the Ewald summation.

 

I am trying to build the QTPIE forces into the "fshift" vector to get the
pressures right, but I don't have any way of telling whether my calculated
pressures are correct.  I suspect there is a problem, since the box keeps on
blowing up if I turn pressure coupling on, but I need some reference method
to compare against; how do you all make sure that the pressure is being
computed properly?  Do you simply repeat energy calculations at different
box sizes and compute the pressure by finite difference?

 

Thanks,

 

-          Lee-Ping

 

From: gmx-developers-bounces at gromacs.org
[mailto:gmx-developers-bounces at gromacs.org] On Behalf Of Gerrit Groenhof
Sent: Wednesday, October 13, 2010 2:34 AM
To: Discussion list for GROMACS development
Subject: Re: [gmx-developers] rlist larger than L/2 in simple neighbor
search

 

On 10/12/2010 06:49 PM, Lee-Ping Wang wrote: 

Hi there,

 

For the past few months, I've been adding a charge equilibration model to
GROMACS (to be specific, QTPIE - a refinement of QEq).  For the model to
work properly, each atom must "see" each other atom in the electrostatic
interaction (it is a screened Coulomb interaction).  This is easy enough to
do with pbc = no; I just set all the cutoffs to zero, and add a new
neighbor-list that includes excluded atom pairs.


Since we'd be interested in that, I would appreciate if you can keep us
informed of your progress.



 

Now I would like to run the simulation with PBC turned on, which means
either implementing Ewald summation or using cutoffs for the electrostatic
interaction.  When cutoffs are turned on (with or without PBC), the model
goes haywire even with very long cutoff lengths.  I haven't implemented
Ewald summation for the screened Coulomb interaction, and I'm not sure if
it's worth the effort if an easier fix is possible.


If I recall correctly, the screened Coulomb here means the overlap integral
between s-type Slater functions. This interaction behave as a normal coulomb
at large distances, so that coudl explain the issues with increasing the
cut-off. But an Ewald expression will be difficult as well, since  the all
pairs-wise interactions are required to compute the charge flow parameters.
I have so far only seen it applied to molecules. For non bonded system like
a box of water, a direct charge flow between atoms far apart is unlikely
though. Could one not simply scale also the charge transfer variables to
zero with distance as well?

 

As a quick fix, I changed the simple neighbor-searching code to set rcut =
1000000 when QTPIE is turned on.

 

I am worried about the other possible consequences of my change; will this
adversely affect the molecular dynamics or create artifacts in some
unexpected way?  

You will get the ordinary artefacts caused by cutoff, like ordering etc.

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20101016/34b657ae/attachment.html>