28 #include "ProxyMgr.decl.h"
46 #include "ComputeQMMgr.decl.h"
53 #define MIN_DEBUG_LEVEL 2
70 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);
74 #define MASS_EPSILON (1.0e-35) //a very small floating point number
78 #if NAMD_SeparateWaters != 0
89 #define IS_HYDROGEN_GROUP_WATER(hgs, mass) \
90 ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))
94 #ifdef TIMER_COLLECTION
95 const char *TimerSet::tlabel[TimerSet::NUMTIMERS] = {
116 settle_initialized = 0;
119 rattleListValid =
false;
124 numGBISP1Arrived = 0;
125 numGBISP2Arrived = 0;
126 numGBISP3Arrived = 0;
127 phase1BoxClosedCalled =
false;
128 phase2BoxClosedCalled =
false;
129 phase3BoxClosedCalled =
false;
137 center = 0.5*(min+max);
142 aAwayDist = (max.x - min.x) * aAway;
149 bAwayDist = (max.y - min.y) * bAway;
156 cAwayDist = (max.z - min.z) * cAway;
161 migrationSuspended =
false;
162 allMigrationIn =
false;
163 marginViolations = 0;
169 flags.maxForceUsed = -1;
171 numAtoms = atom.size();
172 replacementForces = 0;
175 doPairlistCheck_newTolerance =
181 for (
int i = 0; i < numAtoms; ++i ) {
182 numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
186 #ifdef NODEAWARE_PROXY_SPANNINGTREE
195 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
203 #if NAMD_SeparateWaters != 0
206 tempAtom.resize(numAtoms);
221 void HomePatch::write_tip4_props() {
222 printf(
"Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
225 void HomePatch::init_tip4() {
233 void ::HomePatch::init_swm4() {
253 #if NAMD_SeparateWaters != 0
266 { sequencer=sequencerPtr; }
270 { sequencer->
run(); }
272 void HomePatch::readPatchMap() {
277 patchMigrationCounter = numNeighbors
279 DebugM( 1,
"NumNeighbors for pid " <<
patchID<<
" is "<< numNeighbors <<
"\n");
281 for (n=0; n<numNeighbors; n++) {
282 realInfo[n].
destNodeID = p->
node(realInfo[n].destPatchID = nnPatchID[n]);
283 DebugM( 1,
" nnPatchID=" <<nnPatchID[n]<<
" nnNodeID="<< realInfo[n].destNodeID<<
"\n");
288 for (
int i=0; i<3; i++)
289 for (
int j=0; j<3; j++)
290 for (
int k=0; k<3; k++)
295 DebugM(5,
"ERROR, for patchID " <<
patchID <<
" I got neigh pid = " << pid <<
"\n");
301 mInfo[i][j][k] = NULL;
306 for (n = 0; n<numNeighbors; n++) {
307 if (pid == realInfo[n].destPatchID) {
308 mInfo[i][j][k] = &realInfo[n];
312 if (n == numNeighbors) {
313 DebugM(4,
"BAD News, I could not find PID " << pid <<
"\n");
318 DebugM(1,
"Patch("<<
patchID<<
") # of neighbors = " << numNeighbors <<
"\n");
324 #ifdef NODEAWARE_PROXY_SPANNINGTREE
326 #ifdef USE_NODEPATCHMGR
327 delete [] nodeChildren;
337 phase1BoxClosedCalled =
true;
345 }
else if (box == 6) {
347 }
else if (box == 7) {
349 }
else if (box == 8) {
350 phase2BoxClosedCalled =
true;
360 }
else if (box == 9) {
362 }
else if (box == 10) {
370 if ( replacementForces ) {
371 for (
int i = 0; i <
numAtoms; ++i ) {
372 if ( replacementForces[i].replace ) {
377 replacementForces = 0;
379 DebugM(1,
patchID <<
": " << CthSelf() <<
" awakening sequencer "
380 << sequencer->thread <<
"(" <<
patchID <<
") @" << CmiTimer() <<
"\n");
383 phase3BoxClosedCalled =
true;
388 if (numGBISP1Arrived == proxy.
size() &&
389 numGBISP2Arrived == proxy.
size() &&
390 numGBISP3Arrived == proxy.
size()) {
408 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
417 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
422 for ( ; *pe != n; ++pe );
428 #if USE_TOPOMAP && USE_SPANNING_TREE
430 int HomePatch::findSubroots(
int dim,
int* subroots,
int psize,
int* pidscopy){
433 int conesizes[6] = {0,0,0,0,0,0};
434 int conecounters[6] = {0,0,0,0,0,0};
435 int childcounter = 0;
438 for(
int i=0;i<psize;i++){
439 int cone = tmgr.getConeNumberForRank(pidscopy[i]);
440 cones[cone][conesizes[cone]++] = pidscopy[i];
443 while(childcounter<nChild){
444 for(
int i=0;i<6;i++){
445 if(conecounters[i]<conesizes[i]){
446 subroots[childcounter++] = cones[i][conecounters[i]++];
454 #endif // USE_TOPOMAP
458 int d1 = abs(*(
int *)a - CkMyPe());
459 int d2 = abs(*(
int *)b - CkMyPe());
471 CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
479 #ifdef NODEAWARE_PROXY_SPANNINGTREE
480 void HomePatch::buildNodeAwareSpanningTree(
void){
481 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
482 DebugFileTrace *dft = DebugFileTrace::Object();
484 dft->writeTrace(
"HomePatch[%d] has %d proxy on proc[%d] node[%d]\n",
patchID, proxy.
size(), CkMyPe(), CkMyNode());
485 dft->writeTrace(
"Proxies are: ");
486 for(
int i=0; i<proxy.
size(); i++) dft->writeTrace(
"%d(%d), ", proxy[i], CkNodeOf(proxy[i]));
487 dft->writeTrace(
"\n");
498 ProxyMgr::buildSinglePatchNodeAwareSpanningTree(
patchID, proxy, ptnTree);
502 setupChildrenFromProxySpanningTree();
507 void HomePatch::setupChildrenFromProxySpanningTree(){
508 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
511 if(ptnTree.size()==0) {
515 #ifdef USE_NODEPATCHMGR
517 delete [] nodeChildren;
523 CmiAssert(rootnode->
peIDs[0] == CkMyPe());
527 int internalChild = rootnode->
numPes-1;
528 int externalChild = ptnTree.size()-1;
532 if(internalSlots==0) {
536 internalChild = (internalSlots>internalChild)?internalChild:internalSlots;
540 nChild = externalChild+internalChild;
545 child =
new int[nChild];
547 for(
int i=0; i<externalChild; i++) {
548 child[i] = ptnTree.item(i+1).peIDs[0];
550 for(
int i=externalChild, j=1; i<nChild; i++, j++) {
551 child[i] = rootnode->
peIDs[j];
554 #ifdef USE_NODEPATCHMGR
557 CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
567 numNodeChild = externalChild;
568 if(internalChild) numNodeChild++;
569 delete [] nodeChildren;
570 nodeChildren =
new int[numNodeChild];
571 for(
int i=0; i<externalChild; i++) {
572 nodeChildren[i] = ptnTree.item(i+1).nodeID;
577 nodeChildren[numNodeChild-1] = rootnode->
nodeID;
580 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
581 DebugFileTrace *dft = DebugFileTrace::Object();
583 dft->writeTrace(
"HomePatch[%d] has %d children: ",
patchID, nChild);
584 for(
int i=0; i<nChild; i++)
585 dft->writeTrace(
"%d ", child[i]);
586 dft->writeTrace(
"\n");
592 #ifdef NODEAWARE_PROXY_SPANNINGTREE
597 ptnTree.resize(treesize);
598 int *pAllPes = msg->
allPes;
599 for(
int i=0; i<treesize; i++) {
601 delete [] oneNode->
peIDs;
603 oneNode->
nodeID = CkNodeOf(*pAllPes);
605 for(
int j=0; j<oneNode->
numPes; j++) {
606 oneNode->
peIDs[j] = *pAllPes;
611 setupChildrenFromProxySpanningTree();
617 if(ptnTree.size()==0)
return;
621 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
625 CmiAssert(CkMyPe() == msg->
allPes[0]);
634 #ifndef NODEAWARE_PROXY_SPANNINGTREE
641 for (i=0; i<n; i++) {
646 if (tree.
size() <= i)
break;
647 child[i-1] = tree[i];
651 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
661 if (tree.
size() == 0)
return;
664 msg->
node = CkMyPe();
673 #ifndef NODEAWARE_PROXY_SPANNINGTREE
677 int psize = proxy.
size();
678 if (psize == 0)
return;
680 int oldsize = oldtree.
size();
691 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
693 int oldindex = oldtree.
find(*pli);
694 if (oldindex != -1 && oldindex < psize) {
695 tree[oldindex] = *pli;
699 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
701 if (tree.
find(*pli) != -1)
continue;
703 while (tree[e] != -1) { e--;
if (e==-1) e = psize; }
706 while (tree[s] != -1) { s++;
if (s==psize+1) s = 1; }
712 if (oldsize==0 && nNonPatch) {
720 #if USE_TOPOMAP && USE_SPANNING_TREE
723 int *treecopy =
new int [psize];
735 for(i=0;i<psize;i++){
736 treecopy[i] = tree[i+1];
740 tmgr.sortRanksByHops(treecopy,nNonPatch);
741 tmgr.sortRanksByHops(treecopy+nNonPatch,
745 nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
748 for(
int i=1;i<psize+1;i++){
750 for(
int j=0;j<nChild;j++)
751 if(tree[i]==subroots[j]){
755 if(isSubroot)
continue;
758 tmgr.sortIndexByHops(tree[i], subroots,
759 idxes, proxySpanDim);
761 if(subsizes[idxes[j]]<proxySpanDim){
762 subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
767 if( psize > proxySpanDim && ! bAdded ) {
768 NAMD_bug(
"HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
775 if (tree.
size() <= i)
break;
776 child[i-1] = tree[i];
784 CkPrintf(
"[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(),
patchID, psize, nNonPatch);
785 for (
int i=0; i<psize+1; i++) {
786 CkPrintf(
"%d ", tree[i]);
803 char *iszeroPtr = msg->
isZero;
809 for(
int i=0; i<msg->
flLen[k]; i++, rfPtr++, iszeroPtr++) {
810 if((*iszeroPtr)!=1) {
831 for ( ; f_i != f_e; ++f_i, ++
f ) *f += *f_i;
847 int nf = msg->
flLen[k];
849 #pragma disjoint (*f_i, *f)
851 for (
int count = 0; count < nf; count++) {
853 f[count].
x += f_i->
x;
854 f[count].
y += f_i->
y;
855 f[count].
z += f_i->
z;
878 ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
879 CkpvAccess(BOCclass_group).computeQMMgr) ;
887 if ( subP != NULL ) {
893 if (qmAtomGroup[subP->
newID] > 0 && simParams->
PMEOn) {
900 for (
int qmIter=0; qmIter<numQMAtms; qmIter++) {
902 if (qmAtmIndx[qmIter] == subP->
newID) {
931 if (!patchMapRead) { readPatchMap(); }
942 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
948 if (doMigration && simParams->
qmLSSOn)
951 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
958 #if defined(NAMD_CUDA) || defined(NAMD_AVXTILES) || \
959 (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
960 int *ao =
new int[n];
965 for (
int j=0; j<n; ++j ) {
967 if ( a_i[j].atomFixed ) ao[--k2] = j;
973 for (
int j=0; j<n; ++j ) {
980 for (
int i=0; i<n; ++i ) {
985 for (
int i = 0; i < n; ++i) {
994 const BigReal ucenter_x = ucenter.
x;
995 const BigReal ucenter_y = ucenter.
y;
996 const BigReal ucenter_z = ucenter.
z;
998 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
1000 n_16 = (n + 15) & (~15);
1001 cudaAtomList.
resize(n_16);
1003 mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
1004 mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
1005 mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
1006 mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
1007 #elif defined(NAMD_AVXTILES)
1008 int n_avx = (n + 15) & (~15);
1009 cudaAtomList.
resize(n_avx);
1011 tiles.realloc(n, ac);
1017 for (
int k=0; k<n; ++k ) {
1018 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
1020 atom_x[k] = a[j].
position.
x - ucenter_x;
1021 atom_y[k] = a[j].
position.
y - ucenter_y;
1022 atom_z[k] = a[j].
position.
z - ucenter_z;
1023 atom_q[k] = charge_scaling * a[j].
charge;
1029 ac[k].
q = charge_scaling * a[j].
charge;
1032 #ifdef NAMD_AVXTILES
1036 for (
int k=n; k<n_avx; ++k ) {
1045 #ifdef NAMD_AVXTILES
1047 }
else doMigration = doMigration && numNeighbors;
1050 doMigration = doMigration && numNeighbors;
1052 doMigration = doMigration || ! patchMapRead;
1054 doMigration = doMigration || doAtomUpdate;
1055 doAtomUpdate =
false;
1066 #if ! defined(NAMD_CUDA)
1067 #if defined(NAMD_AVXTILES)
1083 for ( i=0; i<n; ++i ) {
1098 numGBISP1Arrived = 0;
1099 phase1BoxClosedCalled =
false;
1100 numGBISP2Arrived = 0;
1101 phase2BoxClosedCalled =
false;
1102 numGBISP3Arrived = 0;
1103 phase3BoxClosedCalled =
false;
1104 if (doMigration || isNewProxyAdded)
1109 if (doMigration || isNewProxyAdded) {
1117 int pidsPreAllocated = 1;
1120 npid = proxy.
size();
1121 pids =
new int[npid];
1122 pidsPreAllocated = 0;
1124 int *pide = pids + proxy.
size();
1125 int patchNodesLast =
1127 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
1137 #ifdef NODEAWARE_PROXY_SPANNINGTREE
1138 #ifdef USE_NODEPATCHMGR
1139 npid = numNodeChild;
1140 pids = nodeChildren;
1147 pidsPreAllocated = 0;
1149 for (
int i=0; i<nChild; i++) pids[i] = child[i];
1156 int pdMsgPLLen = p.size();
1157 int pdMsgAvgPLLen = 0;
1164 pdMsgVLLen =
v.
size();
1169 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1174 int lcpoTypeLen = 0;
1175 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1179 int pdMsgPLExtLen = 0;
1180 if(doMigration || isNewProxyAdded) {
1184 int cudaAtomLen = 0;
1185 #if defined(NAMD_CUDA)
1187 #elif defined(NAMD_AVXTILES)
1189 cudaAtomLen = (
numAtoms + 15) & (~15);
1190 #elif defined(NAMD_MIC)
1191 #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
1194 cudaAtomLen = (
numAtoms + 15) & (~15);
1198 ProxyDataMsg *nmsg =
new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1204 nmsg->
plLen = pdMsgPLLen;
1212 nmsg->
vlLen = pdMsgVLLen;
1218 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1219 for (
int i = 0; i <
numAtoms * 2; i++) {
1224 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1225 for (
int i = 0; i <
numAtoms; i++) {
1231 if(doMigration || isNewProxyAdded){
1236 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
1237 #ifdef NAMD_AVXTILES
1243 #if NAMD_SeparateWaters != 0
1245 nmsg->numWaterAtoms = numWaterAtoms;
1248 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
1249 nmsg->isFromImmMsgCall = 0;
1252 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
1253 DebugFileTrace *dft = DebugFileTrace::Object();
1255 dft->writeTrace(
"HP::posReady: for HomePatch[%d], sending proxy msg to: ",
patchID);
1256 for(
int i=0; i<npid; i++) {
1257 dft->writeTrace(
"%d ", pids[i]);
1259 dft->writeTrace(
"\n");
1263 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1264 if (isProxyChanged || localphs == NULL)
1269 for (
int i=0; i<nphs; i++) {
1270 CmiDestoryPersistent(localphs[i]);
1274 localphs =
new PersistentHandle[npid];
1276 for (
int i=0; i<npid; i++) {
1277 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
1279 localphs[i] = CmiCreateNodePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1282 localphs[i] = CmiCreatePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1286 CmiAssert(nphs == npid && localphs != NULL);
1287 CmiUsePersistentHandle(localphs, nphs);
1289 if(doMigration || isNewProxyAdded) {
1294 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1295 CmiUsePersistentHandle(NULL, 0);
1297 isNewProxyAdded = 0;
1300 if(!pidsPreAllocated)
delete [] pids;
1301 DebugM(4,
"patchID("<<
patchID<<
") doing positions Ready\n");
1303 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
1304 positionPtrBegin = p.begin();
1305 positionPtrEnd = p.end();
1332 replacementForces =
f;
1338 for (
int i = 0; i <
numAtoms; ++i )
1340 f_saved[ftag][i] =
f[ftag][i];
1345 #undef DEBUG_REDISTRIB_FORCE
1346 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
1363 void HomePatch::redistrib_colinear_lp_force(
1367 BigReal distance = distance_f;
1372 distance *= r_jk_rlength;
1373 BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
1374 fj += (1. + scale + distance)*fi - r_jk*fdot;
1375 fk -= (scale + distance)*fi + r_jk*fdot;
1403 #define FIX_FOR_WATER
1404 void HomePatch::redistrib_relative_lp_force(
1407 Tensor *virial,
int midpt) {
1408 #ifdef DEBUG_REDISTRIB_FORCE
1410 foldnet = fi + fj + fk + fl;
1413 Vector fja(0), fka(0), fla(0);
1418 BigReal fdot = (fi*r) * invr2;
1425 s = rj - 0.5*(rk + rl);
1436 #if !defined(FIX_FOR_WATER)
1445 #if !defined(FIX_FOR_WATER)
1448 Vector ft = fi - fr - fp;
1458 fpdot = (fi*
p) / sqrt(p2);
1466 ft =
cross(s,v) * invs2;
1477 BigReal srdot = (s*r) * invs2;
1480 BigReal stdot = (s*t) * invs2;
1483 BigReal fact = rrdot*fpdot*invtt*invq;
1487 fja += fp*(1+srdot) + fq*stdot;
1489 ft = fq*(1+stdot) + fp*srdot;
1501 va +=
outer(fka,rk);
1502 va +=
outer(fla,rl);
1512 #ifdef DEBUG_REDISTRIB_FORCE
1513 #define TOL_REDISTRIB 1e-4
1515 fnewnet = fi + fj + fk + fl;
1517 Vector fdiff = fnewnet - foldnet;
1518 Vector tdiff = tnewnet - toldnet;
1519 if (fdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1520 printf(
"Error: force redistribution for water exceeded tolerance: "
1521 "fdiff=(%f, %f, %f)\n", fdiff.
x, fdiff.
y, fdiff.
z);
1523 if (tdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1524 printf(
"Error: torque redistribution for water exceeded tolerance: "
1525 "tdiff=(%f, %f, %f)\n", tdiff.
x, tdiff.
y, tdiff.
z);
1530 void HomePatch::redistrib_ap_force(
Vector& fi,
Vector& fj) {
1535 fi = fi_old + fj_old;
1536 fj = fi_old + fj_old;
1548 void HomePatch::redistrib_lp_water_force(
1553 #ifdef DEBUG_REDISTRIB_FORCE
1557 Vector totforce, tottorque;
1558 totforce = f_ox + f_h1 + f_h2 + f_lp;
1559 tottorque =
cross(f_ox, p_ox) +
cross(f_h1, p_h1) +
cross(f_h2, p_h2);
1562 tottorque +=
cross(f_lp, p_lp);
1568 Vector fad_ox(0), fad_h(0);
1571 Vector r_ox_lp = p_lp - p_ox;
1573 invlen2_r_ox_lp *= invlen2_r_ox_lp;
1574 BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
1575 Vector f_rad = r_ox_lp * rad_factor;
1580 Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
1581 Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5;
1594 Vector f_ang = f_lp - f_rad;
1598 BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
1599 BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
1600 BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
1602 fad_ox += (f_ang * oxcomp);
1603 fad_h += (f_ang * hydcomp);
1608 vir +=
outer(fad_h, p_h1);
1609 vir +=
outer(fad_h, p_h2);
1610 vir -=
outer(f_lp, p_lp);
1620 #ifdef DEBUG_REDISTRIB_FORCE
1622 Vector newforce, newtorque;
1623 newforce = f_ox + f_h1 + f_h2;
1624 newtorque =
cross(f_ox, p_ox) +
cross(f_h1, p_h1) +
cross(f_h2, p_h2);
1625 Vector fdiff = newforce - totforce;
1626 Vector tdiff = newtorque - tottorque;
1628 if (error > 0.0001) {
1629 printf(
"Error: Force redistribution for water "
1630 "exceeded force tolerance: error=%f\n", error);
1632 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1633 printf(
"Error in net force: %f\n", error);
1637 if (error > 0.0001) {
1638 printf(
"Error: Force redistribution for water "
1639 "exceeded torque tolerance: error=%f\n", error);
1641 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1642 printf(
"Error in net torque: %f\n", error);
1665 void HomePatch::reposition_colinear_lonepair(
1669 BigReal distance = distance_f;
1673 if (r2 < 1e-10 || 100. < r2) {
1674 iout <<
iWARN <<
"Large/small distance between lonepair reference atoms: ("
1675 << rj <<
") (" << rk <<
")\n" <<
endi;
1677 ri = rj + (scale + distance*r_jk.
rlength())*r_jk;
1694 void HomePatch::reposition_relative_lonepair(
1698 if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
1699 iout <<
iWARN <<
"Large distance between lonepair reference atoms: ("
1700 << rj <<
") (" << rk <<
") (" << rl <<
")\n" <<
endi;
1702 BigReal r, t,
p, cst, snt, csp, snp, invlen;
1705 if (distance >= 0) {
1731 ri.
x = rj.x + w.
x*a.
x + w.
y*b.
x + w.
z*c.
x;
1732 ri.
y = rj.y + w.
x*a.
y + w.
y*b.
y + w.
z*c.
y;
1733 ri.
z = rj.z + w.
x*a.
z + w.
y*b.
z + w.
z*c.
z;
1747 ri.
x = (mi * ri_old.
x + mj * rj_old.
x)/(mi + mj);
1748 ri.
y = (mi * ri_old.
y + mj * rj_old.
y)/(mi + mj);
1749 ri.
z = (mi * ri_old.
z + mj * rj_old.
z)/(mi + mj);
1755 void HomePatch::reposition_all_lonepairs(
void) {
1758 if (atom[i].mass < 0.01) {
1764 sprintf(errmsg,
"reposition lone pairs: "
1765 "no Lphost exists for LP %d\n", aid);
1773 sprintf(errmsg,
"reposition lone pairs: "
1774 "LP %d has some Lphost atom off patch\n", aid);
1779 reposition_colinear_lonepair(atom[i].position, atom[j.
index].position,
1783 reposition_relative_lonepair(atom[i].position, atom[j.
index].position,
1784 atom[k.
index].position, atom[l.
index].position,
1791 void HomePatch::reposition_all_alchpairs(
void) {
1795 for (
int i = 0; i <
numAtoms; i++) {
1797 alchPair_id = atom[i].id + numFepInitial;
1799 reposition_alchpair(atom[i].position, atom[j.
index].position, atom[i].mass, atom[j.
index].mass);
1809 pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
1813 vel[2] = (pos[2] - ref[2]) * invdt;
1831 pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
1836 vel[3] = (pos[3] - ref[3]) * invdt;
1843 void HomePatch::redistrib_lonepair_forces(
const int ftag,
Tensor *virial) {
1846 for (
int i = 0; i <
numAtoms; i++) {
1847 if (atom[i].mass < 0.01) {
1853 sprintf(errmsg,
"redistrib lone pair forces: "
1854 "no Lphost exists for LP %d\n", aid);
1862 sprintf(errmsg,
"redistrib lone pair forces: "
1863 "LP %d has some Lphost atom off patch\n", aid);
1868 redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
1869 f_mod[ftag][k.
index], atom[i].position, atom[j.
index].position,
1874 redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
1876 atom[i].position, atom[j.
index].position,
1877 atom[k.
index].position, atom[l.
index].position, virial, midpt);
1883 void HomePatch::redistrib_alchpair_forces(
const int ftag) {
1890 for (
int i = 0; i <
numAtoms; i++) {
1892 alchPair_id = atom[i].id + numFepInitial;
1894 redistrib_ap_force(f_mod[ftag][i],f_mod[ftag][j.
index]);
1899 void HomePatch::redistrib_swm4_forces(
const int ftag,
Tensor *virial) {
1903 for (
int i = 0; i <
numAtoms; i++) {
1904 if (atom[i].mass < 0.01) {
1906 redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
1907 f_mod[ftag][i+2], f_mod[ftag][i],
1908 atom[i-2].position, atom[i+1].position,
1909 atom[i+2].position, atom[i].position, virial);
1914 void HomePatch::redistrib_tip4p_forces(
const int ftag,
Tensor* virial) {
1921 if (atom[i].mass < 0.01) {
1923 redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
1924 f_mod[ftag][i-1], f_mod[ftag][i],
1925 atom[i-3].position, atom[i-2].position,
1926 atom[i-1].position, atom[i].position, virial);
1934 const Force * __restrict force_arr,
1942 if ( atom_arr[i].atomFixed ) {
1943 atom_arr[i].velocity = 0;
1945 BigReal dt_mass = dt * atom_arr[i].recipMass;
1946 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1947 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1948 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1953 BigReal dt_mass = dt * atom_arr[i].recipMass;
1954 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1955 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1956 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1963 const Force * __restrict force_arr1,
1964 const Force * __restrict force_arr2,
1965 const Force * __restrict force_arr3,
1975 if ( atom_arr[i].atomFixed ) {
1976 atom_arr[i].velocity = 0;
1978 BigReal rmass = atom_arr[i].recipMass;
1979 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1980 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1981 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1982 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1983 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1984 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1989 BigReal rmass = atom_arr[i].recipMass;
1990 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1991 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1992 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1993 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1994 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1995 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
2008 if ( ! atom_arr[i].atomFixed ) {
2009 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2010 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2011 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2016 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2017 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2018 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2031 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2033 int dieOnError = simParams->
rigidDie;
2035 BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
2041 double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
2042 Vector v_ab, vb_1, vp_1, vb_2, vp_2, new_vel_a, new_vel_b, new_pos_a, new_pos_b, *new_pos, *new_vel;
2043 double dot_v_r_1, dot_v_r_2;
2044 double vb_cm, dr_a, dr_b;
2056 r_wall_SQ = r_wall*r_wall;
2059 if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) {
2063 v_ab = atom[ib].position - atom[ia].position;
2064 rab_SQ = v_ab.
x*v_ab.
x + v_ab.
y*v_ab.
y + v_ab.
z*v_ab.
z;
2066 if (rab_SQ > r_wall_SQ) {
2068 if ( (rab > (2.0*r_wall)) && dieOnError ) {
2070 <<
"The drude is too far away from atom "
2071 << (atom[ia].id + 1) <<
" d = " << rab <<
"!\n" <<
endi;
2079 if ( fixedAtomsOn && atom[ia].atomFixed ) {
2080 if (atom[ib].atomFixed) {
2084 dot_v_r_2 = atom[ib].velocity.x*v_ab.
x
2085 + atom[ib].velocity.y*v_ab.
y + atom[ib].velocity.z*v_ab.
z;
2086 vb_2 = dot_v_r_2 * v_ab;
2087 vp_2 = atom[ib].velocity - vb_2;
2090 if(dot_v_r_2 == 0.0) {
2094 delta_T = dr/fabs(dot_v_r_2);
2095 if(delta_T > maxtime ) delta_T = maxtime;
2098 dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
2100 vb_2 = dot_v_r_2 * v_ab;
2102 new_vel_a = atom[ia].velocity;
2103 new_vel_b = vp_2 + vb_2;
2105 dr_b = -dr + delta_T*dot_v_r_2;
2107 new_pos_a = atom[ia].position;
2108 new_pos_b = atom[ib].position + dr_b*v_ab;
2112 mass_a = atom[ia].mass;
2113 mass_b = atom[ib].mass;
2114 mass_sum = mass_a+mass_b;
2116 dot_v_r_1 = atom[ia].velocity.x*v_ab.
x
2117 + atom[ia].velocity.y*v_ab.
y + atom[ia].velocity.z*v_ab.
z;
2118 vb_1 = dot_v_r_1 * v_ab;
2119 vp_1 = atom[ia].velocity - vb_1;
2121 dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2122 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2123 vb_2 = dot_v_r_2 * v_ab;
2124 vp_2 = atom[ib].velocity - vb_2;
2126 vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
2133 if(dot_v_r_2 == dot_v_r_1) {
2137 delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);
2138 if(delta_T > maxtime ) delta_T = maxtime;
2142 v_Bond = sqrt(kbt/mass_b);
2145 dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
2146 dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
2148 dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
2149 dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
2151 new_pos_a = atom[ia].position + dr_a*v_ab;
2152 new_pos_b = atom[ib].position + dr_b*v_ab;
2159 vb_1 = dot_v_r_1 * v_ab;
2160 vb_2 = dot_v_r_2 * v_ab;
2162 new_vel_a = vp_1 + vb_1;
2163 new_vel_b = vp_2 + vb_2;
2168 atom[ia].position = new_pos_a;
2169 atom[ib].position = new_pos_b;
2171 else if ( virial == 0 ) {
2172 atom[ia].velocity = new_vel_a;
2173 atom[ib].velocity = new_vel_b;
2176 for ( j = 0; j < 2; j++ ) {
2179 new_pos = &new_pos_a;
2180 new_vel = &new_vel_a;
2184 new_pos = &new_pos_b;
2185 new_vel = &new_vel_b;
2187 Force df = (*new_vel - atom[Idx].velocity) *
2188 ( atom[Idx].mass * invdt );
2191 atom[Idx].velocity = *new_vel;
2192 atom[Idx].position = *new_pos;
2197 int partition = atom[Idx].partition;
2198 int slab = (int)floor((z-zmin)*idz);
2199 if (slab < 0) slab += nslabs;
2200 else if (slab >= nslabs) slab -= nslabs;
2203 ppreduction->
item(ppoffset ) += vir.
xx;
2204 ppreduction->
item(ppoffset+1) += vir.
yy;
2205 ppreduction->
item(ppoffset+2) += vir.
zz;
2224 if ( dt && virial ) *virial += wc;
2233 const int useSettle = simParams->
useSettle;
2241 int watmodel = simParams->
watmodel;
2242 if (watmodel ==
WAT_TIP4) wathgsize = 4;
2243 else if (watmodel ==
WAT_SWM4) wathgsize = 5;
2250 if ( ! settle_initialized ) {
2251 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2253 if (atom[ig].rigidBondLength > 0) {
2263 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2270 atom[ig].rigidBondLength,
2271 atom[oatm].rigidBondLength,
2272 settle_mOrmT, settle_mHrmT, settle_ra,
2273 settle_rb, settle_rc, settle_rra);
2274 settle_initialized = 1;
2289 int posRattleParam = 0;
2296 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2297 int hgs = atom[ig].hydrogenGroupSize;
2304 for (
int i = 0; i < hgs; ++i ) {
2305 ref[i] = atom[ig+i].position;
2306 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2307 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2314 BigReal tmp = atom[ig].rigidBondLength;
2316 if (hgs != wathgsize) {
2318 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms "
2319 "but the specified water model requires %d atoms.\n",
2320 atom[ig].
id+1, hgs, wathgsize);
2324 if (useSettle && !anyfixed) {
2329 if ( !(fixed[1] && fixed[2]) ) {
2330 dsq[icnt] = tmp * tmp;
2336 for (
int i = 1; i < hgs; ++i ) {
2337 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2338 if ( !(fixed[0] && fixed[i]) ) {
2339 dsq[icnt] = tmp * tmp;
2353 rattleListElem.
ig = ig;
2354 rattleListElem.
icnt = icnt;
2356 for (
int i = 0; i < icnt; ++i ) {
2360 rattleParamElem.
ia = a;
2361 rattleParamElem.
ib = b;
2362 rattleParamElem.
dsq = dsq[i];
2363 rattleParamElem.
rma = rmass[a];
2364 rattleParamElem.
rmb = rmass[b];
2372 for (
int ig = 0; ig <
numAtoms; ++ig ) {
2373 Force df = (
velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
2377 atom[ig].velocity =
velNew[ig];
2387 return rattle1old(timestep, virial, ppreduction);
2396 const int useSettle = simParams->
useSettle;
2398 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2401 int dieOnError = simParams->
rigidDie;
2409 for (
int j=0;j < n;j+=2) {
2412 for (
int i = 0; i < 3; ++i ) {
2413 ref[i] = atom[ig+i].position;
2414 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2417 for (
int i = 0; i < 3; ++i ) {
2418 ref[i+3] = atom[ig+i].position;
2419 pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
2422 settle_mOrmT, settle_mHrmT, settle_ra,
2423 settle_rb, settle_rc, settle_rra);
2426 for (
int i = 0; i < 3; ++i ) {
2427 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2431 for (
int i = 0; i < 3; ++i ) {
2432 velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
2440 for (
int i = 0; i < 3; ++i ) {
2441 ref[i] = atom[ig+i].position;
2442 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2445 settle_mOrmT, settle_mHrmT, settle_ra,
2446 settle_rb, settle_rc, settle_rra);
2447 for (
int i = 0; i < 3; ++i ) {
2448 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2466 int hgs = atom[ig].hydrogenGroupSize;
2467 for (
int i = 0; i < hgs; ++i ) {
2468 ref[i] = atom[ig+i].position;
2469 pos[i] = atom[ig+i].position;
2470 if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
2471 pos[i] += atom[ig+i].velocity * dt;
2500 for (
int i = 0; i < hgs; ++i ) {
2506 for (
int i = 0; i < hgs; ++i ) {
2507 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2511 if ( consFailure ) {
2513 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom "
2514 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2517 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom "
2518 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2520 }
else if ( ! done ) {
2522 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom "
2523 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2526 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom "
2527 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2535 int hgs = atom[ig].hydrogenGroupSize;
2536 for (
int i = 0; i < hgs; ++i ) {
2537 velNew[ig+i] = atom[ig+i].velocity;
2538 posNew[ig+i] = atom[ig+i].position;
2543 for (
int ig = 0; ig <
numAtoms; ++ig )
2544 atom[ig].position =
posNew[ig];
2545 }
else if ( virial == 0 ) {
2546 for (
int ig = 0; ig <
numAtoms; ++ig )
2547 atom[ig].velocity =
velNew[ig];
2564 const int useSettle = simParams->
useSettle;
2566 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2569 int dieOnError = simParams->
rigidDie;
2572 int ial[10], ibl[10];
2589 if ( ! settle_initialized ) {
2590 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2592 if (atom[ig].rigidBondLength > 0) {
2602 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2609 atom[ig].rigidBondLength,
2610 atom[oatm].rigidBondLength,
2611 settle_mOrmT, settle_mHrmT, settle_ra,
2612 settle_rb, settle_rc, settle_rra);
2613 settle_initialized = 1;
2627 int watmodel = simParams->
watmodel;
2628 if (watmodel ==
WAT_TIP4) wathgsize = 4;
2629 else if (watmodel ==
WAT_SWM4) wathgsize = 5;
2631 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2632 int hgs = atom[ig].hydrogenGroupSize;
2633 if ( hgs == 1 )
continue;
2636 for ( i = 0; i < hgs; ++i ) {
2637 ref[i] = atom[ig+i].position;
2638 pos[i] = atom[ig+i].position;
2639 vel[i] = atom[ig+i].velocity;
2640 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2642 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2645 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
else pos[i] += vel[i] * dt;
2648 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
2649 if (hgs != wathgsize) {
2651 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms "
2652 "but the specified water model requires %d atoms.\n",
2653 atom[ig].
id+1, hgs, wathgsize);
2657 if (useSettle && !anyfixed) {
2668 settle1(ref+2, pos+2, vel+2, invdt,
2669 settle_mOrmT, settle_mHrmT, settle_ra,
2670 settle_rb, settle_rc, settle_rra);
2678 swm4_omrepos(ref, pos, vel, invdt);
2682 settle_mOrmT, settle_mHrmT, settle_ra,
2683 settle_rb, settle_rc, settle_rra);
2685 tip4_omrepos(ref, pos, vel, invdt);
2692 if ( invdt == 0 )
for ( i = 0; i < wathgsize; ++i ) {
2693 atom[ig+i].position = pos[i];
2694 }
else if ( virial == 0 )
for ( i = 0; i < wathgsize; ++i ) {
2695 atom[ig+i].velocity = vel[i];
2696 }
else for ( i = 0; i < wathgsize; ++i ) {
2697 Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
2701 atom[ig+i].velocity = vel[i];
2707 partition = atom[ig].partition;
2708 int slab = (int)floor((z-zmin)*idz);
2709 if (slab < 0) slab += nslabs;
2710 else if (slab >= nslabs) slab -= nslabs;
2713 ppreduction->
item(ppoffset ) += vir.
xx;
2714 ppreduction->
item(ppoffset+1) += vir.
yy;
2715 ppreduction->
item(ppoffset+2) += vir.
zz;
2720 if ( !(fixed[1] && fixed[2]) ) {
2721 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2724 for ( i = 1; i < hgs; ++i ) {
2725 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2726 if ( !(fixed[0] && fixed[i]) ) {
2727 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
2731 if ( icnt == 0 )
continue;
2732 for ( i = 0; i < icnt; ++i ) {
2733 refab[i] = ref[ial[i]] - ref[ibl[i]];
2735 for ( i = 0; i < hgs; ++i ) {
2740 for ( iter = 0; iter < maxiter; ++iter ) {
2744 for ( i = 0; i < icnt; ++i ) {
2745 int a = ial[i];
int b = ibl[i];
2746 Vector pab = pos[a] - pos[b];
2749 BigReal diffsq = rabsq - pabsq;
2750 if ( fabs(diffsq) > (rabsq * tol2) ) {
2753 if ( rpab < ( rabsq * 1.0e-6 ) ) {
2760 BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
2764 if ( invdt != 0. ) {
2775 if ( consFailure ) {
2777 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom "
2778 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2781 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom "
2782 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2784 }
else if ( ! done ) {
2786 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom "
2787 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2790 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom "
2791 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2797 if ( invdt == 0 )
for ( i = 0; i < hgs; ++i ) {
2798 atom[ig+i].position = pos[i];
2799 }
else if ( virial == 0 )
for ( i = 0; i < hgs; ++i ) {
2800 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2801 }
else for ( i = 0; i < hgs; ++i ) {
2802 Force df = netdp[i] * invdt;
2806 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2810 int partition = atom[ig].partition;
2811 int slab = (int)floor((z-zmin)*idz);
2812 if (slab < 0) slab += nslabs;
2813 else if (slab >= nslabs) slab -= nslabs;
2816 ppreduction->
item(ppoffset ) += vir.
xx;
2817 ppreduction->
item(ppoffset+1) += vir.
yy;
2818 ppreduction->
item(ppoffset+2) += vir.
zz;
2822 if ( dt && virial ) *virial += wc;
2833 const int useSettle = simParams->
useSettle;
2838 int dieOnError = simParams->
rigidDie;
2841 int ial[10], ibl[10];
2855 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2857 int hgs = atom[ig].hydrogenGroupSize;
2858 if ( hgs == 1 )
continue;
2861 for ( i = 0; i < hgs; ++i ) {
2862 ref[i] = atom[ig+i].position;
2863 vel[i] = atom[ig+i].velocity;
2864 rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
2865 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2866 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2869 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
2870 if (hgs != wathgsize) {
2871 NAMD_bug(
"Hydrogen group error caught in rattle2().");
2874 if (useSettle && !anyfixed) {
2883 settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
2890 settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
2895 for (i=0; i<hgs; i++) {
2896 atom[ig+i].velocity = vel[i];
2900 if ( !(fixed[1] && fixed[2]) ) {
2901 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2902 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2906 for ( i = 1; i < hgs; ++i ) {
2907 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2908 if ( !(fixed[0] && fixed[i]) ) {
2909 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2910 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
2911 ibl[icnt] = i; ++icnt;
2915 if ( icnt == 0 )
continue;
2917 for ( i = 0; i < icnt; ++i ) {
2918 refab[i] = ref[ial[i]] - ref[ibl[i]];
2922 for ( iter = 0; iter < maxiter; ++iter ) {
2924 for ( i = 0; i < icnt; ++i ) {
2925 int a = ial[i];
int b = ibl[i];
2926 Vector vab = vel[a] - vel[b];
2930 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2931 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2932 wc +=
outer(dp,rab);
2933 vel[a] += rmass[a] * dp;
2934 vel[b] -= rmass[b] * dp;
2943 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
2946 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
2950 for ( i = 0; i < hgs; ++i ) {
2951 atom[ig+i].velocity = vel[i];
2956 *virial += wc / ( 0.5 * dt );
2968 const int useSettle = simParams->
useSettle;
2973 int dieOnError = simParams->
rigidDie;
2976 int ial[10], ibl[10];
2990 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2992 int hgs = atom[ig].hydrogenGroupSize;
2993 if ( hgs == 1 )
continue;
2996 for ( i = 0; i < hgs; ++i ) {
2997 ref[i] = atom[ig+i].position;
2998 vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
3000 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
3001 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
3004 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
3005 if (hgs != wathgsize) {
3006 NAMD_bug(
"Hydrogen group error caught in rattle2().");
3009 if (useSettle && !anyfixed) {
3018 settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
3025 settle2(1.0, 1.0, ref, vel, dt, virial);
3030 for (i=0; i<hgs; i++) {
3031 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3035 if ( !(fixed[1] && fixed[2]) ) {
3036 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
3037 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3041 for ( i = 1; i < hgs; ++i ) {
3042 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3043 if ( !(fixed[0] && fixed[i]) ) {
3044 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
3045 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
3046 ibl[icnt] = i; ++icnt;
3050 if ( icnt == 0 )
continue;
3052 for ( i = 0; i < icnt; ++i ) {
3053 refab[i] = ref[ial[i]] - ref[ibl[i]];
3057 for ( iter = 0; iter < maxiter; ++iter ) {
3059 for ( i = 0; i < icnt; ++i ) {
3060 int a = ial[i];
int b = ibl[i];
3061 Vector vab = vel[a] - vel[b];
3065 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
3066 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
3067 wc +=
outer(dp,rab);
3068 vel[a] += rmass[a] * dp;
3069 vel[b] -= rmass[b] * dp;
3078 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
3081 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
3085 for ( i = 0; i < hgs; ++i ) {
3086 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3091 *virial += wc / ( 0.5 * dt );
3099 DebugM(2,
"loweAndersenVelocities\n");
3103 for (
int i = 0; i <
numAtoms; ++i) {
3106 v[i].position = atom[i].velocity;
3107 v[i].charge = atom[i].mass;
3109 DebugM(2,
"loweAndersenVelocities\n");
3114 DebugM(2,
"loweAndersenFinish\n");
3125 for (
int i = 0; i <
numAtoms; i++) {
3137 for (
int i = 0; i <
numAtoms; i++) {
3140 intRad[2*i+0] = rad - offset;
3141 intRad[2*i+1] = screen*(rad - offset);
3158 for (
int i = 0; i <
numAtoms; i++) {
3160 rhoi = rhoi0+coulomb_radius_offset;
3186 BigReal fij, expkappa, Dij, dEdai, dedasum;
3187 BigReal rhoi, rhoi0, psii, nbetapsi;
3188 BigReal gammapsi2, tanhi, daidr;
3189 for (
int i = 0; i <
numAtoms; i++) {
3193 expkappa = exp(-kappa*fij);
3194 Dij = epsilon_p_i - expkappa*epsilon_s_i;
3196 dEdai = -0.5*
COULOMB*atom[i].charge*atom[i].charge
3197 *(kappa*epsilon_s_i*expkappa-Dij/fij)/
bornRad[i];
3202 rhoi = rhoi0+coulomb_radius_offset;
3204 nbetapsi = -beta*psii;
3205 gammapsi2 = gamma*psii*psii;
3206 tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
3208 * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
3219 if (proxy.
size() > 0) {
3220 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3221 for (
int i = 0; i < proxy.
size(); i++) {
3222 int node = proxy[i];
3231 cp[node].recvData(msg);
3239 if (proxy.
size() > 0) {
3240 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3243 for (
int i = 0; i < proxy.
size(); i++) {
3244 int node = proxy[i];
3254 cp[node].recvData(msg);
3263 for (
int i = 0; i < msg->
psiSumLen; ++i ) {
3270 if (phase1BoxClosedCalled ==
true &&
3271 numGBISP1Arrived==proxy.
size() ) {
3277 numGBISP1Arrived == proxy.
size() &&
3278 numGBISP2Arrived == proxy.
size() &&
3279 numGBISP3Arrived == proxy.
size()) {
3289 for (
int i = 0; i < msg->
dEdaSumLen; ++i ) {
3296 if (phase2BoxClosedCalled ==
true &&
3297 numGBISP2Arrived==proxy.
size() ) {
3303 numGBISP1Arrived == proxy.
size() &&
3304 numGBISP2Arrived == proxy.
size() &&
3305 numGBISP3Arrived == proxy.
size()) {
3334 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3335 int hgs = atom[ig].hydrogenGroupSize;
3336 if ( hgs == 1 )
continue;
3337 for ( i = 0; i < hgs; ++i ) {
3338 ref[i] = atom[ig+i].position;
3339 rmass[i] = 1. / atom[ig+i].mass;
3340 fixed[i] = ( simParams->
fixedAtomsOn && atom[ig+i].atomFixed );
3341 if ( fixed[i] ) rmass[i] = 0.;
3346 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
3348 NAMD_die(
"Hydrogen group error caught in mollyAverage(). It's a bug!\n");
3350 if ( !(fixed[1] && fixed[2]) ) {
3351 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3354 for ( i = 1; i < hgs; ++i ) {
3355 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3356 if ( !(fixed[0] && fixed[i]) ) {
3357 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3361 if ( icnt == 0 )
continue;
3363 molly_lambda.
resize(numLambdas);
3364 lambda = &(molly_lambda[numLambdas - icnt]);
3365 for ( i = 0; i < icnt; ++i ) {
3366 refab[i] = ref[ial[i]] - ref[ibl[i]];
3369 iter=
average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
3370 if ( iter == maxiter ) {
3371 iout <<
iWARN <<
"Exceeded maximum number of iterations in mollyAverage().\n"<<
endi;
3402 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3403 int hgs = atom[ig].hydrogenGroupSize;
3404 if (hgs == 1 )
continue;
3405 for ( i = 0; i < hgs; ++i ) {
3406 ref[i] = atom[ig+i].position;
3408 rmass[i] = 1. / atom[ig+i].mass;
3409 fixed[i] = ( simParams->
fixedAtomsOn && atom[ig+i].atomFixed );
3410 if ( fixed[i] ) rmass[i] = 0.;
3414 if ( atom[ig].rigidBondLength > 0 ) {
3416 NAMD_die(
"Hydrogen group error caught in mollyMollify(). It's a bug!\n");
3418 if ( !(fixed[1] && fixed[2]) ) {
3419 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3422 for ( i = 1; i < hgs; ++i ) {
3423 if ( atom[ig+i].rigidBondLength > 0 ) {
3424 if ( !(fixed[0] && fixed[i]) ) {
3425 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3430 if ( icnt == 0 )
continue;
3431 lambda = &(molly_lambda[numLambdas]);
3433 for ( i = 0; i < icnt; ++i ) {
3434 refab[i] = ref[ial[i]] - ref[ibl[i]];
3437 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
3439 for ( i = 0; i < hgs; ++i ) {
3450 checkpoint_atom.
copy(atom);
3454 #if NAMD_SeparateWaters != 0
3455 checkpoint_numWaterAtoms = numWaterAtoms;
3462 atom.
copy(checkpoint_atom);
3463 numAtoms = atom.
size();
3466 doAtomUpdate =
true;
3472 #if NAMD_SeparateWaters != 0
3473 numWaterAtoms = checkpoint_numWaterAtoms;
3488 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
3498 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
3510 msg->
replica = CmiMyPartition();
3520 NAMD_bug(
"HomePatch::recvCheckpointLoad called during checkpointFree.");
3522 if ( msg->
replica != remote ) {
3523 NAMD_bug(
"HomePatch::recvCheckpointLoad message from wrong replica.");
3527 strcpy(newmsg->
key,key);
3531 newmsg->
pe = CkMyPe();
3532 newmsg->
replica = CmiMyPartition();
3544 doAtomUpdate =
true;
3569 CkpvAccess(_qd)->process();
3601 CkpvAccess(_qd)->process();
3613 CkpvAccess(_qd)->process();
3614 doAtomUpdate =
true;
3631 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
3655 doPairlistCheck_newTolerance *= (1. - simParams->
pairlistShrink);
3656 doPairlistCheck_lattice =
lattice;
3657 doPairlistCheck_positions.
resize(numAtoms);
3659 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
3666 Lattice &lattice_old = doPairlistCheck_lattice;
3669 Vector center_delta = center_cur - center_old;
3673 for ( i=0; i<2; ++i ) {
3674 for (
int j=0; j<2; ++j ) {
3675 for (
int k=0; k<2; ++k ) {
3678 k ? min.
z : max.
z );
3681 corner_delta -= center_delta;
3683 if ( cd > max_cd ) max_cd = cd;
3687 max_cd = sqrt(max_cd);
3692 for ( i=0; i<n; ++i ) {
3694 p_delta -= center_delta;
3696 if ( pd > max_pd ) max_pd = pd;
3698 max_pd = sqrt(max_pd);
3700 BigReal max_tol = max_pd + max_cd;
3707 doPairlistCheck_newTolerance ) ) {
3708 doPairlistCheck_newTolerance *= (1. + simParams->
pairlistGrow);
3711 if ( max_tol > doPairlistCheck_newTolerance ) {
3712 doPairlistCheck_newTolerance = max_tol / (1. - simParams->
pairlistTrigger);
3730 while ( p_i != p_e ) {
3731 const int hgs = p_i->hydrogenGroupSize;
3734 if ( ngs > 5 ) ngs = 5;
3739 for ( i = 1; i < ngs; ++i ) {
3740 p_i[i].nonbondedGroupSize = 0;
3741 BigReal dx = p_i[i].position.x -
x;
3742 BigReal dy = p_i[i].position.y -
y;
3743 BigReal dz = p_i[i].position.z -
z;
3744 BigReal r2 = dx * dx + dy * dy + dz * dz;
3745 if ( r2 > hgcut )
break;
3746 else if ( r2 > maxrad2 ) maxrad2 = r2;
3748 p_i->nonbondedGroupSize = i;
3749 for ( ; i < hgs; ++i ) {
3750 p_i[i].nonbondedGroupSize = 1;
3756 NAMD_bug(
"hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
3773 if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
3774 ( bAwayDist*sysdimb < minSize*0.9999 ) ||
3775 ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
3777 NAMD_die(
"Periodic cell has become too small for original patch grid!\n"
3778 "Possible solutions are to restart from a recent checkpoint,\n"
3779 "increase margin, or disable useFlexibleCell for liquid simulation.");
3784 BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
3785 BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
3786 BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
3788 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
3789 NAMD_die(
"Periodic cell has become too small for original patch grid!\n"
3790 "There are probably many margin violations already reported.\n"
3791 "Possible solutions are to restart from a recent checkpoint,\n"
3792 "increase margin, or disable useFlexibleCell for liquid simulation.");
3802 int xdev, ydev, zdev;
3803 int problemCount = 0;
3807 for ( ; p_i != p_e; ++p_i ) {
3812 if (s.
x < minx) xdev = 0;
3813 else if (maxx <= s.
x) xdev = 2;
3816 if (s.
y < miny) ydev = 0;
3817 else if (maxy <= s.
y) ydev = 2;
3820 if (s.
z < minz) zdev = 0;
3821 else if (maxz <= s.
z) zdev = 2;
3824 if (mInfo[xdev][ydev][zdev]) {
3844 for (i=0; i<numNeighbors; i++) {
3860 int xdev, ydev, zdev;
3867 #if NAMD_SeparateWaters != 0
3869 int numLostWaterAtoms = 0;
3872 while ( atom_i != atom_e ) {
3879 if ( atom_i->migrationGroupSize ) {
3881 if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
3886 int mgs = atom_i->migrationGroupSize;
3888 for (
int j=atom_i->hydrogenGroupSize; j<mgs;
3889 j+=(atom_i+j)->hydrogenGroupSize ) {
3890 pos += (atom_i+j)->position;
3905 if (s.
x < minx) xdev = 0;
3906 else if (maxx <= s.
x) xdev = 2;
3909 if (s.
y < miny) ydev = 0;
3910 else if (maxy <= s.
y) ydev = 2;
3913 if (s.
z < minz) zdev = 0;
3914 else if (maxz <= s.
z) zdev = 2;
3919 if (mInfo[xdev][ydev][zdev]) {
3924 DebugM(3,
"Migrating atom " << atom_i->id <<
" from patch "
3925 <<
patchID <<
" with position " << atom_i->position <<
"\n");
3932 #if NAMD_SeparateWaters != 0
3941 if (atomIndex < numWaterAtoms)
3942 numLostWaterAtoms++;
3950 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
3958 #if NAMD_SeparateWaters != 0
3959 numWaterAtoms -= numLostWaterAtoms;
3963 int delpos = numAtoms - delnum;
3964 DebugM(4,
"numAtoms " << numAtoms <<
" deleted " << delnum <<
"\n");
3965 atom.
del(delpos,delnum);
3967 numAtoms = atom.
size();
3976 DebugM(1,
"Draining migration buffer ("<<i<<
","<<numMlBuf<<
")\n");
3983 if (!allMigrationIn) {
3984 DebugM(3,
"All Migrations NOT in, we are suspending patch "<<
patchID<<
"\n");
3985 migrationSuspended =
true;
3987 migrationSuspended =
false;
3989 allMigrationIn =
false;
4007 #if NAMD_SeparateWaters != 0
4026 for (mi = migrationList.
begin(); mi != migrationList.
end(); mi++) {
4027 DebugM(1,
"Migrating atom " << mi->id <<
" to patch "
4028 <<
patchID <<
" with position " << mi->position <<
"\n");
4029 if ( mi->migrationGroupSize ) {
4030 if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
4032 int mgs = mi->migrationGroupSize;
4034 for (
int j=mi->hydrogenGroupSize; j<mgs;
4035 j+=(mi+j)->hydrogenGroupSize ) {
4036 pos += (mi+j)->position;
4041 mother_transform = mi->transform;
4045 mi->transform = mother_transform;
4047 mi->position =
lattice.
nearest(mi->position,center,&(mi->transform));
4048 mother_transform = mi->transform;
4053 mi->transform = mother_transform;
4060 #endif // if (NAMD_SeparateWaters != 0)
4063 numAtoms = atom.
size();
4066 DebugM(3,
"Counter on " <<
patchID <<
" = " << patchMigrationCounter <<
"\n");
4067 if (!--patchMigrationCounter) {
4068 DebugM(3,
"All Migrations are in for patch "<<
patchID<<
"\n");
4069 allMigrationIn =
true;
4070 patchMigrationCounter = numNeighbors;
4071 if (migrationSuspended) {
4073 migrationSuspended =
false;
4086 #if NAMD_SeparateWaters != 0
4091 void HomePatch::separateAtoms() {
4103 if (atom.
size() < 0)
return;
4106 tempAtom.resize(numAtoms);
4116 int nonWaterIndex = 0;
4117 while (i < numAtoms) {
4123 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4126 if (waterIndex != i) {
4127 atom[waterIndex ] = atom[i ];
4128 atom[waterIndex + 1] = atom[i + 1];
4129 atom[waterIndex + 2] = atom[i + 2];
4130 if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3];
4131 if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4];
4135 waterIndex += wathgsize;
4141 for (
int j = 0; j < hgs; j++)
4142 tempAtom[nonWaterIndex + j] = atom[i + j];
4144 nonWaterIndex += hgs;
4156 for (i = 0; i < nonWaterIndex; i++)
4157 atom[waterIndex + i] = tempAtom[i];
4160 numWaterAtoms = waterIndex;
4173 if (al.
size() <= 0)
return;
4175 const int orig_atomSize = atom.
size();
4176 const int orig_alSize = al.
size();
4179 atom.
resize(orig_atomSize + orig_alSize);
4182 #if 0 // version where non-waters are moved to scratch first
4204 const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
4205 tempAtom.resize(orig_atom_numNonWaters + al.
size());
4206 for (
int i = 0; i < orig_atom_numNonWaters; i++)
4207 tempAtom[i] = atom[numWaterAtoms + i];
4211 int atom_waterIndex = numWaterAtoms;
4212 int atom_nonWaterIndex = orig_atom_numNonWaters;
4214 while (i < orig_alSize) {
4219 NAMD_bug(
"HomePatch::mergeAtomList() not updated for migration groups!");
4223 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4228 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4229 Transform mother_transform = al[i].transform;
4234 al[i+1].transform = mother_transform;
4239 al[i+2].transform = mother_transform;
4242 atom[atom_waterIndex ] = al[i ];
4243 atom[atom_waterIndex + 1] = al[i + 1];
4244 atom[atom_waterIndex + 2] = al[i + 2];
4246 atom_waterIndex += 3;
4254 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4255 Transform mother_transform = al[i].transform;
4258 for (
int j = 1; j < hgs; j++) {
4261 al[i+j].transform = mother_transform;
4265 for (
int j = 0; j < hgs; j++)
4266 tempAtom[atom_nonWaterIndex + j] = al[i + j];
4268 atom_nonWaterIndex += hgs;
4275 for (
int i = 0; i < atom_nonWaterIndex; i++)
4276 atom[atom_waterIndex + i] = tempAtom[i];
4279 numWaterAtoms = atom_waterIndex;
4280 numAtoms = atom.
size();
4303 int al_numWaterAtoms = 0;
4305 while (i < orig_alSize) {
4311 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4312 al_numWaterAtoms += wathgsize;
4320 if (al_numWaterAtoms > 0) {
4321 for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
4322 atom[i + al_numWaterAtoms] = atom[i];
4329 int atom_waterIndex = numWaterAtoms;
4330 int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
4332 while (i < orig_alSize) {
4337 NAMD_bug(
"HomePatch::mergeAtomList() not updated for migration groups!");
4341 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4346 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4347 Transform mother_transform = al[i].transform;
4352 al[i+1].transform = mother_transform;
4357 al[i+2].transform = mother_transform;
4360 atom[atom_waterIndex ] = al[i ];
4361 atom[atom_waterIndex + 1] = al[i + 1];
4362 atom[atom_waterIndex + 2] = al[i + 2];
4364 if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
4366 atom_waterIndex += wathgsize;
4374 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4375 Transform mother_transform = al[i].transform;
4378 for (
int j = 1; j < hgs; j++) {
4381 al[i+j].transform = mother_transform;
4385 for (
int j = 0; j < hgs; j++)
4386 atom[atom_nonWaterIndex + j] = al[i + j];
4388 atom_nonWaterIndex += hgs;
4395 numWaterAtoms = atom_waterIndex;
4396 numAtoms = atom_nonWaterIndex;
4416 for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
4420 for (i=n-1;i>=0;i--) {
4422 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
4432 double big,dum,sum,temp;
4438 if ((temp=fabs(a[i][j])) > big) big=temp;
4439 if (big == 0.0)
NAMD_die(
"Singular matrix in routine ludcmp\n");
4445 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
4452 sum -= a[i][k]*a[k][j];
4454 if ( (dum=vv[i]*fabs(sum)) >= big) {
4469 if (a[j][j] == 0.0) a[j][j]=
TINY;
4472 for (i=j+1;i<n;i++) a[i][j] *= dum;
4483 gqij[i][ial[i]]=2.0*refab[i];
4484 gqij[i][ibl[i]]=-gqij[i][ial[i]];
4490 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) {
4528 G_q(refab,grhs,n,m,ial,ibl);
4529 for (k=1;k<=ntrial;k++) {
4538 multiplier = lambda[i];
4541 auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
4547 tmp[j]+=auxrhs[i][j];
4557 for ( i = 0; i < m; i++ ) {
4568 G_q(avgab,glhs,n,m,ial,ibl);
4572 iout<<
iINFO << glhs[i*n+0] <<
" " << glhs[i*n+1] <<
" " << glhs[i*n+2] << std::endl<<
endi;
4577 for (j=0; j<n; j++) {
4578 for (i=0; i<m;i++) {
4579 grhs2[i][j] = grhs[i][j]*imass[j];
4588 for (k1=0;k1<n;k1++) {
4589 fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
4617 gij[i]=avgab[i]*avgab[i]-length2[i];
4622 for(i=0;i<m-1;i++) {
4635 for (i=0;i<m;i++) errf += fabs(fvec[i]);
4642 for (i=0;i<m;i++) p[i] = -fvec[i];
4656 <<
" " << lambda[1] <<
" " << lambda[2] << std::endl<<
endi;
4659 if (errx <= tolx)
break;
4666 iout<<
iINFO <<
"LAMBDA:" << lambda[0] <<
" " << lambda[1] <<
" " << lambda[2] << std::endl<<
endi;
4676 Vector zero(0.0,0.0,0.0);
4687 tmpforce[i]=imass[i]*force[i];
4693 for ( i = 0; i < m; i++ ) {
4697 G_q(avgab,glhs,n,m,ial,ibl);
4698 G_q(refab,grhs,n,m,ial,ibl);
4700 for (j=0; j<n; j++) {
4701 for (i=0; i<m;i++) {
4702 grhs2[i][j] = grhs[i][j]*imass[j];
4712 fjac[j][i] += glhs[i][k]*grhs2[j][k];
4723 aux[i]+=grhs[i][j]*tmpforce[j];
4733 y[j] += aux[i]*glhs[i][j];
4742 tmpforce2[i]=imass[i]*y[i];
4751 Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
4752 tmpforce[ial[j]] += tmpf;
4753 tmpforce[ibl[j]] -= tmpf;
4758 force[i]=tmpforce[i]+y[i];
void depositMigration(MigrateAtomsMsg *)
void recvCheckpointLoad(CheckpointAtomsMsg *msg)
template void settle1_SIMD< 1 >(const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
char scriptStringArg1[128]
#define NAMD_EVENT_STOP(eon, id)
std::ostream & iINFO(std::ostream &s)
Data * clientOpen(int count=1)
void addVelocityToPosition(FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
int numNodesWithPatches(void)
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
void registerProxy(RegisterProxyMsg *)
BigReal solvent_dielectric
void copy(ResizeArray< Elem > &ra)
OwnerBox< Patch, Results > forceBox
static ProxyMgr * Object()
unsigned int hydrogenGroupSize
SortedArray< LSSSubsDat > & lssSubs(ComputeQMMgr *mgr)
int flLen[Results::maxNumForces]
const int * get_qmAtmIndx()
static BigReal dielectric_1
static void partition(int *order, const FullAtom *atoms, int begin, int end)
void del(int index, int num=1)
static PatchMap * Object()
void sendProxies(int pid, int *list, int n)
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
BigReal max_c(int pid) const
int find(const Elem &e) const
CompAtom * velocityPtrEnd
int32 numhosts
Either 2 or 3 host atoms, depending on LP type.
BigReal min_a(int pid) const
#define NAMD_SeparateWaters
SimParameters * simParameters
void sendNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *)
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
int index_a(int pid) const
MigrationList migrationList
void addForceToMomentum(FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
void rattle2(const BigReal, Tensor *virial)
static __thread float4 * forces
void recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg)
void gbisComputeAfterP2()
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
void sendExchangeReq(int pid, int src)
CompAtom * avgPositionPtrEnd
void receiveResults(ProxyResultVarsizeMsg *msg)
char const *const NamdProfileEventStr[]
ExchangeAtomsMsg * exchange_msg
void positionsReady(int n=0)
#define PROXY_DATA_PRIORITY
std::ostream & iWARN(std::ostream &s)
virtual void boxClosed(int)
if(ComputeNonbondedUtil::goMethod==2)
BigReal HGMatrixBigReal[MAXHGS][MAXHGS]
CompAtom * avgPositionList
void loweAndersenVelocities()
int numaway_b(void) const
BigReal length(void) const
void sortAtomsForCUDA(int *order, const FullAtom *atoms, int nfree, int n)
void recvCheckpointReq(int task, const char *key, int replica, int pe)
void exchangeCheckpoint(int scriptTask, int &bpc)
Position nearest(Position data, ScaledPosition ref) const
template void settle1_SIMD< 2 >(const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
BigReal min_b(int pid) const
BigReal coulomb_radius_offset
void positionsReady(int doMigration=0)
int settle2(BigReal mO, BigReal mH, const Vector *pos, Vector *vel, BigReal dt, Tensor *virial)
std::vector< RattleList > rattleList
void patchLoad(PatchID id, int nAtoms, int timestep)
void sendCheckpointAck(int pid, int dst, int dstpe)
static ProxyNodeAwareSpanningTreeMsg * getANewMsg(PatchID pid, NodeID nid, proxyTreeNode *tree, int size)
int berendsenPressure_count
void reorder(Elem *a, int n)
int periodic_a(void) const
const Real * get_qmAtomGroup() const
BigReal HGArrayBigReal[MAXHGS]
int const * getLcpoParamType()
Lphost * get_lphost(int atomid) const
static int compDistance(const void *a, const void *b)
void unregisterProxy(UnregisterProxyMsg *)
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const float const PatchPairRecord *__restrict__ const int *__restrict__ atomIndex
void NAMD_bug(const char *err_msg)
ResizeArray< FullAtom > atoms
template void rattlePair< 1 >(const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, bool &consFailure)
int oneAwayNeighbors(int pid, PatchID *neighbor_ids=0)
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)
zVector HGMatrixVector[MAXHGS][MAXHGS]
std::map< std::string, checkpoint_t * > checkpoints
void submitLoadStats(int timestep)
void mollyMollify(Tensor *virial)
Bool staticAtomAssignment
int rattle1old(const BigReal, Tensor *virial, SubmitReduction *)
int index_b(int pid) const
void clientClose(int count=1)
void recvSpanningTree(int *t, int n)
LocalID localID(AtomID id)
void setall(const Elem &elem)
void replaceForces(ExtForce *f)
std::vector< int > settleList
void sendMigrationMsgs(PatchID, MigrationInfo *, int)
void gbisComputeAfterP1()
void NAMD_die(const char *err_msg)
void sendNodeAwareSpanningTree()
void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx, HGArrayBigReal &b)
static LdbCoordinator * Object()
CompAtom * avgPositionPtrBegin
int periodic_b(void) const
int berendsenPressure_count
static AtomMap * Object()
MigrateAtomsMsg * msgbuf[PatchMap::MaxOneAway]
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
Tensor outer(const Vector &v1, const Vector &v2)
void saveForce(const int ftag=Results::normal)
void setGBISIntrinsicRadii()
std::vector< Vector > posNew
int berendsenPressure_count
BigReal max_b(int pid) const
int index_c(int pid) const
int add(const Elem &elem)
std::vector< RattleParam > rattleParam
void sendProxyData(ProxyDataMsg *, int, int *)
BigReal max_a(int pid) const
void sendSpanningTree(ProxySpanningTreeMsg *)
void sendCheckpointLoad(CheckpointAtomsMsg *msg, int dst, int dstpe)
BigReal length2(void) const
int settle1(const Vector *ref, Vector *pos, Vector *vel, BigReal invdt, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
optimized settle1 algorithm, reuses water properties as much as possible
void swap(ResizeArray< Elem > &ra)
void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
void sendProxyList(int pid, int *plist, int size)
#define NAMD_EVENT_START_EX(eon, id, str)
Position apply_transform(Position data, const Transform &t) const
void sendProxyAll(ProxyDataMsg *, int, int *)
Position unscale(ScaledPosition s) const
ForceList * forceList[Results::maxNumForces]
void G_q(const HGArrayVector &refab, HGMatrixVector &gqij, const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl)
Position reverse_transform(Position data, const Transform &t) const
OwnerBox< Patch, GBReal > dEdaSumBox
int pid(int aIndex, int bIndex, int cIndex)
static float MassToRadius(Mass mi)
int flLen[Results::maxNumForces]
static float MassToScreen(Mass mi)
void recvExchangeReq(int req)
void settle1init(BigReal pmO, BigReal pmH, BigReal hhdist, BigReal ohdist, BigReal &mOrmT, BigReal &mHrmT, BigReal &ra, BigReal &rb, BigReal &rc, BigReal &rra)
initialize cached water properties
Elem * find(const Elem &elem)
CompAtomExt * positionExtList
void buildSpanningTree(void)
ScaledPosition scale(Position p) const
int numPatchesOnNode(int node)
BigReal pairlistTolerance
std::vector< Vector > velNew
int numaway_c(void) const
void receiveResult(ProxyGBISP1ResultMsg *msg)
void rattleN(const int icnt, const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
std::ostream & iERROR(std::ostream &s)
BigReal min_c(int pid) const
infostream & endi(infostream &s)
#define SET_PRIORITY(MSG, SEQ, PRIO)
void unregisterIDsFullAtom(const FullAtom *begin, const FullAtom *end)
zVector HGArrayVector[MAXHGS]
void addForceToMomentum3(FullAtom *__restrict atom_arr, const Force *__restrict force_arr1, const Force *__restrict force_arr2, const Force *__restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms)
void registerPatch(int patchID, int numPes, int *pes)
std::vector< int > noconstList
ForceList f[Results::maxNumForces]
#define GB2_PROXY_DATA_PRIORITY
int periodic_c(void) const
void loweAndersenFinish()
int numaway_a(void) const
void addRattleForce(const BigReal invdt, Tensor &wc)
static __thread int num_atoms
void sendCheckpointStore(CheckpointAtomsMsg *msg, int dst, int dstpe)
void sendCheckpointReq(int pid, int remote, const char *key, int task)
void useSequencer(Sequencer *sequencerPtr)
static PatchMgr * Object()
OwnerBox< Patch, GBReal > psiSumBox
#define GB3_PROXY_DATA_PRIORITY
void sendExchangeMsg(ExchangeAtomsMsg *msg, int dst, int dstpe)
void recvCheckpointStore(CheckpointAtomsMsg *msg)
void recvExchangeMsg(ExchangeAtomsMsg *msg)
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)
#define PATCH_PRIORITY(PID)
CompAtom * velocityPtrBegin
void exchangeAtoms(int scriptTask)