[gmx-developers] PBC with p2sub1 from pbc_dx() modifications
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Wed May 31 01:54:40 CEST 2006
Fantastic! Thank you so much for the details. It will take some weeks to test
that everything is working as it should, but I will be sure to post as soon as
everything is tested. In the end, I believe that a modern computer should still
get 0.5-1ns/day on a simple protein-membrane system, which is about 5 times
faster than charmm for me, making it a very worthwhile endevour. Once working, I
will be pulling a palmityl-(oleyl-phosphatidylethanolamine) out of a gram
negative outer menbrane protein (PagP) binding pocket by a derivative of replica
exchange and measuring the PMF making it quite important the the membrane is
happy. I really want to emphasize my appreciation for the rapid and in-depth
comments of Berk and Erik. And I won't hesitate to point out that the code is
generally far more readable than the humble Erik lets on.
Quoting Erik Lindahl <lindahl at sbc.su.se>:
> Hi Chris,
> I've actually given this some thought in more general terms, and the
> appropriate long-term solution is probably to completely separate the
> asymmetric units/symmetry from the periodicity of the unit cell.
> Let's start with the neighborsearching. This code might look
> complicated primarily due to a very deep nested loop over shifts and
> dimensions, but it's not really _that_ hard to understand. Referring
> to src/mdlib/ns.c:
> search_neighbors() is a wrapper function that does some setup,
> assigns grid cell indices to atoms (with fill_grid() in nsgrid.c) and
> then calls ns5_core() to do the real search work. All
> neighborsearching is based on charge group centers-of-geometry rather
> than atoms for reasons of efficiency.
> Still, _exclusions_ from NS need to be defined on an atom-by-atom
> basis, so we solve this in a pretty smart way with bit masks. Charge
> groups are currently allowed to contain up to 32 atoms, so before
> doing advanced grid stuff and locating neighbors for a particular
> group we go through the (usually very short) list of excluded atoms
> for this ("i particule") group and for each of these ("j particle")
> set the bits of the corresponding excluded atoms to 1. You probably
> won't need to alter this, it's just to give you a hint what it does.
> in ns5_core() we go through the neighbor grid cells that might
> contain neighbor particles - the number of cells we need to check
> depends on the box geometry, so we first determine exactly how many
> neighbor cells we should check, and in which directions. However,
> instead of checking periodic boundary conditions for every single
> distance we use "shift vectors" here too. Think of an imaginary 1D-
> box where the current cell is the leftmost, and we want to find
> neighbors both to the right (next box), and to the left (using PBC,
> in the rightmost box). This corresponds to a list of two shift
> vectors. First we search the next box without doing anything special.
> For the second step, we first add a shift vector corresponding to the
> 1D-box length, which you can think of as moving the entire reference
> box to the right of the system, and then search the rightmost box
> without any need for checking PBC.
> You WILL need to understand this "nonbonded shift" concept - there
> are some illustrations in the manual that might help, as well as the
> Gromacs papers!
> Finally, we check distances between charge groups and eventually go
> through atoms and put them in the lists. This is just a bunch of
> special cases to optimize the neighborlists as much as possible, e.g.
> VdW-only, or coulomb-interactions-only between pairs of water molecules.
> As explained in the manual, we need to store the shift indices/
> vectors for each neighborlist to be able to calculate the virial as a
> single sum after calling the nonbonded forces instead of doing it
> inside the kernel. I somewhat suspect that it might be easiest for
> you just to copy coordinates to get a full periodic cell so you don't
> have to worry about rotations in the shift - that could namely affect
> the innerloops too, and you don't want to hack assembly ;-)
> Now for something completely different: shifts.
> There is an related, but slightly different, periodicity problem with
> bonded interactions. We really should have used a different name, but
> unfortunately it's called "molecular shift" or even only "shift" in
> the code. So sorry. Mea culpa...
> The graph code goes through all molecules at setup time and makes
> them continuous, i.e. removes PBC. We then re-apply PBC, but keep
> track of how each atom has been moved/shifted in terms of integer
> multiples of the unit cell vectors.
> Neighborsearching and nonbonded force calculation normally puts all
> charge groups inside the central unit cell, but the molecular shifts
> enable us to rapidly make all molecules continuous again without
> evaluating a single conditional. I'm not sure if you need bother with
> this, but I think we explained this too in more detail in the latest
> JCC Gromacs paper.
> PME is fairly independent from everything else, in src/mdlib/pme.c.
> You provide coordinates and triclinic box vectors and get potential,
> forces on each atom, and long-range virial back.
> My recommendation would be to start with neighborsearch for a trivial
> system without bonded interactions or exclusions, and make sure you
> get the neighborlist correct :-)
More information about the gromacs.org_gmx-developers