GlobalMasterFreeEnergy.C

Go to the documentation of this file.
00001 
00007 /*
00008    Forwards atoms to master node for force evaluation.
00009 */
00010 
00011 #include <string.h>
00012 #include "InfoStream.h"
00013 #include "NamdTypes.h"
00014 #include "FreeEnergyEnums.h"
00015 #include "FreeEnergyAssert.h"
00016 #include "FreeEnergyGroup.h"
00017 #include "Vector.h"
00018 #include "FreeEnergyVector.h"
00019 #include "FreeEnergyRestrain.h"
00020 #include "FreeEnergyRMgr.h"
00021 #include "FreeEnergyLambda.h"
00022 #include "FreeEnergyLambdMgr.h"
00023 #include "GlobalMaster.h"
00024 #include "GlobalMasterFreeEnergy.h"
00025 #include "FreeEnergyParse.h"
00026 
00027 #include "Node.h"
00028 #include "Molecule.h"
00029 #include "ReductionMgr.h"
00030 #include "SimParameters.h"
00031 
00032 #include <stdio.h>
00033 
00034 //#define DEBUGM
00035 #define MIN_DEBUG_LEVEL 1
00036 #include "Debug.h"
00037 
00038 
00039 void GlobalMasterFreeEnergy::update() {
00040 //-----------------------------------------------------------------
00041 // get lambdas from LambdaManager, inform RestraintManager.
00042 // calculate gradients for each center-of-mass of each restraint,
00043 // and apply the forces to the atoms involved
00044 //-----------------------------------------------------------------
00045   double  LambdaKf, LambdaRef;
00046   double  Sum_dU_dLambdas;
00047 
00048   if (m_LambdaManager.GetLambdas(LambdaKf, LambdaRef)) {
00049 
00050     // stuff that's done every time step
00051     m_RestraintManager.SetLambdas(LambdaKf, LambdaRef);
00052     m_RestraintManager.UpdateCOMs(*this);
00053     m_RestraintManager.AddForces(*this);
00054     if (m_LambdaManager.IsTimeToClearAccumulator()) {
00055       m_LambdaManager.ZeroAccumulator();
00056     }
00057     Sum_dU_dLambdas = m_RestraintManager.Sum_dU_dLambdas();
00058     m_LambdaManager.Accumulate(Sum_dU_dLambdas);
00059 
00060     // for integrating all the MCTI averages
00061     if (m_LambdaManager.IsEndOf_MCTI_Step()) {
00062       m_LambdaManager.Integrate_MCTI();
00063     }
00064 
00065     // stuff that's done when it's time to print
00066     if (m_LambdaManager.IsFirstStep()) {
00067       m_LambdaManager.PrintLambdaHeader(simParams->dt);
00068     }
00069     if (m_LambdaManager.IsTimeToPrint()) {
00070       m_LambdaManager.PrintHeader(simParams->dt);
00071       if (m_LambdaManager.IsTimeToPrint_dU_dLambda()) {
00072         m_RestraintManager.Print_dU_dLambda_Info();
00073         if (m_RestraintManager.ThereIsAForcingRestraint()) {
00074           m_LambdaManager.Print_dU_dLambda_Summary(Sum_dU_dLambdas);
00075         }
00076       }
00077       else {
00078         m_LambdaManager.PrintSomeSpaces();
00079       }
00080       m_RestraintManager.PrintEnergyInfo();
00081       m_RestraintManager.PrintRestraintInfo();
00082       if (m_LambdaManager.IsEndOf_MCTI()) {
00083         m_LambdaManager.Print_MCTI_Integration();
00084       }
00085     }
00086   }
00087 }
00088 
00089 
00090 void GlobalMasterFreeEnergy::user_initialize() {
00091 //-----------------------------------------------------------------
00092 // read all the input from config
00093 //-----------------------------------------------------------------
00094 
00095   iout << iINFO << "  FREE ENERGY PERTURBATION CONFIG\n"; 
00096   iout << iINFO << "***********************************\n"; 
00097   int config_len = strlen(config);
00098   if ( config_len < 10000 ) {
00099     iout << config;
00100   } else {
00101     char *new_config = new char[10000 + 10];
00102     strncpy(new_config,config,10000);
00103     new_config[10000] = 0;
00104     strcat(new_config,"\n...\n");
00105     iout << new_config;
00106     delete [] new_config;
00107   }
00108   iout << iINFO << "***********************************\n" << endi; 
00109 
00110   ReadInput(config, m_RestraintManager, m_LambdaManager, *this, simParams->dt);
00111 
00112   // exit if there aren't enough steps to complete all pmf & mcti blocks
00113   int Total = m_LambdaManager.GetTotalNumSteps();
00114   if (Total > simParams->N) {
00115     iout << "FreeEnergy: Not enough steps to complete pfm & mcti blocks" << std::endl;
00116     iout << "FreeEnergy:   Num Steps Needed =    " << Total << std::endl;
00117     iout << "FreeEnergy:   Num Steps Requested = " << simParams->N << std::endl << endi;
00118     NAMD_die("FreeEnergy: Fatal Run-Time Error");
00119   }
00120 }
00121 
00122 
00123 void GlobalMasterFreeEnergy::user_calculate() {
00124 //-----------------------------------------------------------------
00125 // this is what's executed every time-step
00126 //-----------------------------------------------------------------
00127   m_LambdaManager.IncCurrStep();
00128   update();
00129 }
00130 
00131 
00132 int GlobalMasterFreeEnergy::
00133 	getAtomID(const char *segid, int resid, const char *aname)
00134 {
00135   return molecule->get_atom_from_name(segid,resid,aname);
00136 }
00137 
00138 int GlobalMasterFreeEnergy::
00139 	getNumAtoms(const char* segid, int resid) // 0 on error
00140 {
00141   return molecule->get_residue_size(segid,resid);
00142 }
00143 
00144 int GlobalMasterFreeEnergy::
00145 	getAtomID(const char *segid, int resid, int index)
00146 {
00147   return molecule->get_atom_from_index_in_residue(segid,resid,index);
00148 }
00149 
00150 double GlobalMasterFreeEnergy::getMass(int atomid)
00151 {
00152   if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1.;  // failure
00153   return molecule->atommass(atomid);
00154 }
00155 
00156 
00157 int GlobalMasterFreeEnergy::requestAtom(int atomid)
00158 {
00159   if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1;  // failure
00160   modifyRequestedAtoms().add(atomid);
00161   return 0;  // success
00162 }
00163 
00164 int GlobalMasterFreeEnergy::getPosition(int atomid, Position &position)
00165 {
00166   AtomIDList::const_iterator a_i = getAtomIdBegin();
00167   AtomIDList::const_iterator a_e = getAtomIdEnd();
00168   PositionList::const_iterator p_i = getAtomPositionBegin();
00169   for ( ; a_i != a_e; ++a_i, ++p_i ) {
00170     if ( *a_i == atomid ) {
00171       position = *p_i;
00172       return 0;  // success
00173     }
00174   }
00175   return -1;  // failure
00176 }
00177 
00178 int GlobalMasterFreeEnergy::addForce(int atomid, Force force)
00179 {
00180   DebugM(2,"Forcing "<<atomid<<" with "<<force<<"\n");
00181   if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1;  // failure
00182   modifyForcedAtoms().add(atomid);
00183   modifyAppliedForces().add(force);
00184   return 0;  // success
00185 }
00186 
00187 
00188 GlobalMasterFreeEnergy::GlobalMasterFreeEnergy() 
00189 : GlobalMaster() {
00190   DebugM(3,"Constructing GlobalMasterFreeEnergy\n");
00191   molecule = Node::Object()->molecule;
00192   simParams = Node::Object()->simParameters;
00193 
00194   // now set up the free energy stuff
00195   initialize();
00196 }
00197 
00198 GlobalMasterFreeEnergy::~GlobalMasterFreeEnergy() {
00199   DebugM(3,"Destructing GlobalMasterFreeEnergy\n");
00200   delete config;
00201 }
00202 
00203 
00204 void GlobalMasterFreeEnergy::initialize() {
00205   DebugM(4,"Initializing master\n");
00206 
00207   // Get our script
00208   StringList *script = Node::Object()->configList->find("freeEnergyConfig");
00209 
00210   config = new char[1];
00211   config[0] = '\0';
00212 
00213   for ( ; script; script = script->next) {
00214     if ( strstr(script->data,"\n") ) {
00215       size_t add_len = strlen(script->data);
00216       size_t config_len = 0;
00217       config_len = strlen(config);
00218       char *new_config = new char[config_len + add_len + 2];
00219       strcpy(new_config,config);
00220       strcat(new_config,script->data);
00221       strcat(new_config,"\n");  // just to be safe
00222       delete [] config;
00223       config = new_config;
00224     } else {
00225       FILE *infile = fopen(script->data,"r");
00226       if ( ! infile ) {
00227         char errmsg[256];
00228         sprintf(errmsg,"Error trying to read file %s!\n",script->data);
00229         NAMD_die(errmsg);
00230       }
00231       fseek(infile,0,SEEK_END);
00232       size_t add_len = ftell(infile);
00233       size_t config_len = 0;
00234       config_len = strlen(config);
00235       char *new_config = new char[config_len + add_len + 3];
00236       strcpy(new_config,config);
00237       delete [] config;
00238       config = new_config;
00239       new_config += config_len;
00240       rewind(infile);
00241       fread(new_config,sizeof(char),add_len,infile);
00242       new_config += add_len;
00243       new_config[0] = '\n';
00244       new_config[1] = '\0';
00245       fclose(infile);
00246     }
00247   }
00248 
00249   // iout << iDEBUG << "Free energy perturbation - initialize()\n" << endi; 
00250   user_initialize();
00251 }
00252 
00253 
00254 void GlobalMasterFreeEnergy::calculate() {
00255   DebugM(4,"Calculating forces on master\n");
00256   
00257   /* zero out the forces */
00258   modifyForcedAtoms().resize(0);
00259   modifyAppliedForces().resize(0);
00260 
00261   /* XXX is this line needed at all? */
00262   modifyGroupForces().resize(getGroupMassEnd() - getGroupMassBegin());
00263   modifyGroupForces().setall(Vector(0,0,0));
00264 
00265 //  iout << iDEBUG << "Free energy perturbation - calculate()\n" << endi; 
00266   user_calculate();
00267 
00268   // Send results to clients
00269   DebugM(3,"Sending results (" << forcedAtoms().size() << " forces) on master\n");
00270   if ( changedAtoms() || changedGroups() ) {
00271     DebugM(4,"Sending new configuration (" <<
00272                         requestedAtoms().size() << " atoms) on master\n");
00273   }
00274 }

Generated on Tue Sep 26 01:17:13 2017 for NAMD by  doxygen 1.4.7