[gmx-developers] Extracting i-j forces from nb_generic.c

Thu Jan 30 09:54:30 CET 2014

Dear all,

Following your advice I have been adding the heat-flux calculation to
gromacs, running on 1 processor. My test with a simple LJ atom is about 50%
off other literature values, and while I don't rule out a trivial mistake
on my part, I wonder if I am pulling the information about atomic
interactions in the correct manner or not. Would you mind confirming if at
least my intended methodology is correct?

1. In the MD loop I first pull velocity and kinetic energy from update.c
(lines 1069 and 1078 respectively)

https://github.com/jamesresearching/gromacs-heatflux/commit/fb1624e11085422d24c58a2f6118ff590bf12913#diff-6

2. then force from nb_generic.c (~line 457). For LJ atoms, coulomb
interactions don't matter of course, but for charged atoms later, for the
GMX_NBKERNEL_ELEC_EWALD case (line 311) I manually added a short-range
calculation for charge interaction. Is there any reason why this would not
be reasonable?

https://github.com/jamesresearching/gromacs-heatflux/commit/fb1624e11085422d24c58a2f6118ff590bf12913#diff-2

3. Finally, I calculate heat-flux and write out the value. I use leap-frog
integration - is my understanding correct that the energies and forces
collected in steps 1 and 2 are all on the same step?

https://github.com/jamesresearching/gromacs-heatflux/commit/fb1624e11085422d24c58a2f6118ff590bf12913#diff-0

(I realise this last writing-out part could probably be better-placed,
eventually).

In addition to the above questions, if there is any glaring
mis-understanding on my part about this implementation within Gromacs, any
comment is of course greatly appreciated.

Thank you.
With best regards,

James

On 10 January 2014 04:31, Berk Hess <hess at kth.se> wrote:

>  On 01/09/2014 07:30 PM, Mark Abraham wrote:
>
>
> On Jan 9, 2014 7:44 AM, "James" <jamesresearching at gmail.com> wrote:
> >
> > Dear Erik, David and other Gromacs developers,
> >
> > Thank you very much for your replies. They were very helpful! I have run
> with just one thread and I can now see how to extract i-j VDW and
> short-range coulomb potentials (shared equally between pairs) and by
> tracing the temperature evaluation I think I have found the appropriate
> place for extracting kinetic energy and velocities.
> >
> > I wonder if I may follow up about the atomic numbering. I will need to
> move to parallel calculation as soon as possible, and I have been thinking
> about the comments:
> >
> > (Erik) "you can write functionality to move the properties you need
> (just as we move charge, etc) into the local topology when we do
> rebalancing between domains"
> >
> > I'm not sure what exactly is meant by "local topology". Where is
> "rebalancing between domains" performed? Would you mind expanding this
> explanation a little more or giving an example?
>
> See code in runner.c and md.c, probably near init_domain_decomposition or
> domain_decomposition. You'll need to drill into either or both to find how
> to propagate whatever it is you need.
>
> > (David) "There is a function in mtop_util.h that will give you the
> global atom number"
> >
> > Are you referring to the gmx_mtop_atomlookup_init function? I looked at
> other places in the code where this is called; I assume I would have to
> somehow pass *mtop to the force subroutine? If evaluating this is
> expensive, as Erik suggests, I suppose I could do this outside of the force
> loop.
>
> No, that is a global topology, and local atom indices won't help. IIRC
> gmx_local_topology_t is the data type. You should want to build an array at
> domain-decomposition time that can be looked up at each following
> integration step with a local atom index. Or simulate with pen and paper
> for better speed ;-)
>
> The array is already there: cr->dd->gatindex
>
> Cheers,
>
> Berk
>
>  Mark
>
> > Thank you again for your help.
> >
> > With best regards,
> >
> > James
> >
> >
> > On 31 December 2013 18:15, Erik Lindahl <erik.lindahl at scilifelab.se>
> wrote:
> >>
> >> Hi James,
> >>
> >> On 31 Dec 2013, at 08:15, James <jamesresearching at gmail.com> wrote:
> >>
> >> > Dear Gromacs developers,
> >> >
> >> > To measure thermal conductivity using Green-Kubo I am trying to
> extract the i-j atomic forces from nb_generic.c and I wonder if anyone
> could help me with some basic questions?
> >> >
> >> > 1. If I understand correctly, ii and jj are atom numbers for atom i
> and atom j. These don't seem to correspond to the atoms in the order that
> they were read in (PDB file), however. How can I recover this order, or
> what is the new order?
> >>
> >> That would be "ii" and "jnr" in the master branch. If you are only
> using a single thread this should correspond to the atom order in the input
> (although C always starts counting at 0). However, when we run in parallel
> the domain decomposition will mean that the particle indices change
> (frequently) as particules move across node boundaries.
> >>
> >> There are routines you can use to find the global atom indices, but you
> don't want to call those from the innermost kernel. Instead I would
> recommend using a single thread while you develop, and later when you want
> to enable parallelization you can write functionality to move the
> properties you need (just as we move charge, etc) into the local topology
> when we do rebalancing between domains.
> >>
> >> > 2. What units are the positions ix, iy, iz (again in nb_generic.c)?
> >>
> >> nm - same as in the gro files, but it differs from PDB files that use
> Angstrom.
> >> >
> >> > 3. nb_generic is within the loop i = i0 to i1-1 in nonbonded.c. What
> is this loop? Somehow related to different interaction types?
> >>
> >> At one point we experimented with multithreading parallelism over this
> loop (thus the indices), but now the outermost loop is over different
> neighborlists. A certain particle will have different neighborlists for
> neighbors that only have charge, only Lennard-Jones, and both LJ+charge.
> >>
> >> Having said that, for the future we are likely moving entirely to the
> new style "verlet" kernels currently used for GPUs (and that will likely be
> default for CPUs in Gromacs-5.0 too), but we still need to write a generic
> kernel there.
> >>
> >> > 4. Later, I will also want to extract the instantaneous kinetic and
> potential energy of the atom in the same time-step - any advice is very
> welcome.
> >>
> >> For kinetic energy you can just look at how we calculate temperature -
> this is easy. Potential energy is complex if you use lattice summation
> (PME), since you then would need to modify the grid interpolation to get
> potentials of individual particles.
> >> >
> >> > 5. Out of interest, I notice many of the files in the release-4-6
> branch are in include/ but the same files are in legacyheaders/ in the
> master branch. Why is this? Which should I be developing against?
> >>
> >> We are moving pretty fast right now due to a switch from C to partial
> C++, which also involves a lot of cleaning-up and modularization. The
> master branch is the future, but legacyheaders is stuff that still hasn't
> been modularized.
> >>
> >> Cheers,
> >>
> >> Erik
> >> --
> >> Gromacs Developers mailing list
> >>
> >> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
> posting!
> >>
> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>
> >> * For (un)subscribe requests visit
> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developersor send a mail to
> gmx-developers-request at gromacs.org.
> >
> >
> >
> > --
> > Gromacs Developers mailing list
> >
> > * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
> posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developersor send a mail to
> gmx-developers-request at gromacs.org.
>
>
>
>
> --
> Gromacs Developers mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developersor send a mail to
> gmx-developers-request at gromacs.org.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20140130/5f0562a2/attachment-0001.html>