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-2024.
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.algorithms.optimize;
39  
40  import static java.lang.String.format;
41  import static java.lang.System.arraycopy;
42  import static org.apache.commons.math3.util.FastMath.sqrt;
43  
44  import ffx.algorithms.AlgorithmListener;
45  import ffx.algorithms.Terminatable;
46  import ffx.crystal.Crystal;
47  import ffx.numerics.Potential;
48  import ffx.numerics.optimization.OptimizationListener;
49  import ffx.potential.ForceFieldEnergy;
50  import ffx.potential.MolecularAssembly;
51  import ffx.potential.MolecularAssembly.FractionalMode;
52  import ffx.potential.XtalEnergy;
53  import ffx.utilities.Constants;
54  import java.util.logging.Logger;
55  
56  /**
57   * Minimize the energy of a system to an RMS gradient per atom convergence criteria.
58   *
59   * @author Michael J. Schnieders
60   * @since 1.0
61   */
62  public class CrystalMinimize extends Minimize implements OptimizationListener, Terminatable {
63  
64    private static final Logger logger = Logger.getLogger(CrystalMinimize.class.getName());
65  
66    private final ForceFieldEnergy forceFieldEnergy;
67    private Crystal crystal;
68    private final Crystal unitCell;
69  
70    /**
71     * Constructor for CrystalMinimize.
72     *
73     * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
74     * @param xtalEnergy a {@link ffx.potential.XtalEnergy} object.
75     * @param algorithmListener a {@link ffx.algorithms.AlgorithmListener} object.
76     */
77    public CrystalMinimize(MolecularAssembly molecularAssembly, XtalEnergy xtalEnergy,
78        AlgorithmListener algorithmListener) {
79      super(molecularAssembly, xtalEnergy, algorithmListener);
80      crystal = molecularAssembly.getCrystal();
81      unitCell = crystal.getUnitCell();
82      forceFieldEnergy = molecularAssembly.getPotentialEnergy();
83  
84      for (int i = 0; i < n - 6; i += 3) {
85        scaling[i] = 12.0 * crystal.a;
86        scaling[i + 1] = 12.0 * crystal.b;
87        scaling[i + 2] = 12.0 * crystal.c;
88      }
89      scaling[n - 6] = 4.0 * sqrt(crystal.a);
90      scaling[n - 5] = 4.0 * sqrt(crystal.b);
91      scaling[n - 4] = 4.0 * sqrt(crystal.c);
92      scaling[n - 3] = 0.02 * sqrt(crystal.volume);
93      scaling[n - 2] = 0.02 * sqrt(crystal.volume);
94      scaling[n - 1] = 0.02 * sqrt(crystal.volume);
95    }
96  
97    /**
98     * computeStressTensor.
99     *
100    * @param verbose a boolean.
101    */
102   public void computeElasticityTensor(boolean verbose) {
103 
104     double[] xOrig = new double[forceFieldEnergy.getNumberOfVariables()];
105     forceFieldEnergy.getCoordinates(xOrig);
106 
107     forceFieldEnergy.energy(xOrig, verbose);
108 
109     FractionalMode currentFractionalMode = molecularAssembly.getFractionalMode();
110 
111     molecularAssembly.setFractionalMode(FractionalMode.MOLECULE);
112     molecularAssembly.computeFractionalCoordinates();
113 
114     // Conversion from kcal/mol/Ang^3 to Atm.
115     double PRESCON = 6.85684112e4;
116 
117     // Conversion from Atm to GPa.
118     double GPA = 101325 * 1.0e-9;
119 
120     double volume =
121         crystal.getUnitCell().volume / crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
122 
123     boolean currentCheckRestrictions = crystal.getCheckRestrictions();
124 
125     crystal.setCheckRestrictions(false);
126 
127     logger.info(" The Strain and Stress Tensor code is under development.");
128 
129     double delta = 5.0e-3;
130     double eps = 2.0e-3;
131 
132     double c11 = 0.0, c22 = 0.0, c33 = 0.0;
133     double c12 = 0.0, c13 = 0.0, c23 = 0.0;
134     double c14 = 0.0, c15 = 0.0, c16 = 0.0;
135     double c24 = 0.0, c25 = 0.0, c26 = 0.0;
136     double c34 = 0.0, c35 = 0.0, c36 = 0.0;
137     double c44 = 0.0, c55 = 0.0, c66 = 0.0;
138     double c45 = 0.0, c46 = 0.0, c56 = 0.0;
139 
140     // 1 -> XX
141     // 2 -> YY
142     // 3 -> ZZ
143     // 4 -> YZ
144     // 5 -> XZ
145     // 6 -> XY
146 
147     double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
148     arraycopy(xOrig, 0, x, 0, x.length);
149 
150     switch (crystal.spaceGroup.crystalSystem) {
151       case CUBIC:
152         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
153         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
154         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
155         double cavg = (c11 + c22 + c33) / 3.0;
156         c11 = c22 = c33 = cavg;
157 
158         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
159         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
160         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
161         cavg = (c12 + c13 + c23) / 3.0;
162         c12 = c13 = c23 = cavg;
163 
164         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
165         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
166         c66 = dE2dA2(6, 6, delta, x, eps) * PRESCON * GPA / volume;
167         cavg = (c44 + c55 + c66) / 3.0;
168         c44 = c55 = c66 = cavg;
169         break;
170       case HEXAGONAL:
171         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
172         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
173         cavg = (c11 + c22) / 2.0;
174         c11 = c22 = cavg;
175         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
176 
177         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
178         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
179         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
180         cavg = (c13 + c23) / 2.0;
181         c13 = c23 = cavg;
182 
183         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
184         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
185         cavg = (c44 + c55) / 2.0;
186         c44 = c55 = cavg;
187         c66 = 0.5 * (c11 - c12);
188         break;
189       case TRIGONAL:
190         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
191         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
192         cavg = (c11 + c22) / 2.0;
193         c11 = c22 = cavg;
194         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
195 
196         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
197         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
198         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
199         cavg = (c13 + c23) / 2.0;
200         c13 = c23 = cavg;
201 
202         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
203         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
204         cavg = (c44 + c55) / 2.0;
205         c44 = c55 = cavg;
206         c66 = 0.5 * (c11 - c12);
207 
208         // c14 and c24 should be equal and opposite sign.
209         c14 = dE2dA2(1, 4, delta, x, eps) * PRESCON * GPA / volume;
210         c24 = dE2dA2(2, 4, delta, x, eps) * PRESCON * GPA / volume;
211         c56 = c14;
212 
213         String pg = crystal.getUnitCell().spaceGroup.pointGroupName;
214         if (pg.equalsIgnoreCase("PG3") || pg.equalsIgnoreCase("PG3bar")) {
215           // Should be equal in magnitude and opposite sign.
216           c15 = dE2dA2(1, 5, delta, x, eps) * PRESCON * GPA / volume;
217           c25 = dE2dA2(2, 5, delta, x, eps) * PRESCON * GPA / volume;
218           c46 = c25;
219         }
220         break;
221       case TETRAGONAL:
222         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
223         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
224         cavg = (c11 + c22) / 2.0;
225         c11 = c22 = cavg;
226         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
227 
228         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
229         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
230         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
231         cavg = (c13 + c23) / 2.0;
232         c13 = c23 = cavg;
233 
234         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
235         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
236         cavg = (c44 + c55) / 2.0;
237         c44 = c55 = cavg;
238         c66 = dE2dA2(6, 6, delta, x, eps) * PRESCON * GPA / volume;
239         pg = crystal.getUnitCell().spaceGroup.pointGroupName;
240         if (pg.equalsIgnoreCase("PG4") || pg.equalsIgnoreCase("PG4bar") || pg.equalsIgnoreCase(
241             "PG4/m")) {
242           // Should be equal in magnitude and opposite sign.
243           c16 = dE2dA2(1, 6, delta, x, eps) * PRESCON * GPA / volume;
244           c26 = dE2dA2(2, 6, delta, x, eps) * PRESCON * GPA / volume;
245         }
246         break;
247       case ORTHORHOMBIC:
248         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
249         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
250         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
251 
252         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
253         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
254         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
255 
256         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
257         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
258         c66 = dE2dA2(6, 6, delta, x, eps) * PRESCON * GPA / volume;
259         break;
260       case MONOCLINIC:
261         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
262         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
263         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
264 
265         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
266         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
267         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
268 
269         c15 = dE2dA2(1, 5, delta, x, eps) * PRESCON * GPA / volume;
270         c25 = dE2dA2(2, 5, delta, x, eps) * PRESCON * GPA / volume;
271         c35 = dE2dA2(3, 5, delta, x, eps) * PRESCON * GPA / volume;
272         c46 = dE2dA2(4, 6, delta, x, eps) * PRESCON * GPA / volume;
273 
274         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
275         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
276         c66 = dE2dA2(6, 6, delta, x, eps) * PRESCON * GPA / volume;
277         break;
278       case TRICLINIC:
279       default:
280         c11 = dE2dA2(1, 1, delta, x, eps) * PRESCON * GPA / volume;
281         c22 = dE2dA2(2, 2, delta, x, eps) * PRESCON * GPA / volume;
282         c33 = dE2dA2(3, 3, delta, x, eps) * PRESCON * GPA / volume;
283 
284         c12 = dE2dA2(1, 2, delta, x, eps) * PRESCON * GPA / volume;
285         c13 = dE2dA2(1, 3, delta, x, eps) * PRESCON * GPA / volume;
286         c23 = dE2dA2(2, 3, delta, x, eps) * PRESCON * GPA / volume;
287 
288         c14 = dE2dA2(1, 4, delta, x, eps) * PRESCON * GPA / volume;
289         c15 = dE2dA2(1, 5, delta, x, eps) * PRESCON * GPA / volume;
290         c16 = dE2dA2(1, 6, delta, x, eps) * PRESCON * GPA / volume;
291 
292         c24 = dE2dA2(2, 4, delta, x, eps) * PRESCON * GPA / volume;
293         c25 = dE2dA2(2, 5, delta, x, eps) * PRESCON * GPA / volume;
294         c26 = dE2dA2(2, 6, delta, x, eps) * PRESCON * GPA / volume;
295 
296         c34 = dE2dA2(3, 4, delta, x, eps) * PRESCON * GPA / volume;
297         c35 = dE2dA2(3, 5, delta, x, eps) * PRESCON * GPA / volume;
298         c36 = dE2dA2(3, 6, delta, x, eps) * PRESCON * GPA / volume;
299 
300         c44 = dE2dA2(4, 4, delta, x, eps) * PRESCON * GPA / volume;
301         c55 = dE2dA2(5, 5, delta, x, eps) * PRESCON * GPA / volume;
302         c66 = dE2dA2(6, 6, delta, x, eps) * PRESCON * GPA / volume;
303 
304         c45 = dE2dA2(4, 5, delta, x, eps) * PRESCON * GPA / volume;
305         c46 = dE2dA2(4, 6, delta, x, eps) * PRESCON * GPA / volume;
306         c56 = dE2dA2(5, 6, delta, x, eps) * PRESCON * GPA / volume;
307     }
308 
309     if (eps > 0) {
310       logger.info(
311           format(" Elasticity Tensor using minimization and FD step size of %6.3e (GPa) =", delta));
312     } else {
313       logger.info(format(
314           " Elasticity Tensor using rigid body fractional coordinates and FD step size of %6.3e (GPa) =",
315           delta));
316     }
317     logger.info(
318         format(" [ %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f ]", c11, c12, c13, c14, c15, c16));
319     logger.info(format(" [ %12s %12.3f %12.3f %12.3f %12.3f %12.3f ]", "", c22, c23, c24, c25, c26));
320     logger.info(format(" [ %12s %12s %12.3f %12.3f %12.3f %12.3f ]", "", "", c33, c34, c35, c36));
321     logger.info(format(" [ %12s %12s %12s %12.3f %12.3f %12.3f ]", "", "", "", c44, c45, c46));
322     logger.info(format(" [ %12s %12s %12s %12s %12.3f %12.3f ]", "", "", "", "", c55, c56));
323     logger.info(format(" [ %12s %12s %12s %12s %12s %12.3f ]", "", "", "", "", "", c66));
324 
325     // forceFieldEnergy.energy(xOrig, verbose);
326     molecularAssembly.setFractionalMode(currentFractionalMode);
327     crystal.setCheckRestrictions(currentCheckRestrictions);
328   }
329 
330   /**
331    * minimize
332    *
333    * @param m The number of previous steps used to estimate the Hessian.
334    * @param maxIterations The maximum number of minimization steps.
335    * @param eps The minimization convergence criteria.
336    * @return a {@link ffx.numerics.Potential} object.
337    */
338   public Potential minimize(int m, double eps, int maxIterations) {
339     super.minimize(m, eps, maxIterations);
340 
341     crystal = molecularAssembly.getCrystal();
342     logger.info("\n Final lattice parameters" + crystal);
343 
344     return potential;
345   }
346 
347   /**
348    * Print out the partial derivatives of the Energy with respect to components of the 3 vectors that
349    * define the primitive cell.
350    */
351   public void printTensor() {
352     computeStressTensor(true);
353   }
354 
355   /**
356    * Apply a strain delta.
357    *
358    * @param voight Voight index
359    * @param delta Strain delta
360    * @param strain Strain matrix
361    */
362   private void applyStrain(int voight, double delta, double[][] strain) {
363     switch (voight) {
364       case 1 -> // XX
365           strain[0][0] += delta;
366 
367       // strain[1][1] -= 1.0 * delta / 3.0;
368       // strain[2][2] -= 1.0 * delta / 3.0;
369       case 2 -> // YY
370         // strain[0][0] -= 1.0 * delta / 3.0;
371           strain[1][1] += delta;
372 
373       // strain[2][2] -= 1.0 * delta / 3.0;
374       case 3 -> // ZZ
375         // strain[0][0] -= 1.0 * delta / 3.0;
376         // strain[1][1] -= 1.0 * delta / 3.0;
377           strain[2][2] += delta;
378       case 4 -> { // YZ
379         strain[1][2] += delta / 2.0;
380         strain[2][1] += delta / 2.0;
381       }
382       case 5 -> { // XZ
383         strain[0][2] += delta / 2.0;
384         strain[2][0] += delta / 2.0;
385       }
386       case 6 -> { // XY
387         strain[0][1] += delta / 2.0;
388         strain[1][0] += delta / 2.0;
389       }
390     }
391   }
392 
393   private void computeStressTensor(boolean verbose) {
394 
395     double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
396     forceFieldEnergy.getCoordinates(x);
397 
398     forceFieldEnergy.energy(x, verbose);
399 
400     FractionalMode currentFractionalMode = molecularAssembly.getFractionalMode();
401 
402     molecularAssembly.setFractionalMode(FractionalMode.ATOM);
403     molecularAssembly.computeFractionalCoordinates();
404 
405     double delta = 1.0e-5;
406     double[][] dEdV = new double[3][3];
407 
408     dEdV[0][0] = dEdA(0, 0, delta, x);
409     dEdV[0][1] = dEdA(0, 1, delta, x);
410     dEdV[0][2] = dEdA(0, 2, delta, x);
411     dEdV[1][1] = dEdA(1, 1, delta, x);
412     dEdV[1][2] = dEdA(1, 2, delta, x);
413     dEdV[2][2] = dEdA(2, 2, delta, x);
414 
415     logger.info("\n The Stress Tensor code is under development.");
416 
417     logger.info("\n dE/dLvec Derivatives: ");
418     logger.info(format(" [ %16.8f %16.8f %16.8f ]", dEdV[0][0], dEdV[0][1], dEdV[0][2]));
419     logger.info(format(" [ %16.8f %16.8f %16.8f ]", dEdV[1][0], dEdV[1][1], dEdV[1][2]));
420     logger.info(format(" [ %16.8f %16.8f %16.8f ]", dEdV[2][0], dEdV[2][1], dEdV[2][2]));
421 
422     dEdV[1][0] = dEdV[0][1];
423     dEdV[2][0] = dEdV[0][2];
424     dEdV[2][1] = dEdV[1][2];
425 
426     // Compute a numerical virial.
427     double[][] virial = new double[3][3];
428     for (int i = 0; i < 3; i++) {
429       for (int j = 0; j <= i; j++) {
430         virial[j][i] = 0.0;
431         for (int k = 0; k < 3; k++) {
432           virial[j][i] = virial[j][i] + dEdV[k][j] * unitCell.Ai[i][k];
433         }
434         virial[i][j] = virial[j][i];
435       }
436     }
437 
438     logger.info("\n Numerical Virial: ");
439     logger.info(format(" [ %16.8f %16.8f %16.8f ]", virial[0][0], virial[0][1], virial[0][2]));
440     logger.info(format(" [ %16.8f %16.8f %16.8f ]", virial[1][0], virial[1][1], virial[1][2]));
441     logger.info(format(" [ %16.8f %16.8f %16.8f ]", virial[2][0], virial[2][1], virial[2][2]));
442 
443     double dedv = (virial[0][0] + virial[1][1] + virial[2][2]) / (3 * unitCell.volume);
444 
445     // Conversion from kcal/mol/Ang^3 to Atm.
446     double PRESCON = 6.85684112e4;
447     double temperature = 298.15;
448     int nAtoms = molecularAssembly.getAtomArray().length;
449     double pressure = PRESCON * (nAtoms * Constants.R * temperature / unitCell.volume - dedv);
450 
451     logger.info(format("\n Pressure estimate (%6.2f Kelvin): %12.8f (atm)", temperature, pressure));
452 
453     molecularAssembly.setFractionalMode(currentFractionalMode);
454   }
455 
456   private double dEdA(int ii, int jj, double delta, double[] x) {
457 
458     // Store current unit cell parameters.
459     double a = unitCell.a;
460     double b = unitCell.b;
461     double c = unitCell.c;
462     double alpha = unitCell.alpha;
463     double beta = unitCell.beta;
464     double gamma = unitCell.gamma;
465 
466     double[][] cellVectors = new double[3][3];
467     for (int i = 0; i < 3; i++) {
468       System.arraycopy(unitCell.Ai[i], 0, cellVectors[i], 0, 3);
469     }
470 
471     cellVectors[ii][jj] += delta;
472     try {
473       if (!crystal.setCellVectors(cellVectors)) {
474         return 0.0;
475       }
476     } catch (Exception e) {
477       return 0.0;
478     }
479 
480     // Update the ForceFieldEnergy with the new boundary conditions.
481     forceFieldEnergy.setCrystal(crystal);
482 
483     // Update and retrieve coordinates.
484     molecularAssembly.moveToFractionalCoordinates();
485     forceFieldEnergy.getCoordinates(x);
486 
487     // Compute the energy.
488     double eplus = forceFieldEnergy.energy(x);
489 
490     cellVectors[ii][jj] -= 2 * delta;
491     try {
492       if (!crystal.setCellVectors(cellVectors)) {
493         return 0.0;
494       }
495     } catch (Exception e) {
496       return 0.0;
497     }
498 
499     // Update the ForceFieldEnergy with the new boundary conditions.
500     forceFieldEnergy.setCrystal(crystal);
501 
502     // Update and retrieve coordinates.
503     molecularAssembly.moveToFractionalCoordinates();
504     forceFieldEnergy.getCoordinates(x);
505 
506     // Compute the energy.
507     double eminus = forceFieldEnergy.energy(x);
508 
509     // Revert to the original boundary conditions.
510     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
511     forceFieldEnergy.setCrystal(crystal);
512     molecularAssembly.moveToFractionalCoordinates();
513 
514     return (eplus - eminus) / (2 * delta);
515   }
516 
517   private double dE2dA2(int voight1, int voight2, double delta, double[] x, double eps) {
518 
519     // Store current unit cell parameters.
520     double a = unitCell.a;
521     double b = unitCell.b;
522     double c = unitCell.c;
523     double alpha = unitCell.alpha;
524     double beta = unitCell.beta;
525     double gamma = unitCell.gamma;
526 
527     // F(1,1)
528     double[][] dStrain = new double[3][3];
529     applyStrain(voight1, delta, dStrain);
530     applyStrain(voight2, delta, dStrain);
531     try {
532       if (!crystal.perturbCellVectors(dStrain)) {
533         logger.info(" Crystal method perturbCellVectors returned false.");
534         return 0.0;
535       }
536     } catch (Exception e) {
537       logger.info(" Exception from Crystal method perturbCellVectors.");
538       return 0.0;
539     }
540     forceFieldEnergy.setCrystal(crystal);
541     molecularAssembly.moveToFractionalCoordinates();
542     if (eps > 0.0) {
543       Minimize minimize = new Minimize(molecularAssembly, forceFieldEnergy, null);
544       minimize.minimize(eps);
545     }
546     forceFieldEnergy.getCoordinates(x);
547     double e11 = forceFieldEnergy.energy(x);
548 
549     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
550 
551     // F(-1,-1)
552     dStrain = new double[3][3];
553     applyStrain(voight1, -delta, dStrain);
554     applyStrain(voight2, -delta, dStrain);
555     try {
556       if (!crystal.perturbCellVectors(dStrain)) {
557         logger.info(" Crystal method perturbCellVectors returned false.");
558         return 0.0;
559       }
560     } catch (Exception e) {
561       logger.info(" Exception from Crystal method perturbCellVectors.");
562       return 0.0;
563     }
564     forceFieldEnergy.setCrystal(crystal);
565     molecularAssembly.moveToFractionalCoordinates();
566     if (eps > 0.0) {
567       Minimize minimize = new Minimize(molecularAssembly, forceFieldEnergy, null);
568       minimize.minimize(eps);
569     }
570     forceFieldEnergy.getCoordinates(x);
571     double em1m1 = forceFieldEnergy.energy(x);
572 
573     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
574 
575     // F(1,-1)
576     dStrain = new double[3][3];
577     applyStrain(voight1, delta, dStrain);
578     applyStrain(voight2, -delta, dStrain);
579     try {
580       if (!crystal.perturbCellVectors(dStrain)) {
581         logger.info(" Crystal method perturbCellVectors returned false.");
582         return 0.0;
583       }
584     } catch (Exception e) {
585       logger.info(" Exception from Crystal method perturbCellVectors.");
586       return 0.0;
587     }
588     forceFieldEnergy.setCrystal(crystal);
589     molecularAssembly.moveToFractionalCoordinates();
590     if (eps > 0) {
591       Minimize minimize = new Minimize(molecularAssembly, forceFieldEnergy, null);
592       minimize.minimize(eps);
593     }
594     forceFieldEnergy.getCoordinates(x);
595     double e1m1 = forceFieldEnergy.energy(x);
596 
597     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
598 
599     // F(-1,1)
600     dStrain = new double[3][3];
601     applyStrain(voight1, -delta, dStrain);
602     applyStrain(voight2, delta, dStrain);
603     try {
604       if (!crystal.perturbCellVectors(dStrain)) {
605         logger.info(" Crystal method perturbCellVectors returned false.");
606         return 0.0;
607       }
608     } catch (Exception e) {
609       logger.info(" Exception from Crystal method perturbCellVectors.");
610       return 0.0;
611     }
612     forceFieldEnergy.setCrystal(crystal);
613     molecularAssembly.moveToFractionalCoordinates();
614     if (eps > 0) {
615       Minimize minimize = new Minimize(molecularAssembly, forceFieldEnergy, null);
616       minimize.minimize(eps);
617     }
618     forceFieldEnergy.getCoordinates(x);
619     double em11 = forceFieldEnergy.energy(x);
620 
621     // Change back to original unit cell parameters.
622     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
623     forceFieldEnergy.setCrystal(crystal);
624     molecularAssembly.moveToFractionalCoordinates();
625     forceFieldEnergy.getCoordinates(x);
626 
627     return (e11 - e1m1 - em11 + em1m1) / (4 * delta * delta);
628   }
629 }