HomePatch.C

Go to the documentation of this file.

/*
   HomePatch owns the actual atoms of a Patch of space
   Proxy(s) get messages via ProxyMgr from HomePatch(es)
   to update lists of atoms and their coordinates
   HomePatch(es) also have a Sequencer bound to them

   superclass:  Patch           
*/

#include <math.h>
#include "charm++.h"

#include "SimParameters.h"
#include "HomePatch.h"
#include "AtomMap.h"
#include "Node.h"
#include "PatchMap.inl"
#include "main.h"
#include "ProxyMgr.decl.h"
#include "ProxyMgr.h"
#include "Migration.h"
#include "Molecule.h"
#include "PatchMgr.h"
#include "Sequencer.h"
#include "LdbCoordinator.h"
#include "Settle.h"
#include "ReductionMgr.h"
#include "Sync.h"
#include "Random.h"
#include "Priorities.h"

#define TINY 1.0e-20;
#define MAXHGS 10
#define MIN_DEBUG_LEVEL 2
//#define DEBUGM
#include "Debug.h"

typedef int HGArrayInt[MAXHGS];
typedef BigReal HGArrayBigReal[MAXHGS];
typedef zVector HGArrayVector[MAXHGS];
typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS];
typedef zVector HGMatrixVector[MAXHGS][MAXHGS];

int average(CompAtom *qtilde,const HGArrayVector &q,BigReal *lambda,const int n,const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial);

void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGArrayVector &force,const int n, const int m, const HGArrayBigReal &imass,const HGArrayInt &ial,const HGArrayInt &ibl,const HGArrayVector &refab);


// DMK - Atom Separation (water vs. non-water)
#if NAMD_SeparateWaters != 0

// Macro to test if a hydrogen group represents a water molecule.
// NOTE: This test is the same test in Molecule.C for setting the
//   OxygenAtom flag in status.
// hgtype should be the number of atoms in a water hydrogen group
// It must now be set based on simulation parameters because we might
// be using tip4p

// DJH: This will give false positive for full Drude model,
//      e.g. O D H is not water but has hgs==3
#define IS_HYDROGEN_GROUP_WATER(hgs, mass)                 \
  ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))

#endif


HomePatch::HomePatch(PatchID pd, int atomCnt) : Patch(pd)
// DMK - Atom Separation (water vs. non-water)
#if NAMD_SeparateWaters != 0
  ,tempAtom()
#endif
{
  min.x = PatchMap::Object()->min_a(patchID);
  min.y = PatchMap::Object()->min_b(patchID);
  min.z = PatchMap::Object()->min_c(patchID);
  max.x = PatchMap::Object()->max_a(patchID);
  max.y = PatchMap::Object()->max_b(patchID);
  max.z = PatchMap::Object()->max_c(patchID);
  center = 0.5*(min+max);

  aAway = PatchMap::Object()->numaway_a();
  bAway = PatchMap::Object()->numaway_b();
  cAway = PatchMap::Object()->numaway_c();

  migrationSuspended = false;
  allMigrationIn = false;
  marginViolations = 0;
  patchMapRead = 0; // We delay read of PatchMap data
                    // to make sure it is really valid
  inMigration = false;
  numMlBuf = 0;
  flags.sequence = -1;

  numAtoms = atomCnt;
  replacementForces = 0;

  SimParameters *simParams = Node::Object()->simParameters;
  doPairlistCheck_newTolerance = 
        0.5 * ( simParams->pairlistDist - simParams->cutoff );


  numFixedAtoms = 0;
  //if ( simParams->fixedAtomsOn ) {
  //  for ( int i = 0; i < numAtoms; ++i ) {
  //    numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
  //  }
  //}

#ifdef NODEAWARE_PROXY_SPANNINGTREE
  ptnTree.resize(0);
  /*children = NULL;
  numChild = 0;*/
#else
  child =  new int[proxySpanDim];
  nChild = 0;   // number of proxy spanning tree children
#endif

#if CMK_PERSISTENT_COMM
  phsReady = 0;
  nphs = 0;
  localphs = new PersistentHandle[CkNumPes()];
  for (int i=0; i<CkNumPes(); i++) localphs[i] = 0;
#endif


  // DMK - Atom Separation (water vs. non-water)
  #if NAMD_SeparateWaters != 0

    // Create the scratch memory for separating atoms
    tempAtom.resize(numAtoms);
    numWaterAtoms = -1;

  #endif
  // Handle unusual water models here
  if (simParams->watmodel == WAT_TIP4) init_tip4();
  else if (simParams->watmodel == WAT_SWM4) init_swm4();

#ifdef MEM_OPT_VERSION
  isNewProxyAdded = 0;
#endif
}

HomePatch::HomePatch(PatchID pd, FullAtomList al) : Patch(pd), atom(al)
// DMK - Atom Separation (water vs. non-water)
#if NAMD_SeparateWaters != 0
  ,tempAtom()
#endif
{ 
  min.x = PatchMap::Object()->min_a(patchID);
  min.y = PatchMap::Object()->min_b(patchID);
  min.z = PatchMap::Object()->min_c(patchID);
  max.x = PatchMap::Object()->max_a(patchID);
  max.y = PatchMap::Object()->max_b(patchID);
  max.z = PatchMap::Object()->max_c(patchID);
  center = 0.5*(min+max);

  aAway = PatchMap::Object()->numaway_a();
  bAway = PatchMap::Object()->numaway_b();
  cAway = PatchMap::Object()->numaway_c();

  migrationSuspended = false;
  allMigrationIn = false;
  marginViolations = 0;
  patchMapRead = 0; // We delay read of PatchMap data
                    // to make sure it is really valid
  inMigration = false;
  numMlBuf = 0;
  flags.sequence = -1;

  numAtoms = atom.size();
  replacementForces = 0;

  SimParameters *simParams = Node::Object()->simParameters;
  doPairlistCheck_newTolerance = 
        0.5 * ( simParams->pairlistDist - simParams->cutoff );


  numFixedAtoms = 0;
  if ( simParams->fixedAtomsOn ) {
    for ( int i = 0; i < numAtoms; ++i ) {
      numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
    }
  }

#ifdef NODEAWARE_PROXY_SPANNINGTREE
  ptnTree.resize(0);
  /*children = NULL;
  numChild = 0;*/
#else
  child =  new int[proxySpanDim];
  nChild = 0;   // number of proxy spanning tree children
#endif

#if CMK_PERSISTENT_COMM
  phsReady = 0;
  nphs = 0;
  localphs = new PersistentHandle[CkNumPes()];
  for (int i=0; i<CkNumPes(); i++) localphs[i] = 0;
#endif


  // DMK - Atom Separation (water vs. non-water)
  #if NAMD_SeparateWaters != 0

    // Create the scratch memory for separating atoms
    tempAtom.resize(numAtoms);
    numWaterAtoms = -1;

    // Separate the current list of atoms
    separateAtoms();

  #endif
    
  // Handle unusual water models here
  if (simParams->watmodel == WAT_TIP4) init_tip4();
  else if (simParams->watmodel == WAT_SWM4) init_swm4();

#ifdef MEM_OPT_VERSION
  isNewProxyAdded = 0;
#endif

}

