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 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
58
59
60
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
72
73
74
75
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
99
100
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
115 double PRESCON = 6.85684112e4;
116
117
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
141
142
143
144
145
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
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
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
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
326 molecularAssembly.setFractionalMode(currentFractionalMode);
327 crystal.setCheckRestrictions(currentCheckRestrictions);
328 }
329
330
331
332
333
334
335
336
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
349
350
351 public void printTensor() {
352 computeStressTensor(true);
353 }
354
355
356
357
358
359
360
361
362 private void applyStrain(int voight, double delta, double[][] strain) {
363 switch (voight) {
364 case 1 ->
365 strain[0][0] += delta;
366
367
368
369 case 2 ->
370
371 strain[1][1] += delta;
372
373
374 case 3 ->
375
376
377 strain[2][2] += delta;
378 case 4 -> {
379 strain[1][2] += delta / 2.0;
380 strain[2][1] += delta / 2.0;
381 }
382 case 5 -> {
383 strain[0][2] += delta / 2.0;
384 strain[2][0] += delta / 2.0;
385 }
386 case 6 -> {
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
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
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
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
481 forceFieldEnergy.setCrystal(crystal);
482
483
484 molecularAssembly.moveToFractionalCoordinates();
485 forceFieldEnergy.getCoordinates(x);
486
487
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
500 forceFieldEnergy.setCrystal(crystal);
501
502
503 molecularAssembly.moveToFractionalCoordinates();
504 forceFieldEnergy.getCoordinates(x);
505
506
507 double eminus = forceFieldEnergy.energy(x);
508
509
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
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
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
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
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
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
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 }