1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
59
60
61
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
73
74
75
76
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
100
101
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
116 double PRESCON = 6.85684112e4;
117
118
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
142
143
144
145
146
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
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
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
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
327 molecularAssembly.setFractionalMode(currentFractionalMode);
328 crystal.setCheckRestrictions(currentCheckRestrictions);
329 }
330
331
332
333
334
335
336
337
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
350
351
352 public void printTensor() {
353 computeStressTensor(true);
354 }
355
356
357
358
359
360
361
362
363 private void applyStrain(int voight, double delta, double[][] strain) {
364 switch (voight) {
365 case 1 ->
366 strain[0][0] += delta;
367
368
369
370 case 2 ->
371
372 strain[1][1] += delta;
373
374
375 case 3 ->
376
377
378 strain[2][2] += delta;
379 case 4 -> {
380 strain[1][2] += delta / 2.0;
381 strain[2][1] += delta / 2.0;
382 }
383 case 5 -> {
384 strain[0][2] += delta / 2.0;
385 strain[2][0] += delta / 2.0;
386 }
387 case 6 -> {
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
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
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
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
482 forceFieldEnergy.setCrystal(crystal);
483
484
485 molecularAssembly.moveToFractionalCoordinates();
486 forceFieldEnergy.getCoordinates(x);
487
488
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
501 forceFieldEnergy.setCrystal(crystal);
502
503
504 molecularAssembly.moveToFractionalCoordinates();
505 forceFieldEnergy.getCoordinates(x);
506
507
508 double eminus = forceFieldEnergy.energy(x);
509
510
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
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
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
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
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
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
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 }