86926343eb76e067c52f8d1bdda5a4231117c19c
[namd.git] / src / ComputeBonds.C
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
7 #include "InfoStream.h"
8 #include "ComputeBonds.h"
9 #include "Molecule.h"
10 #include "Parameters.h"
11 #include "Node.h"
12 #include "ReductionMgr.h"
13 #include "Lattice.h"
14 #include "PressureProfile.h"
15 #include "Debug.h"
16
17
18 // static initialization
19 int BondElem::pressureProfileSlabs = 0;
20 int BondElem::pressureProfileAtomTypes = 1;
21 BigReal BondElem::pressureProfileThickness = 0;
22 BigReal BondElem::pressureProfileMin = 0;
23
24 void BondElem::getMoleculePointers
25     (Molecule* mol, int* count, int32*** byatom, Bond** structarray)
26 {
27 #ifdef MEM_OPT_VERSION
28   NAMD_die("Should not be called in BondElem::getMoleculePointers in memory optimized version!");
29 #else
30   *count = mol->numBonds;
31   *byatom = mol->bondsByAtom;
32   *structarray = mol->bonds;
33 #endif
34 }
35
36 void BondElem::getParameterPointers(Parameters *p, const BondValue **v) {
37   *v = p->bond_array;
38 }
39
40 void BondElem::computeForce(BondElem *tuples, int ntuple, BigReal *reduction, 
41                             BigReal *pressureProfileData)
42 {
43  SimParameters *const simParams = Node::Object()->simParameters;
44  const Lattice & lattice = tuples[0].p[0]->p->lattice;
45
46  //fepb BKR
47  const int step = tuples[0].p[0]->p->flags.step;
48  const BigReal alchLambda = simParams->getCurrentLambda(step);
49  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
50  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
51  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
52  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
53  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
54  Molecule *const mol = Node::Object()->molecule;
55  //fepe
56
57  for ( int ituple=0; ituple<ntuple; ++ituple ) {
58   const BondElem &tup = tuples[ituple];
59   enum { size = 2 };
60   const AtomID (&atomID)[size](tup.atomID);
61   const int    (&localIndex)[size](tup.localIndex);
62   TuplePatchElem * const(&p)[size](tup.p);
63   const Real (&scale)(tup.scale);
64   const BondValue * const(&value)(tup.value);
65
66   DebugM(1, "::computeForce() localIndex = " << localIndex[0] << " "
67                << localIndex[1] << std::endl);
68
69   // skip Lonepair bonds (other k=0. bonds have been filtered out)
70   if (0. == value->k) continue;
71
72   BigReal scal;         // force scaling
73   BigReal energy;       // energy from the bond
74
75   //DebugM(3, "::computeForce() -- starting with bond type " << bondType << std::endl);
76
77   // get the bond information
78   Real k = value->k * scale;
79   Real x0 = value->x0;
80
81   // compute vectors between atoms and their distances
82   const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
83                                         p[1]->x[localIndex[1]].position);
84
85   if (0. == x0) {
86     BigReal r2 = r12.length2();
87     scal = -2.0*k;    // bond interaction for equilibrium length 0
88     energy = k*r2;
89     // TODO: Instead flag by mass for Drude particles, otherwise mixing and 
90     // matching Drude and alchemical "twin" particles will likely not work as 
91     // expected.
92     if(simParams->drudeOn) { // for Drude bonds
93       BigReal drudeBondLen = simParams->drudeBondLen;
94       BigReal drudeBondConst = simParams->drudeBondConst;
95       if ( (drudeBondConst > 0) && ( ! simParams->drudeHardWallOn ) &&
96           (r2 > drudeBondLen*drudeBondLen) ) {
97         // add a quartic restraining potential to keep Drude bond short
98         BigReal r = sqrt(r2);
99         BigReal diff = r - drudeBondLen;
100         BigReal diff2 = diff*diff;
101
102         scal += -4*drudeBondConst * diff2 * diff / r;
103         energy += drudeBondConst * diff2 * diff2;
104       }
105     }
106   }
107   else {
108     BigReal r = r12.length();  // Distance between atoms
109     BigReal diff = r - x0;     // Compare it to the rest bond
110
111     //  Add the energy from this bond to the total energy
112     energy = k*diff*diff;
113
114     //  Determine the magnitude of the force
115     diff *= -2.0*k;
116
117     //  Scale the force vector accordingly
118     scal = (diff/r);
119   }
120
121   //fepb - BKR scaling of alchemical bonded terms
122   //       NB: TI derivative is the _unscaled_ energy.
123   if ( simParams->alchOn ) {
124     switch ( mol->get_fep_bonded_type(atomID, 2) ) {
125     case 1:
126       reduction[bondEnergyIndex_ti_1] += energy;
127       reduction[bondEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
128       energy *= bond_lambda_1;
129       scal *= bond_lambda_1;
130       break;
131     case 2:
132       reduction[bondEnergyIndex_ti_2] += energy;
133       reduction[bondEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
134       energy *= bond_lambda_2;
135       scal *= bond_lambda_2;
136       break;
137     }
138   }
139   //fepe
140   const Force f12 = scal * r12;
141
142
143   //  Now add the forces to each force vector
144   p[0]->f[localIndex[0]] += f12;
145   p[1]->f[localIndex[1]] -= f12;
146
147   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
148   reduction[bondEnergyIndex] += energy;
149   reduction[virialIndex_XX] += f12.x * r12.x;
150   reduction[virialIndex_XY] += f12.x * r12.y;
151   reduction[virialIndex_XZ] += f12.x * r12.z;
152   reduction[virialIndex_YX] += f12.y * r12.x;
153   reduction[virialIndex_YY] += f12.y * r12.y;
154   reduction[virialIndex_YZ] += f12.y * r12.z;
155   reduction[virialIndex_ZX] += f12.z * r12.x;
156   reduction[virialIndex_ZY] += f12.z * r12.y;
157   reduction[virialIndex_ZZ] += f12.z * r12.z;
158
159   if (pressureProfileData) {
160     BigReal z1 = p[0]->x[localIndex[0]].position.z;
161     BigReal z2 = p[1]->x[localIndex[1]].position.z;
162     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
163     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
164     pp_clamp(n1, pressureProfileSlabs);
165     pp_clamp(n2, pressureProfileSlabs);
166     int p1 = p[0]->x[localIndex[0]].partition;
167     int p2 = p[1]->x[localIndex[1]].partition;
168     int pn = pressureProfileAtomTypes;
169     pp_reduction(pressureProfileSlabs,
170                 n1, n2, 
171                 p1, p2, pn, 
172                 f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
173                 pressureProfileData);
174   } 
175
176  }
177 }
178
179 void BondElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
180 {
181   reduction->item(REDUCTION_BOND_ENERGY) += data[bondEnergyIndex];
182   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[bondEnergyIndex_f];
183   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[bondEnergyIndex_ti_1];
184   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[bondEnergyIndex_ti_2];
185   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
186 }
187