[gmx-developers] PBC with p2sub1 from pbc_dx() modifications

[Erik Lindahl]

Hi Chris,

On May 30, 2006, at 9:57 PM, chris.neale at utoronto.ca wrote:
>
> Thanks for the pointers. It would have taken me forever to figure  
> that out. It
> looked easy enough in pbc_dx(), so I will give it a try in the  
> proper places. My
> coding is good, but it may take some time to actually find the (i) 
> grid neighbour
> search, (ii)pressure, (iii) graph code, (iv) PME, (v) anything else  
> you can
> think of, routines as this is my first adventure into the gromacs  
> code. If
> anybody knows the names of the functions off-hand it would be a  
> great help. I
> will post the modifications as soon as I have it running = well  
> before any
> publication. I realize that I am fishing for some rather detailed  
> direction. If
> this is not an appropriate type of question then I sincerly  
> appologize.

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

Cheers,

Erik