Difference for src/Controller.C from version 1.1320 to 1.1321

version 1.1320version 1.1321
Line 7
Line 7
 /***************************************************************************** /*****************************************************************************
  * $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"
Line 1261
Line 1261
 //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;
Line 1276
Line 1277
 }  } 
 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 
Line 2250
Line 2251
  
     // 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);
Line 2273
Line 2274
  
     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);
Line 2785
Line 2800
         } 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);
         }          } 
       }       }
Line 2894
Line 2909
    //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");
Line 2910
Line 2926
         }         }
  
  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.
Line 2967
Line 2984
       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));
Line 3065
Line 3084
 } }
  
 //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;
Line 3103
Line 3108
     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);
Line 3137
Line 3143
             << "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();
   }   }
Line 3151
Line 3156
   }   }
  }  }
 } }
  }
  
 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;
Line 3167
Line 3173
     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)
Line 3179
Line 3198
     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 
Line 3206
Line 3213
     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
Line 3240
Line 3247
              << 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;    
Line 3254
Line 3259
       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) {
Line 3270
Line 3273
            << "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.
 */ */
Line 3384
Line 3389
     }     }
   }   }
 } }
 //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);
Line 3407
Line 3412
   }   }
   tiFile << std::endl;   tiFile << std::endl;
 } }
  //fepe
  
 void Controller::outputExtendedSystem(int step) void Controller::outputExtendedSystem(int step)
 { {

Legend:
Removed in v.1.1320 
changed lines
 Added in v.1.1321


Made by using version 1.53 of cvs2html