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) || defined(NAMD_ROCTX_ENABLED)
948 if (doMigration && simParams->
qmLSSOn)
951 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
958 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || 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;
1067 #if ! (defined(NAMD_CUDA) || defined(NAMD_HIP))
1068 #if defined(NAMD_AVXTILES)
1084 for ( i=0; i<n; ++i ) {
1099 numGBISP1Arrived = 0;
1100 phase1BoxClosedCalled =
false;
1101 numGBISP2Arrived = 0;
1102 phase2BoxClosedCalled =
false;
1103 numGBISP3Arrived = 0;
1104 phase3BoxClosedCalled =
false;
1105 if (doMigration || isNewProxyAdded)
1110 if (doMigration || isNewProxyAdded) {
1118 int pidsPreAllocated = 1;
1121 npid = proxy.
size();
1122 pids =
new int[npid];
1123 pidsPreAllocated = 0;
1125 int *pide = pids + proxy.
size();
1126 int patchNodesLast =
1128 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
1138 #ifdef NODEAWARE_PROXY_SPANNINGTREE
1139 #ifdef USE_NODEPATCHMGR
1140 npid = numNodeChild;
1141 pids = nodeChildren;
1148 pidsPreAllocated = 0;
1150 for (
int i=0; i<nChild; i++) pids[i] = child[i];
1157 int pdMsgPLLen = p.size();
1158 int pdMsgAvgPLLen = 0;
1165 pdMsgVLLen =
v.
size();
1170 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1175 int lcpoTypeLen = 0;
1176 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1180 int pdMsgPLExtLen = 0;
1181 if(doMigration || isNewProxyAdded) {
1185 int cudaAtomLen = 0;
1186 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
1188 #elif defined(NAMD_AVXTILES)
1190 cudaAtomLen = (
numAtoms + 15) & (~15);
1191 #elif defined(NAMD_MIC)
1192 #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
1195 cudaAtomLen = (
numAtoms + 15) & (~15);
1199 ProxyDataMsg *nmsg =
new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1205 nmsg->
plLen = pdMsgPLLen;
1213 nmsg->
vlLen = pdMsgVLLen;
1219 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1220 for (
int i = 0; i <
numAtoms * 2; i++) {
1225 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1226 for (
int i = 0; i <
numAtoms; i++) {
1232 if(doMigration || isNewProxyAdded){
1237 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
1238 #ifdef NAMD_AVXTILES
1244 #if NAMD_SeparateWaters != 0
1246 nmsg->numWaterAtoms = numWaterAtoms;
1249 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
1250 nmsg->isFromImmMsgCall = 0;
1253 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
1254 DebugFileTrace *dft = DebugFileTrace::Object();
1256 dft->writeTrace(
"HP::posReady: for HomePatch[%d], sending proxy msg to: ",
patchID);
1257 for(
int i=0; i<npid; i++) {
1258 dft->writeTrace(
"%d ", pids[i]);
1260 dft->writeTrace(
"\n");
1264 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1265 if (isProxyChanged || localphs == NULL)
1270 for (
int i=0; i<nphs; i++) {
1271 CmiDestoryPersistent(localphs[i]);
1275 localphs =
new PersistentHandle[npid];
1277 for (
int i=0; i<npid; i++) {
1278 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
1280 localphs[i] = CmiCreateNodePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1283 localphs[i] = CmiCreatePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1287 CmiAssert(nphs == npid && localphs != NULL);
1288 CmiUsePersistentHandle(localphs, nphs);
1290 if(doMigration || isNewProxyAdded) {
1295 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1296 CmiUsePersistentHandle(NULL, 0);
1298 isNewProxyAdded = 0;
1301 if(!pidsPreAllocated)
delete [] pids;
1302 DebugM(4,
"patchID("<<
patchID<<
") doing positions Ready\n");
1304 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
1305 positionPtrBegin = p.begin();
1306 positionPtrEnd = p.end();
1333 replacementForces =
f;
1339 for (
int i = 0; i <
numAtoms; ++i )
1341 f_saved[ftag][i] =
f[ftag][i];
1346 #undef DEBUG_REDISTRIB_FORCE
1347 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
1364 void HomePatch::redistrib_colinear_lp_force(
1368 BigReal distance = distance_f;
1373 distance *= r_jk_rlength;
1374 BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
1375 fj += (1. + scale + distance)*fi - r_jk*fdot;
1376 fk -= (scale + distance)*fi + r_jk*fdot;
1404 #define FIX_FOR_WATER
1405 void HomePatch::redistrib_relative_lp_force(
1408 Tensor *virial,
int midpt) {
1409 #ifdef DEBUG_REDISTRIB_FORCE
1411 foldnet = fi + fj + fk + fl;
1414 Vector fja(0), fka(0), fla(0);
1419 BigReal fdot = (fi*r) * invr2;
1426 s = rj - 0.5*(rk + rl);
1437 #if !defined(FIX_FOR_WATER)
1446 #if !defined(FIX_FOR_WATER)
1449 Vector ft = fi - fr - fp;
1459 fpdot = (fi*
p) / sqrt(p2);
1467 ft =
cross(s,v) * invs2;
1478 BigReal srdot = (s*r) * invs2;
1481 BigReal stdot = (s*t) * invs2;
1484 BigReal fact = rrdot*fpdot*invtt*invq;
1488 fja += fp*(1+srdot) + fq*stdot;
1490 ft = fq*(1+stdot) + fp*srdot;
1502 va +=
outer(fka,rk);
1503 va +=
outer(fla,rl);
1513 #ifdef DEBUG_REDISTRIB_FORCE
1514 #define TOL_REDISTRIB 1e-4
1516 fnewnet = fi + fj + fk + fl;
1518 Vector fdiff = fnewnet - foldnet;
1519 Vector tdiff = tnewnet - toldnet;
1520 if (fdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1521 printf(
"Error: force redistribution for water exceeded tolerance: "
1522 "fdiff=(%f, %f, %f)\n", fdiff.
x, fdiff.
y, fdiff.
z);
1524 if (tdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1525 printf(
"Error: torque redistribution for water exceeded tolerance: "
1526 "tdiff=(%f, %f, %f)\n", tdiff.
x, tdiff.
y, tdiff.
z);
1531 void HomePatch::redistrib_ap_force(
Vector& fi,
Vector& fj) {
1536 fi = fi_old + fj_old;
1537 fj = fi_old + fj_old;
1549 void HomePatch::redistrib_lp_water_force(
1554 #ifdef DEBUG_REDISTRIB_FORCE
1558 Vector totforce, tottorque;
1559 totforce = f_ox + f_h1 + f_h2 + f_lp;
1560 tottorque =
cross(f_ox, p_ox) +
cross(f_h1, p_h1) +
cross(f_h2, p_h2);
1563 tottorque +=
cross(f_lp, p_lp);
1569 Vector fad_ox(0), fad_h(0);
1572 Vector r_ox_lp = p_lp - p_ox;
1574 invlen2_r_ox_lp *= invlen2_r_ox_lp;
1575 BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
1576 Vector f_rad = r_ox_lp * rad_factor;
1581 Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
1582 Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5;
1595 Vector f_ang = f_lp - f_rad;
1599 BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
1600 BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
1601 BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
1603 fad_ox += (f_ang * oxcomp);
1604 fad_h += (f_ang * hydcomp);
1609 vir +=
outer(fad_h, p_h1);
1610 vir +=
outer(fad_h, p_h2);
1611 vir -=
outer(f_lp, p_lp);
1621 #ifdef DEBUG_REDISTRIB_FORCE
1623 Vector newforce, newtorque;
1624 newforce = f_ox + f_h1 + f_h2;
1625 newtorque =
cross(f_ox, p_ox) +
cross(f_h1, p_h1) +
cross(f_h2, p_h2);
1626 Vector fdiff = newforce - totforce;
1627 Vector tdiff = newtorque - tottorque;
1629 if (error > 0.0001) {
1630 printf(
"Error: Force redistribution for water "
1631 "exceeded force tolerance: error=%f\n", error);
1633 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1634 printf(
"Error in net force: %f\n", error);
1638 if (error > 0.0001) {
1639 printf(
"Error: Force redistribution for water "
1640 "exceeded torque tolerance: error=%f\n", error);
1642 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1643 printf(
"Error in net torque: %f\n", error);
1666 void HomePatch::reposition_colinear_lonepair(
1670 BigReal distance = distance_f;
1674 if (r2 < 1e-10 || 100. < r2) {
1675 iout <<
iWARN <<
"Large/small distance between lonepair reference atoms: ("
1676 << rj <<
") (" << rk <<
")\n" <<
endi;
1678 ri = rj + (scale + distance*r_jk.
rlength())*r_jk;
1695 void HomePatch::reposition_relative_lonepair(
1699 if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
1700 iout <<
iWARN <<
"Large distance between lonepair reference atoms: ("
1701 << rj <<
") (" << rk <<
") (" << rl <<
")\n" <<
endi;
1703 BigReal r, t,
p, cst, snt, csp, snp, invlen;
1706 if (distance >= 0) {
1732 ri.
x = rj.x + w.
x*a.
x + w.
y*b.
x + w.
z*c.
x;
1733 ri.
y = rj.y + w.
x*a.
y + w.
y*b.
y + w.
z*c.
y;
1734 ri.
z = rj.z + w.
x*a.
z + w.
y*b.
z + w.
z*c.
z;
1748 ri.
x = (mi * ri_old.
x + mj * rj_old.
x)/(mi + mj);
1749 ri.
y = (mi * ri_old.
y + mj * rj_old.
y)/(mi + mj);
1750 ri.
z = (mi * ri_old.
z + mj * rj_old.
z)/(mi + mj);
1756 void HomePatch::reposition_all_lonepairs(
void) {
1759 if (atom[i].mass < 0.01) {
1765 sprintf(errmsg,
"reposition lone pairs: "
1766 "no Lphost exists for LP %d\n", aid);
1774 sprintf(errmsg,
"reposition lone pairs: "
1775 "LP %d has some Lphost atom off patch\n", aid);
1780 reposition_colinear_lonepair(atom[i].position, atom[j.
index].position,
1784 reposition_relative_lonepair(atom[i].position, atom[j.
index].position,
1785 atom[k.
index].position, atom[l.
index].position,
1792 void HomePatch::reposition_all_alchpairs(
void) {
1796 for (
int i = 0; i <
numAtoms; i++) {
1798 alchPair_id = atom[i].id + numFepInitial;
1800 reposition_alchpair(atom[i].position, atom[j.
index].position, atom[i].mass, atom[j.
index].mass);
1810 pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
1814 vel[2] = (pos[2] - ref[2]) * invdt;
1832 pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
1837 vel[3] = (pos[3] - ref[3]) * invdt;
1844 void HomePatch::redistrib_lonepair_forces(
const int ftag,
Tensor *virial) {
1847 for (
int i = 0; i <
numAtoms; i++) {
1848 if (atom[i].mass < 0.01) {
1854 sprintf(errmsg,
"redistrib lone pair forces: "
1855 "no Lphost exists for LP %d\n", aid);
1863 sprintf(errmsg,
"redistrib lone pair forces: "
1864 "LP %d has some Lphost atom off patch\n", aid);
1869 redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
1870 f_mod[ftag][k.
index], atom[i].position, atom[j.
index].position,
1875 redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
1877 atom[i].position, atom[j.
index].position,
1878 atom[k.
index].position, atom[l.
index].position, virial, midpt);
1884 void HomePatch::redistrib_alchpair_forces(
const int ftag) {
1891 for (
int i = 0; i <
numAtoms; i++) {
1893 alchPair_id = atom[i].id + numFepInitial;
1895 redistrib_ap_force(f_mod[ftag][i],f_mod[ftag][j.
index]);
1900 void HomePatch::redistrib_swm4_forces(
const int ftag,
Tensor *virial) {
1904 for (
int i = 0; i <
numAtoms; i++) {
1905 if (atom[i].mass < 0.01) {
1907 redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
1908 f_mod[ftag][i+2], f_mod[ftag][i],
1909 atom[i-2].position, atom[i+1].position,
1910 atom[i+2].position, atom[i].position, virial);
1915 void HomePatch::redistrib_tip4p_forces(
const int ftag,
Tensor* virial) {
1922 if (atom[i].mass < 0.01) {
1924 redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
1925 f_mod[ftag][i-1], f_mod[ftag][i],
1926 atom[i-3].position, atom[i-2].position,
1927 atom[i-1].position, atom[i].position, virial);
1935 const Force * __restrict force_arr,
1943 if ( atom_arr[i].atomFixed ) {
1944 atom_arr[i].velocity = 0;
1946 BigReal dt_mass = dt * atom_arr[i].recipMass;
1947 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1948 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1949 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1954 BigReal dt_mass = dt * atom_arr[i].recipMass;
1955 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1956 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1957 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1964 const Force * __restrict force_arr1,
1965 const Force * __restrict force_arr2,
1966 const Force * __restrict force_arr3,
1976 if ( atom_arr[i].atomFixed ) {
1977 atom_arr[i].velocity = 0;
1979 BigReal rmass = atom_arr[i].recipMass;
1980 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1981 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1982 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1983 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1984 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1985 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1990 BigReal rmass = atom_arr[i].recipMass;
1991 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1992 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1993 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1994 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1995 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1996 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
2009 if ( ! atom_arr[i].atomFixed ) {
2010 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2011 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2012 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2017 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2018 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2019 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2032 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2034 int dieOnError = simParams->
rigidDie;
2036 BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
2042 double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
2043 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;
2044 double dot_v_r_1, dot_v_r_2;
2045 double vb_cm, dr_a, dr_b;
2057 r_wall_SQ = r_wall*r_wall;
2060 if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) {
2064 v_ab = atom[ib].position - atom[ia].position;
2065 rab_SQ = v_ab.
x*v_ab.
x + v_ab.
y*v_ab.
y + v_ab.
z*v_ab.
z;
2067 if (rab_SQ > r_wall_SQ) {
2069 if ( (rab > (2.0*r_wall)) && dieOnError ) {
2071 <<
"The drude is too far away from atom "
2072 << (atom[ia].id + 1) <<
" d = " << rab <<
"!\n" <<
endi;
2080 if ( fixedAtomsOn && atom[ia].atomFixed ) {
2081 if (atom[ib].atomFixed) {
2085 dot_v_r_2 = atom[ib].velocity.x*v_ab.
x
2086 + atom[ib].velocity.y*v_ab.
y + atom[ib].velocity.z*v_ab.
z;
2087 vb_2 = dot_v_r_2 * v_ab;
2088 vp_2 = atom[ib].velocity - vb_2;
2091 if(dot_v_r_2 == 0.0) {
2095 delta_T = dr/fabs(dot_v_r_2);
2096 if(delta_T > maxtime ) delta_T = maxtime;
2099 dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
2101 vb_2 = dot_v_r_2 * v_ab;
2103 new_vel_a = atom[ia].velocity;
2104 new_vel_b = vp_2 + vb_2;
2106 dr_b = -dr + delta_T*dot_v_r_2;
2108 new_pos_a = atom[ia].position;
2109 new_pos_b = atom[ib].position + dr_b*v_ab;
2113 mass_a = atom[ia].mass;
2114 mass_b = atom[ib].mass;
2115 mass_sum = mass_a+mass_b;
2117 dot_v_r_1 = atom[ia].velocity.x*v_ab.
x
2118 + atom[ia].velocity.y*v_ab.
y + atom[ia].velocity.z*v_ab.
z;
2119 vb_1 = dot_v_r_1 * v_ab;
2120 vp_1 = atom[ia].velocity - vb_1;
2122 dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2123 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2124 vb_2 = dot_v_r_2 * v_ab;
2125 vp_2 = atom[ib].velocity - vb_2;
2127 vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
2134 if(dot_v_r_2 == dot_v_r_1) {
2138 delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);
2139 if(delta_T > maxtime ) delta_T = maxtime;
2143 v_Bond = sqrt(kbt/mass_b);
2146 dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
2147 dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
2149 dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
2150 dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
2152 new_pos_a = atom[ia].position + dr_a*v_ab;
2153 new_pos_b = atom[ib].position + dr_b*v_ab;
2160 vb_1 = dot_v_r_1 * v_ab;
2161 vb_2 = dot_v_r_2 * v_ab;
2163 new_vel_a = vp_1 + vb_1;
2164 new_vel_b = vp_2 + vb_2;
2169 atom[ia].position = new_pos_a;
2170 atom[ib].position = new_pos_b;
2172 else if ( virial == 0 ) {
2173 atom[ia].velocity = new_vel_a;
2174 atom[ib].velocity = new_vel_b;
2177 for ( j = 0; j < 2; j++ ) {
2180 new_pos = &new_pos_a;
2181 new_vel = &new_vel_a;
2185 new_pos = &new_pos_b;
2186 new_vel = &new_vel_b;
2188 Force df = (*new_vel - atom[Idx].velocity) *
2189 ( atom[Idx].mass * invdt );
2192 atom[Idx].velocity = *new_vel;
2193 atom[Idx].position = *new_pos;
2198 int partition = atom[Idx].partition;
2199 int slab = (int)floor((z-zmin)*idz);
2200 if (slab < 0) slab += nslabs;
2201 else if (slab >= nslabs) slab -= nslabs;
2204 ppreduction->
item(ppoffset ) += vir.
xx;
2205 ppreduction->
item(ppoffset+1) += vir.
yy;
2206 ppreduction->
item(ppoffset+2) += vir.
zz;
2225 if ( dt && virial ) *virial += wc;
2234 const int useSettle = simParams->
useSettle;
2242 int watmodel = simParams->
watmodel;
2243 if (watmodel ==
WAT_TIP4) wathgsize = 4;
2244 else if (watmodel ==
WAT_SWM4) wathgsize = 5;
2251 if ( ! settle_initialized ) {
2252 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2254 if (atom[ig].rigidBondLength > 0) {
2264 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2271 atom[ig].rigidBondLength,
2272 atom[oatm].rigidBondLength,
2273 settle_mOrmT, settle_mHrmT, settle_ra,
2274 settle_rb, settle_rc, settle_rra);
2275 settle_initialized = 1;
2290 int posRattleParam = 0;
2297 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2298 int hgs = atom[ig].hydrogenGroupSize;
2305 for (
int i = 0; i < hgs; ++i ) {
2306 ref[i] = atom[ig+i].position;
2307 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2308 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2315 BigReal tmp = atom[ig].rigidBondLength;
2317 if (hgs != wathgsize) {
2319 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms "
2320 "but the specified water model requires %d atoms.\n",
2321 atom[ig].
id+1, hgs, wathgsize);
2325 if (useSettle && !anyfixed) {
2330 if ( !(fixed[1] && fixed[2]) ) {
2331 dsq[icnt] = tmp * tmp;
2337 for (
int i = 1; i < hgs; ++i ) {
2338 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2339 if ( !(fixed[0] && fixed[i]) ) {
2340 dsq[icnt] = tmp * tmp;
2354 rattleListElem.
ig = ig;
2355 rattleListElem.
icnt = icnt;
2357 for (
int i = 0; i < icnt; ++i ) {
2361 rattleParamElem.
ia = a;
2362 rattleParamElem.
ib = b;
2363 rattleParamElem.
dsq = dsq[i];
2364 rattleParamElem.
rma = rmass[a];
2365 rattleParamElem.
rmb = rmass[b];
2373 for (
int ig = 0; ig <
numAtoms; ++ig ) {
2374 Force df = (
velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
2378 atom[ig].velocity =
velNew[ig];
2388 return rattle1old(timestep, virial, ppreduction);
2397 const int useSettle = simParams->
useSettle;
2399 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2402 int dieOnError = simParams->
rigidDie;
2410 for (
int j=0;j < n;j+=2) {
2413 for (
int i = 0; i < 3; ++i ) {
2414 ref[i] = atom[ig+i].position;
2415 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2418 for (
int i = 0; i < 3; ++i ) {
2419 ref[i+3] = atom[ig+i].position;
2420 pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
2423 settle_mOrmT, settle_mHrmT, settle_ra,
2424 settle_rb, settle_rc, settle_rra);
2427 for (
int i = 0; i < 3; ++i ) {
2428 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2432 for (
int i = 0; i < 3; ++i ) {
2433 velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
2441 for (
int i = 0; i < 3; ++i ) {
2442 ref[i] = atom[ig+i].position;
2443 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2446 settle_mOrmT, settle_mHrmT, settle_ra,
2447 settle_rb, settle_rc, settle_rra);
2448 for (
int i = 0; i < 3; ++i ) {
2449 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2467 int hgs = atom[ig].hydrogenGroupSize;
2468 for (
int i = 0; i < hgs; ++i ) {
2469 ref[i] = atom[ig+i].position;
2470 pos[i] = atom[ig+i].position;
2471 if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
2472 pos[i] += atom[ig+i].velocity * dt;
2501 for (
int i = 0; i < hgs; ++i ) {
2507 for (
int i = 0; i < hgs; ++i ) {
2508 velNew[ig+i] = (pos[i] - ref[i])*invdt;
2512 if ( consFailure ) {
2514 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom "
2515 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2518 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom "
2519 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2521 }
else if ( ! done ) {
2523 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom "
2524 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2527 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom "
2528 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2536 int hgs = atom[ig].hydrogenGroupSize;
2537 for (
int i = 0; i < hgs; ++i ) {
2538 velNew[ig+i] = atom[ig+i].velocity;
2539 posNew[ig+i] = atom[ig+i].position;
2544 for (
int ig = 0; ig <
numAtoms; ++ig )
2545 atom[ig].position =
posNew[ig];
2546 }
else if ( virial == 0 ) {
2547 for (
int ig = 0; ig <
numAtoms; ++ig )
2548 atom[ig].velocity =
velNew[ig];
2565 const int useSettle = simParams->
useSettle;
2567 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
2570 int dieOnError = simParams->
rigidDie;
2573 int ial[10], ibl[10];
2590 if ( ! settle_initialized ) {
2591 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2593 if (atom[ig].rigidBondLength > 0) {
2603 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2610 atom[ig].rigidBondLength,
2611 atom[oatm].rigidBondLength,
2612 settle_mOrmT, settle_mHrmT, settle_ra,
2613 settle_rb, settle_rc, settle_rra);
2614 settle_initialized = 1;
2628 int watmodel = simParams->
watmodel;
2629 if (watmodel ==
WAT_TIP4) wathgsize = 4;
2630 else if (watmodel ==
WAT_SWM4) wathgsize = 5;
2632 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2633 int hgs = atom[ig].hydrogenGroupSize;
2634 if ( hgs == 1 )
continue;
2637 for ( i = 0; i < hgs; ++i ) {
2638 ref[i] = atom[ig+i].position;
2639 pos[i] = atom[ig+i].position;
2640 vel[i] = atom[ig+i].velocity;
2641 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2643 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2646 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
else pos[i] += vel[i] * dt;
2649 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
2650 if (hgs != wathgsize) {
2652 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms "
2653 "but the specified water model requires %d atoms.\n",
2654 atom[ig].
id+1, hgs, wathgsize);
2658 if (useSettle && !anyfixed) {
2669 settle1(ref+2, pos+2, vel+2, invdt,
2670 settle_mOrmT, settle_mHrmT, settle_ra,
2671 settle_rb, settle_rc, settle_rra);
2679 swm4_omrepos(ref, pos, vel, invdt);
2683 settle_mOrmT, settle_mHrmT, settle_ra,
2684 settle_rb, settle_rc, settle_rra);
2686 tip4_omrepos(ref, pos, vel, invdt);
2693 if ( invdt == 0 )
for ( i = 0; i < wathgsize; ++i ) {
2694 atom[ig+i].position = pos[i];
2695 }
else if ( virial == 0 )
for ( i = 0; i < wathgsize; ++i ) {
2696 atom[ig+i].velocity = vel[i];
2697 }
else for ( i = 0; i < wathgsize; ++i ) {
2698 Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
2702 atom[ig+i].velocity = vel[i];
2708 partition = atom[ig].partition;
2709 int slab = (int)floor((z-zmin)*idz);
2710 if (slab < 0) slab += nslabs;
2711 else if (slab >= nslabs) slab -= nslabs;
2714 ppreduction->
item(ppoffset ) += vir.
xx;
2715 ppreduction->
item(ppoffset+1) += vir.
yy;
2716 ppreduction->
item(ppoffset+2) += vir.
zz;
2721 if ( !(fixed[1] && fixed[2]) ) {
2722 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2725 for ( i = 1; i < hgs; ++i ) {
2726 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2727 if ( !(fixed[0] && fixed[i]) ) {
2728 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
2732 if ( icnt == 0 )
continue;
2733 for ( i = 0; i < icnt; ++i ) {
2734 refab[i] = ref[ial[i]] - ref[ibl[i]];
2736 for ( i = 0; i < hgs; ++i ) {
2741 for ( iter = 0; iter < maxiter; ++iter ) {
2745 for ( i = 0; i < icnt; ++i ) {
2746 int a = ial[i];
int b = ibl[i];
2747 Vector pab = pos[a] - pos[b];
2750 BigReal diffsq = rabsq - pabsq;
2751 if ( fabs(diffsq) > (rabsq * tol2) ) {
2754 if ( rpab < ( rabsq * 1.0e-6 ) ) {
2761 BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
2765 if ( invdt != 0. ) {
2776 if ( consFailure ) {
2778 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom "
2779 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2782 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom "
2783 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2785 }
else if ( ! done ) {
2787 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom "
2788 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2791 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom "
2792 << (atom[ig].id + 1) <<
"!\n" <<
endi;
2798 if ( invdt == 0 )
for ( i = 0; i < hgs; ++i ) {
2799 atom[ig+i].position = pos[i];
2800 }
else if ( virial == 0 )
for ( i = 0; i < hgs; ++i ) {
2801 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2802 }
else for ( i = 0; i < hgs; ++i ) {
2803 Force df = netdp[i] * invdt;
2807 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2811 int partition = atom[ig].partition;
2812 int slab = (int)floor((z-zmin)*idz);
2813 if (slab < 0) slab += nslabs;
2814 else if (slab >= nslabs) slab -= nslabs;
2817 ppreduction->
item(ppoffset ) += vir.
xx;
2818 ppreduction->
item(ppoffset+1) += vir.
yy;
2819 ppreduction->
item(ppoffset+2) += vir.
zz;
2823 if ( dt && virial ) *virial += wc;
2834 const int useSettle = simParams->
useSettle;
2839 int dieOnError = simParams->
rigidDie;
2842 int ial[10], ibl[10];
2856 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2858 int hgs = atom[ig].hydrogenGroupSize;
2859 if ( hgs == 1 )
continue;
2862 for ( i = 0; i < hgs; ++i ) {
2863 ref[i] = atom[ig+i].position;
2864 vel[i] = atom[ig+i].velocity;
2865 rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
2866 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2867 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2870 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
2871 if (hgs != wathgsize) {
2872 NAMD_bug(
"Hydrogen group error caught in rattle2().");
2875 if (useSettle && !anyfixed) {
2884 settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
2891 settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
2896 for (i=0; i<hgs; i++) {
2897 atom[ig+i].velocity = vel[i];
2901 if ( !(fixed[1] && fixed[2]) ) {
2902 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2903 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2907 for ( i = 1; i < hgs; ++i ) {
2908 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2909 if ( !(fixed[0] && fixed[i]) ) {
2910 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2911 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
2912 ibl[icnt] = i; ++icnt;
2916 if ( icnt == 0 )
continue;
2918 for ( i = 0; i < icnt; ++i ) {
2919 refab[i] = ref[ial[i]] - ref[ibl[i]];
2923 for ( iter = 0; iter < maxiter; ++iter ) {
2925 for ( i = 0; i < icnt; ++i ) {
2926 int a = ial[i];
int b = ibl[i];
2927 Vector vab = vel[a] - vel[b];
2931 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2932 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2933 wc +=
outer(dp,rab);
2934 vel[a] += rmass[a] * dp;
2935 vel[b] -= rmass[b] * dp;
2944 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
2947 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
2951 for ( i = 0; i < hgs; ++i ) {
2952 atom[ig+i].velocity = vel[i];
2957 *virial += wc / ( 0.5 * dt );
2969 const int useSettle = simParams->
useSettle;
2974 int dieOnError = simParams->
rigidDie;
2977 int ial[10], ibl[10];
2991 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2993 int hgs = atom[ig].hydrogenGroupSize;
2994 if ( hgs == 1 )
continue;
2997 for ( i = 0; i < hgs; ++i ) {
2998 ref[i] = atom[ig+i].position;
2999 vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
3001 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
3002 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
3005 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
3006 if (hgs != wathgsize) {
3007 NAMD_bug(
"Hydrogen group error caught in rattle2().");
3010 if (useSettle && !anyfixed) {
3019 settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
3026 settle2(1.0, 1.0, ref, vel, dt, virial);
3031 for (i=0; i<hgs; i++) {
3032 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3036 if ( !(fixed[1] && fixed[2]) ) {
3037 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
3038 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3042 for ( i = 1; i < hgs; ++i ) {
3043 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3044 if ( !(fixed[0] && fixed[i]) ) {
3045 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
3046 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
3047 ibl[icnt] = i; ++icnt;
3051 if ( icnt == 0 )
continue;
3053 for ( i = 0; i < icnt; ++i ) {
3054 refab[i] = ref[ial[i]] - ref[ibl[i]];
3058 for ( iter = 0; iter < maxiter; ++iter ) {
3060 for ( i = 0; i < icnt; ++i ) {
3061 int a = ial[i];
int b = ibl[i];
3062 Vector vab = vel[a] - vel[b];
3066 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
3067 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
3068 wc +=
outer(dp,rab);
3069 vel[a] += rmass[a] * dp;
3070 vel[b] -= rmass[b] * dp;
3079 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
3082 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
3086 for ( i = 0; i < hgs; ++i ) {
3087 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3092 *virial += wc / ( 0.5 * dt );
3100 DebugM(2,
"loweAndersenVelocities\n");
3104 for (
int i = 0; i <
numAtoms; ++i) {
3107 v[i].position = atom[i].velocity;
3108 v[i].charge = atom[i].mass;
3110 DebugM(2,
"loweAndersenVelocities\n");
3115 DebugM(2,
"loweAndersenFinish\n");
3126 for (
int i = 0; i <
numAtoms; i++) {
3138 for (
int i = 0; i <
numAtoms; i++) {
3141 intRad[2*i+0] = rad - offset;
3142 intRad[2*i+1] = screen*(rad - offset);
3159 for (
int i = 0; i <
numAtoms; i++) {
3161 rhoi = rhoi0+coulomb_radius_offset;
3187 BigReal fij, expkappa, Dij, dEdai, dedasum;
3188 BigReal rhoi, rhoi0, psii, nbetapsi;
3189 BigReal gammapsi2, tanhi, daidr;
3190 for (
int i = 0; i <
numAtoms; i++) {
3194 expkappa = exp(-kappa*fij);
3195 Dij = epsilon_p_i - expkappa*epsilon_s_i;
3197 dEdai = -0.5*
COULOMB*atom[i].charge*atom[i].charge
3198 *(kappa*epsilon_s_i*expkappa-Dij/fij)/
bornRad[i];
3203 rhoi = rhoi0+coulomb_radius_offset;
3205 nbetapsi = -beta*psii;
3206 gammapsi2 = gamma*psii*psii;
3207 tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
3209 * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
3220 if (proxy.
size() > 0) {
3221 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3222 for (
int i = 0; i < proxy.
size(); i++) {
3223 int node = proxy[i];
3232 cp[node].recvData(msg);
3240 if (proxy.
size() > 0) {
3241 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3244 for (
int i = 0; i < proxy.
size(); i++) {
3245 int node = proxy[i];
3255 cp[node].recvData(msg);
3264 for (
int i = 0; i < msg->
psiSumLen; ++i ) {
3271 if (phase1BoxClosedCalled ==
true &&
3272 numGBISP1Arrived==proxy.
size() ) {
3278 numGBISP1Arrived == proxy.
size() &&
3279 numGBISP2Arrived == proxy.
size() &&
3280 numGBISP3Arrived == proxy.
size()) {
3290 for (
int i = 0; i < msg->
dEdaSumLen; ++i ) {
3297 if (phase2BoxClosedCalled ==
true &&
3298 numGBISP2Arrived==proxy.
size() ) {
3304 numGBISP1Arrived == proxy.
size() &&
3305 numGBISP2Arrived == proxy.
size() &&
3306 numGBISP3Arrived == proxy.
size()) {
3335 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3336 int hgs = atom[ig].hydrogenGroupSize;
3337 if ( hgs == 1 )
continue;
3338 for ( i = 0; i < hgs; ++i ) {
3339 ref[i] = atom[ig+i].position;
3340 rmass[i] = 1. / atom[ig+i].mass;
3341 fixed[i] = ( simParams->
fixedAtomsOn && atom[ig+i].atomFixed );
3342 if ( fixed[i] ) rmass[i] = 0.;
3347 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
3349 NAMD_die(
"Hydrogen group error caught in mollyAverage(). It's a bug!\n");
3351 if ( !(fixed[1] && fixed[2]) ) {
3352 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3355 for ( i = 1; i < hgs; ++i ) {
3356 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3357 if ( !(fixed[0] && fixed[i]) ) {
3358 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3362 if ( icnt == 0 )
continue;
3364 molly_lambda.
resize(numLambdas);
3365 lambda = &(molly_lambda[numLambdas - icnt]);
3366 for ( i = 0; i < icnt; ++i ) {
3367 refab[i] = ref[ial[i]] - ref[ibl[i]];
3370 iter=
average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
3371 if ( iter == maxiter ) {
3372 iout <<
iWARN <<
"Exceeded maximum number of iterations in mollyAverage().\n"<<
endi;
3403 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3404 int hgs = atom[ig].hydrogenGroupSize;
3405 if (hgs == 1 )
continue;
3406 for ( i = 0; i < hgs; ++i ) {
3407 ref[i] = atom[ig+i].position;
3409 rmass[i] = 1. / atom[ig+i].mass;
3410 fixed[i] = ( simParams->
fixedAtomsOn && atom[ig+i].atomFixed );
3411 if ( fixed[i] ) rmass[i] = 0.;
3415 if ( atom[ig].rigidBondLength > 0 ) {
3417 NAMD_die(
"Hydrogen group error caught in mollyMollify(). It's a bug!\n");
3419 if ( !(fixed[1] && fixed[2]) ) {
3420 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3423 for ( i = 1; i < hgs; ++i ) {
3424 if ( atom[ig+i].rigidBondLength > 0 ) {
3425 if ( !(fixed[0] && fixed[i]) ) {
3426 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3431 if ( icnt == 0 )
continue;
3432 lambda = &(molly_lambda[numLambdas]);
3434 for ( i = 0; i < icnt; ++i ) {
3435 refab[i] = ref[ial[i]] - ref[ibl[i]];
3438 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
3440 for ( i = 0; i < hgs; ++i ) {
3451 checkpoint_atom.
copy(atom);
3455 #if NAMD_SeparateWaters != 0
3456 checkpoint_numWaterAtoms = numWaterAtoms;
3463 atom.
copy(checkpoint_atom);
3464 numAtoms = atom.
size();
3467 doAtomUpdate =
true;
3473 #if NAMD_SeparateWaters != 0
3474 numWaterAtoms = checkpoint_numWaterAtoms;
3489 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
3499 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
3511 msg->
replica = CmiMyPartition();
3521 NAMD_bug(
"HomePatch::recvCheckpointLoad called during checkpointFree.");
3523 if ( msg->
replica != remote ) {
3524 NAMD_bug(
"HomePatch::recvCheckpointLoad message from wrong replica.");
3528 strcpy(newmsg->
key,key);
3532 newmsg->
pe = CkMyPe();
3533 newmsg->
replica = CmiMyPartition();
3545 doAtomUpdate =
true;
3570 CkpvAccess(_qd)->process();
3602 CkpvAccess(_qd)->process();
3614 CkpvAccess(_qd)->process();
3615 doAtomUpdate =
true;
3631 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3642 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3654 doPairlistCheck_newTolerance *= (1. - simParams->
pairlistShrink);
3655 doPairlistCheck_lattice =
lattice;
3656 doPairlistCheck_positions.
resize(numAtoms);
3658 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
3659 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3665 Lattice &lattice_old = doPairlistCheck_lattice;
3668 Vector center_delta = center_cur - center_old;
3672 for ( i=0; i<2; ++i ) {
3673 for (
int j=0; j<2; ++j ) {
3674 for (
int k=0; k<2; ++k ) {
3677 k ? min.
z : max.
z );
3680 corner_delta -= center_delta;
3682 if ( cd > max_cd ) max_cd = cd;
3686 max_cd = sqrt(max_cd);
3691 for ( i=0; i<n; ++i ) {
3693 p_delta -= center_delta;
3695 if ( pd > max_pd ) max_pd = pd;
3697 max_pd = sqrt(max_pd);
3699 BigReal max_tol = max_pd + max_cd;
3706 doPairlistCheck_newTolerance ) ) {
3707 doPairlistCheck_newTolerance *= (1. + simParams->
pairlistGrow);
3710 if ( max_tol > doPairlistCheck_newTolerance ) {
3711 doPairlistCheck_newTolerance = max_tol / (1. - simParams->
pairlistTrigger);
3713 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3729 while ( p_i != p_e ) {
3730 const int hgs = p_i->hydrogenGroupSize;
3733 if ( ngs > 5 ) ngs = 5;
3738 for ( i = 1; i < ngs; ++i ) {
3739 p_i[i].nonbondedGroupSize = 0;
3740 BigReal dx = p_i[i].position.x -
x;
3741 BigReal dy = p_i[i].position.y -
y;
3742 BigReal dz = p_i[i].position.z -
z;
3743 BigReal r2 = dx * dx + dy * dy + dz * dz;
3744 if ( r2 > hgcut )
break;
3745 else if ( r2 > maxrad2 ) maxrad2 = r2;
3747 p_i->nonbondedGroupSize = i;
3748 for ( ; i < hgs; ++i ) {
3749 p_i[i].nonbondedGroupSize = 1;
3755 NAMD_bug(
"hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
3772 if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
3773 ( bAwayDist*sysdimb < minSize*0.9999 ) ||
3774 ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
3776 NAMD_die(
"Periodic cell has become too small for original patch grid!\n"
3777 "Possible solutions are to restart from a recent checkpoint,\n"
3778 "increase margin, or disable useFlexibleCell for liquid simulation.");
3783 BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
3784 BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
3785 BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
3787 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
3788 NAMD_die(
"Periodic cell has become too small for original patch grid!\n"
3789 "There are probably many margin violations already reported.\n"
3790 "Possible solutions are to restart from a recent checkpoint,\n"
3791 "increase margin, or disable useFlexibleCell for liquid simulation.");
3801 int xdev, ydev, zdev;
3802 int problemCount = 0;
3806 for ( ; p_i != p_e; ++p_i ) {
3811 if (s.
x < minx) xdev = 0;
3812 else if (maxx <= s.
x) xdev = 2;
3815 if (s.
y < miny) ydev = 0;
3816 else if (maxy <= s.
y) ydev = 2;
3819 if (s.
z < minz) zdev = 0;
3820 else if (maxz <= s.
z) zdev = 2;
3823 if (mInfo[xdev][ydev][zdev]) {
3842 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3845 for (i=0; i<numNeighbors; i++) {
3861 int xdev, ydev, zdev;
3868 #if NAMD_SeparateWaters != 0
3870 int numLostWaterAtoms = 0;
3873 while ( atom_i != atom_e ) {
3880 if ( atom_i->migrationGroupSize ) {
3882 if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
3887 int mgs = atom_i->migrationGroupSize;
3889 for (
int j=atom_i->hydrogenGroupSize; j<mgs;
3890 j+=(atom_i+j)->hydrogenGroupSize ) {
3891 pos += (atom_i+j)->position;
3906 if (s.
x < minx) xdev = 0;
3907 else if (maxx <= s.
x) xdev = 2;
3910 if (s.
y < miny) ydev = 0;
3911 else if (maxy <= s.
y) ydev = 2;
3914 if (s.
z < minz) zdev = 0;
3915 else if (maxz <= s.
z) zdev = 2;
3920 if (mInfo[xdev][ydev][zdev]) {
3925 DebugM(3,
"Migrating atom " << atom_i->id <<
" from patch "
3926 <<
patchID <<
" with position " << atom_i->position <<
"\n");
3933 #if NAMD_SeparateWaters != 0
3942 if (atomIndex < numWaterAtoms)
3943 numLostWaterAtoms++;
3951 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
3959 #if NAMD_SeparateWaters != 0
3960 numWaterAtoms -= numLostWaterAtoms;
3964 int delpos = numAtoms - delnum;
3965 DebugM(4,
"numAtoms " << numAtoms <<
" deleted " << delnum <<
"\n");
3966 atom.
del(delpos,delnum);
3968 numAtoms = atom.
size();
3977 DebugM(1,
"Draining migration buffer ("<<i<<
","<<numMlBuf<<
")\n");
3982 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3986 if (!allMigrationIn) {
3987 DebugM(3,
"All Migrations NOT in, we are suspending patch "<<
patchID<<
"\n");
3988 migrationSuspended =
true;
3990 migrationSuspended =
false;
3992 allMigrationIn =
false;
4010 #if NAMD_SeparateWaters != 0
4029 for (mi = migrationList.
begin(); mi != migrationList.
end(); mi++) {
4030 DebugM(1,
"Migrating atom " << mi->id <<
" to patch "
4031 <<
patchID <<
" with position " << mi->position <<
"\n");
4032 if ( mi->migrationGroupSize ) {
4033 if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
4035 int mgs = mi->migrationGroupSize;
4037 for (
int j=mi->hydrogenGroupSize; j<mgs;
4038 j+=(mi+j)->hydrogenGroupSize ) {
4039 pos += (mi+j)->position;
4044 mother_transform = mi->transform;
4048 mi->transform = mother_transform;
4050 mi->position =
lattice.
nearest(mi->position,center,&(mi->transform));
4051 mother_transform = mi->transform;
4056 mi->transform = mother_transform;
4063 #endif // if (NAMD_SeparateWaters != 0)
4066 numAtoms = atom.
size();
4069 DebugM(3,
"Counter on " <<
patchID <<
" = " << patchMigrationCounter <<
"\n");
4070 if (!--patchMigrationCounter) {
4071 DebugM(3,
"All Migrations are in for patch "<<
patchID<<
"\n");
4072 allMigrationIn =
true;
4073 patchMigrationCounter = numNeighbors;
4074 if (migrationSuspended) {
4076 migrationSuspended =
false;
4089 #if NAMD_SeparateWaters != 0
4094 void HomePatch::separateAtoms() {
4106 if (atom.
size() < 0)
return;
4109 tempAtom.resize(numAtoms);
4119 int nonWaterIndex = 0;
4120 while (i < numAtoms) {
4126 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4129 if (waterIndex != i) {
4130 atom[waterIndex ] = atom[i ];
4131 atom[waterIndex + 1] = atom[i + 1];
4132 atom[waterIndex + 2] = atom[i + 2];
4133 if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3];
4134 if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4];
4138 waterIndex += wathgsize;
4144 for (
int j = 0; j < hgs; j++)
4145 tempAtom[nonWaterIndex + j] = atom[i + j];
4147 nonWaterIndex += hgs;
4159 for (i = 0; i < nonWaterIndex; i++)
4160 atom[waterIndex + i] = tempAtom[i];
4163 numWaterAtoms = waterIndex;
4176 if (al.
size() <= 0)
return;
4178 const int orig_atomSize = atom.
size();
4179 const int orig_alSize = al.
size();
4182 atom.
resize(orig_atomSize + orig_alSize);
4185 #if 0 // version where non-waters are moved to scratch first
4207 const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
4208 tempAtom.resize(orig_atom_numNonWaters + al.
size());
4209 for (
int i = 0; i < orig_atom_numNonWaters; i++)
4210 tempAtom[i] = atom[numWaterAtoms + i];
4214 int atom_waterIndex = numWaterAtoms;
4215 int atom_nonWaterIndex = orig_atom_numNonWaters;
4217 while (i < orig_alSize) {
4222 NAMD_bug(
"HomePatch::mergeAtomList() not updated for migration groups!");
4226 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4231 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4232 Transform mother_transform = al[i].transform;
4237 al[i+1].transform = mother_transform;
4242 al[i+2].transform = mother_transform;
4245 atom[atom_waterIndex ] = al[i ];
4246 atom[atom_waterIndex + 1] = al[i + 1];
4247 atom[atom_waterIndex + 2] = al[i + 2];
4249 atom_waterIndex += 3;
4257 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4258 Transform mother_transform = al[i].transform;
4261 for (
int j = 1; j < hgs; j++) {
4264 al[i+j].transform = mother_transform;
4268 for (
int j = 0; j < hgs; j++)
4269 tempAtom[atom_nonWaterIndex + j] = al[i + j];
4271 atom_nonWaterIndex += hgs;
4278 for (
int i = 0; i < atom_nonWaterIndex; i++)
4279 atom[atom_waterIndex + i] = tempAtom[i];
4282 numWaterAtoms = atom_waterIndex;
4283 numAtoms = atom.
size();
4306 int al_numWaterAtoms = 0;
4308 while (i < orig_alSize) {
4314 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4315 al_numWaterAtoms += wathgsize;
4323 if (al_numWaterAtoms > 0) {
4324 for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
4325 atom[i + al_numWaterAtoms] = atom[i];
4332 int atom_waterIndex = numWaterAtoms;
4333 int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
4335 while (i < orig_alSize) {
4340 NAMD_bug(
"HomePatch::mergeAtomList() not updated for migration groups!");
4344 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4349 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4350 Transform mother_transform = al[i].transform;
4355 al[i+1].transform = mother_transform;
4360 al[i+2].transform = mother_transform;
4363 atom[atom_waterIndex ] = al[i ];
4364 atom[atom_waterIndex + 1] = al[i + 1];
4365 atom[atom_waterIndex + 2] = al[i + 2];
4367 if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
4369 atom_waterIndex += wathgsize;
4377 al[i].position =
lattice.
nearest(al[i].position, center, &(al[i].transform));
4378 Transform mother_transform = al[i].transform;
4381 for (
int j = 1; j < hgs; j++) {
4384 al[i+j].transform = mother_transform;
4388 for (
int j = 0; j < hgs; j++)
4389 atom[atom_nonWaterIndex + j] = al[i + j];
4391 atom_nonWaterIndex += hgs;
4398 numWaterAtoms = atom_waterIndex;
4399 numAtoms = atom_nonWaterIndex;
4419 for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
4423 for (i=n-1;i>=0;i--) {
4425 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
4435 double big,dum,sum,temp;
4441 if ((temp=fabs(a[i][j])) > big) big=temp;
4442 if (big == 0.0)
NAMD_die(
"Singular matrix in routine ludcmp\n");
4448 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
4455 sum -= a[i][k]*a[k][j];
4457 if ( (dum=vv[i]*fabs(sum)) >= big) {
4472 if (a[j][j] == 0.0) a[j][j]=
TINY;
4475 for (i=j+1;i<n;i++) a[i][j] *= dum;
4486 gqij[i][ial[i]]=2.0*refab[i];
4487 gqij[i][ibl[i]]=-gqij[i][ial[i]];
4493 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) {
4531 G_q(refab,grhs,n,m,ial,ibl);
4532 for (k=1;k<=ntrial;k++) {
4541 multiplier = lambda[i];
4544 auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
4550 tmp[j]+=auxrhs[i][j];
4560 for ( i = 0; i < m; i++ ) {
4571 G_q(avgab,glhs,n,m,ial,ibl);
4575 iout<<
iINFO << glhs[i*n+0] <<
" " << glhs[i*n+1] <<
" " << glhs[i*n+2] << std::endl<<
endi;
4580 for (j=0; j<n; j++) {
4581 for (i=0; i<m;i++) {
4582 grhs2[i][j] = grhs[i][j]*imass[j];
4591 for (k1=0;k1<n;k1++) {
4592 fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
4620 gij[i]=avgab[i]*avgab[i]-length2[i];
4625 for(i=0;i<m-1;i++) {
4638 for (i=0;i<m;i++) errf += fabs(fvec[i]);
4645 for (i=0;i<m;i++) p[i] = -fvec[i];
4659 <<
" " << lambda[1] <<
" " << lambda[2] << std::endl<<
endi;
4662 if (errx <= tolx)
break;
4669 iout<<
iINFO <<
"LAMBDA:" << lambda[0] <<
" " << lambda[1] <<
" " << lambda[2] << std::endl<<
endi;
4679 Vector zero(0.0,0.0,0.0);
4690 tmpforce[i]=imass[i]*force[i];
4696 for ( i = 0; i < m; i++ ) {
4700 G_q(avgab,glhs,n,m,ial,ibl);
4701 G_q(refab,grhs,n,m,ial,ibl);
4703 for (j=0; j<n; j++) {
4704 for (i=0; i<m;i++) {
4705 grhs2[i][j] = grhs[i][j]*imass[j];
4715 fjac[j][i] += glhs[i][k]*grhs2[j][k];
4726 aux[i]+=grhs[i][j]*tmpforce[j];
4736 y[j] += aux[i]*glhs[i][j];
4745 tmpforce2[i]=imass[i]*y[i];
4754 Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
4755 tmpforce[ial[j]] += tmpf;
4756 tmpforce[ibl[j]] -= tmpf;
4761 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)
std::ostream & endi(std::ostream &s)
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()
#define NAMD_EVENT_START(eon, id)
Lphost * get_lphost(int atomid) const
static int compDistance(const void *a, const void *b)
void unregisterProxy(UnregisterProxyMsg *)
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
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t float const PatchPairRecord *__restrict__ const int *__restrict__ atomIndex
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
#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)