Re: pair interaction energy for 2 water molecules

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Thu Dec 02 2010 - 21:35:23 CST

hugh,

On Thu, Dec 2, 2010 at 9:21 PM, Hugh Heldenbrand <helde010_at_umn.edu> wrote:
> Hello-
>
> My understanding of PME is pretty basic, but it seems like you could access
> the information you are looking for if you account for the self-energy of
> the charges within PME.  I am trying to do something similar at the moment,
> so I have a tcl script in VMD that goes something like this:

no. this cannot work. the kspace part of any ewald sum is not
pairwise additive, i.e. there is no way to divide this up in pairs.
what you compute is the interaction of each individual atom
in the (infinite!) lattice of all other atoms. removing the self-energy
contribution really doesn't help with that.

you can compute the explicit coulomb energy between each pair,
but then you have to stick with minimum image and that will not
add up to the total energy of the total system.

sorry,
   axel.

> 1. load the trajectory frame I want to calculate pairwise electrostatics for
>
> 2. write an individual coordinate file for each atom in the trajectory, to
> be used in the self-energy calculation
>
> 3. loop over all atom indices
>
>    4. loop over all atom indices less than the outer loop (to avoid double
> counting)
>
>    5. write a coordinate file for each pair of atoms
>
>    6. call a script that writes a parameter file for the individual atoms as
> well as the pair of atoms together
>
>    7. call a script that writes namd input files to do a single energy
> calculation for the individual atoms as well as the pair of atoms (using the
> same PME parameters as the calculation that originally created the
> trajectory).
>
>    8. call NAMD to do the 3 energy calculations (the two atoms together and
> each individually)
>
>    9. open each output file from the 3 NAMD calculations and parse them with
> a regular expression to get the electrostatic energy in each case
>
>    10.  The pairwise energy is the energy of the two atoms together minus
> each individual atom's self-energy
>
> I've only tried this on small systems of a few waters and a few ions, but if
> I total all the pairwise interactions I get within less than 1 kcal of the
> reported total electrostatic energy of the original trajectory frame.  I
> don't know how useful my actual scripts would be to you though, Divya, since
> I am using the AMBER forcefield and their parameter file generation tools.
>
> -Hugh Heldenbrand
> PhD candidate
> University of Minnesota
>
>
>
> On 11/29/2010 04:56 AM, Axel Kohlmeyer wrote:
>>
>> On Mon, Nov 29, 2010 at 5:00 AM, divya nayar<divya.alchemist_at_gmail.com>
>>  wrote:
>>>
>>> Hi,
>>> I want to calculate the pair interaction energy between 2 water molecules
>>> in
>>> a system of 256 waters. My conf file is below.
>>> I am able to reproduce the van der Waal's component of pair interaction
>>> energy by manual calculations using LJ potential defined in User's guide.
>>> I
>>> want to ask how should I change my PME parameters to tell NAMD to
>>> calculate
>>> the electrostatic energy between only those two tagged water molecules in
>>> a
>>> system of 256 water molecules because the electrostatic energy  for 2
>>> molecules that I am getting is not very much different from that of
>>> between
>>> one molecule interacting with all other 255 molecules.
>>
>> you cannot do what you describe _at all_ with PME enabled
>> or any other type of ewald sum.  you have a component of the
>> electrostatic interaction that is not pair-wise additive and thus
>> cannot be attributed to arbitrary pairs of atoms or molecules.
>>
>> axel.
>>
>>
>>> Also, is there a way I can check the values tht I am getting for
>>> electrostatic energy are correct?
>>> structure          ./../tip4p.psf
>>> coordinates        ../../tip4p-prod260.coor
>>>
>>> set outputname     4tip.tpe
>>>
>>> temperature         260
>>> seed                4010
>>>
>>> paraTypeCharmm      on
>>> #parameters          ../par_all27_prot_lipid.inp
>>> parameters          ../../tip4p.par
>>>
>>> # Force-Field Parameters
>>> exclude             scaled1-4
>>> 1-4scaling          1.0
>>> dielectric              1.0
>>> cutoff              9.5
>>> switching           on
>>> switchdist          7.5
>>> pairlistdist        11.5
>>> margin              0.0
>>> # Integrator Parameters
>>> timestep            1.0  ;# 2fs/step
>>> rigidBonds          all  ;# needed for 2fs steps
>>> nonbondedFreq       1
>>> fullElectFrequency  2
>>> stepspercycle       20
>>> rigidTolerance      0.00001
>>> rigidIterations     100
>>> wrapAll             on
>>> wrapWater           on
>>>
>>> waterModel        tip4
>>> # Periodic Boundary Conditions
>>> PME                 yes
>>> PMETolerance        0.000001
>>> PMEGridSizeX        40
>>> PMEGridSizeY        40
>>> PMEGridSizeZ        40
>>>
>>>
>>> # Output
>>> outputName          $outputname
>>> outputTiming        500
>>> extendedSystem      ../../tip4p-prod260.xsc
>>> pairInteraction     on
>>> pairInteractionFile  ../pdb/abc.pdb
>>>
>>> pairInteractionCol B
>>> pairInteractionGroup1 1
>>> pairInteractionGroup2 2
>>> #############################################################
>>> ## EXECUTION SCRIPT                                        ##
>>> #############################################################
>>> set ts 100
>>>
>>> coorfile open dcd   ../8000.dcd
>>>
>>> while {! [coorfile read] } {
>>>    firstTimestep $ts
>>>    run 0
>>>    incr ts 100
>>> }
>>>
>>> coorfile close
>>>
>>>
>>>
>>>
>>> Please guide me. Thanks in advance
>>>
>>> Divya
>>
>>
>
>

-- 
Dr. Axel Kohlmeyer
akohlmey_at_gmail.com  http://goo.gl/1wk0
Institute for Computational Molecular Science
Temple University, Philadelphia PA, USA.

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:54:50 CST