void HomePatch::write_tip4_props() {
  printf("Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
}

void HomePatch::init_tip4() {
  // initialize the distances needed for the tip4p water model

  Molecule *mol = Node::Object()->molecule;
  r_om = mol->r_om;
  r_ohc = mol->r_ohc;
}


void ::HomePatch::init_swm4() {
  // initialize the distances needed for the SWM4 water model
  SimParameters *simParams = Node::Object()->simParameters;
  Molecule *mol = Node::Object()->molecule;
  int ig;
  if (RIGID_NONE == simParams->rigidBonds) return;
  if (WAT_SWM4 != simParams->watmodel) return;
  for (ig = 0;  ig < numAtoms;  ig += atom[ig].hydrogenGroupSize ) {
    // find a water
    if (mol->rigid_bond_length(atom[ig].id) > 0) {
      // water is guaranteed by Molecule to have order:  O D LP H1 H2
      BigReal r_hh = mol->rigid_bond_length(atom[ig].id);
      BigReal r_oh = mol->rigid_bond_length(atom[ig+3].id);
      r_om = mol->rigid_bond_length(atom[ig+2].id);
      r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
      //printf("r_om and r_ohc initialized to %f and %f\n", r_om, r_ohc);
      break;
    }
  }
}

void HomePatch::reinitAtoms(FullAtomList al) {
  atom = al;
  numAtoms = atom.size();

  // DMK - Atom Separation (water vs. non-water)
  #if NAMD_SeparateWaters != 0

    // Reset the numWaterAtoms value
    numWaterAtoms = -1;

    // Separate the atoms
    separateAtoms();

  #endif
}

// Bind a Sequencer to this HomePatch
void HomePatch::useSequencer(Sequencer *sequencerPtr)
{ sequencer=sequencerPtr; }

// start simulation over this Patch of atoms
void HomePatch::runSequencer(void)
{ sequencer->run(); }

void HomePatch::readPatchMap() {
  // iout << "Patch " << patchID << " has " << proxy.size() << " proxies.\n" << endi;
  PatchMap *p = PatchMap::Object();
  PatchID nnPatchID[PatchMap::MaxOneAway];

  patchMigrationCounter = numNeighbors 
    = PatchMap::Object()->oneAwayNeighbors(patchID, nnPatchID);
  DebugM( 1, "NumNeighbors for pid " <<patchID<<" is "<< numNeighbors << "\n");
  int n;
  for (n=0; n<numNeighbors; n++) {
    realInfo[n].destNodeID = p->node(realInfo[n].destPatchID = nnPatchID[n]);
     DebugM( 1, " nnPatchID=" <<nnPatchID[n]<<" nnNodeID="<< realInfo[n].destNodeID<< "\n");
    realInfo[n].mList.resize(0);
  }

  // Make mapping from the 3x3x3 cube of pointers to real migration info
  for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
      for (int k=0; k<3; k++)
      {
        int pid =  p->pid(p->index_a(patchID)+i-1, 
            p->index_b(patchID)+j-1, p->index_c(patchID)+k-1);
        if (pid < 0) {
           DebugM(5, "ERROR, for patchID " << patchID <<" I got neigh pid = " << pid << "\n");
        }
        if (pid == patchID && ! (
                ( (i-1) && p->periodic_a() ) ||
                ( (j-1) && p->periodic_b() ) ||
                ( (k-1) && p->periodic_c() ) )) {
          mInfo[i][j][k] = NULL;
        }
        else {
          // Does not work as expected for periodic with only two patches.
          // Also need to check which image we want, but OK for now.  -JCP
          for (n = 0; n<numNeighbors; n++) {
            if (pid == realInfo[n].destPatchID) {
              mInfo[i][j][k] = &realInfo[n];
              break;
            }
          }
          if (n == numNeighbors) { // disaster! 
            DebugM(4,"BAD News, I could not find PID " << pid << "\n");
          }
        }
      }

  DebugM(1,"Patch("<<patchID<<") # of neighbors = " << numNeighbors << "\n");
}

HomePatch::~HomePatch()
{
#ifdef NODEAWARE_PROXY_SPANNINGTREE
    ptnTree.resize(0);
    delete [] children;
    #ifdef USE_NODEPATCHMGR
    delete [] nodeChildren;    
    #endif
#else
  delete [] child;
#endif
}


void HomePatch::boxClosed(int)
{
  if ( ! --boxesOpen )
  {
    if ( replacementForces ) {
      for ( int i = 0; i < numAtoms; ++i ) {
        if ( replacementForces[i].replace ) {
          for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
          f[Results::normal][i] = replacementForces[i].force;
        }
      }
      replacementForces = 0;
    }
    DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
        << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
    // only awaken suspended threads.  Then say it is suspended.
    sequencer->awaken();
    return;
  }
  else
  {
    DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
  }
}

void HomePatch::registerProxy(RegisterProxyMsg *msg) {
  DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
  proxy.add(ProxyListElem(msg->node,forceBox.checkOut()));

#ifdef MEM_OPT_VERSION
  isNewProxyAdded = 1;
#endif

  Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
  delete msg;
}

void HomePatch::unregisterProxy(UnregisterProxyMsg *msg) {
  int n = msg->node;
  ProxyListElem *pe = proxy.begin();
  for ( ; pe->node != n; ++pe );
  forceBox.checkIn(pe->forceBox);
  proxy.del(pe - proxy.begin());
  delete msg;
}

#if USE_TOPOMAP && USE_SPANNING_TREE

int HomePatch::findSubroots(int dim, int* subroots, int psize, int* pidscopy){
  int nChild = 0;
  int cones[6][proxySpanDim*proxySpanDim+proxySpanDim];
  int conesizes[6] = {0,0,0,0,0,0};
  int conecounters[6] = {0,0,0,0,0,0};
  int childcounter = 0;
  nChild = (psize>proxySpanDim)?proxySpanDim:psize;
  TopoManager tmgr;
  for(int i=0;i<psize;i++){
    int cone = tmgr.getConeNumberForRank(pidscopy[i]);
    cones[cone][conesizes[cone]++] = pidscopy[i];
  }

  while(childcounter<nChild){
    for(int i=0;i<6;i++){
      if(conecounters[i]<conesizes[i]){
        subroots[childcounter++] = cones[i][conecounters[i]++];
      }
    }
  }
  for(int i=nChild;i<proxySpanDim;i++)
    subroots[i] = -1;
  return nChild;
}
#endif // USE_TOPOMAP 

static int compDistance(const void *a, const void *b)
{
  int d1 = abs(*(int *)a - CkMyPe());
  int d2 = abs(*(int *)b - CkMyPe());
  if (d1 < d2) 
    return -1;
  else if (d1 == d2) 
    return 0;
  else 
    return 1;
}

void HomePatch::sendProxies()
{
  ProxyListIter pli(proxy);
  NodeIDList list;
  for ( pli = pli.begin(); pli != pli.end(); ++pli )
  {
    list.add(pli->node);
  }
  ProxyMgr::Object()->sendProxies(patchID, &list[0], list.size());  
}

#ifdef NODEAWARE_PROXY_SPANNINGTREE
void HomePatch::buildNodeAwareSpanningTree(void){
    //build the naive spanning tree for this home patch    
    int *proxyNodeMap = new int[CkNumNodes()]; //each element indiates the number of proxies residing on this node 
    NodeIDList proxyPeList;
    int proxyCnt = proxy.size();
    if(proxyCnt==0) {
        //this case will not happen in practice.
        //In debugging state where spanning tree is enforced, then this could happen
        //Chao Mei        
        #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
        DebugFileTrace *dft = DebugFileTrace::Object();
        dft->openTrace();
        dft->writeTrace("HomePatch[%d] has 0 proxy on proc[%d] node[%d]\n", patchID, CkMyPe(), CkMyNode());
        dft->closeTrace();
        #endif
        return;
    }
    proxyPeList.resize(proxyCnt);
    ProxyListIter pli(proxy);
    int i=0;
    for ( pli = pli.begin(); pli != pli.end(); ++pli, i++ )
        proxyPeList[i] = pli->node;    
    ProxyMgr::buildSinglePatchNodeAwareSpanningTree(patchID, proxyPeList, ptnTree, proxyNodeMap);    
    delete [] proxyNodeMap;
    proxyPeList.resize(0);

    //optimize on the naive spanning tree

    //setup the children
    setupChildrenFromProxySpanningTree();
    //send down to children
    sendNodeAwareSpanningTree();
}

void HomePatch::setupChildrenFromProxySpanningTree(){
    if(ptnTree.size()==0) {
        numChild = 0;
        delete [] children;
        children = NULL;
        #ifdef USE_NODEPATCHMGR
        numNodeChild = 0;
        delete [] nodeChildren;
        nodeChildren = NULL;        
        #endif
        return;
    }
    proxyTreeNode *rootnode = &ptnTree.item(0);
    CmiAssert(rootnode->peIDs[0] == CkMyPe());
    //set up children
    //1. add external children (the first proc inside the proxy tree node)    
    //2. add internal children (with threshold that would enable spanning    
    int internalChild;
    int externalChild;
    internalChild = rootnode->numPes-1;
    numChild = internalChild;
    if(numChild > inNodeProxySpanDim) {        
        //tree construction within a node)
        CmiAbort("Enabling in-node spanning tree construction has not been implemented yet!\n");
    }else{
        //exclude the root node
        int treesize = ptnTree.size();
        externalChild = (proxySpanDim>(treesize-1))?(treesize-1):proxySpanDim;
        numChild += externalChild;        

        delete [] children;
        #ifdef USE_NODEPATCHMGR
        delete [] nodeChildren;
        #endif
        if(numChild==0){
            children = NULL;
            #ifdef USE_NODEPATCHMGR
            nodeChildren = NULL;
            numNodeChild = 0;
            #endif
            return;
        }
        children = new int[numChild];    
        for(int i=0; i<externalChild; i++) {
            children[i] = ptnTree.item(i+1).peIDs[0];
        }
        for(int i=externalChild, j=1; i<numChild; i++, j++) {
            children[i] = rootnode->peIDs[j];
        }
    }

    #ifdef USE_NODEPATCHMGR
    //only register the cores that have proxy patches. The HomePach's core
    //doesn't need to be registered.
    CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
    NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
    if(rootnode->numPes==1){
        npm->registerPatch(patchID, 0, NULL);        
    }
    else{
        npm->registerPatch(patchID, rootnode->numPes-1, &rootnode->peIDs[1]);        
    }

    //set up childrens in terms of node ids
    numNodeChild = externalChild;
    if(internalChild) numNodeChild++;
    nodeChildren = new int[numNodeChild];    
    for(int i=0; i<externalChild; i++) {
        nodeChildren[i] = ptnTree.item(i+1).nodeID;        
    }
    //the last entry always stores this node id if there are proxies
    //on other cores of the same node for this patch.
    if(internalChild)
        nodeChildren[numNodeChild-1] = rootnode->nodeID;
    #endif
    
    #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
    DebugFileTrace *dft = DebugFileTrace::Object();
    dft->openTrace();
    dft->writeTrace("HomePatch[%d] has %d children: ", patchID, numChild);
    for(int i=0; i<numChild; i++)
        dft->writeTrace("%d ", children[i]);
    dft->writeTrace("\n");
    dft->closeTrace();
    #endif
    
}
#endif

#ifdef NODEAWARE_PROXY_SPANNINGTREE
//This is not an entry method, but takes an argument of message type
void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){
    //set up the whole tree ptnTree
    int treesize = msg->numNodesWithProxies;    
    ptnTree.resize(treesize);    
    int *pAllPes = msg->allPes;
    for(int i=0; i<treesize; i++) {
        proxyTreeNode *oneNode = &ptnTree.item(i);
        delete oneNode->peIDs;
        oneNode->numPes = msg->numPesOfNode[i];
        oneNode->nodeID = CkNodeOf(*pAllPes);
        oneNode->peIDs = new int[oneNode->numPes];
        for(int j=0; j<oneNode->numPes; j++) {
            oneNode->peIDs[j] = *pAllPes;
            pAllPes++;
        }
    }
    //setup children
    setupChildrenFromProxySpanningTree();
    //send down to children
    sendNodeAwareSpanningTree();
}

