1 /* -*- c++ -*- */
2 #ifndef ISGLENNARDJONESFORCE_H
3 #define ISGLENNARDJONESFORCE_H
4 
5 #include "GenericTopology.h"
6 #include "LennardJonesParameters.h"
7 #include "ScalarStructure.h"
8 #include "ExclusionTable.h"
9 #include "Parameter.h"
10 #include <string>
11 
12 // uncomment for debugging purposes
13 //#define DEBUG_NORMAL_LJ
14 //#define DEBUG_SOFTCORE
15 //#define DEBUG_LJ_DISTANCE_CHECK
16 
17 namespace ProtoMol {
18 
19   //_________________________________________________________________ ISGLennardJonesForce
20   class iSGLennardJonesForce {
21 
22     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
23     // Constructors, destructors, assignment
24     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
25   public:
iSGLennardJonesForce()26     iSGLennardJonesForce() {};
27 
28     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29     // New methods of class iSGLennardJonesForce
30     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31   public:
operator()32     void operator()(Real &energy,
33 		    Real &force,
34 		    Real &DeltaMu,
35 		    Real DistSquared,
36 		    Real rDistSquared,
37 		    const Vector3D &,
38 		    const GenericTopology* topo,
39 		    int atom1, int atom2,
40 		    ExclusionClass excl) const {
41 
42 #ifdef DEBUG_LJ_DISTANCE_CHECK
43       if(rDistSquared > 1.0/0.25)
44 	Report::report << Report::warning << "iSGLennardJonesForce::operator(): atom "<<atom1<<" and "<<atom2<<" get closer"
45 	       << " than 0.5 angstroms (="<<sqrt(1.0/rDistSquared)<<")." << Report::endr;
46 #endif
47 
48       // Fast LJ computations
49       Real r2   = rDistSquared;
50       Real r6   = r2*r2*r2;
51       Real r12, rij_6;
52 
53       // soft-core interaction variables
54       Real Vr6_old = 0.0;
55       Real Vr6_new = 0.0;
56       Real Vr12_old = 0.0;
57       Real Vr12_new = 0.0;
58       Real Scale_new, Scale_old, Lambda, alpha;
59       Real SoftDist_old, SoftDist_new, LambdaMu;
60 
61       // energy and size parameter variables
62       LennardJonesParameters params, params_old, params_new;
63       Real A_old = 0.0;
64       Real A_new = 0.0;
65       Real B_old = 0.0;
66       Real B_new = 0.0;
67 
68       unsigned int type1_old, type1_new;
69 
70       // LJ parameter bank indices
71       unsigned int type2, i, j;
72 
73       // transforming atom indicators
74       bool atom1_scaled, atom2_scaled;
75       atom1_scaled = atom2_scaled = false;
76       int myStage, OldStage, atom_stage;
77       myStage = OldStage = atom_stage = 0;
78 
79       // Use the molecule types to determine if the atoms
80       // are being transformed or not
81       int M1 = topo->atoms[atom1].molecule;
82       int M2 = topo->atoms[atom2].molecule;
83 
84       // if lambda for either molecule is nonzero then the molecule is being transformed
85       if (topo->molecules[M1].lambda != 0.0) {atom1_scaled = true;}
86       if (topo->molecules[M2].lambda != 0.0) {atom2_scaled = true;}
87 
88       // integer used to select which interaction type to compute
89       int Choice;
90 
91       // determine the value of Choice
92       // neither atom is being transformed
93       if ( !(atom1_scaled) && !(atom2_scaled) ) {Choice = 0;}
94       // intramolecular interaction on a transforming molecule
95       else if ( (atom1_scaled) && (atom2_scaled) ) {
96 
97         // determine the current transformation stage of the molecule
98         myStage = static_cast<int>(floor(topo->molecules[M1].lambda)) + 1;
99 
100         // get the transformation stage #'s for these atoms
101         int atom1_stage = topo->atoms[atom1].stageNumber;
102         int atom2_stage = topo->atoms[atom2].stageNumber;
103 
104         // if the current stage is not the same as either atom1_stage or atom2_stage
105         // then neither atom is currently being transformed
106         if (myStage != atom1_stage && myStage != atom2_stage) Choice = 0;
107         else Choice = 2;
108       }
109       // one of the atoms is being transformed
110       else {
111         // determine the current transformation stage of the molecule
112         // and get the transformation stage # for this atoms
113         if (atom1_scaled) {
114           myStage = static_cast<int>(floor(topo->molecules[M1].lambda)) + 1;
115           atom_stage = topo->atoms[atom1].stageNumber;}
116         else {
117           myStage = static_cast<int>(floor(topo->molecules[M2].lambda)) + 1;
118           atom_stage = topo->atoms[atom2].stageNumber;}
119 
120         // if the current stage is not the same as atom_stage
121         // then neither atom is currently being transformed
122         if (myStage != atom_stage) Choice = 0;
123         else Choice = 1;
124       }
125 
126       //************************************************************************
127       // Compute the appropriate LJ interaction
128       DeltaMu = 0.0;
129 
130       switch (Choice) {
131         //----------------------------------------------------------
132       case 0:
133 	// interaction between two non-transforming atoms
134 
135 	// get the identities of the two atoms
136 	// i, j and are used to retrieve the correct set of
137 	// parameters from the LJ parameter bank
138 	// to access the proper LJ bank, atomI must always be the
139         // atom with the LOWER atom_type number since that is how
140         // the bank was constructed in the buildISGTopology function
141         if (topo->atoms[atom1].type <= topo->atoms[atom2].type) {
142           i = topo->molecules[M1].type;
143           j = topo->molecules[M2].type;}
144         else {
145           i = topo->molecules[M2].type;
146           j = topo->molecules[M1].type;}
147 
148         // check to see if atom1 may have already been successfully transformed in a previous stage
149         if (atom1_scaled) {
150           // this atom has already been transformed if myStage > atom_stage
151           if (myStage > atom_stage) {
152             if (topo->atoms[atom1].type <= topo->atoms[atom2].type) i = topo->molecules[M1].newtype;
153             else j = topo->molecules[M1].newtype;
154           }
155         } // end if statement
156         // check to see if atom2 may have already been successfully transformed in a previous stage
157         if (atom2_scaled) {
158           // this atom has already been transformed if myStage > atom_stage
159           if (myStage > atom_stage) {
160             if (topo->atoms[atom1].type <= topo->atoms[atom2].type) j = topo->molecules[M2].newtype;
161             else i = topo->molecules[M2].newtype;
162           }
163         } // end if statement
164 
165 	// get the LJ parameters for this pair
166 	params = topo->isgLJParms(i,j,topo->atoms[atom1].type,topo->atoms[atom2].type);
167 	A_old = (excl != EXCLUSION_MODIFIED ? params.A : params.A14);
168 	B_old = (excl != EXCLUSION_MODIFIED ? params.B : params.B14);
169 
170 	// if any of these atoms are dummy atoms (epsilon = 0), then skip
171 	if (A_old == 0.0) break;
172 
173 	// attractive LJ term
174 	Vr6_old = B_old * r6;
175 
176 	// repulsive LJ term
177 	r12 = r6 * r6;
178 	Vr12_old = A_old * r12;
179 
180 	// energy and force
181 	energy = Vr12_old - Vr6_old;
182 	force = 12.0*Vr12_old*r2 - 6.0*r2*Vr6_old;
183 
184 #ifdef DEBUG_NORMAL_LJ
185     	Report::report.precision(6);
186 	Report::report << "i = " << atom1+1 << ", j = " << atom2+1 << ", " << i << ", " << j
187                        << ", eps = " << B_old*B_old / (4.0*A_old) << ", energy = "
188                        << energy << Report::endr;
189 #endif
190 	break;
191         //----------------------------------------------------------
192       case 1:
193 	// interaction between non-transforming atom and a transforming atom
194 	// 'old' indicates OldType, 'new' indicates NewType
195 
196 	if (atom1_scaled) {
197 	  // get the lambda factor
198 	  Lambda = topo->molecules[M1].lambda;
199 
200 	  // get the identities of the two atoms
201 	  type1_old = topo->molecules[M1].type;
202 	  type1_new = topo->molecules[M1].newtype;
203 	  type2 = topo->molecules[M2].type;
204 
205           // get the alphaLJ parameter for this atom
206           alpha = topo->atoms[atom1].alphaLJ;
207 	}
208 	else {
209 	  // get the lambda factor
210 	  Lambda = topo->molecules[M2].lambda;
211 
212 	  // get the identities of the two atoms
213 	  type1_old = topo->molecules[M2].type;
214 	  type1_new = topo->molecules[M2].newtype;
215 	  type2 = topo->molecules[M1].type;
216 
217           // get the alphaLJ parameter for this atom
218           alpha = topo->atoms[atom2].alphaLJ;
219 	}
220 
221         // i, j and are used to retrieve the correct set of
222         // parameters from the LJ parameter bank
223         // if atom 1 is the transforming atom...
224         if (atom1_scaled) {
225 
226           // to access the proper LJ bank, atomI must always be the
227           // atom with the LOWER atom_type number -- TIM
228           if (topo->atoms[atom1].type <= topo->atoms[atom2].type) {
229             i = type1_old;
230             j = type2;}
231           else {
232             i = type2;
233             j = type1_old;
234           }
235         }
236         // if atom2 is the transforming atom then the procedure
237         // below is the opposite of the procedure above.  AtomI must still
238         // be the atom with the LOWER atom_type number, but remember that
239         // if atom2 is the transforming atom then the variable type2 is equal
240         // to the molecule_type of atom 1!!!
241         else {
242 
243           // if atom1 has a lower atom_type, then i = type2 since now type2
244           // is equal to the molecule_type of atom1
245           if (topo->atoms[atom1].type <= topo->atoms[atom2].type) {
246             i = type2;
247             j = type1_old;}
248           // if atom2 has the lower atom_type, then i = type1A
249           else {
250             i = type1_old;
251             j = type2;
252           }
253         }
254 
255 
256         // compute the A-B and B-B epsilon and sigmas
257         // retrieve the A & B parameters for these interaction types
258         // get the LJ parameters for this pair when the transforming atom is in the old state
259         params_old = topo->isgLJParms(i,j,topo->atoms[atom1].type,topo->atoms[atom2].type);
260 
261         // determine the new i and j when the transforming atom is in the new state
262         if (atom1_scaled) {
263           if (topo->atoms[atom1].type <= topo->atoms[atom2].type) i = type1_new;
264           else j = type1_new;}
265         else {
266           if (topo->atoms[atom1].type <= topo->atoms[atom2].type) j = type1_new;
267           else i = type1_new;
268         }
269 
270         // get the LJ parameters for this pair when the transforming atom is in the new state
271         params_new = topo->isgLJParms(i,j,topo->atoms[atom1].type,topo->atoms[atom2].type);
272         A_old = (params_old.A);
273         B_old = (params_old.B);
274         A_new = (params_new.A);
275         B_new = (params_new.B);
276 
277         // if the type2 atom is a dummy atom then skip (epsilon = 0)
278         if (A_old == 0.0 && A_new == 0.0) break;
279 
280         // determine the current transformation stage of the molecule
281         OldStage = myStage - 1;
282 
283         // LJ soft-core energy scaling factors
284         LambdaMu = 2.0 * alpha * (Lambda - OldStage) * (myStage - Lambda);
285         Scale_old = alpha * (Lambda - OldStage) * (Lambda - OldStage);
286         Scale_new = alpha * (myStage - Lambda) * (myStage - Lambda);
287 
288         // soft-core interaction distances
289         // if the transforming atom is a dummy atom in one of its identities then
290         // we must make sure not to divide by zero
291         rij_6 = DistSquared * DistSquared * DistSquared;
292         if (A_old != 0.0) SoftDist_old = rij_6 + (A_old / B_old) * Scale_old;
293         else SoftDist_old = 1.0;
294         if (A_new != 0.0) SoftDist_new = rij_6 + (A_new / B_new) * Scale_new;
295         else SoftDist_new = 1.0;
296 
297         // attractive LJ terms
298         Vr6_old = B_old / SoftDist_old;
299         Vr6_new = B_new / SoftDist_new;
300 
301         // repulsive LJ terms
302         Vr12_old = A_old / (SoftDist_old * SoftDist_old);
303         Vr12_new = A_new / (SoftDist_new * SoftDist_new);
304 
305         // energy, force, and chemical potential difference
306         // if the transforming atom is a dummy atom in one of its identities then
307         // we must make sure not to divide by zero
308         if (A_old == 0.0) {
309           A_old = 1.0;
310           B_old = 1.0;}
311         if (A_new == 0.0) {
312           A_new = 1.0;
313           B_new = 1.0;}
314         energy = (myStage - Lambda) * (Vr12_old - Vr6_old) + (Lambda - OldStage) * (Vr12_new - Vr6_new);
315         force = 6.0 * r2 / r6 * ((myStage - Lambda) * Vr12_old * (2.0 * Vr6_old/B_old - B_old/A_old)
316           + (Lambda - OldStage) * Vr12_new * (2.0 * Vr6_new/B_new - B_new/A_new));
317         DeltaMu = (Vr12_new - Vr6_new) - (Vr12_old - Vr6_old)
318          + LambdaMu * (Vr12_new * (2.0 * Vr6_new*(A_new/(B_new*B_new)) - 1.0)
319          - Vr12_old * (2.0 * Vr6_old*(A_old/(B_old*B_old)) - 1.0));
320 
321 #ifdef DEBUG_SOFTCORE
322         Report::report.precision(6);
323         Report::report << Report::warning << "i = " << atom1+1 << ", j = " << atom2+1 << ", " << type1_old << "->" << type1_new
324                        << ", " << type2 << ", rij = " << 1.0/sqrt(r2) << ", Energy = " << energy
325                        << ", force = " << force << ", DeltaMu = " << DeltaMu << Report::endr;
326 #endif //DEBUG_SOFTCORE
327 	break;
328         //----------------------------------------------------------
329       case 2:
330 	// intramolecular interaction on a transforming molecule
331 
332         // Get the lambda factor for each atom
333 	Lambda = topo->molecules[M1].lambda;
334 
335         // Get the alphaLJ parameter for this interaction
336         Real tempAlpha1 = topo->atoms[atom1].alphaLJ;
337         Real tempAlpha2 = topo->atoms[atom2].alphaLJ;
338         if (tempAlpha1 > tempAlpha2) alpha = tempAlpha1;
339         else alpha = tempAlpha2;
340 
341         // retrieve the A-A and B-B energy and size parameters
342         // j is used to retrieve the correct set of
343         // parameters from the LJ parameter bank
344         j = topo->molecules[M1].type;
345 
346         // get the LJ parameters for this pair when the atoms are in the old state
347         params_old = topo->isgLJParms(j,j,topo->atoms[atom1].type,topo->atoms[atom2].type);
348 
349         // determine the proper index # into to the bank when the atoms are in the new state
350         j = topo->molecules[M1].newtype;
351 
352         // get the LJ parameters for this pair when the atoms are in the new state
353         params_new = topo->isgLJParms(j,j,topo->atoms[atom1].type,topo->atoms[atom2].type);
354         A_old = (excl != EXCLUSION_MODIFIED ? params_old.A : params_old.A14);
355         B_old = (excl != EXCLUSION_MODIFIED ? params_old.B : params_old.B14);
356         A_new = (excl != EXCLUSION_MODIFIED ? params_new.A : params_new.A14);
357         B_new = (excl != EXCLUSION_MODIFIED ? params_new.B : params_new.B14);
358 
359         // if both of these atoms are dummy atoms (epsilon = 0), then skip
360         if (A_old == 0.0 && A_new == 0.0) break;
361 
362         // determine the current transformation stage of the molecule
363         OldStage = myStage - 1;
364 
365         // LJ soft-core energy scaling factors
366         LambdaMu = 2.0 * alpha * (Lambda - OldStage) * (myStage - Lambda);
367         Scale_old = alpha * (Lambda - OldStage) * (Lambda - OldStage);
368         Scale_new = alpha * (myStage - Lambda) * (myStage - Lambda);
369 
370         // soft-core interaction distances
371         // if the transforming atom is a dummy atom in one of its identities then
372         // we must make sure not to divide by zero
373         rij_6 = DistSquared * DistSquared * DistSquared;
374         if (A_old != 0.0) SoftDist_old = rij_6 + (A_old / B_old) * Scale_old;
375         else SoftDist_old = 1.0;
376         if (A_new != 0.0) SoftDist_new = rij_6 + (A_new / B_new) * Scale_new;
377         else SoftDist_new = 1.0;
378 
379         // attractive LJ terms
380         Vr6_old = B_old / SoftDist_old;
381         Vr6_new = B_new / SoftDist_new;
382 
383         // repulsive LJ terms
384         Vr12_old = A_old / (SoftDist_old * SoftDist_old);
385         Vr12_new = A_new / (SoftDist_new * SoftDist_new);
386 
387         // get the transformation stage #'s for these atoms
388         int atom1_stage = topo->atoms[atom1].stageNumber;
389         int atom2_stage = topo->atoms[atom2].stageNumber;
390         int max_stage = max(atom1_stage,atom2_stage);
391 
392         // if the transforming atom is a dummy atom in one of its identities then
393         // we must make sure not to divide by zer
394         if (A_old == 0.0) {
395           A_old = 1.0;
396           B_old = 1.0;}
397         if (A_new == 0.0) {
398           A_new = 1.0;
399           B_new = 1.0;}
400 
401         // if the current stage is not the same as the largest atom stage number then
402         // we do not compute any deltaMu, and if the current stage is less than the largest
403         // atom stage number then we only compute interactions in the old identity, if the
404         // current stage is larger than both atom stage numbers then we use only the new identity
405         if (myStage == max_stage) {
406           // energy, force, and chemical potential difference
407           energy = (myStage - Lambda) * (Vr12_old - Vr6_old) + (Lambda - OldStage) * (Vr12_new - Vr6_new);
408           force = 6.0 * r2 * rij_6 * ((myStage - Lambda) * Vr12_old * (2.0 * Vr6_old/B_old - B_old/A_old)
409             + (Lambda - OldStage) * Vr12_new * (2.0 * Vr6_new/B_new - B_new/A_new));
410           DeltaMu = (Vr12_new - Vr6_new) - (Vr12_old - Vr6_old)
411             + LambdaMu * (Vr12_new * (2.0 * Vr6_new*(A_new/(B_new*B_new)) - 1.0)
412             - Vr12_old * (2.0 * Vr6_old*(A_old/(B_old*B_old)) - 1.0));
413         }
414         else if (myStage > max_stage) {
415           // energy, force, and chemical potential difference (new identity)
416           energy = (Lambda - OldStage) * (Vr12_new - Vr6_new);
417           force = 6.0 * r2 * rij_6 * (Lambda - OldStage) * Vr12_new * (2.0 * Vr6_new/B_new - B_new/A_new);
418         }
419         else {
420           // energy, force, and chemical potential difference (old identity)
421           energy = (myStage - Lambda) * (Vr12_old - Vr6_old);
422           force = 6.0 * r2 * rij_6 * (myStage - Lambda) * Vr12_old * (2.0 * Vr6_old/B_old - B_old/A_old);
423         }
424         break;
425 	//************************************************************************
426       } // end switch structure
427 
428     } // end operator()
429 
430 
accumulateEnergy(ScalarStructure * energies,Real energy,Real deltaMu)431     static void accumulateEnergy(ScalarStructure* energies, Real energy, Real deltaMu) {
432       (*energies)[ScalarStructure::LENNARDJONES] += energy;
433       (*energies)[ScalarStructure::LENNARDJONES_DELTAMU] += deltaMu;
434     }
getEnergy(const ScalarStructure * energies)435     static Real getEnergy(const ScalarStructure* energies) {
436       return  (*energies)[ScalarStructure::LENNARDJONES];}
437 
438     // Parsing
getId()439     static std::string getId() {return keyword;}
getParameterSize()440     static unsigned int getParameterSize() {return 0;}
getParameters(std::vector<Parameter> &)441     void getParameters(std::vector<Parameter>&) const{}
make(std::string &,const std::vector<Value> &)442     static iSGLennardJonesForce make(std::string& , const std::vector<Value>&) {
443       return iSGLennardJonesForce();
444     }
445 
446     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
447     // My data members
448     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
449   public:
450     static const std::string keyword;
451 
452   private:
453 
454   };
455   //______________________________________________________________________ INLINES
456 }
457 
458 #endif /* ISGLENNARDJONESFORCE_H */
459