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