void HomePatch::sendNodeAwareSpanningTree(){
    if(ptnTree.size()==0) return;    
    ProxyNodeAwareSpanningTreeMsg *msg = 
        ProxyNodeAwareSpanningTreeMsg::getANewMsg(patchID, CkMyPe(), ptnTree.begin(), ptnTree.size());
   
    #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
    msg->printOut("HP::sendST");
    #endif

    CmiAssert(CkMyPe() == msg->allPes[0]);
    ProxyMgr::Object()->sendNodeAwareSpanningTree(msg);

}
#else
void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){}
void HomePatch::sendNodeAwareSpanningTree(){}
#endif

#ifndef NODEAWARE_PROXY_SPANNINGTREE
// recv a spanning tree from load balancer
void HomePatch::recvSpanningTree(int *t, int n)
{
  int i;
  nChild=0;
  tree.resize(n);
  for (i=0; i<n; i++) {
    tree[i] = t[i];
  }

  for (i=1; i<=proxySpanDim; i++) {
    if (tree.size() <= i) break;
    child[i-1] = tree[i];
    nChild++;
  }

  // send down to children
  sendSpanningTree();
}

void HomePatch::sendSpanningTree()
{
  if (tree.size() == 0) return;
  ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
  msg->patch = patchID;
  msg->node = CkMyPe();
  msg->tree = tree;
  ProxyMgr::Object()->sendSpanningTree(msg);  
}
#else
void HomePatch::recvSpanningTree(int *t, int n){}
void HomePatch::sendSpanningTree(){}
#endif

