version 1.1320 | version 1.1321 |
---|
| |
/***************************************************************************** | /***************************************************************************** |
* $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $ | * $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $ |
* $Author: jim $ | * $Author: jim $ |
* $Date: 2016/09/05 18:25:42 $ | * $Date: 2016/09/14 20:22:03 $ |
* $Revision: 1.1320 $ | * $Revision: 1.1321 $ |
*****************************************************************************/ | *****************************************************************************/ |
| |
#include "InfoStream.h" | #include "InfoStream.h" |
| |
//Modifications for alchemical fep | //Modifications for alchemical fep |
void Controller::printFepMessage(int step) | void Controller::printFepMessage(int step) |
{ | { |
if (simParams->alchFepOn) { | if (simParams->alchOn && simParams->alchFepOn |
| && !simParams->alchLambdaFreq) { |
const BigReal alchLambda = simParams->alchLambda; | const BigReal alchLambda = simParams->alchLambda; |
const BigReal alchLambda2 = simParams->alchLambda2; | const BigReal alchLambda2 = simParams->alchLambda2; |
const BigReal alchTemp = simParams->alchTemp; | const BigReal alchTemp = simParams->alchTemp; |
| |
} | } |
void Controller::printTiMessage(int step) | void Controller::printTiMessage(int step) |
{ | { |
if (simParams->alchThermIntOn) { | if (simParams->alchOn && simParams->alchThermIntOn |
| && !simParams->alchLambdaFreq) { |
const BigReal alchLambda = simParams->alchLambda; | const BigReal alchLambda = simParams->alchLambda; |
const BigReal alchTemp = simParams->alchTemp; | |
const int alchEquilSteps = simParams->alchEquilSteps; | const int alchEquilSteps = simParams->alchEquilSteps; |
iout << "TI: RESETTING FOR NEW WINDOW " | iout << "TI: RESETTING FOR NEW WINDOW " |
<< "LAMBDA SET TO " << alchLambda | << "LAMBDA SET TO " << alchLambda |
| |
| |
// Some basic correctness checking | // Some basic correctness checking |
BigReal checksum, checksum_b; | BigReal checksum, checksum_b; |
| char errmsg[256]; |
| |
checksum = reduction->item(REDUCTION_ATOM_CHECKSUM); | checksum = reduction->item(REDUCTION_ATOM_CHECKSUM); |
if ( ((int)checksum) != molecule->numAtoms ) { | if ( ((int)checksum) != molecule->numAtoms ) { |
char errmsg[256]; | |
sprintf(errmsg, "Bad global atom count! (%d vs %d)\n", | sprintf(errmsg, "Bad global atom count! (%d vs %d)\n", |
(int)checksum, molecule->numAtoms); | (int)checksum, molecule->numAtoms); |
NAMD_bug(errmsg); | NAMD_bug(errmsg); |
| |
| |
checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM); | checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) { |
| sprintf(errmsg, "Bad global bond count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcBonds); |
if ( forgiving && (((int)checksum) < molecule->numCalcBonds) ) | if ( forgiving && (((int)checksum) < molecule->numCalcBonds) ) |
iout << iWARN << "Bad global bond count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global bond count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM); | checksum = reduction->item(REDUCTION_ANGLE_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) { |
| sprintf(errmsg, "Bad global angle count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcAngles); |
if ( forgiving && (((int)checksum) < molecule->numCalcAngles) ) | if ( forgiving && (((int)checksum) < molecule->numCalcAngles) ) |
iout << iWARN << "Bad global angle count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global angle count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM); | checksum = reduction->item(REDUCTION_DIHEDRAL_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) { |
| sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcDihedrals); |
if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) ) | if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) ) |
iout << iWARN << "Bad global dihedral count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global dihedral count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM); | checksum = reduction->item(REDUCTION_IMPROPER_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) { |
| sprintf(errmsg, "Bad global improper count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcImpropers); |
if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) ) | if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) ) |
iout << iWARN << "Bad global improper count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global improper count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_THOLE_CHECKSUM); | checksum = reduction->item(REDUCTION_THOLE_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) { |
| sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcTholes); |
if ( forgiving && (((int)checksum) < molecule->numCalcTholes) ) | if ( forgiving && (((int)checksum) < molecule->numCalcTholes) ) |
iout << iWARN << "Bad global Thole count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global Thole count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_ANISO_CHECKSUM); | checksum = reduction->item(REDUCTION_ANISO_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) { |
| sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcAnisos); |
if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) ) | if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) ) |
iout << iWARN << "Bad global Aniso count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global Aniso count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM); | checksum = reduction->item(REDUCTION_CROSSTERM_CHECKSUM); |
if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) { | if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) { |
| sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n", |
| (int)checksum, molecule->numCalcCrossterms); |
if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) ) | if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) ) |
iout << iWARN << "Bad global crossterm count!\n" << endi; | iout << iWARN << errmsg << endi; |
else NAMD_bug("Bad global crossterm count!\n"); | else NAMD_bug(errmsg); |
} | } |
| |
checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM); | checksum = reduction->item(REDUCTION_EXCLUSION_CHECKSUM); |
| |
} else if (simParameters->alchFepOn) { | } else if (simParameters->alchFepOn) { |
CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy + | CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy + |
improperEnergy + bondedEnergyDiff_f); | improperEnergy + bondedEnergyDiff_f); |
CALLBACKLIST("ELEC2", electEnergy_f); | CALLBACKLIST("ELEC2", electEnergy_f + electEnergySlow_f); |
CALLBACKLIST("VDW2", ljEnergy_f); | CALLBACKLIST("VDW2", ljEnergy_f); |
} | } |
} | } |
| |
//iout << FORMAT("GOAVG"); | //iout << FORMAT("GOAVG"); |
} | } |
// End of port -- JLai | // End of port -- JLai |
| |
if (simParameters->alchOn && simParameters->alchThermIntOn) { | if (simParameters->alchOn && simParameters->alchThermIntOn) { |
iout << " "; | iout << "\n#TITITLE: TS"; |
iout << FORMAT("BOND1"); | iout << FORMAT("BOND1"); |
iout << FORMAT("ELECT1"); | iout << FORMAT("ELECT1"); |
iout << FORMAT("VDW1"); | iout << FORMAT("VDW1"); |
iout << FORMAT("BOND2"); | iout << FORMAT("BOND2"); |
iout << FORMAT("ELECT2"); | |
iout << " "; | iout << " "; |
| iout << FORMAT("ELECT2"); |
iout << FORMAT("VDW2"); | iout << FORMAT("VDW2"); |
if (simParameters->alchLambdaFreq > 0) { | if (simParameters->alchLambdaFreq > 0) { |
iout << FORMAT("LAMBDA"); | iout << FORMAT("LAMBDA"); |
| |
} | } |
| |
iout << "\n\n" << endi; | iout << "\n\n" << endi; |
| |
} | } |
| |
// N.B. HP's aCC compiler merges FORMAT calls in the same expression. | // N.B. HP's aCC compiler merges FORMAT calls in the same expression. |
| |
iout << FORMAT(goTotalEnergy); | iout << FORMAT(goTotalEnergy); |
//iout << FORMAT("not implemented"); | //iout << FORMAT("not implemented"); |
} // End of port -- JLai | } // End of port -- JLai |
| |
if (simParameters->alchOn && simParameters->alchThermIntOn) { | if (simParameters->alchOn && simParameters->alchThermIntOn) { |
iout << " "; | iout << "\n"; |
| iout << TITITLE(step); |
iout << FORMAT(bondedEnergy_ti_1); | iout << FORMAT(bondedEnergy_ti_1); |
iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + | iout << FORMAT(electEnergy_ti_1 + electEnergySlow_ti_1 + |
electEnergyPME_ti_1); | electEnergyPME_ti_1); |
iout << FORMAT(ljEnergy_ti_1); | iout << FORMAT(ljEnergy_ti_1); |
iout << FORMAT(bondedEnergy_ti_2); | iout << FORMAT(bondedEnergy_ti_2); |
| iout << " "; |
iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 + | iout << FORMAT(electEnergy_ti_2 + electEnergySlow_ti_2 + |
electEnergyPME_ti_2); | electEnergyPME_ti_2); |
iout << " "; | |
iout << FORMAT(ljEnergy_ti_2); | iout << FORMAT(ljEnergy_ti_2); |
if (simParameters->alchLambdaFreq > 0) { | if (simParameters->alchLambdaFreq > 0) { |
iout << FORMAT(simParameters->getCurrentLambda(step)); | iout << FORMAT(simParameters->getCurrentLambda(step)); |
| |
} | } |
| |
//Modifications for alchemical fep | //Modifications for alchemical fep |
static char *FEPTITLE(int X) | |
{ | |
static char tmp_string[21]; | |
sprintf(tmp_string, "FepEnergy: %6d ",X); | |
return tmp_string; | |
} | |
| |
static char *TITITLE(int X) | |
{ | |
static char tmp_string[21]; | |
sprintf(tmp_string, "TI: %7d",X); | |
return tmp_string; | |
} | |
| |
void Controller::outputFepEnergy(int step) { | void Controller::outputFepEnergy(int step) { |
if (simParams->alchFepOn) { | if (simParams->alchOn && simParams->alchFepOn) { |
const int stepInRun = step - simParams->firstTimestep; | const int stepInRun = step - simParams->firstTimestep; |
const int alchEquilSteps = simParams->alchEquilSteps; | const int alchEquilSteps = simParams->alchEquilSteps; |
const BigReal alchLambda = simParams->alchLambda; | const BigReal alchLambda = simParams->alchLambda; |
| |
net_dE += dE; | net_dE += dE; |
} | } |
| |
| if (simParams->alchOutFreq) { |
if (stepInRun == 0) { | if (stepInRun == 0) { |
if (!fepFile.is_open()) { | if (!fepFile.is_open()) { |
NAMD_backup_file(simParams->alchOutFile); | NAMD_backup_file(simParams->alchOutFile); |
| |
<< "LAMBDA " << simParams->alchLambda << " COMPLETED\n" | << "LAMBDA " << simParams->alchLambda << " COMPLETED\n" |
<< "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl; | << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl; |
} | } |
if ((simParams->N) && | if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) { |
(simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0))) { | |
writeFepEnergyData(step, fepFile); | writeFepEnergyData(step, fepFile); |
fepFile.flush(); | fepFile.flush(); |
} | } |
| |
} | } |
} | } |
} | } |
| } |
| |
void Controller::outputTiEnergy(int step) { | void Controller::outputTiEnergy(int step) { |
if (simParams->alchThermIntOn) { | if (simParams->alchOn && simParams->alchThermIntOn) { |
const int stepInRun = step - simParams->firstTimestep; | const int stepInRun = step - simParams->firstTimestep; |
const int alchEquilSteps = simParams->alchEquilSteps; | const int alchEquilSteps = simParams->alchEquilSteps; |
const int alchLambdaFreq = simParams->alchLambdaFreq; | const int alchLambdaFreq = simParams->alchLambdaFreq; |
| |
net_dEdl_lj_1 = 0; | net_dEdl_lj_1 = 0; |
net_dEdl_lj_2 = 0; | net_dEdl_lj_2 = 0; |
} | } |
| TiNo++; |
| net_dEdl_bond_1 += bondedEnergy_ti_1; |
| net_dEdl_bond_2 += bondedEnergy_ti_2; |
| net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 |
| + electEnergyPME_ti_1); |
| net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 |
| + electEnergyPME_ti_2); |
| net_dEdl_lj_1 += ljEnergy_ti_1; |
| net_dEdl_lj_2 += ljEnergy_ti_2; |
| |
| // Don't accumulate block averages (you'll get a divide by zero!) or even |
| // open the TI output if alchOutFreq is zero. |
| if (simParams->alchOutFreq) { |
if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) { | if (stepInRun == 0 || (! ((step - 1) % simParams->alchOutFreq))) { |
// output of instantaneous dU/dl now replaced with running average | // output of instantaneous dU/dl now replaced with running average |
// over last alchOutFreq steps (except for step 0) | // over last alchOutFreq steps (except for step 0) |
| |
recent_dEdl_lj_2 = 0; | recent_dEdl_lj_2 = 0; |
recent_alchWork = 0; | recent_alchWork = 0; |
} | } |
TiNo++; | |
recent_TiNo++; | recent_TiNo++; |
net_dEdl_bond_1 += bondedEnergy_ti_1; | |
net_dEdl_bond_2 += bondedEnergy_ti_2; | |
// FB - PME is no longer scaled by global lambda, but by the respective | |
// lambda as dictated by elecLambdaStart. All electrostatics now go together. | |
net_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 | |
+ electEnergyPME_ti_1); | |
net_dEdl_elec_2 += (electEnergy_ti_2 + electEnergySlow_ti_2 | |
+ electEnergyPME_ti_2); | |
net_dEdl_lj_1 += ljEnergy_ti_1; | |
net_dEdl_lj_2 += ljEnergy_ti_2; | |
| |
recent_dEdl_bond_1 += bondedEnergy_ti_1; | recent_dEdl_bond_1 += bondedEnergy_ti_1; |
recent_dEdl_bond_2 += bondedEnergy_ti_2; | recent_dEdl_bond_2 += bondedEnergy_ti_2; |
recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 | recent_dEdl_elec_1 += (electEnergy_ti_1 + electEnergySlow_ti_1 |
| |
if (!tiFile.is_open()) { | if (!tiFile.is_open()) { |
NAMD_backup_file(simParams->alchOutFile); | NAMD_backup_file(simParams->alchOutFile); |
tiFile.open(simParams->alchOutFile); | tiFile.open(simParams->alchOutFile); |
/* BKR - This has been rather drastically updated to better match stdout. | /* BKR - This has been rather drastically updated to better match |
This was necessary for several reasons: | stdout. This was necessary for several reasons: |
1) PME global scaling is obsolete (now removed) | 1) PME global scaling is obsolete (now removed) |
2) scaling of bonded terms was added | 2) scaling of bonded terms was added |
3) alchemical work is now accumulated when switching is active | 3) alchemical work is now accumulated when switching is active |
| |
<< simParams->alchLambda << " --> " << simParams->alchLambda2 | << simParams->alchLambda << " --> " << simParams->alchLambda2 |
<< "\n#LAMBDA SCHEDULE: " | << "\n#LAMBDA SCHEDULE: " |
<< "dL: " << simParams->getLambdaDelta() | << "dL: " << simParams->getLambdaDelta() |
<< " Freq: " << alchLambdaFreq | << " Freq: " << alchLambdaFreq; |
<< "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K" | |
<< std::endl; | |
} | } |
else { | else { |
const BigReal alchLambda = simParams->alchLambda; | const BigReal alchLambda = simParams->alchLambda; |
| |
const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda); | const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda); |
tiFile << "#NEW TI WINDOW: " | tiFile << "#NEW TI WINDOW: " |
<< "LAMBDA " << alchLambda | << "LAMBDA " << alchLambda |
<< "\n#PARTITION 1 BOND LAMBDA " << bond_lambda_1 | << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1 |
<< "\n#PARTITION 1 VDW LAMBDA " << vdw_lambda_1 | << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1 |
<< "\n#PARTITION 1 ELEC LAMBDA " << elec_lambda_1 | << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2 |
<< "\n#PARTITION 2 BOND LAMBDA " << bond_lambda_2 | << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2; |
<< "\n#PARTITION 2 VDW LAMBDA " << vdw_lambda_2 | |
<< "\n#PARTITION 2 ELEC LAMBDA " << elec_lambda_2 | |
<< "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K" | |
<< std::endl; | |
} | } |
| tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K" |
| << std::endl; |
} | } |
| |
if (stepInRun == alchEquilSteps) { | if (stepInRun == alchEquilSteps) { |
| |
<< "LAMBDA " << simParams->alchLambda << " COMPLETED\n" | << "LAMBDA " << simParams->alchLambda << " COMPLETED\n" |
<< "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl; | << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl; |
} | } |
if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) { | if ((step%simParams->alchOutFreq) == 0) { |
writeTiEnergyData(step, tiFile); | writeTiEnergyData(step, tiFile); |
tiFile.flush(); | tiFile.flush(); |
} | } |
} | } |
} | } |
| } |
| |
/* Work is accumulated whenever alchLambda changes. In the current scheme we | /* |
| Work is accumulated whenever alchLambda changes. In the current scheme we |
always increment lambda _first_, then integrate in time. Therefore the work | always increment lambda _first_, then integrate in time. Therefore the work |
is wrt the "old" lambda before the increment. | is wrt the "old" lambda before the increment. |
*/ | */ |
| |
} | } |
} | } |
} | } |
//fepe | |
void Controller::writeTiEnergyData(int step, ofstream_namd &file) { | void Controller::writeTiEnergyData(int step, ofstream_namd &file) { |
tiFile << TITITLE(step); | tiFile << TITITLE(step); |
tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo); | tiFile << FORMAT(recent_dEdl_bond_1 / recent_TiNo); |
| |
} | } |
tiFile << std::endl; | tiFile << std::endl; |
} | } |
| //fepe |
| |
void Controller::outputExtendedSystem(int step) | void Controller::outputExtendedSystem(int step) |
{ | { |