[gmx-users] selected forces

[Mark Abraham]

----- Original Message -----
From: Ali Hassanali <ahassana at chemistry.ohio-state.edu>
Date: Friday, September 10, 2010 2:21
Subject: [gmx-users] selected forces
To: gmx-users at gromacs.org

> Dear Gromacs Users,
> 
> I've read through various posts on a related issue but am hoping 
> that someone will have some ideas on this or at least point me 
> in the right direction. I am also hoping that some discussion on 
> this will help develop some routines that would allow users to 
> approach this problem as a black box:
> 
> The system I intend to simulate consists of a polymer surrounded 
> by water, lets call them as energygroups Protein and SOL. At 
> each time step during a molecular dynamics simulation, I need 
> (or would like to) have access to the force on each solvent atom 
> exerted only by the protein atoms and not the other solvent 
> atoms. For now lets assume that the protein is uncharged and so 
> only interacts with the solvent via LJ interactions.
> I've been ploughing through the code and have gotten to the following:
> 
> 1) do_force(f) has a force vector array f that is invoked in md.c
> 2) do_force(f) in sim_util.c calls do_force_lowlevel(f) to 
> determine the bonded and non-bonded interactions.
> 3) do_force_lowlevel(f) in force.f calls do_nonbonded(f) which I 
> assume will determine the LJ interactions for all particles.
> 4) do_nonbonded(f) in nonbonded.c seems to have different kernel 
> routines for water and other atoms (as the manual states, forces 
> on water molecules are determined with special kernel loops)
> 5) do_nonbonded(f) calls a specific nonbonded kernel function 
> and this is where its getting a bit tricky. Any help would be 
> appreciated. I'd like to access the solvent atoms and sum up the 
> forces coming from just the protein/polymer and not the other 
> solvent atoms.

You've certainly done your homework - that's a fairly good picture of how things work. If you define a suitable energy monitor groups (see 3.3 of manual) then GROMACS will calculate energies and forces within and between those groups, e.g. SOL-SOL, SOL-Protein, Protein-Protein. So it occurs to me that what you want to do is arrange the ordering of your groups so that GROMACS computes SOL-Protein first (that might require some coding to reorder the looping - step through with a debugger to see how this is working). Then the contents of the f vector will be the LJ + charge forces that you want, so you just need to interrupt the looping over energy groups at that point and do whatever it is you need to do. (Check that bonded forces haven't been computed yet.)

If you really need LJ only, and charges are present in both groups, then you will need to rework the nblist construction so that this pair of groups is sent first to an LJ-only kernel, then to a charge-only kernel, so that you can interrupt when you need to. Whatever you do, don't touch the code inside any kernel or your simulations will still be running when your grandchildren retire :-)

Later when you consider charge interactions, be aware that this kind of group-wise decomposition is impossible for the long-range component of PME, if you use that.

Mark

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100910/4a6afab7/attachment.html>