NAMD
FreeEnergyRestrain.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 <iomanip.h>
11 #include "common.h"
12 #include "InfoStream.h"
13 #include "FreeEnergyAssert.h"
14 #include "FreeEnergyEnums.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 
27 // initialize static member variables
28 // (lambda is the same for all forcing restraints)
29 double ARestraint::m_LambdaKf = 1.0;
30 double ARestraint::m_LambdaRef = 0.0;
31 
32 #if defined(_DEBUG)
33 void ARestraint::SetCOM(int Index, AVector& Pos) {
34 //-----------------------------------------------------------------
35 // only for testing!!!!!!
36 //-----------------------------------------------------------------
37  ASSERT( (Index>=0) && (Index<m_NumGroups) );
38  m_pCOMs[Index] = Pos;
39 }
40 #endif
41 
43 //-----------------------------------------------------------------
44 // constructor for base class
45 //-----------------------------------------------------------------
46  m_pGroups = NULL;
47  m_pCOMs = NULL;
48  m_NumGroups = 0;
49 }
50 
51 
53 //-----------------------------------------------------------------
54 // free space that may have been allocated for Groups and COM's
55 //-----------------------------------------------------------------
56  if (m_pGroups != NULL) {
57  ASSERT(m_pCOMs != NULL);
58  delete []m_pGroups;
59  delete []m_pCOMs;
60  }
61 }
62 
63 
64 void ARestraint::EarlyExit(const char* Str, int AtomID) {
65 //-----------------------------------------------------------------
66 // unrecoverable error
67 //-----------------------------------------------------------------
68  iout << "FreeEnergy: " << std::endl << endi;
69  iout << "FreeEnergy: " << Str;
70  iout << " for AtomID: " << AtomID;
71  iout << std::endl << endi;
72  NAMD_die("FreeEnergy: Fatal Error with Fixed or Forcing Restraints");
73 }
74 
75 
76 void ARestraint::SetGroup(AGroup& Group, int GroupIndex) {
77 //-----------------------------------------------------------------
78 // set one group of atoms
79 //-----------------------------------------------------------------
80  ASSERT( (GroupIndex>=0) && (GroupIndex<m_NumGroups) );
81  m_pGroups[GroupIndex] = Group;
82 }
83 
84 
86 //-----------------------------------------------------------------
87 // set one group of atoms
88 //-----------------------------------------------------------------
89  ASSERT(m_NumGroups >= 1);
90  m_pGroups[0] = Group1;
91 }
92 
93 
94 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2) {
95 //-----------------------------------------------------------------
96 // set two groups of atoms
97 //-----------------------------------------------------------------
98  ASSERT(m_NumGroups >= 2);
99  m_pGroups[0] = Group1;
100  m_pGroups[1] = Group2;
101 }
102 
103 
104 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2, AGroup& Group3) {
105 //-----------------------------------------------------------------
106 // set three groups of atoms
107 //-----------------------------------------------------------------
108  ASSERT(m_NumGroups >= 3);
109  m_pGroups[0] = Group1;
110  m_pGroups[1] = Group2;
111  m_pGroups[2] = Group3;
112 }
113 
114 
115 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2, AGroup& Group3, AGroup& Group4) {
116 //-----------------------------------------------------------------
117 // set four groups of atoms
118 //-----------------------------------------------------------------
119  ASSERT(m_NumGroups >= 4);
120  m_pGroups[0] = Group1;
121  m_pGroups[1] = Group2;
122  m_pGroups[2] = Group3;
123  m_pGroups[3] = Group4;
124 }
125 
126 
128 //-----------------------------------------------------------------
129 // determine the angle formed by the points A-B-C
130 //-----------------------------------------------------------------
131  double u;
132  double a = B.Dist(C);
133  double b = A.Dist(C);
134  double c = A.Dist(B);
135 
136  u = (a*a + c*c - b*b) / (2.0*a*c);
137  // protect against acos(<-1.0) and acos(>1.0)
138  if (u < -1.0) {u = -1.0;}
139  if (u > 1.0) {u = 1.0;}
140  return(acos(u));
141 }
142 
143 
145 //-----------------------------------------------------------------
146 // determine the dihedral angle formed by the points A-B-C-D
147 //-----------------------------------------------------------------
148  AVector CD(D - C);
149  AVector CB(B - C);
150  AVector BC(C - B);
151  AVector BA(A - B);
152  AVector CDxCB, BCxBA;
153  double top, bot, cos_u, sin_u, Angle;
154  AVector topVec;
155 
156  CDxCB = CD.cross(CB);
157  BCxBA = BC.cross(BA);
158 
159  top = CDxCB.dot(BCxBA);
160  bot = CDxCB.Dist() * BCxBA.Dist();
161  cos_u = top/bot;
162 
163  // protect against acos(<-1.0) and acos(>1.0)
164  if (cos_u < -1.0) {cos_u = -1.0;}
165  if (cos_u > 1.0) {cos_u = 1.0;}
166 
167  topVec = CDxCB.cross(BCxBA);
168  sin_u = (topVec/bot).dot(CB/CB.Dist());
169 
170  // protect against asin(<-1.0) and asin(>1.0)
171  if (sin_u < -1.0) {sin_u = -1.0;}
172  if (sin_u > 1.0) {sin_u = 1.0;}
173 
174  Angle = atan2(sin_u, cos_u);
175  return(Angle);
176 }
177 
178 
180  GlobalMasterFreeEnergy& CFE) {
181 //----------------------------------------------------------------------
182 // Distribute Force among the group of atoms specified by WhichGroup
183 //
184 // note: m_pGroups points to an array of Groups
185 // m_pGroups[WhichGroup] references one of the Groups
186 // m_pGroups[WhichGroup][i] returns an AtomID from the Group
187 // (operator[] is defined to return an item from the Group)
188 //----------------------------------------------------------------------
189  int i, AtomID, NumAtoms, RetVal;
190  double Mass, TotalMass=0;
191  AVector SmallForce;
192  Vector NAMD_Vector;
193 
194  ASSERT( (WhichGroup>=0) && (WhichGroup<m_NumGroups) );
195 
196  // calculate the total mass for the group
197  NumAtoms = m_pGroups[WhichGroup].GetNumInGroup();
198  for (i=0; i<NumAtoms; i++) {
199  AtomID = m_pGroups[WhichGroup][i];
200  Mass = CFE.getMass(AtomID);
201  if (Mass < 0) {EarlyExit("Negative Mass", AtomID);};
202  TotalMass += Mass;
203  }
204 
205  // distribute Force according to mass of each atom in the group
206  for (i=0; i<NumAtoms; i++) {
207  AtomID = m_pGroups[WhichGroup][i];
208  Mass = CFE.getMass(AtomID);
209  if (Mass < 0) {EarlyExit("Negative Mass", AtomID);}
210  SmallForce = Force * (Mass/TotalMass);
211  // cast SmallForce to a NAMD-type vector (addForce uses Vector class)
212  SetEqual(NAMD_Vector, SmallForce);
213  RetVal = CFE.addForce(AtomID, NAMD_Vector);
214  if (RetVal < 0) {EarlyExit("Can't add Force", AtomID);}
215  }
216 }
217 
218 
220 //-----------------------------------------------------------------
221 // calculate the center-of-mass of each group of atoms
222 //
223 // note: m_pGroups points to an array of Groups
224 // m_pGroups[i] references one of the Groups
225 // m_pGroups[i][j] returns an AtomID from the Group
226 // (operator[] is defined to return an item from the Group)
227 //-----------------------------------------------------------------
228  int i, j, AtomID, RetVal;
229  Vector NAMD_Vector;
230  AVector COM, Pos;
231  double Mass, TotalMass;
232 
233  ASSERT(m_NumGroups > 0);
234  // for each group of atoms
235  for (i=0; i<m_NumGroups; i++) {
236  TotalMass = 0;
237  COM.Set(0,0,0);
238  // for each atom in the group
239  for (j=0; j<m_pGroups[i].GetNumInGroup(); j++) {
240  AtomID = m_pGroups[i][j];
241  // get its position, weight position with atom's mass
242  RetVal = CFE.getPosition(AtomID, NAMD_Vector);
243  if (RetVal < 0) {EarlyExit("Can't get Position", AtomID);}
244  // cast NAMD_Vector to AVector (getPosition uses Vector class)
245  SetEqual(Pos, NAMD_Vector);
246  Mass = CFE.getMass(AtomID);
247  if (Mass < 0) {EarlyExit("Negative Mass", AtomID);}
248  TotalMass += Mass;
249  COM += Pos * Mass;
250  }
251  m_pCOMs[i] = COM / TotalMass;
252  }
253 }
254 
255 
257 //-----------------------------------------------------------------
258 // each APosRestraint restrains 1 group of atoms to a location
259 //-----------------------------------------------------------------
260  m_NumGroups = 1;
262  m_pCOMs = new AVector[m_NumGroups];
263 }
264 
265 
267 //--------------------------------------------------------------------
268 // print the position for this position restraint
269 //--------------------------------------------------------------------
270  double Distance;
271  char Str[20];
272 
273  Distance = GetDistance();
274  sprintf(Str, "%7.3f", Distance);
275 
276 #if defined(_VERBOSE_PMF)
277  iout << "Position = ";
278  m_pCOMs[0].Out();
279  iout << " Target = ";
280  GetPosTarget().Out();
281  iout << " Distance = ";
282  iout << Str;
283  iout << std::endl << endi;
284 #else
285  m_pCOMs[0].Out();
286  iout << " ";
287  GetPosTarget().Out();
288  iout << " ";
289  iout << Str;
290  iout << " | ";
291 #endif
292 }
293 
294 
295 double APosRestraint::GetE(AVector RefPos, double LambdaKf) {
296 //--------------------------------------------------------------------
297 // calculate and return the Energy for this position restraint.
298 //
299 // E = (Kf/2) * (|ri - rref|)**2
300 //
301 // where |ri - rref| is the distance between a) the center-of-mass
302 // of the restrained atoms and b) the reference position.
303 //
304 // Note: COM is calculated before this routine is called.
305 //--------------------------------------------------------------------
306  return( ((m_Kf*LambdaKf)/2.0) * m_pCOMs[0].DistSqr(RefPos) );
307 }
308 
309 
310 AVector APosRestraint::GetGrad(int /* WhichGroup */,
311  AVector RefPos, double LambdaKf) {
312 //-------------------------------------------------------------------------
313 // calculate and return the gradient for this position restraint.
314 //
315 // E = (Kf/2) * (|ri - rref|)**2
316 //
317 // return: grad(E)
318 //
319 // Notes: COM is calculated before this routine is called.
320 // m_pCOMs points to an array of vectors
321 // m_pCOMS[0] references the only COM for a position restraint
322 // m_pCOMS[0][0] returns the x value from the vector (x,y,z)
323 // (operator[] is defined to return an item from the vector)
324 //-------------------------------------------------------------------------
325  // WhichGroup = 0; // don't care -- there's only 1 atom restrained
326  AVector Vec(m_pCOMs[0][0] - RefPos[0],
327  m_pCOMs[0][1] - RefPos[1],
328  m_pCOMs[0][2] - RefPos[2]);
329  return(Vec*m_Kf*LambdaKf);
330 }
331 
332 
334 //-----------------------------------------------------------------------
335 // each ADistRestraint restrains the distance between 2 groups of atoms
336 //-----------------------------------------------------------------------
337  m_NumGroups = 2;
339  m_pCOMs = new AVector[m_NumGroups];
340 }
341 
342 
344 //--------------------------------------------------------------------
345 // print the distance for this distance restraint
346 //--------------------------------------------------------------------
347  double Distance;
348  char Str1[20], Str2[20];
349 
350  Distance = m_pCOMs[0].Dist(m_pCOMs[1]);
351  sprintf(Str1, "%7.3f", Distance);
352  Distance = GetDistTarget();
353  sprintf(Str2, "%7.3f", Distance);
354 
355 #if defined(_VERBOSE_PMF)
356  iout << "Distance = ";
357  iout << Str1;
358  iout << " Target = ";
359  iout << Str2;
360  iout << std::endl << endi;
361 #else
362  iout << Str1;
363  iout << " ";
364  iout << Str2;
365  iout << " | ";
366 #endif
367 }
368 
369 
370 double ADistRestraint::GetE(double RefDist, double LambdaKf) {
371 //---------------------------------------------------------------------------
372 // calculate and return the Energy for this distance restraint.
373 //
374 // E = (Kf/2) * (di-dref)**2
375 //
376 // where di is the distance between 2 centers-of-mass of restrained atoms,
377 // and dref is the reference distance.
378 //
379 // Note: COM's are calculated before this routine is called.
380 //---------------------------------------------------------------------------
381  double Dist, Diff;
382 
383  Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
384  Diff = Dist - RefDist;
385  return( ((m_Kf*LambdaKf)/2.0) * (Diff*Diff) );
386 }
387 
388 
390  double RefDist, double LambdaKf) {
391 //---------------------------------------------------------------------------
392 // calculate and return the gradient for this distance restraint.
393 //
394 // E = (Kf/2) * (di-dref)**2
395 //
396 // return: grad(E)
397 //
398 // Notes: COM is calculated before this routine is called.
399 // m_pCOMS[0 & 1] reference the COM's of each group of atoms
400 //---------------------------------------------------------------------------
401  double Dist;
402  AVector Vec, A, B;
403 
404  ASSERT( (WhichGroup==0) || (WhichGroup==1) );
405  if (WhichGroup == 0) {
406  A = m_pCOMs[0];
407  B = m_pCOMs[1];
408  }
409  else {
410  A = m_pCOMs[1];
411  B = m_pCOMs[0];
412  }
413 
414  Dist = A.Dist(B);
415  Vec.Set(A[0]-B[0], A[1]-B[1], A[2]-B[2]);
416  Vec *= m_Kf * LambdaKf * (Dist - RefDist) / Dist;
417  return(Vec);
418 }
419 
420 
422 //-----------------------------------------------------------------------
423 // each AnAngleRestraint restrains the angle between 3 groups of atoms
424 //-----------------------------------------------------------------------
425  m_NumGroups = 3;
427  m_pCOMs = new AVector[m_NumGroups];
428 }
429 
430 
432 //--------------------------------------------------------------------
433 // print the angle for this angle restraint
434 //--------------------------------------------------------------------
435  double Angle;
436  char Str1[20], Str2[20];
437 
438  Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]) * (180/kPi);
439  sprintf(Str1, "%8.3f", Angle);
440  Angle = GetAngleTarget() * (180/kPi);
441  sprintf(Str2, "%8.3f", Angle);
442 
443 #if defined(_VERBOSE_PMF)
444  iout << "Angle = ";
445  iout << Str1;
446  iout << " degrees";
447  iout << " Target = ";
448  iout << Str2;
449  iout << " degrees";
450  iout << std::endl << endi;
451 #else
452  iout << Str1;
453  iout << " ";
454  iout << Str2;
455  iout << " | ";
456 #endif
457 }
458 
459 
460 double AnAngleRestraint::GetE(double RefAngle, double LambdaKf) {
461 //---------------------------------------------------------------------------
462 // calculate and return the Energy for this angle restraint.
463 //
464 // E = (Kf/2) * (Theta-ThetaRef)**2
465 //
466 // where Theta is the angle between 3 centers-of-mass of restrained atoms,
467 // m_pCOMs[0] -- m_pCOMs[1] -- m_pCOMs[2].
468 // ThetaRef is the reference angle.
469 //
470 // Note: COM's are calculated before this routine is called.
471 //---------------------------------------------------------------------------
472  double Angle, Diff;
473 
474  Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
475  Diff = Angle - RefAngle;
476  return( ((m_Kf*LambdaKf)/2.0) * (Diff*Diff) );
477 }
478 
479 
481  double RefAngle, double LambdaKf) {
482 //---------------------------------------------------------------------------
483 // calculate and return the gradient for this angle restraint.
484 //
485 // E = (Kf/2) * (Theta-ThetaRef)**2
486 //
487 // return: grad(E)
488 //
489 // Notes: COM's are calculated before this routine is called.
490 //---------------------------------------------------------------------------
491  AVector A, B, C;
492  double Angle;
493  double a, b, c;
494  double u;
495  double Const1, Const2, Const3, Const4, Const5;
496  AVector Vec1, Vec2;
497 
498  ASSERT( (WhichGroup==0) || (WhichGroup==1) || (WhichGroup==2) );
499  A = m_pCOMs[0];
500  B = m_pCOMs[1];
501  C = m_pCOMs[2];
502 
503  a = B.Dist(C);
504  b = A.Dist(C);
505  c = A.Dist(B);
506 
507  u = (a*a + c*c - b*b) / (2.0*a*c);
508 
509  // protect against acos(<-1.0), acos(>1.0), sqrt(<0), and divide by 0
510  if (u < -1.0) {u = -1.0;}
511  if (u > kAlmostOne) {u = kAlmostOne;}
512  Angle = acos(u);
513 
514  Const1 = -1.0 / sqrt(1.0 - u*u);
515  Const2 = m_Kf * LambdaKf * (Angle-RefAngle) * Const1;
516 
517  Const3 = -a/(2.0*c*c) + 1.0/(a+a) + (b*b)/(2.0*a*c*c);
518  Const4 = -c/(2.0*a*a) + 1.0/(c+c) + (b*b)/(2.0*c*a*a);
519  Const5 = -b/(a*c);
520 
521  if (WhichGroup == 0) {
522  Vec1 = (A-C) * (Const5/b);
523  Vec2 = (A-B) * (Const3/c);
524  }
525  else if (WhichGroup == 1) {
526  Vec1 = (B-A) * (Const3/c);
527  Vec2 = (B-C) * (Const4/a);
528  }
529  else {
530  Vec1 = (C-A) * (Const5/b);
531  Vec2 = (C-B) * (Const4/a);
532  }
533  return( (Vec1+Vec2)*Const2);
534 }
535 
536 
538 //----------------------------------------------------------------------------
539 // each ADiheRestraint restrains the dihedral angle between 4 groups of atoms
540 //----------------------------------------------------------------------------
541  m_NumGroups = 4;
543  m_pCOMs = new AVector[m_NumGroups];
544 }
545 
546 
548 //--------------------------------------------------------------------
549 // print the dihedral angle for this dihedral restraint
550 //--------------------------------------------------------------------
551  double Dihedral;
552  char Str1[20], Str2[20], Str3[20];
553 
554  Dihedral = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]) * (180/kPi);
555  sprintf(Str1, "%8.3f", Dihedral);
556  Dihedral = GetDiheTarget1() * (180/kPi);
557  sprintf(Str2, "%8.3f", Dihedral);
558  Dihedral = GetDiheTarget2() * (180/kPi);
559  sprintf(Str3, "%8.3f", Dihedral);
560 
561 #if defined(_VERBOSE_PMF)
562  iout << "Dihedral = ";
563  iout << Str1;
564  iout << " degrees";
565  iout << " Target = ";
566  iout << Str2;
567  iout << " degrees";
568  if (TwoTargets()) {
569  iout << " to ";
570  iout << Str3;
571  iout << " degrees";
572  }
573  iout << std::endl << endi;
574 #else
575  iout << Str1;
576  iout << " ";
577  iout << Str2;
578  if (TwoTargets()) {
579  iout << ", ";
580  iout << Str3;
581  }
582  iout << " | ";
583 #endif
584 }
585 
586 
587 double ADiheRestraint::GetE(double RefAngle, double Const) {
588 //---------------------------------------------------------------------------
589 // calculate and return the Energy for this angle restraint.
590 //
591 // E = (E0/2) * (1 - cos(Chi - ChiRef))
592 //
593 // where Chi is the dihedral angle between 4 centers-of-mass of restrained atoms,
594 // m_pCOMs[0] -- m_pCOMs[1] -- m_pCOMs[2] -- m_pCOMs[3].
595 // ChiRef is the reference angle.
596 //
597 // Note: COM's are calculated before this routine is called.
598 //---------------------------------------------------------------------------
599  double Angle;
600 
601  Angle = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
602  return( (Const/2.0) * (1.0 - cos(Angle-RefAngle)) );
603 }
604 
605 
607  double RefAngle, double Const) {
608 //---------------------------------------------------------------------------
609 // calculate and return the gradient for this dihedral angle restraint.
610 //
611 // E = (E0/2) * (1 - cos(Chi - ChiRef))
612 //
613 // return: grad(E)
614 //
615 // Notes: COM's are calculated before this routine is called.
616 //---------------------------------------------------------------------------
617  AVector A, B, C, D;
618 
619  ASSERT((WhichGroup==0)||(WhichGroup==1)||(WhichGroup==2)||(WhichGroup==3));
620 
621  if ((WhichGroup==0) || (WhichGroup==1)) {
622  A = m_pCOMs[0];
623  B = m_pCOMs[1];
624  C = m_pCOMs[2];
625  D = m_pCOMs[3];
626  }
627  // re-state the problem so the gradient is solved for either atoms 0 or 1
628  else {
629  A = m_pCOMs[3];
630  B = m_pCOMs[2];
631  C = m_pCOMs[1];
632  D = m_pCOMs[0];
633  if (WhichGroup==3) {WhichGroup=0;}
634  if (WhichGroup==2) {WhichGroup=1;}
635  }
636 
637  AVector CD(D - C);
638  AVector CB(B - C);
639  AVector BC(C - B);
640  AVector BA(A - B);
641  AVector AC(C - A);
642  AVector CDxCB, BCxBA;
643  AVector Vec;
644  double phi;
645  double top, bot, u;
646 
647  CDxCB = CD.cross(CB);
648  BCxBA = BC.cross(BA);
649 
650  top = CDxCB.dot(BCxBA);
651  bot = CDxCB.Dist() * BCxBA.Dist();
652 
653  u = top/bot;
654  // protect against acos(<-1.0), acos(>1.0), sqrt(<0), and divide by 0
655  if (u < kAlmostMinusOne) {u = kAlmostMinusOne;}
656  if (u > kAlmostOne) {u = kAlmostOne;}
657 
658  // get dihedral using atan
659  phi = GetDihe(A,B,C,D);
660 
661  AVector dP1, dP2, dP3;
662  AVector dP4, dP5, dP6;
663  ASSERT((WhichGroup==0) || (WhichGroup==1));
664  if (WhichGroup==0) {
665  dP1.Set( 0, 0, 0 );
666  dP2.Set( 0, 0, 0 );
667  dP3.Set( 0, 0, 0 );
668  dP4.Set( 0, -BC[2], BC[1]);
669  dP5.Set( BC[2], 0, -BC[0]);
670  dP6.Set(-BC[1], BC[0], 0 );
671  }
672  else {
673  dP1.Set( 0, -CD[2], CD[1]);
674  dP2.Set( CD[2], 0, -CD[0]);
675  dP3.Set(-CD[1], CD[0], 0 );
676  dP4.Set( 0, AC[2], -AC[1]);
677  dP5.Set(-AC[2], 0, AC[0]);
678  dP6.Set( AC[1], -AC[0], 0 );
679  }
680 
681  Vec = gradU(CDxCB, BCxBA, dP1, dP2, dP3, dP4, dP5, dP6);
682  Vec *= (Const/2.0) * sin(phi-RefAngle) * (-1.0/sqrt(1.0 - u*u));
683 
684  // flip gradient for negative angles
685  if (phi < 0) {
686  Vec *= -1.0;
687  }
688 
689  return(Vec);
690 }
691 
692 
694  AVector& dP1, AVector& dP2, AVector& dP3,
695  AVector& dP4, AVector& dP5, AVector& dP6) {
696 //----------------------------------------------------------------
697 // calculate the gradient for ((P1P2P3.dot.P4P5P6)/(mag(P1P2P3)*mag(P4P5P6)))
698 // P1P2P3 = (P1)i + (P2)j + (P3)k
699 // P4P5P6 = (P4)i + (P5)j + (P6)k
700 // dP1 = (d(P1)/dx)i + (d(P1)/dy)j +(d(P1)/dz)k
701 // dP2 = (d(P2)/dx)i + (d(P2)/dy)j +(d(P2)/dz)k
702 // dP3 = (d(P3)/dx)i + (d(P3)/dy)j +(d(P3)/dz)k
703 // dP4 = (d(P4)/dx)i + (d(P4)/dy)j +(d(P4)/dz)k
704 // dP5 = (d(P5)/dx)i + (d(P5)/dy)j +(d(P5)/dz)k
705 // dP6 = (d(P6)/dx)i + (d(P6)/dy)j +(d(P6)/dz)k
706 //----------------------------------------------------------------
707  double Mag123, Mag456, Dot;
708  double Const1, Const2, Const3;
709  double P1, P2, P3, P4, P5, P6;
710  AVector RetVec;
711 
712  P1 = P1P2P3[0]; P2 = P1P2P3[1]; P3 = P1P2P3[2];
713  P4 = P4P5P6[0]; P5 = P4P5P6[1]; P6 = P4P5P6[2];
714 
715  Mag123 = P1P2P3.Dist();
716  Mag456 = P4P5P6.Dist();
717  Dot = P1P2P3.dot(P4P5P6);
718 
719  Const1 = 1.0 / (Mag123*Mag456);
720  Const2 = -Dot * (1.0 / (Mag123*Mag456*Mag456*Mag456));
721  Const3 = -Dot * (1.0 / (Mag456*Mag123*Mag123*Mag123));
722 
723  RetVec = (dP4*P1 + dP1*P4 + dP5*P2 + dP2*P5 + dP6*P3 + dP3*P6) * Const1 +
724  (dP4*P4 + dP5*P5 + dP6*P6) * Const2 +
725  (dP1*P1 + dP2*P2 + dP3*P3) * Const3;
726 
727  return(RetVec);
728 }
729 
730 
732 //--------------------------------------------------------------------
733 // return the Energy for this fixed position restraint.
734 //--------------------------------------------------------------------
735  return(GetE(m_RefPos));
736 }
737 
738 
740 //-------------------------------------------------------------------------
741 // return the Gradient for this fixed position restraint.
742 //-------------------------------------------------------------------------
743  return(GetGrad(WhichGroup, m_RefPos));
744 }
745 
746 
748 //--------------------------------------------------------------------
749 // calculate and return the Energy for this bound position restraint.
750 //
751 // Note: This is an exception because the form for the E term is
752 // different from the other postion restraints.
753 //--------------------------------------------------------------------
754  double E, Dist, Diff;
755 
756  E = 0.0;
757  Dist = m_pCOMs[0].Dist(m_RefPos);
758  if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
759  ((m_Bound==kLower) && (Dist<m_RefDist))) {
760  Diff = Dist - m_RefDist;
761  E = (m_Kf/2.0) * (Diff*Diff);
762  }
763  return(E);
764 }
765 
766 
768 //---------------------------------------------------------------------------
769 // calculate and return the gradient for this bound position restraint.
770 //
771 // Note: This is an exception because the form for the E term is
772 // different from the other postion restraints.
773 //---------------------------------------------------------------------------
774  double Dist;
775  AVector Vec; // Vec is initialized to (0,0,0)
776 
777  // WhichGroup = 0; // don't care -- there's only 1 atom restrained
778  Dist = m_pCOMs[0].Dist(m_RefPos);
779  if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
780  ((m_Bound==kLower) && (Dist<m_RefDist))) {
781  Vec.Set(m_pCOMs[0][0] - m_RefPos[0],
782  m_pCOMs[0][1] - m_RefPos[1],
783  m_pCOMs[0][2] - m_RefPos[2]);
784  Vec *= m_Kf * (Dist - m_RefDist) / Dist;
785  }
786  return(Vec);
787 }
788 
789 
791 //--------------------------------------------------------------------
792 // return the Energy for this forcing position restraint.
793 //
794 // rref = lambda*r1 + (1-lambda)*r0.
795 // where r0 is the starting position and r1 is the final position
796 //--------------------------------------------------------------------
797  AVector RefPos;
798 
799  RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
800  return(GetE(RefPos, m_LambdaKf));
801 }
802 
803 
805 //---------------------------------------------------------------------------
806 // return the gradient for this forcing position restraint.
807 //
808 // rref = lambda*r1 + (1-lambda)*r0.
809 // where r0 is the starting position and r1 is the final position
810 //---------------------------------------------------------------------------
811  AVector RefPos;
812 
813  RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
814  return(GetGrad(WhichGroup, RefPos, m_LambdaKf));
815 }
816 
817 
819 //---------------------------------------------------------------------------
820 // return dU/dLambda for this forcing position restraint
821 //---------------------------------------------------------------------------
822  AVector RefPos;
823  double T1, T2, T3;
824 
825  RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
826  T1 = (m_pCOMs[0][0] - RefPos[0]) * (m_StartPos[0] - m_StopPos[0]);
827  T2 = (m_pCOMs[0][1] - RefPos[1]) * (m_StartPos[1] - m_StopPos[1]);
828  T3 = (m_pCOMs[0][2] - RefPos[2]) * (m_StartPos[2] - m_StopPos[2]);
829  return( m_Kf * m_LambdaKf * (T1+T2+T3) );
830 }
831 
832 
834 //---------------------------------------------------------------------------
835 // return the Energy for this fixed distance restraint.
836 //---------------------------------------------------------------------------
837  return(GetE(m_RefDist));
838 }
839 
840 
842 //---------------------------------------------------------------------------
843 // return the gradient for this fixed distance restraint.
844 //---------------------------------------------------------------------------
845  return(GetGrad(WhichGroup, m_RefDist));
846 }
847 
848 
850 //---------------------------------------------------------------------------
851 // return the Energy for this bound distance restraint.
852 //---------------------------------------------------------------------------
853  double Dist, E;
854 
855  E = 0.0;
856  Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
857  if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
858  ((m_Bound==kLower) && (Dist<m_RefDist))) {
859  E = GetE(m_RefDist);
860  }
861  return(E);
862 }
863 
864 
866 //---------------------------------------------------------------------------
867 // return the gradient for this bound distance restraint.
868 //---------------------------------------------------------------------------
869  double Dist;
870  AVector Vec;
871 
872  Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
873  if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
874  ((m_Bound==kLower) && (Dist<m_RefDist))) {
875  Vec = GetGrad(WhichGroup, m_RefDist);
876  }
877  return(Vec);
878 }
879 
880 
882 //---------------------------------------------------------------------------
883 // return the Energy for this forcing distance restraint.
884 //---------------------------------------------------------------------------
885  double RefDist;
886 
887  RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
888  return(GetE(RefDist, m_LambdaKf));
889 }
890 
891 
893 //---------------------------------------------------------------------------
894 // return the gradient for this forcing distance restraint.
895 //---------------------------------------------------------------------------
896  double RefDist;
897 
898  RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
899  return(GetGrad(WhichGroup, RefDist, m_LambdaKf));
900 }
901 
902 
904 //---------------------------------------------------------------------------
905 // return dU/dLambda for this forcing distance restraint
906 //---------------------------------------------------------------------------
907  double Dist;
908  double RefDist;
909 
910  Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
911  RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
912  return( m_Kf * m_LambdaKf * (Dist-RefDist)*(m_StartDist-m_StopDist) );
913 }
914 
915 
917 //---------------------------------------------------------------------------
918 // return the Energy for this fixed angle restraint.
919 //---------------------------------------------------------------------------
920  return(GetE(m_RefAngle));
921 }
922 
923 
925 //---------------------------------------------------------------------------
926 // return the gradient for this fixed angle restraint.
927 //---------------------------------------------------------------------------
928  return(GetGrad(WhichGroup, m_RefAngle));
929 }
930 
931 
933 //---------------------------------------------------------------------------
934 // return the Energy for this bound angle restraint.
935 //---------------------------------------------------------------------------
936  double E, Angle;
937 
938  E = 0.0;
939  Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
940  if (((m_Bound==kUpper) && (Angle>m_RefAngle)) ||
941  ((m_Bound==kLower) && (Angle<m_RefAngle))) {
942  E = GetE(m_RefAngle);
943  }
944  return(E);
945 }
946 
947 
949 //---------------------------------------------------------------------------
950 // return the gradient for this bound angle restraint
951 //---------------------------------------------------------------------------
952  double Angle;
953  AVector Vec;
954 
955  Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
956  if (((m_Bound==kUpper) && (Angle>m_RefAngle)) ||
957  ((m_Bound==kLower) && (Angle<m_RefAngle))) {
958  Vec = GetGrad(WhichGroup, m_RefAngle);
959  }
960  return(Vec);
961 }
962 
963 
965 //---------------------------------------------------------------------------
966 // return the Energy for this forcing angle restraint.
967 //---------------------------------------------------------------------------
968  double RefAngle;
969 
970  RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
971  return(GetE(RefAngle, m_LambdaKf));
972 }
973 
974 
976 //---------------------------------------------------------------------------
977 // return the gradient for this forcing angle restraint.
978 //---------------------------------------------------------------------------
979  double RefAngle;
980 
981  RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
982  return(GetGrad(WhichGroup, RefAngle, m_LambdaKf));
983 }
984 
985 
987 //---------------------------------------------------------------------------
988 // return dU/dLambda for this forcing angle restraint
989 //---------------------------------------------------------------------------
990  double Angle;
991  double RefAngle;
992 
993  Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
994  RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
995  return( m_Kf * m_LambdaKf * (Angle-RefAngle)*(m_StartAngle-m_StopAngle) );
996 }
997 
998 
1000 //---------------------------------------------------------------------------
1001 // return the Energy for this fixed dihedral angle restraint.
1002 //---------------------------------------------------------------------------
1003  return(GetE(m_RefAngle, m_Kf));
1004 }
1005 
1006 
1008 //---------------------------------------------------------------------------
1009 // return the gradient for this fixed dihedral angle restraint.
1010 //---------------------------------------------------------------------------
1011  return(GetGrad(WhichGroup, m_RefAngle, m_Kf));
1012 }
1013 
1014 
1016 //---------------------------------------------------------------------------
1017 // return the Energy for this bound dihedral angle restraint.
1018 //---------------------------------------------------------------------------
1019  double E, Dihe, Const;
1020 
1021  Const = m_Kf / (1.0 - cos(m_IntervalAngle));
1022  Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
1023  // dihedral angle is between LowerAngle and UpperAngle
1024  if ( (Dihe>m_LowerAngle) && (Dihe<m_UpperAngle) ) {
1025  E = 0.0;
1026  }
1027  // dihedral angle is between LowerAngle and LowerAngle-IntervalAngle
1028  else if ( (Dihe<m_LowerAngle) && (Dihe>(m_LowerAngle-m_IntervalAngle)) ) {
1029  E = GetE(m_LowerAngle, Const);
1030  }
1031  // dihedral angle is between UpperAngle and UpperAngle+IntervalAngle
1032  else if ( (Dihe>m_UpperAngle) && (Dihe<(m_UpperAngle+m_IntervalAngle)) ) {
1033  E = GetE(m_UpperAngle, Const);
1034  }
1035  // dihedral angle is more than UpperAngle or less than LowerAngle
1036  else {
1037  E = Const;
1038  }
1039  return(E);
1040 }
1041 
1042 
1044 //---------------------------------------------------------------------------
1045 // return the gradient for this bound dihedral angle restraint.
1046 //---------------------------------------------------------------------------
1047  AVector Vec;
1048  double Dihe, Const;
1049 
1050  Const = m_Kf / (1.0 - cos(m_IntervalAngle));
1051  Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
1052  // dihedral angle is between LowerAngle and LowerAngle-IntervalAngle
1053  if ( (Dihe<m_LowerAngle) && (Dihe>(m_LowerAngle-m_IntervalAngle)) ) {
1054  Vec = GetGrad(WhichGroup, m_LowerAngle, Const);
1055  }
1056  // dihedral angle is between UpperAngle and UpperAngle+IntervalAngle
1057  else if ( (Dihe>m_UpperAngle) && (Dihe<(m_UpperAngle+m_IntervalAngle)) ) {
1058  Vec = GetGrad(WhichGroup, m_UpperAngle, Const);
1059  }
1060  return(Vec);
1061 }
1062 
1063 
1065 //---------------------------------------------------------------------------
1066 // return the Energy for this forcing dihedral angle restraint.
1067 //---------------------------------------------------------------------------
1068  double RefDihe;
1069 
1070  RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
1071  return(GetE(RefDihe, m_Kf*m_LambdaKf));
1072 }
1073 
1074 
1076 //---------------------------------------------------------------------------
1077 // return the gradient for this forcing dihedral angle restraint.
1078 //---------------------------------------------------------------------------
1079  double RefDihe;
1080 
1081  RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
1082  return(GetGrad(WhichGroup, RefDihe, m_Kf*m_LambdaKf));
1083 }
1084 
1085 
1087 //---------------------------------------------------------------------------
1088 // return dU/dLambda for this forcing dihedral angle restraint
1089 //---------------------------------------------------------------------------
1090  double Dihe;
1091  double RefDihe;
1092 
1093  Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
1094  RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
1095  return((m_Kf/2)*m_LambdaKf * sin(Dihe-RefDihe) * (m_StartAngle-m_StopAngle));
1096 }
1097 
struct angle Angle
Bool_t Set(double x, double y, double z)
int getPosition(int atomid, Position &position)
const double kAlmostMinusOne
void SetGroup(AGroup &Group, int GroupIndex)
AVector cross(const AVector &Vector)
AVector GetGrad(int WhichGroup, double RefDihe, double Const)
int AtomID
Definition: NamdTypes.h:29
AVector gradU(AVector &P1P2P3, AVector &P4P5P6, AVector &dP1, AVector &dP2, AVector &dP3, AVector &dP4, AVector &dP5, AVector &dP6)
double GetE(double RefAngle, double LambdaKf=1.0)
AVector GetGradient(int WhichGroup)
const BigReal A
Definition: Vector.h:64
int GetNumInGroup()
struct dihedral Dihedral
double GetDihe(AVector &A, AVector &B, AVector &C, AVector &D)
#define iout
Definition: InfoStream.h:87
virtual ~ARestraint()
int addForce(int atomid, Force force)
double GetE(double RefDist, double LambdaKf=1.0)
double GetE(AVector RefPos, double LambdaKf=1.0)
virtual double GetDiheTarget1()=0
int Index
Definition: structures.h:26
virtual double GetDistTarget()=0
double dot(const AVector &Vector)
AVector GetGradient(int WhichGroup)
double GetAngle(AVector &A, AVector &B, AVector &C)
#define ASSERT(E)
AVector GetGrad(int WhichGroup, double RefAngle, double LambdaKf=1.0)
AVector GetGradient(int WhichGroup)
virtual double GetDiheTarget2()=0
void DistributeForce(int WhichGroup, AVector Force, GlobalMasterFreeEnergy &CFE)
AVector GetGradient(int WhichGroup)
virtual Bool_t TwoTargets()=0
AVector GetGradient(int WhichGroup)
virtual AVector GetPosTarget()=0
const double kPi
AVector GetGradient(int WhichGroup)
void NAMD_die(const char *err_msg)
Definition: common.C:83
static double m_LambdaRef
AVector GetGrad(int WhichGroup, AVector RefPos, double LambdaKf=1.0)
AVector GetGradient(int WhichGroup)
AVector GetGrad(int WhichGroup, double RefDist, double LambdaKf=1.0)
void UpdateCOMs(GlobalMasterFreeEnergy &CFE)
void SetGroups(AGroup &Group1)
AVector GetGradient(int WhichGroup)
AVector GetGradient(int WhichGroup)
__device__ __forceinline__ float dot(const float3 v1, const float3 v2)
static double m_LambdaKf
void EarlyExit(const char *Str, int AtomID)
AVector GetGradient(int WhichGroup)
void SetEqual(AVector &Vec1, const Vector &Vec2)
double GetE(double RefDihe, double Const)
const BigReal B
double Dist()
float Mass
Definition: ComputeGBIS.inl:20
infostream & endi(infostream &s)
Definition: InfoStream.C:38
AVector GetGradient(int WhichGroup)
AVector GetGradient(int WhichGroup)
const double kAlmostOne
virtual double GetDistance()=0
AGroup * m_pGroups
virtual double GetAngleTarget()=0
AVector * m_pCOMs