Loading...
Searching...
No Matches
ccl.cpp
Go to the documentation of this file.
1/*
2** This file is part of eOn.
3**
4** SPDX-License-Identifier: BSD-3-Clause
5**
6** Copyright (c) 2010--present, eOn Development Team
7** All rights reserved.
8**
9** Repo:
10** https://github.com/TheochemUI/eOn
11*/
18
19#include "eon/potentials/Water/ccl.hpp"
20#include <cassert>
21#include <cmath>
22// #include "unit_system.hpp"
23
24using namespace forcefields;
25namespace {
26#if defined(FORCEFIELDS_UNIT_SYSTEM_HPP) && \
27 (FORCEFIELDS_UNIT_SYSTEM_HPP != \
28 FORCEFIELDS_UNIT_SYSTEM_ELECTRONVOLT_ANGSTROM_FEMTOSECOND_ECHARGE)
29using namespace unit_system;
30
31double const re_ = 0.9572 * ANGSTROM;
32double const thetae_ = 104.52 * DEGREE;
33
34double const re2_ = re_ * re_;
35// double const ERGS_PER_ANGSTROM2=ERGS/ANGSTROM2;
36
37// ------------------------ Quadratic ---------------------------
38double const ro_2_ = 84.54e-12 * ERGS_PER_ANGSTROM2;
39double const ro1_ro2_ = -1.01e-12 * ERGS_PER_ANGSTROM2;
40double const ro_theta_ = 2.288e-12 * ERGS_PER_ANGSTROM2;
41double const theta_2_ = 7.607e-12 * ERGS_PER_ANGSTROM2;
42
43// --------------------------------------------- Cubic
44// ----------------------------------------
45double const ro_3_ = -10.168e-12 * ERGS_PER_ANGSTROM2;
46double const ro_ro1_ro2_ = 0.201e-12 * ERGS_PER_ANGSTROM2;
47double const ro_2_theta_ = 4.308e-12 * ERGS_PER_ANGSTROM2;
48double const ro1_ro2_theta_ = -4.020e-12 * ERGS_PER_ANGSTROM2;
49double const ro_theta_2_ = -1.175e-12 * ERGS_PER_ANGSTROM2;
50double const thetat_3_ = -1.595e-12 * ERGS_PER_ANGSTROM2;
51
52// --------------------------------------------- Quartic
53// ----------------------------------------
54double const ro_4_ = -10.684e-12 * ERGS_PER_ANGSTROM2;
55double const ro1_ro2_ro_2_ = -6.162e-12 * ERGS_PER_ANGSTROM2;
56double const ro1_2_ro2_2_ = 2.717e-12 * ERGS_PER_ANGSTROM2;
57
58double const ro_3_theta_ = 6.328e-12 * ERGS_PER_ANGSTROM2;
59double const ro_ro1_ro2_theta_ = -4.020e-12 * ERGS_PER_ANGSTROM2;
60
61double const ro_2_theta_2_ = -4.70e-12 * ERGS_PER_ANGSTROM2;
62double const ro1_ro2_theta_2_ = 3.05e-12 * ERGS_PER_ANGSTROM2;
63
64// ro_theta_3 = 0
65double const theta_4_ = -0.0318e-12 * ERGS_PER_ANGSTROM2;
66#else
67double const re_ = 0.9572; // ANGSTROM
68double const thetae_ = 1.82421813418447321; // RADIANS
69
70double const re2_ = 0.91623184; // ANGSTROM^2
71// double const ERGS_PER_ANGSTROM2 = 624150947960.771851; // eV / Angstrom^2
72
73// ------------------------ Quadratic ---------------------------
74double const ro_2_ = 52.7657211406036524; // eV / Angstrom^2
75double const ro1_ro2_ = -0.630392457440379528; // eV / Angstrom^2
76double const ro_theta_ = 1.42805736893424595; // eV / Angstrom^2
77double const theta_2_ = 4.74791626113759158; // eV / Angstrom^2
78
79// --------------------------------------------- Cubic
80// ----------------------------------------
81double const ro_3_ = -6.3463668388651282; // eV / Angstrom^2
82double const ro_ro1_ro2_ = 0.125454340540115145; // eV / Angstrom^2
83double const ro_2_theta_ = 2.68884228381500501; // eV / Angstrom^2
84double const ro1_ro2_theta_ = -2.509086810802303; // eV / Angstrom^2
85double const ro_theta_2_ = -0.733377363853906838; // eV / Angstrom^2
86double const thetat_3_ = -0.995520761997431114; // eV / Angstrom^2
87
88// --------------------------------------------- Quartic
89// ----------------------------------------
90double const ro_4_ = -6.668428728012886; // eV / Angstrom^2
91double const ro1_ro2_ro_2_ = -3.84601814133427622; // eV / Angstrom^2
92double const ro1_2_ro2_2_ = 1.69581812560941692; // eV / Angstrom^2
93double const ro_3_theta_ = 3.94962719869576429; // eV / Angstrom^2
94double const ro_ro1_ro2_theta_ = -2.509086810802303; // eV / Angstrom^2
95double const ro_2_theta_2_ = -2.93350945541562735; // eV / Angstrom^2
96double const ro1_ro2_theta_2_ = 1.90366039128035425; // eV / Angstrom^2
97double const theta_4_ = -0.0198480001451525438; // eV / Angstrom^2
98#endif
99} // namespace
100
111
113double const Ccl::re_ = ::re_;
115double const Ccl::thetae_ = ::thetae_;
116
130void Ccl::computeHH_O_(const int nAtoms, const double R[], double F[],
131 double &U, const double b[]) {
132 computeHH_O_(nAtoms, R, F, U, b, 0);
133}
134
148void Ccl::computeHH_O_(const int nAtoms, const double R[], double F[],
149 double &U, const double b[], const bool fixed[]) {
150 int const nMolecules = nAtoms / 3;
151 const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(R);
152 const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&R[3]);
153 const double (*const ro)[3] =
154 reinterpret_cast<const double (*)[3]>(&R[nMolecules * 6]);
155 double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(F);
156 double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&F[3]);
157 double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&F[nMolecules * 6]);
158 bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
159 bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
160 bool const(*const xo)[1] =
161 reinterpret_cast<bool const(*)[1]>(&fixed[nMolecules * 2]);
162 computeTemplate(nMolecules, rh1, rh2, ro, fh1, fh2, fo, U, b, xh1, xh2, xo);
163}
164
166char const *Ccl::getName() const { return "Ccl"; }
167
172Ccl::Ccl(double cutoff, double switchingWidth)
173 : PotentialBase(cutoff, switchingWidth) {}
174
182
197
199void Ccl::initialiseRho(Vector3 const &v, Rho &ro) {
200 ro._1 = (v._1 - re_) / v._1;
201 ro._2 = ro._1 * ro._1;
202 ro._3 = ro._2 * ro._1;
203 double const d = re_ / v._2 / v._1;
204 for (int i = 0; i < 3; ++i)
205 ro.n[i] = v.v[i] * d;
206}
207
209void Ccl::initialiseDtheta(Vector3 const &v1, Vector3 const &v2, Dtheta &dth,
210 double const thetaEquilibrium) {
211 double const cos_theta = dotProduct(v1.v, v2.v) / v1._1 / v2._1;
212 double const theta = std::acos(cos_theta);
213 dth._1 = theta - thetaEquilibrium;
214 dth._2 = dth._1 * dth._1;
215 dth._3 = dth._2 * dth._1;
216 double const d_theta = -1.0 / std::sqrt(1.0 - cos_theta * cos_theta);
217 for (int k = 0; k < 3; ++k) {
218 dth.n1[k] = 0.0;
219 dth.n2[k] = 0.0;
220 for (int j = 0; j < 3; ++j) {
221 dth.n1[k] -= v2.v[j] / v2._1 / v1._1 * v1.v[j] * v1.v[k] / v1._2;
222 dth.n2[k] -= v1.v[j] / v1._1 / v2._1 * v2.v[j] * v2.v[k] / v2._2;
223 };
224 dth.n1[k] += v2.v[k] / v2._1 / v1._1;
225 dth.n2[k] += v1.v[k] / v1._1 / v2._1;
226 dth.n1[k] *= d_theta;
227 dth.n2[k] *= d_theta;
228 };
229}
230
238void Ccl::intramolecular(double const rh1[], double const rh2[],
239 double const ro[], double fh1[], double fh2[],
240 double fo[], double &energy) {
241 Vector3 v1, v2;
242 distance(rh1, ro, v1);
243 distance(rh2, ro, v2);
244 // ------prepare _
245 Rho ro1, ro2;
246 initialiseRho(v1, ro1);
247 initialiseRho(v2, ro2);
248
249 // ------- prepare Delta theta
250 Dtheta dth;
251 initialiseDtheta(v1, v2, dth, thetae_);
252 intramolecular(ro1, ro2, dth, energy, fh1, fh2, fo);
253}
254
260void Ccl::intramolecular(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
261 double &energy, double fh1[], double fh2[],
262 double fo[]) {
263 double d1 = 0, d2 = 0,
264 d3 = 0; // d1=d(E)/d(rho1), d2=d(E)/d(rho2), d(E)/d(theta)
265
266 // -------- quadratic interaction
267 ro_2(ro1, ro2, energy, d1, d2);
268 ro1_ro2(ro1, ro2, energy, d1, d2);
269 ro_theta(ro1, ro2, dth, energy, d1, d2, d3);
270 theta_2(dth, energy, d3);
271 //*/
272
273 // -------- Cubic interaction
274 ro_3(ro1, ro2, energy, d1, d2);
275 ro_ro1_ro2(ro1, ro2, energy, d1, d2);
276 ro_2_theta(ro1, ro2, dth, energy, d1, d2, d3);
277 ro1_ro2_theta(ro1, ro2, dth, energy, d1, d2, d3);
278 ro_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
279 theta_3(dth, energy, d3);
280 //*/
281
282 // --------------------------------------------- Quartic
283 // -------------------------------------
284 ro_4(ro1, ro2, energy, d1, d2);
285 ro1_ro2_ro_2(ro1, ro2, energy, d1, d2);
286 ro1_2_ro2_2(ro1, ro2, energy, d1, d2);
287
288 ro_3_theta(ro1, ro2, dth, energy, d1, d2, d3);
289 ro_ro1_ro2_theta(ro1, ro2, dth, energy, d1, d2, d3);
290
291 ro_2_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
292 ro1_ro2_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
293
294 // ro_theta_3 = 0
295 theta_4(dth, energy, d3);
296 //*/
297 for (int i = 0; i < 3; ++i) {
298 fh1[i] -= d1 * ro1.n[i] + d3 * dth.n1[i];
299 fh2[i] -= d2 * ro2.n[i] + d3 * dth.n2[i];
300 fo[i] += d1 * ro1.n[i] + d2 * ro2.n[i] + d3 * dth.n1[i] + d3 * dth.n2[i];
301 };
302}
303
322template <int H, int O>
323void Ccl::computeTemplate(
324 const int nMolecules, const double (*const rh1)[H * 3],
325 const double (*const rh2)[H * 3], const double (*const ro)[O * 3],
326 double (*const fh1)[H * 3], double (*const fh2)[H * 3],
327 double (*const fo)[O * 3], double &energy, double const b[],
328 bool const (*const xh1)[H], bool const (*const xh2)[H],
329 bool const (*const xo)[O]) {
330 for (int i = 0; i < nMolecules; ++i) {
331 for (int a = 0; a < 3; a++) {
332 fh1[i][a] = 0.0;
333 fh2[i][a] = 0.0;
334 fo[i][a] = 0.0;
335 };
336 };
337 energy = 0.0;
338 setPeriodicity(b);
339
340 for (int i = nMolecules - 1; i >= 0; --i) {
341 bool const areFixed = xh1[i][0] and xh2[i][0] and xo[i][0];
342 if (not areFixed) {
343 intramolecular(rh1[i], rh2[i], ro[i], fh1[i], fh2[i], fo[i], energy);
344 };
345 };
346 assert(not std::isnan(energy) and not std::isinf(energy));
347}
348
349// --------------------------------------------- Square
350// ----------------------------------------
351
361void Ccl::ro_2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
362 double &d2) {
363 energy += ro_2_ * re2_ * (ro1._2 + ro2._2) / 2.0;
364 d1 += ro_2_ * re2_ * ro1._1;
365 d2 += ro_2_ * re2_ * ro2._1;
366}
367
368void Ccl::ro1_ro2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
369 double &d2) {
370 energy += ro1_ro2_ * re2_ * ro1._1 * ro2._1;
371 d1 += ro1_ro2_ * re2_ * ro2._1;
372 d2 += ro1_ro2_ * re2_ * ro1._1;
373}
374
375void Ccl::ro_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
376 double &energy, double &d1, double &d2, double &d3) {
377 energy += ro_theta_ * re2_ * (ro1._1 + ro2._1) * dth._1;
378 d1 += ro_theta_ * re2_ * dth._1; // dE/drho1
379 d2 += ro_theta_ * re2_ * dth._1; // dE/drho2
380 d3 += ro_theta_ * re2_ * (ro1._1 + ro2._1); // dE/d(theta)
381}
382
392void Ccl::theta_2(Dtheta const &dth, double &energy, double &d3) {
393 energy += theta_2_ * re2_ * dth._2 / 2.0;
394 d3 += theta_2_ * re2_ * dth._1; // dE/d(theta)= - moment
395}
396
397// --------------------------------------------- Cubic
398// ----------------------------------------
399
400void Ccl::ro_3(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
401 double &d2) {
402 energy += ro_3_ * re2_ * (ro1._3 + ro2._3);
403 d1 += 3.0 * ro_3_ * re2_ * ro1._2;
404 d2 += 3.0 * ro_3_ * re2_ * ro2._2;
405}
406
407void Ccl::ro_ro1_ro2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
408 double &d2) {
409 energy += ro_ro1_ro2_ * re2_ * (ro1._1 + ro2._1) * ro1._1 * ro2._1;
410 d1 += ro_ro1_ro2_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro2._2);
411 d2 += ro_ro1_ro2_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro1._2);
412}
413
414void Ccl::ro_2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
415 double &energy, double &d1, double &d2, double &d3) {
416 energy += ro_2_theta_ * re2_ * (ro1._2 + ro2._2) * dth._1;
417 d1 += ro_2_theta_ * re2_ * 2.0 * ro1._1 * dth._1; // dE/drho1
418 d2 += ro_2_theta_ * re2_ * 2.0 * ro2._1 * dth._1; // dE/drho2
419 d3 += ro_2_theta_ * re2_ * (ro1._2 + ro2._2); // dE/d(theta)
420}
421
422void Ccl::ro1_ro2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
423 double &energy, double &d1, double &d2, double &d3) {
424 energy += ro1_ro2_theta_ * re2_ * ro1._1 * ro2._1 * dth._1;
425 d1 += ro1_ro2_theta_ * re2_ * ro2._1 * dth._1; // dE/drho1
426 d2 += ro1_ro2_theta_ * re2_ * ro1._1 * dth._1; // dE/drho2
427 d3 += ro1_ro2_theta_ * re2_ * ro1._1 * ro2._1; // dE/d(theta)
428}
429
430void Ccl::ro_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
431 double &energy, double &d1, double &d2, double &d3) {
432 energy += ro_theta_2_ * re2_ * (ro1._1 + ro2._1) * dth._2;
433 d1 += ro_theta_2_ * re2_ * dth._2; // dE/drho1
434 d2 += ro_theta_2_ * re2_ * dth._2; // dE/drho2
435 d3 += ro_theta_2_ * re2_ * (ro1._1 + ro2._1) * 2.0 * dth._1; // dE/d(theta)
436}
437
438void Ccl::theta_3(Dtheta const &dth, double &energy, double &d3) {
439 energy += thetat_3_ * re2_ * dth._3;
440 d3 += thetat_3_ * re2_ * 3.0 * dth._2; // dE/d(theta)
441}
442
443// --------------------------------------------- Quartic
444// ----------------------------------------
445
446void Ccl::ro_4(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
447 double &d2) {
448 energy += ro_4_ * re2_ * (ro1._2 * ro1._2 + ro2._2 * ro2._2);
449 d1 += ro_4_ * re2_ * 4.0 * ro1._3; // dE/drho1
450 d2 += ro_4_ * re2_ * 4.0 * ro2._3; // dE/drho2
451}
452
453void Ccl::ro1_ro2_ro_2(Rho const &ro1, Rho const &ro2, double &energy,
454 double &d1, double &d2) {
455 energy += ro1_ro2_ro_2_ * re2_ * ro1._1 * ro2._1 * (ro1._2 + ro2._2);
456 d1 += ro1_ro2_ro_2_ * re2_ * (3.0 * ro2._1 * ro1._2 + ro2._3); // dE/drho1
457 d2 += ro1_ro2_ro_2_ * re2_ * (3.0 * ro1._1 * ro2._2 + ro1._3); // dE/drho2
458}
459
460void Ccl::ro1_2_ro2_2(Rho const &ro1, Rho const &ro2, double &energy,
461 double &d1, double &d2) {
462 energy += ro1_2_ro2_2_ * re2_ * ro1._2 * ro2._2;
463 d1 += ro1_2_ro2_2_ * re2_ * 2.0 * ro1._1 * ro2._2; // dE/drho1
464 d2 += ro1_2_ro2_2_ * re2_ * 2.0 * ro1._2 * ro2._1; // dE/drho2
465}
466
467void Ccl::ro_3_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
468 double &energy, double &d1, double &d2, double &d3) {
469 energy += ro_3_theta_ * re2_ * (ro1._3 + ro2._3) * dth._1;
470 d1 += ro_3_theta_ * re2_ * 3.0 * ro1._2 * dth._1; // dE/drho1
471 d2 += ro_3_theta_ * re2_ * 3.0 * ro2._2 * dth._1; // dE/drho2
472 d3 += ro_3_theta_ * re2_ * (ro1._3 + ro2._3); // dE/d(theta)
473}
474
475void Ccl::ro_ro1_ro2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
476 double &energy, double &d1, double &d2, double &d3) {
477 energy +=
478 ro_ro1_ro2_theta_ * re2_ * (ro1._1 + ro2._1) * ro1._1 * ro2._1 * dth._1;
479 d1 += ro_ro1_ro2_theta_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro2._2) *
480 dth._1; // dE/drho1
481 d2 += ro_ro1_ro2_theta_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro1._2) *
482 dth._1; // dE/drho2
483 d3 += ro_ro1_ro2_theta_ * re2_ * (ro1._1 + ro2._1) * ro1._1 *
484 ro2._1; // dE/d(theta)
485}
486
487void Ccl::ro_2_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
488 double &energy, double &d1, double &d2, double &d3) {
489 energy += ro_2_theta_2_ * re2_ * (ro1._2 + ro2._2) * dth._2;
490 d1 += ro_2_theta_2_ * re2_ * 2.0 * ro1._1 * dth._2; // dE/drho1
491 d2 += ro_2_theta_2_ * re2_ * 2.0 * ro2._1 * dth._2; // dE/drho2
492 d3 += ro_2_theta_2_ * re2_ * 2.0 * (ro1._2 + ro2._2) * dth._1; // dE/d(theta)
493}
494
495void Ccl::ro1_ro2_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
496 double &energy, double &d1, double &d2, double &d3) {
497 energy += ro1_ro2_theta_2_ * re2_ * ro1._1 * ro2._1 * dth._2;
498 d1 += ro1_ro2_theta_2_ * re2_ * ro2._1 * dth._2; // dE/drho1
499 d2 += ro1_ro2_theta_2_ * re2_ * ro1._1 * dth._2; // dE/drho2
500 d3 += ro1_ro2_theta_2_ * re2_ * ro1._1 * ro2._1 * 2.0 * dth._1; // dE/d(theta)
501}
502
503// rCcl::o_theta_3 = 0
504void Ccl::theta_4(Dtheta const &dth, double &energy, double &d3) {
505 energy += theta_4_ * re2_ * dth._2 * dth._2;
506 d3 += theta_4_ * re2_ * 4.0 * dth._3; // dE/d(theta)
507}
ccl.hpp
Potential CCL Table II for intramolecular interactions in water.
forcefields::Ccl::re_
static double const re_
Distance OH at equilibrium.
Definition ccl.hpp:56
forcefields::Ccl::thetae_
static double const thetae_
Distance OH at equilibrium.
Definition ccl.hpp:57
forcefields::Ccl::theta_4
void theta_4(Dtheta const &s3, double &energy, double &d3)
Definition ccl.cpp:504
forcefields::Ccl::ro_3
void ro_3(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:400
forcefields::Ccl::theta_3
void theta_3(Dtheta const &s3, double &energy, double &d3)
Definition ccl.cpp:438
forcefields::Ccl::ro_ro1_ro2_theta
void ro_ro1_ro2_theta(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:475
forcefields::Ccl::ro_theta
void ro_theta(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:375
forcefields::Ccl::ro_3_theta
void ro_3_theta(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:467
forcefields::Ccl::theta_2
void theta_2(Dtheta const &s3, double &energy, double &d3)
Energy and derivatives for Term 4.
Definition ccl.cpp:392
forcefields::Ccl::ro1_ro2_theta
void ro1_ro2_theta(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:422
forcefields::Ccl::ro1_ro2
void ro1_ro2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:368
forcefields::Ccl::getName
char const * getName() const
Name of the potential.
Definition ccl.cpp:166
forcefields::Ccl::initialiseRho
void initialiseRho(Vector3 const &v, Rho &r)
Initialise Rho.
Definition ccl.cpp:199
forcefields::Ccl::initialiseDtheta
void initialiseDtheta(Vector3 const &v1, Vector3 const &v2, Dtheta &dth, double const thetaEquilibrium)
Initialise Dtheta.
Definition ccl.cpp:209
forcefields::Ccl::intramolecular
void intramolecular(double const rh1[], double const rh2[], double const ro[], double fh1[], double fh2[], double fo[], double &energy)
Interactions inside one molecules.
Definition ccl.cpp:238
forcefields::Ccl::ro_2
void ro_2(Rho const &ro1, Rho const &ro2, double &energy, double &d1, double &d2)
Energy and derivatives for Term 1.
Definition ccl.cpp:361
forcefields::Ccl::ro_4
void ro_4(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:446
forcefields::Ccl::ro_2_theta_2
void ro_2_theta_2(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:487
forcefields::Ccl::ro_2_theta
void ro_2_theta(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:414
forcefields::Ccl::computeHH_O_
void computeHH_O_(const int nAtoms, const double R[], double F[], double &U, const double b[])
Compute the forces and the energy.
Definition ccl.cpp:130
forcefields::Ccl::ro1_ro2_ro_2
void ro1_ro2_ro_2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:453
forcefields::Ccl::computeTemplate
void computeTemplate(const int nMolecules, const double(*const rh1)[H *3], const double(*const rh2)[H *3], const double(*const ro)[O *3], double(*const fh1)[H *3], double(*const fh2)[H *3], double(*const fo)[O *3], double &energy, double const b[], bool const (*const xh1)[H]=0, bool const (*const xh2)[H]=0, bool const (*const xo)[O]=0)
Compute forces and energy.
Definition ccl.cpp:323
forcefields::Ccl::ro_theta_2
void ro_theta_2(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:430
forcefields::Ccl::ro1_ro2_theta_2
void ro1_ro2_theta_2(Rho const &s1, Rho const &s2, Dtheta const &s3, double &energy, double &d1, double &d2, double &d3)
Definition ccl.cpp:495
forcefields::Ccl::ro_ro1_ro2
void ro_ro1_ro2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:407
forcefields::Ccl::ro1_2_ro2_2
void ro1_2_ro2_2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:460
forcefields::PotentialBase::PotentialBase
PotentialBase()
Non bond interaction cutoff.
Definition potential_base.cpp:40
forcefields::PotentialBase::setPeriodicity
void setPeriodicity(const double periods[])
Set periodicity.
Definition potential_base.cpp:603
forcefields::PotentialBase::distance
void distance(const double x[], const double y[], double z[], double &z1, double &z2)
Distance vector, norm and norm square.
Definition potential_base.cpp:555
forcefields::PotentialBase::dotProduct
static double dotProduct(double const v[], double const w[])
Dot product.
Definition potential_base.hpp:216
forcefields::unit_system
Physical constants and unit conversion.
Definition tip4p_unit_system.hpp:105
forcefields::unit_system::DEGREE
const double DEGREE
Definition tip4p_unit_system.hpp:125
forcefields::unit_system::ANGSTROM
const double ANGSTROM
Definition tip4p_unit_system.hpp:265
forcefields::Ccl::Dtheta
parameter in paper.
Definition ccl.hpp:43
forcefields::Ccl::Dtheta::n2
double n2[3]
Same as n1 but for hydrogen 2.
Definition ccl.hpp:49
forcefields::Ccl::Dtheta::_2
double _2
(Delta theta)^2
Definition ccl.hpp:45
forcefields::Ccl::Dtheta::n1
double n1[3]
Convert derivative to force.
Definition ccl.hpp:47
forcefields::Ccl::Dtheta::_3
double _3
(Delta theta)^3
Definition ccl.hpp:46
forcefields::Ccl::Dtheta::_1
double _1
Delta theta.
Definition ccl.hpp:44
forcefields::Ccl::Rho
parameter in paper.
Definition ccl.hpp:36
forcefields::Ccl::Rho::_2
double _2
rho^2
Definition ccl.hpp:38
forcefields::Ccl::Rho::_3
double _3
rho^3
Definition ccl.hpp:39
forcefields::Ccl::Rho::n
double n[3]
Convert derivative to force.
Definition ccl.hpp:40
forcefields::Ccl::Rho::_1
double _1
rho^1
Definition ccl.hpp:37
forcefields::PotentialBase::Vector3::_2
double _2
Vector's norm square.
Definition potential_base.hpp:112