30 #include "ProxyMgr.decl.h" 49 #include "ComputeQMMgr.decl.h" 56 #define MIN_DEBUG_LEVEL 2 71 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);
78 const BigReal tol2,
const int maxiter,
79 bool& done,
bool& consFailure);
80 #define MASS_EPSILON (1.0e-35) //a very small floating point number 84 #if NAMD_SeparateWaters != 0 95 #define IS_HYDROGEN_GROUP_WATER(hgs, mass) \ 96 ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0))) 100 #ifdef TIMER_COLLECTION 101 const char *TimerSet::tlabel[TimerSet::NUMTIMERS] = {
122 settle_initialized = 0;
125 rattleListValid =
false;
126 rattleListValid_SOA =
false;
131 numGBISP1Arrived = 0;
132 numGBISP2Arrived = 0;
133 numGBISP3Arrived = 0;
134 phase1BoxClosedCalled =
false;
135 phase2BoxClosedCalled =
false;
136 phase3BoxClosedCalled =
false;
144 center = 0.5*(min+max);
149 aAwayDist = (max.x - min.x) * aAway;
156 bAwayDist = (max.y - min.y) * bAway;
163 cAwayDist = (max.z - min.z) * cAway;
168 migrationSuspended =
false;
169 allMigrationIn =
false;
170 marginViolations = 0;
176 flags.maxForceUsed = -1;
178 numAtoms = atom.size();
179 replacementForces = 0;
182 doPairlistCheck_newTolerance =
188 for (
int i = 0; i < numAtoms; ++i ) {
189 numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
195 sizeCudaAtomList = 0;
198 #ifdef NODEAWARE_PROXY_SPANNINGTREE 207 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 215 sort_solvent_atoms();
219 buildRattleList_SOA();
220 rattleListValid_SOA =
true;
226 #if NAMD_SeparateWaters != 0 229 tempAtom.resize(numAtoms);
239 gridForceIdxChecked=
false;
244 void HomePatch::write_tip4_props() {
245 printf(
"Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
248 void HomePatch::init_tip4() {
256 void ::HomePatch::init_swm4() {
278 sort_solvent_atoms();
289 #if NAMD_SeparateWaters != 0 302 { sequencer=sequencerPtr; }
306 { sequencer->
run(); }
308 void HomePatch::readPatchMap() {
313 patchMigrationCounter = numNeighbors
315 DebugM( 1,
"NumNeighbors for pid " <<
patchID<<
" is "<< numNeighbors <<
"\n");
317 for (n=0; n<numNeighbors; n++) {
318 realInfo[n].
destNodeID =
p->node(realInfo[n].destPatchID = nnPatchID[n]);
319 DebugM( 1,
" nnPatchID=" <<nnPatchID[n]<<
" nnNodeID="<< realInfo[n].destNodeID<<
"\n");
324 for (
int i=0; i<3; i++)
325 for (
int j=0; j<3; j++)
326 for (
int k=0; k<3; k++)
328 int pid =
p->pid(
p->index_a(
patchID)+i-1,
331 DebugM(5,
"ERROR, for patchID " <<
patchID <<
" I got neigh pid = " << pid <<
"\n");
334 ( (i-1) &&
p->periodic_a() ) ||
335 ( (j-1) &&
p->periodic_b() ) ||
336 ( (k-1) &&
p->periodic_c() ) )) {
337 mInfo[i][j][k] = NULL;
342 for (n = 0; n<numNeighbors; n++) {
343 if (pid == realInfo[n].destPatchID) {
344 mInfo[i][j][k] = &realInfo[n];
348 if (n == numNeighbors) {
349 DebugM(4,
"BAD News, I could not find PID " << pid <<
"\n");
354 DebugM(1,
"Patch("<<
patchID<<
") # of neighbors = " << numNeighbors <<
"\n");
360 #ifdef NODEAWARE_PROXY_SPANNINGTREE 362 #ifdef USE_NODEPATCHMGR 363 delete [] nodeChildren;
373 phase1BoxClosedCalled =
true;
382 }
else if (box == 6) {
384 }
else if (box == 7) {
386 }
else if (box == 8) {
387 phase2BoxClosedCalled =
true;
398 }
else if (box == 9) {
400 }
else if (box == 10) {
408 if ( replacementForces ) {
409 for (
int i = 0; i <
numAtoms; ++i ) {
410 if ( replacementForces[i].replace ) {
415 replacementForces = 0;
417 DebugM(1,
patchID <<
": " << CthSelf() <<
" awakening sequencer " 418 << sequencer->thread <<
"(" <<
patchID <<
") @" << CmiTimer() <<
"\n");
421 phase3BoxClosedCalled =
true;
426 if (numGBISP1Arrived == proxy.
size() &&
427 numGBISP2Arrived == proxy.
size() &&
428 numGBISP3Arrived == proxy.
size()) {
449 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 458 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 463 for ( ; *pe != n; ++pe );
469 #if USE_TOPOMAP && USE_SPANNING_TREE 471 int HomePatch::findSubroots(
int dim,
int* subroots,
int psize,
int* pidscopy){
474 int conesizes[6] = {0,0,0,0,0,0};
475 int conecounters[6] = {0,0,0,0,0,0};
476 int childcounter = 0;
479 for(
int i=0;i<psize;i++){
480 int cone = tmgr.getConeNumberForRank(pidscopy[i]);
481 cones[cone][conesizes[cone]++] = pidscopy[i];
484 while(childcounter<nChild){
485 for(
int i=0;i<6;i++){
486 if(conecounters[i]<conesizes[i]){
487 subroots[childcounter++] = cones[i][conecounters[i]++];
495 #endif // USE_TOPOMAP 499 int d1 = abs(*(
int *)a - CkMyPe());
500 int d2 = abs(*(
int *)b - CkMyPe());
512 CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
520 #ifdef NODEAWARE_PROXY_SPANNINGTREE 521 void HomePatch::buildNodeAwareSpanningTree(
void){
522 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 523 DebugFileTrace *dft = DebugFileTrace::Object();
525 dft->writeTrace(
"HomePatch[%d] has %d proxy on proc[%d] node[%d]\n",
patchID, proxy.
size(), CkMyPe(), CkMyNode());
526 dft->writeTrace(
"Proxies are: ");
527 for(
int i=0; i<proxy.
size(); i++) dft->writeTrace(
"%d(%d), ", proxy[i], CkNodeOf(proxy[i]));
528 dft->writeTrace(
"\n");
539 ProxyMgr::buildSinglePatchNodeAwareSpanningTree(
patchID, proxy, ptnTree);
543 setupChildrenFromProxySpanningTree();
548 void HomePatch::setupChildrenFromProxySpanningTree(){
549 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 552 if(ptnTree.size()==0) {
556 #ifdef USE_NODEPATCHMGR 558 delete [] nodeChildren;
564 CmiAssert(rootnode->
peIDs[0] == CkMyPe());
568 int internalChild = rootnode->
numPes-1;
569 int externalChild = ptnTree.size()-1;
573 if(internalSlots==0) {
577 internalChild = (internalSlots>internalChild)?internalChild:internalSlots;
581 nChild = externalChild+internalChild;
586 child =
new int[nChild];
588 for(
int i=0; i<externalChild; i++) {
589 child[i] = ptnTree.item(i+1).peIDs[0];
591 for(
int i=externalChild, j=1; i<nChild; i++, j++) {
592 child[i] = rootnode->
peIDs[j];
595 #ifdef USE_NODEPATCHMGR 598 CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
608 numNodeChild = externalChild;
609 if(internalChild) numNodeChild++;
610 delete [] nodeChildren;
611 nodeChildren =
new int[numNodeChild];
612 for(
int i=0; i<externalChild; i++) {
613 nodeChildren[i] = ptnTree.item(i+1).nodeID;
618 nodeChildren[numNodeChild-1] = rootnode->
nodeID;
621 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 622 DebugFileTrace *dft = DebugFileTrace::Object();
624 dft->writeTrace(
"HomePatch[%d] has %d children: ",
patchID, nChild);
625 for(
int i=0; i<nChild; i++)
626 dft->writeTrace(
"%d ", child[i]);
627 dft->writeTrace(
"\n");
633 #ifdef NODEAWARE_PROXY_SPANNINGTREE 638 ptnTree.resize(treesize);
639 int *pAllPes = msg->
allPes;
640 for(
int i=0; i<treesize; i++) {
642 delete [] oneNode->
peIDs;
644 oneNode->
nodeID = CkNodeOf(*pAllPes);
646 for(
int j=0; j<oneNode->
numPes; j++) {
647 oneNode->
peIDs[j] = *pAllPes;
652 setupChildrenFromProxySpanningTree();
658 if(ptnTree.size()==0)
return;
662 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 666 CmiAssert(CkMyPe() == msg->
allPes[0]);
675 #ifndef NODEAWARE_PROXY_SPANNINGTREE 682 for (i=0; i<n; i++) {
687 if (tree.
size() <= i)
break;
688 child[i-1] = tree[i];
692 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 702 if (tree.
size() == 0)
return;
705 msg->
node = CkMyPe();
714 #ifndef NODEAWARE_PROXY_SPANNINGTREE 718 int psize = proxy.
size();
719 if (psize == 0)
return;
721 int oldsize = oldtree.
size();
732 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
734 int oldindex = oldtree.
find(*pli);
735 if (oldindex != -1 && oldindex < psize) {
736 tree[oldindex] = *pli;
740 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
742 if (tree.
find(*pli) != -1)
continue;
744 while (tree[e] != -1) { e--;
if (e==-1) e = psize; }
747 while (tree[s] != -1) { s++;
if (s==psize+1) s = 1; }
753 if (oldsize==0 && nNonPatch) {
761 #if USE_TOPOMAP && USE_SPANNING_TREE 764 int *treecopy =
new int [psize];
776 for(i=0;i<psize;i++){
777 treecopy[i] = tree[i+1];
781 tmgr.sortRanksByHops(treecopy,nNonPatch);
782 tmgr.sortRanksByHops(treecopy+nNonPatch,
786 nChild = findSubroots(
proxySpanDim,subroots,psize,treecopy);
789 for(
int i=1;i<psize+1;i++){
791 for(
int j=0;j<nChild;j++)
792 if(tree[i]==subroots[j]){
796 if(isSubroot)
continue;
799 tmgr.sortIndexByHops(tree[i], subroots,
803 subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
809 NAMD_bug(
"HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
816 if (tree.
size() <= i)
break;
817 child[i-1] = tree[i];
825 CkPrintf(
"[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(),
patchID, psize, nNonPatch);
826 for (
int i=0; i<psize+1; i++) {
827 CkPrintf(
"%d ", tree[i]);
844 char *iszeroPtr = msg->
isZero;
850 for(
int i=0; i<msg->
flLen[k]; i++, rfPtr++, iszeroPtr++) {
851 if((*iszeroPtr)!=1) {
872 for ( ; f_i != f_e; ++f_i, ++
f ) *
f += *f_i;
888 int nf = msg->
flLen[k];
890 #pragma disjoint (*f_i, *f) 892 for (
int count = 0; count < nf; count++) {
894 f[count].x += f_i->
x;
895 f[count].y += f_i->
y;
896 f[count].z += f_i->
z;
919 ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
920 CkpvAccess(BOCclass_group).computeQMMgr) ;
928 if ( subP != NULL ) {
941 for (
int qmIter=0; qmIter<numQMAtms; qmIter++) {
943 if (qmAtmIndx[qmIter] == subP->
newID) {
979 if (!patchMapRead) { readPatchMap(); }
980 if (numNeighbors && !
simParams->staticAtomAssignment) {
984 copy_updates_to_AOS();
997 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED) 1008 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1009 if ( doMigration ) {
1012 int * __restrict ao_SOA = patchDataSOA.
sortOrder;
1013 int * __restrict unsort_SOA = patchDataSOA.
unsortOrder;
1014 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0) 1020 int * __restrict atomFixed = patchDataSOA.
atomFixed;
1021 for (
int j=0; j<n; ++j ) {
1023 if ( atomFixed[j] ) ao_SOA[--k2] = j;
1024 else ao_SOA[k++] = j;
1029 for (
int j=0; j<n; ++j ) {
1039 for (
int i = 0; i < n; ++i) {
1041 printf(
"ao_SOA[%d] = %d\n", i, ao_SOA[i]);
1046 for (
int i = 0; i < n; ++i) {
1054 const BigReal ucenter_x = ucenter.
x;
1055 const BigReal ucenter_y = ucenter.
y;
1056 const BigReal ucenter_z = ucenter.
z;
1058 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0) 1060 n_16 = (n + 15) & (~15);
1061 cudaAtomList.
resize(n_16);
1063 mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
1064 mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
1065 mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
1066 mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
1067 #elif defined(NAMD_AVXTILES) 1068 int n_avx = (n + 15) & (~15);
1069 cudaAtomList.
resize(n_avx);
1071 tiles.realloc(n, ac);
1073 if(doMigration) cudaAtomList.
resize(n);
1076 double * __restrict pos_x = patchDataSOA.
pos_x;
1077 double * __restrict pos_y = patchDataSOA.
pos_y;
1078 double * __restrict pos_z = patchDataSOA.
pos_z;
1079 float * __restrict charge = patchDataSOA.
charge;
1080 int * __restrict ao_SOA = patchDataSOA.
sortOrder;
1081 #ifndef NODEGROUP_FORCE_REGISTER 1083 for (
int k=0; k<n; ++k ) {
1084 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0) 1086 atom_x[k] = pos_x[j] - ucenter_x;
1087 atom_y[k] = pos_y[j] - ucenter_y;
1088 atom_z[k] = pos_z[j] - ucenter_z;
1089 atom_q[k] = charge_scaling * charge[j];
1097 ac[k].
x = pos_x[j] - ucenter_x;
1098 ac[k].
y = pos_y[j] - ucenter_y;
1099 ac[k].
z = pos_z[j] - ucenter_z;
1104 ac[k].
q = charge_scaling * charge[j];
1108 if(!
simParams->CUDASOAintegrate || doMigration){
1109 for (
int k=0; k<n; ++k ) {
1110 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0) 1112 atom_x[k] = pos_x[j] - ucenter_x;
1113 atom_y[k] = pos_y[j] - ucenter_y;
1114 atom_z[k] = pos_z[j] - ucenter_z;
1115 atom_q[k] = charge_scaling * charge[j];
1122 ac[k].
x = pos_x[j] - ucenter_x;
1123 ac[k].
y = pos_y[j] - ucenter_y;
1124 ac[k].
z = pos_z[j] - ucenter_z;
1129 ac[k].
q = charge_scaling * charge[j];
1136 doMigration = doMigration && numNeighbors;
1138 doMigration = doMigration || ! patchMapRead;
1140 doMigration = doMigration || doAtomUpdate;
1141 doAtomUpdate =
false;
1152 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1164 const double * __restrict pos_x = patchDataSOA.
pos_x;
1165 const double * __restrict pos_y = patchDataSOA.
pos_y;
1166 const double * __restrict pos_z = patchDataSOA.
pos_z;
1167 const float * __restrict charge = patchDataSOA.
charge;
1168 const int * __restrict vdwType = patchDataSOA.
vdwType;
1170 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1176 for (
int i=0; i < n; i++) {
1177 p_i[i].position.x = pos_x[i];
1178 p_i[i].position.y = pos_y[i];
1179 p_i[i].position.z = pos_z[i];
1180 p_i[i].charge = charge[i];
1181 p_i[i].vdwType = vdwType[i];
1183 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1184 p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
1186 p_i[i].hydrogenGroupSize = hydrogenGroupSize[i];
1189 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1190 const int * __restrict sortOrder = patchDataSOA.
sortOrder;
1192 const int * __restrict
id = patchDataSOA.
id;
1193 #if defined(MEM_OPT_VERSION) 1194 const int * __restrict exclId = patchDataSOA.
exclId;
1195 const int * __restrict sigId = patchDataSOA.
sigId;
1197 const int * __restrict atomFixed = patchDataSOA.
atomFixed;
1198 const int * __restrict groupFixed = patchDataSOA.
groupFixed;
1201 #ifndef USE_NO_BITFIELDS 1203 #endif // USE_NO_BITFIELDS 1204 for (
int i=0; i < n; i++) {
1205 #ifndef USE_NO_BITFIELDS 1206 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1209 #if defined(MEM_OPT_VERSION) 1210 pExtCopy_i[i].
id =
id[i];
1211 pExtCopy_i[i].exclId = exclId[i];
1212 ASSERT(atomFixed[i] == 0 || atomFixed[i] == 1);
1213 ASSERT(groupFixed[i] == 0 || groupFixed[i] == 1);
1216 pExtCopy_i[i].sigId = (sigId[i] | atomFixedBit | groupFixedBit);
1218 ASSERT(atomFixed[i] == 0 || atomFixed[i] == 1);
1219 ASSERT(groupFixed[i] == 0 || groupFixed[i] == 1);
1222 pExtCopy_i[i].
id = (
id[i] | atomFixedBit | groupFixedBit);
1223 #endif // if defined(MEM_OPT_VERSION) 1225 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1226 pExt_i[i].sortOrder = sortOrder[i];
1228 pExt_i[i].id =
id[i];
1229 #if defined(MEM_OPT_VERSION) 1230 pExt_i[i].exclId = exclId[i];
1231 pExt_i[i].sigId = sigId[i];
1233 pExt_i[i].atomFixed = atomFixed[i];
1234 pExt_i[i].groupFixed = groupFixed[i];
1235 #endif // USE_NO_BITFIELDS 1244 const double * __restrict pos_x = patchDataSOA.
pos_x;
1245 const double * __restrict pos_y = patchDataSOA.
pos_y;
1246 const double * __restrict pos_z = patchDataSOA.
pos_z;
1247 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1251 #ifndef NODEGROUP_FORCE_REGISTER 1253 for (
int i=0; i < n; i++) {
1254 p_i[i].position.x = pos_x[i];
1255 p_i[i].position.y = pos_y[i];
1256 p_i[i].position.z = pos_z[i];
1257 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1258 p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
1262 if(!
simParams->CUDASOAintegrate || doMigration){
1263 for (
int i=0; i < n; i++) {
1264 p_i[i].position.x = pos_x[i];
1265 p_i[i].position.y = pos_y[i];
1266 p_i[i].position.z = pos_z[i];
1267 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP) 1268 p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
1273 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1274 const int * __restrict sortOrder = patchDataSOA.
sortOrder;
1275 #ifdef NODEGROUP_FORCE_REGISTER 1277 if(!
simParams->CUDASOAintegrate || doMigration){
1278 for (
int i=0; i < n; i++) {
1279 pExt_i[i].sortOrder = sortOrder[i];
1284 for (
int i=0; i < n; i++) {
1285 pExt_i[i].sortOrder = sortOrder[i];
1293 #ifdef NODEGROUP_FORCE_REGISTER 1309 numGBISP1Arrived = 0;
1310 phase1BoxClosedCalled =
false;
1311 numGBISP2Arrived = 0;
1312 phase2BoxClosedCalled =
false;
1313 numGBISP3Arrived = 0;
1314 phase3BoxClosedCalled =
false;
1315 if (doMigration || isNewProxyAdded)
1320 if (doMigration || isNewProxyAdded) {
1329 int pidsPreAllocated = 1;
1332 npid = proxy.
size();
1333 pids =
new int[npid];
1334 pidsPreAllocated = 0;
1336 int *pide = pids + proxy.
size();
1337 int patchNodesLast =
1339 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
1349 #ifdef NODEAWARE_PROXY_SPANNINGTREE 1350 #ifdef USE_NODEPATCHMGR 1351 npid = numNodeChild;
1352 pids = nodeChildren;
1359 pidsPreAllocated = 0;
1361 for (
int i=0; i<nChild; i++) pids[i] = child[i];
1368 int pdMsgPLLen =
p.
size();
1369 int pdMsgAvgPLLen = 0;
1379 pdMsgVLLen =
v.
size();
1386 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1392 int lcpoTypeLen = 0;
1394 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1399 int pdMsgPLExtLen = 0;
1400 if(doMigration || isNewProxyAdded) {
1404 int cudaAtomLen = 0;
1405 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 1410 #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0 1413 cudaAtomLen = (
numAtoms + 15) & (~15);
1416 ProxyDataMsg *nmsg =
new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1422 nmsg->
plLen = pdMsgPLLen;
1434 nmsg->
vlLen = pdMsgVLLen;
1443 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
1444 for (
int i = 0; i <
numAtoms * 2; i++) {
1451 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
1452 for (
int i = 0; i <
numAtoms; i++) {
1458 if(doMigration || isNewProxyAdded){
1463 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) 1469 #if NAMD_SeparateWaters != 0 1471 nmsg->numWaterAtoms = numWaterAtoms;
1474 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK) 1475 nmsg->isFromImmMsgCall = 0;
1478 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 1479 DebugFileTrace *dft = DebugFileTrace::Object();
1481 dft->writeTrace(
"HP::posReady: for HomePatch[%d], sending proxy msg to: ",
patchID);
1482 for(
int i=0; i<npid; i++) {
1483 dft->writeTrace(
"%d ", pids[i]);
1485 dft->writeTrace(
"\n");
1489 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 1490 if (isProxyChanged || localphs == NULL)
1495 for (
int i=0; i<nphs; i++) {
1496 CmiDestoryPersistent(localphs[i]);
1500 localphs =
new PersistentHandle[npid];
1502 for (
int i=0; i<npid; i++) {
1503 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) 1505 localphs[i] = CmiCreateNodePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1508 localphs[i] = CmiCreatePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1512 CmiAssert(nphs == npid && localphs != NULL);
1513 CmiUsePersistentHandle(localphs, nphs);
1515 if(doMigration || isNewProxyAdded) {
1520 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 1521 CmiUsePersistentHandle(NULL, 0);
1523 isNewProxyAdded = 0;
1526 if(!pidsPreAllocated)
delete [] pids;
1527 DebugM(4,
"patchID("<<
patchID<<
") doing positions Ready\n");
1529 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY 1530 positionPtrBegin =
p.
begin();
1531 positionPtrEnd =
p.
end();
1558 fprintf(stderr,
"Patch %d atom IDS\n", this->
patchID);
1560 for(
int i = 0 ; i <
numAtoms; i++){
1561 fprintf(stderr,
"atom[%d] = %d %d %d\n", i, atom[i].
id,
1577 #ifdef NODEGROUP_FORCE_REGISTER 1582 if (n > MigrationCUDAKernel::kMaxAtomsPerPatch) {
1583 NAMD_die(
"Device migration does not currently support patches with greater than 2048 atoms.\n" 1584 "Please run with a smaller margin or without device migration.\n");
1594 size_t nbytes = patchDataSOA.
numBytes;
1596 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 1606 void HomePatch::clearAtomMap() {
1621 const int numAtomAlloc = 2048;
1623 cudaAtomList.
reserve(numAtomAlloc);
1638 size_t nbytes = patchDataSOA.
numBytes;
1640 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 1650 doMigration = doMigration || ! patchMapRead;
1652 doMigration = doMigration || doAtomUpdate;
1653 doAtomUpdate =
false;
1657 const bool updateAtomMap =
simParams->updateAtomMap;
1661 if (startup || updateAtomMap) {
1664 #ifndef USE_NO_BITFIELDS 1671 #ifndef USE_NO_BITFIELDS 1672 ASSERT(atomFixed == 0 || atomFixed == 1);
1673 ASSERT(groupFixed == 0 || groupFixed == 1);
1676 pExtCopy_i[i].
id = (
id | atomFixedBit | groupFixedBit);
1679 pExt_i[i].atomFixed = atomFixed;
1680 pExt_i[i].groupFixed = groupFixed;
1681 #endif // USE_NO_BITFIELDS 1688 p_i[i].
vdwType = atom[i].vdwType;
1692 p_i[i].partition = atom[i].partition;
1697 patchDataSOA.
id[i] = int(atom[i].
id);
1698 patchDataSOA.
atomFixed[i] = int(atom[i].atomFixed);
1699 patchDataSOA.
groupFixed[i] = int(atom[i].groupFixed);
1701 }
else if (updateAtomMap) {
1703 copy_atoms_to_SOA();
1708 #ifdef NODEGROUP_FORCE_REGISTER 1719 int pidsPreAllocated = 1;
1722 npid = proxy.
size();
1723 pids =
new int[npid];
1724 pidsPreAllocated = 0;
1726 int *pide = pids + proxy.
size();
1727 int patchNodesLast =
1729 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
1739 #ifdef NODEAWARE_PROXY_SPANNINGTREE 1740 #ifdef USE_NODEPATCHMGR 1741 npid = numNodeChild;
1742 pids = nodeChildren;
1749 pidsPreAllocated = 0;
1751 for (
int i=0; i<nChild; i++) pids[i] = child[i];
1758 int pdMsgPLLen =
p.
size();
1759 int pdMsgAvgPLLen = 0;
1767 int lcpoTypeLen = 0;
1769 int pdMsgPLExtLen = 0;
1770 if(doMigration || isNewProxyAdded) {
1774 int cudaAtomLen = 0;
1780 #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0 1783 cudaAtomLen = (
numAtoms + 15) & (~15);
1786 ProxyDataMsg *nmsg =
new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1792 nmsg->
plLen = pdMsgPLLen;
1799 nmsg->
vlLen = pdMsgVLLen;
1803 if(doMigration || isNewProxyAdded){
1808 #if defined(NAMD_CUDA) || defined(NAMD_MIC) 1814 #if NAMD_SeparateWaters != 0 1816 nmsg->numWaterAtoms = numWaterAtoms;
1819 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK) 1820 nmsg->isFromImmMsgCall = 0;
1823 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 1824 DebugFileTrace *dft = DebugFileTrace::Object();
1826 dft->writeTrace(
"HP::posReady: for HomePatch[%d], sending proxy msg to: ",
patchID);
1827 for(
int i=0; i<npid; i++) {
1828 dft->writeTrace(
"%d ", pids[i]);
1830 dft->writeTrace(
"\n");
1834 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 1835 if (isProxyChanged || localphs == NULL)
1838 for (
int i=0; i<nphs; i++) {
1839 CmiDestoryPersistent(localphs[i]);
1843 localphs =
new PersistentHandle[npid];
1845 for (
int i=0; i<npid; i++) {
1846 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) 1848 localphs[i] = CmiCreateNodePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1851 localphs[i] = CmiCreatePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
1855 CmiAssert(nphs == npid && localphs != NULL);
1856 CmiUsePersistentHandle(localphs, nphs);
1858 if(doMigration || isNewProxyAdded) {
1863 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 1864 CmiUsePersistentHandle(NULL, 0);
1866 isNewProxyAdded = 0;
1869 if(!pidsPreAllocated)
delete [] pids;
1870 DebugM(4,
"patchID("<<
patchID<<
") doing positions Ready\n");
1872 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY 1873 positionPtrBegin =
p.
begin();
1874 positionPtrEnd =
p.
end();
1887 #endif // NODEGROUP_FORCE_REGISTER 1902 if (!patchMapRead) { readPatchMap(); }
1904 if (numNeighbors && !
simParams->staticAtomAssignment) {
1913 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED) 1922 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES) || defined(NAMD_HIP) 1923 #ifdef NAMD_AVXTILES 1926 if ( doMigration ) {
1929 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_AVXTILES) || \ 1930 (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0) 1931 int *ao =
new int[n];
1936 for (
int j=0; j<n; ++j ) {
1938 if ( a_i[j].atomFixed ) ao[--k2] = j;
1944 for (
int j=0; j<n; ++j ) {
1949 for (
int i=0; i<n; ++i ) {
1954 for (
int i = 0; i < n; ++i) {
1963 const BigReal ucenter_x = ucenter.
x;
1964 const BigReal ucenter_y = ucenter.
y;
1965 const BigReal ucenter_z = ucenter.
z;
1967 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0) 1969 n_16 = (n + 15) & (~15);
1970 cudaAtomList.
resize(n_16);
1972 mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
1973 mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
1974 mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
1975 mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
1976 #elif defined(NAMD_AVXTILES) 1977 int n_avx = (n + 15) & (~15);
1978 cudaAtomList.
resize(n_avx);
1980 tiles.realloc(n, ac);
1986 reallocate_host<CudaAtom>(&cudaAtomList, &sizeCudaAtomList, n);
1991 for (
int k=0; k<n; ++k ) {
1992 #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0) 1994 atom_x[k] = a[j].
position.
x - ucenter_x;
1995 atom_y[k] = a[j].
position.
y - ucenter_y;
1996 atom_z[k] = a[j].
position.
z - ucenter_z;
1997 atom_q[k] = charge_scaling * a[j].
charge;
2003 ac[k].q = charge_scaling * a[j].
charge;
2006 #ifdef NAMD_AVXTILES 2010 for (
int k=n; k<n_avx; ++k ) {
2019 #ifdef NAMD_AVXTILES 2021 }
else doMigration = doMigration && numNeighbors;
2024 doMigration = doMigration && numNeighbors;
2026 doMigration = doMigration || ! patchMapRead;
2028 doMigration = doMigration || doAtomUpdate;
2029 doAtomUpdate =
false;
2040 #if ! (defined(NAMD_CUDA) || defined(NAMD_HIP)) 2041 #if defined(NAMD_AVXTILES) 2057 for ( i=0; i<n; ++i ) {
2072 numGBISP1Arrived = 0;
2073 phase1BoxClosedCalled =
false;
2074 numGBISP2Arrived = 0;
2075 phase2BoxClosedCalled =
false;
2076 numGBISP3Arrived = 0;
2077 phase3BoxClosedCalled =
false;
2078 if (doMigration || isNewProxyAdded)
2083 if (doMigration || isNewProxyAdded) {
2091 int pidsPreAllocated = 1;
2094 npid = proxy.
size();
2095 pids =
new int[npid];
2096 pidsPreAllocated = 0;
2098 int *pide = pids + proxy.
size();
2099 int patchNodesLast =
2101 for ( pli = proxy.
begin(); pli != proxy.
end(); ++pli )
2111 #ifdef NODEAWARE_PROXY_SPANNINGTREE 2112 #ifdef USE_NODEPATCHMGR 2113 npid = numNodeChild;
2114 pids = nodeChildren;
2121 pidsPreAllocated = 0;
2123 for (
int i=0; i<nChild; i++) pids[i] = child[i];
2130 int pdMsgPLLen =
p.
size();
2131 int pdMsgAvgPLLen = 0;
2138 pdMsgVLLen =
v.
size();
2143 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
2148 int lcpoTypeLen = 0;
2149 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
2153 int pdMsgPLExtLen = 0;
2154 if(doMigration || isNewProxyAdded) {
2158 int cudaAtomLen = 0;
2159 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 2161 #elif defined(NAMD_AVXTILES) 2163 cudaAtomLen = (
numAtoms + 15) & (~15);
2164 #elif defined(NAMD_MIC) 2165 #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0 2168 cudaAtomLen = (
numAtoms + 15) & (~15);
2172 ProxyDataMsg *nmsg =
new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
2178 nmsg->
plLen = pdMsgPLLen;
2188 nmsg->
vlLen = pdMsgVLLen;
2194 if (
flags.
doGBIS && (doMigration || isNewProxyAdded)) {
2195 for (
int i = 0; i <
numAtoms * 2; i++) {
2200 if (
flags.
doLCPO && (doMigration || isNewProxyAdded)) {
2201 for (
int i = 0; i <
numAtoms; i++) {
2207 if(doMigration || isNewProxyAdded){
2212 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES) || defined(NAMD_HIP) 2213 #ifdef NAMD_AVXTILES 2223 #if NAMD_SeparateWaters != 0 2225 nmsg->numWaterAtoms = numWaterAtoms;
2228 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK) 2229 nmsg->isFromImmMsgCall = 0;
2232 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG) 2233 DebugFileTrace *dft = DebugFileTrace::Object();
2235 dft->writeTrace(
"HP::posReady: for HomePatch[%d], sending proxy msg to: ",
patchID);
2236 for(
int i=0; i<npid; i++) {
2237 dft->writeTrace(
"%d ", pids[i]);
2239 dft->writeTrace(
"\n");
2243 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 2244 if (isProxyChanged || localphs == NULL)
2249 for (
int i=0; i<nphs; i++) {
2250 CmiDestoryPersistent(localphs[i]);
2254 localphs =
new PersistentHandle[npid];
2256 for (
int i=0; i<npid; i++) {
2257 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) 2259 localphs[i] = CmiCreateNodePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
2262 localphs[i] = CmiCreatePersistent(pids[i], persist_size,
sizeof(envelope)+
sizeof(
ProxyDataMsg));
2266 CmiAssert(nphs == npid && localphs != NULL);
2267 CmiUsePersistentHandle(localphs, nphs);
2269 if(doMigration || isNewProxyAdded) {
2274 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE 2275 CmiUsePersistentHandle(NULL, 0);
2277 isNewProxyAdded = 0;
2280 if(!pidsPreAllocated)
delete [] pids;
2281 DebugM(4,
"patchID("<<
patchID<<
") doing positions Ready\n");
2283 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY 2284 positionPtrBegin =
p.
begin();
2285 positionPtrEnd =
p.
end();
2312 replacementForces =
f;
2318 for (
int i = 0; i <
numAtoms; ++i )
2320 f_saved[ftag][i] =
f[ftag][i];
2325 void HomePatch::sort_solvent_atoms() {
2326 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED) 2338 numSolventAtoms = 0;
2342 if (a[i].
isWater) numSolventAtoms++;
2345 numSoluteAtoms =
numAtoms - numSolventAtoms;
2356 atom.
swap(sortatom);
2359 numWaters = numSolventAtoms / 3;
2380 p->
maxAtoms = ((natoms + padding - 1) / padding) * padding;
2387 (numdoubles *
sizeof(double) + 9*
sizeof(
float) + 17*
sizeof(int));
2398 p->
buffer = (
unsigned char *) mybuffer;
2399 unsigned char *t = p->
buffer;
2400 p->
pos_x = (
double *) t;
2402 p->
pos_y = (
double *) t;
2404 p->
pos_z = (
double *) t;
2426 p->
sigId = (
int *) t;
2432 p->
vel_x = (
double *) t;
2434 p->
vel_y = (
double *) t;
2436 p->
vel_z = (
double *) t;
2446 p->
mass = (
float *) t;
2524 void HomePatch::copy_atoms_to_SOA() {
2528 size_t nbytes = patchDataSOA.
numBytes;
2530 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 2541 patchDataSOA.
pos_x[i] = atom[i].position.x;
2542 patchDataSOA.
pos_y[i] = atom[i].position.y;
2543 patchDataSOA.
pos_z[i] = atom[i].position.z;
2544 patchDataSOA.
charge[i] = atom[i].charge;
2545 patchDataSOA.
vdwType[i] = int(atom[i].vdwType);
2554 patchDataSOA.
id[i] = int(atom[i].
id);
2555 #ifdef MEM_OPT_VERSION 2556 patchDataSOA.
exclId[i] = int(atom[i].exclId);
2557 patchDataSOA.
sigId[i] = int(atom[i].sigId);
2559 patchDataSOA.
atomFixed[i] = int(atom[i].atomFixed);
2560 patchDataSOA.
groupFixed[i] = int(atom[i].groupFixed);
2562 patchDataSOA.
vel_x[i] = atom[i].velocity.x;
2563 patchDataSOA.
vel_y[i] = atom[i].velocity.y;
2564 patchDataSOA.
vel_z[i] = atom[i].velocity.z;
2568 patchDataSOA.
mass[i] = atom[i].mass;
2570 patchDataSOA.
status[i] = atom[i].status;
2571 patchDataSOA.
transform_i[i] = atom[i].transform.i;
2572 patchDataSOA.
transform_j[i] = atom[i].transform.j;
2573 patchDataSOA.
transform_k[i] = atom[i].transform.k;
2579 calculate_derived_SOA();
2583 void HomePatch::calculate_derived_SOA() {
2586 patchDataSOA.
recipMass[i] = (atom[i].mass > 0 ? 1.f / atom[i].mass : 0);
2604 void HomePatch::copy_forces_to_SOA() {
2646 for (
int i = 0; i <
numAtoms; i++) {
2656 for (
int i = 0; i <
numAtoms; i++) {
2664 for (
int i = 0; i <
numAtoms; i++) {
2675 void HomePatch::copy_updates_to_AOS() {
2677 atom[i].velocity.x = patchDataSOA.
vel_x[i];
2678 atom[i].velocity.y = patchDataSOA.
vel_y[i];
2679 atom[i].velocity.z = patchDataSOA.
vel_z[i];
2680 atom[i].position.x = patchDataSOA.
pos_x[i];
2681 atom[i].position.y = patchDataSOA.
pos_y[i];
2682 atom[i].position.z = patchDataSOA.
pos_z[i];
2687 void HomePatch::copy_forces_to_AOS() {
2692 current_force[i].x = patchDataSOA.
f_slow_x[i];
2693 current_force[i].y = patchDataSOA.
f_slow_y[i];
2694 current_force[i].z = patchDataSOA.
f_slow_z[i];
2700 current_force[i].x = patchDataSOA.
f_nbond_x[i];
2701 current_force[i].y = patchDataSOA.
f_nbond_y[i];
2702 current_force[i].z = patchDataSOA.
f_nbond_z[i];
2708 current_force[i].x = patchDataSOA.
f_normal_x[i];
2709 current_force[i].y = patchDataSOA.
f_normal_y[i];
2710 current_force[i].z = patchDataSOA.
f_normal_z[i];
2717 void HomePatch::zero_global_forces_SOA() {
2723 #undef DEBUG_REDISTRIB_FORCE 2724 #undef DEBUG_REDISTRIB_FORCE_VERBOSE 2741 void HomePatch::redistrib_colinear_lp_force(
2745 BigReal distance = distance_f;
2750 distance *= r_jk_rlength;
2751 BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
2752 const Vector fja = (1. + scale + distance)*fi - r_jk*fdot;
2753 const Vector fka = fi - fja;
2787 #define FIX_FOR_WATER 2788 void HomePatch::redistrib_relative_lp_force(
2791 Tensor *virial,
int midpt) {
2792 #ifdef DEBUG_REDISTRIB_FORCE 2794 foldnet = fi + fj + fk + fl;
2795 toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
2797 Vector fja(0), fka(0), fla(0);
2802 BigReal fdot = (fi*r) * invr2;
2809 s = rj - 0.5*(rk + rl);
2820 #if !defined(FIX_FOR_WATER) 2829 #if !defined(FIX_FOR_WATER) 2832 Vector ft = fi - fr - fp;
2842 fpdot = (fi*
p) / sqrt(p2);
2850 ft = cross(s,
v) * invs2;
2861 BigReal srdot = (s*r) * invs2;
2864 BigReal stdot = (s*t) * invs2;
2867 BigReal fact = rrdot*fpdot*invtt*invq;
2871 fja += fp*(1+srdot) + fq*stdot;
2873 ft = fq*(1+stdot) + fp*srdot;
2885 va +=
outer(fka,rk);
2886 va +=
outer(fla,rl);
2896 #ifdef DEBUG_REDISTRIB_FORCE 2897 #define TOL_REDISTRIB 1e-4 2899 fnewnet = fi + fj + fk + fl;
2900 tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
2901 Vector fdiff = fnewnet - foldnet;
2902 Vector tdiff = tnewnet - toldnet;
2903 if (fdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
2904 printf(
"Error: force redistribution for water exceeded tolerance: " 2905 "fdiff=(%f, %f, %f)\n", fdiff.
x, fdiff.
y, fdiff.
z);
2907 if (tdiff.
length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
2908 printf(
"Error: torque redistribution for water exceeded tolerance: " 2909 "tdiff=(%f, %f, %f)\n", tdiff.
x, tdiff.
y, tdiff.
z);
2914 void HomePatch::redistrib_ap_force(
Vector& fi,
Vector& fj) {
2919 fi = fi_old + fj_old;
2920 fj = fi_old + fj_old;
2932 void HomePatch::redistrib_lp_water_force(
2937 #ifdef DEBUG_REDISTRIB_FORCE 2941 Vector totforce, tottorque;
2942 totforce = f_ox + f_h1 + f_h2 + f_lp;
2943 tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
2946 tottorque += cross(f_lp, p_lp);
2952 Vector fad_ox(0), fad_h(0);
2955 Vector r_ox_lp = p_lp - p_ox;
2957 invlen2_r_ox_lp *= invlen2_r_ox_lp;
2958 BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
2959 Vector f_rad = r_ox_lp * rad_factor;
2964 Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
2978 Vector f_ang = f_lp - f_rad;
2982 BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
2983 BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
2984 BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
2986 fad_ox += (f_ang * oxcomp);
2987 fad_h += (f_ang * hydcomp);
2992 vir +=
outer(fad_h, p_h1);
2993 vir +=
outer(fad_h, p_h2);
2994 vir -=
outer(f_lp, p_lp);
3004 #ifdef DEBUG_REDISTRIB_FORCE 3006 Vector newforce, newtorque;
3007 newforce = f_ox + f_h1 + f_h2;
3008 newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
3009 Vector fdiff = newforce - totforce;
3010 Vector tdiff = newtorque - tottorque;
3012 if (error > 0.0001) {
3013 printf(
"Error: Force redistribution for water " 3014 "exceeded force tolerance: error=%f\n", error);
3016 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE 3017 printf(
"Error in net force: %f\n", error);
3021 if (error > 0.0001) {
3022 printf(
"Error: Force redistribution for water " 3023 "exceeded torque tolerance: error=%f\n", error);
3025 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE 3026 printf(
"Error in net torque: %f\n", error);
3049 void HomePatch::reposition_colinear_lonepair(
3053 BigReal distance = distance_f;
3057 if (r2 < 1e-10 || 100. < r2) {
3058 iout <<
iWARN <<
"Large/small distance between lonepair reference atoms: (" 3059 << rj <<
") (" << rk <<
")\n" <<
endi;
3061 ri = rj + (scale + distance*r_jk.
rlength())*r_jk;
3078 void HomePatch::reposition_relative_lonepair(
3082 if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
3083 iout <<
iWARN <<
"Large distance between lonepair reference atoms: (" 3084 << rj <<
") (" << rk <<
") (" << rl <<
")\n" <<
endi;
3086 BigReal r, t,
p, cst, snt, csp, snp, invlen;
3089 if (distance >= 0) {
3115 ri.
x = rj.x + w.
x*a.
x + w.
y*b.
x + w.
z*c.
x;
3116 ri.
y = rj.y + w.
x*a.
y + w.
y*b.
y + w.
z*c.
y;
3117 ri.
z = rj.z + w.
x*a.
z + w.
y*b.
z + w.
z*c.
z;
3131 ri.
x = (mi * ri_old.
x + mj * rj_old.
x)/(mi + mj);
3132 ri.
y = (mi * ri_old.
y + mj * rj_old.
y)/(mi + mj);
3133 ri.
z = (mi * ri_old.
z + mj * rj_old.
z)/(mi + mj);
3139 void HomePatch::reposition_all_lonepairs(
void) {
3142 if (atom[i].mass < 0.01) {
3148 sprintf(errmsg,
"reposition lone pairs: " 3149 "no Lphost exists for LP %d\n", aid);
3157 sprintf(errmsg,
"reposition lone pairs: " 3158 "LP %d has some Lphost atom off patch\n", aid);
3163 reposition_colinear_lonepair(atom[i].position, atom[j.
index].position,
3167 reposition_relative_lonepair(atom[i].position, atom[j.
index].position,
3168 atom[k.
index].position, atom[l.
index].position,
3175 void HomePatch::reposition_all_alchpairs(
void) {
3179 for (
int i = 0; i <
numAtoms; i++) {
3181 alchPair_id = atom[i].id + numFepInitial;
3183 reposition_alchpair(atom[i].position, atom[j.
index].position, atom[i].mass, atom[j.
index].mass);
3193 pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
3197 vel[2] = (pos[2] - ref[2]) * invdt;
3215 pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
3220 vel[3] = (pos[3] - ref[3]) * invdt;
3227 void HomePatch::redistrib_lonepair_forces(
const int ftag,
Tensor *virial) {
3230 for (
int i = 0; i <
numAtoms; i++) {
3231 if (atom[i].mass < 0.01) {
3237 sprintf(errmsg,
"redistrib lone pair forces: " 3238 "no Lphost exists for LP %d\n", aid);
3246 sprintf(errmsg,
"redistrib lone pair forces: " 3247 "LP %d has some Lphost atom off patch\n", aid);
3252 redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
3253 f_mod[ftag][k.
index], atom[i].position, atom[j.
index].position,
3258 redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.
index],
3260 atom[i].position, atom[j.
index].position,
3261 atom[k.
index].position, atom[l.
index].position, virial, midpt);
3267 void HomePatch::redistrib_alchpair_forces(
const int ftag) {
3274 for (
int i = 0; i <
numAtoms; i++) {
3276 alchPair_id = atom[i].id + numFepInitial;
3278 redistrib_ap_force(f_mod[ftag][i],f_mod[ftag][j.
index]);
3283 void HomePatch::redistrib_swm4_forces(
const int ftag,
Tensor *virial) {
3287 for (
int i = 0; i <
numAtoms; i++) {
3288 if (atom[i].mass < 0.01) {
3290 redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
3291 f_mod[ftag][i+2], f_mod[ftag][i],
3292 atom[i-2].position, atom[i+1].position,
3293 atom[i+2].position, atom[i].position, virial);
3298 void HomePatch::redistrib_tip4p_forces(
const int ftag,
Tensor* virial) {
3304 if (atom[i].mass < 0.01) {
3306 redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
3307 f_mod[ftag][i-1], f_mod[ftag][i],
3308 atom[i-3].position, atom[i-2].position,
3309 atom[i-1].position, atom[i].position, virial);
3317 const Force * __restrict force_arr,
3324 for (
int i = 0; i < num_atoms; ++i ) {
3325 if ( atom_arr[i].atomFixed ) {
3326 atom_arr[i].velocity = 0;
3328 BigReal dt_mass = dt * atom_arr[i].recipMass;
3329 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
3330 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
3331 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
3335 for (
int i = 0; i < num_atoms; ++i ) {
3336 BigReal dt_mass = dt * atom_arr[i].recipMass;
3337 atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
3338 atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
3339 atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
3346 const Force * __restrict force_arr1,
3347 const Force * __restrict force_arr2,
3348 const Force * __restrict force_arr3,
3357 for (
int i = 0; i < num_atoms; ++i ) {
3358 if ( atom_arr[i].atomFixed ) {
3359 atom_arr[i].velocity = 0;
3361 BigReal rmass = atom_arr[i].recipMass;
3362 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
3363 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
3364 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
3365 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
3366 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
3367 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
3371 for (
int i = 0; i < num_atoms; ++i ) {
3372 BigReal rmass = atom_arr[i].recipMass;
3373 atom_arr[i].velocity.x += (force_arr1[i].x*dt1
3374 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
3375 atom_arr[i].velocity.y += (force_arr1[i].y*dt1
3376 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
3377 atom_arr[i].velocity.z += (force_arr1[i].z*dt1
3378 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
3390 for (
int i = 0; i < num_atoms; ++i ) {
3391 if ( ! atom_arr[i].atomFixed ) {
3392 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
3393 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
3394 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
3398 for (
int i = 0; i < num_atoms; ++i ) {
3399 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
3400 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
3401 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
3412 const int fixedAtomsOn =
simParams->fixedAtomsOn;
3414 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
3418 BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
3424 double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
3425 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;
3426 double dot_v_r_1, dot_v_r_2;
3427 double vb_cm, dr_a, dr_b;
3433 nslabs =
simParams->pressureProfileSlabs;
3439 r_wall_SQ = r_wall*r_wall;
3442 if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) {
3446 v_ab = atom[ib].position - atom[ia].position;
3447 rab_SQ = v_ab.
x*v_ab.
x + v_ab.
y*v_ab.
y + v_ab.
z*v_ab.
z;
3449 if (rab_SQ > r_wall_SQ) {
3451 if ( (rab > (2.0*r_wall)) && dieOnError ) {
3453 <<
"The drude is too far away from atom " 3454 << (atom[ia].id + 1) <<
" d = " << rab <<
"!\n" <<
endi;
3462 if ( fixedAtomsOn && atom[ia].atomFixed ) {
3463 if (atom[ib].atomFixed) {
3467 dot_v_r_2 = atom[ib].velocity.x*v_ab.
x 3468 + atom[ib].velocity.y*v_ab.
y + atom[ib].velocity.z*v_ab.
z;
3469 vb_2 = dot_v_r_2 * v_ab;
3470 vp_2 = atom[ib].velocity - vb_2;
3473 if(dot_v_r_2 == 0.0) {
3477 delta_T = dr/fabs(dot_v_r_2);
3478 if(delta_T > maxtime ) delta_T = maxtime;
3481 dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
3483 vb_2 = dot_v_r_2 * v_ab;
3485 new_vel_a = atom[ia].velocity;
3486 new_vel_b = vp_2 + vb_2;
3488 dr_b = -dr + delta_T*dot_v_r_2;
3490 new_pos_a = atom[ia].position;
3491 new_pos_b = atom[ib].position + dr_b*v_ab;
3495 mass_a = atom[ia].mass;
3496 mass_b = atom[ib].mass;
3497 mass_sum = mass_a+mass_b;
3499 dot_v_r_1 = atom[ia].velocity.x*v_ab.
x 3500 + atom[ia].velocity.y*v_ab.
y + atom[ia].velocity.z*v_ab.
z;
3501 vb_1 = dot_v_r_1 * v_ab;
3502 vp_1 = atom[ia].velocity - vb_1;
3504 dot_v_r_2 = atom[ib].velocity.x*v_ab.x
3505 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
3506 vb_2 = dot_v_r_2 * v_ab;
3507 vp_2 = atom[ib].velocity - vb_2;
3509 vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
3516 if(dot_v_r_2 == dot_v_r_1) {
3520 delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);
3521 if(delta_T > maxtime ) delta_T = maxtime;
3525 v_Bond = sqrt(kbt/mass_b);
3528 dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
3529 dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
3531 dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
3532 dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
3534 new_pos_a = atom[ia].position + dr_a*v_ab;
3535 new_pos_b = atom[ib].position + dr_b*v_ab;
3542 vb_1 = dot_v_r_1 * v_ab;
3543 vb_2 = dot_v_r_2 * v_ab;
3545 new_vel_a = vp_1 + vb_1;
3546 new_vel_b = vp_2 + vb_2;
3551 atom[ia].position = new_pos_a;
3552 atom[ib].position = new_pos_b;
3554 else if ( virial == 0 ) {
3555 atom[ia].velocity = new_vel_a;
3556 atom[ib].velocity = new_vel_b;
3559 for ( j = 0; j < 2; j++ ) {
3562 new_pos = &new_pos_a;
3563 new_vel = &new_vel_a;
3567 new_pos = &new_pos_b;
3568 new_vel = &new_vel_b;
3570 Force df = (*new_vel - atom[Idx].velocity) *
3571 ( atom[Idx].mass * invdt );
3574 atom[Idx].velocity = *new_vel;
3575 atom[Idx].position = *new_pos;
3581 int slab = (int)floor((z-zmin)*idz);
3582 if (slab < 0) slab += nslabs;
3583 else if (slab >= nslabs) slab -= nslabs;
3586 ppreduction->
item(ppoffset ) += vir.
xx;
3587 ppreduction->
item(ppoffset+1) += vir.
yy;
3588 ppreduction->
item(ppoffset+2) += vir.
zz;
3607 if ( dt && virial ) *virial += wc;
3613 #ifdef DEBUG_MINIMIZE 3615 printf(
"Step %d, patch %d: buildRattleList()\n",
3620 const int fixedAtomsOn =
simParams->fixedAtomsOn;
3621 const int useSettle =
simParams->useSettle;
3636 if ( ! settle_initialized ) {
3637 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3639 if (atom[ig].rigidBondLength > 0) {
3649 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
3656 atom[ig].rigidBondLength,
3657 atom[oatm].rigidBondLength,
3658 settle_mO, settle_mH,
3659 settle_mOrmT, settle_mHrmT, settle_ra,
3660 settle_rb, settle_rc, settle_rra);
3661 settle_initialized = 1;
3676 int posRattleParam = 0;
3683 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3684 int hgs = atom[ig].hydrogenGroupSize;
3691 for (
int i = 0; i < hgs; ++i ) {
3692 ref[i] = atom[ig+i].position;
3693 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
3694 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
3701 BigReal tmp = atom[ig].rigidBondLength;
3703 if (hgs != wathgsize) {
3705 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms " 3706 "but the specified water model requires %d atoms.\n",
3707 atom[ig].
id+1, hgs, wathgsize);
3711 if (useSettle && !anyfixed) {
3716 if ( !(fixed[1] && fixed[2]) ) {
3717 dsq[icnt] = tmp * tmp;
3723 for (
int i = 1; i < hgs; ++i ) {
3724 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3725 if ( !(fixed[0] && fixed[i]) ) {
3726 dsq[icnt] = tmp * tmp;
3740 rattleListElem.
ig = ig;
3741 rattleListElem.
icnt = icnt;
3743 for (
int i = 0; i < icnt; ++i ) {
3747 rattleParamElem.
ia = a;
3748 rattleParamElem.
ib = b;
3749 rattleParamElem.
dsq = dsq[i];
3750 rattleParamElem.
rma = rmass[a];
3751 rattleParamElem.
rmb = rmass[b];
3755 for (
int i = icnt; i < 4; ++i ) {
3757 rattleParamElem.
ia = 0;
3758 rattleParamElem.
ib = 0;
3759 rattleParamElem.
dsq = 0;
3760 rattleParamElem.
rma = 0;
3761 rattleParamElem.
rmb = 0;
3765 for(
int i = 0; i < 4; ++i) {
3768 std::cout << std::endl;
3775 for (
int ig = 0; ig <
numAtoms; ++ig ) {
3776 Force df = (
velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
3780 atom[ig].velocity =
velNew[ig];
3790 return rattle1old(timestep, virial, ppreduction);
3798 const int fixedAtomsOn =
simParams->fixedAtomsOn;
3799 const int useSettle =
simParams->useSettle;
3801 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
3812 for (
int j=0;j < n;j+=2) {
3815 for (
int i = 0; i < 3; ++i ) {
3816 ref[i] = atom[ig+i].position;
3817 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
3820 for (
int i = 0; i < 3; ++i ) {
3821 ref[i+3] = atom[ig+i].position;
3822 pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
3825 settle_mOrmT, settle_mHrmT, settle_ra,
3826 settle_rb, settle_rc, settle_rra);
3829 for (
int i = 0; i < 3; ++i ) {
3830 velNew[ig+i] = (pos[i] - ref[i])*invdt;
3834 for (
int i = 0; i < 3; ++i ) {
3835 velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
3843 for (
int i = 0; i < 3; ++i ) {
3844 ref[i] = atom[ig+i].position;
3845 pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
3848 settle_mOrmT, settle_mHrmT, settle_ra,
3849 settle_rb, settle_rc, settle_rra);
3850 for (
int i = 0; i < 3; ++i ) {
3851 velNew[ig+i] = (pos[i] - ref[i])*invdt;
3869 int hgs = atom[ig].hydrogenGroupSize;
3870 for (
int i = 0; i < hgs; ++i ) {
3871 ref[i] = atom[ig+i].position;
3872 pos[i] = atom[ig+i].position;
3873 if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
3874 pos[i] += atom[ig+i].velocity * dt;
3904 tol2, maxiter, done, consFailure);
3917 for (
int i = 0; i < hgs; ++i ) {
3923 for (
int i = 0; i < hgs; ++i ) {
3924 velNew[ig+i] = (pos[i] - ref[i])*invdt;
3928 if ( consFailure ) {
3930 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom " 3931 << (atom[ig].id + 1) <<
"!\n" <<
endi;
3934 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom " 3935 << (atom[ig].id + 1) <<
"!\n" <<
endi;
3937 }
else if ( ! done ) {
3939 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom " 3940 << (atom[ig].id + 1) <<
"!\n" <<
endi;
3943 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom " 3944 << (atom[ig].id + 1) <<
"!\n" <<
endi;
3952 int hgs = atom[ig].hydrogenGroupSize;
3953 for (
int i = 0; i < hgs; ++i ) {
3954 velNew[ig+i] = atom[ig+i].velocity;
3955 posNew[ig+i] = atom[ig+i].position;
3960 for (
int ig = 0; ig <
numAtoms; ++ig )
3961 atom[ig].position =
posNew[ig];
3962 }
else if ( virial == 0 ) {
3963 for (
int ig = 0; ig <
numAtoms; ++ig )
3964 atom[ig].velocity =
velNew[ig];
3981 const int fixedAtomsOn =
simParams->fixedAtomsOn;
3982 const int useSettle =
simParams->useSettle;
3984 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
3990 int ial[10], ibl[10];
4011 if ( ! settle_initialized ) {
4012 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
4014 if (atom[ig].rigidBondLength > 0) {
4024 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
4031 atom[ig].rigidBondLength,
4032 atom[oatm].rigidBondLength,
4033 settle_mO, settle_mH,
4034 settle_mOrmT, settle_mHrmT, settle_ra,
4035 settle_rb, settle_rc, settle_rra);
4036 settle_initialized = 1;
4043 nslabs =
simParams->pressureProfileSlabs;
4048 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
4049 int hgs = atom[ig].hydrogenGroupSize;
4050 if ( hgs == 1 )
continue;
4053 for ( i = 0; i < hgs; ++i ) {
4054 ref[i] = atom[ig+i].position;
4055 pos[i] = atom[ig+i].position;
4056 vel[i] = atom[ig+i].velocity;
4057 rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
4059 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
4062 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
else pos[i] += vel[i] * dt;
4065 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
4066 if (hgs != wathgsize) {
4068 sprintf(errmsg,
"Water molecule starting with atom %d contains %d atoms " 4069 "but the specified water model requires %d atoms.\n",
4070 atom[ig].
id+1, hgs, wathgsize);
4074 if (useSettle && !anyfixed) {
4085 settle1(ref+2, pos+2, vel+2, invdt,
4086 settle_mOrmT, settle_mHrmT, settle_ra,
4087 settle_rb, settle_rc, settle_rra);
4095 swm4_omrepos(ref, pos, vel, invdt);
4099 settle_mOrmT, settle_mHrmT, settle_ra,
4100 settle_rb, settle_rc, settle_rra);
4102 tip4_omrepos(ref, pos, vel, invdt);
4109 if ( invdt == 0 )
for ( i = 0; i < wathgsize; ++i ) {
4110 atom[ig+i].position = pos[i];
4111 }
else if ( virial == 0 )
for ( i = 0; i < wathgsize; ++i ) {
4112 atom[ig+i].velocity = vel[i];
4113 }
else for ( i = 0; i < wathgsize; ++i ) {
4114 Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
4118 atom[ig+i].velocity = vel[i];
4125 int slab = (int)floor((z-zmin)*idz);
4126 if (slab < 0) slab += nslabs;
4127 else if (slab >= nslabs) slab -= nslabs;
4130 ppreduction->
item(ppoffset ) += vir.
xx;
4131 ppreduction->
item(ppoffset+1) += vir.
yy;
4132 ppreduction->
item(ppoffset+2) += vir.
zz;
4137 if ( !(fixed[1] && fixed[2]) ) {
4138 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
4141 for ( i = 1; i < hgs; ++i ) {
4142 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
4143 if ( !(fixed[0] && fixed[i]) ) {
4144 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
4148 if ( icnt == 0 )
continue;
4149 for ( i = 0; i < icnt; ++i ) {
4150 refab[i] = ref[ial[i]] - ref[ibl[i]];
4152 for ( i = 0; i < hgs; ++i ) {
4157 for ( iter = 0; iter < maxiter; ++iter ) {
4161 for ( i = 0; i < icnt; ++i ) {
4162 int a = ial[i];
int b = ibl[i];
4163 Vector pab = pos[a] - pos[b];
4166 BigReal diffsq = rabsq - pabsq;
4167 if ( fabs(diffsq) > (rabsq * tol2) ) {
4170 if ( rpab < ( rabsq * 1.0e-6 ) ) {
4177 BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
4181 if ( invdt != 0. ) {
4192 if ( consFailure ) {
4194 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom " 4195 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4198 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom " 4199 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4201 }
else if ( ! done ) {
4203 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom " 4204 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4207 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom " 4208 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4214 if ( invdt == 0 )
for ( i = 0; i < hgs; ++i ) {
4215 atom[ig+i].position = pos[i];
4216 }
else if ( virial == 0 )
for ( i = 0; i < hgs; ++i ) {
4217 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
4218 }
else for ( i = 0; i < hgs; ++i ) {
4219 Force df = netdp[i] * invdt;
4223 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
4228 int slab = (int)floor((z-zmin)*idz);
4229 if (slab < 0) slab += nslabs;
4230 else if (slab >= nslabs) slab -= nslabs;
4233 ppreduction->
item(ppoffset ) += vir.
xx;
4234 ppreduction->
item(ppoffset+1) += vir.
yy;
4235 ppreduction->
item(ppoffset+2) += vir.
zz;
4239 if ( dt && virial ) *virial += wc;
4249 const int fixedAtomsOn =
simParams->fixedAtomsOn;
4250 const int useSettle =
simParams->useSettle;
4258 int ial[10], ibl[10];
4271 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
4273 int hgs = atom[ig].hydrogenGroupSize;
4274 if ( hgs == 1 )
continue;
4277 for ( i = 0; i < hgs; ++i ) {
4278 ref[i] = atom[ig+i].position;
4279 vel[i] = atom[ig+i].velocity;
4280 rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
4281 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
4282 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
4285 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
4286 if (hgs != wathgsize) {
4287 NAMD_bug(
"Hydrogen group error caught in rattle2().");
4290 if (useSettle && !anyfixed) {
4299 settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
4306 settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
4311 for (i=0; i<hgs; i++) {
4312 atom[ig+i].velocity = vel[i];
4316 if ( !(fixed[1] && fixed[2]) ) {
4317 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
4318 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
4322 for ( i = 1; i < hgs; ++i ) {
4323 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
4324 if ( !(fixed[0] && fixed[i]) ) {
4325 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
4326 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
4327 ibl[icnt] = i; ++icnt;
4331 if ( icnt == 0 )
continue;
4333 for ( i = 0; i < icnt; ++i ) {
4334 refab[i] = ref[ial[i]] - ref[ibl[i]];
4338 for ( iter = 0; iter < maxiter; ++iter ) {
4340 for ( i = 0; i < icnt; ++i ) {
4341 int a = ial[i];
int b = ibl[i];
4342 Vector vab = vel[a] - vel[b];
4346 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
4347 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
4348 wc +=
outer(dp,rab);
4349 vel[a] += rmass[a] * dp;
4350 vel[b] -= rmass[b] * dp;
4359 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
4362 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
4366 for ( i = 0; i < hgs; ++i ) {
4367 atom[ig+i].velocity = vel[i];
4372 *virial += wc / ( 0.5 * dt );
4383 const int fixedAtomsOn =
simParams->fixedAtomsOn;
4384 const int useSettle =
simParams->useSettle;
4392 int ial[10], ibl[10];
4405 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
4407 int hgs = atom[ig].hydrogenGroupSize;
4408 if ( hgs == 1 )
continue;
4411 for ( i = 0; i < hgs; ++i ) {
4412 ref[i] = atom[ig+i].position;
4413 vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
4415 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
4416 if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
4419 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
4420 if (hgs != wathgsize) {
4421 NAMD_bug(
"Hydrogen group error caught in rattle2().");
4424 if (useSettle && !anyfixed) {
4433 settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
4440 settle2(1.0, 1.0, ref, vel, dt, virial);
4445 for (i=0; i<hgs; i++) {
4446 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
4450 if ( !(fixed[1] && fixed[2]) ) {
4451 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
4452 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
4456 for ( i = 1; i < hgs; ++i ) {
4457 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
4458 if ( !(fixed[0] && fixed[i]) ) {
4459 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
4460 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
4461 ibl[icnt] = i; ++icnt;
4465 if ( icnt == 0 )
continue;
4467 for ( i = 0; i < icnt; ++i ) {
4468 refab[i] = ref[ial[i]] - ref[ibl[i]];
4472 for ( iter = 0; iter < maxiter; ++iter ) {
4474 for ( i = 0; i < icnt; ++i ) {
4475 int a = ial[i];
int b = ibl[i];
4476 Vector vab = vel[a] - vel[b];
4480 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
4481 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
4482 wc +=
outer(dp,rab);
4483 vel[a] += rmass[a] * dp;
4484 vel[b] -= rmass[b] * dp;
4493 NAMD_die(
"Exceeded maximum number of iterations in rattle2().");
4496 "Exceeded maximum number of iterations in rattle2().\n" <<
endi;
4500 for ( i = 0; i < hgs; ++i ) {
4501 ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
4506 *virial += wc / ( 0.5 * dt );
4517 #ifdef DEBUG_MINIMIZE 4519 printf(
"Step %d, patch %d: buildRattleList_SOA()\n",
4528 const double * __restrict pos_x = patchDataSOA.
pos_x;
4529 const double * __restrict pos_y = patchDataSOA.
pos_y;
4530 const double * __restrict pos_z = patchDataSOA.
pos_z;
4531 const float * __restrict mass = patchDataSOA.
mass;
4532 const double * __restrict recipMass = patchDataSOA.
recipMass;
4533 const float * __restrict rigidBondLength = patchDataSOA.
rigidBondLength;
4537 const int fixedAtomsOn =
simParams->fixedAtomsOn;
4538 const int useSettle =
simParams->useSettle;
4549 if ( ! settle_initialized ) {
4550 for (
int ig = numSoluteAtoms;
4552 ig += hydrogenGroupSize[ig]) {
4554 if (rigidBondLength[ig] > 0) {
4564 if (mass[ig] < 0.5 || mass[ig+1] < 0.5) {
4571 rigidBondLength[ig],
4572 rigidBondLength[oatm],
4573 settle_mO, settle_mH,
4574 settle_mOrmT, settle_mHrmT, settle_ra,
4575 settle_rb, settle_rc, settle_rra);
4576 settle_initialized = 1;
4590 for (
int ig = 0; ig < numSoluteAtoms; ig += hydrogenGroupSize[ig]) {
4591 int hgs = hydrogenGroupSize[ig];
4599 BigReal tmp = rigidBondLength[ig];
4601 dsq[icnt] = tmp * tmp;
4606 for (
int i = 1; i < hgs; ++i ) {
4607 if ( ( tmp = rigidBondLength[ig+i] ) > 0 ) {
4608 dsq[icnt] = tmp * tmp;
4621 rattleListElem.
ig = ig;
4622 rattleListElem.
icnt = icnt;
4624 for (
int i = 0; i < icnt; ++i ) {
4628 rattleParamElem.
ia = a;
4629 rattleParamElem.
ib = b;
4630 rattleParamElem.
dsq = dsq[i];
4631 rattleParamElem.
rma = recipMass[ig+a];
4632 rattleParamElem.
rmb = recipMass[ig+b];
4636 for (
int i = icnt; i < 4; ++i )
4639 rattleParamElem.
ia = 0;
4640 rattleParamElem.
ib = 0;
4641 rattleParamElem.
dsq = 0;
4642 rattleParamElem.
rma = 0;
4643 rattleParamElem.
rmb = 0;
4665 double * __restrict pos_x = patchDataSOA.
pos_x;
4666 double * __restrict pos_y = patchDataSOA.
pos_y;
4667 double * __restrict pos_z = patchDataSOA.
pos_z;
4668 double * __restrict vel_x = patchDataSOA.
vel_x;
4669 double * __restrict vel_y = patchDataSOA.
vel_y;
4670 double * __restrict vel_z = patchDataSOA.
vel_z;
4671 double * __restrict posNew_x = patchDataSOA.
posNew_x;
4672 double * __restrict posNew_y = patchDataSOA.
posNew_y;
4673 double * __restrict posNew_z = patchDataSOA.
posNew_z;
4674 double * __restrict velNew_x = patchDataSOA.
velNew_x;
4675 double * __restrict velNew_y = patchDataSOA.
velNew_y;
4676 double * __restrict velNew_z = patchDataSOA.
velNew_z;
4678 #ifdef __INTEL_COMPILER 4679 __assume_aligned(pos_x,64);
4680 __assume_aligned(pos_y,64);
4681 __assume_aligned(pos_z,64);
4682 __assume_aligned(vel_x,64);
4683 __assume_aligned(vel_y,64);
4684 __assume_aligned(vel_z,64);
4685 __assume_aligned(posNew_x,64);
4686 __assume_aligned(posNew_y,64);
4687 __assume_aligned(posNew_z,64);
4688 __assume_aligned(velNew_x,64);
4689 __assume_aligned(velNew_y,64);
4690 __assume_aligned(velNew_z,64);
4691 __assume_aligned(hydrogenGroupSize,64);
4694 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
4701 posNew_x[i] = pos_x[i] + vel_x[i] * dt;
4702 posNew_y[i] = pos_y[i] + vel_y[i] * dt;
4703 posNew_z[i] = pos_z[i] + vel_z[i] * dt;
4709 if (numSolventAtoms > 0) {
4710 int n = numSoluteAtoms;
4712 &posNew_x[n], &posNew_y[n], &posNew_z[n],
4714 settle_mOrmT, settle_mHrmT, settle_ra,
4715 settle_rb, settle_rc, settle_rra);
4725 &pos_x[ig], &pos_y[ig], &pos_z[ig],
4726 &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
4734 &pos_x[ig], &pos_y[ig], &pos_z[ig],
4735 &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
4741 &pos_x[ig], &pos_y[ig], &pos_z[ig],
4742 &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
4743 tol2, maxiter, done, consFailure);
4748 &pos_x[ig], &pos_y[ig], &pos_z[ig],
4749 &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
4757 if ( consFailure ) {
4759 iout <<
iERROR <<
"Constraint failure in RATTLE algorithm for atom " 4760 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4763 iout <<
iWARN <<
"Constraint failure in RATTLE algorithm for atom " 4764 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4766 }
else if ( ! done ) {
4768 iout <<
iERROR <<
"Exceeded RATTLE iteration limit for atom " 4769 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4772 iout <<
iWARN <<
"Exceeded RATTLE iteration limit for atom " 4773 << (atom[ig].id + 1) <<
"!\n" <<
endi;
4781 velNew_x[i] = (posNew_x[i] - pos_x[i]) * invdt;
4782 velNew_y[i] = (posNew_y[i] - pos_y[i]) * invdt;
4783 velNew_z[i] = (posNew_z[i] - pos_z[i]) * invdt;
4789 for (
int j=0; j < numNoconst; j++) {
4791 posNew_x[ig] = pos_x[ig];
4792 posNew_y[ig] = pos_y[ig];
4793 posNew_z[ig] = pos_z[ig];
4794 velNew_x[ig] = vel_x[ig];
4795 velNew_y[ig] = vel_y[ig];
4796 velNew_z[ig] = vel_z[ig];
4800 for (
int ig = 0; ig <
numAtoms; ++ig ) {
4801 pos_x[ig] = posNew_x[ig];
4802 pos_y[ig] = posNew_y[ig];
4803 pos_z[ig] = posNew_z[ig];
4806 else if ( virial == 0 ) {
4807 for (
int ig = 0; ig <
numAtoms; ++ig ) {
4808 vel_x[ig] = velNew_x[ig];
4809 vel_y[ig] = velNew_y[ig];
4810 vel_z[ig] = velNew_z[ig];
4814 const float * __restrict mass = patchDataSOA.
mass;
4815 double * __restrict f_normal_x = patchDataSOA.
f_normal_x;
4816 double * __restrict f_normal_y = patchDataSOA.
f_normal_y;
4817 double * __restrict f_normal_z = patchDataSOA.
f_normal_z;
4818 #ifdef __INTEL_COMPILER 4819 __assume_aligned(mass,64);
4820 __assume_aligned(f_normal_x,64);
4821 __assume_aligned(f_normal_y,64);
4822 __assume_aligned(f_normal_z,64);
4825 for (
int ig = 0; ig <
numAtoms; ig++) {
4826 BigReal df_x = (velNew_x[ig] - vel_x[ig]) * ( mass[ig] * invdt );
4827 BigReal df_y = (velNew_y[ig] - vel_y[ig]) * ( mass[ig] * invdt );
4828 BigReal df_z = (velNew_z[ig] - vel_z[ig]) * ( mass[ig] * invdt );
4829 f_normal_x[ig] += df_x;
4830 f_normal_y[ig] += df_y;
4831 f_normal_z[ig] += df_z;
4832 vir.
xx += df_x * pos_x[ig];
4833 vir.
xy += df_x * pos_y[ig];
4834 vir.
xz += df_x * pos_z[ig];
4835 vir.
yx += df_y * pos_x[ig];
4836 vir.
yy += df_y * pos_y[ig];
4837 vir.
yz += df_y * pos_z[ig];
4838 vir.
zx += df_z * pos_x[ig];
4839 vir.
zy += df_z * pos_y[ig];
4840 vir.
zz += df_z * pos_z[ig];
4843 for (
int ig = 0; ig <
numAtoms; ig++) {
4844 vel_x[ig] = velNew_x[ig];
4845 vel_y[ig] = velNew_y[ig];
4846 vel_z[ig] = velNew_z[ig];
4862 DebugM(2,
"loweAndersenVelocities\n");
4866 for (
int i = 0; i <
numAtoms; ++i) {
4869 v[i].position = atom[i].velocity;
4870 v[i].charge = atom[i].mass;
4872 DebugM(2,
"loweAndersenVelocities\n");
4877 DebugM(2,
"loweAndersenFinish\n");
4888 for (
int i = 0; i <
numAtoms; i++) {
4900 for (
int i = 0; i <
numAtoms; i++) {
4903 intRad[2*i+0] = rad - offset;
4904 intRad[2*i+1] = screen*(rad - offset);
4921 for (
int i = 0; i <
numAtoms; i++) {
4923 rhoi = rhoi0+coulomb_radius_offset;
4949 BigReal fij, expkappa, Dij, dEdai, dedasum;
4950 BigReal rhoi, rhoi0, psii, nbetapsi;
4951 BigReal gammapsi2, tanhi, daidr;
4952 for (
int i = 0; i <
numAtoms; i++) {
4956 expkappa = exp(-kappa*fij);
4957 Dij = epsilon_p_i - expkappa*epsilon_s_i;
4959 dEdai = -0.5*
COULOMB*atom[i].charge*atom[i].charge
4960 *(kappa*epsilon_s_i*expkappa-Dij/fij)/
bornRad[i];
4965 rhoi = rhoi0+coulomb_radius_offset;
4967 nbetapsi = -beta*psii;
4968 gammapsi2 = gamma*psii*psii;
4969 tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
4971 * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
4982 if (proxy.
size() > 0) {
4983 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
4984 for (
int i = 0; i < proxy.
size(); i++) {
4985 int node = proxy[i];
4994 cp[node].recvData(msg);
5002 if (proxy.
size() > 0) {
5003 CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
5006 for (
int i = 0; i < proxy.
size(); i++) {
5007 int node = proxy[i];
5017 cp[node].recvData(msg);
5026 for (
int i = 0; i < msg->
psiSumLen; ++i ) {
5033 if (phase1BoxClosedCalled ==
true &&
5034 numGBISP1Arrived==proxy.
size() ) {
5041 numGBISP1Arrived == proxy.
size() &&
5042 numGBISP2Arrived == proxy.
size() &&
5043 numGBISP3Arrived == proxy.
size()) {
5054 for (
int i = 0; i < msg->
dEdaSumLen; ++i ) {
5061 if (phase2BoxClosedCalled ==
true &&
5062 numGBISP2Arrived==proxy.
size() ) {
5069 numGBISP1Arrived == proxy.
size() &&
5070 numGBISP2Arrived == proxy.
size() &&
5071 numGBISP3Arrived == proxy.
size()) {
5101 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
5102 int hgs = atom[ig].hydrogenGroupSize;
5103 if ( hgs == 1 )
continue;
5104 for ( i = 0; i < hgs; ++i ) {
5105 ref[i] = atom[ig+i].position;
5106 rmass[i] = 1. / atom[ig+i].mass;
5107 fixed[i] = (
simParams->fixedAtomsOn && atom[ig+i].atomFixed );
5108 if ( fixed[i] ) rmass[i] = 0.;
5113 if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {
5115 NAMD_die(
"Hydrogen group error caught in mollyAverage(). It's a bug!\n");
5117 if ( !(fixed[1] && fixed[2]) ) {
5118 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
5121 for ( i = 1; i < hgs; ++i ) {
5122 if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
5123 if ( !(fixed[0] && fixed[i]) ) {
5124 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
5128 if ( icnt == 0 )
continue;
5130 molly_lambda.
resize(numLambdas);
5131 lambda = &(molly_lambda[numLambdas - icnt]);
5132 for ( i = 0; i < icnt; ++i ) {
5133 refab[i] = ref[ial[i]] - ref[ibl[i]];
5136 iter=
average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
5137 if ( iter == maxiter ) {
5138 iout <<
iWARN <<
"Exceeded maximum number of iterations in mollyAverage().\n"<<
endi;
5169 for (
int ig = 0; ig <
numAtoms; ig += atom[ig].hydrogenGroupSize ) {
5170 int hgs = atom[ig].hydrogenGroupSize;
5171 if (hgs == 1 )
continue;
5172 for ( i = 0; i < hgs; ++i ) {
5173 ref[i] = atom[ig+i].position;
5175 rmass[i] = 1. / atom[ig+i].mass;
5176 fixed[i] = (
simParams->fixedAtomsOn && atom[ig+i].atomFixed );
5177 if ( fixed[i] ) rmass[i] = 0.;
5181 if ( atom[ig].rigidBondLength > 0 ) {
5183 NAMD_die(
"Hydrogen group error caught in mollyMollify(). It's a bug!\n");
5185 if ( !(fixed[1] && fixed[2]) ) {
5186 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
5189 for ( i = 1; i < hgs; ++i ) {
5190 if ( atom[ig+i].rigidBondLength > 0 ) {
5191 if ( !(fixed[0] && fixed[i]) ) {
5192 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
5197 if ( icnt == 0 )
continue;
5198 lambda = &(molly_lambda[numLambdas]);
5200 for ( i = 0; i < icnt; ++i ) {
5201 refab[i] = ref[ial[i]] - ref[ibl[i]];
5204 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
5206 for ( i = 0; i < hgs; ++i ) {
5217 checkpoint_atom.
copy(atom);
5221 #if NAMD_SeparateWaters != 0 5222 checkpoint_numWaterAtoms = numWaterAtoms;
5229 atom.
copy(checkpoint_atom);
5233 doAtomUpdate =
true;
5241 sort_solvent_atoms();
5242 copy_atoms_to_SOA();
5252 #if NAMD_SeparateWaters != 0 5253 numWaterAtoms = checkpoint_numWaterAtoms;
5260 const int remote =
simParams->scriptIntArg1;
5261 const char *key =
simParams->scriptStringArg1;
5268 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
5278 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
5290 msg->
replica = CmiMyPartition();
5297 const int remote =
simParams->scriptIntArg1;
5298 const char *key =
simParams->scriptStringArg1;
5300 NAMD_bug(
"HomePatch::recvCheckpointLoad called during checkpointFree.");
5302 if ( msg->
replica != remote ) {
5303 NAMD_bug(
"HomePatch::recvCheckpointLoad message from wrong replica.");
5307 strcpy(newmsg->
key,key);
5311 newmsg->
pe = CkMyPe();
5312 newmsg->
replica = CmiMyPartition();
5324 doAtomUpdate =
true;
5329 sort_solvent_atoms();
5330 copy_atoms_to_SOA();
5360 CkpvAccess(_qd)->process();
5392 CkpvAccess(_qd)->process();
5404 CkpvAccess(_qd)->process();
5405 doAtomUpdate =
true;
5411 sort_solvent_atoms();
5412 copy_atoms_to_SOA();
5435 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED) 5459 doPairlistCheck_newTolerance *= (1. -
simParams->pairlistShrink);
5460 doPairlistCheck_lattice =
lattice;
5463 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
5470 Lattice &lattice_old = doPairlistCheck_lattice;
5473 Vector center_delta = center_cur - center_old;
5477 for ( i=0; i<2; ++i ) {
5478 for (
int j=0; j<2; ++j ) {
5479 for (
int k=0; k<2; ++k ) {
5482 corner_delta -= center_delta;
5484 if ( cd > max_cd ) max_cd = cd;
5488 max_cd = sqrt(max_cd);
5493 for ( i=0; i<n; ++i ) {
5496 p_delta -= center_delta;
5498 if ( pd > max_pd ) max_pd = pd;
5500 max_pd = sqrt(max_pd);
5502 BigReal max_tol = max_pd + max_cd;
5508 if ( max_tol > ( (1. -
simParams->pairlistTrigger) *
5509 doPairlistCheck_newTolerance ) ) {
5512 doPairlistCheck_newTolerance *= (1. +
simParams->pairlistGrow);
5515 if ( max_tol > doPairlistCheck_newTolerance ) {
5518 doPairlistCheck_newTolerance = max_tol / (1. -
simParams->pairlistTrigger);
5539 double * __restrict pos_x = patchDataSOA.
pos_x;
5540 double * __restrict pos_y = patchDataSOA.
pos_y;
5541 double * __restrict pos_z = patchDataSOA.
pos_z;
5547 const int hgs = hydrogenGroupSize[j];
5550 if ( ngs > 5 ) ngs = 5;
5555 for ( i = 1; i < ngs; ++i ) {
5556 nonbondedGroupSize[j+i] = 0;
5560 BigReal r2 = dx * dx + dy * dy + dz * dz;
5561 if ( r2 > hgcut )
break;
5562 else if ( r2 > maxrad2 ) maxrad2 = r2;
5564 nonbondedGroupSize[j] = i;
5565 for ( ; i < hgs; ++i ) {
5566 nonbondedGroupSize[j+i] = 1;
5572 NAMD_bug(
"hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
5595 while ( p_i != p_e ) {
5599 if ( ngs > 5 ) ngs = 5;
5604 for ( i = 1; i < ngs; ++i ) {
5609 BigReal r2 = dx * dx + dy * dy + dz * dz;
5610 if ( r2 > hgcut )
break;
5611 else if ( r2 > maxrad2 ) maxrad2 = r2;
5614 for ( ; i < hgs; ++i ) {
5621 NAMD_bug(
"hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
5643 if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
5644 ( bAwayDist*sysdimb < minSize*0.9999 ) ||
5645 ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
5647 NAMD_die(
"Periodic cell has become too small for original patch grid!\n" 5648 "Possible solutions are to restart from a recent checkpoint,\n" 5649 "increase margin, or disable useFlexibleCell for liquid simulation.");
5654 BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
5655 BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
5656 BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
5658 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
5659 NAMD_die(
"Periodic cell has become too small for original patch grid!\n" 5660 "There are probably many margin violations already reported.\n" 5661 "Possible solutions are to restart from a recent checkpoint,\n" 5662 "increase margin, or disable useFlexibleCell for liquid simulation.");
5672 int xdev, ydev, zdev;
5673 int problemCount = 0;
5675 double * __restrict pos_x = patchDataSOA.
pos_x;
5676 double * __restrict pos_y = patchDataSOA.
pos_y;
5677 double * __restrict pos_z = patchDataSOA.
pos_z;
5679 Vector pos(pos_x[i],pos_y[i],pos_z[i]);
5684 if (s.
x < minx) xdev = 0;
5685 else if (maxx <= s.
x) xdev = 2;
5688 if (s.
y < miny) ydev = 0;
5689 else if (maxy <= s.
y) ydev = 2;
5692 if (s.
z < minz) zdev = 0;
5693 else if (maxz <= s.
z) zdev = 2;
5696 if (mInfo[xdev][ydev][zdev]) {
5725 if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
5726 ( bAwayDist*sysdimb < minSize*0.9999 ) ||
5727 ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
5729 NAMD_die(
"Periodic cell has become too small for original patch grid!\n" 5730 "Possible solutions are to restart from a recent checkpoint,\n" 5731 "increase margin, or disable useFlexibleCell for liquid simulation.");
5736 BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
5737 BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
5738 BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
5740 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
5741 NAMD_die(
"Periodic cell has become too small for original patch grid!\n" 5742 "There are probably many margin violations already reported.\n" 5743 "Possible solutions are to restart from a recent checkpoint,\n" 5744 "increase margin, or disable useFlexibleCell for liquid simulation.");
5754 int xdev, ydev, zdev;
5755 int problemCount = 0;
5759 for ( ; p_i != p_e; ++p_i ) {
5764 if (s.
x < minx) xdev = 0;
5765 else if (maxx <= s.
x) xdev = 2;
5768 if (s.
y < miny) ydev = 0;
5769 else if (maxy <= s.
y) ydev = 2;
5772 if (s.
z < minz) zdev = 0;
5773 else if (maxz <= s.
z) zdev = 2;
5776 if (mInfo[xdev][ydev][zdev]) {
5794 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED) 5803 for (i=0; i<numNeighbors; i++) {
5819 int xdev, ydev, zdev;
5826 #if NAMD_SeparateWaters != 0 5828 int numLostWaterAtoms = 0;
5831 while ( atom_i != atom_e ) {
5848 j+=(atom_i+j)->hydrogenGroupSize ) {
5849 pos += (atom_i+j)->position;
5864 if (s.
x < minx) xdev = 0;
5865 else if (maxx <= s.
x) xdev = 2;
5868 if (s.
y < miny) ydev = 0;
5869 else if (maxy <= s.
y) ydev = 2;
5872 if (s.
z < minz) zdev = 0;
5873 else if (maxz <= s.
z) zdev = 2;
5878 if (mInfo[xdev][ydev][zdev]) {
5883 DebugM(3,
"Migrating atom " << atom_i->
id <<
" from patch " 5890 #if NAMD_SeparateWaters != 0 5898 int atomIndex = atom_i - atom_first;
5899 if (atomIndex < numWaterAtoms)
5900 numLostWaterAtoms++;
5908 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
5916 #if NAMD_SeparateWaters != 0 5917 numWaterAtoms -= numLostWaterAtoms;
<