OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
elasticConstants.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include <cmath>
46#include <cstdio>
47#include <cstdlib>
48#include <cstring>
49#include <fstream>
50#include <iostream>
51#include <map>
52#include <memory>
53#include <string>
54
55#include "brains/ForceManager.hpp"
56#include "brains/Register.hpp"
57#include "brains/SimCreator.hpp"
58#include "brains/SimInfo.hpp"
59#include "brains/Thermo.hpp"
60#include "brains/Velocitizer.hpp"
61#include "constraints/Shake.hpp"
62#include "elasticConstantsCmd.hpp"
63#include "flucq/FluctuatingChargeConstraints.hpp"
64#include "flucq/FluctuatingChargeDamped.hpp"
65#include "io/DumpReader.hpp"
66#include "math/LU.hpp"
67#include "math/QR.hpp"
68#include "math/SquareMatrix3.hpp"
69#include "optimization/BoxObjectiveFunction.hpp"
70#include "optimization/Constraint.hpp"
71#include "optimization/Method.hpp"
72#include "optimization/OptimizationFactory.hpp"
73#include "optimization/Problem.hpp"
74#include "utils/StringUtils.hpp"
75
76using namespace OpenMD;
77using namespace JAMA;
78#include <algorithm>
79#include <iomanip>
80#include <iostream>
81#include <numeric>
82#include <vector>
83
84template<typename Real>
85class Vector6 : public Vector<Real, 6> {
86public:
87 using ElemType = Real;
88 using ElemPoinerType = Real*;
89
90 Vector6() : Vector<Real, 6>() {}
91
92 /** Constructs and initializes a Vector6d from individual coordinates */
93 inline Vector6(RealType v0, RealType v1, RealType v2, RealType v3,
94 RealType v4, RealType v5) {
95 this->data_[0] = v0;
96 this->data_[1] = v1;
97 this->data_[2] = v2;
98 this->data_[3] = v3;
99 this->data_[4] = v4;
100 this->data_[5] = v5;
101 }
102 /** Constructs and initializes from an array*/
103 inline Vector6(Real* array) : Vector<Real, 6>(array) {}
104
105 inline Vector6(const Vector<Real, 6>& v) : Vector<Real, 6>(v) {}
106
107 inline Vector6<Real>& operator=(const Vector<Real, 6>& v) {
108 if (this == &v) { return *this; }
109 Vector<Real, 6>::operator=(v);
110 return *this;
111 }
112};
113
114using Vector6d = Vector6<RealType>;
115
116RealType slope(const std::vector<RealType>& x, const std::vector<RealType>& y) {
117 // for (size_t i = 0; i < x.size(); i++) {
118 // std::cerr << x[i] << "\t" << y[i] << "\n";
119 // }
120 // std::cerr << "&\n";
121
122 const size_t n = x.size();
123 const RealType s_x = std::accumulate(x.begin(), x.end(), 0.0);
124 const RealType s_y = std::accumulate(y.begin(), y.end(), 0.0);
125 const RealType s_xx = std::inner_product(x.begin(), x.end(), x.begin(), 0.0);
126 const RealType s_xy = std::inner_product(x.begin(), x.end(), y.begin(), 0.0);
127 const RealType a = (n * s_xy - s_x * s_y) / (n * s_xx - s_x * s_x);
128 return a;
129}
130
131void quadraticFit(const std::vector<RealType>& x,
132 const std::vector<RealType>& y, RealType& a, RealType& b,
133 RealType& c) {
134 RealType s00 = RealType(x.size());
135 RealType s10(0.0), s20(0.0), s30(0.0), s40(0.0);
136 RealType s01(0.0), s11(0.0), s21(0.0);
137
138 for (size_t i = 0; i < x.size(); i++) {
139 // std::cerr << x[i] << "\t" << y[i] << "\n";
140 s10 += x[i];
141 s20 += pow(x[i], 2);
142 s30 += pow(x[i], 3);
143 s40 += pow(x[i], 4);
144 s01 += y[i];
145 s11 += x[i] * y[i];
146 s21 += pow(x[i], 2) * y[i];
147 }
148 // std::cerr << "&\n";
149
150 RealType D = (s40 * (s20 * s00 - s10 * s10) - s30 * (s30 * s00 - s10 * s20) +
151 s20 * (s30 * s10 - s20 * s20));
152
153 a = (s21 * (s20 * s00 - s10 * s10) - s11 * (s30 * s00 - s10 * s20) +
154 s01 * (s30 * s10 - s20 * s20)) /
155 D;
156
157 b = (s40 * (s11 * s00 - s01 * s10) - s30 * (s21 * s00 - s01 * s20) +
158 s20 * (s21 * s10 - s11 * s20)) /
159 D;
160
161 c = (s40 * (s20 * s01 - s10 * s11) - s30 * (s30 * s01 - s10 * s21) +
162 s20 * (s30 * s11 - s20 * s21)) /
163 D;
164
165 return;
166}
167
168void writeMatrix(DynamicRectMatrix<RealType> M, std::string title,
169 std::string units) {
170 std::cout << left << title << " (" << units << "):" << std::endl;
171 std::cout << std::endl;
172
173 std::cout << " ";
174 std::cout << right << setw(12) << M(0, 0) << " ";
175 std::cout << right << setw(12) << M(0, 1) << " ";
176 std::cout << right << setw(12) << M(0, 2) << " ";
177 std::cout << right << setw(12) << M(0, 3) << " ";
178 std::cout << right << setw(12) << M(0, 4) << " ";
179 std::cout << right << setw(12) << M(0, 5);
180 std::cout << std::endl;
181
182 std::cout << " ";
183 std::cout << right << setw(12) << M(1, 0) << " ";
184 std::cout << right << setw(12) << M(1, 1) << " ";
185 std::cout << right << setw(12) << M(1, 2) << " ";
186 std::cout << right << setw(12) << M(1, 3) << " ";
187 std::cout << right << setw(12) << M(1, 4) << " ";
188 std::cout << right << setw(12) << M(1, 5);
189 std::cout << std::endl;
190
191 std::cout << " ";
192 std::cout << right << setw(12) << M(2, 0) << " ";
193 std::cout << right << setw(12) << M(2, 1) << " ";
194 std::cout << right << setw(12) << M(2, 2) << " ";
195 std::cout << right << setw(12) << M(2, 3) << " ";
196 std::cout << right << setw(12) << M(2, 4) << " ";
197 std::cout << right << setw(12) << M(2, 5);
198 std::cout << std::endl;
199
200 std::cout << " ";
201 std::cout << right << setw(12) << M(3, 0) << " ";
202 std::cout << right << setw(12) << M(3, 1) << " ";
203 std::cout << right << setw(12) << M(3, 2) << " ";
204 std::cout << right << setw(12) << M(3, 3) << " ";
205 std::cout << right << setw(12) << M(3, 4) << " ";
206 std::cout << right << setw(12) << M(3, 5);
207 std::cout << std::endl;
208
209 std::cout << " ";
210 std::cout << right << setw(12) << M(4, 0) << " ";
211 std::cout << right << setw(12) << M(4, 1) << " ";
212 std::cout << right << setw(12) << M(4, 2) << " ";
213 std::cout << right << setw(12) << M(4, 3) << " ";
214 std::cout << right << setw(12) << M(4, 4) << " ";
215 std::cout << right << setw(12) << M(4, 5);
216 std::cout << std::endl;
217
218 std::cout << " ";
219 std::cout << right << setw(12) << M(5, 0) << " ";
220 std::cout << right << setw(12) << M(5, 1) << " ";
221 std::cout << right << setw(12) << M(5, 2) << " ";
222 std::cout << right << setw(12) << M(5, 3) << " ";
223 std::cout << right << setw(12) << M(5, 4) << " ";
224 std::cout << right << setw(12) << M(5, 5);
225 std::cout << std::endl;
226 std::cout << std::endl;
227}
228
229void writeBoxGeometries(Mat3x3d org, Mat3x3d opt, std::string title1,
230 std::string title2, std::string units) {
231 std::cout << left << title1 << " (" << units << "):" << std::endl;
232 std::cout << std::endl;
233 std::cout << " ";
234 std::cout << right << setw(12) << org(0, 0) << " ";
235 std::cout << right << setw(12) << org(0, 1) << " ";
236 std::cout << right << setw(12) << org(0, 2) << " ";
237 std::cout << std::endl;
238 std::cout << " ";
239 std::cout << right << setw(12) << org(1, 0) << " ";
240 std::cout << right << setw(12) << org(1, 1) << " ";
241 std::cout << right << setw(12) << org(1, 2) << " ";
242 std::cout << std::endl;
243 std::cout << " ";
244 std::cout << right << setw(12) << org(2, 0) << " ";
245 std::cout << right << setw(12) << org(2, 1) << " ";
246 std::cout << right << setw(12) << org(2, 2) << " ";
247 std::cout << std::endl;
248 std::cout << std::endl;
249 std::cout << left << title2 << " (" << units << "):" << std::endl;
250 std::cout << std::endl;
251 std::cout << " ";
252 std::cout << right << setw(12) << opt(0, 0) << " ";
253 std::cout << right << setw(12) << opt(0, 1) << " ";
254 std::cout << right << setw(12) << opt(0, 2) << " ";
255 std::cout << std::endl;
256 std::cout << " ";
257 std::cout << right << setw(12) << opt(1, 0) << " ";
258 std::cout << right << setw(12) << opt(1, 1) << " ";
259 std::cout << right << setw(12) << opt(1, 2) << " ";
260 std::cout << std::endl;
261 std::cout << " ";
262 std::cout << right << setw(12) << opt(2, 0) << " ";
263 std::cout << right << setw(12) << opt(2, 1) << " ";
264 std::cout << right << setw(12) << opt(2, 2) << " ";
265 std::cout << std::endl;
266 std::cout << std::endl;
267}
268
269void writeMaterialProperties(DynamicRectMatrix<RealType> C,
270 DynamicRectMatrix<RealType> S) {
271 RealType C11 = C(0, 0);
272 RealType C22 = C(1, 1);
273 RealType C33 = C(2, 2);
274 RealType C12 = C(0, 1);
275 RealType C23 = C(1, 2);
276 RealType C31 = C(2, 0);
277 RealType C44 = C(3, 3);
278 RealType C55 = C(4, 4);
279 RealType C66 = C(5, 5);
280
281 RealType S11 = S(0, 0);
282 RealType S22 = S(1, 1);
283 RealType S33 = S(2, 2);
284 RealType S12 = S(0, 1);
285 RealType S23 = S(1, 2);
286 RealType S31 = S(2, 0);
287 RealType S44 = S(3, 3);
288 RealType S55 = S(4, 4);
289 RealType S66 = S(5, 5);
290
291 // Material Properties defined in:
292 //
293 // "Charting the complete elastic properties of inorganic
294 // crystalline compounds," Maarten de Jong, Wei Chen, Thomas
295 // Angsten, Anubhav Jain, Randy Notestine, Anthony Gamst, Marcel
296 // Sluiter, Chaitanya Krishna Ande, Sybrand van der Zwaag, Jose J
297 // Plata, Cormac Toher, Stefano Curtarolo, Gerbrand Ceder, Kristin
298 // A. Persson & Mark Asta,
299 //
300 // Scientific Data volume 2, Article number: 150009 (2015)
301 // doi:10.1038/sdata.2015.9
302 //
303 // And in
304 //
305 // "ELATE: an open-source online application for analysis and
306 // visualization of elastic tensors," Romain Gaillac, Pluton
307 // Pullumbi and François-Xavier Coudert,
308 //
309 // J. Phys.: Condens. Matter 28 275201 (2016)
310 // doi:10.1088/0953-8984/28/27/275201
311
312 // Bulk modulus (Voigt average):
313 RealType Kv = ((C11 + C22 + C33) + 2.0 * (C12 + C23 + C31)) / 9.0;
314 // Bulk modulus (Reuss average):
315 RealType Kr = 1.0 / ((S11 + S22 + S33) + 2.0 * (S12 + S23 + S31));
316 // Shear modulus (Voigt average):
317 RealType Gv =
318 ((C11 + C22 + C33) - (C12 + C23 + C31) + 3.0 * (C44 + C55 + C66)) / 15.0;
319 // Shear modulus (Reuss average):
320 RealType Gr = 15.0 / (4.0 * (S11 + S22 + S33) - 4.0 * (S12 + S23 + S31) +
321 3.0 * (S44 + S55 + S66));
322 // Bulk modulus (Hill average):
323 RealType Kh = (Kv + Kr) / 2.0;
324 // Shear modulus (Hill average):
325 RealType Gh = (Gv + Gr) / 2.0;
326 // Universal elastic anisotropy
327 RealType Au = 5.0 * (Gv / Gr) + (Kv / Kr) - 6.0;
328 // Isotropic Poisson ratio
329 RealType muv = (1.0 - 3.0 * Gv / (3.0 * Kv + Gv)) / 2.0;
330 RealType mur = (1.0 - 3.0 * Gr / (3.0 * Kr + Gr)) / 2.0;
331 RealType muh = (1.0 - 3.0 * Gh / (3.0 * Kh + Gh)) / 2.0;
332 // Isotropic Young's modulus
333 RealType Ev = 1.0 / (1.0 / (3.0 * Gv) + 1.0 / (9.0 * Kv));
334 RealType Er = 1.0 / (1.0 / (3.0 * Gr) + 1.0 / (9.0 * Kr));
335 RealType Eh = 1.0 / (1.0 / (3.0 * Gh) + 1.0 / (9.0 * Kh));
336
337 std::cout << " " << setw(12) << "Voigt"
338 << " " << setw(12) << "Reuss"
339 << " " << setw(12) << "Hill\n";
340 std::cout << "Bulk modulus: " << setw(12) << Kv << " "
341 << setw(12) << Kr << " " << setw(12) << Kh << " (GPa)\n";
342
343 std::cout << "Shear modulus: " << setw(12) << Gv << " "
344 << setw(12) << Gr << " " << setw(12) << Gh << " (GPa)\n";
345
346 std::cout << "Young\'s modulus (isotropic): " << setw(12) << Ev << " "
347 << setw(12) << Er << " " << setw(12) << Eh << " (GPa)\n";
348
349 std::cout << "Poisson\'s Ratio: " << setw(12) << muv << " "
350 << setw(12) << mur << " " << setw(12) << muh << "\n";
351
352 std::cout << "Universal elastic Anisotropy: " << setw(12) << Au << "\n";
353
354 // Assume a cubic crystal, and use symmetries:
355
356 // C11 = (C(0,0) + C(1,1) + C(2,2)) / 3.0;
357 // C12 = (C(0,1) + C(0,2) + C(1,2)) / 3.0;
358 // C44 = (C(3,3) + C(4,4) + C(5,5)) / 3.0;
359 // S11 = (S(0,0) + S(1,1) + S(2,2)) / 3.0;
360 // S12 = (S(0,1) + S(0,2) + S(1,2)) / 3.0;
361 // S44 = (S(3,3) + S(4,4) + S(5,5)) / 3.0;
362
363 // Anisotropy factor
364 // RealType A1 = 2.0*C44 / (C11 - C12);
365 // RealType A2 = 2.0*(S11 - S12) / S44;
366
367 // std::cout << "Anisotropy factor = " << A1 << " " << A2 << "\n";
368
369 // Effective Elastic constants for propagation in Cubic Crystals
370 // RealType kL_100 = C11;
371 // RealType kT_100 = C44;
372 // RealType kL_110 = 0.5 * (C11 + C12 + 2.0*C44);
373 // RealType kT1_110 = C44;
374 // RealType kT2_110 = 0.5*(C11 - C12);
375 // RealType kL_111 = (C11 + 2*C12 + 4*C44) / 3.0;
376 // RealType kT_111 = (C11 - C12 + C44) / 3.0;
377}
378
379int main(int argc, char* argv[]) {
380 std::string method;
381 std::string inputFileName;
382 std::string outputFileName;
383
384 gengetopt_args_info args_info;
385
386 // parse command line arguments
387 if (cmdline_parser(argc, argv, &args_info) != 0) {
388 cmdline_parser_print_help();
389 exit(1);
390 }
391
392 // get input file name
393 if (args_info.input_given) {
394 inputFileName = args_info.input_arg;
395 } else {
396 if (args_info.inputs_num) {
397 inputFileName = args_info.inputs[0];
398 } else {
399 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
400 "No input file name was specified on the command line");
401 painCave.severity = OPENMD_ERROR;
402 painCave.isFatal = 1;
403 simError();
404 }
405 }
406
407 method = "energy";
408 if (args_info.method_given) {
409 method = args_info.method_arg;
410 toLower(method);
411 }
412
413 int nMax = args_info.npoints_arg;
414
415 RealType dmax;
416 if (args_info.delta_given) {
417 dmax = args_info.delta_arg;
418 } else {
419 if (!method.compare("energy")) {
420 dmax = 0.09;
421 } else {
422 dmax = 0.003;
423 }
424 }
425
426 // The strain basis sets are originally from:
427
428 // "Calculations of single-crystal elastic constants made simple,"
429 // R. Yu, J. Zhu, H.Q. Ye,
430 // Computer Physics Communications 181 (2010) 671–675,
431 // DOI: 10.1016/j.cpc.2009.11.017
432 //
433 // and
434 //
435 // "ElaStic: A tool for calculating second-order elastic
436 // constants from first principles," Rostam Golesorkhtabara,
437 // Pasquale Pavonea, Jürgen Spitalera, Peter Puschniga, Claudia Draxl,
438 // Computer Physics Communications 184 (2013) 1861–1873,
439 // DOI: 10.1016/j.cpc.2013.03.010
440 //
441 // Note that this our version assumes the worst about the crystal
442 // system present the box (e.g. a Triclinic box in the N Laue
443 // group).
444
445 std::vector<Vector6d> eStrains;
446 // Only the strains for the "N" Laue group are used (1-21, skipping 0)
447 // eStrains.push_back(Vector6d( 1., 1., 1., 0., 0., 0.));
448 eStrains.push_back(Vector6d(1., 0., 0., 0., 0., 0.));
449 eStrains.push_back(Vector6d(0., 1., 0., 0., 0., 0.));
450 eStrains.push_back(Vector6d(0., 0., 1., 0., 0., 0.));
451 eStrains.push_back(Vector6d(0., 0., 0., 2., 0., 0.));
452 eStrains.push_back(Vector6d(0., 0., 0., 0., 2., 0.));
453 eStrains.push_back(Vector6d(0., 0., 0., 0., 0., 2.));
454 eStrains.push_back(Vector6d(1., 1., 0., 0., 0., 0.));
455 eStrains.push_back(Vector6d(1., 0., 1., 0., 0., 0.));
456 eStrains.push_back(Vector6d(1., 0., 0., 2., 0., 0.));
457 eStrains.push_back(Vector6d(1., 0., 0., 0., 2., 0.));
458 eStrains.push_back(Vector6d(1., 0., 0., 0., 0., 2.));
459 eStrains.push_back(Vector6d(0., 1., 1., 0., 0., 0.));
460 eStrains.push_back(Vector6d(0., 1., 0., 2., 0., 0.));
461 eStrains.push_back(Vector6d(0., 1., 0., 0., 2., 0.));
462 eStrains.push_back(Vector6d(0., 1., 0., 0., 0., 2.));
463 eStrains.push_back(Vector6d(0., 0., 1., 2., 0., 0.));
464 eStrains.push_back(Vector6d(0., 0., 1., 0., 2., 0.));
465 eStrains.push_back(Vector6d(0., 0., 1., 0., 0., 2.));
466 eStrains.push_back(Vector6d(0., 0., 0., 2., 2., 0.));
467 eStrains.push_back(Vector6d(0., 0., 0., 2., 0., 2.));
468 eStrains.push_back(Vector6d(0., 0., 0., 0., 2., 2.));
469
470 // The rest (22-28) are only used for crystals of higher symmetry,
471 // and are not utilized in this code:
472 //
473 // eStrains.push_back(Vector6d( 0., 0., 0., 2., 2., 2.));
474 // eStrains.push_back(Vector6d(-1., .5, .5, 0., 0., 0.));
475 // eStrains.push_back(Vector6d( .5,-1., .5, 0., 0., 0.));
476 // eStrains.push_back(Vector6d( .5, .5,-1., 0., 0., 0.));
477 // eStrains.push_back(Vector6d( 1.,-1., 0., 0., 0., 0.));
478 // eStrains.push_back(Vector6d( 1.,-1., 0., 0., 0., 2.));
479 // eStrains.push_back(Vector6d( 0., 1.,-1., 0., 0., 2.));
480 // eStrains.push_back(Vector6d( .5, .5,-1., 0., 0., 2.));
481 // eStrains.push_back(Vector6d( 1., 0., 0., 2., 2., 0.));
482 // eStrains.push_back(Vector6d( 1., 1.,-1., 0., 0., 0.));
483 // eStrains.push_back(Vector6d( 1., 1., 1.,-2.,-2.,-2.));
484 // eStrains.push_back(Vector6d( .5, .5,-1., 2., 2., 2.));
485 // eStrains.push_back(Vector6d( 0., 0., 0., 2., 2., 4.));
486
487 // The Universal Linear-Independent Coupling Strains (ULICS) are from:
488 //
489 // "Calculations of single-crystal elastic constants made simple,"
490 // R. Yu, J. Zhu, H.Q. Ye,
491 // Computer Physics Communications 181 (2010) 671–675,
492 // DOI: 10.1016/j.cpc.2009.11.017
493
494 std::vector<Vector6d> sStrains;
495 sStrains.push_back(Vector6d(1., 2., 3., 4., 5., 6.));
496 sStrains.push_back(Vector6d(-2., 1., 4., -3., 6., -5.));
497 sStrains.push_back(Vector6d(3., -5., -1., 6., 2., -4.));
498 sStrains.push_back(Vector6d(-4., -6., 5., 1., -3., 2.));
499 sStrains.push_back(Vector6d(5., 4., 6., -2., -1., -3.));
500 sStrains.push_back(Vector6d(-6., 3., -2., 5., -4., 1.));
501
502 std::vector<Vector6d> strainBasis;
503
504 if (!method.compare("energy")) {
505 strainBasis = eStrains;
506 } else {
507 strainBasis = sStrains;
508 }
509
510 // The matrix to perform the linear least squares fits from the
511 // ULICS to find the elastic constants were derived for the ElaStic
512 // code, Golesorkhtabara, et al. (cited above):
513
514 RealType mat[36][21] = {
515 {1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
516 {0, 1, 0, 0, 0, 0, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
517 {0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0},
518 {0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 4, 5, 6, 0, 0, 0},
519 {0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 4, 0, 5, 6, 0},
520 {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 4, 0, 5, 6},
521 {-2, 1, 4, -3, 6, -5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
522 {0, -2, 0, 0, 0, 0, 1, 4, -3, 6, -5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
523 {0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 4, -3, 6, -5, 0, 0, 0, 0, 0, 0},
524 {0, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 4, 0, 0, -3, 6, -5, 0, 0, 0},
525 {0, 0, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 4, 0, 0, -3, 0, 6, -5, 0},
526 {0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 0, 4, 0, 0, -3, 0, 6, -5},
527 {3, -5, -1, 6, 2, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
528 {0, 3, 0, 0, 0, 0, -5, -1, 6, 2, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
529 {0, 0, 3, 0, 0, 0, 0, -5, 0, 0, 0, -1, 6, 2, -4, 0, 0, 0, 0, 0, 0},
530 {0, 0, 0, 3, 0, 0, 0, 0, -5, 0, 0, 0, -1, 0, 0, 6, 2, -4, 0, 0, 0},
531 {0, 0, 0, 0, 3, 0, 0, 0, 0, -5, 0, 0, 0, -1, 0, 0, 6, 0, 2, -4, 0},
532 {0, 0, 0, 0, 0, 3, 0, 0, 0, 0, -5, 0, 0, 0, -1, 0, 0, 6, 0, 2, -4},
533 {-4, -6, 5, 1, -3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
534 {0, -4, 0, 0, 0, 0, -6, 5, 1, -3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
535 {0, 0, -4, 0, 0, 0, 0, -6, 0, 0, 0, 5, 1, -3, 2, 0, 0, 0, 0, 0, 0},
536 {0, 0, 0, -4, 0, 0, 0, 0, -6, 0, 0, 0, 5, 0, 0, 1, -3, 2, 0, 0, 0},
537 {0, 0, 0, 0, -4, 0, 0, 0, 0, -6, 0, 0, 0, 5, 0, 0, 1, 0, -3, 2, 0},
538 {0, 0, 0, 0, 0, -4, 0, 0, 0, 0, -6, 0, 0, 0, 5, 0, 0, 1, 0, -3, 2},
539 {5, 4, 6, -2, -1, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
540 {0, 5, 0, 0, 0, 0, 4, 6, -2, -1, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
541 {0, 0, 5, 0, 0, 0, 0, 4, 0, 0, 0, 6, -2, -1, -3, 0, 0, 0, 0, 0, 0},
542 {0, 0, 0, 5, 0, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, -2, -1, -3, 0, 0, 0},
543 {0, 0, 0, 0, 5, 0, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, -2, 0, -1, -3, 0},
544 {0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, -2, 0, -1, -3},
545 {-6, 3, -2, 5, -4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
546 {0, -6, 0, 0, 0, 0, 3, -2, 5, -4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
547 {0, 0, -6, 0, 0, 0, 0, 3, 0, 0, 0, -2, 5, -4, 1, 0, 0, 0, 0, 0, 0},
548 {0, 0, 0, -6, 0, 0, 0, 0, 3, 0, 0, 0, -2, 0, 0, 5, -4, 1, 0, 0, 0},
549 {0, 0, 0, 0, -6, 0, 0, 0, 0, 3, 0, 0, 0, -2, 0, 0, 5, 0, -4, 1, 0},
550 {0, 0, 0, 0, 0, -6, 0, 0, 0, 0, 3, 0, 0, 0, -2, 0, 0, 5, 0, -4, 1}};
551
552 DynamicRectMatrix<RealType> drMat = DynamicRectMatrix<RealType>(36, 21, *mat);
553
554 // register forcefields, integrators and minimizers
555 registerAll();
556
557 // Parse the input file, set up the system, and read the last frame:
558 SimCreator creator;
559 SimInfo* info = creator.createSim(inputFileName, true);
560 Globals* simParams = info->getSimParams();
561 ForceManager* forceMan = new ForceManager(info);
562
563 std::unique_ptr<Velocitizer> veloSet {std::make_unique<Velocitizer>(info)};
564
565 forceMan->initialize();
566 info->update();
567
568 Shake* shake = new Shake(info);
569 bool hasFlucQ = false;
570 FluctuatingChargePropagator* flucQ = new FluctuatingChargeDamped(info);
571
572 if (info->usesFluctuatingCharges()) {
573 if (info->getNFluctuatingCharges() > 0) {
574 hasFlucQ = true;
575 flucQ->setForceManager(forceMan);
576 flucQ->initialize();
577 }
578 }
579
580 // Important utility classes for computing system properties:
581 Thermo thermo(info);
582
583 // Just in case we were passed a system that is on the move:
584 veloSet->removeComDrift();
585
586 if (args_info.box_flag) {
587 std::cout << "Doing box optimization\n\n";
588 Mat3x3d oldHmat =
589 info->getSnapshotManager()->getCurrentSnapshot()->getHmat();
590
591 MinimizerParameters* miniPars = simParams->getMinimizerParameters();
592 // OptimizationMethod* minim =
593 // OptimizationFactory::getInstance().createOptimization(toUpperCopy(miniPars->getMethod()),
594 // info);
595 OptimizationMethod* minim =
596 OptimizationFactory::getInstance().createOptimization("CG", info);
597
598 if (minim == NULL) {
599 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
600 "Optimization Factory can not create %s OptimizationMethod\n",
601 miniPars->getMethod().c_str());
602 painCave.isFatal = 1;
603 simError();
604 }
605
606 BoxObjectiveFunction boxObjf(info, forceMan);
607 NoConstraint noConstraint {};
608 DumpStatusFunction dsf(info);
609 DynamicVector<RealType> initCoords = boxObjf.setInitialCoords();
610 Problem problem(boxObjf, noConstraint, dsf, initCoords);
611
612 int maxIter = miniPars->getMaxIterations();
613 int mssIter = miniPars->getMaxStationaryStateIterations();
614 RealType rEps = miniPars->getRootEpsilon();
615 RealType fEps = miniPars->getFunctionEpsilon();
616 RealType gnEps = miniPars->getGradientNormEpsilon();
617 RealType initialStepSize = miniPars->getInitialStepSize();
618
619 EndCriteria endCriteria(maxIter, mssIter, rEps, fEps, gnEps);
620
621 minim->minimize(problem, endCriteria, initialStepSize);
622 delete minim;
623
624 Mat3x3d newHmat =
625 info->getSnapshotManager()->getCurrentSnapshot()->getHmat();
626 writeBoxGeometries(oldHmat, newHmat, "Original Box Geometry",
627 "Optimized Box Geometry", "Angstroms");
628 }
629
630 Snapshot* snap = info->getSnapshotManager()->getCurrentSnapshot();
631 Mat3x3d refHmat = snap->getHmat();
632
633 forceMan->calcForces();
634 Mat3x3d ptRef = thermo.getPressureTensor();
635 RealType V0 = thermo.getVolume();
636 ptRef.negate();
637 ptRef *= Constants::elasticConvert;
638
639 Vector6d stress(0.0);
640 Vector6d strain(0.0);
641 Vector6d lstress(0.0);
642 Mat3x3d epsilon(0.0);
643
644 std::vector<std::vector<RealType>> stressStrain;
645 std::vector<RealType> strainValues;
646 std::vector<RealType> energyValues;
647 DynamicRectMatrix<RealType> C(6, 6, 0.0);
648 DynamicRectMatrix<RealType> S(6, 6, 0.0);
649 DynamicVector<RealType> ci(21, 0.0);
650 DynamicVector<RealType> sigma(36, 0.0);
651
652 Vector3d pos;
653 SimInfo::MoleculeIterator miter;
654 Molecule* mol;
655 RealType de;
656 Vector3d delta;
657
658 Vector6d L(0.0);
659 Mat3x3d eta(0.0);
660 Mat3x3d eps(0.0);
661 Mat3x3d x(0.0);
662 Mat3x3d test(0.0);
663 Mat3x3d deformation(0.0);
664 std::vector<RealType> A2;
665 RealType norm;
666 RealType a, b, c;
667 RealType energy;
668 Mat3x3d pressureTensor;
669
670 for (std::vector<Vector6d>::iterator it = strainBasis.begin();
671 it != strainBasis.end(); ++it) {
672 strain = *it;
673 int ii = std::distance(strainBasis.begin(), it);
674
675 strainValues.clear();
676 energyValues.clear();
677 stressStrain.clear();
678 stressStrain.resize(6);
679
680 for (int n = 0; n < nMax; n++) {
681 // First, set up the deformation of the box and coodinates:
682 de = -0.5 * dmax + dmax * RealType(n) / RealType(nMax - 1);
683 L = strain * de;
684
685 // η is the Lagrangian strain tensor:
686 eta.setupVoigtTensor(L[0], L[1], L[2], L[3] / 2., L[4] / 2., L[5] / 2.);
687
688 // Make sure the deformation isn't too large:
689 if (eta.frobeniusNorm() > 0.7) {
690 std::cerr << "Deformation is too large!\n";
691 }
692
693 // Find the physical strain tensor, ε, from the Lagrangian strain, η:
694 // η = ε + 0.5 * ε^2
695 norm = 1.0;
696 eps = eta;
697 while (norm > 1.0e-10) {
698 x = eta - eps * eps / 2.0;
699 test = x - eps;
700 norm = test.frobeniusNorm();
701 eps = x;
702 }
703 deformation = SquareMatrix3<RealType>::identity() + eps;
704
705 // Second, do the deformation and compute the energy or stress
706 // tensor for this deformation:
707 info->getSnapshotManager()->advance();
708 for (mol = info->beginMolecule(miter); mol != NULL;
709 mol = info->nextMolecule(miter)) {
710 pos = mol->getCom();
711 delta = deformation * pos;
712 mol->moveCom(delta - pos);
713 }
714 Mat3x3d Hmat = deformation * refHmat;
715 snap->setHmat(Hmat);
716 shake->constraintR();
717 forceMan->calcForces();
718 if (hasFlucQ) flucQ->applyConstraints();
719 shake->constraintF();
720
721 // Third, record the energy or the stress:
722 if (!method.compare("energy")) {
723 energy = thermo.getPotential();
724 energyValues.push_back(energy);
725
726 } else {
727 // Find the Lagragian stress tensor, τ, from the physical
728 // stress tensor, σ, that was computed from the pressureTensor
729 // in this code.
730 // τ = det(1+ε) (1+ε)^−1 · σ · (1+ε)^−1
731 // (Note that 1+ε is the deformation tensor computed above.)
732
733 Mat3x3d idm = deformation.inverse();
734 RealType ddm = deformation.determinant();
735
736 pressureTensor = thermo.getPressureTensor();
737 pressureTensor.negate();
738 pressureTensor *= Constants::elasticConvert;
739
740 Mat3x3d tao = idm * (pressureTensor * idm);
741 tao *= ddm;
742
743 lstress = tao.toVoigtTensor();
744
745 for (int j = 0; j < 6; j++) {
746 stressStrain[j].push_back(lstress[j]);
747 }
748 }
749
750 strainValues.push_back(de);
751 info->getSnapshotManager()->resetToPrevious();
752 }
753
754 // Fit the energy vs. strain (quadratic) or stress vs. strain (linear)
755 if (!method.compare("energy")) {
756 quadraticFit(strainValues, energyValues, a, b, c);
757 A2.push_back(a * Constants::energyElasticConvert / V0);
758 } else {
759 for (int j = 0; j < 6; j++) {
760 quadraticFit(strainValues, stressStrain[j], a, b, c);
761 sigma(6 * ii + j) = b;
762 // sigma(6*ii + j) = slope(strainValues, stressStrain[j]);
763 }
764 }
765 }
766
767 if (!method.compare("energy")) {
768 C(0, 0) = 2. * A2[0];
769 C(0, 1) = 1. * (-A2[0] - A2[1] + A2[6]);
770 C(0, 2) = 1. * (-A2[0] - A2[2] + A2[7]);
771 C(0, 3) = .5 * (-A2[0] - A2[3] + A2[8]);
772 C(0, 4) = .5 * (-A2[0] + A2[9] - A2[4]);
773 C(0, 5) = .5 * (-A2[0] + A2[10] - A2[5]);
774 C(1, 1) = 2. * A2[1];
775 C(1, 2) = 1. * (A2[11] - A2[1] - A2[2]);
776 C(1, 3) = .5 * (A2[12] - A2[1] - A2[3]);
777 C(1, 4) = .5 * (A2[13] - A2[1] - A2[4]);
778 C(1, 5) = .5 * (A2[14] - A2[1] - A2[5]);
779 C(2, 2) = 2. * A2[2];
780 C(2, 3) = .5 * (A2[15] - A2[2] - A2[3]);
781 C(2, 4) = .5 * (A2[16] - A2[2] - A2[4]);
782 C(2, 5) = .5 * (A2[17] - A2[2] - A2[5]);
783 C(3, 3) = .5 * A2[3];
784 C(3, 4) = .25 * (A2[18] - A2[3] - A2[4]);
785 C(3, 5) = .25 * (A2[19] - A2[3] - A2[5]);
786 C(4, 4) = .5 * A2[4];
787 C(4, 5) = .25 * (A2[20] - A2[4] - A2[5]);
788 C(5, 5) = .5 * A2[5];
789 } else {
790 // Least squares to map fits of stress-strain relationships onto
791 // elastic matrix:
792
793 QR<RealType> qr(drMat);
794 ci = qr.solve(sigma);
795
796 C(0, 0) = ci(0);
797 C(0, 1) = ci(1);
798 C(0, 2) = ci(2);
799 C(0, 3) = ci(3);
800 C(0, 4) = ci(4);
801 C(0, 5) = ci(5);
802 C(1, 1) = ci(6);
803 C(1, 2) = ci(7);
804 C(1, 3) = ci(8);
805 C(1, 4) = ci(9);
806 C(1, 5) = ci(10);
807 C(2, 2) = ci(11);
808 C(2, 3) = ci(12);
809 C(2, 4) = ci(13);
810 C(2, 5) = ci(14);
811 C(3, 3) = ci(15);
812 C(3, 4) = ci(16);
813 C(3, 5) = ci(17);
814 C(4, 4) = ci(18);
815 C(4, 5) = ci(19);
816 C(5, 5) = ci(20);
817 }
818
819 // Symmetrize C:
820
821 for (int i = 0; i < 5; i++) {
822 for (int j = i + 1; j < 6; j++) {
823 C(j, i) = C(i, j);
824 }
825 }
826
827 // matrix is destroyed during inversion, so make a working copy:
828 DynamicRectMatrix<RealType> tmpMat(C);
829 invertMatrix(tmpMat, S);
830
831 writeMatrix(C, "Elastic Tensor", "GPa");
832 writeMatrix(S, "Compliance Tensor", "GPa^-1");
833
834 writeMaterialProperties(C, S);
835
836 return 0;
837}
Constraint.hpp
Abstract constraint class.
Method.hpp
Abstract optimization method class.
Problem.hpp
Abstract optimization problem class.
OpenMD::FluctuatingChargePropagator
abstract class for propagating fluctuating charge variables
Definition FluctuatingChargePropagator.hpp:59
OpenMD::ForceManager
ForceManager is responsible for calculating both the short range (bonded) interactions and long range...
Definition ForceManager.hpp:80
OpenMD::Molecule::moveCom
void moveCom(const Vector3d &delta)
Moves the center of this molecule.
Definition Molecule.cpp:341
OpenMD::Molecule::getCom
Vector3d getCom()
Returns the current center of mass position of this molecule.
Definition Molecule.cpp:298
OpenMD::OptimizationFactory::getInstance
static OptimizationFactory & getInstance()
Returns an instance of Optimization factory.
Definition OptimizationFactory.hpp:78
OpenMD::RectMatrix::negate
void negate()
Negates the value of this matrix in place.
Definition RectMatrix.hpp:242
OpenMD::SimCreator
The only responsibility of SimCreator is to parse the meta-data file and create a SimInfo instance ba...
Definition SimCreator.hpp:71
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
OpenMD::SimInfo::beginMolecule
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
OpenMD::SimInfo::update
void update()
update
Definition SimInfo.cpp:697
OpenMD::SimInfo::getNFluctuatingCharges
int getNFluctuatingCharges()
Returns the total number of fluctuating charges that are present.
Definition SimInfo.hpp:217
OpenMD::SimInfo::nextMolecule
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
OpenMD::SimInfo::getSnapshotManager
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
OpenMD::Snapshot
The Snapshot class is a repository storing dynamic data during a Simulation.
Definition Snapshot.hpp:147
OpenMD::Snapshot::getHmat
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
OpenMD::Snapshot::setHmat
void setHmat(const Mat3x3d &m)
Sets the H-Matrix.
Definition Snapshot.cpp:217
OpenMD::SnapshotManager::getCurrentSnapshot
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Definition SnapshotManager.hpp:89
OpenMD::SquareMatrix3::inverse
SquareMatrix3< Real > inverse() const
Sets the value of this matrix to the inverse of itself.
Definition SquareMatrix3.hpp:422
OpenMD::SquareMatrix3::determinant
Real determinant() const
Returns the determinant of this matrix.
Definition SquareMatrix3.hpp:399
OpenMD::Vector
Fix length vector class.
Definition Vector.hpp:78
OpenMD::Vector::operator=
Vector< Real, Dim > & operator=(const Vector< Real, Dim > &v)
copy assignment operator
Definition Vector.hpp:93
QuantLib::EndCriteria
Criteria to end optimization process:
Definition EndCriteria.hpp:42
QuantLib::OptimizationMethod
Abstract class for constrained optimization method.
Definition Method.hpp:36
QuantLib::Problem
Constrained optimization problem.
Definition Problem.hpp:37
elasticConstantsCmd.hpp
The header file for the command line option parser generated by GNU Gengetopt version 2....
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60
OpenMD::registerAll
void registerAll()
register force fields, integrators and optimizers
Definition Register.cpp:140
OpenMD::invertMatrix
bool invertMatrix(MatrixType &A, MatrixType &AI)
Invert input square matrix A into matrix AI.
Definition LU.hpp:98
gengetopt_args_info
Where the command line options are stored.
Definition Dump2XYZCmd.hpp:39