#include <Controller.h>
|
|
|
Definition at line 203 of file Controller.C. 00204 {
00205 delete broadcast;
00206 delete reduction;
00207 delete amd_reduction;
00208 delete ppbonded;
00209 delete ppnonbonded;
00210 delete ppint;
00211 delete [] pressureProfileAverage;
00212 delete random;
00213 }
|
|
|
Definition at line 1424 of file Controller.C. References SimParameters::adaptTempAutoDt, adaptTempAutoDt, adaptTempBetaMax, adaptTempBetaMin, adaptTempBetaN, adaptTempBins, SimParameters::adaptTempBins, adaptTempCg, SimParameters::adaptTempCgamma, adaptTempDBeta, SimParameters::adaptTempDt, adaptTempDt, adaptTempDTave, adaptTempDTavenum, adaptTempDtMax, adaptTempDtMin, SimParameters::adaptTempInFile, SimParameters::adaptTempOn, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, adaptTempRestartFile, SimParameters::adaptTempRestartFile, SimParameters::adaptTempRestartFreq, adaptTempT, SimParameters::adaptTempTmax, SimParameters::adaptTempTmin, BigReal, iINFO(), SimParameters::initialTemp, iout, j, SimParameters::langevinOn, SimParameters::langevinTemp, NAMD_backup_file(), NAMD_die(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, and simParams. Referenced by adaptTempUpdate(). 01424 {
01425 if (!simParams->adaptTempOn) return;
01426 iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
01427 adaptTempDtMin = 0;
01428 adaptTempDtMax = 0;
01429 adaptTempAutoDt = false;
01430 if (simParams->adaptTempBins == 0) {
01431 iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
01432 std::ifstream adaptTempRead(simParams->adaptTempInFile);
01433 if (adaptTempRead) {
01434 int readInt;
01435 BigReal readReal;
01436 bool readBool;
01437 // step
01438 adaptTempRead >> readInt;
01439 // Start with min and max temperatures
01440 adaptTempRead >> adaptTempT; // KELVIN
01441 adaptTempRead >> adaptTempBetaMin; // KELVIN
01442 adaptTempRead >> adaptTempBetaMax; // KELVIN
01443 adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
01444 adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
01445 // In case file is manually edited
01446 if (adaptTempBetaMin > adaptTempBetaMax){
01447 readReal = adaptTempBetaMax;
01448 adaptTempBetaMax = adaptTempBetaMin;
01449 adaptTempBetaMin = adaptTempBetaMax;
01450 }
01451 adaptTempRead >> adaptTempBins;
01452 adaptTempRead >> adaptTempCg;
01453 adaptTempRead >> adaptTempDt;
01454 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
01455 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
01456 adaptTempPotEnergySamples = new int[adaptTempBins];
01457 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
01458 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
01459 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
01460 adaptTempBetaN = new BigReal[adaptTempBins];
01461 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
01462 for(int j = 0; j < adaptTempBins; ++j) {
01463 adaptTempRead >> adaptTempPotEnergyAve[j];
01464 adaptTempRead >> adaptTempPotEnergyVar[j];
01465 adaptTempRead >> adaptTempPotEnergySamples[j];
01466 adaptTempRead >> adaptTempPotEnergyAveNum[j];
01467 adaptTempRead >> adaptTempPotEnergyVarNum[j];
01468 adaptTempRead >> adaptTempPotEnergyAveDen[j];
01469 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
01470 }
01471 adaptTempRead.close();
01472 }
01473 else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
01474 }
01475 else {
01476 adaptTempBins = simParams->adaptTempBins;
01477 adaptTempPotEnergyAveNum = new BigReal[adaptTempBins];
01478 adaptTempPotEnergyAveDen = new BigReal[adaptTempBins];
01479 adaptTempPotEnergySamples = new int[adaptTempBins];
01480 adaptTempPotEnergyVarNum = new BigReal[adaptTempBins];
01481 adaptTempPotEnergyVar = new BigReal[adaptTempBins];
01482 adaptTempPotEnergyAve = new BigReal[adaptTempBins];
01483 adaptTempBetaN = new BigReal[adaptTempBins];
01484 adaptTempBetaMax = 1./simParams->adaptTempTmin;
01485 adaptTempBetaMin = 1./simParams->adaptTempTmax;
01486 adaptTempCg = simParams->adaptTempCgamma;
01487 adaptTempDt = simParams->adaptTempDt;
01488 adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
01489 adaptTempT = simParams->initialTemp;
01490 if (simParams->langevinOn)
01491 adaptTempT = simParams->langevinTemp;
01492 else if (simParams->rescaleFreq > 0)
01493 adaptTempT = simParams->rescaleTemp;
01494 for(int j = 0; j < adaptTempBins; ++j){
01495 adaptTempPotEnergyAveNum[j] = 0.;
01496 adaptTempPotEnergyAveDen[j] = 0.;
01497 adaptTempPotEnergySamples[j] = 0;
01498 adaptTempPotEnergyVarNum[j] = 0.;
01499 adaptTempPotEnergyVar[j] = 0.;
01500 adaptTempPotEnergyAve[j] = 0.;
01501 adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
01502 }
01503 }
01504 if (simParams->adaptTempAutoDt > 0.0) {
01505 adaptTempAutoDt = true;
01506 adaptTempDtMin = simParams->adaptTempAutoDt - 0.01;
01507 adaptTempDtMax = simParams->adaptTempAutoDt + 0.01;
01508 }
01509 adaptTempDTave = 0;
01510 adaptTempDTavenum = 0;
01511 iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
01512 iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
01513 iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
01514 iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
01515 iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
01516 if (simParams->adaptTempRestartFreq > 0) {
01517 NAMD_backup_file(simParams->adaptTempRestartFile);
01518 adaptTempRestartFile.open(simParams->adaptTempRestartFile);
01519 if(!adaptTempRestartFile)
01520 NAMD_die("Error opening restart file");
01521 }
01522 }
|
|
||||||||||||
|
Definition at line 1550 of file Controller.C. References adaptTempBetaMax, adaptTempBetaMin, adaptTempBetaN, adaptTempBin, adaptTempBins, adaptTempCg, adaptTempDBeta, SimParameters::adaptTempDebug, adaptTempDt, adaptTempDTave, adaptTempDTavenum, adaptTempDtMax, ControllerBroadcasts::adaptTemperature, SimParameters::adaptTempFirstStep, SimParameters::adaptTempFreq, adaptTempInit(), SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, SimParameters::adaptTempOutFreq, adaptTempPotEnergyAve, adaptTempPotEnergyAveDen, adaptTempPotEnergyAveNum, adaptTempPotEnergySamples, adaptTempPotEnergyVar, adaptTempPotEnergyVarNum, SimParameters::adaptTempRandom, adaptTempT, adaptTempWriteRestart(), BigReal, BOLTZMANN, broadcast, SimParameters::firstTimestep, Random::gaussian(), iout, iWARN(), j, SimpleBroadcastObject< T >::publish(), random, simParams, totalEnergy, and Random::uniform(). Referenced by integrate(). 01551 {
01552 //Beta = 1./T
01553 if ( !simParams->adaptTempOn ) return;
01554 int j = 0;
01555 if (step == simParams->firstTimestep) {
01556 adaptTempInit(step);
01557 return;
01558 }
01559 if ( minimize || (step < simParams->adaptTempFirstStep ) ||
01560 ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
01561 const int adaptTempOutFreq = simParams->adaptTempOutFreq;
01562 const bool adaptTempDebug = simParams->adaptTempDebug;
01563 //Calculate Current inverse temperature and bin
01564 BigReal adaptTempBeta = 1./adaptTempT;
01565 adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
01566
01567 if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
01568 iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
01569 << " adaptTempBeta: " << adaptTempBeta
01570 << " adaptTempDBeta: " << adaptTempDBeta
01571 << " betaMin:" << adaptTempBetaMin
01572 << " betaMax: " << adaptTempBetaMax << "\n";
01573 adaptTempPotEnergySamples[adaptTempBin] += 1;
01574 BigReal gammaAve = 1.-adaptTempCg/adaptTempPotEnergySamples[adaptTempBin];
01575
01576 BigReal potentialEnergy;
01577 BigReal potEnergyAverage;
01578 BigReal potEnergyVariance;
01579 potentialEnergy = totalEnergy-kineticEnergy;
01580
01581 //calculate new bin average and variance using adaptive averaging
01582 adaptTempPotEnergyAveNum[adaptTempBin] = adaptTempPotEnergyAveNum[adaptTempBin]*gammaAve + potentialEnergy;
01583 adaptTempPotEnergyAveDen[adaptTempBin] = adaptTempPotEnergyAveDen[adaptTempBin]*gammaAve + 1;
01584 adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
01585
01586 potEnergyAverage = adaptTempPotEnergyAveNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
01587 potEnergyVariance = adaptTempPotEnergyVarNum[adaptTempBin]/adaptTempPotEnergyAveDen[adaptTempBin];
01588 potEnergyVariance -= potEnergyAverage*potEnergyAverage;
01589
01590 adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
01591 adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
01592
01593 // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
01594 // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
01595 // is constant for each bin. This is to estimate <E(beta)> where beta \in
01596 // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
01597 if ( ! ( step % simParams->adaptTempFreq ) ) {
01598 // If adaptTempBin not at the edge, calculate improved average:
01599 if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
01600 // Get Averaging Limits:
01601 BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
01602 BigReal betaPlus;
01603 BigReal betaMinus;
01604 int nMinus =0;
01605 int nPlus = 0;
01606 if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
01607 if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
01608 betaMinus = adaptTempBeta - deltaBeta;
01609 betaPlus = adaptTempBeta + deltaBeta;
01610 BigReal betaMinus2 = betaMinus*betaMinus;
01611 BigReal betaPlus2 = betaPlus*betaPlus;
01612 nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
01613 nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
01614 // Variables for <E(beta)> estimate:
01615 BigReal potEnergyAve0 = 0.0;
01616 BigReal potEnergyAve1 = 0.0;
01617 // Integral terms
01618 BigReal A0 = 0;
01619 BigReal A1 = 0;
01620 BigReal A2 = 0;
01621 //A0 phi_s integral for beta_minus < beta < beta_{i+1}
01622 BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
01623 BigReal DeltaE2Ave;
01624 BigReal DeltaE2AveSum = 0;
01625 BigReal betaj;
01626 for (j = nMinus+1; j <= adaptTempBin; ++j) {
01627 potEnergyAve0 += adaptTempPotEnergyAve[j];
01628 DeltaE2Ave = adaptTempPotEnergyVar[j];
01629 DeltaE2AveSum += DeltaE2Ave;
01630 A0 += j*DeltaE2Ave;
01631 }
01632 A0 *= adaptTempDBeta;
01633 A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
01634 A0 *= adaptTempDBeta;
01635 betaj = adaptTempBetaN[nMinus+1]-betaMinus;
01636 betaj *= betaj;
01637 A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
01638 A0 /= (betaNp1-betaMinus);
01639
01640 //A1 phi_s integral for beta_{i+1} < beta < beta_plus
01641 DeltaE2AveSum = 0;
01642 for (j = adaptTempBin+1; j < nPlus; j++) {
01643 potEnergyAve1 += adaptTempPotEnergyAve[j];
01644 DeltaE2Ave = adaptTempPotEnergyVar[j];
01645 DeltaE2AveSum += DeltaE2Ave;
01646 A1 += j*DeltaE2Ave;
01647 }
01648 A1 *= adaptTempDBeta;
01649 A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
01650 A1 *= adaptTempDBeta;
01651 if ( nPlus < adaptTempBins ) {
01652 betaj = betaPlus-adaptTempBetaN[nPlus];
01653 betaj *= betaj;
01654 A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
01655 }
01656 A1 /= (betaPlus-betaNp1);
01657
01658 //A2 phi_t integral for beta_i
01659 A2 = 0.5*adaptTempDBeta*potEnergyVariance;
01660
01661 // Now calculate a+ and a-
01662 DeltaE2Ave = A0-A1;
01663 BigReal aplus = 0;
01664 BigReal aminus = 0;
01665 if (DeltaE2Ave != 0) {
01666 aplus = (A0-A2)/(A0-A1);
01667 if (aplus < 0) {
01668 aplus = 0;
01669 }
01670 if (aplus > 1) {
01671 aplus = 1;
01672 }
01673 aminus = 1-aplus;
01674 potEnergyAve0 *= adaptTempDBeta;
01675 potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
01676 potEnergyAve0 /= (betaNp1-betaMinus);
01677 potEnergyAve1 *= adaptTempDBeta;
01678 if ( nPlus < adaptTempBins ) {
01679 potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
01680 }
01681 potEnergyAve1 /= (betaPlus-betaNp1);
01682 potEnergyAverage = aminus*potEnergyAve0;
01683 potEnergyAverage += aplus*potEnergyAve1;
01684 }
01685 if (simParams->adaptTempDebug) {
01686 iout << "ADAPTEMP DEBUG:" << "\n"
01687 << " adaptTempBin: " << adaptTempBin << "\n"
01688 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
01689 << " adaptTempBeta: " << adaptTempBeta << "\n"
01690 << " adaptTemp: " << adaptTempT<< "\n"
01691 << " betaMin: " << adaptTempBetaMin << "\n"
01692 << " betaMax: " << adaptTempBetaMax << "\n"
01693 << " gammaAve: " << gammaAve << "\n"
01694 << " deltaBeta: " << deltaBeta << "\n"
01695 << " betaMinus: " << betaMinus << "\n"
01696 << " betaPlus: " << betaPlus << "\n"
01697 << " nMinus: " << nMinus << "\n"
01698 << " nPlus: " << nPlus << "\n"
01699 << " A0: " << A0 << "\n"
01700 << " A1: " << A1 << "\n"
01701 << " A2: " << A2 << "\n"
01702 << " a+: " << aplus << "\n"
01703 << " a-: " << aminus << "\n"
01704 << endi;
01705 }
01706 }
01707 else {
01708 if (simParams->adaptTempDebug) {
01709 iout << "ADAPTEMP DEBUG:" << "\n"
01710 << " adaptTempBin: " << adaptTempBin << "\n"
01711 << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
01712 << " adaptTempBeta: " << adaptTempBeta << "\n"
01713 << " adaptTemp: " << adaptTempT<< "\n"
01714 << " betaMin: " << adaptTempBetaMin << "\n"
01715 << " betaMax: " << adaptTempBetaMax << "\n"
01716 << " gammaAve: " << gammaAve << "\n"
01717 << " deltaBeta: N/A\n"
01718 << " betaMinus: N/A\n"
01719 << " betaPlus: N/A\n"
01720 << " nMinus: N/A\n"
01721 << " nPlus: N/A\n"
01722 << " A0: N/A\n"
01723 << " A1: N/A\n"
01724 << " A2: N/A\n"
01725 << " a+: N/A\n"
01726 << " a-: N/A\n"
01727 << endi;
01728 }
01729 }
01730
01731 //dT is new temperature
01732 BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
01733 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
01734 dT += adaptTempT;
01735
01736 // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
01737 // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
01738 if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
01739 dT = ((potentialEnergy-adaptTempPotEnergyAve[adaptTempBin])/BOLTZMANN+adaptTempT)*adaptTempDt;
01740 dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
01741 dT += adaptTempT;
01742 // Check again, if not then keep original adaptTempTor assign random.
01743 if ( dT > 1./adaptTempBetaMin ) {
01744 if (!simParams->adaptTempRandom) {
01745 //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
01746 // << " K higher than adaptTempTmax."
01747 // << " Keeping temperature at "
01748 // << adaptTempT<< "\n"<< endi;
01749 dT = adaptTempT;
01750 }
01751 else {
01752 //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
01753 // << " K higher than adaptTempTmax."
01754 // << " Assigning random temperature in range\n" << endi;
01755 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
01756 dT = 1./dT;
01757 }
01758 }
01759 else if ( dT < 1./adaptTempBetaMax ) {
01760 if (!simParams->adaptTempRandom) {
01761 //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
01762 // << " K lower than adaptTempTmin."
01763 // << " Keeping temperature at " << adaptTempT<< "\n" << endi;
01764 dT = adaptTempT;
01765 }
01766 else {
01767 //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
01768 // << " K lower than adaptTempTmin."
01769 // << " Assigning random temperature in range\n" << endi;
01770 dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
01771 dT = 1./dT;
01772 }
01773 }
01774 else if (adaptTempAutoDt) {
01775 //update temperature step size counter
01776 //FOR "TRUE" ADAPTIVE TEMPERING
01777 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
01778 if (adaptTempTdiff > 0) {
01779 adaptTempDTave += adaptTempTdiff;
01780 adaptTempDTavenum++;
01781 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
01782 }
01783 if(adaptTempDTavenum == 100){
01784 BigReal Frac;
01785 adaptTempDTave /= adaptTempDTavenum;
01786 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
01787 Frac = adaptTempDTave/Frac;
01788 //if average temperature jump is > 3% of temperature range,
01789 //modify jump size to match 3%
01790 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
01791 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
01792 Frac = adaptTempDtMax/Frac;
01793 iout << "ADAPTEMP: Updating adaptTempDt to ";
01794 adaptTempDt *= Frac;
01795 iout << adaptTempDt << "\n" << endi;
01796 }
01797 adaptTempDTave = 0;
01798 adaptTempDTavenum = 0;
01799 }
01800 }
01801 }
01802 else if (adaptTempAutoDt) {
01803 //update temperature step size counter
01804 // FOR "TRUE" ADAPTIVE TEMPERING
01805 BigReal adaptTempTdiff = fabs(dT-adaptTempT);
01806 if (adaptTempTdiff > 0) {
01807 adaptTempDTave += adaptTempTdiff;
01808 adaptTempDTavenum++;
01809 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
01810 }
01811 if(adaptTempDTavenum == 100){
01812 BigReal Frac;
01813 adaptTempDTave /= adaptTempDTavenum;
01814 Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
01815 Frac = adaptTempDTave/Frac;
01816 //if average temperature jump is > 3% of temperature range,
01817 //modify jump size to match 3%
01818 iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
01819 if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
01820 Frac = adaptTempDtMax/Frac;
01821 iout << "ADAPTEMP: Updating adaptTempDt to ";
01822 adaptTempDt *= Frac;
01823 iout << adaptTempDt << "\n" << endi;
01824 }
01825 adaptTempDTave = 0;
01826 adaptTempDTavenum = 0;
01827
01828 }
01829
01830 }
01831
01832 adaptTempT = dT;
01833 broadcast->adaptTemperature.publish(step,adaptTempT);
01834 }
01835 adaptTempWriteRestart(step);
01836 if ( ! (step % adaptTempOutFreq) ) {
01837 iout << "ADAPTEMP: STEP " << step
01838 << " TEMP " << adaptTempT
01839 << " ENERGY " << std::setprecision(10) << potentialEnergy
01840 << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
01841 << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
01842 iout << "\n" << endi;
01843 }
01844
01845 }
|
|
|
|
Definition at line 233 of file Controller.C. References BackEnd::awaken(), berendsenPressure_avg, berendsenPressure_count, broadcast, checkpoint_berendsenPressure_avg, checkpoint_berendsenPressure_count, checkpoint_langevinPiston_strainRate, checkpoint_lattice, checkpoint_smooth2_avg, checkpoint_stored, END_OF_RUN, enqueueCollections(), EVAL_MEASURE, FILE_OUTPUT, SimParameters::firstTimestep, FORCE_OUTPUT, SimpleBroadcastObject< T >::get(), SimParameters::initialTemp, integrate(), iout, langevinPiston_strainRate, NamdState::lattice, minimize(), NAMD_die(), outputExtendedSystem(), SCRIPT_CHECKPOINT, SCRIPT_FORCEOUTPUT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, simParams, smooth2_avg, state, and terminate(). 00234 {
00235 int scriptTask;
00236 int scriptSeq = 0;
00237 BackEnd::awaken();
00238 while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
00239 switch ( scriptTask ) {
00240 case SCRIPT_OUTPUT:
00241 enqueueCollections(FILE_OUTPUT);
00242 outputExtendedSystem(FILE_OUTPUT);
00243 break;
00244 case SCRIPT_FORCEOUTPUT:
00245 enqueueCollections(FORCE_OUTPUT);
00246 break;
00247 case SCRIPT_MEASURE:
00248 enqueueCollections(EVAL_MEASURE);
00249 break;
00250 case SCRIPT_REINITVELS:
00251 iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
00252 << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
00253 break;
00254 case SCRIPT_RESCALEVELS:
00255 iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
00256 << " BY " << simParams->scriptArg1 << "\n" << endi;
00257 break;
00258 case SCRIPT_CHECKPOINT:
00259 iout << "CHECKPOINTING POSITIONS AT STEP " << simParams->firstTimestep
00260 << "\n" << endi;
00261 checkpoint_stored = 1;
00262 checkpoint_lattice = state->lattice;
00263 checkpoint_langevinPiston_strainRate = langevinPiston_strainRate;
00264 checkpoint_berendsenPressure_avg = berendsenPressure_avg;
00265 checkpoint_berendsenPressure_count = berendsenPressure_count;
00266 checkpoint_smooth2_avg = smooth2_avg;
00267 break;
00268 case SCRIPT_REVERT:
00269 iout << "REVERTING POSITIONS AT STEP " << simParams->firstTimestep
00270 << "\n" << endi;
00271 if ( ! checkpoint_stored )
00272 NAMD_die("Unable to revert, checkpoint was never called!");
00273 state->lattice = checkpoint_lattice;
00274 langevinPiston_strainRate = checkpoint_langevinPiston_strainRate;
00275 berendsenPressure_avg = checkpoint_berendsenPressure_avg;
00276 berendsenPressure_count = checkpoint_berendsenPressure_count;
00277 smooth2_avg = checkpoint_smooth2_avg;
00278 break;
00279 case SCRIPT_MINIMIZE:
00280 minimize();
00281 break;
00282 case SCRIPT_RUN:
00283 integrate();
00284 break;
00285 }
00286 BackEnd::awaken();
00287 }
00288 enqueueCollections(END_OF_RUN);
00289 outputExtendedSystem(END_OF_RUN);
00290 // note: this is a Converse interface call that gets translated to a
00291 // null method call if CMK_OPTIMIZE is set.
00292 if(traceAvailable())
00293 traceClose();
00294 terminate();
00295 }
|
|
|
Definition at line 35 of file Controller.h. Referenced by LdbCoordinator::awakenSequencers(), resumeAfterTraceBarrier(), and run(). 00035 { CthAwaken(thread); };
|
|
|
Definition at line 647 of file Controller.C. References berendsenPressure_avg, berendsenPressure_count, SimParameters::berendsenPressureCompressibility, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, SimParameters::berendsenPressureRelaxationTime, SimParameters::berendsenPressureTarget, broadcast, Tensor::diagonal(), SimParameters::dt, Tensor::identity(), iERROR(), iout, NamdState::lattice, LIMIT_SCALING, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), Lattice::rescale(), simParams, state, Tensor::xx, Tensor::yy, and Tensor::zz. Referenced by integrate(). 00648 {
00649 if ( simParams->berendsenPressureOn ) {
00650 berendsenPressure_count += 1;
00651 berendsenPressure_avg += controlPressure;
00652 const int freq = simParams->berendsenPressureFreq;
00653 if ( ! (berendsenPressure_count % freq) ) {
00654 Tensor factor = berendsenPressure_avg / berendsenPressure_count;
00655 berendsenPressure_avg = 0;
00656 berendsenPressure_count = 0;
00657 // We only use on-diagonal terms (for now)
00658 factor = Tensor::diagonal(diagonal(factor));
00659 factor -= Tensor::identity(simParams->berendsenPressureTarget);
00660 factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
00661 simParams->dt * freq / simParams->berendsenPressureRelaxationTime );
00662 factor += Tensor::identity(1.0);
00663 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
00664 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
00665 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
00666 int limited = 0;
00667 LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
00668 LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
00669 LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
00670 #undef LIMIT_SCALING
00671 if ( limited ) {
00672 iout << iERROR << "Step " << step <<
00673 " cell rescaling factor limited.\n" << endi;
00674 }
00675 broadcast->positionRescaleFactor.publish(step,factor);
00676 state->lattice.rescale(factor);
00677 }
00678 } else {
00679 berendsenPressure_avg = 0;
00680 berendsenPressure_count = 0;
00681 }
00682 }
|
|
||||||||||||
|
Definition at line 1848 of file Controller.C. References BigReal, computeChecksum, SimParameters::fullElectFrequency, iERROR(), iINFO(), iout, RequireReduction::item(), iWARN(), marginViolations, Node::molecule, SimParameters::mollyOn, NAMD_bug(), Molecule::numAtoms, Molecule::numCalcAngles, Molecule::numCalcAnisos, Molecule::numCalcBonds, Molecule::numCalcCrossterms, Molecule::numCalcDihedrals, Molecule::numCalcExclusions, Molecule::numCalcImpropers, Molecule::numCalcTholes, Node::Object(), SimParameters::outputPairlists, pairlistWarnings, reduction, REDUCTION_ANGLE_CHECKSUM, REDUCTION_ANISO_CHECKSUM, REDUCTION_ATOM_CHECKSUM, REDUCTION_BOND_CHECKSUM, REDUCTION_COMPUTE_CHECKSUM, REDUCTION_CROSSTERM_CHECKSUM, REDUCTION_DIHEDRAL_CHECKSUM, REDUCTION_EXCLUSION_CHECKSUM, REDUCTION_IMPROPER_CHECKSUM, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_PAIRLIST_WARNINGS, REDUCTION_STRAY_CHARGE_ERRORS, REDUCTION_THOLE_CHECKSUM, and simParams. Referenced by printDynamicsEnergies(), and printMinimizeEnergies(). 01848 {
01849 Node *node = Node::Object();
01850 Molecule *molecule = node->molecule;
01851
01852 // Some basic correctness checking
01853 BigReal checksum, checksum_b;
01854
01855 checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
01856 if ( ((int)checksum) != molecule->numAtoms ) {
01857 char errmsg[256];
01858 sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
01859 (int)checksum, molecule->numAtoms);
01860 NAMD_bug(errmsg);
01861 }
01862
01863 checksum = reduction->item(REDUCTION_COMPUTE_CHECKSUM);
01864 if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
01865 if ( ! computeChecksum ) {
01866 computesPartitioned = 0;
01867 } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
01868 computesPartitioned = 1;
01869 } else {
01870 NAMD_bug("Bad global compute count!\n");
01871 }
01872 computeChecksum = ((int)checksum);
01873 }
01874
01875 checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
01876 if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
01877 if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
01878 iout << iWARN << "Bad global bond count!\n" << endi;
01879 else NAMD_bug("Bad global bond count!\n");
01880 }
01881
01882 checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM);
01883 if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
01884 if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
01885 iout << iWARN << "Bad global angle count!\n" << endi;
01886 else NAMD_bug("Bad global angle count!\n");
01887 }
01888
01889 checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
01890 if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
01891 if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
01892 iout << iWARN << "Bad global dihedral count!\n" << endi;
01893 else NAMD_bug("Bad global dihedral count!\n");
01894 }
01895
01896 checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM);
01897 if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
01898 if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
01899 iout << iWARN << "Bad global improper count!\n" << endi;
01900 else NAMD_bug("Bad global improper count!\n");
01901 }
01902
01903 checksum = reduction->item(REDUCTION_THOLE_CHECKSUM);
01904 if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
01905 if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
01906 iout << iWARN << "Bad global Thole count!\n" << endi;
01907 else NAMD_bug("Bad global Thole count!\n");
01908 }
01909
01910 checksum = reduction->item(REDUCTION_ANISO_CHECKSUM);
01911 if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
01912 if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
01913 iout << iWARN << "Bad global Aniso count!\n" << endi;
01914 else NAMD_bug("Bad global Aniso count!\n");
01915 }
01916
01917 checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM);
01918 if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
01919 if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
01920 iout << iWARN << "Bad global crossterm count!\n" << endi;
01921 else NAMD_bug("Bad global crossterm count!\n");
01922 }
01923
01924 checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM);
01925 if ( ((int)checksum) > molecule->numCalcExclusions &&
01926 ( ! simParams->mollyOn || step % slowFreq ) ) {
01927 if ( forgiving )
01928 iout << iWARN << "High global exclusion count ("
01929 << ((int)checksum) << " vs "
01930 << molecule->numCalcExclusions << "), possible error!\n"
01931 << iWARN << "This warning is not unusual during minimization.\n"
01932 << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
01933 else {
01934 char errmsg[256];
01935 sprintf(errmsg, "High global exclusion count! (%d vs %d) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
01936 (int)checksum, molecule->numCalcExclusions);
01937 NAMD_bug(errmsg);
01938 }
01939 }
01940 if ( ((int)checksum) && ((int)checksum) < molecule->numCalcExclusions ) {
01941 if ( forgiving || ! simParams->fullElectFrequency ) {
01942 iout << iWARN << "Low global exclusion count! ("
01943 << ((int)checksum) << " vs " << molecule->numCalcExclusions << ")";
01944 if ( forgiving ) {
01945 iout << "\n"
01946 << iWARN << "This warning is not unusual during minimization.\n"
01947 << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
01948 } else {
01949 iout << " System unstable or pairlistdist or cutoff too small.\n";
01950 }
01951 iout << endi;
01952 } else {
01953 char errmsg[256];
01954 sprintf(errmsg, "Low global exclusion count! (%d vs %d) System unstable or pairlistdist or cutoff too small.\n",
01955 (int)checksum, molecule->numCalcExclusions);
01956 NAMD_bug(errmsg);
01957 }
01958 }
01959
01960 checksum = reduction->item(REDUCTION_MARGIN_VIOLATIONS);
01961 if ( ((int)checksum) && ! marginViolations ) {
01962 iout << iERROR << "Margin is too small for " << ((int)checksum) <<
01963 " atoms during timestep " << step << ".\n" << iERROR <<
01964 "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
01965 }
01966 marginViolations += (int)checksum;
01967
01968 checksum = reduction->item(REDUCTION_PAIRLIST_WARNINGS);
01969 if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
01970 iout << iINFO <<
01971 "Pairlistdist is too small for " << ((int)checksum) <<
01972 " computes during timestep " << step << ".\n" << endi;
01973 }
01974 if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
01975
01976 checksum = reduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
01977 if ( checksum ) {
01978 if ( forgiving )
01979 iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
01980 else NAMD_bug("Stray PME grid charges detected!\n");
01981 }
01982 }
|
|
|
Definition at line 909 of file Controller.C. References BigReal, broadcast, iERROR(), iout, RequireReduction::item(), Vector::length2(), ControllerBroadcasts::momentumCorrection, NAMD_die(), SimpleBroadcastObject< T >::publish(), reduction, REDUCTION_MOMENTUM_MASS, slowFreq, Vector::x, Vector::y, and Vector::z. Referenced by integrate(). 00909 {
00910
00911 Vector momentum;
00912 momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
00913 momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
00914 momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
00915 const BigReal mass = reduction->item(REDUCTION_MOMENTUM_MASS);
00916
00917 if ( momentum.length2() == 0. )
00918 iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
00919 if ( mass == 0. )
00920 NAMD_die("Total mass is zero in Controller::correctMomentum");
00921
00922 momentum *= (-1./mass);
00923
00924 broadcast->momentumCorrection.publish(step+slowFreq,momentum);
00925 }
|
|
||||||||||||
|
Definition at line 2854 of file Controller.C. References broadcast. Referenced by integrate(). 02854 {
02855 #if USE_BARRIER
02856 if (doBarrier) {
02857 broadcast->cycleBarrier.publish(step,1);
02858 CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
02859 CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02860 }
02861 #endif
02862 }
|
|
|
Definition at line 2519 of file Controller.C. References collection, Output::coordinateNeeded(), CollectionMaster::enqueueForces(), CollectionMaster::enqueuePositions(), CollectionMaster::enqueueVelocities(), Output::forceNeeded(), NamdState::lattice, state, and Output::velocityNeeded(). Referenced by algorithm(), integrate(), and minimize(). 02520 {
02521 if ( Output::coordinateNeeded(timestep) ){
02522 collection->enqueuePositions(timestep,state->lattice);
02523 }
02524 if ( Output::velocityNeeded(timestep) )
02525 collection->enqueueVelocities(timestep);
02526 if ( Output::forceNeeded(timestep) )
02527 collection->enqueueForces(timestep);
02528 }
|
|
|
Definition at line 313 of file Controller.C. References adaptTempUpdate(), berendsenPressure(), correctMomentum(), cycleBarrier(), enqueueCollections(), eventEndOfTimeStep, SimParameters::firstTimestep, SimParameters::fullElectFrequency, langevinPiston1(), langevinPiston2(), SimParameters::N, nbondFreq, SimParameters::nonbondedFrequency, SimParameters::numTraceSteps, Node::Object(), outputExtendedSystem(), outputFepEnergy(), outputTiEnergy(), printDynamicsEnergies(), printFepMessage(), printTiMessage(), reassignVelocities(), rebalanceLoad(), receivePressure(), rescaleaccelMD(), rescaleVelocities(), simParams, slowFreq, Node::specialTracing, SimParameters::statsOn, SimParameters::stepsPerCycle, tcoupleVelocities(), traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum. Referenced by algorithm(). 00313 {
00314 char traceNote[24];
00315
00316 int step = simParams->firstTimestep;
00317
00318 const int numberOfSteps = simParams->N;
00319 const int stepsPerCycle = simParams->stepsPerCycle;
00320
00321 const int zeroMomentum = simParams->zeroMomentum;
00322
00323 nbondFreq = simParams->nonbondedFrequency;
00324 const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00325 if (dofull)
00326 slowFreq = simParams->fullElectFrequency;
00327 else
00328 slowFreq = simParams->nonbondedFrequency;
00329 if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
00330
00331 reassignVelocities(step); // only for full-step velecities
00332 rescaleaccelMD(step);
00333 adaptTempUpdate(step); // Init adaptive tempering;
00334
00335 receivePressure(step);
00336 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00337 correctMomentum(step);
00338 printFepMessage(step);
00339 printTiMessage(step);
00340 printDynamicsEnergies(step);
00341 outputFepEnergy(step);
00342 outputTiEnergy(step);
00343 if(traceIsOn()){
00344 traceUserEvent(eventEndOfTimeStep);
00345 sprintf(traceNote, "s:%d", step);
00346 traceUserSuppliedNote(traceNote);
00347 }
00348 outputExtendedSystem(step);
00349 rebalanceLoad(step);
00350
00351 // Handling SIGINT doesn't seem to be working on Lemieux, and it
00352 // sometimes causes the net-xxx versions of NAMD to segfault on exit,
00353 // so disable it for now.
00354 // namd_sighandler_t oldhandler = signal(SIGINT,
00355 // (namd_sighandler_t)my_sigint_handler);
00356 for ( ++step ; step <= numberOfSteps; ++step )
00357 {
00358 adaptTempUpdate(step);
00359 rescaleVelocities(step);
00360 tcoupleVelocities(step);
00361 berendsenPressure(step);
00362 langevinPiston1(step);
00363 rescaleaccelMD(step);
00364 enqueueCollections(step); // after lattice scaling!
00365 receivePressure(step);
00366 if ( zeroMomentum && dofull && ! (step % slowFreq) )
00367 correctMomentum(step);
00368 langevinPiston2(step);
00369 reassignVelocities(step);
00370 printDynamicsEnergies(step);
00371 outputFepEnergy(step);
00372 outputTiEnergy(step);
00373 if(traceIsOn()){
00374 traceUserEvent(eventEndOfTimeStep);
00375 sprintf(traceNote, "s:%d", step);
00376 traceUserSuppliedNote(traceNote);
00377 }
00378 // if (gotsigint) {
00379 // iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
00380 // NAMD_quit();
00381 // }
00382 outputExtendedSystem(step);
00383 #if CYCLE_BARRIER
00384 cycleBarrier(!((step+1) % stepsPerCycle),step);
00385 #elif PME_BARRIER
00386 cycleBarrier(dofull && !(step%slowFreq),step);
00387 #elif STEP_BARRIER
00388 cycleBarrier(1, step);
00389 #endif
00390
00391 if(Node::Object()->specialTracing || simParams->statsOn){
00392 int bstep = simParams->traceStartStep;
00393 int estep = bstep + simParams->numTraceSteps;
00394 if(step == bstep){
00395 traceBarrier(1, step);
00396 }
00397 if(step == estep){
00398 traceBarrier(0, step);
00399 }
00400 }
00401
00402 #ifdef MEASURE_NAMD_WITH_PAPI
00403 if(simParams->papiMeasure) {
00404 int bstep = simParams->papiMeasureStartStep;
00405 int estep = bstep + simParams->numPapiMeasureSteps;
00406 if(step == bstep) {
00407 papiMeasureBarrier(1, step);
00408 }
00409 if(step == estep) {
00410 papiMeasureBarrier(0, step);
00411 }
00412 }
00413 #endif
00414
00415 rebalanceLoad(step);
00416
00417 #if PME_BARRIER
00418 cycleBarrier(dofull && !((step+1)%slowFreq),step); // step before PME
00419 #endif
00420 }
00421 // signal(SIGINT, oldhandler);
00422 }
|
|
|
Definition at line 684 of file Controller.C. References BigReal, BOLTZMANN, broadcast, Lattice::c(), controlNumDegFreedom, controlPressure, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, NamdState::lattice, nbondFreq, ControllerBroadcasts::positionRescaleFactor, SimpleBroadcastObject< T >::publish(), random, Lattice::rescale(), simParams, slowFreq, state, SimParameters::surfaceTensionTarget, SimParameters::useConstantArea, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz. Referenced by integrate(). 00685 {
00686 if ( simParams->langevinPistonOn )
00687 {
00688 Tensor &strainRate = langevinPiston_strainRate;
00689 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00690 BigReal dt = simParams->dt;
00691 BigReal dt_long = slowFreq * dt;
00692 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00693 BigReal tau = simParams->langevinPistonPeriod;
00694 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00695
00696 #ifdef DEBUG_PRESSURE
00697 iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
00698 #endif
00699
00700 if ( ! ( (step-1) % slowFreq ) )
00701 {
00702 BigReal gamma = 1 / simParams->langevinPistonDecay;
00703 BigReal f1 = exp( -0.5 * dt_long * gamma );
00704 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00705 strainRate *= f1;
00706 if ( simParams->useFlexibleCell ) {
00707 // We only use on-diagonal terms (for now)
00708 if ( simParams->useConstantRatio ) {
00709 BigReal r = f2 * random->gaussian();
00710 strainRate.xx += r;
00711 strainRate.yy += r;
00712 strainRate.zz += f2 * random->gaussian();
00713 } else {
00714 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00715 }
00716 } else {
00717 strainRate += f2 * Tensor::identity(random->gaussian());
00718 }
00719
00720 #ifdef DEBUG_PRESSURE
00721 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00722 #endif
00723 }
00724
00725 // Apply surface tension. If surfaceTensionTarget is zero, we get
00726 // the default (isotropic pressure) case.
00727
00728 Tensor ptarget;
00729 ptarget.zz = simParams->langevinPistonTarget;
00730 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00731 simParams->surfaceTensionTarget / state->lattice.c().z;
00732
00733 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00734 ( controlPressure - ptarget );
00735
00736 if (simParams->fixCellDims) {
00737 if (simParams->fixCellDimX) strainRate.xx = 0;
00738 if (simParams->fixCellDimY) strainRate.yy = 0;
00739 if (simParams->fixCellDimZ) strainRate.zz = 0;
00740 }
00741
00742
00743 #ifdef DEBUG_PRESSURE
00744 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00745 #endif
00746
00747 if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
00748 {
00749 // We only use on-diagonal terms (for now)
00750 Tensor factor;
00751 if ( !simParams->useConstantArea ) {
00752 factor.xx = exp( dt_long * strainRate.xx );
00753 factor.yy = exp( dt_long * strainRate.yy );
00754 } else {
00755 factor.xx = factor.yy = 1;
00756 }
00757 factor.zz = exp( dt_long * strainRate.zz );
00758 broadcast->positionRescaleFactor.publish(step,factor);
00759 state->lattice.rescale(factor);
00760 #ifdef DEBUG_PRESSURE
00761 iout << iINFO << "rescaling by: " << factor << "\n";
00762 #endif
00763 }
00764
00765 // corrections to integrator
00766 if ( ! ( step % nbondFreq ) )
00767 {
00768 #ifdef DEBUG_PRESSURE
00769 iout << iINFO << "correcting strain rate for nbond, ";
00770 #endif
00771 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00772 ( (nbondFreq - 1) * controlPressure_nbond );
00773 #ifdef DEBUG_PRESSURE
00774 iout << "strain rate: " << strainRate << "\n";
00775 #endif
00776 }
00777 if ( ! ( step % slowFreq ) )
00778 {
00779 #ifdef DEBUG_PRESSURE
00780 iout << iINFO << "correcting strain rate for slow, ";
00781 #endif
00782 strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00783 ( (slowFreq - 1) * controlPressure_slow );
00784 #ifdef DEBUG_PRESSURE
00785 iout << "strain rate: " << strainRate << "\n";
00786 #endif
00787 }
00788 if (simParams->fixCellDims) {
00789 if (simParams->fixCellDimX) strainRate.xx = 0;
00790 if (simParams->fixCellDimY) strainRate.yy = 0;
00791 if (simParams->fixCellDimZ) strainRate.zz = 0;
00792 }
00793
00794 }
00795
00796 }
|
|
|
Definition at line 798 of file Controller.C. References BigReal, BOLTZMANN, Lattice::c(), controlNumDegFreedom, controlPressure, Tensor::diagonal(), SimParameters::dt, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, Random::gaussian(), Random::gaussian_vector(), Tensor::identity(), iINFO(), iout, SimParameters::langevinPistonDecay, SimParameters::langevinPistonOn, SimParameters::langevinPistonPeriod, SimParameters::langevinPistonTarget, SimParameters::langevinPistonTemp, NamdState::lattice, nbondFreq, random, simParams, slowFreq, state, SimParameters::surfaceTensionTarget, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, Lattice::volume(), Tensor::xx, Tensor::yy, Vector::z, and Tensor::zz. Referenced by integrate(). 00799 {
00800 if ( simParams->langevinPistonOn )
00801 {
00802 Tensor &strainRate = langevinPiston_strainRate;
00803 int cellDims = simParams->useFlexibleCell ? 1 : 3;
00804 BigReal dt = simParams->dt;
00805 BigReal dt_long = slowFreq * dt;
00806 BigReal kT = BOLTZMANN * simParams->langevinPistonTemp;
00807 BigReal tau = simParams->langevinPistonPeriod;
00808 BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
00809
00810 // corrections to integrator
00811 if ( ! ( step % nbondFreq ) )
00812 {
00813 #ifdef DEBUG_PRESSURE
00814 iout << iINFO << "correcting strain rate for nbond, ";
00815 #endif
00816 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00817 ( (nbondFreq - 1) * controlPressure_nbond );
00818 #ifdef DEBUG_PRESSURE
00819 iout << "strain rate: " << strainRate << "\n";
00820 #endif
00821 }
00822 if ( ! ( step % slowFreq ) )
00823 {
00824 #ifdef DEBUG_PRESSURE
00825 iout << iINFO << "correcting strain rate for slow, ";
00826 #endif
00827 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00828 ( (slowFreq - 1) * controlPressure_slow );
00829 #ifdef DEBUG_PRESSURE
00830 iout << "strain rate: " << strainRate << "\n";
00831 #endif
00832 }
00833
00834 // Apply surface tension. If surfaceTensionTarget is zero, we get
00835 // the default (isotropic pressure) case.
00836
00837 Tensor ptarget;
00838 ptarget.zz = simParams->langevinPistonTarget;
00839 ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
00840 simParams->surfaceTensionTarget / state->lattice.c().z;
00841
00842 strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
00843 ( controlPressure - ptarget );
00844
00845
00846 #ifdef DEBUG_PRESSURE
00847 iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
00848 #endif
00849
00850 if ( ! ( step % slowFreq ) )
00851 {
00852 BigReal gamma = 1 / simParams->langevinPistonDecay;
00853 BigReal f1 = exp( -0.5 * dt_long * gamma );
00854 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
00855 strainRate *= f1;
00856 if ( simParams->useFlexibleCell ) {
00857 // We only use on-diagonal terms (for now)
00858 if ( simParams->useConstantRatio ) {
00859 BigReal r = f2 * random->gaussian();
00860 strainRate.xx += r;
00861 strainRate.yy += r;
00862 strainRate.zz += f2 * random->gaussian();
00863 } else {
00864 strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
00865 }
00866 } else {
00867 strainRate += f2 * Tensor::identity(random->gaussian());
00868 }
00869 #ifdef DEBUG_PRESSURE
00870 iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
00871 #endif
00872 }
00873
00874 #ifdef DEBUG_PRESSURE
00875 iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
00876 #endif
00877 if (simParams->fixCellDims) {
00878 if (simParams->fixCellDimX) strainRate.xx = 0;
00879 if (simParams->fixCellDimY) strainRate.yy = 0;
00880 if (simParams->fixCellDimZ) strainRate.zz = 0;
00881 }
00882 }
00883
00884
00885 }
|
|
|
Definition at line 465 of file Controller.C. References BigReal, broadcast, CALCULATE, minpoint::dudx, enqueueCollections(), SimParameters::firstTimestep, iout, min_f_dot_f, min_f_dot_v, min_huge_count, min_v_dot_v, SimParameters::minBabyStep, ControllerBroadcasts::minimizeCoefficient, SimParameters::minLineGoal, SimParameters::minTinyStep, SimParameters::minVerbose, MOVETO, SimParameters::N, nbondFreq, minpoint::noGradient, PRINT_BRACKET, SimpleBroadcastObject< T >::publish(), simParams, slowFreq, minpoint::u, and minpoint::x. Referenced by algorithm(). 00465 {
00466 // iout << "Controller::minimize() called.\n" << endi;
00467
00468 const int minVerbose = simParams->minVerbose;
00469 const int numberOfSteps = simParams->N;
00470 int step = simParams->firstTimestep;
00471 slowFreq = nbondFreq = 1;
00472 BigReal tinystep = simParams->minTinyStep; // 1.0e-6
00473 BigReal babystep = simParams->minBabyStep; // 1.0e-2
00474 BigReal linegoal = simParams->minLineGoal; // 1.0e-3
00475 BigReal initstep = tinystep;
00476 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
00477
00478 CALCULATE
00479
00480 int minSeq = 0;
00481
00482 // just move downhill until initial bad contacts go away or 100 steps
00483 int old_min_huge_count = 2000000000;
00484 int steps_at_same_min_huge_count = 0;
00485 for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
00486 if ( min_huge_count >= old_min_huge_count ) {
00487 if ( ++steps_at_same_min_huge_count > 10 ) break;
00488 } else {
00489 old_min_huge_count = min_huge_count;
00490 steps_at_same_min_huge_count = 0;
00491 }
00492 iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
00493 broadcast->minimizeCoefficient.publish(minSeq++,1.);
00494 enqueueCollections(step);
00495 CALCULATE
00496 }
00497 if ( min_huge_count ) {
00498 iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
00499 }
00500 iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00501
00502 int atStart = 2;
00503 int errorFactor = 10;
00504 BigReal old_f_dot_f = min_f_dot_f;
00505 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00506 broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
00507 int newDir = 1;
00508 min_f_dot_v = min_f_dot_f;
00509 min_v_dot_v = min_f_dot_f;
00510 while ( 1 ) {
00511 // line minimization
00512 // bracket minimum on line
00513 minpoint lo,hi,mid,last;
00514 BigReal x = 0;
00515 lo.x = x;
00516 lo.u = min_energy;
00517 lo.dudx = -1. * min_f_dot_v;
00518 lo.noGradient = min_huge_count;
00519 mid = lo;
00520 last = mid;
00521 if ( minVerbose ) {
00522 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
00523 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
00524 iout << "\n" << endi;
00525 }
00526 BigReal tol = fabs( linegoal * min_f_dot_v );
00527 if ( initstep > babystep ) initstep = babystep;
00528 if ( initstep < 1.0e-300 ) initstep = 1.0e-300;
00529 iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
00530 fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
00531 int start_with_huge = last.noGradient;
00532 x = initstep;
00533 x *= sqrt( min_f_dot_f / min_v_dot_v ); MOVETO(x)
00534 if ( ! start_with_huge ) {
00535 for ( int i_huge = 0; last.noGradient && i_huge < 10; ++i_huge ) {
00536 x *= 0.25; MOVETO(x);
00537 initstep *= 0.25;
00538 }
00539 }
00540 // bracket minimum on line
00541 initstep *= 0.25;
00542 BigReal maxinitstep = initstep * 16.0;
00543 while ( last.u < mid.u ) {
00544 initstep *= 2.0;
00545 lo = mid; mid = last;
00546 x *= 2.0; MOVETO(x)
00547 }
00548 if ( initstep > maxinitstep ) initstep = maxinitstep;
00549 hi = last;
00550 #define PRINT_BRACKET \
00551 iout << "LINE MINIMIZER BRACKET: DX " \
00552 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
00553 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
00554 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
00555 PRINT_BRACKET
00556 // converge on minimum on line
00557 int itcnt;
00558 for ( itcnt = 10; itcnt > 0; --itcnt ) {
00559 int progress = 1;
00560 // select new position
00561 if ( mid.noGradient ) {
00562 if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) { // subdivide left side
00563 x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
00564 MOVETO(x)
00565 if ( last.u <= mid.u ) {
00566 hi = mid; mid = last;
00567 } else {
00568 lo = last;
00569 }
00570 } else { // subdivide right side
00571 x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
00572 MOVETO(x)
00573 if ( last.u <= mid.u ) {
00574 lo = mid; mid = last;
00575 } else {
00576 hi = last;
00577 }
00578 }
00579 } else {
00580 if ( mid.dudx > 0. ) { // subdivide left side
00581 BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
00582 BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
00583 x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
00584 BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
00585 if ( xdiv ) x /= xdiv; else x = last.x;
00586 if ( x > altxhi ) x = altxhi;
00587 if ( x < altxlo ) x = altxlo;
00588 if ( x-last.x == 0 ) break;
00589 MOVETO(x)
00590 if ( last.u <= mid.u ) {
00591 hi = mid; mid = last;
00592 } else {
00593 if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
00594 lo = last;
00595 }
00596 } else { // subdivide right side
00597 BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
00598 BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
00599 x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
00600 BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
00601 if ( xdiv ) x /= xdiv; else x = last.x;
00602 if ( x < altxlo ) x = altxlo;
00603 if ( x > altxhi ) x = altxhi;
00604 if ( x-last.x == 0 ) break;
00605 MOVETO(x)
00606 if ( last.u <= mid.u ) {
00607 lo = mid; mid = last;
00608 } else {
00609 if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
00610 hi = last;
00611 }
00612 }
00613 }
00614 PRINT_BRACKET
00615 if ( ! progress ) {
00616 MOVETO(mid.x);
00617 break;
00618 }
00619 if ( fabs(last.dudx) < tol ) break;
00620 if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
00621 if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
00622 }
00623 // new direction
00624 broadcast->minimizeCoefficient.publish(minSeq++,0.);
00625 BigReal c = min_f_dot_f / old_f_dot_f;
00626 c = ( c > 1.5 ? 1.5 : c );
00627 if ( atStart ) { c = 0; --atStart; }
00628 if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
00629 c = 0;
00630 if ( errorFactor < 100 ) errorFactor += 10;
00631 }
00632 if ( c == 0 ) {
00633 iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
00634 }
00635 broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
00636 newDir = 1;
00637 old_f_dot_f = min_f_dot_f;
00638 min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
00639 min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
00640 }
00641
00642 }
|
|
|
Definition at line 2724 of file Controller.C. References FILE_OUTPUT, SimParameters::firstTimestep, iout, SimParameters::N, NAMD_backup_file(), NAMD_err(), SimParameters::outputFilename, SimParameters::restartFilename, SimParameters::restartFrequency, SimParameters::restartSave, simParams, writeExtendedSystemData(), writeExtendedSystemLabels(), xstFile, SimParameters::xstFilename, and SimParameters::xstFrequency. Referenced by algorithm(), and integrate(). 02725 {
02726
02727 if ( step >= 0 ) {
02728
02729 // Write out eXtended System Trajectory (XST) file
02730 if ( simParams->xstFrequency &&
02731 ((step % simParams->xstFrequency) == 0) )
02732 {
02733 if ( ! xstFile.rdbuf()->is_open() )
02734 {
02735 iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02736 NAMD_backup_file(simParams->xstFilename);
02737 xstFile.open(simParams->xstFilename);
02738 while (!xstFile) {
02739 if ( errno == EINTR ) {
02740 CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
02741 xstFile.clear();
02742 xstFile.open(simParams->xstFilename);
02743 continue;
02744 }
02745 char err_msg[257];
02746 sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
02747 NAMD_err(err_msg);
02748 }
02749 xstFile << "# NAMD extended system trajectory file" << std::endl;
02750 writeExtendedSystemLabels(xstFile);
02751 }
02752 writeExtendedSystemData(step,xstFile);
02753 xstFile.flush();
02754 }
02755
02756 // Write out eXtended System Configuration (XSC) files
02757 // Output a restart file
02758 if ( simParams->restartFrequency &&
02759 ((step % simParams->restartFrequency) == 0) &&
02760 (step != simParams->firstTimestep) )
02761 {
02762 iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
02763 << step << "\n" << endi;
02764 char fname[140];
02765 strcpy(fname, simParams->restartFilename);
02766 if ( simParams->restartSave ) {
02767 char timestepstr[20];
02768 sprintf(timestepstr,".%d",step);
02769 strcat(fname, timestepstr);
02770 }
02771 strcat(fname, ".xsc");
02772 NAMD_backup_file(fname,".old");
02773 std::ofstream xscFile(fname);
02774 while (!xscFile) {
02775 if ( errno == EINTR ) {
02776 CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
02777 xscFile.clear();
02778 xscFile.open(fname);
02779 continue;
02780 }
02781 char err_msg[257];
02782 sprintf(err_msg, "Error opening XSC restart file %s",fname);
02783 NAMD_err(err_msg);
02784 }
02785 xscFile << "# NAMD extended system configuration restart file" << std::endl;
02786 writeExtendedSystemLabels(xscFile);
02787 writeExtendedSystemData(step,xscFile);
02788 if (!xscFile) {
02789 char err_msg[257];
02790 sprintf(err_msg, "Error writing XSC restart file %s",fname);
02791 NAMD_err(err_msg);
02792 }
02793 }
02794
02795 }
02796
02797 // Output final coordinates
02798 if (step == FILE_OUTPUT || step == END_OF_RUN)
02799 {
02800 int realstep = ( step == FILE_OUTPUT ?
02801 simParams->firstTimestep : simParams->N );
02802 iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
02803 << realstep << "\n" << endi;
02804 static char fname[140];
02805 strcpy(fname, simParams->outputFilename);
02806 strcat(fname, ".xsc");
02807 NAMD_backup_file(fname);
02808 std::ofstream xscFile(fname);
02809 while (!xscFile) {
02810 if ( errno == EINTR ) {
02811 CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
02812 xscFile.clear();
02813 xscFile.open(fname);
02814 continue;
02815 }
02816 char err_msg[257];
02817 sprintf(err_msg, "Error opening XSC output file %s",fname);
02818 NAMD_err(err_msg);
02819 }
02820 xscFile << "# NAMD extended system configuration output file" << std::endl;
02821 writeExtendedSystemLabels(xscFile);
02822 writeExtendedSystemData(realstep,xscFile);
02823 if (!xscFile) {
02824 char err_msg[257];
02825 sprintf(err_msg, "Error writing XSC output file %s",fname);
02826 NAMD_err(err_msg);
02827 }
02828 }
02829
02830 // Close trajectory file
02831 if (step == END_OF_RUN) {
02832 if ( xstFile.rdbuf()->is_open() ) {
02833 xstFile.close();
02834 iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
02835 }
02836 }
02837
02838 }
|
|
|
Definition at line 2545 of file Controller.C. References SimParameters::alchEnsembleAvg, SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchTemp, BigReal, BOLTZMANN, electEnergy, electEnergy_f, electEnergySlow, electEnergySlow_f, exp_dE_ByRT, fepFile, FepNo, fepSum, SimParameters::firstTimestep, iout, ljEnergy_f, SimParameters::N, NAMD_backup_file(), net_dE, simParams, and writeFepEnergyData(). Referenced by integrate(). 02545 {
02546 if (simParams->alchFepOn) {
02547 const int stepInRun = step - simParams->firstTimestep;
02548 const int alchEquilSteps = simParams->alchEquilSteps;
02549 const BigReal alchLambda = simParams->alchLambda;
02550 const BigReal alchLambda2 = simParams->alchLambda2;
02551 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
02552 if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
02553 FepNo = 0;
02554 exp_dE_ByRT = 0.0;
02555 net_dE = 0.0;
02556 }
02557 BigReal dE = electEnergy_f + electEnergySlow_f + ljEnergy_f
02558 - (electEnergy + electEnergySlow + ljEnergy);
02559 BigReal RT = BOLTZMANN * simParams->alchTemp;
02560
02561 if (alchEnsembleAvg){
02562 FepNo++;
02563 exp_dE_ByRT += exp(-dE/RT);
02564 net_dE += dE;
02565 }
02566
02567 if (stepInRun == 0) {
02568 if (!fepFile.rdbuf()->is_open()) {
02569 NAMD_backup_file(simParams->alchOutFile);
02570 fepFile.open(simParams->alchOutFile);
02571 iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
02572 if(alchEnsembleAvg){
02573 fepSum = 0.0;
02574 fepFile << "# STEP Elec "
02575 << "vdW dE dE_avg Temp dG\n"
02576 << "# l l+dl "
02577 << " l l+dl E(l+dl)-E(l)" << std::endl;
02578 }
02579 else{
02580 fepFile << "# STEP Elec "
02581 << "vdW dE Temp\n"
02582 << "# l l+dl "
02583 << " l l+dl E(l+dl)-E(l)" << std::endl;
02584 }
02585 }
02586 if(!step){
02587 fepFile << "#NEW FEP WINDOW: "
02588 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
02589 << alchLambda2 << std::endl;
02590 }
02591 }
02592 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
02593 fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
02594 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
02595 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
02596 }
02597 if ((simParams->N) && (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0))) {
02598 writeFepEnergyData(step, fepFile);
02599 fepFile.flush();
02600 }
02601 if (alchEnsembleAvg && (step == simParams->N)) {
02602 fepSum = fepSum + dG;
02603 fepFile << "#Free energy change for lambda window [ " << alchLambda
02604 << " " << alchLambda2 << " ] is " << dG << " ; net change until now is " << fepSum << std::endl;
02605 }
02606 }
02607 }
|
|
|
Definition at line 2609 of file Controller.C. References SimParameters::alchElecLambdaStart, SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchOutFile, SimParameters::alchOutFreq, SimParameters::alchThermIntOn, SimParameters::alchVdwLambdaEnd, BigReal, electEnergy_ti_1, electEnergy_ti_2, electEnergySlow_ti_1, electEnergySlow_ti_2, SimParameters::firstTimestep, iout, NAMD_backup_file(), net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, simParams, tiFile, TiNo, and writeTiEnergyData(). Referenced by integrate(). 02609 {
02610 if (simParams->alchThermIntOn) {
02611 const int stepInRun = step - simParams->firstTimestep;
02612 const int alchEquilSteps = simParams->alchEquilSteps;
02613 const BigReal alchLambda = simParams->alchLambda;
02614
02615 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
02616 TiNo = 0;
02617 net_dEdl_elec_1 = 0;
02618 net_dEdl_elec_2 = 0;
02619 net_dEdl_lj_1 = 0;
02620 net_dEdl_lj_2 = 0;
02621 }
02622 if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) {
02623 // output of instantaneous dU/dl now replaced with running average
02624 // over last alchOutFreq steps (except for step 0)
02625 recent_TiNo = 0;
02626 recent_dEdl_elec_1 = 0;
02627 recent_dEdl_elec_2 = 0;
02628 recent_dEdl_lj_1 = 0;
02629 recent_dEdl_lj_2 = 0;
02630 }
02631 TiNo++;
02632 recent_TiNo++;
02633 // FB - PME is no longer scaled by global lambda, but by the respective
02634 // lambda as dictated by elecLambdaStart. All electrostatics now go together.
02635 net_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
02636 net_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
02637 net_dEdl_lj_1 += ljEnergy_ti_1;
02638 net_dEdl_lj_2 += ljEnergy_ti_2;
02639 recent_dEdl_elec_1 += electEnergy_ti_1 + electEnergySlow_ti_1 + electEnergyPME_ti_1;
02640 recent_dEdl_elec_2 += electEnergy_ti_2 + electEnergySlow_ti_2 + electEnergyPME_ti_2;
02641 recent_dEdl_lj_1 += ljEnergy_ti_1;
02642 recent_dEdl_lj_2 += ljEnergy_ti_2;
02643
02644 if (stepInRun == 0) {
02645 BigReal alchElecLambdaStart = simParams->alchElecLambdaStart;
02646 BigReal alchVdwLambdaEnd = simParams->alchVdwLambdaEnd;
02647 BigReal elec_lambda_1 = (alchLambda <= alchElecLambdaStart)? 0. : \
02648 (alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02649 BigReal elec_lambda_2 = ((1-alchLambda) <= alchElecLambdaStart)? 0. : \
02650 ((1-alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02651 BigReal vdw_lambda_1 = (alchLambda >= alchVdwLambdaEnd)? 1. : \
02652 alchLambda / alchVdwLambdaEnd;
02653 BigReal vdw_lambda_2 = ((1-alchLambda) >= alchVdwLambdaEnd)? 1. : \
02654 (1-alchLambda) / alchVdwLambdaEnd;
02655 if (!tiFile.rdbuf()->is_open()) {
02656 //tiSum = 0.0;
02657 NAMD_backup_file(simParams->alchOutFile);
02658 tiFile.open(simParams->alchOutFile);
02659 iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
02660 tiFile << "# STEP Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg Elec_dU/dl Elec_avg vdW_dU/dl vdw_avg PME_dU/dl PME_avg\n"
02661 << "# <---------------------PARTITION 1------------------------> <---------------------PARTITION 2--------------------->"
02662 << std::endl;
02663 }
02664 tiFile << "#NEW TI WINDOW: "
02665 << "LAMBDA " << alchLambda
02666 << "\n#PARTITION 1 VDW LAMBDA " << vdw_lambda_1
02667 << "\n#PARTITION 1 ELEC LAMBDA " << elec_lambda_1
02668 << "\n#PARTITION 2 VDW LAMBDA " << vdw_lambda_2
02669 << "\n#PARTITION 2 ELEC LAMBDA " << elec_lambda_2
02670 << "\n" << std::endl;
02671 }
02672 if (stepInRun == alchEquilSteps) {
02673 tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
02674 << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
02675 << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
02676 }
02677 if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) {
02678 writeTiEnergyData(step, tiFile);
02679 tiFile.flush();
02680 }
02681 }
02682 }
|
|
|
Definition at line 2040 of file Controller.C. References compareChecksums(), NamdState::lattice, Node::molecule, Node::Object(), printEnergies(), Node::simParameters, and state. Referenced by integrate(). 02040 {
02041
02042 Node *node = Node::Object();
02043 Molecule *molecule = node->molecule;
02044 SimParameters *simParameters = node->simParameters;
02045 Lattice &lattice = state->lattice;
02046
02047 compareChecksums(step);
02048
02049 printEnergies(step,0);
02050 }
|
|
||||||||||||
|
Definition at line 2052 of file Controller.C. References Lattice::a(), Lattice::a_p(), avg_count, Lattice::b(), Lattice::b_p(), BigReal, Lattice::c(), Lattice::c_p(), NamdState::callback_labelstring, NamdState::callback_valuestring, CALLBACKDATA, CALLBACKLIST, ScriptTcl::doCallback(), drudeBondTemp, drudeBondTempAvg, drudeComTemp, drudeComTempAvg, SimParameters::drudeOn, SimParameters::dt, IMDEnergies::Eangle, IMDEnergies::Ebond, IMDEnergies::Edihe, IMDEnergies::Eelec, IMDEnergies::Eimpr, electEnergy, electEnergy_f, electEnergy_ti_1, electEnergy_ti_2, electEnergyPME_ti_1, electEnergyPME_ti_2, electEnergySlow, electEnergySlow_f, electEnergySlow_ti_1, electEnergySlow_ti_2, IMDEnergies::Epot, ETITLE(), IMDEnergies::Etot, IMDEnergies::Evdw, SimParameters::firstLdbStep, SimParameters::firstTimestep, FORMAT(), IMDOutput::gather_energies(), GET_VECTOR, PressureProfileReduction::getData(), Node::getScript(), SimParameters::goForcesOn, goNativeEnergy, goNonnativeEnergy, goTotalEnergy, groupPressure, groupPressure_avg, groupPressure_tavg, iERROR(), iINFO(), Node::imd, SimParameters::IMDfreq, SimParameters::IMDon, iout, RequireReduction::item(), iWARN(), kineticEnergyHalfstep, NamdState::lattice, SimParameters::LJcorrection, ljEnergy, ljEnergy_f, ljEnergy_ti_1, ljEnergy_ti_2, marginViolations, memusage_MB(), SimParameters::mergeCrossterms, Node::molecule, SimParameters::N, Node::Object(), Lattice::origin(), SimParameters::outputEnergies, SimParameters::outputMomenta, SimParameters::outputPairlists, SimParameters::outputPressure, SimParameters::pairInteractionOn, pairlistWarnings, ppbonded, ppint, ppnonbonded, pressure, pressure_avg, pressure_tavg, PRESSUREFACTOR, pressureProfileAverage, pressureProfileCount, SimParameters::pressureProfileFreq, printTiming(), reduction, REDUCTION_ANGLE_ENERGY, REDUCTION_BC_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_F, REDUCTION_ELECT_ENERGY_PME_TI_1, REDUCTION_ELECT_ENERGY_PME_TI_2, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_ELECT_ENERGY_SLOW_F, REDUCTION_ELECT_ENERGY_SLOW_TI_1, REDUCTION_ELECT_ENERGY_SLOW_TI_2, REDUCTION_ELECT_ENERGY_TI_1, REDUCTION_ELECT_ENERGY_TI_2, REDUCTION_GO_NATIVE_ENERGY, REDUCTION_GO_NONNATIVE_ENERGY, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_LJ_ENERGY_F, REDUCTION_LJ_ENERGY_TI_1, REDUCTION_LJ_ENERGY_TI_2, REDUCTION_MISC_ENERGY, REDUCTION_PAIR_ELECT_FORCE, REDUCTION_PAIR_VDW_FORCE, Node::simParameters, simParams, smooth2_avg, state, IMDEnergies::T, Molecule::tail_corr_ener, tavg_count, temp_avg, temperature, totalEnergy, IMDEnergies::tstep, Lattice::volume(), Vector::x, Vector::y, and Vector::z. Referenced by printDynamicsEnergies(), and printMinimizeEnergies(). 02053 {
02054 Node *node = Node::Object();
02055 Molecule *molecule = node->molecule;
02056 SimParameters *simParameters = node->simParameters;
02057 Lattice &lattice = state->lattice;
02058
02059 // Drude model ANISO energy is added into BOND energy
02060 // and THOLE energy is added into ELECT energy
02061
02062 BigReal bondEnergy;
02063 BigReal angleEnergy;
02064 BigReal dihedralEnergy;
02065 BigReal improperEnergy;
02066 BigReal crosstermEnergy;
02067 BigReal boundaryEnergy;
02068 BigReal miscEnergy;
02069 BigReal potentialEnergy;
02070 BigReal flatEnergy;
02071 BigReal smoothEnergy;
02072
02073 Vector momentum;
02074 Vector angularMomentum;
02075 BigReal volume = lattice.volume();
02076
02077 bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
02078 angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
02079 dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
02080 improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
02081 crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
02082 boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
02083 miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
02084
02085 if ( minimize || ! ( step % nbondFreq ) )
02086 {
02087 electEnergy = reduction->item(REDUCTION_ELECT_ENERGY);
02088 ljEnergy = reduction->item(REDUCTION_LJ_ENERGY);
02089
02090 // JLai
02091 goNativeEnergy = reduction->item(REDUCTION_GO_NATIVE_ENERGY);
02092 goNonnativeEnergy = reduction->item(REDUCTION_GO_NONNATIVE_ENERGY);
02093 goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
02094
02095 //fepb
02096 electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
02097 ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
02098
02099 electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
02100 ljEnergy_ti_1 = reduction->item(REDUCTION_LJ_ENERGY_TI_1);
02101 electEnergy_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_TI_2);
02102 ljEnergy_ti_2 = reduction->item(REDUCTION_LJ_ENERGY_TI_2);
02103 //fepe
02104 }
02105
02106 if ( minimize || ! ( step % slowFreq ) )
02107 {
02108 electEnergySlow = reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
02109 //fepb
02110 electEnergySlow_f = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F);
02111
02112 electEnergySlow_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1);
02113 electEnergySlow_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_2);
02114 electEnergyPME_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1);
02115 electEnergyPME_ti_2 = reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2);
02116 //fepe
02117 }
02118
02119 if (simParameters->LJcorrection && volume) {
02120 // Apply tail correction to energy
02121 //printf("Volume is %f\n", volume);
02122 //printf("Applying tail correction of %f to energy\n", molecule->tail_corr_ener / volume);
02123 ljEnergy += molecule->tail_corr_ener / volume;
02124 }
02125
02126
02127 momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
02128 momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
02129 momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
02130 angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
02131 angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
02132 angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
02133
02134 // Ported by JLai
02135 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
02136 improperEnergy + electEnergy + electEnergySlow + ljEnergy +
02137 crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy;
02138 // End of port
02139 totalEnergy = potentialEnergy + kineticEnergy;
02140 flatEnergy = totalEnergy +
02141 (1.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02142 if ( !(step%slowFreq) ) {
02143 // only adjust based on most accurate energies
02144 BigReal s = (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02145 if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
02146 if ( step != simParams->firstTimestep ) {
02147 smooth2_avg *= 0.9375;
02148 smooth2_avg += 0.0625 * s;
02149 }
02150 }
02151 smoothEnergy = flatEnergy + smooth2_avg -
02152 (4.0/3.0)*( kineticEnergyHalfstep - kineticEnergyCentered);
02153
02154 if ( simParameters->outputMomenta && ! minimize &&
02155 ! ( step % simParameters->outputMomenta ) )
02156 {
02157 iout << "MOMENTUM: " << step
02158 << " P: " << momentum
02159 << " L: " << angularMomentum
02160 << "\n" << endi;
02161 }
02162
02163 if ( simParameters->outputPressure ) {
02164 pressure_tavg += pressure;
02165 groupPressure_tavg += groupPressure;
02166 tavg_count += 1;
02167 if ( minimize || ! ( step % simParameters->outputPressure ) ) {
02168 iout << "PRESSURE: " << step << " "
02169 << PRESSUREFACTOR * pressure << "\n"
02170 << "GPRESSURE: " << step << " "
02171 << PRESSUREFACTOR * groupPressure << "\n";
02172 if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
02173 << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
02174 << "GPRESSAVG: " << step << " "
02175 << (PRESSUREFACTOR/tavg_count) * groupPressure_tavg << "\n";
02176 iout << endi;
02177 pressure_tavg = 0;
02178 groupPressure_tavg = 0;
02179 tavg_count = 0;
02180 }
02181 }
02182
02183 // pressure profile reductions
02184 if (pressureProfileSlabs) {
02185 const int freq = simParams->pressureProfileFreq;
02186 const int arraysize = 3*pressureProfileSlabs;
02187
02188 BigReal *total = new BigReal[arraysize];
02189 memset(total, 0, arraysize*sizeof(BigReal));
02190 const int first = simParams->firstTimestep;
02191
02192 if (ppbonded) ppbonded->getData(first, step, lattice, total);
02193 if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
02194 if (ppint) ppint->getData(first, step, lattice, total);
02195 for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
02196 pressureProfileCount++;
02197
02198 if (!(step % freq)) {
02199 // convert NAMD internal virial to pressure in units of bar
02200 BigReal scalefac = PRESSUREFACTOR * pressureProfileSlabs;
02201
02202 iout << "PRESSUREPROFILE: " << step << " ";
02203 if (step == first) {
02204 // output pressure profile for this step
02205 for (int i=0; i<arraysize; i++) {
02206 iout << total[i] * scalefac << " ";
02207 }
02208 } else {
02209 // output pressure profile averaged over the last count steps.
02210 scalefac /= pressureProfileCount;
02211 for (int i=0; i<arraysize; i++)
02212 iout << pressureProfileAverage[i]*scalefac << " ";
02213 }
02214 iout << "\n" << endi;
02215
02216 // Clear the average for the next block
02217 memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
02218 pressureProfileCount = 0;
02219 }
02220 delete [] total;
02221 }
02222
02223 int stepInRun = step - simParams->firstTimestep;
02224 if ( stepInRun % simParams->firstLdbStep == 0 ) {
02225 int benchPhase = stepInRun / simParams->firstLdbStep;
02226 if ( benchPhase > 0 && benchPhase < 10 ) {
02227 #ifdef NAMD_CUDA
02228 if ( simParams->outputEnergies < 60 ) {
02229 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
02230 }
02231 #endif
02232 iout << iINFO;
02233 if ( benchPhase < 4 ) iout << "Initial time: ";
02234 else iout << "Benchmark time: ";
02235 iout << CkNumPes() << " CPUs ";
02236 {
02237 BigReal wallPerStep =
02238 (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
02239 BigReal ns = simParams->dt / 1000000.0;
02240 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
02241 BigReal daysPerNano = wallPerStep * days / ns;
02242 iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
02243 iout << memusage_MB() << " MB memory\n" << endi;
02244 }
02245 }
02246 startBenchTime = CmiWallTimer();
02247 }
02248
02249 printTiming(step);
02250
02251 Vector pairVDWForce, pairElectForce;
02252 if ( simParameters->pairInteractionOn ) {
02253 GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
02254 GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
02255 }
02256
02257 // callback to Tcl with whatever we can
02258 #ifdef NAMD_TCL
02259 #define CALLBACKDATA(LABEL,VALUE) \
02260 labels << (LABEL) << " "; values << (VALUE) << " ";
02261 #define CALLBACKLIST(LABEL,VALUE) \
02262 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
02263 if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
02264 std::ostringstream labels, values;
02265 #if CMK_BLUEGENEL
02266 // the normal version below gives a compiler error
02267 values.precision(16);
02268 #else
02269 values << std::setprecision(16);
02270 #endif
02271 CALLBACKDATA("TS",step);
02272 CALLBACKDATA("BOND",bondEnergy);
02273 CALLBACKDATA("ANGLE",angleEnergy);
02274 CALLBACKDATA("DIHED",dihedralEnergy);
02275 CALLBACKDATA("CROSS",crosstermEnergy);
02276 CALLBACKDATA("IMPRP",improperEnergy);
02277 CALLBACKDATA("ELECT",electEnergy+electEnergySlow);
02278 CALLBACKDATA("VDW",ljEnergy);
02279 CALLBACKDATA("BOUNDARY",boundaryEnergy);
02280 CALLBACKDATA("MISC",miscEnergy);
02281 CALLBACKDATA("POTENTIAL",potentialEnergy);
02282 CALLBACKDATA("KINETIC",kineticEnergy);
02283 CALLBACKDATA("TOTAL",totalEnergy);
02284 CALLBACKDATA("TEMP",temperature);
02285 CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
02286 CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
02287 CALLBACKDATA("VOLUME",lattice.volume());
02288 CALLBACKLIST("CELL_A",lattice.a());
02289 CALLBACKLIST("CELL_B",lattice.b());
02290 CALLBACKLIST("CELL_C",lattice.c());
02291 CALLBACKLIST("CELL_O",lattice.origin());
02292 labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
02293 << lattice.b_p() << " " << lattice.c_p() << "} ";
02294 if ( simParameters->pairInteractionOn ) {
02295 CALLBACKLIST("VDW_FORCE",pairVDWForce);
02296 CALLBACKLIST("ELECT_FORCE",pairElectForce);
02297 }
02298
02299 labels << '\0'; values << '\0'; // insane but makes Linux work
02300 state->callback_labelstring = labels.str();
02301 state->callback_valuestring = values.str();
02302 // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
02303 }
02304 #undef CALLBACKDATA
02305 #endif
02306
02307 drudeComTempAvg += drudeComTemp;
02308 drudeBondTempAvg += drudeBondTemp;
02309
02310 temp_avg += temperature;
02311 pressure_avg += trace(pressure)/3.;
02312 groupPressure_avg += trace(groupPressure)/3.;
02313 avg_count += 1;
02314
02315 if ( simParams->outputPairlists && pairlistWarnings &&
02316 ! (step % simParams->outputPairlists) ) {
02317 iout << iINFO << pairlistWarnings <<
02318 " pairlist warnings in past " << simParams->outputPairlists <<
02319 " steps.\n" << endi;
02320 pairlistWarnings = 0;
02321 }
02322
02323 // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
02324 if ( ! minimize && step % simParameters->outputEnergies ) return;
02325 // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
02326
02327 if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
02328 IMDEnergies energies;
02329 energies.tstep = step;
02330 energies.T = temp_avg/avg_count;
02331 energies.Etot = totalEnergy;
02332 energies.Epot = potentialEnergy;
02333 energies.Evdw = ljEnergy;
02334 energies.Eelec = electEnergy + electEnergySlow;
02335 energies.Ebond = bondEnergy;
02336 energies.Eangle = angleEnergy;
02337 energies.Edihe = dihedralEnergy + crosstermEnergy;
02338 energies.Eimpr = improperEnergy;
02339 Node::Object()->imd->gather_energies(&energies);
02340 }
02341
02342 if ( marginViolations ) {
02343 iout << iERROR << marginViolations <<
02344 " margin violations detected since previous energy output.\n" << endi;
02345 }
02346 marginViolations = 0;
02347
02348 if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
02349 {
02350 iout << "ETITLE: TS";
02351 iout << FORMAT("BOND");
02352 iout << FORMAT("ANGLE");
02353 iout << FORMAT("DIHED");
02354 if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
02355 iout << FORMAT("IMPRP");
02356 iout << " ";
02357 iout << FORMAT("ELECT");
02358 iout << FORMAT("VDW");
02359 iout << FORMAT("BOUNDARY");
02360 iout << FORMAT("MISC");
02361 iout << FORMAT("KINETIC");
02362 iout << " ";
02363 iout << FORMAT("TOTAL");
02364 iout << FORMAT("TEMP");
02365 iout << FORMAT("POTENTIAL");
02366 // iout << FORMAT("TOTAL2");
02367 iout << FORMAT("TOTAL3");
02368 iout << FORMAT("TEMPAVG");
02369 if ( volume != 0. ) {
02370 iout << " ";
02371 iout << FORMAT("PRESSURE");
02372 iout << FORMAT("GPRESSURE");
02373 iout << FORMAT("VOLUME");
02374 iout << FORMAT("PRESSAVG");
02375 iout << FORMAT("GPRESSAVG");
02376 }
02377 if (simParameters->drudeOn) {
02378 iout << " ";
02379 iout << FORMAT("DRUDECOM");
02380 iout << FORMAT("DRUDEBOND");
02381 iout << FORMAT("DRCOMAVG");
02382 iout << FORMAT("DRBONDAVG");
02383 }
02384 // Ported by JLai
02385 if (simParameters->goForcesOn) {
02386 iout << " ";
02387 iout << FORMAT("NATIVE");
02388 iout << FORMAT("NONNATIVE");
02389 //iout << FORMAT("REL_NATIVE");
02390 //iout << FORMAT("REL_NONNATIVE");
02391 iout << FORMAT("GOTOTAL");
02392 //iout << FORMAT("GOAVG");
02393 }
02394 // End of port -- JLai
02395 iout << "\n\n" << endi;
02396 }
02397
02398 // N.B. HP's aCC compiler merges FORMAT calls in the same expression.
02399 // Need separate statements because data returned in static array.
02400
02401 iout << ETITLE(step);
02402 iout << FORMAT(bondEnergy);
02403 iout << FORMAT(angleEnergy);
02404 if ( simParameters->mergeCrossterms ) {
02405 iout << FORMAT(dihedralEnergy+crosstermEnergy);
02406 } else {
02407 iout << FORMAT(dihedralEnergy);
02408 iout << FORMAT(crosstermEnergy);
02409 }
02410 iout << FORMAT(improperEnergy);
02411 iout << " ";
02412 iout << FORMAT(electEnergy+electEnergySlow);
02413 iout << FORMAT(ljEnergy);
02414 iout << FORMAT(boundaryEnergy);
02415 iout << FORMAT(miscEnergy);
02416 iout << FORMAT(kineticEnergy);
02417 iout << " ";
02418 iout << FORMAT(totalEnergy);
02419 iout << FORMAT(temperature);
02420 iout << FORMAT(potentialEnergy);
02421 // iout << FORMAT(flatEnergy);
02422 iout << FORMAT(smoothEnergy);
02423 iout << FORMAT(temp_avg/avg_count);
02424 if ( volume != 0. )
02425 {
02426 iout << " ";
02427 iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
02428 iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
02429 iout << FORMAT(volume);
02430 iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
02431 iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
02432 }
02433 if (simParameters->drudeOn) {
02434 iout << " ";
02435 iout << FORMAT(drudeComTemp);
02436 iout << FORMAT(drudeBondTemp);
02437 iout << FORMAT(drudeComTempAvg/avg_count);
02438 iout << FORMAT(drudeBondTempAvg/avg_count);
02439 }
02440 // Ported by JLai
02441 if (simParameters->goForcesOn) {
02442 iout << " ";
02443 iout << FORMAT(goNativeEnergy);
02444 iout << FORMAT(goNonnativeEnergy);
02445 //iout << FORMAT(relgoNativeEnergy);
02446 //iout << FORMAT(relgoNonnativeEnergy);
02447 iout << FORMAT(goTotalEnergy);
02448 //iout << FORMAT("not implemented");
02449 } // End of port -- JLai
02450 iout << "\n\n" << endi;
02451
02452 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
02453 char webout[80];
02454 sprintf(webout,"%d %d %d %d",(int)totalEnergy,
02455 (int)(potentialEnergy),
02456 (int)kineticEnergy,(int)temperature);
02457 CApplicationDepositNode0Data(webout);
02458 #endif
02459
02460 if (simParameters->pairInteractionOn) {
02461 iout << "PAIR INTERACTION:";
02462 iout << " STEP: " << step;
02463 iout << " VDW_FORCE: ";
02464 iout << FORMAT(pairVDWForce.x);
02465 iout << FORMAT(pairVDWForce.y);
02466 iout << FORMAT(pairVDWForce.z);
02467 iout << " ELECT_FORCE: ";
02468 iout << FORMAT(pairElectForce.x);
02469 iout << FORMAT(pairElectForce.y);
02470 iout << FORMAT(pairElectForce.z);
02471 iout << "\n" << endi;
02472 }
02473 drudeComTempAvg = 0;
02474 drudeBondTempAvg = 0;
02475 temp_avg = 0;
02476 pressure_avg = 0;
02477 groupPressure_avg = 0;
02478 avg_count = 0;
02479
02480 if ( fflush_count ) {
02481 --fflush_count;
02482 fflush(stdout);
02483 }
02484 }
|
|
|
Definition at line 928 of file Controller.C. References SimParameters::alchEquilSteps, SimParameters::alchFepOn, SimParameters::alchLambda, SimParameters::alchLambda2, SimParameters::alchTemp, BigReal, iout, and simParams. Referenced by integrate(). 00929 {
00930 if (simParams->alchFepOn) {
00931 const BigReal alchLambda = simParams->alchLambda;
00932 const BigReal alchLambda2 = simParams->alchLambda2;
00933 const BigReal alchTemp = simParams->alchTemp;
00934 const int alchEquilSteps = simParams->alchEquilSteps;
00935 iout << "FEP: RESETTING FOR NEW FEP WINDOW "
00936 << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
00937 << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
00938 << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
00939 << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
00940 << " K FOR FEP CALCULATION\n" << endi;
00941 }
00942 }
|
|
|
Definition at line 2021 of file Controller.C. References compareChecksums(), RequireReduction::item(), min_energy, min_f_dot_f, min_f_dot_v, min_huge_count, min_v_dot_v, Node::molecule, Node::Object(), printEnergies(), receivePressure(), reduction, REDUCTION_MIN_F_DOT_F, REDUCTION_MIN_F_DOT_V, REDUCTION_MIN_HUGE_COUNT, REDUCTION_MIN_V_DOT_V, and rescaleaccelMD(). 02021 {
02022
02023 rescaleaccelMD(step,1);
02024 receivePressure(step,1);
02025
02026 Node *node = Node::Object();
02027 Molecule *molecule = node->molecule;
02028 compareChecksums(step,1);
02029
02030 printEnergies(step,1);
02031
02032 min_energy = totalEnergy;
02033 min_f_dot_f = reduction->item(REDUCTION_MIN_F_DOT_F);
02034 min_f_dot_v = reduction->item(REDUCTION_MIN_F_DOT_V);
02035 min_v_dot_v = reduction->item(REDUCTION_MIN_V_DOT_V);
02036 min_huge_count = (int) (reduction->item(REDUCTION_MIN_HUGE_COUNT));
02037
02038 }
|
|
|
Definition at line 943 of file Controller.C. References SimParameters::alchEquilSteps, SimParameters::alchLambda, SimParameters::alchTemp, SimParameters::alchThermIntOn, BigReal, iout, and simParams. Referenced by integrate(). 00944 {
00945 if (simParams->alchThermIntOn) {
00946 const BigReal alchLambda = simParams->alchLambda;
00947 const BigReal alchTemp = simParams->alchTemp;
00948 const int alchEquilSteps = simParams->alchEquilSteps;
00949 iout << "TI: RESETTING FOR NEW WINDOW "
00950 << "LAMBDA SET TO " << alchLambda
00951 << "\nTI: WINDOW TO HAVE " << alchEquilSteps
00952 << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
00953 }
00954 }
|
|
|
Definition at line 1984 of file Controller.C. References fflush_count, SimParameters::firstTimestep, iout, iWARN(), memusage_MB(), SimParameters::N, SimParameters::outputEnergies, SimParameters::outputTiming, and simParams. Referenced by printEnergies(). 01984 {
01985
01986 if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
01987 {
01988 const double endWTime = CmiWallTimer() - firstWTime;
01989 const double endCTime = CmiTimer() - firstCTime;
01990
01991 // fflush about once per minute
01992 if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
01993
01994 const double elapsedW =
01995 (endWTime - startWTime) / simParams->outputTiming;
01996 const double elapsedC =
01997 (endCTime - startCTime) / simParams->outputTiming;
01998
01999 const double remainingW = elapsedW * (simParams->N - step);
02000 const double remainingW_hours = remainingW / 3600;
02001
02002 startWTime = endWTime;
02003 startCTime = endCTime;
02004
02005 #ifdef NAMD_CUDA
02006 if ( simParams->outputEnergies < 60 &&
02007 step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
02008 iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
02009 }
02010 #endif
02011 if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
02012 CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
02013 ", %g hours remaining, %f MB of memory in use.\n",
02014 step, endCTime, elapsedC, endWTime, elapsedW,
02015 remainingW_hours, memusage_MB());
02016 if ( fflush_count ) { --fflush_count; fflush(stdout); }
02017 }
02018 }
02019 }
|
|
|
Definition at line 957 of file Controller.C. References BigReal, iout, SimParameters::reassignFreq, SimParameters::reassignHold, SimParameters::reassignIncr, SimParameters::reassignTemp, and simParams. Referenced by integrate(). 00958 {
00959 const int reassignFreq = simParams->reassignFreq;
00960 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00961 BigReal newTemp = simParams->reassignTemp;
00962 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00963 if ( simParams->reassignIncr > 0.0 ) {
00964 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00965 newTemp = simParams->reassignHold;
00966 } else {
00967 if ( newTemp < simParams->reassignHold )
00968 newTemp = simParams->reassignHold;
00969 }
00970 iout << "REASSIGNING VELOCITIES AT STEP " << step
00971 << " TO " << newTemp << " KELVIN.\n" << endi;
00972 }
00973 }
|
|
|
Definition at line 2840 of file Controller.C. References fflush_count, LdbCoordinator::getNumStepsToRun(), ldbSteps, Node::Object(), LdbCoordinator::Object(), Node::outputPatchComputeMaps(), and LdbCoordinator::rebalance(). Referenced by integrate(). 02841 {
02842 if ( ! ldbSteps ) {
02843 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02844 }
02845 if ( ! --ldbSteps ) {
02846 startBenchTime -= CmiWallTimer();
02847 Node::Object()->outputPatchComputeMaps("before_ldb", step);
02848 LdbCoordinator::Object()->rebalance(this);
02849 startBenchTime += CmiWallTimer();
02850 fflush_count = 3;
02851 }
02852 }
|
|
||||||||||||
|
Definition at line 1010 of file Controller.C. References SimParameters::accelMDDebugOn, SimParameters::accelMDOn, SimParameters::accelMDOutFreq, AVGXY, BigReal, SimParameters::comMove, controlNumDegFreedom, controlPressure, controlPressure_nbond, controlPressure_normal, controlPressure_slow, Tensor::diagonal(), drudeBondTemp, drudeComTemp, SimParameters::drudeOn, SimParameters::fixCellDims, SimParameters::fixCellDimX, SimParameters::fixCellDimY, SimParameters::fixCellDimZ, SimParameters::fixedAtomsOn, GET_TENSOR, GET_VECTOR, groupPressure, groupPressure_nbond, groupPressure_normal, groupPressure_slow, Tensor::identity(), iINFO(), iout, RequireReduction::item(), kineticEnergy, kineticEnergyCentered, kineticEnergyHalfstep, SimParameters::langevinOn, NamdState::lattice, SimParameters::LJcorrection, Node::molecule, Molecule::num_deg_freedom(), Molecule::num_fixed_atoms(), Molecule::num_fixed_groups(), Molecule::num_group_deg_freedom(), Molecule::numAtoms, Molecule::numConstraints, numDegFreedom, Molecule::numDrudeAtoms, Molecule::numFepInitial, Molecule::numFixedAtoms, Molecule::numFixedGroups, Molecule::numFixedRigidBonds, Molecule::numHydrogenGroups, Molecule::numLonepairs, Molecule::numRigidBonds, Node::Object(), Lattice::origin(), outer(), SimParameters::pairInteractionOn, pressure, pressure_nbond, pressure_normal, pressure_slow, PRESSUREFACTOR, reduction, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, RequireReduction::require(), Node::simParameters, simParams, state, Molecule::tail_corr_virial, temperature, SimParameters::useConstantRatio, SimParameters::useFlexibleCell, SimParameters::useGroupPressure, virial_amd, and Lattice::volume(). Referenced by integrate(), and printMinimizeEnergies(). 01011 {
01012 Node *node = Node::Object();
01013 Molecule *molecule = node->molecule;
01014 SimParameters *simParameters = node->simParameters;
01015 Lattice &lattice = state->lattice;
01016
01017 reduction->require();
01018
01019 Tensor virial;
01020 Tensor virial_normal;
01021 Tensor virial_nbond;
01022 Tensor virial_slow;
01023 Tensor pressure_amd;
01024 #ifdef ALTVIRIAL
01025 Tensor altVirial_normal;
01026 Tensor altVirial_nbond;
01027 Tensor altVirial_slow;
01028 #endif
01029 Tensor intVirial;
01030 Tensor intVirial_normal;
01031 Tensor intVirial_nbond;
01032 Tensor intVirial_slow;
01033 Vector extForce_normal;
01034 Vector extForce_nbond;
01035 Vector extForce_slow;
01036 BigReal volume;
01037
01038 #if 1
01039 numDegFreedom = molecule->num_deg_freedom();
01040 int numGroupDegFreedom = molecule->num_group_deg_freedom();
01041 int numFixedGroups = molecule->num_fixed_groups();
01042 int numFixedAtoms = molecule->num_fixed_atoms();
01043 #endif
01044 #if 0
01045 int numAtoms = molecule->numAtoms;
01046 numDegFreedom = 3 * numAtoms;
01047 int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
01048 int numFixedAtoms =
01049 ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
01050 int numLonepairs = molecule->numLonepairs;
01051 int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
01052 if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
01053 if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
01054 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
01055 if ( ! ( numFixedAtoms || molecule->numConstraints
01056 || simParameters->comMove || simParameters->langevinOn ) ) {
01057 numDegFreedom -= 3;
01058 numGroupDegFreedom -= 3;
01059 }
01060 if (simParameters->pairInteractionOn) {
01061 // this doesn't attempt to deal with fixed atoms or constraints
01062 numDegFreedom = 3 * molecule->numFepInitial;
01063 }
01064 int numRigidBonds = molecule->numRigidBonds;
01065 int numFixedRigidBonds =
01066 ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
01067
01068 // numLonepairs is subtracted here because all lonepairs have a rigid bond
01069 // to oxygen, but all of the LP degrees of freedom are dealt with above
01070 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
01071 #endif
01072
01073 kineticEnergyHalfstep = reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY);
01074 kineticEnergyCentered = reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY);
01075
01076 BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
01077 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY);
01078 BigReal groupKineticEnergyCentered = kineticEnergyCentered -
01079 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY);
01080
01081 BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
01082 / ( numDegFreedom * BOLTZMANN );
01083 BigReal atomTempCentered = 2.0 * kineticEnergyCentered
01084 / ( numDegFreedom * BOLTZMANN );
01085 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
01086 / ( numGroupDegFreedom * BOLTZMANN );
01087 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
01088 / ( numGroupDegFreedom * BOLTZMANN );
01089
01090 /* test code for comparing different temperatures
01091 iout << "TEMPTEST: " << step << " " <<
01092 atomTempHalfstep << " " <<
01093 atomTempCentered << " " <<
01094 groupTempHalfstep << " " <<
01095 groupTempCentered << "\n" << endi;
01096 iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
01097 numGroupDegFreedom << "\n" << endi;
01098 */
01099
01100 GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
01101 GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
01102 GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
01103
01104 #ifdef ALTVIRIAL
01105 GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
01106 GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
01107 GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
01108 #endif
01109
01110 GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
01111 GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
01112 GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
01113
01114 GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
01115 GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
01116 GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
01117 Vector extPosition = lattice.origin();
01118 virial_normal -= outer(extForce_normal,extPosition);
01119 virial_nbond -= outer(extForce_nbond,extPosition);
01120 virial_slow -= outer(extForce_slow,extPosition);
01121
01122
01123 kineticEnergy = kineticEnergyCentered;
01124 temperature = 2.0 * kineticEnergyCentered / ( numDegFreedom * BOLTZMANN );
01125
01126 if (simParameters->drudeOn) {
01127 BigReal drudeComKE
01128 = reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY);
01129 BigReal drudeBondKE
01130 = reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY);
01131 int g_bond = 3 * molecule->numDrudeAtoms;
01132 int g_com = numDegFreedom - g_bond;
01133 drudeComTemp = 2.0 * drudeComKE / (g_com * BOLTZMANN);
01134 drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
01135 }
01136
01137 if ( (volume=lattice.volume()) != 0. )
01138 {
01139
01140 if (simParameters->LJcorrection && volume) {
01141 // Apply tail correction to pressure
01142 //printf("Volume is %f\n", volume);
01143 //printf("Applying tail correction of %f to virial\n", molecule->tail_corr_virial / volume);
01144 virial_normal += Tensor::identity(molecule->tail_corr_virial / volume);
01145 }
01146
01147 // kinetic energy component included in virials
01148 pressure_normal = virial_normal / volume;
01149 groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
01150
01151 if (simParameters->accelMDOn) {
01152 pressure_amd = virial_amd / volume;
01153 pressure_normal += pressure_amd;
01154 groupPressure_normal += pressure_amd;
01155 }
01156
01157 if ( minimize || ! ( step % nbondFreq ) )
01158 {
01159 pressure_nbond = virial_nbond / volume;
01160 groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
01161 }
01162
01163 if ( minimize || ! ( step % slowFreq ) )
01164 {
01165 pressure_slow = virial_slow / volume;
01166 groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
01167 }
01168
01169 /*
01170 iout << "VIRIALS: " << virial_normal << " " << virial_nbond << " " <<
01171 virial_slow << " " << ( virial_normal - intVirial_normal ) << " " <<
01172 ( virial_nbond - intVirial_nbond ) << " " <<
01173 ( virial_slow - intVirial_slow ) << "\n";
01174 */
01175
01176 pressure = pressure_normal + pressure_nbond + pressure_slow;
01177 groupPressure = groupPressure_normal + groupPressure_nbond +
01178 groupPressure_slow;
01179 }
01180 else
01181 {
01182 pressure = Tensor();
01183 groupPressure = Tensor();
01184 }
01185
01186 if ( simParameters->useGroupPressure )
01187 {
01188 controlPressure_normal = groupPressure_normal;
01189 controlPressure_nbond = groupPressure_nbond;
01190 controlPressure_slow = groupPressure_slow;
01191 controlPressure = groupPressure;
01192 controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
01193 if ( ! ( numFixedAtoms || molecule->numConstraints
01194 || simParameters->comMove || simParameters->langevinOn ) ) {
01195 controlNumDegFreedom -= 1;
01196 }
01197 }
01198 else
01199 {
01200 controlPressure_normal = pressure_normal;
01201 controlPressure_nbond = pressure_nbond;
01202 controlPressure_slow = pressure_slow;
01203 controlPressure = pressure;
01204 controlNumDegFreedom = numDegFreedom / 3;
01205 }
01206
01207 if (simParameters->fixCellDims) {
01208 if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
01209 if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
01210 if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
01211 }
01212
01213 if ( simParameters->useFlexibleCell ) {
01214 // use symmetric pressure to control rotation
01215 // controlPressure_normal = symmetric(controlPressure_normal);
01216 // controlPressure_nbond = symmetric(controlPressure_nbond);
01217 // controlPressure_slow = symmetric(controlPressure_slow);
01218 // controlPressure = symmetric(controlPressure);
01219 // only use on-diagonal components for now
01220 controlPressure_normal = Tensor::diagonal(diagonal(controlPressure_normal));
01221 controlPressure_nbond = Tensor::diagonal(diagonal(controlPressure_nbond));
01222 controlPressure_slow = Tensor::diagonal(diagonal(controlPressure_slow));
01223 controlPressure = Tensor::diagonal(diagonal(controlPressure));
01224 if ( simParameters->useConstantRatio ) {
01225 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
01226 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
01227 AVGXY(controlPressure_normal);
01228 AVGXY(controlPressure_nbond);
01229 AVGXY(controlPressure_slow);
01230 AVGXY(controlPressure);
01231 #undef AVGXY
01232 }
01233 } else {
01234 controlPressure_normal =
01235 Tensor::identity(trace(controlPressure_normal)/3.);
01236 controlPressure_nbond =
01237 Tensor::identity(trace(controlPressure_nbond)/3.);
01238 controlPressure_slow =
01239 Tensor::identity(trace(controlPressure_slow)/3.);
01240 controlPressure =
01241 Tensor::identity(trace(controlPressure)/3.);
01242 }
01243
01244 #ifdef DEBUG_PRESSURE
01245 iout << iINFO << "Control pressure = " << controlPressure <<
01246 " = " << controlPressure_normal << " + " <<
01247 controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
01248 #endif
01249 if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
01250 iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
01251 << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
01252 << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
01253 << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
01254 << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
01255 << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
01256 << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
01257 << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
01258 << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
01259 << endi;
01260 }
01261 }
|
|
||||||||||||
|
Definition at line 1263 of file Controller.C. References SimParameters::accelMDalpha, SimParameters::accelMDDebugOn, SimParameters::accelMDdihe, SimParameters::accelMDdual, accelMDdVAverage, SimParameters::accelMDE, SimParameters::accelMDFirstStep, SimParameters::accelMDLastStep, SimParameters::accelMDOn, SimParameters::accelMDOutFreq, ControllerBroadcasts::accelMDRescaleFactor, SimParameters::accelMDTalpha, SimParameters::accelMDTE, amd_reduction, BigReal, broadcast, SimParameters::firstTimestep, GET_TENSOR, iout, RequireReduction::item(), iWARN(), NamdState::lattice, Node::molecule, Node::Object(), PRESSUREFACTOR, SimpleBroadcastObject< T >::publish(), REDUCTION_ANGLE_ENERGY, REDUCTION_BOND_ENERGY, REDUCTION_CROSSTERM_ENERGY, REDUCTION_DIHEDRAL_ENERGY, REDUCTION_ELECT_ENERGY, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_IMPROPER_ENERGY, REDUCTION_LJ_ENERGY, REDUCTION_VIRIAL_AMD_DIHE, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, RequireReduction::require(), simParams, state, virial_amd, and Lattice::volume(). Referenced by integrate(), and printMinimizeEnergies(). 01264 {
01265 if ( !simParams->accelMDOn ) return;
01266
01267 amd_reduction->require();
01268
01269 if (step == simParams->firstTimestep) accelMDdVAverage = 0;
01270 // if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
01271 if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
01272
01273 Node *node = Node::Object();
01274 Molecule *molecule = node->molecule;
01275 Lattice &lattice = state->lattice;
01276
01277 const BigReal accelMDE = simParams->accelMDE;
01278 const BigReal accelMDalpha = simParams->accelMDalpha;
01279 const BigReal accelMDTE = simParams->accelMDTE;
01280 const BigReal accelMDTalpha = simParams->accelMDTalpha;
01281 const int accelMDOutFreq = simParams->accelMDOutFreq;
01282
01283 BigReal bondEnergy;
01284 BigReal angleEnergy;
01285 BigReal dihedralEnergy;
01286 BigReal improperEnergy;
01287 BigReal crosstermEnergy;
01288 BigReal amd_electEnergy;
01289 BigReal amd_ljEnergy;
01290 BigReal amd_electEnergySlow;
01291 BigReal potentialEnergy;
01292 BigReal factor_dihe = 1;
01293 BigReal factor_tot = 1;
01294 BigReal testV;
01295 BigReal dV = 0;
01296 Vector accelMDfactor;
01297 Tensor vir;
01298 Tensor vir_dihe;
01299 Tensor vir_normal;
01300 Tensor vir_nbond;
01301 Tensor vir_slow;
01302
01303 bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
01304 angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
01305 dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
01306 improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
01307 crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
01308
01309 GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
01310 GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
01311 GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
01312 GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
01313
01314 if ( !( step % nbondFreq ) ) {
01315 amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
01316 amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
01317 } else {
01318 amd_electEnergy = electEnergy;
01319 amd_ljEnergy = ljEnergy;
01320 }
01321
01322 if ( !( step % slowFreq ) ) {
01323 amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
01324 } else {
01325 amd_electEnergySlow = electEnergySlow;
01326 }
01327
01328 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy + improperEnergy
01329 + crosstermEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy;
01330
01331 if (simParams->accelMDdihe) {
01332
01333 testV = dihedralEnergy + crosstermEnergy;
01334 if ( testV < accelMDE ) {
01335 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
01336 factor_dihe *= factor_dihe;
01337 vir = vir_dihe * (factor_dihe - 1.0);
01338 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01339 accelMDdVAverage += dV;
01340 }
01341
01342 } else if (simParams->accelMDdual) {
01343
01344 testV = dihedralEnergy + crosstermEnergy;
01345 if ( testV < accelMDE ) {
01346 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
01347 factor_dihe *= factor_dihe;
01348 vir = vir_dihe * (factor_dihe - 1.0) ;
01349 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01350 }
01351
01352 testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
01353 if ( testV < accelMDTE ) {
01354 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
01355 factor_tot *= factor_tot;
01356 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
01357 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
01358 }
01359 accelMDdVAverage += dV;
01360
01361 } else {
01362
01363 testV = potentialEnergy;
01364 if ( testV < accelMDE ) {
01365 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
01366 factor_tot *= factor_tot;
01367 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
01368 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
01369 accelMDdVAverage += dV;
01370 }
01371 }
01372
01373 accelMDfactor[0]=factor_dihe;
01374 accelMDfactor[1]=factor_tot;
01375 accelMDfactor[2]=1;
01376 broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
01377 virial_amd = vir;
01378
01379 if ( factor_tot < 0.001 ) {
01380 iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
01381 << "\n" << endi;
01382 }
01383
01384 if ( ! (step % accelMDOutFreq) ) {
01385 if ( !(step == simParams->firstTimestep) ) {
01386 accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
01387 }
01388 iout << "ACCELERATED MD: STEP " << step
01389 << " dV " << dV
01390 << " dVAVG " << accelMDdVAverage
01391 << " BOND " << bondEnergy
01392 << " ANGLE " << angleEnergy
01393 << " DIHED " << dihedralEnergy+crosstermEnergy
01394 << " IMPRP " << improperEnergy
01395 << " ELECT " << amd_electEnergy+amd_electEnergySlow
01396 << " VDW " << amd_ljEnergy
01397 << " POTENTIAL " << potentialEnergy << "\n"
01398 << endi;
01399
01400 accelMDdVAverage = 0;
01401
01402 if (simParams->accelMDDebugOn) {
01403 BigReal volume;
01404 Tensor p_normal;
01405 Tensor p_nbond;
01406 Tensor p_slow;
01407 Tensor p;
01408 if ( (volume=lattice.volume()) != 0. ) {
01409 p_normal = vir_normal/volume;
01410 p_nbond = vir_nbond/volume;
01411 p_slow = vir_slow/volume;
01412 p = vir/volume;
01413 }
01414 iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
01415 << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
01416 << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
01417 << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
01418 << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
01419 << endi;
01420 }
01421 }
01422 }
|
|
|
|
Definition at line 2871 of file Controller.C. References awaken(), broadcast, SimpleBroadcastObject< T >::publish(), and ControllerBroadcasts::traceBarrier. Referenced by Node::resumeAfterTraceBarrier(). 02871 {
02872 broadcast->traceBarrier.publish(step,1);
02873 CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02874 awaken();
02875 }
|
|
|
Definition at line 220 of file Controller.C. References awaken(), CTRL_STK_SZ, and DebugM. Referenced by NamdState::runController(). 00221 {
00222 // create a Thread and invoke it
00223 DebugM(4, "Starting thread in controller on this=" << this << "\n");
00224 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
00225 CthSetStrategyDefault(thread);
00226 #if CMK_BLUEGENE_CHARM
00227 BgAttach(thread);
00228 #endif
00229 awaken();
00230 }
|
|
|
Definition at line 975 of file Controller.C. References BigReal, broadcast, SimpleBroadcastObject< T >::publish(), simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, SimParameters::tCoupleTemp, and temperature. Referenced by integrate(). 00976 {
00977 if ( simParams->tCoupleOn )
00978 {
00979 const BigReal tCoupleTemp = simParams->tCoupleTemp;
00980 BigReal coefficient = 1.;
00981 if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
00982 broadcast->tcoupleCoefficient.publish(step,coefficient);
00983 }
00984 }
|
|
|
Definition at line 2892 of file Controller.C. References BackEnd::awaken(). Referenced by algorithm(). 02892 {
02893 BackEnd::awaken();
02894 CthFree(thread);
02895 CthSuspend();
02896 }
|
|
||||||||||||
|
Definition at line 2864 of file Controller.C. References ControllerBroadcasts::traceBarrier. Referenced by integrate(). 02864 {
02865 CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
02866 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
02867 nd.traceBarrier(turnOnTrace, step);
02868 CthSuspend();
02869 }
|
|
||||||||||||
|
Definition at line 2499 of file Controller.C. References Lattice::a(), Lattice::a_p(), Lattice::b(), Lattice::b_p(), Lattice::c(), Lattice::c_p(), langevinPiston_strainRate, SimParameters::langevinPistonOn, NamdState::lattice, Lattice::origin(), simParams, state, Vector::x, Vector::y, and Vector::z. Referenced by outputExtendedSystem(). 02499 {
02500 Lattice &lattice = state->lattice;
02501 file << step;
02502 if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
02503 if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
02504 if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
02505 file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
02506 if ( simParams->langevinPistonOn ) {
02507 Vector strainRate = diagonal(langevinPiston_strainRate);
02508 Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
02509 file << " " << strainRate.x;
02510 file << " " << strainRate.y;
02511 file << " " << strainRate.z;
02512 file << " " << strainRate2.x;
02513 file << " " << strainRate2.y;
02514 file << " " << strainRate2.z;
02515 }
02516 file << std::endl;
02517 }
|
|
|
Definition at line 2486 of file Controller.C. References Lattice::a_p(), Lattice::b_p(), Lattice::c_p(), SimParameters::langevinPistonOn, NamdState::lattice, simParams, and state. Referenced by outputExtendedSystem(). 02486 {
02487 Lattice &lattice = state->lattice;
02488 file << "#$LABELS step";
02489 if ( lattice.a_p() ) file << " a_x a_y a_z";
02490 if ( lattice.b_p() ) file << " b_x b_y b_z";
02491 if ( lattice.c_p() ) file << " c_x c_y c_z";
02492 file << " o_x o_y o_z";
02493 if ( simParams->langevinPistonOn ) {
02494 file << " s_x s_y s_z s_u s_v s_w";
02495 }
02496 file << std::endl;
02497 }
|
|
||||||||||||
|
Definition at line 2684 of file Controller.C. References SimParameters::alchEnsembleAvg, SimParameters::alchTemp, BigReal, BOLTZMANN, dG, electEnergy, electEnergy_f, exp_dE_ByRT, fepFile, FepNo, FEPTITLE(), SimParameters::firstTimestep, FORMAT(), ljEnergy_f, net_dE, simParams, and temperature. Referenced by outputFepEnergy(). 02684 {
02685 BigReal eeng = electEnergy+electEnergySlow;
02686 BigReal eeng_f = electEnergy_f + electEnergySlow_f;
02687 BigReal dE = eeng_f + ljEnergy_f - eeng - ljEnergy;
02688 BigReal RT = BOLTZMANN * simParams->alchTemp;
02689 const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
02690 const int stepInRun = step - simParams->firstTimestep;
02691 if(stepInRun){
02692 fepFile << FEPTITLE(step);
02693 fepFile << FORMAT(eeng);
02694 fepFile << FORMAT(eeng_f);
02695 fepFile << FORMAT(ljEnergy);
02696 fepFile << FORMAT(ljEnergy_f);
02697 fepFile << FORMAT(dE);
02698 if(alchEnsembleAvg){
02699 BigReal dE_avg = net_dE/FepNo;
02700 fepFile << FORMAT(dE_avg);
02701 }
02702 fepFile << FORMAT(temperature);
02703 if(alchEnsembleAvg){
02704 dG = -(RT * log(exp_dE_ByRT/FepNo));
02705 fepFile << FORMAT(dG);
02706 }
02707 fepFile << std::endl;
02708 }
02709 }
|
|
||||||||||||
|
Definition at line 2711 of file Controller.C. References FORMAT(), net_dEdl_elec_1, net_dEdl_elec_2, net_dEdl_lj_1, net_dEdl_lj_2, recent_dEdl_elec_1, recent_dEdl_elec_2, recent_dEdl_lj_1, recent_dEdl_lj_2, recent_TiNo, tiFile, TiNo, and TITITLE(). Referenced by outputTiEnergy(). 02711 {
02712 tiFile << TITITLE(step);
02713 tiFile << FORMAT(recent_dEdl_elec_1 / recent_TiNo);
02714 tiFile << FORMAT(net_dEdl_elec_1/TiNo);
02715 tiFile << FORMAT(recent_dEdl_lj_1 / recent_TiNo);
02716 tiFile << FORMAT(net_dEdl_lj_1/TiNo);
02717 tiFile << FORMAT(recent_dEdl_elec_2 / recent_TiNo);
02718 tiFile << FORMAT(net_dEdl_elec_2/TiNo);
02719 tiFile << FORMAT(recent_dEdl_lj_2 / recent_TiNo);
02720 tiFile << FORMAT(net_dEdl_lj_2/TiNo);
02721 tiFile << std::endl;
02722 }
|
|
|
Definition at line 42 of file Controller.h. |
|
|
Definition at line 205 of file Controller.h. Referenced by rescaleaccelMD(). |
|
|
Definition at line 228 of file Controller.h. Referenced by adaptTempInit(). |
|
|
Definition at line 222 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 221 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 217 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 223 of file Controller.h. Referenced by adaptTempUpdate(). |
|
|
Definition at line 224 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 226 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 225 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 227 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 219 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 220 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 230 of file Controller.h. Referenced by adaptTempInit(), and adaptTempUpdate(). |
|
|
Definition at line 229 of file Controller.h. Referenced by adaptTempInit(). |
|
|
Definition at line 214 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 212 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 211 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 216 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 215 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 213 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 231 of file Controller.h. Referenced by adaptTempInit(), and adaptTempWriteRestart(). |
|
|
Definition at line 218 of file Controller.h. Referenced by adaptTempInit(), adaptTempUpdate(), and adaptTempWriteRestart(). |
|
|
Definition at line 168 of file Controller.h. Referenced by Controller(), and rescaleaccelMD(). |
|
|
Definition at line 64 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 144 of file Controller.h. Referenced by algorithm(), berendsenPressure(), and Controller(). |
|
|
Definition at line 145 of file Controller.h. Referenced by algorithm(), berendsenPressure(), and Controller(). |
|
|
Definition at line 180 of file Controller.h. Referenced by adaptTempUpdate(), algorithm(), berendsenPressure(), Controller(), correctMomentum(), cycleBarrier(), langevinPiston1(), minimize(), rescaleaccelMD(), rescaleVelocities(), resumeAfterTraceBarrier(), and tcoupleVelocities(). |
|
|
Definition at line 199 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 200 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 198 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 197 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 201 of file Controller.h. Referenced by algorithm(). |
|
|
Definition at line 196 of file Controller.h. Referenced by algorithm(), and Controller(). |
|
|
Definition at line 178 of file Controller.h. Referenced by enqueueCollections(). |
|
|
Definition at line 69 of file Controller.h. Referenced by compareChecksums(). |
|
|
Definition at line 134 of file Controller.h. Referenced by langevinPiston1(), langevinPiston2(), and receivePressure(). |
|
|
Definition at line 135 of file Controller.h. Referenced by langevinPiston1(), langevinPiston2(), and receivePressure(). |
|
|
Definition at line 57 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 56 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 58 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 95 of file Controller.h. Referenced by writeFepEnergyData(). |
|
|
Definition at line 122 of file Controller.h. Referenced by Controller(), printEnergies(), and receivePressure(). |
|
|
Definition at line 124 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 121 of file Controller.h. Referenced by Controller(), printEnergies(), and receivePressure(). |
|
|
Definition at line 123 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 83 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 90 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 101 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 104 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 111 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 112 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 84 of file Controller.h. Referenced by outputFepEnergy(), and printEnergies(). |
|
|
Definition at line 91 of file Controller.h. Referenced by outputFepEnergy(), and printEnergies(). |
|
|
Definition at line 102 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 105 of file Controller.h. Referenced by outputTiEnergy(), and printEnergies(). |
|
|
Definition at line 93 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 187 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 96 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 98 of file Controller.h. Referenced by outputFepEnergy(). |
|
|
Definition at line 152 of file Controller.h. Referenced by printTiming(), and rebalanceLoad(). |
|
|
Definition at line 86 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 87 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 88 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 133 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 63 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 54 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 53 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 55 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 66 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 126 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 128 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 127 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 148 of file Controller.h. Referenced by algorithm(), Controller(), and writeExtendedSystemData(). |
|
|
Definition at line 150 of file Controller.h. Referenced by rebalanceLoad(). |
|
|
Definition at line 85 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 92 of file Controller.h. Referenced by outputFepEnergy(), printEnergies(), and writeFepEnergyData(). |
|
|
Definition at line 103 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 106 of file Controller.h. Referenced by printEnergies(). |
|
|
Definition at line 70 of file Controller.h. Referenced by compareChecksums(), and printEnergies(). |
|
|
Definition at line 74 of file Controller.h. Referenced by printMinimizeEnergies(). |
|
|
Definition at line 75 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 76 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 78 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 77 of file Controller.h. Referenced by minimize(), and printMinimizeEnergies(). |
|
|
Definition at line 59 of file Controller.h. Referenced by integrate(), langevinPiston1(), langevinPiston2(), and minimize(). |
|
|
Definition at line 94 of file Controller.h. Referenced by outputFepEnergy(), and writeFepEnergyData(). |
|
|
Definition at line 107 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 108 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 109 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 110 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 81 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 71 of file Controller.h. Referenced by compareChecksums(), and printEnergies(). |
|
|
Definition at line 171 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 173 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 172 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 132 of file Controller.h. Referenced by printEnergies(), and receivePressure(). |
|
|
Definition at line 62 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 50 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 49 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 51 of file Controller.h. Referenced by receivePressure(). |
|
|
Definition at line 65 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 176 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 175 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 174 of file Controller.h. Referenced by Controller(). |
|
|
Definition at line 164 of file Controller.h. Referenced by adaptTempUpdate(), Controller(), langevinPiston1(), and langevinPiston2(). |
|
|
Definition at line 114 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 115 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 116 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 117 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 118 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 167 of file Controller.h. Referenced by compareChecksums(), Controller(), correctMomentum(), printEnergies(), printMinimizeEnergies(), and receivePressure(). |
|
|
Definition at line 140 of file Controller.h. Referenced by Controller(), and rescaleVelocities(). |
|
|
Definition at line 139 of file Controller.h. Referenced by Controller(), and rescaleVelocities(). |
|
|
|
Definition at line 60 of file Controller.h. Referenced by correctMomentum(), integrate(), langevinPiston1(), langevinPiston2(), and minimize(). |
|
|
Definition at line 130 of file Controller.h. Referenced by algorithm(), Controller(), and printEnergies(). |
|
|
Definition at line 131 of file Controller.h. |
|
|
Definition at line 166 of file Controller.h. Referenced by algorithm(), berendsenPressure(), enqueueCollections(), langevinPiston1(), langevinPiston2(), printDynamicsEnergies(), printEnergies(), receivePressure(), rescaleaccelMD(), writeExtendedSystemData(), and writeExtendedSystemLabels(). |
|
|
Definition at line 67 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 61 of file Controller.h. Referenced by Controller(), and printEnergies(). |
|
|
Definition at line 129 of file Controller.h. Referenced by printEnergies(), receivePressure(), tcoupleVelocities(), and writeFepEnergyData(). |
|
|
Definition at line 191 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 113 of file Controller.h. Referenced by outputTiEnergy(), and writeTiEnergyData(). |
|
|
Definition at line 82 of file Controller.h. Referenced by adaptTempUpdate(), and printEnergies(). |
|
|
Definition at line 52 of file Controller.h. Referenced by receivePressure(), and rescaleaccelMD(). |
|
|
Definition at line 181 of file Controller.h. Referenced by outputExtendedSystem(). |
1.3.9.1