Update default behavior and documentation for Drude
[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 && !simParams->drudeHardWallOn ) { // for Drude bonds
93       BigReal drudeBondLen = simParams->drudeBondLen;
94       BigReal drudeBondConst = simParams->drudeBondConst;
95       if ( drudeBondConst > 0 && r2 > drudeBondLen*drudeBondLen ) {
96         // add a quartic restraining potential to keep Drude bond short
97         BigReal r = sqrt(r2);
98         BigReal diff = r - drudeBondLen;
99         BigReal diff2 = diff*diff;
100
101         scal += -4*drudeBondConst * diff2 * diff / r;
102         energy += drudeBondConst * diff2 * diff2;
103       }
104     }
105   }
106   else {
107     BigReal r = r12.length();  // Distance between atoms
108     BigReal diff = r - x0;     // Compare it to the rest bond
109
110     //  Add the energy from this bond to the total energy
111     energy = k*diff*diff;
112
113     //  Determine the magnitude of the force
114     diff *= -2.0*k;
115
116     //  Scale the force vector accordingly
117     scal = (diff/r);
118   }
119
120   //fepb - BKR scaling of alchemical bonded terms
121   //       NB: TI derivative is the _unscaled_ energy.
122   if ( simParams->alchOn ) {
123     switch ( mol->get_fep_bonded_type(atomID, 2) ) {
124     case 1:
125       reduction[bondEnergyIndex_ti_1] += energy;
126       reduction[bondEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
127       energy *= bond_lambda_1;
128       scal *= bond_lambda_1;
129       break;
130     case 2:
131       reduction[bondEnergyIndex_ti_2] += energy;
132       reduction[bondEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
133       energy *= bond_lambda_2;
134       scal *= bond_lambda_2;
135       break;
136     }
137   }
138   //fepe
139   const Force f12 = scal * r12;
140
141
142   //  Now add the forces to each force vector
143   p[0]->f[localIndex[0]] += f12;
144   p[1]->f[localIndex[1]] -= f12;
145
146   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
147   reduction[bondEnergyIndex] += energy;
148   reduction[virialIndex_XX] += f12.x * r12.x;
149   reduction[virialIndex_XY] += f12.x * r12.y;
150   reduction[virialIndex_XZ] += f12.x * r12.z;
151   reduction[virialIndex_YX] += f12.y * r12.x;
152   reduction[virialIndex_YY] += f12.y * r12.y;
153   reduction[virialIndex_YZ] += f12.y * r12.z;
154   reduction[virialIndex_ZX] += f12.z * r12.x;
155   reduction[virialIndex_ZY] += f12.z * r12.y;
156   reduction[virialIndex_ZZ] += f12.z * r12.z;
157
158   if (pressureProfileData) {
159     BigReal z1 = p[0]->x[localIndex[0]].position.z;
160     BigReal z2 = p[1]->x[localIndex[1]].position.z;
161     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
162     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
163     pp_clamp(n1, pressureProfileSlabs);
164     pp_clamp(n2, pressureProfileSlabs);
165     int p1 = p[0]->x[localIndex[0]].partition;
166     int p2 = p[1]->x[localIndex[1]].partition;
167     int pn = pressureProfileAtomTypes;
168     pp_reduction(pressureProfileSlabs,
169                 n1, n2, 
170                 p1, p2, pn, 
171                 f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
172                 pressureProfileData);
173   } 
174
175  }
176 }
177
178 void BondElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
179 {
180   reduction->item(REDUCTION_BOND_ENERGY) += data[bondEnergyIndex];
181   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[bondEnergyIndex_f];
182   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[bondEnergyIndex_ti_1];
183   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[bondEnergyIndex_ti_2];
184   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
185 }
186