NAMD
FreeEnergyLambda.C
Go to the documentation of this file.
1 
7 // written by David Hurwitz, March to May 1998.
8 
9 #include <string.h>
10 #include "InfoStream.h"
11 #include "FreeEnergyEnums.h"
12 #include "FreeEnergyAssert.h"
13 #include "FreeEnergyLambda.h"
14 #include "Vector.h"
15 #include "FreeEnergyVector.h"
16 
17 // initialize static member variables
18 int ALambdaControl::m_CurrStep = 0;
19 
21 //------------------------------------------------------------------------
22 // initialize member variables
23 //------------------------------------------------------------------------
24  // initialize these for the first block.
25  m_StartStep = 0;
26  m_LambdaKf = 1.0; // applied to Kf
27  m_LambdaRef = 0.0; // applied to reference (pos, dist, angle, dihe)
28  m_Sum_dU_dLambda = 0.0; // for free energy measurement
29  m_Num_dU_dLambda = 0; // "
30  m_MCTI_Integration = 0.0; // "
31 
32  // set the rest to illegal values.
33  m_Task = kUnknownTask;
34  m_NumSteps = -1;
35  m_NumEquilSteps = -1;
36  m_NumAccumSteps = -1;
37  m_NumPrintSteps = -1;
38  m_NumRepeats = -1;
39  m_StopStep = -1;
40 }
41 
42 
44 //------------------------------------------------------------------------
45 // initialize this object using the settings for the prior object.
46 //------------------------------------------------------------------------
47  m_Task = PriorBlock.m_Task;
48  m_NumSteps = PriorBlock.m_NumSteps;
49  m_NumEquilSteps = PriorBlock.m_NumEquilSteps;
50  m_NumAccumSteps = PriorBlock.m_NumAccumSteps;
51  m_NumPrintSteps = PriorBlock.m_NumPrintSteps;
52  m_NumRepeats = PriorBlock.m_NumRepeats;
53  switch (PriorBlock.m_Task) {
54 
55  case kUp: m_LambdaKf=1.0; m_LambdaRef=1.0; break;
56  case kStepUp: m_LambdaKf=1.0; m_LambdaRef=1.0; break;
57  case kDown: m_LambdaKf=1.0; m_LambdaRef=0.0; break;
58  case kStepDown: m_LambdaKf=1.0; m_LambdaRef=0.0; break;
59  case kStop: m_LambdaKf=1.0; m_LambdaRef=PriorBlock.m_LambdaRef; break;
60 
61  case kGrow: m_LambdaKf=1.0; m_LambdaRef=PriorBlock.m_LambdaRef; break;
62  case kStepGrow: m_LambdaKf=1.0; m_LambdaRef=PriorBlock.m_LambdaRef; break;
63  case kFade: m_LambdaKf=0.0; m_LambdaRef=PriorBlock.m_LambdaRef; break;
64  case kStepFade: m_LambdaKf=0.0; m_LambdaRef=PriorBlock.m_LambdaRef; break;
65  case kNoGrow: m_LambdaKf = PriorBlock.m_LambdaKf;
66  m_LambdaRef = PriorBlock.m_LambdaRef; break;
67  default: ASSERT(kFalse); break; //should never get here
68 
69  }
70  m_StartStep = PriorBlock.GetLastStep();
71 }
72 
73 
75 //------------------------------------------------------------------------
76 // copy everything from PmfBlock to this block.
77 //------------------------------------------------------------------------
78  m_NumSteps = PmfBlock.m_NumSteps;
79  m_NumEquilSteps = PmfBlock.m_NumEquilSteps;
80  m_NumAccumSteps = PmfBlock.m_NumAccumSteps;
81  m_NumRepeats = PmfBlock.m_NumRepeats;
82  m_NumPrintSteps = PmfBlock.m_NumPrintSteps;
83  m_StartStep = PmfBlock.m_StartStep;
84  m_StopStep = PmfBlock.m_StopStep;
85  m_LambdaKf = PmfBlock.m_LambdaKf;
86  m_LambdaRef = PmfBlock.m_LambdaRef;
87  m_Task = PmfBlock.m_Task;
88  return(*this);
89 }
90 
91 
92 void ALambdaControl::Accumulate(double dU_dLambda) {
93 //------------------------------------------------------------------------
94 // if lambda is constant, sum dU/dLambda
95 // if lambda is changing, sum dU/dLambda * dLambda
96 //------------------------------------------------------------------------
97 
98  m_Num_dU_dLambda++;
99  switch (m_Task) {
100  // lambda is constant
101  case kStop: case kNoGrow: case kStepUp:
102  case kStepDown: case kStepGrow: case kStepFade:
103  m_Sum_dU_dLambda += dU_dLambda;
104  break;
105  // lambda is increasing
106  case kUp:
107  case kGrow:
108  m_Sum_dU_dLambda += dU_dLambda / (double)m_NumSteps;
109  break;
110  // lambda is decreasing
111  case kDown:
112  case kFade:
113  m_Sum_dU_dLambda -= dU_dLambda / (double)m_NumSteps;
114  break;
115  // should never get here
116  default:
117  ASSERT(kFalse);
118  break;
119  }
120 }
121 
122 
124 //------------------------------------------------------------------------
125 // integrate MCTI: <dU/dLambda> * dLambda
126 //------------------------------------------------------------------------
127  ASSERT(m_Task==kStepUp || m_Task==kStepDown ||
128  m_Task==kStepGrow || m_Task==kStepFade);
129 
130  switch (m_Task) {
131  // lambda is increasing
132  case kStepUp:
133  case kStepGrow:
134  m_MCTI_Integration += GetAccumulation() / (double)m_NumRepeats;
135  break;
136  // lambda is decreasing
137  case kStepDown:
138  case kStepFade:
139  m_MCTI_Integration -= GetAccumulation() / (double)m_NumRepeats;
140  break;
141  // should never get here
142  default:
143  ASSERT(kFalse);
144  break;
145  }
146 }
147 
148 
150 //------------------------------------------------------------------------
151 // get MCTI integral. integral(dU/dLambda> * dLambda
152 //------------------------------------------------------------------------
153  ASSERT(m_Task==kStepUp || m_Task==kStepDown ||
154  m_Task==kStepGrow || m_Task==kStepFade);
155 
156  return(m_MCTI_Integration);
157 }
158 
159 
161 //------------------------------------------------------------------------
162 // if lambda is constant, return (sum/N)
163 // if lambda is changing, return (sum)
164 //------------------------------------------------------------------------
165  switch (m_Task) {
166  // lambda is constant
167  case kStop: case kNoGrow:
168  case kStepUp: case kStepDown: case kStepGrow: case kStepFade:
169  ASSERT(m_Num_dU_dLambda != 0);
170  return(m_Sum_dU_dLambda / (double)m_Num_dU_dLambda);
171  // lambda is changing
172  case kUp: case kDown: case kGrow: case kFade:
173  return(m_Sum_dU_dLambda);
174  // should never get here
175  default:
176  ASSERT(kFalse);
177  return(0);
178  }
179 }
180 
181 
183 //------------------------------------------------------------------------
184 // return the total number of steps dU/dLambda has been accumulated
185 // this is only called during MCTI
186 //------------------------------------------------------------------------
187  ASSERT(m_Task==kStepUp || m_Task==kStepDown ||
188  m_Task==kStepGrow || m_Task==kStepFade);
189 
190  int Count, Maybe;
191 
192  Count = (m_CurrStep-m_StartStep) / (m_NumAccumSteps+m_NumEquilSteps);
193  Count *= m_NumAccumSteps;
194  Maybe = (m_CurrStep-m_StartStep) % (m_NumAccumSteps+m_NumEquilSteps);
195  Maybe -= m_NumEquilSteps;
196  if (Maybe > 0) {
197  Count += Maybe;
198  }
199  return(Count);
200 }
201 
202 
204 //------------------------------------------------------------------------
205 // get the number of steps needed for this pmf or mcti block
206 //------------------------------------------------------------------------
207  // make sure m_StopStep is calculated
208  GetLastStep();
209 
210  return( (m_StopStep - m_StartStep) + 1 );
211 }
212 
213 
214 Bool_t ALambdaControl::IsLastStep() {
215 //------------------------------------------------------------------------
216 // return true if we're on the last step of this pmf or mcti block
217 //------------------------------------------------------------------------
218  if (m_CurrStep == GetLastStep()) {
219  return(kTrue);
220  }
221  else {
222  return(kFalse);
223  }
224 }
225 
226 
227 int ALambdaControl::GetLastStep() {
228 //------------------------------------------------------------------------
229 // get the last step of this task
230 //------------------------------------------------------------------------
231  // if it's already calculated, just return it
232  if (m_StopStep > 0) {
233  return(m_StopStep);
234  }
235  // otherwise calculate it
236  switch (m_Task) {
237  case kStepUp:
238  case kStepDown:
239  case kStepGrow:
240  case kStepFade:
241  m_StopStep = m_StartStep +
242  (m_NumAccumSteps+m_NumEquilSteps) * m_NumRepeats;
243  break;
244  default:
245  m_StopStep = m_StartStep + m_NumSteps;
246  break;
247  }
248  // and return it
249  return(m_StopStep);
250 }
251 
252 
254 //------------------------------------------------------------------------
255 // get a string that describes this task
256 //------------------------------------------------------------------------
257  switch (m_Task) {
258  case kUp: strcpy(Str, " Up"); break;
259  case kDown: strcpy(Str, " Down"); break;
260  case kStop: strcpy(Str, " Stop"); break;
261  case kGrow: strcpy(Str, " Grow"); break;
262  case kFade: strcpy(Str, " Fade"); break;
263  case kNoGrow: strcpy(Str, " NoGrow"); break;
264  case kStepUp: strcpy(Str, " StepUp"); break;
265  case kStepDown: strcpy(Str, "StepDown"); break;
266  case kStepGrow: strcpy(Str, "StepGrow"); break;
267  case kStepFade: strcpy(Str, "StepFade"); break;
268  default: strcpy(Str, "Bug Alert!!!"); break;
269  }
270 }
271 
272 
273 void ALambdaControl::GetTaskStr(char* Str) {
274 //------------------------------------------------------------------------
275 // get a string that describes this task
276 //------------------------------------------------------------------------
277  switch (m_Task) {
278  case kUp: strcpy(Str, "Up"); break;
279  case kDown: strcpy(Str, "Down"); break;
280  case kStop: strcpy(Str, "Stop"); break;
281  case kGrow: strcpy(Str, "Grow"); break;
282  case kFade: strcpy(Str, "Fade"); break;
283  case kNoGrow: strcpy(Str, "NoGrow"); break;
284  case kStepUp: strcpy(Str, "StepUp"); break;
285  case kStepDown: strcpy(Str, "StepDown"); break;
286  case kStepGrow: strcpy(Str, "StepGrow"); break;
287  case kStepFade: strcpy(Str, "StepFade"); break;
288  default: strcpy(Str, "Bug Alert!!!"); break;
289  }
290 }
291 
292 
294 //------------------------------------------------------------------------
295 // ASSUMING that this object is currently active, decide if it's time
296 // to start accumulating dU/dLambda from zero.
297 // (clear the accumulator on the 1st step)
298 //------------------------------------------------------------------------
299  Bool_t RetVal=kFalse;
300 
301  switch (m_Task) {
302  // for pmf blocks, clear accumulator on first step
303  case kUp: case kDown: case kStop:
304  case kGrow: case kFade: case kNoGrow:
305  if (m_CurrStep == m_StartStep) {
306  RetVal = kTrue;
307  }
308  break;
309  // for mcti blocks, clear accumulator after each equilibration
310  case kStepUp: case kStepDown: case kStepGrow: case kStepFade:
311  if ( (m_CurrStep-(m_StartStep+m_NumEquilSteps)) %
312  (m_NumEquilSteps+m_NumAccumSteps) == 1) {
313  RetVal = kTrue;
314  }
315  break;
316  // should never get here
317  default:
318  ASSERT(kFalse);
319  break;
320  }
321  return(RetVal);
322 }
323 
324 
326 //------------------------------------------------------------------------
327 // ASSUMING that this object is currently active, decide if it's the
328 // first step of the control object.
329 //------------------------------------------------------------------------
330  Bool_t RetVal=kFalse;
331 
332  if (m_CurrStep == (m_StartStep+1)) {
333  RetVal = kTrue;
334  }
335  return(RetVal);
336 }
337 
338 
340 //------------------------------------------------------------------------
341 // ASSUMING that this object is currently active, decide if it's time
342 // to print out restraint information.
343 //------------------------------------------------------------------------
344  Bool_t RetVal=kFalse;
345 
346  // if printing is required
347  if (m_NumPrintSteps > 0) {
348  // if number-of-steps from StartStep is an even multiple of NumPrintSteps
349  // or if it's the last step of this pmf or mcti block,
350  // then it's time to print
351  if ( IsLastStep() || (((m_CurrStep-m_StartStep)%m_NumPrintSteps)==0) ) {
352  RetVal = kTrue;
353  }
354  }
355  return(RetVal);
356 }
357 
358 
360 //------------------------------------------------------------------------
361 // ASSUMING that this object is currently active, decide if this is
362 // the last time step of an mcti block.
363 //------------------------------------------------------------------------
364  Bool_t RetVal=kFalse;
365 
366  // if an MCTI block is currently active
367  if ( (m_Task==kStepUp) || (m_Task==kStepDown) ||
368  (m_Task==kStepGrow) || (m_Task==kStepFade) ) {
369  // if this is the last step
370  if (m_CurrStep == GetLastStep()) {
371  RetVal = kTrue;
372  }
373  }
374  return(RetVal);
375 }
376 
377 
379 //------------------------------------------------------------------------
380 // ASSUMING that this object is currently active, decide if this is
381 // the last time step of an mcti step.
382 //------------------------------------------------------------------------
383  Bool_t RetVal=kFalse;
384 
385  // if an MCTI block is currently active
386  if ( (m_Task==kStepUp) || (m_Task==kStepDown) ||
387  (m_Task==kStepGrow) || (m_Task==kStepFade) ) {
388  // if this is the last step of accumulation
389  if (((m_CurrStep-m_StartStep)%(m_NumEquilSteps+m_NumAccumSteps))==0) {
390  // then this is the last time step of this mcti step
391  RetVal = kTrue;
392  }
393  }
394  return(RetVal);
395 }
396 
397 
399 //------------------------------------------------------------------------
400 // ASSUMING that this object is currently active, decide if it's time
401 // to print out dU/dLambda information.
402 //------------------------------------------------------------------------
403  Bool_t RetVal=kFalse;
404 
405  // if printing is required
406  if (m_NumPrintSteps > 0) {
407  // if number-of-steps from StartStep is an even multiple of NumPrintSteps
408  // or if it's the last step of this pmf or mcti block,
409  // then it might be time to print du/dLambda
410  if ( IsLastStep() || (((m_CurrStep-m_StartStep)%m_NumPrintSteps)==0) ) {
411  // for mcti blocks
412  if ((m_Task==kStepUp) || (m_Task==kStepDown) ||
413  (m_Task==kStepGrow) || (m_Task==kStepFade)) {
414  // it's only time to print if we're no longer equilibrating
415  if ( ((m_CurrStep-m_StartStep-1) % (m_NumEquilSteps+m_NumAccumSteps))
416  > m_NumEquilSteps) {
417  RetVal = kTrue;
418  }
419  }
420  // for other blocks (up, down, grow, fade) it's time
421  else {
422  RetVal = kTrue;
423  }
424  }
425  }
426  return(RetVal);
427 }
428 
429 
431 //----------------------------------------------------------------------------
432 // ASSUMING that this object is currently active, write out the header for
433 // a new lambda control object
434 //----------------------------------------------------------------------------
435  // double Time;
436 
437  // calculate current time in femto-seconds
438  // Time = (double)m_CurrStep * dT;
439 
440 iout << "FreeEnergy: ";
441 #if !defined(_VERBOSE_PMF)
442  iout << "nstep time(ps) ";
443  iout << " task lambdaKf lambdaRef delta-G #steps n*{value target |}" << std::endl;
444  iout << "FreeEnergy: ----- -------- ";
445  iout << "-------- -------- --------- ---------- ------ -------------------" << std::endl;
446  iout << endi;
447 #endif
448 }
449 
450 
452 //----------------------------------------------------------------------------
453 // ASSUMING that this object is currently active, write out the current time
454 //----------------------------------------------------------------------------
455  double Time;
456  char Str[100], Str2[100];
457 
458  // calculate current time in femto-seconds
459  Time = (double)m_CurrStep * dT;
460 
461 #if defined(_VERBOSE_PMF)
462  iout << "FreeEnergy: " << std::endl << endi;
463  iout << "FreeEnergy: ";
464  iout << "Time Step = " << m_CurrStep << ", ";
465  iout << "Time = ";
466  // write out time in ps
467  iout << Time/1000.0 << " ps, ";
468  iout << "Lambda_Kf = " << m_LambdaKf << ", ";
469  iout << "Lambda_Ref = " << m_LambdaRef << " ";
470  GetTaskStr(Str);
471  iout << "(" << Str << ")";
472  iout << std::endl << endi;
473  iout << "FreeEnergy: ";
474  iout << "------------------------------------------------";
475  iout << "-------------------";
476  iout << std::endl << endi;
477 #else
478  sprintf(Str, "%5d", m_CurrStep);
479  // write out time in ps
480  sprintf(Str2, "%8.3f", Time/1000.0);
481  iout << "FreeEnergy: ";
482  iout << Str << " " << Str2 << " ";
483  GetPaddedTaskStr(Str);
484  iout << Str << " ";
485  sprintf(Str, "%8.5f", m_LambdaKf);
486  iout << Str << " ";
487  sprintf(Str, "%9.5f", m_LambdaRef);
488  iout << Str << " ";
489 #endif
490 }
491 
492 
494 //------------------------------------------------------------------------
495 // determine if this object is currently active
496 //------------------------------------------------------------------------
497  if ( (m_CurrStep>=m_StartStep) && (m_CurrStep<=GetLastStep()) ) {
498  return(kTrue);
499  }
500  else {
501  return(kFalse);
502  }
503 }
504 
505 
507 //------------------------------------------------------------------------
508 // calculate LambdaKf for grow, fade, stepgrow, and stepfade.
509 // LambdaKf=1.0, for up, down, stepup, stepdown, and stop.
510 // for nogrow, LambdaKf is Lambda from the config file.
511 //------------------------------------------------------------------------
512  int N;
513  double CurrStep = m_CurrStep;
514  double StartStep = m_StartStep;
515  double NumSteps = m_NumSteps;
516  double NumEquilSteps = m_NumEquilSteps;
517  double NumAccumSteps = m_NumAccumSteps;
518  double NumRepeats = m_NumRepeats;
519 
520  if (IsActive()) {
521  switch (m_Task) {
522  case kGrow:
523  m_LambdaKf = (CurrStep-StartStep)/NumSteps;
524  break;
525  case kFade:
526  m_LambdaKf = 1.0-(CurrStep-StartStep)/NumSteps;
527  break;
528  case kStepGrow:
529  N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
530  m_LambdaKf = (N+1)/NumRepeats;
531  break;
532  case kStepFade:
533  N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
534  m_LambdaKf = 1.0 - (N+1)/NumRepeats;
535  break;
536  case kNoGrow:
537  break; // return prior setting of m_LambdaKf
538  default:
539  m_LambdaKf=1.0;
540  }
541  }
542  else {
543  m_LambdaKf=1.0;
544  }
545  return(m_LambdaKf);
546 }
547 
548 
550 //------------------------------------------------------------------------
551 // calculate LambdaRef for up, down, stepup, and stepdown.
552 // for stop, LambdaRef is Lambda from the config file.
553 // for grow, fade, stepgrow, stepfade, and nogrow,
554 // LambdaRef is LambdaT from the config file.
555 //------------------------------------------------------------------------
556  int N;
557  double CurrStep = m_CurrStep;
558  double StartStep = m_StartStep;
559  double NumSteps = m_NumSteps;
560  double NumEquilSteps = m_NumEquilSteps;
561  double NumAccumSteps = m_NumAccumSteps;
562  double NumRepeats = m_NumRepeats;
563 
564  if (IsActive()) {
565  switch (m_Task) {
566  case kUp:
567  m_LambdaRef = (CurrStep-StartStep)/NumSteps;
568  break;
569  case kDown:
570  m_LambdaRef = 1.0-(CurrStep-StartStep)/NumSteps;
571  break;
572  case kStepUp:
573  N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
574  m_LambdaRef = (N+1)/NumRepeats;
575  break;
576  case kStepDown:
577  N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
578  m_LambdaRef = 1.0 - (N+1)/NumRepeats;
579  default:
580  break; // return prior setting of m_LambdaRef
581  }
582  }
583  else {
584  m_LambdaRef=0.0;
585  }
586  return(m_LambdaRef);
587 }
588 
void GetPaddedTaskStr(char *Str)
void GetTaskStr(char *Str)
double GetIntegration()
void PrintHeader(double dT)
Bool_t IsTimeToPrint_dU_dLambda()
double GetAccumulation()
#define iout
Definition: InfoStream.h:87
Bool_t
void PrintLambdaHeader(double dT)
#define ASSERT(E)
void Accumulate(double dU_dLambda)
ALambdaControl & operator=(ALambdaControl &PmfBlock)
Bool_t IsTimeToPrint()
void Init(ALambdaControl &PriorBlock)
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Bool_t IsTimeToClearAccumulator()
Bool_t IsEndOf_MCTI_Step()