#ifndef NODEAWARE_PROXY_SPANNINGTREE
void HomePatch::buildSpanningTree(void)
{
  nChild = 0;
  int psize = proxy.size();
  if (psize == 0) return;
  NodeIDList oldtree = tree;
  int oldsize = oldtree.size();
  tree.resize(psize + 1);
  tree.setall(-1);
  tree[0] = CkMyPe();
  int s=1, e=psize+1;
  ProxyListIter pli(proxy);
  int patchNodesLast =
    ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
  int nNonPatch = 0;
#if 0
  int* pelists = new int[psize];
  for (int i=0; i<psize; i++) pelists[i] = -1;
  for ( pli = pli.begin(); pli != pli.end(); ++pli ) {
    int idx = rand()%psize;
    while (pelists[idx] != -1) { idx++; if (idx == psize) idx = 0; }
    pelists[idx] = pli->node;
  }
  for ( int i=0; i<psize; i++)
  {
    if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pelists[i]) ) {
      tree[--e] = pelists[i];
    } else {
      tree[s++] = pelists[i];
      nNonPatch++;
    }
  }
  delete [] pelists;
#else
    // try to put it to the same old tree
  for ( pli = pli.begin(); pli != pli.end(); ++pli )
  {
    int oldindex = oldtree.find(pli->node);
    if (oldindex != -1 && oldindex < psize) {
      tree[oldindex] = pli->node;
    }
  }
  s=1; e=psize;
  for ( pli = pli.begin(); pli != pli.end(); ++pli )
  {
    if (tree.find(pli->node) != -1) continue;    // already assigned
    if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
      while (tree[e] != -1) { e--; if (e==-1) e = psize; }
      tree[e] = pli->node;
    } else {
      while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
      tree[s] = pli->node;
      nNonPatch++;
    }
  }
#if 1
  if (oldsize==0 && nNonPatch) {
    // first time, sort by distance
    qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
  }
#endif

  //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);

#if USE_TOPOMAP && USE_SPANNING_TREE
  
  //Right now only works for spanning trees with two levels
  int *treecopy = new int [psize];
  int subroots[proxySpanDim];
  int subsizes[proxySpanDim];
  int subtrees[proxySpanDim][proxySpanDim];
  int idxes[proxySpanDim];
  int i = 0;

  for(i=0;i<proxySpanDim;i++){
    subsizes[i] = 0;
    idxes[i] = i;
  }
  
  for(i=0;i<psize;i++){
    treecopy[i] = tree[i+1];
  }
  
  TopoManager tmgr;
  tmgr.sortRanksByHops(treecopy,nNonPatch);
  tmgr.sortRanksByHops(treecopy+nNonPatch,
                                                psize-nNonPatch);  
  
  /* build tree and subtrees */
  nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
  delete [] treecopy;
  
  for(int i=1;i<psize+1;i++){
    int isSubroot=0;
    for(int j=0;j<nChild;j++)
      if(tree[i]==subroots[j]){
        isSubroot=1;
        break;
      }
    if(isSubroot) continue;
    
    int bAdded = 0;
    tmgr.sortIndexByHops(tree[i], subroots,
                                                  idxes, proxySpanDim);
    for(int j=0;j<proxySpanDim;j++){
      if(subsizes[idxes[j]]<proxySpanDim){
        subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
        bAdded = 1; 
        break;
      }
    }
    if( psize > proxySpanDim && ! bAdded ) {
      NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
    }
  }

#else /* USE_TOPOMAP && USE_SPANNING_TREE */
  
  for (int i=1; i<=proxySpanDim; i++) {
    if (tree.size() <= i) break;
    child[i-1] = tree[i];
    nChild++;
  }
#endif
#endif
  
#if 0
  // for debugging
  CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
  for (int i=0; i<psize+1; i++) {
    CkPrintf("%d ", tree[i]);
  }
  CkPrintf("\n");
#endif
  // send to children nodes
  sendSpanningTree();
}
#endif


void HomePatch::receiveResults(ProxyResultVarsizeMsg *msg){
    DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
    int n = msg->node;
    ProxyListElem *pe = proxy.begin();
    for ( ; pe->node != n; ++pe );
    Results *r = pe->forceBox->open();

    char *iszeroPtr = msg->isZero;
    Force *msgFPtr = msg->forceArr;

    for ( int k = 0; k < Results::maxNumForces; ++k )
    {
      Force *rfPtr = r->f[k];
      for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
          if((*iszeroPtr)!=1) {
              *rfPtr += *msgFPtr;
              msgFPtr++;
          }
      }      
    }
    pe->forceBox->close(&r);
    delete msg;
}

