Re: pair interaction energy for 2 water molecules

From: Hugh Heldenbrand (helde010_at_umn.edu)
Date: Thu Dec 02 2010 - 20:21:01 CST

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:

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

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:56:26 CST