version 1.1318 | version 1.1319 |
---|
| |
/***************************************************************************** | /***************************************************************************** |
* $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/05/24 19:41:49 $ | * $Date: 2016/06/08 22:06:21 $ |
* $Revision: 1.1318 $ | * $Revision: 1.1319 $ |
*****************************************************************************/ | *****************************************************************************/ |
| |
#include "InfoStream.h" | #include "InfoStream.h" |
| |
const BigReal accelMDTalpha = simParams->accelMDTalpha; | const BigReal accelMDTalpha = simParams->accelMDTalpha; |
const int accelMDOutFreq = simParams->accelMDOutFreq; | const int accelMDOutFreq = simParams->accelMDOutFreq; |
| |
//BigReal bondEnergy; | BigReal bondEnergy; |
//BigReal angleEnergy; | BigReal angleEnergy; |
//BigReal dihedralEnergy; | BigReal dihedralEnergy; |
//BigReal improperEnergy; | BigReal improperEnergy; |
// BKR - bonded terms are now declared globally for easy FEP differencing | |
BigReal crosstermEnergy; | BigReal crosstermEnergy; |
BigReal boundaryEnergy; | BigReal boundaryEnergy; |
BigReal miscEnergy; | BigReal miscEnergy; |
| |
// Drude model ANISO energy is added into BOND energy | // Drude model ANISO energy is added into BOND energy |
// and THOLE energy is added into ELECT energy | // and THOLE energy is added into ELECT energy |
| |
//BigReal bondEnergy; | BigReal bondEnergy; |
//BigReal angleEnergy; | BigReal angleEnergy; |
//BigReal dihedralEnergy; | BigReal dihedralEnergy; |
//BigReal improperEnergy; | BigReal improperEnergy; |
// BKR - bonded terms are now declared globally for easy FEP differencing | |
BigReal crosstermEnergy; | BigReal crosstermEnergy; |
BigReal boundaryEnergy; | BigReal boundaryEnergy; |
BigReal miscEnergy; | BigReal miscEnergy; |
| |
goTotalEnergy = goNativeEnergy + goNonnativeEnergy; | goTotalEnergy = goNativeEnergy + goNonnativeEnergy; |
| |
//fepb | //fepb |
bondedEnergy_f = reduction->item(REDUCTION_BONDED_ENERGY_F); | bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F); |
electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F); | electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F); |
ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F); | ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F); |
ljEnergy_f_left = reduction->item(REDUCTION_LJ_ENERGY_F_LEFT); | ljEnergy_f_left = reduction->item(REDUCTION_LJ_ENERGY_F_LEFT); |
| |
//fepb BKR - Compute alchemical work if using dynamic lambda. This is here so | //fepb BKR - Compute alchemical work if using dynamic lambda. This is here so |
// that the cumulative work can be given during a callback. | // that the cumulative work can be given during a callback. |
if (simParameters->alchLambdaFreq > 0) { | if (simParameters->alchLambdaFreq > 0) { |
| if (step <= |
| simParameters->firstTimestep + simParameters->alchEquilSteps) { |
| cumAlchWork = 0.0; |
| } |
alchWork = computeAlchWork(step); | alchWork = computeAlchWork(step); |
cumAlchWork += alchWork; | cumAlchWork += alchWork; |
} | } |
| |
CALLBACKLIST("CUMALCHWORK", cumAlchWork); | CALLBACKLIST("CUMALCHWORK", cumAlchWork); |
} | } |
} else if (simParameters->alchFepOn) { | } else if (simParameters->alchFepOn) { |
CALLBACKLIST("BOND2", bondedEnergy_f); | CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy + |
| improperEnergy + bondedEnergyDiff_f); |
CALLBACKLIST("ELEC2", electEnergy_f); | CALLBACKLIST("ELEC2", electEnergy_f); |
CALLBACKLIST("VDW2", ljEnergy_f); | CALLBACKLIST("VDW2", ljEnergy_f); |
} | } |
| |
exp_dE_ByRT = 0.0; | exp_dE_ByRT = 0.0; |
net_dE = 0.0; | net_dE = 0.0; |
} | } |
BigReal dE = bondedEnergy_f + electEnergy_f + electEnergySlow_f + ljEnergy_f | dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f - |
- (bondEnergy + angleEnergy + dihedralEnergy + improperEnergy | electEnergy - electEnergySlow - ljEnergy; |
+ electEnergy + electEnergySlow + ljEnergy); | |
BigReal RT = BOLTZMANN * simParams->alchTemp; | BigReal RT = BOLTZMANN * simParams->alchTemp; |
| iout << step << " bondedEnergyDiff: " << bondedEnergyDiff_f << "\n" << endi; |
| |
if (alchEnsembleAvg){ | if (alchEnsembleAvg){ |
FepNo++; | FepNo++; |
| |
<< "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) && (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0))) { | if ((simParams->N) && |
| (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0))) { |
writeFepEnergyData(step, fepFile); | writeFepEnergyData(step, fepFile); |
fepFile.flush(); | fepFile.flush(); |
} | } |
if (alchEnsembleAvg && (step == simParams->N)) { | if (alchEnsembleAvg && (step == simParams->N)) { |
fepSum = fepSum + dG; | fepSum = fepSum + dG; |
fepFile << "#Free energy change for lambda window [ " << alchLambda | fepFile << "#Free energy change for lambda window [ " << alchLambda |
<< " " << alchLambda2 << " ] is " << dG << " ; net change until now is " << fepSum << std::endl; | << " " << alchLambda2 << " ] is " << dG |
| << " ; net change until now is " << fepSum << std::endl; |
fepFile.flush(); | fepFile.flush(); |
} | } |
} | } |
| |
if (simParams->alchThermIntOn) { | if (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 stepInSwitch = stepInRun - alchEquilSteps; | |
const int alchLambdaFreq = simParams->alchLambdaFreq; | const int alchLambdaFreq = simParams->alchLambdaFreq; |
| |
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 (alchLambdaFreq > 0 && stepInSwitch >= 0 && step != simParams->N ) { | |
// Work is accumulated whenever alchLambda changes. In the current | |
// scheme we always increment lambda _first_, then integrate in time. | |
// Therefore the work is wrt the "old" lambda before the increment. | |
alchWork = computeAlchWork(step); | |
cumAlchWork += alchWork; | |
} | |
if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) { | if (simParams->alchOutFreq && ((step%simParams->alchOutFreq)==0)) { |
writeTiEnergyData(step, tiFile); | writeTiEnergyData(step, tiFile); |
tiFile.flush(); | tiFile.flush(); |
| |
void Controller::writeFepEnergyData(int step, ofstream_namd &file) { | void Controller::writeFepEnergyData(int step, ofstream_namd &file) { |
BigReal eeng = electEnergy+electEnergySlow; | BigReal eeng = electEnergy+electEnergySlow; |
BigReal eeng_f = electEnergy_f + electEnergySlow_f; | BigReal eeng_f = electEnergy_f + electEnergySlow_f; |
BigReal dE = eeng_f + ljEnergy_f - eeng - ljEnergy; | |
BigReal dE_Left = eeng_f + ljEnergy_f_left - eeng - ljEnergy; | BigReal dE_Left = eeng_f + ljEnergy_f_left - eeng - ljEnergy; |
BigReal RT = BOLTZMANN * simParams->alchTemp; | BigReal RT = BOLTZMANN * simParams->alchTemp; |
| |
const bool alchEnsembleAvg = simParams->alchEnsembleAvg; | const bool alchEnsembleAvg = simParams->alchEnsembleAvg; |
const int stepInRun = step - simParams->firstTimestep; | const int stepInRun = step - simParams->firstTimestep; |
const bool FepWhamOn = simParams->alchFepWhamOn; | const bool FepWhamOn = simParams->alchFepWhamOn; |
| |
const BigReal alchRepLambda = simParams->alchRepLambda; | const BigReal alchRepLambda = simParams->alchRepLambda; |
const BigReal alchDispLambda = simParams->alchDispLambda; | const BigReal alchDispLambda = simParams->alchDispLambda; |
const BigReal alchElecLambda = simParams->alchElecLambda; | const BigReal alchElecLambda = simParams->alchElecLambda; |
| |
if(stepInRun){ | if(stepInRun){ |
if(!FepWhamOn){ | if(!FepWhamOn){ |
fepFile << FEPTITLE(step); | fepFile << FEPTITLE(step); |