[gmx-developers] Extracting i-j forces from nb_generic.c
James
jamesresearching at gmail.com
Sun Feb 2 03:35:41 CET 2014
Further to my previous email, I tried it for water and got values about 10x
literature values. After a detailed search I'm unable to identify a trivial
error, so I'm wondering if my understanding of how gromacs works is somehow
wrong.
I realise everyone is busy with the upcoming version 5.0 and that heat flux
is not a priority for biological simulation software, but if an experienced
gromacs developer is able to spot something that I am unable to do, this
would be immensely helpful.
For want of better ideas, I guessed maybe the short-range charge
contribution to heat flux (qkqf and gkqv) which I calculate independently
when coming through GMX_NBKERNEL_ELEC_EWALD was somehow not being cut off,
so I changed
if (bExactElecCutoff)
{
felec = (rsq <= rcoulomb2) ? felec : 0.0;
velec = (rsq <= rcoulomb2) ? velec : 0.0;
gkqf = (rsq <= rcoulomb2) ? gkqf : 0.0;
gkqv = (rsq <= rcoulomb2) ? gkqv : 0.0;
}
to
if (bExactElecCutoff)
{
felec = (rsq <= rcoulomb2) ? felec : 0.0;
velec = (rsq <= rcoulomb2) ? velec : 0.0;
}
gkqf = (rsq <= rcoulomb2) ? gkqf : 0.0;
gkqv = (rsq <= rcoulomb2) ? gkqv : 0.0;
https://github.com/jamesresearching/gromacs-heatflux/commit/fb1624e11085422d24c58a2f6118ff590bf12913#diff-2
(lines 311 and 338)
but this actually increases the potential (gkqv) contribution rather than
diminishing it! So it appears I am not understanding something fundamental
here.
Thank you in advance for your time.
With best regards,
James
On 30 January 2014 17:54, James <jamesresearching at gmail.com> wrote:
> 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/20140202/13490b70/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list