void HomePatch::receiveResults(ProxyResultMsg *msg)
{
  DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
  int n = msg->node;
  ProxyListElem *pe = proxy.begin();
  for ( ; pe->node != n; ++pe );
  Results *r = pe->forceBox->open();
  for ( int k = 0; k < Results::maxNumForces; ++k )
  {
    Force *f = r->f[k];
    register ForceList::iterator f_i, f_e;
    f_i = msg->forceList[k].begin();
    f_e = msg->forceList[k].end();
    for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
  }
  pe->forceBox->close(&r);
  delete msg;
}

void HomePatch::receiveResults(ProxyCombinedResultMsg *msg)
{
  DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodes.size()<<")\n");
//CkPrintf("[%d] Homepatch: %d receiveResults from %d nodes\n", CkMyPe(), patchID, n);
  for (int i=0; i<msg->nodes.size(); i++) {
    int node = msg->nodes[i];
    ProxyListElem *pe = proxy.begin();
    for ( ; pe->node != node; ++pe );
    Results *r = pe->forceBox->open();
    if (i==0) {
      for ( int k = 0; k < Results::maxNumForces; ++k )
      {
        Force *f = r->f[k];
        register ForceList::iterator f_i, f_e;
        f_i = msg->forceList[k].begin();
        f_e = msg->forceList[k].end();
        //for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;

        int nf = f_e - f_i;
#ifdef ARCH_POWERPC
#pragma disjoint (*f_i, *f)
#pragma unroll(4)
#endif
        for (int count = 0; count < nf; count++) {
          f[count].x += f_i[count].x;      
          f[count].y += f_i[count].y;      
          f[count].z += f_i[count].z;
        }
      }
    }
    pe->forceBox->close(&r);
  }
  delete msg;
}

void HomePatch::positionsReady(int doMigration)
{    
  flags.sequence++;

  if (!patchMapRead) { readPatchMap(); }
      
  if (numNeighbors) {
    if (doMigration) {
      doAtomMigration();
    } else {
      doMarginCheck();
    }
  }
  doMigration = (doMigration && numNeighbors) || ! patchMapRead;

  // Workaround for oversize groups
  doGroupSizeCheck();

  // Copy information needed by computes and proxys to Patch::p.
  p.resize(numAtoms);
  CompAtom *p_i = p.begin();
#ifdef MEM_OPT_VERSION
  pExt.resize(numAtoms);
  CompAtomExt *pExt_i = pExt.begin();
#endif
  FullAtom *a_i = atom.begin();
  int i; int n = numAtoms;
  for ( i=0; i<n; ++i ) { 
    p_i[i] = a_i[i]; 
    #ifdef MEM_OPT_VERSION
    pExt_i[i] = a_i[i];
    #endif
  }

  // Measure atom movement to test pairlist validity
  doPairlistCheck();

  if (flags.doMolly) mollyAverage();

  // Must Add Proxy Changes when migration completed!
  ProxyListIter pli(proxy);
  int *pids = NULL;
  int npid;
  if (proxySendSpanning == 0) {
    npid = proxy.size();
    pids = new int[npid];
    int *pidi = pids;
    int *pide = pids + proxy.size();
    int patchNodesLast =
      ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
    for ( pli = pli.begin(); pli != pli.end(); ++pli )
    {
      if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
        *(--pide) = pli->node;
      } else {
        *(pidi++) = pli->node;
      }
    }
  }
  else {
#ifdef NODEAWARE_PROXY_SPANNINGTREE
    #ifdef USE_NODEPATCHMGR
    npid = numNodeChild;
    if(numNodeChild>0) {
        pids = new int[npid];
        memcpy(pids, nodeChildren, sizeof(int)*npid);
    }
    #else
    npid = numChild;
    if(numChild>0) {
        pids = new int[numChild];
        memcpy(pids, children, sizeof(int)*numChild);
    }
    #endif
#else
    npid = nChild;
    pids = new int[proxySpanDim];
    for (int i=0; i<nChild; i++) pids[i] = child[i];
#endif
  }
  if (npid) {
#if CMK_PERSISTENT_COMM
    if (phsReady == 0)
      {
//CmiPrintf("Build on %d phs0:%d\n", CkMyPe(), localphs[0]);
     for (int i=0; i<npid; i++) {
       localphs[i] = CmiCreatePersistent(pids[i], 30000);
     }
     nphs = npid;
     phsReady = 1;
    }
#endif
    int seq = flags.sequence;
    int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
    //begin to prepare proxy msg and send it
    int pdMsgPLLen = p.size();
    int pdMsgAvgPLLen = 0;
    if(flags.doMolly) {
        pdMsgAvgPLLen = p_avg.size();
    }
    int pdMsgPLExtLen = 0;
#ifdef MEM_OPT_VERSION
    if(doMigration || isNewProxyAdded) {
        pdMsgPLExtLen = pExt.size();
    }
#endif
    ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgPLExtLen, PRIORITY_SIZE) ProxyDataMsg;
    SET_PRIORITY(nmsg,seq,priority);
    nmsg->patch = patchID;
    nmsg->flags = flags;
    nmsg->plLen = pdMsgPLLen;                
    //copying data to the newly created msg
    memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
    nmsg->avgPlLen = pdMsgAvgPLLen;        
    if(flags.doMolly) {
        memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
    }
    nmsg->plExtLen = pdMsgPLExtLen;
#ifdef MEM_OPT_VERSION
    if(doMigration || isNewProxyAdded){     
        memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
    }
#endif
    
#if NAMD_SeparateWaters != 0
    //DMK - Atom Separation (water vs. non-water)
    nmsg->numWaterAtoms = numWaterAtoms;
#endif

    
    #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
    DebugFileTrace *dft = DebugFileTrace::Object();
    dft->openTrace();
    dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
    for(int i=0; i<npid; i++) {
        dft->writeTrace("%d ", pids[i]);
    }
    dft->writeTrace("\n");
    dft->closeTrace();
    #endif

    if(doMigration) {
        ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
    }else{
        ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
    }
#if CMK_PERSISTENT_COMM
    CmiUsePersistentHandle(NULL, 0);
#endif
#ifdef MEM_OPT_VERSION
    isNewProxyAdded = 0;
#endif
  }
  delete [] pids;
  DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");

#ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
  positionPtrBegin = p.begin();
  positionPtrEnd = p.end();
