27 #define PRINT_FORCES 0
50 #define MIN_DEBUG_LEVEL 4
60 #define START_HPM_STEP 1000
61 #define STOP_HPM_STEP 1500
64 #define SPECIAL_PATCH_ID 91
119 void Sequencer::threadRun(
Sequencer* arg)
129 DebugM(4,
"::run() - this = " <<
this <<
"\n" );
130 thread = CthCreate((CthVoidFn)&(threadRun),(
void*)(
this),
SEQ_STK_SZ);
131 CthSetStrategyDefault(thread);
152 switch ( scriptTask ) {
208 NAMD_bug(
"Unknown task in Sequencer::algorithm");
220 char tracePrefix[20];
230 #ifdef TIMER_COLLECTION
231 TimerSet& t =
patch->timerSet;
267 const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
269 doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
276 if ( dofull )
slowFreq = fullElectFrequency;
277 const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
279 doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
288 if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed =
Results::amdf;
346 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
351 doEnergy = ! ( step % energyFrequency );
353 if ( accelMDOn && !accelMDdihe ) doEnergy=1;
355 if ( adaptTempOn ) doEnergy=1;
357 D_MSG(
"runComputeObjects()");
365 if ( staleForces || doTcl || doColvars ) {
370 D_MSG(
"newtonianVelocities()");
382 D_MSG(
"submitHalfstep()");
387 D_MSG(
"newtonianVelocities()");
396 if (doTcl || doColvars)
398 D_MSG(
"submitHalfstep()");
404 D_MSG(
"newtonianVelocities()");
410 D_MSG(
"submitReductions()");
416 traceUserEvent(eventEndOfTimeStep);
417 sprintf(traceNote,
"%s%d",tracePrefix,step);
418 traceUserSuppliedNote(traceNote);
426 bool doMultigratorRattle =
false;
440 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
444 int beginStep =
simParams->beginEventStep;
449 for ( ++step; step <= numberOfSteps; ++
step )
451 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
452 eon = epid && (beginStep < step && step <= endStep);
454 if (controlProfiling && step == beginStep) {
457 if (controlProfiling && step == endStep) {
473 newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
478 if (doMultigratorRattle)
rattle1(timestep, 1);
531 #endif // UPPER_BOUND
533 doNonbonded = !(step%nonbondedFrequency);
534 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
537 if ( zeroMomentum && doFullElectrostatics ) {
556 if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed =
Results::amdf;
559 doEnergy = ! ( step % energyFrequency );
560 if ( accelMDOn && !accelMDdihe ) doEnergy=1;
561 if ( adaptTempOn ) doEnergy=1;
585 if ( staleForces || doTcl || doColvars ) {
591 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
593 newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
602 newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
616 if (doTcl || doColvars)
626 newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
665 if(step == bstep || step == estep){
670 #ifdef MEASURE_NAMD_WITH_PAPI
672 int bstep =
simParams->papiMeasureStartStep;
673 int estep = bstep +
simParams->numPapiMeasureSteps;
674 if(step == bstep || step==estep) {
675 papiMeasureBarrier(step);
681 traceUserEvent(eventEndOfTimeStep);
682 sprintf(traceNote,
"%s%d",tracePrefix,step);
683 traceUserSuppliedNote(traceNote);
685 #endif // UPPER_BOUND
690 cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
694 if(step == START_HPM_STEP)
695 (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
697 if(step == STOP_HPM_STEP)
698 (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
704 #ifdef TIMER_COLLECTION
706 printf(
"Timer collection reporting in microseconds for "
710 #endif // TIMER_COLLECTION
724 Vector movDragVel, dragIncrement;
725 for (
int i = 0; i < numAtoms; ++i )
731 dragIncrement = movDragGlobVel * movDragVel *
dt;
745 Vector rotDragAxis, rotDragPivot, dragIncrement;
746 for (
int i = 0; i < numAtoms; ++i )
752 dAngle = rotDragGlobVel * rotDragVel *
dt;
753 rotDragAxis /= rotDragAxis.
length();
754 atomRadius = atom[i].
position - rotDragPivot;
755 dragIncrement =
cross(rotDragAxis, atomRadius) * dAngle;
782 doFullElectrostatics = dofull;
816 if ( doTcl || doColvars ) {
829 for ( ++step; step <= numberOfSteps; ++
step ) {
857 if ( doTcl || doColvars ) {
881 for (
int i = 0; i < numAtoms; ++i ) {
887 for (
int j=1; j<hgs; ++j ) {
905 for (
int i = 0; i < numAtoms; ++i ) {
908 if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) {
911 if ( fixedAtomsOn && a[i].atomFixed ) a[i].
velocity = 0;
913 if ( v2 > maxv2 ) maxv2 = v2;
919 for (
int i = 0; i < numAtoms; ++i ) {
920 if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) {
923 if ( fixedAtomsOn && a[i].atomFixed ) a[i].
velocity = 0;
925 if ( v2 > maxv2 ) maxv2 = v2;
935 for (
int i = 0; i < numAtoms; i += hgs ) {
938 if ( minChildVel < fmax2 )
continue;
939 int adjustChildren = 1;
940 for (
int j = i+1; j < (i+hgs); ++j ) {
941 if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
943 if ( adjustChildren ) {
945 for (
int j = i+1; j < (i+hgs); ++j ) {
946 if (a[i].mass < 0.01)
continue;
960 for (
int i = 0; i < numAtoms; ++i ) {
965 for (
int i = 1; i < numAtoms; ++i ) {
966 if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) {
975 for (
int i = 1; i < numAtoms; ++i ) {
976 if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) {
988 for (
int i = 0; i < numAtoms; ++i ) {
1001 for (
int i = 0; i < numAtoms; ++i ) {
1006 for (
int i = 0; i < numAtoms; ++i ) {
1022 NAMD_die(
"Cannot zero momentum when fixed atoms are present.");
1033 for (
int i = 0; i < numAtoms; ++i ) {
1034 a[i].velocity += dv * a[i].recipMass;
1035 a[i].position += dx * a[i].recipMass;
1038 for (
int i = 0; i < numAtoms; ++i ) {
1039 a[i].velocity += dv;
1040 a[i].position += dx;
1052 NAMD_bug(
"Sequencer::scalePositionsVelocities, fixed atoms not implemented");
1056 for (
int i = 0; i < numAtoms; i += hgs ) {
1061 for (
int j=0;j < hgs;++j) {
1062 m_cm += a[i+j].
mass;
1071 for (
int j=0;j < hgs;++j) {
1077 for (
int i = 0; i < numAtoms; i++) {
1094 Tensor posScaleTensor = scaleTensor;
1123 for (
int i = 0; i < numAtoms; ++i ) {
1126 virialNormal.
outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1130 for (
int i = 0; i < numAtoms; ++i ) {
1131 if (a[i].mass < 0.01)
continue;
1133 virialNormal.
outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1137 kineticEnergy *= 0.5;
1145 Vector fixForceNormal = 0;
1146 Vector fixForceNbond = 0;
1149 calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
1165 for (
int i = 0; i < numAtoms; i += hgs ) {
1171 for ( j = i; j < (i+hgs); ++j ) {
1182 NAMD_bug(
"Sequencer::multigratorPressure, this part needs to be implemented correctly");
1183 for ( j = i; j < (i+hgs); ++j ) {
1188 intVirialNormal2.
outerAdd (mass, v, dv);
1196 for ( j = i; j < (i+hgs); ++j ) {
1200 intVirialNormal2.
outerAdd(mass, v, dv);
1222 for (
int i = 0; i < numAtoms; i++) {
1233 for (
int i = 0; i < numAtoms; ++i ) {
1239 for (
int i = 0; i < numAtoms; ++i ) {
1243 kineticEnergy *= 0.5;
1244 return kineticEnergy;
1262 for (
int i = 0; i < numAtoms; i += hgs ) {
1268 for ( j = i; j < (i+hgs); ++j ) {
1273 momentumSqrSum.
outerAdd(1.0/m_cm, v_cm, v_cm);
1276 for (
int i = 0; i < numAtoms; i++) {
1277 momentumSqrSum.
outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1296 const int staleForces,
1297 const int doNonbonded,
1298 const int doFullElectrostatics)
1301 NamdProfileEvent::NEWTONIAN_VELOCITIES);
1304 if (staleForces || (doNonbonded && doFullElectrostatics)) {
1310 if (staleForces || doNonbonded)
1312 if (staleForces || doFullElectrostatics)
1339 for (
int i = 0; i < numAtoms; ++i )
1342 if ( ! dt_gamma )
continue;
1344 BigReal f1 = exp( -dt_gamma );
1345 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
1357 NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
1367 for (i = 0; i < numAtoms; i++) {
1369 if (i < numAtoms-1 &&
1370 a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1382 if (dt_gamma != 0.0) {
1383 v_com *= ( 1. - 0.5 * dt_gamma );
1388 if (dt_gamma != 0.0) {
1389 v_bnd *= ( 1. - 0.5 * dt_gamma );
1400 if ( ! dt_gamma )
continue;
1402 a[i].
velocity *= ( 1. - 0.5 * dt_gamma );
1413 for ( i = 0; i < numAtoms; ++i )
1416 if ( ! dt_gamma )
continue;
1418 a[i].
velocity *= ( 1. - 0.5 * dt_gamma );
1430 NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
1457 for (i = 0; i < numAtoms; i++) {
1459 if (i < numAtoms-1 &&
1460 a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1472 if (dt_gamma != 0.0) {
1475 sqrt( 2 * dt_gamma * kbT *
1476 ( a[i].
partition ? tempFactor : 1.0 ) / mass );
1477 v_com /= ( 1. + 0.5 * dt_gamma );
1482 if (dt_gamma != 0.0) {
1485 sqrt( 2 * dt_gamma * kbT_bnd *
1486 ( a[i+1].
partition ? tempFactor : 1.0 ) / mass );
1487 v_bnd /= ( 1. + 0.5 * dt_gamma );
1498 if ( ! dt_gamma )
continue;
1501 sqrt( 2 * dt_gamma * kbT *
1502 ( a[i].
partition ? tempFactor : 1.0 ) * a[i].recipMass );
1503 a[i].
velocity /= ( 1. + 0.5 * dt_gamma );
1517 for ( i = 0; i < numAtoms; ++i )
1520 if ( ! dt_gamma )
continue;
1523 sqrt( 2 * dt_gamma * kbT *
1524 ( a[i].
partition ? tempFactor : 1.0 ) * a[i].recipMass );
1525 a[i].
velocity /= ( 1. + 0.5 * dt_gamma );
1549 for (
int i = 0; i < numAtoms; i += hgs ) {
1553 for ( j = i; j < (i+hgs); ++j ) {
1555 a[j].fixedPosition,a[j].transform);
1561 for ( j = i; j < (i+hgs); ++j ) {
1569 Position delta_x_cm = new_x_cm - x_cm;
1570 for ( j = i; j < (i+hgs); ++j ) {
1573 a[j].fixedPosition,a[j].transform);
1582 for (
int i = 0; i < numAtoms; ++i )
1586 a[i].fixedPosition,a[i].transform);
1618 for (
int i = 0; i < numAtoms; i += hgs ) {
1622 for ( j = i; j < (i+hgs); ++j ) {
1624 a[j].fixedPosition,a[j].transform);
1631 for ( j = i; j < (i+hgs); ++j ) {
1640 Position delta_x_cm = new_x_cm - x_cm;
1643 delta_v_cm.
x = ( velFactor.
x - 1 ) * v_cm.
x;
1644 delta_v_cm.
y = ( velFactor.
y - 1 ) * v_cm.
y;
1645 delta_v_cm.
z = ( velFactor.
z - 1 ) * v_cm.
z;
1646 for ( j = i; j < (i+hgs); ++j ) {
1649 a[j].fixedPosition,a[j].transform);
1660 for (
int i = 0; i < numAtoms; ++i )
1664 a[i].fixedPosition,a[i].transform);
1681 if ( rescaleFreq > 0 ) {
1688 for (
int i = 0; i < numAtoms; ++i )
1704 const BigReal factor_dihe = accelMDfactor[0];
1705 const BigReal factor_tot = accelMDfactor[1];
1709 NAMD_die(
"accelMD broadcasting error!\n");
1711 NAMD_die(
"accelMD broadcasting error!\n");
1714 for (
int i = 0; i < numAtoms; ++i)
1720 for (
int i = 0; i < numAtoms; ++i)
1723 for (
int i = 0; i < numAtoms; ++i)
1726 if (doFullElectrostatics) {
1727 for (
int i = 0; i < numAtoms; ++i)
1733 for (
int i = 0; i < numAtoms; ++i)
1744 if ( (step < simParams->adaptTempFirstStep ) ||
1759 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1768 if ( newTemp < simParams->reassignHold )
1776 for (
int i = 0; i < numAtoms; ++i )
1780 sqrt(kbT * (a[i].
partition ? tempFactor : 1.0) * a[i].recipMass) *
1784 NAMD_bug(
"Sequencer::reassignVelocities called improperly!");
1798 for (
int i = 0; i < numAtoms; ++i )
1801 a[i].mass <= 0.) ?
Vector(0,0,0) :
1802 sqrt(kbT * (a[i].
partition ? tempFactor : 1.0) * a[i].recipMass) *
1815 for (
int i = 0; i < numAtoms; ++i )
1826 for (
int i = 0; i < numAtoms; ++i )
1838 BigReal sqrt_factor = sqrt(factor);
1841 for (
int i = 0; i < numAtoms; ++i ) {
1859 for (
int i = 0; i < numAtoms; ++i )
1861 BigReal f1 = exp( coefficient * a[i].langevinParam );
1883 DebugM(4,
"stochastically rescaling velocities at step " << step <<
" by " << coefficient <<
"\n");
1884 for (
int i = 0; i < numAtoms; ++i )
1905 BigReal timestep,
const int ftag,
const int useSaved
1908 NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1910 CmiNetworkProgressAfter (0);
1920 const BigReal timestep1,
const int ftag1,
const int useSaved1,
1921 const BigReal timestep2,
const int ftag2,
const int useSaved2,
1922 const BigReal timestep3,
const int ftag3,
const int useSaved3
1925 NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1927 CmiNetworkProgressAfter (0);
1946 NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
1948 CmiNetworkProgressAfter (0);
1959 Tensor *vp = ( pressure ? &virial : 0 );
1961 iout <<
iERROR <<
"Constraint failure in HardWallDrude(); "
1962 <<
"simulation may become unstable.\n" <<
endi;
1975 Tensor *vp = ( pressure ? &virial : 0 );
1978 "Constraint failure; simulation has become unstable.\n" <<
endi;
1983 printf(
"virial = %g %g %g %g %g %g %g %g %g\n",
1984 virial.
xx, virial.
xy, virial.
xz,
1985 virial.
yx, virial.
yy, virial.
yz,
1986 virial.
zx, virial.
zy, virial.
zz);
1993 printf(
"pos[%d] = %g %g %g\n", n,
1997 printf(
"vel[%d] = %g %g %g\n", n,
2002 printf(
"force[%d] = %g %g %g\n", n,
2036 const BigReal maxvel2 = maxvel * maxvel;
2037 for (
int i=0; i<numAtoms; ++i ) {
2038 if ( a[i].velocity.length2() > maxvel2 ) {
2045 const BigReal maxvel2 = maxvel * maxvel;
2047 for (
int i=0; i<numAtoms; ++i ) {
2052 for (
int i=0; i<numAtoms; ++i ) {
2053 if ( a[i].velocity.length2() > maxvel2 ) {
2055 iout <<
iERROR <<
"Atom " << (a[i].
id + 1) <<
" velocity is "
2058 << i <<
" of " << numAtoms <<
" on patch "
2063 "Atoms moving too fast; simulation has become unstable ("
2065 <<
" pe " << CkMyPe() <<
").\n" <<
endi;
2077 for (
int i=0; i<numAtoms; ++i ) {
2093 CmiNetworkProgressAfter (0);
2104 for (
int i = 0; i < numAtoms; ++i ) {
2107 virial.
outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2111 for (
int i = 0; i < numAtoms; ++i ) {
2112 if (a[i].mass < 0.01)
continue;
2114 virial.
outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2119 momentumSqrSum = virial;
2121 kineticEnergy *= 0.5 * 0.5;
2145 for (
int i=0; i<numAtoms; i += hgs) {
2152 for (j=i; j< i+hgs; ++j) {
2157 for (j=i; j < i+hgs; ++j) {
2159 if (! (useGroupPressure && j != i)) {
2161 int slab = (int)floor((z-zmin)*idz);
2162 if (slab < 0) slab += nslabs;
2163 else if (slab >= nslabs) slab -= nslabs;
2164 ppoffset = 3*(slab + partition*nslabs);
2167 if (useGroupPressure) {
2190 for (
int i = 0; i < numAtoms; i += hgs ) {
2193 CmiNetworkProgress ();
2200 for ( j = i; j < (i+hgs); ++j ) {
2205 momentumSqrSum.
outerAdd(1.0/m_cm, v_cm, v_cm);
2210 for ( j = i; j < (i+hgs); ++j ) {
2215 intKineticEnergy += mass * (v * dv);
2216 intVirialNormal.
outerAdd (mass, v, dv);
2220 for ( j = i; j < (i+hgs); ++j ) {
2224 intKineticEnergy += mass * (v * dv);
2225 intVirialNormal.
outerAdd(mass, v, dv);
2230 intKineticEnergy *= 0.5 * 0.5;
2232 intVirialNormal *= 0.5;
2235 momentumSqrSum *= 0.5;
2248 for (
int j = 0; j < numAtoms; j++ ) {
2266 NamdProfileEvent::SUBMIT_REDUCTIONS);
2272 CmiNetworkProgressAfter(0);
2284 Vector angularMomentum = 0;
2289 for (i = 0; i < numAtoms; ++i ) {
2293 angularMomentum +=
cross(a[i].mass,a[i].position-o,a[i].velocity);
2297 for (i = 0; i < numAtoms; ++i ) {
2300 angularMomentum +=
cross(a[i].mass,a[i].position-o,a[i].velocity);
2306 for (i = 0; i < numAtoms; i++) {
2307 if (i < numAtoms-1 &&
2308 a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
2318 drudeComKE += m_com * v_com.
length2();
2319 drudeBondKE += m_bond * v_bond.length2();
2338 kineticEnergy *= 0.5;
2348 for (
int i = 0; i < numAtoms; ++i ) {
2355 for (
int i = 0; i < numAtoms; ++i ) {
2362 for (
int i = 0; i < numAtoms; ++i ) {
2378 for (
int i = 0; i < numAtoms; i += hgs ) {
2380 CmiNetworkProgress();
2387 for ( j = i; j < (i+hgs); ++j ) {
2397 for ( j = i; j < (i+hgs); ++j ) {
2399 ( pairInteractionSelf || a[j].
partition != 2 ) )
continue;
2401 if ( fixedAtomsOn && a[j].atomFixed )
continue;
2405 intKineticEnergy += mass * (v * dv);
2412 for ( j = i; j < (i+hgs); ++j ) {
2414 if ( fixedAtomsOn && a[j].atomFixed )
continue;
2418 intKineticEnergy += mass * (v * dv);
2427 intKineticEnergy *= 0.5;
2443 for (
int i=0; i<numAtoms; i += hgs) {
2448 for (j=i; j< i+hgs; ++j) {
2455 int slab = (int)floor((z-zmin)*idz);
2456 if (slab < 0) slab += nslabs;
2457 else if (slab >= nslabs) slab -= nslabs;
2459 int ppoffset = 3*(slab + nslabs*
partition);
2460 for (j=i; j < i+hgs; ++j) {
2466 BigReal wxx = (fnormal.
x + fnbond.
x + fslow.
x) * dx.
x;
2467 BigReal wyy = (fnormal.
y + fnbond.
y + fslow.
y) * dx.
y;
2468 BigReal wzz = (fnormal.
z + fnbond.
z + fslow.
z) * dx.
z;
2483 Vector fixForceNormal = 0;
2484 Vector fixForceNbond = 0;
2487 calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
2497 #endif // UPPER_BOUND
2514 const double drudeBondLen2 = drudeBondLen * drudeBondLen;
2516 const double drudeMove = 0.01;
2517 const double drudeStep2 = drudeStep * drudeStep;
2518 const double drudeMove2 = drudeMove * drudeMove;
2523 for (
int i = 0; i < numAtoms; ++i ) {
2524 f1[i] += f2[i] + f3[i];
2525 if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) {
2526 if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
2527 if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
2530 a[i].
position += drudeStep * f1[i];
2532 if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
2536 Vector netf = f1[i-1] + f1[i];
2537 if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
2541 if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
2548 for (
int i = 0; i < numAtoms; ++i ) {
2551 if ( v2 > maxv2 ) maxv2 = v2;
2554 if ( v2 > maxv2 ) maxv2 = v2;
2565 for (
int i = 0; i < numAtoms; ++i ) {
2567 if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) )
continue;
2576 BigReal fmult = 1.01 * sqrt(fmax2/ff);
2577 f *= fmult; ff = f * f;
2596 for (
int i = 0; i < numAtoms; i += hgs ) {
2601 for ( j = i; j < (i+hgs); ++j ) {
2606 for ( j = i; j < (i+hgs); ++j ) {
2626 Vector fixForceNormal = 0;
2627 Vector fixForceNbond = 0;
2630 calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
2653 NamdProfileEvent::SUBMIT_COLLECTIONS);
2671 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
2681 if ( !
simParams->usePairlists ) pairlists = 0;
2682 patch->flags.usePairlists = pairlists || pairlistsAreValid;
2683 patch->flags.savePairlists =
2684 pairlists && ! pairlistsAreValid;
2686 if (
simParams->singleTopology ) patch->reposition_all_alchpairs();
2687 if (
simParams->lonepairs ) patch->reposition_all_lonepairs();
2699 patch->positionsReady(migration);
2701 int seq = patch->flags.sequence;
2702 int basePriority = ( (seq & 0xffff) << 15 )
2704 if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
2707 patch->gbisComputeAfterP1();
2710 patch->gbisComputeAfterP2();
2734 if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
2735 pairlistsAreValid = 1;
2740 if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
2748 if (patch->flags.doNonbonded) {
2753 if (patch->flags.doFullElectrostatics) {
2764 if (patch->flags.doNonbonded) {
2769 if (patch->flags.doFullElectrostatics) {
2780 if (patch->flags.doNonbonded) {
2785 if (patch->flags.doFullElectrostatics) {
2794 if (patch->flags.doNonbonded) {
2797 if (patch->flags.doFullElectrostatics) {
2802 if ( patch->flags.doMolly ) {
2804 patch->mollyMollify(&virial);
2810 if (patch->flags.doLoweAndersen) {
2811 patch->loweAndersenFinish();
2814 #ifdef NAMD_CUDA_XXX
2815 int numAtoms = patch->numAtoms;
2817 for (
int i=0; i<numAtoms; ++i ) {
2818 CkPrintf(
"%d %g %g %g\n", a[i].
id,
2828 CkPrintf(
"%d %g %g %g\n", a[i].
id,
2832 CkPrintf(
"%d %g %g %g\n", a[i].
id,
2836 CkPrintf(
"%d %g %g %g\n", a[i].
id,
2844 int numAtoms = patch->numAtoms;
2846 for (
int i=0; i<numAtoms; ++i ) {
2856 float fx = fxNo+fxNb+fxSl;
2857 float fy = fyNo+fyNb+fySl;
2858 float fz = fzNo+fzNb+fzSl;
2860 float f = sqrt(fx*fx+fy*fy+fz*fz);
2861 int id = patch->pExt[i].id;
2862 int seq = patch->flags.sequence;
2863 float x = patch->p[i].position.x;
2864 float y = patch->p[i].position.y;
2865 float z = patch->p[i].position.z;
2867 CkPrintf(
"FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,
id,
2902 #ifdef MEASURE_NAMD_WITH_PAPI
2903 void Sequencer::papiMeasureBarrier(
int step){
2905 broadcast->papiMeasureBarrier.get(step);
SubmitReduction * multigratorReduction
Vector gaussian_vector(void)
void rescaleVelocities(int)
void minimizationQuenchVelocity(void)
void max(int i, BigReal v)
void tcoupleVelocities(BigReal, int)
void addMovDragToPosition(BigReal)
BigReal soluteScalingFactorCharge
virtual void algorithm(void)
#define NAMD_EVENT_STOP(eon, id)
SubmitReduction * pressureProfileReduction
void addVelocityToPosition(FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
static int velocityNeeded(int)
void scaleVelocities(const BigReal velScale)
#define GB1_COMPUTE_HOME_PRIORITY
void addVelocityToPosition(BigReal)
SubmitReduction * reduction
SubmitReduction * min_reduction
SimpleBroadcastObject< int > traceBarrier
void maximumMove(BigReal)
unsigned int hydrogenGroupSize
Bool is_atom_movdragged(int atomnum) const
static int coordinateNeeded(int)
void cycleBarrier(int, int)
static void partition(int *order, const FullAtom *atoms, int begin, int end)
SimpleBroadcastObject< Vector > momentumCorrection
void addRotDragToPosition(BigReal)
static PatchMap * Object()
void saveForce(const int ftag=Results::normal)
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
void langevinVelocitiesBBK2(BigReal)
#define ADD_TENSOR_OBJECT(R, RL, D)
void newMinimizeDirection(BigReal)
void newMinimizePosition(BigReal)
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
void addForceToMomentum(FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
void startWork(const LDObjHandle &handle)
void langevinVelocitiesBBK1(BigReal)
std::ostream & endi(std::ostream &s)
char const *const NamdProfileEventStr[]
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
int berendsenPressureFreq
SubmitReduction * willSubmit(int setID, int size=-1)
void rattle1(BigReal, int)
void saveTotalForces(HomePatch *)
if(ComputeNonbondedUtil::goMethod==2)
void submitVelocities(int seq, int zero, FullAtomList &a)
SimpleBroadcastObject< BigReal > adaptTemperature
void rebalanceLoad(int timestep)
static ReductionMgr * Object(void)
void minimizeMoveDownhill(BigReal fmax2)
void addForceToMomentum(BigReal, const int ftag=Results::normal, const int useSaved=0)
BigReal length(void) const
void pauseWork(const LDObjHandle &handle)
SimpleBroadcastObject< BigReal > tcoupleCoefficient
static int forceNeeded(int)
void exchangeCheckpoint(int scriptTask, int &bpc)
void split(int iStream, int numStreams)
void addForceToMomentum3(const BigReal timestep1, const int ftag1, const int useSaved1, const BigReal timestep2, const int ftag2, const int useSaved2, const BigReal timestep3, const int ftag3, const int useSaved3)
SimpleBroadcastObject< BigReal > stochRescaleCoefficient
void submitCollections(int step, int zeroVel=0)
BigReal calcKineticEnergy()
void adaptTempUpdate(int)
void rescale(Tensor factor)
#define TIMER_START(T, TYPE)
void calcFixVirial(Tensor &fixVirialNormal, Tensor &fixVirialNbond, Tensor &fixVirialSlow, Vector &fixForceNormal, Vector &fixForceNbond, Vector &fixForceSlow)
Bool is_atom_rotdragged(int atomnum) const
#define NAMD_PROFILE_START()
#define NAMD_EVENT_START(eon, id)
#define COMPUTE_HOME_PRIORITY
void NAMD_bug(const char *err_msg)
void multigratorPressure(int step, int callNumber)
void berendsenPressure(int)
void submitMomentum(int step)
int rescaleVelocities_numTemps
void submitLoadStats(int timestep)
void runComputeObjects(int migration=1, int pairlists=0, int pressureStep=0)
void rebalance(Sequencer *seq, PatchID id)
void rescaleaccelMD(int, int, int)
SimpleBroadcastObject< Tensor > velocityRescaleTensor2
unsigned char get_ss_type(int anum) const
#define NAMD_EVENT_RANGE_2(eon, id)
SimpleBroadcastObject< int > scriptBarrier
void scalePositionsVelocities(const Tensor &posScale, const Tensor &velScale)
SimpleBroadcastObject< BigReal > velocityRescaleFactor2
void NAMD_die(const char *err_msg)
static LdbCoordinator * Object()
static Tensor identity(BigReal v1=1.0)
#define TIMER_INIT_WIDTH(T, TYPE, WIDTH)
void submitPositions(int seq, FullAtomList &a, Lattice l, int prec)
int berendsenPressure_count
SimpleBroadcastObject< BigReal > velocityRescaleFactor
SimpleBroadcastObject< BigReal > minimizeCoefficient
void reassignVelocities(BigReal, int)
SimpleBroadcastObject< Vector > accelMDRescaleFactor
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
ComputeGlobal * computeGlobalObject
void saveForce(const int ftag=Results::normal)
SimpleBroadcastObject< Tensor > positionRescaleFactor
void langevinVelocities(BigReal)
void submitForces(int seq, FullAtomList &a, int maxForceUsed, ForceList *f)
void hardWallDrude(BigReal, int)
#define TIMER_STOP(T, TYPE)
void multigratorTemperature(int step, int callNumber)
BigReal length2(void) const
void reinitVelocities(void)
int pressureProfileAtomTypes
int checkpoint_berendsenPressure_count
int numPatches(void) const
ControllerBroadcasts * broadcast
#define NAMD_EVENT_START_EX(eon, id, str)
Position apply_transform(Position data, const Transform &t) const
CollectionMgr *const collection
void stochRescaleVelocities(BigReal, int)
#define GB2_COMPUTE_HOME_PRIORITY
void get_movdrag_params(Vector &v, int atomnum) const
#define NAMD_PROFILE_STOP()
void get_rotdrag_params(BigReal &v, Vector &a, Vector &p, int atomnum) const
int getNumStepsToRun(void)
void newtonianVelocities(BigReal, const BigReal, const BigReal, const BigReal, const int, const int, const int)
void rescaleSoluteCharges(BigReal)
Real atomcharge(int anum) const
void submitMinimizeReductions(int, BigReal fmax2)
#define ADD_VECTOR_OBJECT(R, RL, D)
const_iterator const_begin(void) const
int multigratorPressureFreq
void correctMomentum(int step, BigReal drifttime)
Bool is_atom_exPressure(int atomnum) const
std::ostream & iERROR(std::ostream &s)
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)
SubmitReduction * reduction
Used to submit restraint energy as MISC.
ForceList f[Results::maxNumForces]
void enableEarlyExit(void)
void submitReductions(int)
SimpleBroadcastObject< Tensor > positionRescaleFactor2
SimParameters *const simParams
SimpleBroadcastObject< Tensor > velocityRescaleTensor
void rescaleVelocitiesByFactor(BigReal)
int multigratorTemperatureFreq
#define PATCH_PRIORITY(PID)
void exchangeAtoms(int scriptTask)