NAMD
FreeEnergyRMgr.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 <string.h>
11 // #include <iomanip.h>
12 #include "InfoStream.h"
13 #include "FreeEnergyEnums.h"
14 #include "FreeEnergyAssert.h"
15 #include "Vector.h"
16 #include "FreeEnergyVector.h"
17 #include "FreeEnergyGroup.h"
18 #include "FreeEnergyRestrain.h"
19 #include "FreeEnergyRMgr.h"
20 #include "FreeEnergyLambda.h"
21 #include "FreeEnergyLambdMgr.h"
22 
23 #include "NamdTypes.h"
24 #include "GlobalMaster.h"
25 #include "GlobalMasterFreeEnergy.h"
26 
28 //------------------------------------------------------------------------
29 // allocate space for restraint POINTERS
30 //------------------------------------------------------------------------
31  m_ppRestraints = new pRestr[kNumToStart];
32  m_NumRestraints = 0;
33  m_MaxNum = kNumToStart;
34 }
35 
36 
38 //------------------------------------------------------------------------
39 // when the RestraintManager goes out of scope,
40 // free the space that was allocated for each restraint
41 // (see: ARestraint* GetRestraint(char*, int&)
42 // then, free the space that was allocated for the pointers
43 //------------------------------------------------------------------------
44  for (int i=0; i<m_NumRestraints; i++) {
45  delete m_ppRestraints[i];
46  }
47  delete []m_ppRestraints;
48 }
49 
50 
52 //------------------------------------------------------------------------
53 // get a pointer
54 //------------------------------------------------------------------------
55  ASSERT( (Index>=0) && (Index<m_NumRestraints) );
56  return(m_ppRestraints[Index]);
57 }
58 
59 
61 //------------------------------------------------------------------------
62 // add a pointer to the list. if there's not enough room, make room.
63 //------------------------------------------------------------------------
64  ARestraint** ppRestraints;
65 
66  // if there's no room for a new pointer
67  if (m_NumRestraints == m_MaxNum) {
68  // create an array with more space
69  m_MaxNum *= kMultiplier;
70  ppRestraints = new pRestr[m_MaxNum];
71  // fast copy from the full array to the new one (memcpy(dest, src, bytes))
72  memcpy(ppRestraints, m_ppRestraints, sizeof(ARestraint*)*m_NumRestraints);
73  // return the space used for the full array
74  delete []m_ppRestraints;
75  // point to the bigger array
76  m_ppRestraints = ppRestraints;
77  }
78 
79  // add the int to the int array
80  m_ppRestraints[m_NumRestraints] = pRestraint;
81  m_NumRestraints++;
82 }
83 
84 
86 //------------------------------------------------------------------------
87 // update the centers-of-mass in each restraint in the list
88 // note: m_ppRestraints points to a list of restraint pointers
89 // m_ppRestraints[i] references one of the pointers
90 // m_ppRestraints[i]-> accesses a restraint function
91 //------------------------------------------------------------------------
92  for (int i=0; i<m_NumRestraints; i++) {
93  m_ppRestraints[i]->UpdateCOMs(CFE);
94  }
95 }
96 
97 
99 //---------------------------------------------------------------------------
100 // for each restraint, apply restraining force to each COM (center-of-mass).
101 //---------------------------------------------------------------------------
102  int i, j, NumCOMs;
103  AVector Force;
104 
105  // for each restraint
106  for (i=0; i<m_NumRestraints; i++) {
107  // for each center-of-mass
108  NumCOMs = m_ppRestraints[i]->GetNumGroups();
109  for (j=0; j<NumCOMs; j++) {
110  Force = m_ppRestraints[i]->GetGradient(j);
111  // apply restraining force in opposite direction from gradient
112  Force *= -1.0;
113  m_ppRestraints[i]->DistributeForce(j, Force, CFE);
114  }
115  }
116 }
117 
118 
120 //---------------------------------------------------------------------------
121 // sum up dU/dLambda from each forcing restraint
122 //---------------------------------------------------------------------------
123  double Sum=0;
124 
125  for (int i=0; i<m_NumRestraints; i++) {
126  Sum += m_ppRestraints[i]->Get_dU_dLambda();
127  }
128  return(Sum);
129 }
130 
131 
133 //---------------------------------------------------------------------------
134 // return kTrue if there's at least one forcing restraint
135 //---------------------------------------------------------------------------
136  for (int i=0; i<m_NumRestraints; i++) {
137  if (m_ppRestraints[i]->IsForcing()) {
138  return(kTrue);
139  }
140  }
141  return(kFalse);
142 }
143 
144 
146 //---------------------------------------------------------------------------
147 // for a restraint, print restraint type and Energy.
148 //---------------------------------------------------------------------------
149 #if defined(_VERBOSE_PMF)
150  for (int i=0; i<m_NumRestraints; i++) {
151  PrintPreInfo(i);
152  iout << "Energy = ";
153  iout << m_ppRestraints[i]->GetEnergy() << std::endl << endi;
154  }
155 #endif
156 }
157 
158 
160 //---------------------------------------------------------------------------
161 // for a restraint, print its position, distance, angle, or dihedral angle.
162 //---------------------------------------------------------------------------
163  for (int i=0; i<m_NumRestraints; i++) {
164 #if defined(_VERBOSE_PMF)
165  PrintPreInfo(i);
166 #endif
167  m_ppRestraints[i]->PrintInfo();
168  }
169 #if !defined(_VERBOSE_PMF)
170  iout << std::endl << endi;
171 #endif
172 }
173 
174 
176 //---------------------------------------------------------------------------
177 // if restraint is a forcing restraint, print dU/dLambda.
178 //---------------------------------------------------------------------------
179 #if defined(_VERBOSE_PMF)
180  for (int i=0; i<m_NumRestraints; i++) {
181  if (m_ppRestraints[i]->IsForcing()) {
182  PrintPreInfo(i);
183  iout << "dU/dLambda = ";
184  iout << m_ppRestraints[i]->Get_dU_dLambda() << std::endl << endi;
185  }
186  }
187 #endif
188 }
189 
190 
192 //---------------------------------------------------------------------------
193 // print "Restraint xxx: Type of Restraint: "
194 //---------------------------------------------------------------------------
195  char Str[100];
196  char NumStr[20];
197 
198  sprintf(NumStr, "%3d", Index+1);
199  iout << "FreeEnergy: " << "Restraint " << NumStr << ": ";
200  m_ppRestraints[Index]->GetStr(Str);
201  iout << Str << ": ";
202 }
203 
Vector Force
Definition: NamdTypes.h:26
void PrintPreInfo(int Index)
double Sum_dU_dLambdas()
#define iout
Definition: InfoStream.h:87
const int kMultiplier
virtual void PrintInfo()=0
void Print_dU_dLambda_Info()
void Add(ARestraint *pRestraint)
int Index
Definition: structures.h:26
virtual double GetEnergy()=0
Bool_t
#define ASSERT(E)
virtual double Get_dU_dLambda()
void DistributeForce(int WhichGroup, AVector Force, GlobalMasterFreeEnergy &CFE)
void UpdateCOMs(GlobalMasterFreeEnergy &CFE)
const int kNumToStart
void UpdateCOMs(GlobalMasterFreeEnergy &CFE)
void AddForces(GlobalMasterFreeEnergy &CFE)
Bool_t ThereIsAForcingRestraint()
infostream & endi(infostream &s)
Definition: InfoStream.C:38
ARestraint * operator[](int Index)
virtual void GetStr(char *Str)=0
virtual AVector GetGradient(int WhichGroup)=0