#endif

  if(flags.doMolly) {
      avgPositionPtrBegin = p_avg.begin();
      avgPositionPtrEnd = p_avg.end();
  }
  Patch::positionsReady(doMigration);

  patchMapRead = 1;

  // gzheng
  if (useSync) Sync::Object()->PatchReady();
}

void HomePatch::replaceForces(ExtForce *f)
{
  replacementForces = f;
}

void HomePatch::saveForce(const int ftag)
{
  f_saved[ftag].resize(numAtoms);
  for ( int i = 0; i < numAtoms; ++i )
  {
    f_saved[ftag][i] = f[ftag][i];
  }
}

#undef DEBUG_REDISTRIB_FORCE 
#undef DEBUG_REDISTRIB_FORCE_VERBOSE
/* Redistribute forces from the massless lonepair charge particle onto
 * the other atoms of the water.
 *
 * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
 *
 * Pass by reference the forces (O H1 H2 LP) to be modified,
 * pass by constant reference the corresponding positions,
 * and a pointer to virial.
 */
void HomePatch::redistrib_lp_force(
    Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
    const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
    const Vector& p_lp, Tensor *virial) {

  Tensor wc;  // accumulate virial contribution from force redistribution

#ifdef DEBUG_REDISTRIB_FORCE 
  // Debug information to check against results at end

  // total force and torque relative to origin
  Vector totforce(0.0, 0.0, 0.0);
  Vector tottorque(0.0, 0.0, 0.0);

  totforce = f_ox + f_h1 + f_h2 + f_lp;

  tottorque += cross(f_ox, p_ox);
  tottorque += cross(f_h1, p_h1);
  tottorque += cross(f_h2, p_h2);
  //printf("Torque without LP is %f/%f/%f\n",
  //    tottorque.x, tottorque.y, tottorque.z);
  tottorque += cross(f_lp, p_lp);
  //printf("Torque with LP is %f/%f/%f\n",
  //    tottorque.x, tottorque.y, tottorque.z);
#endif

  // Calculate the radial component of the force and add it to the oxygen
  Vector r_ox_lp = p_lp - p_ox;
  BigReal rad_factor = (f_lp * r_ox_lp) * r_ox_lp.rlength() * r_ox_lp.rlength();
  Vector f_rad = r_ox_lp * rad_factor;

  Tensor vir = outer(f_rad, p_ox);
  wc += vir;

  f_ox = f_ox + f_rad;

  // Calculate the angular component
  Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
  Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1

  // deviation from collinearity of charge site
  Vector r_oop = cross(r_ox_lp, r_hcom_ox);
  //
  // vector out of o-h-h plane
  Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);

  // Here we assume that Ox/Lp/Hcom are linear
  // If you want to correct for deviations, this is the place

//  printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);

  Vector f_ang = f_lp - f_rad; // leave the angular component

  // now split this component onto the other atoms
  BigReal oxcomp = (r_hcom_ox.length() - r_ox_lp.length()) *
    r_hcom_ox.rlength();
  BigReal hydcomp = 0.5 * r_ox_lp.length() * r_hcom_ox.rlength();

  f_ox = f_ox + (f_ang * oxcomp);
  f_h1 = f_h1 + (f_ang * hydcomp);
  f_h2 = f_h2 + (f_ang * hydcomp);

  // Add virial contributions
  vir = outer(f_ang * oxcomp, p_ox);
  wc += vir;
  vir = outer(f_ang * hydcomp, p_h1);
  wc += vir;
  vir = outer(f_ang * hydcomp, p_h2);
  wc += vir;
  vir = outer(-1.0 * f_lp, p_lp);
  wc += vir;

  if ( virial ) *virial += wc;

  //Vector zerovec(0.0, 0.0, 0.0);
  f_lp = Vector(0.0, 0.0, 0.0);

#ifdef DEBUG_REDISTRIB_FORCE 
  // Check that the total force and torque come out right
  Vector newforce(0.0, 0.0, 0.0);
  Vector newtorque(0.0, 0.0, 0.0);
  BigReal error = 0.0;

  newforce = f_ox + f_h1 + f_h2;

  newtorque += cross(f_ox, p_ox);
  newtorque += cross(f_h1, p_h1);
  newtorque += cross(f_h2, p_h2);

  error = fabs(newforce.length() - totforce.length());
  if (error > 0.0001) {
     printf("Error:  Force redistribution for water "
         "exceeded force tolerance (%f vs. %f)\n",
         newforce.length(), totforce.length());
  }
#ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
  printf("Error in force length:  %f\n", error);
#endif

  error = fabs(newtorque.length() - tottorque.length());
  if (error > 0.0001) {
     printf("Error:  Force redistribution for water "
         "exceeded torque tolerance (%f vs. %f)\n",
         newtorque.length(), tottorque.length());
  }
#ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
  printf("Error in torque:  %f\n", error);
#endif
#endif /* DEBUG */
}

void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
    BigReal invdt) {
  // Reposition lonepair (Om) particle of Drude SWM4 water.
  // Same comments apply as to tip4_omrepos(), but the ordering of atoms
  // is different: O, D, LP, H1, H2.
  pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
  // Now, adjust velocity of particle to get it to appropriate place
  if (invdt != 0) {
    vel[2] = (pos[2] - ref[2]) * invdt;
  }
  // No virial correction needed since lonepair is massless
}

void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
  /* Reposition the om particle of a tip4p water
   * A little geometry shows that the appropriate position is given by
   * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O ) 
   * Here r_om is the distance from the oxygen to Om site, and r_ohc
   * is the altitude from the oxygen to the hydrogen center of mass
   * Those quantities are precalculated upon initialization of HomePatch
   *
   * Ordering of TIP4P atoms: O, H1, H2, LP.
   */

//  printf("rom/rohc are %f %f\n", r_om, r_ohc);
//  printf("Other positions are: \n  0: %f %f %f\n  1: %f %f %f\n  2: %f %f %f\n", pos[0].x, pos[0].y, pos[0].z, pos[1].x, pos[1].y, pos[1].z, pos[2].x, pos[2].y, pos[2].z);
  pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc); 
//  printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);

  // Now, adjust the velocity of the particle to get it to the appropriate place
  if (invdt != 0) {
    vel[3] = (pos[3] - ref[3]) * invdt;
  }

  // No virial correction needed, since this is a massless particle
  return;
}

