NAMD
FreeEnergyLambdMgr.C
Go to the documentation of this file.
1 
7 // written by David Hurwitz, March to May 1998.
8 
9 #include <memory.h>
10 //#include <iomanip.h>
11 #include <stdio.h>
12 #include "InfoStream.h"
13 #include "FreeEnergyEnums.h"
14 #include "FreeEnergyAssert.h"
15 #include "FreeEnergyLambda.h"
16 #include "FreeEnergyLambdMgr.h"
17 #include "FreeEnergyGroup.h"
18 
20 //------------------------------------------------------------------------
21 // make room for some LambdaControl objects.
22 //------------------------------------------------------------------------
23  m_ActiveIndex = 0;
24  m_NumObjects = 0;
25  m_pPmfBlocks = new ALambdaControl[kLambdaNumToStart];
26  m_MaxNum = kLambdaNumToStart;
27 }
28 
29 
31 //------------------------------------------------------------------------
32 // return borrowed memory to the free store.
33 //------------------------------------------------------------------------
34  delete []m_pPmfBlocks;
35 }
36 
37 
39 //------------------------------------------------------------------------
40 // get number times dU_dLambda accumulated from active lambda object.
41 //------------------------------------------------------------------------
42  ASSERT((*this)[m_ActiveIndex].IsActive());
43  return((*this)[m_ActiveIndex].GetNum_dU_dLambda());
44 }
45 
46 
48 //------------------------------------------------------------------------
49 // integrate <dU/dLambda> for MCTI
50 //------------------------------------------------------------------------
51  ASSERT((*this)[m_ActiveIndex].IsActive());
52  (*this)[m_ActiveIndex].Integrate_MCTI();
53 }
54 
55 
56 void ALambdaManager::Accumulate(double dU_dLambda) {
57 //------------------------------------------------------------------------
58 // add dU_dLambda to the active accumulator (in one of the lambda objects)
59 //------------------------------------------------------------------------
60  ASSERT((*this)[m_ActiveIndex].IsActive());
61  (*this)[m_ActiveIndex].Accumulate(dU_dLambda);
62 }
63 
64 
66 //------------------------------------------------------------------------
67 // get the accumulation of dU_dLambda from active lambda object.
68 //------------------------------------------------------------------------
69  ASSERT((*this)[m_ActiveIndex].IsActive());
70  return((*this)[m_ActiveIndex].GetAccumulation());
71 }
72 
73 
75 //------------------------------------------------------------------------
76 // get accumulation of <dU_dLambda> * dLambda from active lambda object.
77 //------------------------------------------------------------------------
78  ASSERT((*this)[m_ActiveIndex].IsActive());
79  return((*this)[m_ActiveIndex].GetIntegration());
80 }
81 
82 
84 //------------------------------------------------------------------------
85 // zero accumulation of dU_dLambda in the active lambda object.
86 //------------------------------------------------------------------------
87  ASSERT((*this)[m_ActiveIndex].IsActive());
88  (*this)[m_ActiveIndex].ZeroAccumulator();
89 }
90 
91 
93 //------------------------------------------------------------------------
94 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
95 // decide if it's time to print restraint information
96 //------------------------------------------------------------------------
97  ASSERT((*this)[m_ActiveIndex].IsActive());
98  return((*this)[m_ActiveIndex].IsFirstStep());
99 }
100 
101 
103 //------------------------------------------------------------------------
104 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
105 // decide if it's time to print restraint information
106 //------------------------------------------------------------------------
107  ASSERT((*this)[m_ActiveIndex].IsActive());
108  return((*this)[m_ActiveIndex].IsTimeToPrint());
109 }
110 
111 
113 //------------------------------------------------------------------------
114 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
115 // decide if it's time to print du/dLambda information
116 //------------------------------------------------------------------------
117  ASSERT((*this)[m_ActiveIndex].IsActive());
118  return((*this)[m_ActiveIndex].IsTimeToPrint_dU_dLambda());
119 }
120 
121 
123 //------------------------------------------------------------------------
124 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
125 // decide if it's time to start accumulating dU/dLambda from zero
126 //------------------------------------------------------------------------
127  ASSERT((*this)[m_ActiveIndex].IsActive());
128  return((*this)[m_ActiveIndex].IsTimeToClearAccumulator());
129 }
130 
131 
133 //------------------------------------------------------------------------
134 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
135 // decide if this is the last time step of an MCTI step
136 //------------------------------------------------------------------------
137  ASSERT((*this)[m_ActiveIndex].IsActive());
138  return((*this)[m_ActiveIndex].IsEndOf_MCTI_Step());
139 }
140 
141 
143 //------------------------------------------------------------------------
144 // ASSUMING that the m_ActiveIndex'th LambdaControl is active,
145 // decide if this is the last time step of an MCTI block
146 //------------------------------------------------------------------------
147  ASSERT((*this)[m_ActiveIndex].IsActive());
148  return((*this)[m_ActiveIndex].IsEndOf_MCTI());
149 }
150 
151 
153 //------------------------------------------------------------------------
154 // print header for a new lambda control object
155 //------------------------------------------------------------------------
156  ASSERT((*this)[m_ActiveIndex].IsActive());
157  (*this)[m_ActiveIndex].PrintLambdaHeader(dT);
158 }
159 
160 
162 //------------------------------------------------------------------------
163 // print information about current time step
164 //------------------------------------------------------------------------
165  ASSERT((*this)[m_ActiveIndex].IsActive());
166  (*this)[m_ActiveIndex].PrintHeader(dT);
167 }
168 
169 
171 //------------------------------------------------------------------------
172 // return the number of steps taken in the active LambdaControl block
173 //------------------------------------------------------------------------
174  ASSERT((*this)[m_ActiveIndex].IsActive());
175  return((*this)[m_ActiveIndex].GetNumStepsSoFar());
176 }
177 
178 
180 //------------------------------------------------------------------------
181 // return the total number of steps dU/dLambda has been accumulated
182 //------------------------------------------------------------------------
183  ASSERT((*this)[m_ActiveIndex].IsActive());
184  return((*this)[m_ActiveIndex].GetNumAccumStepsSoFar());
185 }
186 
187 
189 //------------------------------------------------------------------------
190 // some stuff to make the output look nice
191 //------------------------------------------------------------------------
192 #if !defined(_VERBOSE_PMF)
193  iout << " ";
194 #endif
195 }
196 
197 
198 void ALambdaManager::Print_dU_dLambda_Summary(double Sum_dU_dLambdas) {
199 //------------------------------------------------------------------------
200 // print sum of dU/dLambda's for current time-step and
201 // the accumulation of the above.
202 //------------------------------------------------------------------------
203  char Str[100];
204 
205 #if defined(_VERBOSE_PMF)
206  iout << "FreeEnergy: ";
207  iout << "For all forcing restraints, dU/dLambda = ";
208  iout << Sum_dU_dLambdas << std::endl << endi;
209  iout << "FreeEnergy: ";
210  iout << "For all forcing restraints, Free Energy = ";
211  iout << GetAccumulation();
212  iout << " for " << GetNum_dU_dLambda() << " steps" << std::endl << endi;
213 #else
214  sprintf(Str, "%10.2e", GetAccumulation());
215  iout << Str << " ";
216  sprintf(Str, "%6d", GetNum_dU_dLambda());
217  iout << Str << " ";
218 #endif
219 }
220 
221 
223 //------------------------------------------------------------------------
224 // print the integral of: <dU/dLambda> * dLambda
225 //------------------------------------------------------------------------
226  iout << "FreeEnergy: ";
227  iout << "For MCTI, Free Energy Integral = ";
228  iout << GetIntegration();
229  iout << " for " << GetNumAccumStepsSoFar() << " steps" << std::endl << endi;
230 }
231 
232 
233 Bool_t ALambdaManager::GetLambdas(double& LambdaKf, double& LambdaRef) {
234 //------------------------------------------------------------------------
235 // get LambdaKf and LambdaRef from the active LambdaControl
236 //
237 // return(kTrue) if an active LambdaControl is found
238 // return(kFalse) if all LambdaControls have expired
239 //------------------------------------------------------------------------
240  // don't continue if all LamdaControl's have expired
241  if (m_ActiveIndex == m_NumObjects) {
242  return(kFalse);
243  }
244 
245  // if the m_ActiveIndex'th LambdaControl is no longer active
246  if ( !(*this)[m_ActiveIndex].IsActive()) {
247  // move on to the next one
248  m_ActiveIndex++;
249  // if there is no next object, return KFalse
250  if (m_ActiveIndex == m_NumObjects) {
251  return(kFalse);
252  }
253  // otherwise, make sure the next one's active
254  else {
255  ASSERT( (*this)[m_ActiveIndex].IsActive() );
256  }
257  }
258  // return LambdaKf and LambdaRef from the active LambdaControl
259  LambdaKf = (*this)[m_ActiveIndex].GetLambdaKf();
260  LambdaRef = (*this)[m_ActiveIndex].GetLambdaRef();
261  return(kTrue);
262 }
263 
264 
266 //------------------------------------------------------------------------
267 // leave memory allocation alone.
268 //------------------------------------------------------------------------
269  m_NumObjects = 0;
270 }
271 
272 
274 //------------------------------------------------------------------------
275 // add an object to the list. if there's not enough room, make room.
276 // return an index to the added oject.
277 //------------------------------------------------------------------------
278  ALambdaControl* pPmfBlocks;
279 
280  // if there's no room for another object
281  if (m_NumObjects == m_MaxNum) {
282  // create an array with more space
283  m_MaxNum *= kLambdaMultiplier;
284  pPmfBlocks = new ALambdaControl[m_MaxNum];
285  // copy from the full array to the new one
286  for (int i=0; i<m_NumObjects; i++) {
287  pPmfBlocks[i] = m_pPmfBlocks[i];
288  }
289  // return the space used for the full array
290  delete []m_pPmfBlocks;
291  // point to the bigger array
292  m_pPmfBlocks = pPmfBlocks;
293  }
294  // add the object to the array
295  m_pPmfBlocks[m_NumObjects] = PmfBlock;
296  m_NumObjects++;
297  return(m_NumObjects-1);
298 }
299 
300 
302 //------------------------------------------------------------------------
303 // return an object from this group of objects.
304 //------------------------------------------------------------------------
305  ASSERT((Index>=0) && (Index<m_NumObjects));
306  return(m_pPmfBlocks[Index]);
307 }
308 
309 
311 //------------------------------------------------------------------------
312 // calculate and return the total number of steps needed for all
313 // pmf and mcti blocks
314 //------------------------------------------------------------------------
315  int Total, i;
316 
317  Total = 0;
318  for (i=0; i<m_NumObjects; i++) {
319  Total += m_pPmfBlocks[i].GetNumSteps();
320  }
321  return(Total);
322 }
323 
Bool_t IsTimeToPrint_dU_dLambda()
int Add(ALambdaControl &PmfBlock)
Bool_t IsTimeToClearAccumulator()
void PrintLambdaHeader(double dT)
void Accumulate(double dU_dLambda)
const int kLambdaMultiplier
#define iout
Definition: InfoStream.h:87
void PrintHeader(double dT)
int Index
Definition: structures.h:26
void Print_dU_dLambda_Summary(double Sum_dU_dLambdas)
Bool_t
#define ASSERT(E)
Bool_t GetLambdas(double &LambdaKf, double &LambdaRef)
Bool_t IsEndOf_MCTI_Step()
infostream & endi(infostream &s)
Definition: InfoStream.C:38
const int kLambdaNumToStart
double GetAccumulation()
ALambdaControl & operator[](int index)