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-2025.
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.realspace;
39  
40  import ffx.crystal.Crystal;
41  import ffx.crystal.CrystalPotential;
42  import ffx.potential.bonded.Atom;
43  import ffx.potential.bonded.LambdaInterface;
44  import ffx.xray.RefinementMinimize.RefinementMode;
45  import ffx.xray.RefinementModel;
46  
47  import java.util.logging.Logger;
48  
49  import static java.lang.String.format;
50  
51  /**
52   * Combine the Real Space target and chemical potential energy.
53   *
54   * @author Timothy D. Fenn
55   * @author Michael J. Schnieders
56   * @since 1.0
57   */
58  public class RealSpaceEnergy implements LambdaInterface, CrystalPotential {
59  
60    private static final Logger logger = Logger.getLogger(RealSpaceEnergy.class.getName());
61  
62    /** The Real Space data to refine against. */
63    private final RealSpaceData realSpaceData;
64    /** The Refinement Model that contains info on mapping between alternate conformers. */
65    private final RefinementModel refinementModel;
66    /** Value of the lambda state variable. */
67    protected double lambda = 1.0;
68    /** The number of parameters that are being refined. */
69    private int nXYZ;
70    /** The refinement mode. */
71    private RefinementMode refinementMode;
72    /** If true, the XYZ coordinates will be refined. */
73    private boolean refineXYZ = false;
74    /** Optimization scaling used to improve convergence. */
75    private double[] optimizationScaling = null;
76    /** Total energy of the refinement. */
77    private double totalEnergy;
78    /** The RealSpaceEnergy is updated with FAST varying energy terms. */
79    private STATE state = STATE.BOTH;
80  
81    /**
82     * Diffraction data energy target
83     *
84     * @param realSpaceData {@link ffx.realspace.RealSpaceData} object to associate with the target
85     * @param nxyz number of xyz parameters
86     * @param nb number of b factor parameters
87     * @param nocc number of occupancy parameters
88     * @param refinementMode the {@link ffx.xray.RefinementMinimize.RefinementMode} type of refinement
89     *     requested
90     */
91    public RealSpaceEnergy(
92        RealSpaceData realSpaceData, int nxyz, int nb, int nocc, RefinementMode refinementMode) {
93      this.realSpaceData = realSpaceData;
94      this.refinementModel = realSpaceData.getRefinementModel();
95      this.refinementMode = refinementMode;
96      this.nXYZ = nxyz;
97      setRefinementBooleans();
98    }
99  
100   /** {@inheritDoc} */
101   @Override
102   public boolean destroy() {
103     return realSpaceData.destroy();
104   }
105 
106   /**
107    * The parameters passed in are only for "active" atoms.
108    *
109    * <p>{@inheritDoc}
110    */
111   @Override
112   public double energy(double[] x) {
113 
114     double e = 0.0;
115 
116     // Unscale the coordinates.
117     unscaleCoordinates(x);
118 
119     if (refineXYZ) {
120       int index = 0;
121       double[] xyz = new double[3];
122       for (Atom a : refinementModel.getTotalAtomArray()) {
123         if (a.isActive()) {
124           int i = index * 3;
125           xyz[0] = x[i];
126           xyz[1] = x[i + 1];
127           xyz[2] = x[i + 2];
128           a.setXYZ(xyz);
129           a.setXYZGradient(0.0, 0.0, 0.0);
130           a.setLambdaXYZGradient(0.0, 0.0, 0.0);
131           index++;
132         }
133       }
134     }
135 
136     // Target function for real space refinement.
137     if (refineXYZ) {
138       e = realSpaceData.computeRealSpaceTarget();
139     }
140 
141     // Scale the coordinates.
142     scaleCoordinates(x);
143 
144     totalEnergy = e;
145     return e;
146   }
147 
148   /** {@inheritDoc} */
149   @Override
150   public double energyAndGradient(double[] x, double[] g) {
151     double e = 0.0;
152     // Unscale the coordinates.
153     unscaleCoordinates(x);
154 
155     if (refineXYZ) {
156       int index = 0;
157       double[] xyz = new double[3];
158       for (Atom a : refinementModel.getTotalAtomArray()) {
159         if (a.isActive()) {
160           int i = index * 3;
161           xyz[0] = x[i];
162           xyz[1] = x[i + 1];
163           xyz[2] = x[i + 2];
164           a.setXYZ(xyz);
165           a.setXYZGradient(0.0, 0.0, 0.0);
166           a.setLambdaXYZGradient(0.0, 0.0, 0.0);
167           index++;
168         }
169       }
170     }
171 
172     // Target function for real space refinement
173     if (refineXYZ) {
174       e = realSpaceData.computeRealSpaceTarget();
175       // Pack gradients into gradient array
176       getXYZGradients(g);
177     }
178 
179     // Scale the coordinates and gradients.
180     scaleCoordinatesAndGradient(x, g);
181 
182     totalEnergy = e;
183     return e;
184   }
185 
186   /** {@inheritDoc} */
187   @Override
188   public double[] getAcceleration(double[] acceleration) {
189     int n = getNumberOfVariables();
190     if (acceleration == null || acceleration.length < n) {
191       acceleration = new double[n];
192     }
193     int index = 0;
194     double[] acc = new double[3];
195     for (Atom a : refinementModel.getTotalAtomArray()) {
196       if (a.isActive()) {
197         a.getAcceleration(acc);
198         acceleration[index++] = acc[0];
199         acceleration[index++] = acc[1];
200         acceleration[index++] = acc[2];
201       }
202     }
203     return acceleration;
204   }
205 
206   /** {@inheritDoc} */
207   @Override
208   public double[] getCoordinates(double[] x) {
209     int n = getNumberOfVariables();
210     if (x == null || x.length < n) {
211       x = new double[n];
212     }
213     int index = 0;
214     for (Atom a : refinementModel.getTotalAtomArray()) {
215       if (a.isActive()) {
216         x[index++] = a.getX();
217         x[index++] = a.getY();
218         x[index++] = a.getZ();
219       }
220     }
221     return x;
222   }
223 
224   /** {@inheritDoc} */
225   @Override
226   public ffx.crystal.Crystal getCrystal() {
227     realSpaceData.getCrystal();
228     return null;
229   }
230 
231   /** {@inheritDoc} */
232   @Override
233   public void setCrystal(Crystal crystal) {
234     logger.severe(" RealSpaceEnergy does implement setCrystal yet.");
235   }
236 
237   /** {@inheritDoc} */
238   @Override
239   public STATE getEnergyTermState() {
240     return state;
241   }
242 
243   /** {@inheritDoc} */
244   @Override
245   public void setEnergyTermState(STATE state) {
246     this.state = state;
247   }
248 
249   /** {@inheritDoc} */
250   @Override
251   public double getLambda() {
252     return lambda;
253   }
254 
255   /** {@inheritDoc} */
256   @Override
257   public void setLambda(double lambda) {
258     if (lambda <= 1.0 && lambda >= 0.0) {
259       this.lambda = lambda;
260       realSpaceData.setLambda(lambda);
261     } else {
262       String message = format(" Lambda value %8.3f is not in the range [0..1].", lambda);
263       logger.warning(message);
264     }
265   }
266 
267   /** {@inheritDoc} */
268   @Override
269   public double[] getMass() {
270     double[] mass = new double[nXYZ];
271     int i = 0;
272     if (refineXYZ) {
273       for (Atom a : refinementModel.getTotalAtomArray()) {
274         if (a.isActive()) {
275           double m = a.getMass();
276           mass[i++] = m;
277           mass[i++] = m;
278           mass[i++] = m;
279         }
280       }
281     }
282     return mass;
283   }
284 
285   /** {@inheritDoc} */
286   @Override
287   public int getNumberOfVariables() {
288     return nXYZ;
289   }
290 
291   /** {@inheritDoc} */
292   @Override
293   public double[] getPreviousAcceleration(double[] previousAcceleration) {
294     int n = getNumberOfVariables();
295     if (previousAcceleration == null || previousAcceleration.length < n) {
296       previousAcceleration = new double[n];
297     }
298     int index = 0;
299     double[] prev = new double[3];
300     for (Atom a : refinementModel.getTotalAtomArray()) {
301       if (a.isActive()) {
302         a.getPreviousAcceleration(prev);
303         previousAcceleration[index++] = prev[0];
304         previousAcceleration[index++] = prev[1];
305         previousAcceleration[index++] = prev[2];
306       }
307     }
308     return previousAcceleration;
309   }
310 
311   /**
312    * Getter for the field <code>refinementMode</code>.
313    *
314    * @return a {@link ffx.xray.RefinementMinimize.RefinementMode} object.
315    */
316   public RefinementMode getRefinementMode() {
317     return refinementMode;
318   }
319 
320   /**
321    * Setter for the field <code>refinementMode</code>.
322    *
323    * @param refinementMode a {@link ffx.xray.RefinementMinimize.RefinementMode} object.
324    */
325   public void setRefinementMode(RefinementMode refinementMode) {
326     this.refinementMode = refinementMode;
327     setRefinementBooleans();
328   }
329 
330   /** {@inheritDoc} */
331   @Override
332   public double[] getScaling() {
333     return optimizationScaling;
334   }
335 
336   /** {@inheritDoc} */
337   @Override
338   public void setScaling(double[] scaling) {
339     optimizationScaling = scaling;
340   }
341 
342   /** {@inheritDoc} */
343   @Override
344   public double getTotalEnergy() {
345     return totalEnergy;
346   }
347 
348   /**
349    * {@inheritDoc}
350    *
351    * <p>Return a reference to each variables type.
352    */
353   @Override
354   public VARIABLE_TYPE[] getVariableTypes() {
355 
356     int nActive = 0;
357     for (Atom a : refinementModel.getTotalAtomArray()) {
358       if (a.isActive()) {
359         nActive++;
360       }
361     }
362 
363     VARIABLE_TYPE[] type = new VARIABLE_TYPE[nActive * 3];
364 
365     int index = 0;
366     for (int i = 0; i < nActive; i++) {
367       type[index++] = VARIABLE_TYPE.X;
368       type[index++] = VARIABLE_TYPE.Y;
369       type[index++] = VARIABLE_TYPE.Z;
370     }
371     return type;
372   }
373 
374   /** {@inheritDoc} */
375   @Override
376   public double[] getVelocity(double[] velocity) {
377     int n = getNumberOfVariables();
378     if (velocity == null || velocity.length < n) {
379       velocity = new double[n];
380     }
381     int index = 0;
382     double[] v = new double[3];
383     for (Atom a : refinementModel.getTotalAtomArray()) {
384       if (a.isActive()) {
385         a.getVelocity(v);
386         velocity[index++] = v[0];
387         velocity[index++] = v[1];
388         velocity[index++] = v[2];
389       }
390     }
391     return velocity;
392   }
393 
394   /** {@inheritDoc} */
395   @Override
396   public double getd2EdL2() {
397     return 0.0;
398   }
399 
400   /** {@inheritDoc} */
401   @Override
402   public double getdEdL() {
403     return realSpaceData.getdEdL();
404   }
405 
406   /** {@inheritDoc} */
407   @Override
408   public void getdEdXdL(double[] gradient) {
409     realSpaceData.getdEdXdL(gradient);
410   }
411 
412   /** {@inheritDoc} */
413   @Override
414   public void setAcceleration(double[] acceleration) {
415     if (acceleration == null) {
416       return;
417     }
418     int index = 0;
419     double[] accel = new double[3];
420     for (Atom a : refinementModel.getTotalAtomArray()) {
421       if (a.isActive()) {
422         accel[0] = acceleration[index++];
423         accel[1] = acceleration[index++];
424         accel[2] = acceleration[index++];
425         a.setAcceleration(accel);
426       }
427     }
428   }
429 
430   /**
431    * Set atomic coordinates positions.
432    *
433    * @param x an array of coordinates for active atoms.
434    */
435   public void setCoordinates(double[] x) {
436     assert (x != null);
437     double[] xyz = new double[3];
438     int index = 0;
439     for (Atom a : refinementModel.getTotalAtomArray()) {
440       if (a.isActive()) {
441         xyz[0] = x[index++];
442         xyz[1] = x[index++];
443         xyz[2] = x[index++];
444         a.moveTo(xyz);
445       }
446     }
447   }
448 
449   /** {@inheritDoc} */
450   @Override
451   public void setPreviousAcceleration(double[] previousAcceleration) {
452     if (previousAcceleration == null) {
453       return;
454     }
455     int index = 0;
456     double[] prev = new double[3];
457     for (Atom a : refinementModel.getTotalAtomArray()) {
458       if (a.isActive()) {
459         prev[0] = previousAcceleration[index++];
460         prev[1] = previousAcceleration[index++];
461         prev[2] = previousAcceleration[index++];
462         a.setPreviousAcceleration(prev);
463       }
464     }
465   }
466 
467   /** {@inheritDoc} */
468   @Override
469   public void setVelocity(double[] velocity) {
470     if (velocity == null) {
471       return;
472     }
473     int index = 0;
474     double[] vel = new double[3];
475     for (Atom a : refinementModel.getTotalAtomArray()) {
476       if (a.isActive()) {
477         vel[0] = velocity[index++];
478         vel[1] = velocity[index++];
479         vel[2] = velocity[index++];
480         a.setVelocity(vel);
481       }
482     }
483   }
484 
485   /**
486    * If the refinement mode has changed, this should be called to update which parameters are being
487    * fit.
488    */
489   private void setRefinementBooleans() {
490     refineXYZ =
491         (refinementMode == RefinementMode.COORDINATES
492             || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS
493             || refinementMode == RefinementMode.COORDINATES_AND_OCCUPANCIES
494             || refinementMode == RefinementMode.COORDINATES_AND_BFACTORS_AND_OCCUPANCIES);
495   }
496 
497   /**
498    * Get the number of xyz parameters being fit.
499    *
500    * @return the number of xyz parameters
501    */
502   private int getNXYZ() {
503     return nXYZ;
504   }
505 
506   /**
507    * Set the number of xyz parameters.
508    *
509    * @param nxyz requested number of xyz parameters
510    */
511   private void setNXYZ(int nxyz) {
512     this.nXYZ = nxyz;
513   }
514 
515   /**
516    * Fill gradient array with xyz gradient.
517    *
518    * @param g array to add gradient to
519    */
520   private void getXYZGradients(double[] g) {
521     assert (g != null);
522     double[] grad = new double[3];
523     int index = 0;
524     for (Atom a : refinementModel.getTotalAtomArray()) {
525       if (a.isActive()) {
526         a.getXYZGradient(grad);
527         g[index++] = grad[0];
528         g[index++] = grad[1];
529         g[index++] = grad[2];
530       }
531     }
532   }
533 }