void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
  // Loop over the patch's atoms and apply the appropriate corrections
  // to get all forces off of lone pairs
  ForceList *f_mod = f;
  for (int i = 0;  i < numAtoms;  i++) {
    if (atom[i].mass < 0.01) {
      // found lonepair
      redistrib_lp_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
          f_mod[ftag][i+2], f_mod[ftag][i],
          atom[i-2].position, atom[i+1].position,
          atom[i+2].position, atom[i].position, virial);
    }
  }
}

void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
  // Loop over the patch's atoms and apply the appropriate corrections
  // to get all forces off of lone pairs
  // Atom ordering:  O H1 H2 LP

  ForceList *f_mod =f;
  for (int i=0; i<numAtoms; i++) {
    if (atom[i].mass < 0.01) {
      // found lonepair
      redistrib_lp_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
          f_mod[ftag][i-1], f_mod[ftag][i],
          atom[i-3].position, atom[i-2].position,
          atom[i-1].position, atom[i].position, virial);
    }
  }
}


void HomePatch::addForceToMomentum(const BigReal timestep, const int ftag,
                                                        const int useSaved)
{
  SimParameters *simParams = Node::Object()->simParameters;
  const BigReal dt = timestep / TIMEFACTOR;
  ForceList *f_use = (useSaved ? f_saved : f);

  if ( simParams->fixedAtomsOn ) {
    for ( int i = 0; i < numAtoms; ++i ) {
      if ( atom[i].atomFixed || atom[i].mass <= 0. ) atom[i].velocity = 0;
      else atom[i].velocity += f_use[ftag][i] * ( dt * namd_reciprocal (atom[i].mass) );
    }
  } else {
    FullAtom *atom_arr  = atom.begin();
    const Force    *force_arr = f_use[ftag].const_begin();
#ifdef ARCH_POWERPC
#pragma disjoint (*force_arr, *atom_arr)
#endif
    for ( int i = 0; i < numAtoms; ++i ) {
      if (atom[i].mass == 0.) continue;
      BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.); 
      //printf("Taking reciprocal of mass %f\n", atom[i].mass);
      atom_arr[i].velocity.x += force_arr[i].x * recip_val;
      atom_arr[i].velocity.y += force_arr[i].y * recip_val;
      atom_arr[i].velocity.z += force_arr[i].z * recip_val;
    }
  }
}

void HomePatch::addVelocityToPosition(const BigReal timestep)
{
  SimParameters *simParams = Node::Object()->simParameters;
  const BigReal dt = timestep / TIMEFACTOR;
  if ( simParams->fixedAtomsOn ) {
    for ( int i = 0; i < numAtoms; ++i ) {
      if ( ! atom[i].atomFixed ) atom[i].position += atom[i].velocity * dt;
    }
  } else {
    FullAtom *atom_arr  = atom.begin();
    for ( int i = 0; i < numAtoms; ++i ) {
      atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
      atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
      atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
    }
  }
}

