View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2021.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.potential.nonbonded;
39  
40  import static ffx.numerics.math.DoubleMath.length2;
41  import static java.lang.String.format;
42  import static java.lang.System.arraycopy;
43  import static java.util.Arrays.fill;
44  import static org.apache.commons.math3.util.FastMath.pow;
45  
46  import ffx.potential.bonded.Atom;
47  import ffx.potential.bonded.LambdaInterface;
48  import ffx.potential.parameters.AtomType;
49  import ffx.potential.parameters.ForceField;
50  import java.util.Arrays;
51  import java.util.Comparator;
52  import java.util.logging.Logger;
53  import org.apache.commons.configuration2.CompositeConfiguration;
54  
55  /**
56   * Restrain atoms to their initial coordinates.
57   *
58   * @author Michael J. Schnieders
59   * @since 1.0
60   */
61  public class CoordRestraint implements LambdaInterface {
62  
63    private static final Logger logger = Logger.getLogger(CoordRestraint.class.getName());
64  
65    private final Atom[] atoms;
66    private final double[][] initialCoordinates;
67    /** Force constant variable stores K/2 in Kcal/mol/A. E = K/2 * dx^2. */
68    private final double forceConstant;
69  
70    private final double[] a1 = new double[3];
71    private final double[] dx = new double[3];
72    private final double lambdaExp = 1.0;
73    private final double[] lambdaGradient;
74    private final String atomIndexRestraints;
75    private final int atom1Index;
76    private final int atom2Index;
77    private final int atom3Index;
78    private final String atomTypeRestraints;
79    private final AtomType atom1Type;
80    private final AtomType atom2Type;
81    private final AtomType atom3Type;
82    private int nAtoms = 0;
83    private double lambda = 1.0;
84    private double lambdaPow = pow(lambda, lambdaExp);
85    private double dLambdaPow = lambdaExp * pow(lambda, lambdaExp - 1.0);
86    private double d2LambdaPow = 0;
87    private double dEdL = 0.0;
88    private double d2EdL2 = 0.0;
89    private boolean lambdaTerm;
90    private boolean ignoreHydrogen = true;
91  
92    /**
93     * This CoordRestraint is based on the unit cell parameters and symmetry operators of the supplied
94     * crystal.
95     *
96     * @param atoms the Atom array to base this CoordRestraint on.
97     * @param forceField the ForceField to apply.
98     */
99    public CoordRestraint(Atom[] atoms, ForceField forceField) {
100     this(atoms, forceField, forceField.getBoolean("RESTRAIN_WITH_LAMBDA", true));
101   }
102 
103   /**
104    * This CoordRestraint is based on the unit cell parameters and symmetry operators of the supplied
105    * crystal.
106    *
107    * @param atoms the Atom array to base this CoordRestraint on.
108    * @param forceField the ForceField to apply.
109    * @param useLambda If false, do not apply lambda term to this restraint.
110    */
111   public CoordRestraint(Atom[] atoms, ForceField forceField, boolean useLambda) {
112     this(atoms, forceField, useLambda, forceField.getDouble("RESTRAINT_K", 10.0));
113   }
114 
115   /**
116    * This CoordRestraint is based on the unit cell parameters and symmetry operators of the supplied
117    * crystal.
118    *
119    * @param atoms the Atom array to base this CoordRestraint on.
120    * @param forceField the ForceField to apply.
121    * @param useLambda If false, do not apply lambda term to this restraint.
122    * @param forceConst Force constant to apply
123    */
124   public CoordRestraint(Atom[] atoms, ForceField forceField, boolean useLambda, double forceConst) {
125     this.atoms = atoms;
126     nAtoms = atoms.length;
127 
128     if (useLambda) {
129       lambdaTerm = forceField.getBoolean("LAMBDATERM", false);
130     } else {
131       lambdaTerm = false;
132     }
133 
134     if (lambdaTerm) {
135       lambdaGradient = new double[nAtoms * 3];
136     } else {
137       lambdaGradient = null;
138       this.lambda = 1.0;
139       lambdaPow = 1.0;
140       dLambdaPow = 0.0;
141       d2LambdaPow = 0.0;
142     }
143 
144     forceConstant = forceConst;
145 
146     logger.info(
147         format(
148             "\n Coordinate Restraint Atoms (k = %6.3f, lambdaTerm=%s):",
149             forceConstant * 2.0, lambdaTerm));
150 
151     initialCoordinates = new double[3][nAtoms];
152     for (int i = 0; i < nAtoms; i++) {
153       Atom a = atoms[i];
154       initialCoordinates[0][i] = a.getX();
155       initialCoordinates[1][i] = a.getY();
156       initialCoordinates[2][i] = a.getZ();
157     }
158 
159     CompositeConfiguration properties = forceField.getProperties();
160 
161     // pdbAtomRestraints uses atom indexes to restrain specific atoms
162     atomIndexRestraints = properties.getString("pdbAtomRestraints");
163     if (atomIndexRestraints != null) {
164       String[] tokens = atomIndexRestraints.split(",");
165       // Pick up atom index for reference when looking at multiple molecules.
166       atom1Index = Integer.parseInt(tokens[0]);
167       atom2Index = Integer.parseInt(tokens[1]);
168       atom3Index = Integer.parseInt(tokens[2]);
169     } else {
170       atom1Index = -1;
171       atom2Index = -1;
172       atom3Index = -1;
173     }
174 
175     /*
176      xyzAtomRestraints uses atom types to restrain specific atoms. This can result in more
177      atoms being restrained than desired since atom types are not unique to each atom.
178     */
179     atomTypeRestraints = properties.getString("xyzAtomRestraints");
180     if (atomTypeRestraints != null) {
181       String[] tokens = atomTypeRestraints.split(",");
182       // Pick up atom type for reference when looking at multiple molecules
183       atom1Type = atoms[Integer.parseInt(tokens[0])].getAtomType();
184       atom2Type = atoms[Integer.parseInt(tokens[1])].getAtomType();
185       atom3Type = atoms[Integer.parseInt(tokens[2])].getAtomType();
186     } else {
187       atom1Type = null;
188       atom2Type = null;
189       atom3Type = null;
190     }
191 
192     Arrays.stream(atoms).sorted(Comparator.comparingInt(Atom::getIndex)).forEach(Atom::print);
193   }
194 
195   /**
196    * Returns a copy of the atoms array.
197    *
198    * @return Copy of the atom array.
199    */
200   public Atom[] getAtoms() {
201     Atom[] retArray = new Atom[nAtoms];
202     arraycopy(atoms, 0, retArray, 0, nAtoms);
203     return retArray;
204   }
205 
206   /**
207    * getCoordinatePin.
208    *
209    * @param xyz an array of {@link double} objects.
210    * @return an array of {@link double} objects.
211    */
212   public double[][] getCoordinatePin(double[][] xyz) {
213     if (xyz == null) {
214       xyz = new double[nAtoms][3];
215     } else if (xyz.length != nAtoms) {
216       throw new IllegalArgumentException(" Incorrect number of atoms!");
217     }
218     for (int i = 0; i < nAtoms; i++) {
219       arraycopy(initialCoordinates[i], 0, xyz[i], 0, 3);
220     }
221     return xyz;
222   }
223 
224   /**
225    * Returns the force constant in kcal/mol/Angstrom^2.
226    *
227    * @return a double.
228    */
229   public double getForceConstant() {
230     return forceConstant;
231   }
232 
233   /** {@inheritDoc} */
234   @Override
235   public double getLambda() {
236     return lambda;
237   }
238 
239   /** {@inheritDoc} */
240   @Override
241   public void setLambda(double lambda) {
242     if (lambdaTerm) {
243       this.lambda = lambda;
244 
245       double lambdaWindow = 1.0;
246       if (this.lambda <= lambdaWindow) {
247         double dldgl = 1.0 / lambdaWindow;
248         double l = dldgl * this.lambda;
249         double l2 = l * l;
250         double l3 = l2 * l;
251         double l4 = l2 * l2;
252         double l5 = l4 * l;
253         double c3 = 10.0;
254         double c4 = -15.0;
255         double c5 = 6.0;
256         double threeC3 = 3.0 * c3;
257         double sixC3 = 6.0 * c3;
258         double fourC4 = 4.0 * c4;
259         double twelveC4 = 12.0 * c4;
260         double fiveC5 = 5.0 * c5;
261         double twentyC5 = 20.0 * c5;
262         lambdaPow = c3 * l3 + c4 * l4 + c5 * l5;
263         dLambdaPow = (threeC3 * l2 + fourC4 * l3 + fiveC5 * l4) * dldgl;
264         d2LambdaPow = (sixC3 * l + twelveC4 * l2 + twentyC5 * l3) * dldgl * dldgl;
265       } else {
266         lambdaPow = 1.0;
267         dLambdaPow = 0.0;
268         d2LambdaPow = 0.0;
269       }
270     } else {
271       this.lambda = 1.0;
272       lambdaPow = 1.0;
273       dLambdaPow = 0.0;
274       d2LambdaPow = 0.0;
275     }
276   }
277 
278   /**
279    * getNumAtoms.
280    *
281    * @return a int.
282    */
283   public int getNumAtoms() {
284     return nAtoms;
285   }
286 
287   /**
288    * Returns the original coordinates of this restraint, indexed by atoms then x,y,z. This is the
289    * opposite order of the internal storage.
290    *
291    * @return Original coordinates [atoms][xyz]
292    */
293   public double[][] getOriginalCoordinates() {
294     double[][] retArray = new double[nAtoms][3];
295     for (int i = 0; i < nAtoms; i++) {
296       for (int j = 0; j < 3; j++) {
297         // Mild subtlety here: it is stored internally as [xyz][atoms], but returned as [atoms][xyz]
298         retArray[i][j] = initialCoordinates[j][i];
299       }
300     }
301     return retArray;
302   }
303 
304   /** {@inheritDoc} */
305   @Override
306   public double getd2EdL2() {
307     if (lambdaTerm) {
308       return d2EdL2;
309     } else {
310       return 0.0;
311     }
312   }
313 
314   /** {@inheritDoc} */
315   @Override
316   public double getdEdL() {
317     if (lambdaTerm) {
318       return dEdL;
319     } else {
320       return 0.0;
321     }
322   }
323 
324   /** {@inheritDoc} */
325   @Override
326   public void getdEdXdL(double[] gradient) {
327     if (lambdaTerm) {
328       int n3 = nAtoms * 3;
329       for (int i = 0; i < n3; i++) {
330         gradient[i] += lambdaGradient[i];
331       }
332     }
333   }
334 
335   /** resetCoordinatePin. */
336   public void resetCoordinatePin() {
337     double[] aixyz = new double[3];
338     for (int i = 0; i < nAtoms; i++) {
339       atoms[i].getXYZ(aixyz);
340       arraycopy(aixyz, 0, initialCoordinates[i], 0, 3);
341     }
342   }
343 
344   /**
345    * Calculates energy and gradients for this coordinate restraint.
346    *
347    * @param gradient Calculate gradients
348    * @param print Unused
349    * @return Energy in the coordinate restraint
350    */
351   public double residual(boolean gradient, boolean print) {
352 
353     if (lambdaTerm) {
354       dEdL = 0.0;
355       d2EdL2 = 0.0;
356       fill(lambdaGradient, 0.0);
357     }
358 
359     int atomIndex = 0;
360 
361     // Assuming that the first molecule in pdb is labeled as 1.
362     int moleculeNumber = 1;
363     double residual = 0.0;
364     double fx2 = forceConstant * 2.0;
365     for (int i = 0; i < nAtoms; i++) {
366       // Current atomic coordinates.
367       Atom atom = atoms[i];
368       if (atomTypeRestraints != null) {
369         if (atom.getAtomType() == atom1Type
370             || atom.getAtomType() == atom2Type
371             || atom.getAtomType() == atom3Type) {
372           atom.getXYZ(a1);
373           // Compute their separation.
374           dx[0] = a1[0] - initialCoordinates[0][i];
375           dx[1] = a1[1] - initialCoordinates[1][i];
376           dx[2] = a1[2] - initialCoordinates[2][i];
377         } else {
378           continue;
379         }
380         double r2 = length2(dx);
381         residual += r2;
382         if (gradient || lambdaTerm) {
383           final double dedx = dx[0] * fx2;
384           final double dedy = dx[1] * fx2;
385           final double dedz = dx[2] * fx2;
386           if (gradient) {
387             atom.addToXYZGradient(lambdaPow * dedx, lambdaPow * dedy, lambdaPow * dedz);
388           }
389           if (lambdaTerm) {
390             int j3 = i * 3;
391             lambdaGradient[j3] = dLambdaPow * dedx;
392             lambdaGradient[j3 + 1] = dLambdaPow * dedy;
393             lambdaGradient[j3 + 2] = dLambdaPow * dedz;
394           }
395         }
396       } else if (atomIndexRestraints != null) {
397         atomIndex += 1;
398         // Following if statement works only for .pdb files since they give residueNumber
399         if (atom.getResidueNumber() != moleculeNumber) {
400           atomIndex = 1;
401           moleculeNumber = atom.getResidueNumber();
402         }
403         if (atomIndex == atom1Index || atomIndex == atom2Index || atomIndex == atom3Index) {
404           atom.getXYZ(a1);
405           // Compute their separation.
406           dx[0] = a1[0] - initialCoordinates[0][i];
407           dx[1] = a1[1] - initialCoordinates[1][i];
408           dx[2] = a1[2] - initialCoordinates[2][i];
409         } else {
410           continue;
411         }
412         double r2 = length2(dx);
413         residual += r2;
414         if (gradient || lambdaTerm) {
415           final double dedx = dx[0] * fx2;
416           final double dedy = dx[1] * fx2;
417           final double dedz = dx[2] * fx2;
418           if (gradient) {
419             atom.addToXYZGradient(lambdaPow * dedx, lambdaPow * dedy, lambdaPow * dedz);
420           }
421           if (lambdaTerm) {
422             int j3 = i * 3;
423             lambdaGradient[j3] = dLambdaPow * dedx;
424             lambdaGradient[j3 + 1] = dLambdaPow * dedy;
425             lambdaGradient[j3 + 2] = dLambdaPow * dedz;
426           }
427         }
428       } else {
429 
430         if (ignoreHydrogen && atom.isHydrogen()) {
431           continue;
432         }
433 
434         atom.getXYZ(a1);
435         // Compute their separation.
436         dx[0] = a1[0] - initialCoordinates[0][i];
437         dx[1] = a1[1] - initialCoordinates[1][i];
438         dx[2] = a1[2] - initialCoordinates[2][i];
439         double r2 = length2(dx);
440         residual += r2;
441         if (gradient || lambdaTerm) {
442           final double dedx = dx[0] * fx2;
443           final double dedy = dx[1] * fx2;
444           final double dedz = dx[2] * fx2;
445           if (gradient) {
446             atom.addToXYZGradient(lambdaPow * dedx, lambdaPow * dedy, lambdaPow * dedz);
447           }
448           if (lambdaTerm) {
449             int j3 = i * 3;
450             lambdaGradient[j3] = dLambdaPow * dedx;
451             lambdaGradient[j3 + 1] = dLambdaPow * dedy;
452             lambdaGradient[j3 + 2] = dLambdaPow * dedz;
453           }
454         }
455       }
456     }
457     if (lambdaTerm) {
458       dEdL = dLambdaPow * forceConstant * residual;
459       d2EdL2 = d2LambdaPow * forceConstant * residual;
460     }
461 
462     return forceConstant * residual * lambdaPow;
463   }
464 
465   /**
466    * Sets the initial coordinates to new values.
467    *
468    * @param newInitialCoordinates Indexed xyz,atom number.
469    */
470   public void setCoordinatePin(double[][] newInitialCoordinates) {
471     if (newInitialCoordinates.length != initialCoordinates.length) {
472       throw new IllegalArgumentException(" Incorrect number of atoms!");
473     }
474     for (int i = 0; i < 3; i++) {
475       arraycopy(
476           newInitialCoordinates[i], 0, initialCoordinates[i], 0, initialCoordinates[0].length);
477     }
478   }
479 
480   /**
481    * Setter for the field <code>ignoreHydrogen</code>.
482    *
483    * @param ignoreHydrogen a boolean.
484    */
485   public void setIgnoreHydrogen(boolean ignoreHydrogen) {
486     this.ignoreHydrogen = ignoreHydrogen;
487   }
488 }