NAMD
GlobalMasterFreeEnergy.C
Go to the documentation of this file.
1 
7 /*
8  Forwards atoms to master node for force evaluation.
9 */
10 
11 #include <string.h>
12 #include "InfoStream.h"
13 #include "NamdTypes.h"
14 #include "FreeEnergyEnums.h"
15 #include "FreeEnergyAssert.h"
16 #include "FreeEnergyGroup.h"
17 #include "Vector.h"
18 #include "FreeEnergyVector.h"
19 #include "FreeEnergyRestrain.h"
20 #include "FreeEnergyRMgr.h"
21 #include "FreeEnergyLambda.h"
22 #include "FreeEnergyLambdMgr.h"
23 #include "GlobalMaster.h"
24 #include "GlobalMasterFreeEnergy.h"
25 #include "FreeEnergyParse.h"
26 
27 #include "Node.h"
28 #include "Molecule.h"
29 #include "ReductionMgr.h"
30 #include "SimParameters.h"
31 
32 #include <stdio.h>
33 
34 //#define DEBUGM
35 #define MIN_DEBUG_LEVEL 1
36 #include "Debug.h"
37 
38 
39 void GlobalMasterFreeEnergy::update() {
40 //-----------------------------------------------------------------
41 // get lambdas from LambdaManager, inform RestraintManager.
42 // calculate gradients for each center-of-mass of each restraint,
43 // and apply the forces to the atoms involved
44 //-----------------------------------------------------------------
45  double LambdaKf, LambdaRef;
46  double Sum_dU_dLambdas;
47 
48  if (m_LambdaManager.GetLambdas(LambdaKf, LambdaRef)) {
49 
50  // stuff that's done every time step
51  m_RestraintManager.SetLambdas(LambdaKf, LambdaRef);
52  m_RestraintManager.UpdateCOMs(*this);
53  m_RestraintManager.AddForces(*this);
54  if (m_LambdaManager.IsTimeToClearAccumulator()) {
55  m_LambdaManager.ZeroAccumulator();
56  }
57  Sum_dU_dLambdas = m_RestraintManager.Sum_dU_dLambdas();
58  m_LambdaManager.Accumulate(Sum_dU_dLambdas);
59 
60  // for integrating all the MCTI averages
61  if (m_LambdaManager.IsEndOf_MCTI_Step()) {
62  m_LambdaManager.Integrate_MCTI();
63  }
64 
65  // stuff that's done when it's time to print
66  if (m_LambdaManager.IsFirstStep()) {
67  m_LambdaManager.PrintLambdaHeader(simParams->dt);
68  }
69  if (m_LambdaManager.IsTimeToPrint()) {
70  m_LambdaManager.PrintHeader(simParams->dt);
71  if (m_LambdaManager.IsTimeToPrint_dU_dLambda()) {
72  m_RestraintManager.Print_dU_dLambda_Info();
73  if (m_RestraintManager.ThereIsAForcingRestraint()) {
74  m_LambdaManager.Print_dU_dLambda_Summary(Sum_dU_dLambdas);
75  }
76  }
77  else {
78  m_LambdaManager.PrintSomeSpaces();
79  }
80  m_RestraintManager.PrintEnergyInfo();
81  m_RestraintManager.PrintRestraintInfo();
82  if (m_LambdaManager.IsEndOf_MCTI()) {
83  m_LambdaManager.Print_MCTI_Integration();
84  }
85  }
86  }
87 }
88 
89 
90 void GlobalMasterFreeEnergy::user_initialize() {
91 //-----------------------------------------------------------------
92 // read all the input from config
93 //-----------------------------------------------------------------
94 
95  iout << iINFO << " FREE ENERGY PERTURBATION CONFIG\n";
96  iout << iINFO << "***********************************\n";
97  int config_len = strlen(config);
98  if ( config_len < 10000 ) {
99  iout << config;
100  } else {
101  char *new_config = new char[10000 + 10];
102  strncpy(new_config,config,10000);
103  new_config[10000] = 0;
104  strcat(new_config,"\n...\n");
105  iout << new_config;
106  delete [] new_config;
107  }
108  iout << iINFO << "***********************************\n" << endi;
109 
110  ReadInput(config, m_RestraintManager, m_LambdaManager, *this, simParams->dt);
111 
112  // exit if there aren't enough steps to complete all pmf & mcti blocks
113  int Total = m_LambdaManager.GetTotalNumSteps();
114  if (Total > simParams->N) {
115  iout << "FreeEnergy: Not enough steps to complete pfm & mcti blocks" << std::endl;
116  iout << "FreeEnergy: Num Steps Needed = " << Total << std::endl;
117  iout << "FreeEnergy: Num Steps Requested = " << simParams->N << std::endl << endi;
118  NAMD_die("FreeEnergy: Fatal Run-Time Error");
119  }
120 }
121 
122 
123 void GlobalMasterFreeEnergy::user_calculate() {
124 //-----------------------------------------------------------------
125 // this is what's executed every time-step
126 //-----------------------------------------------------------------
127  m_LambdaManager.IncCurrStep();
128  update();
129 }
130 
131 
133  getAtomID(const char *segid, int resid, const char *aname)
134 {
135  return molecule->get_atom_from_name(segid,resid,aname);
136 }
137 
139  getNumAtoms(const char* segid, int resid) // 0 on error
140 {
141  return molecule->get_residue_size(segid,resid);
142 }
143 
145  getAtomID(const char *segid, int resid, int index)
146 {
147  return molecule->get_atom_from_index_in_residue(segid,resid,index);
148 }
149 
151 {
152  if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1.; // failure
153  return molecule->atommass(atomid);
154 }
155 
156 
158 {
159  if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1; // failure
160  modifyRequestedAtoms().add(atomid);
161  return 0; // success
162 }
163 
165 {
169  for ( ; a_i != a_e; ++a_i, ++p_i ) {
170  if ( *a_i == atomid ) {
171  position = *p_i;
172  return 0; // success
173  }
174  }
175  return -1; // failure
176 }
177 
179 {
180  DebugM(2,"Forcing "<<atomid<<" with "<<force<<"\n");
181  if ( atomid < 0 || atomid >= molecule->numAtoms ) return -1; // failure
182  modifyForcedAtoms().add(atomid);
183  modifyAppliedForces().add(force);
184  return 0; // success
185 }
186 
187 
189 : GlobalMaster() {
190  DebugM(3,"Constructing GlobalMasterFreeEnergy\n");
191  molecule = Node::Object()->molecule;
192  simParams = Node::Object()->simParameters;
193 
194  // now set up the free energy stuff
195  initialize();
196 }
197 
199  DebugM(3,"Destructing GlobalMasterFreeEnergy\n");
200  delete config;
201 }
202 
203 
204 void GlobalMasterFreeEnergy::initialize() {
205  DebugM(4,"Initializing master\n");
206 
207  // Get our script
208  StringList *script = Node::Object()->configList->find("freeEnergyConfig");
209 
210  config = new char[1];
211  config[0] = '\0';
212 
213  for ( ; script; script = script->next) {
214  if ( strstr(script->data,"\n") ) {
215  size_t add_len = strlen(script->data);
216  size_t config_len = 0;
217  config_len = strlen(config);
218  char *new_config = new char[config_len + add_len + 2];
219  strcpy(new_config,config);
220  strcat(new_config,script->data);
221  strcat(new_config,"\n"); // just to be safe
222  delete [] config;
223  config = new_config;
224  } else {
225  FILE *infile = fopen(script->data,"r");
226  if ( ! infile ) {
227  char errmsg[256];
228  sprintf(errmsg,"Error trying to read file %s!\n",script->data);
229  NAMD_die(errmsg);
230  }
231  fseek(infile,0,SEEK_END);
232  size_t add_len = ftell(infile);
233  size_t config_len = 0;
234  config_len = strlen(config);
235  char *new_config = new char[config_len + add_len + 3];
236  strcpy(new_config,config);
237  delete [] config;
238  config = new_config;
239  new_config += config_len;
240  rewind(infile);
241  fread(new_config,sizeof(char),add_len,infile);
242  new_config += add_len;
243  new_config[0] = '\n';
244  new_config[1] = '\0';
245  fclose(infile);
246  }
247  }
248 
249  // iout << iDEBUG << "Free energy perturbation - initialize()\n" << endi;
250  user_initialize();
251 }
252 
253 
254 void GlobalMasterFreeEnergy::calculate() {
255  DebugM(4,"Calculating forces on master\n");
256 
257  /* zero out the forces */
260 
261  /* XXX is this line needed at all? */
263  modifyGroupForces().setall(Vector(0,0,0));
264 
265 // iout << iDEBUG << "Free energy perturbation - calculate()\n" << endi;
266  user_calculate();
267 
268  // Send results to clients
269  DebugM(3,"Sending results (" << forcedAtoms().size() << " forces) on master\n");
270  if ( changedAtoms() || changedGroups() ) {
271  DebugM(4,"Sending new configuration (" <<
272  requestedAtoms().size() << " atoms) on master\n");
273  }
274 }
static Node * Object()
Definition: Node.h:86
Bool_t IsTimeToPrint_dU_dLambda()
int get_atom_from_name(const char *segid, int resid, const char *aname) const
Definition: Molecule.C:121
int get_residue_size(const char *segid, int resid) const
Definition: Molecule.C:144
int getPosition(int atomid, Position &position)
ForceList & modifyAppliedForces()
Definition: GlobalMaster.C:162
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
AtomIDList & modifyRequestedAtoms()
Definition: GlobalMaster.C:127
Bool_t IsTimeToClearAccumulator()
void SetLambdas(double LambdaKf, double LambdaRef)
int getAtomID(const char *segid, int resid, const char *aname)
const Elem * const_iterator
Definition: ResizeArray.h:38
void PrintLambdaHeader(double dT)
void Accumulate(double dU_dLambda)
AtomIDList::const_iterator getAtomIdBegin()
Definition: GlobalMaster.C:190
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
double Sum_dU_dLambdas()
bool changedAtoms()
Definition: GlobalMaster.C:107
#define DebugM(x, y)
Definition: Debug.h:59
int get_atom_from_index_in_residue(const char *segid, int resid, int index) const
Definition: Molecule.C:158
#define iout
Definition: InfoStream.h:87
void PrintHeader(double dT)
PositionList::const_iterator getAtomPositionBegin()
Definition: GlobalMaster.C:198
int addForce(int atomid, Force force)
BigRealList::const_iterator getGroupMassEnd()
Definition: GlobalMaster.C:239
void Print_dU_dLambda_Info()
void Print_dU_dLambda_Summary(double Sum_dU_dLambdas)
void ReadInput(char *Str, ARestraintManager &RMgr, ALambdaManager &LMgr, GlobalMasterFreeEnergy &CFE, double dT)
void setall(const Elem &elem)
Definition: ResizeArray.h:90
bool changedGroups()
Definition: GlobalMaster.C:115
AtomIDList & modifyForcedAtoms()
Definition: GlobalMaster.C:157
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
Bool_t GetLambdas(double &LambdaKf, double &LambdaRef)
BigRealList::const_iterator getGroupMassBegin()
Definition: GlobalMaster.C:234
ForceList & modifyGroupForces()
Definition: GlobalMaster.C:167
ConfigList * configList
Definition: Node.h:179
void UpdateCOMs(GlobalMasterFreeEnergy &CFE)
int add(const Elem &elem)
Definition: ResizeArray.h:97
StringList * next
Definition: ConfigList.h:49
void resize(int i)
Definition: ResizeArray.h:84
AtomIDList::const_iterator getAtomIdEnd()
Definition: GlobalMaster.C:194
char * data
Definition: ConfigList.h:48
StringList * find(const char *name) const
Definition: ConfigList.C:341
Bool_t IsEndOf_MCTI_Step()
void AddForces(GlobalMasterFreeEnergy &CFE)
Real atommass(int anum) const
Definition: Molecule.h:1038
Bool_t ThereIsAForcingRestraint()
infostream & endi(infostream &s)
Definition: InfoStream.C:38
const AtomIDList & forcedAtoms()
Definition: GlobalMaster.C:133
Molecule * molecule
Definition: Node.h:176
int getNumAtoms(const char *segid, int resid)
const AtomIDList & requestedAtoms()
Definition: GlobalMaster.C:123