//  RATTLE algorithm from Allen & Tildesley
int HomePatch::rattle1(const BigReal timestep, Tensor *virial, 
    SubmitReduction *ppreduction)
{
  Molecule *mol = Node::Object()->molecule;
  SimParameters *simParams = Node::Object()->simParameters;
  const int fixedAtomsOn = simParams->fixedAtomsOn;
  const int useSettle = simParams->useSettle;
  const BigReal dt = timestep / TIMEFACTOR;
  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
  BigReal tol2 = 2.0 * simParams->rigidTol;
  int maxiter = simParams->rigidIter;
  int dieOnError = simParams->rigidDie;
  int i, iter;
  BigReal dsq[10], tmp;
  int ial[10], ibl[10];
  Vector ref[10];  // reference position
  Vector refab[10];  // reference vector
  Vector pos[10];  // new position
  Vector vel[10];  // new velocity
  Vector netdp[10];  // total momentum change from constraint
  BigReal rmass[10];  // 1 / mass
  int fixed[10];  // is atom fixed?
  Tensor wc;  // constraint virial
  BigReal idz, zmin;
  int nslabs;

  // Initialize the settle algorithm with water parameters
  // settle1() assumes all waters are identical,
  // and will generate bad results if they are not.
  // XXX this will move to Molecule::build_atom_status when that 
  // version is debugged
  if (!settle1isinitted()) {
    for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
      // find a water
      if (mol->rigid_bond_length(atom[ig].id) > 0) {
        int oatm;
        if (simParams->watmodel == WAT_SWM4) {
          oatm = ig+3;  // skip over Drude and Lonepair
          //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
          //    ig, atom[ig].mass, oatm, atom[oatm].mass);
        }
        else {
          oatm = ig+1;
          // Avoid using the Om site to set this by mistake
          if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
            oatm += 1;
          }
        }

        // initialize settle water parameters
        settle1init(atom[ig].mass, atom[oatm].mass, 
                    mol->rigid_bond_length(atom[ig].id), 
                    mol->rigid_bond_length(atom[oatm].id));
        break; // done with init
      }
    }
  }

  if (ppreduction) {
    nslabs = simParams->pressureProfileSlabs;
    idz = nslabs/lattice.c().z;
    zmin = lattice.origin().z - 0.5*lattice.c().z;
  }

  // Size of a hydrogen group for water
  int wathgsize = 3;
  int watmodel = simParams->watmodel;
  if (watmodel == WAT_TIP4) wathgsize = 4;
  else if (watmodel == WAT_SWM4) wathgsize = 5;
  
  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
    int hgs = atom[ig].hydrogenGroupSize;
    if ( hgs == 1 ) continue;  // only one atom in group
    // cache data in local arrays and integrate positions normally
    for ( i = 0; i < hgs; ++i ) {
      ref[i] = atom[ig+i].position;
      pos[i] = atom[ig+i].position;
      vel[i] = atom[ig+i].velocity;
      rmass[i] = (atom[ig+1].mass > 0. ? 1. / atom[ig+i].mass : 0.);
      //printf("rmass of %i is %f\n", ig+i, rmass[i]);
      fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
      // undo addVelocityToPosition to get proper reference coordinates
      if ( fixed[i] ) rmass[i] = 0.; else pos[i] += vel[i] * dt;
    }
    int icnt = 0;
    if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) {  // for water
      if (hgs != wathgsize) {
        NAMD_bug("Hydrogen group error caught in rattle1().");
      }
      // Use SETTLE for water unless some of the water atoms are fixed,
      // for speed we test groupFixed rather than the individual atoms
      if (useSettle && !atom[ig].groupFixed) {
        if (simParams->watmodel == WAT_SWM4) {
          // SWM4 ordering:  O D LP H1 H2
          // do swap(O,LP) and call settle with subarray O H1 H2
          // swap back after we return
          Vector lp_ref = ref[2];
          Vector lp_pos = pos[2];
          Vector lp_vel = vel[2];
          ref[2] = ref[0];
          pos[2] = pos[0];
          vel[2] = vel[0];
          settle1(ref+2, pos+2, vel+2, invdt);
          ref[0] = ref[2];
          pos[0] = pos[2];
          vel[0] = vel[2];
          ref[2] = lp_ref;
          pos[2] = lp_pos;
          vel[2] = lp_vel;
          // determine for LP updated pos and vel
          swm4_omrepos(ref, pos, vel, invdt);
        }
        else {
          settle1(ref, pos, vel, invdt);
          if (simParams->watmodel == WAT_TIP4) {
            tip4_omrepos(ref, pos, vel, invdt);
          }
        }

        // which slab the hydrogen group will belong to
        // for pprofile calculations.
        int ppoffset, partition;
        if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
          atom[ig+i].position = pos[i];
        } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
          atom[ig+i].velocity = vel[i];
        } else for ( i = 0; i < wathgsize; ++i ) {
          Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
          Tensor vir = outer(df, ref[i]);
          wc += vir;
          f[Results::normal][ig+i] += df;
          atom[ig+i].velocity = vel[i];
          if (ppreduction) {
            // put all the atoms from a water in the same slab.  Atom 0
            // should be the parent atom.
            if (!i) {
              BigReal z = pos[i].z;
              partition = atom[ig].partition;
              int slab = (int)floor((z-zmin)*idz);
              if (slab < 0) slab += nslabs;
              else if (slab >= nslabs) slab -= nslabs;
              ppoffset = 3*(slab + nslabs*partition);
            }
            ppreduction->item(ppoffset  ) += vir.xx;
            ppreduction->item(ppoffset+1) += vir.yy;
            ppreduction->item(ppoffset+2) += vir.zz;
          }
        }
        continue;
      }
      if ( !(fixed[1] && fixed[2]) ) {
        dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
      }
    }
    for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
      if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) > 0 ) {
        if ( !(fixed[0] && fixed[i]) ) {
          dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
        }
      }
    }
    if ( icnt == 0 ) continue;  // no constraints
    for ( i = 0; i < icnt; ++i ) {
      refab[i] = ref[ial[i]] - ref[ibl[i]];
    }
    for ( i = 0; i < hgs; ++i ) {
      netdp[i] = 0.;
    }
    int done;
    int consFailure;
    for ( iter = 0; iter < maxiter; ++iter ) {
//if (iter > 0) CkPrintf("iteration %d\n", iter);
      done = 1;
      consFailure = 0;
      for ( i = 0; i < icnt; ++i ) {
        int a = ial[i];  int b = ibl[i];
        Vector pab = pos[a] - pos[b];
        BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
        BigReal rabsq = dsq[i];
        BigReal diffsq = rabsq - pabsq;
        if ( fabs(diffsq) > (rabsq * tol2) ) {
          Vector &rab = refab[i];
          BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
          if ( rpab < ( rabsq * 1.0e-6 ) ) {
            done = 0;
            consFailure = 1;
            continue;
          }
          BigReal rma = rmass[a];
          BigReal rmb = rmass[b];
          BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
          Vector dp = rab * gab;
          pos[a] += rma * dp;
          pos[b] -= rmb * dp;
          if ( invdt != 0. ) {
            dp *= invdt;
            netdp[a] += dp;
            netdp[b] -= dp;
          }
          done = 0;
        }
      }
      if ( done ) break;
    }

    if ( consFailure ) {
      if ( dieOnError ) {
        iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
                        << (atom[ig].id + 1) << "!\n" << endi;
        return -1;  // triggers early exit
      } else {
        iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
                        << (atom[ig].id + 1) << "!\n" << endi;
      }
    } else if ( ! done ) {
      if ( dieOnError ) {
        iout << iERROR << "Exceeded RATTLE iteration limit for atom "
                        << (atom[ig].id + 1) << "!\n" << endi;
        return -1;  // triggers early exit
      } else {
        iout << iWARN << "Exceeded RATTLE iteration limit for atom "
                        << (atom[ig].id + 1) << "!\n" << endi;
      }
    }

    // store data back to patch
    int ppoffset, partition;
    if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
      atom[ig+i].position = pos[i];
    } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
      atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
    } else for ( i = 0; i < hgs; ++i ) {
      Force df = netdp[i] * invdt;
      Tensor vir = outer(df, ref[i]);
      wc += vir;
      f[Results::normal][ig+i] += df;
      atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
      if (ppreduction) {
        if (!i) {
          BigReal z = pos[i].z;
          int partition = atom[ig].partition;
          int slab = (int)floor((z-zmin)*idz);
          if (slab < 0) slab += nslabs;
          else if (slab >= nslabs) slab -= nslabs;
          ppoffset = 3*(slab + nslabs*partition);
        }
        ppreduction->item(ppoffset  ) += vir.xx;
        ppreduction->item(ppoffset+1) += vir.yy;
        ppreduction->item(ppoffset+2) += vir.zz;
      }
    }
  }
  if ( dt && virial ) *virial += wc;

  return 0;
}

//  RATTLE algorithm from Allen & Tildesley
void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
{
  Molecule *mol = Node::Object()->molecule;
  SimParameters *simParams = Node::Object()->simParameters;
  const int fixedAtomsOn = simParams->fixedAtomsOn;
  const int useSettle = simParams->useSettle;
  const BigReal dt = timestep / TIMEFACTOR;
  Tensor wc;  // constraint virial
  BigReal tol = simParams->rigidTol;
  int maxiter = simParams->rigidIter;
  int dieOnError = simParams->rigidDie;
  int i, iter;
  BigReal dsqi[10], tmp;
  int ial[10], ibl[10];
  Vector ref[10];  // reference position
  Vector refab[10];  // reference vector
  Vector vel[10];  // new velocity
  BigReal rmass[10];  // 1 / mass
  BigReal redmass[10];  // reduced mass
  int fixed[10];  // is atom fixed?
  
  // Size of a hydrogen group for water
  int wathgsize = 3;
  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;

  //  CkPrintf("In rattle2!\n");
  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
    //    CkPrintf("ig=%d\n",ig);
    int hgs = atom[ig].hydrogenGroupSize;
    if ( hgs == 1 ) continue;  // only one atom in group
    // cache data in lo