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

[chris.neale at utoronto.ca]

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 :-)
> 
> Cheers,
> 
> Erik
>