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
23namespace forcefields {
24namespace {
25#if defined(FORCEFIELDS_UNIT_SYSTEM_HPP) && \
26 (FORCEFIELDS_UNIT_SYSTEM_HPP != \
27 FORCEFIELDS_UNIT_SYSTEM_ELECTRONVOLT_ANGSTROM_FEMTOSECOND_ECHARGE)
28double const re_ = 0.9572 * unit_system::ANGSTROM;
29double const thetae_ = 104.52 * unit_system::DEGREE;
30
31double const re2_ = re_ * re_;
32
33// ------------------------ Quadratic ---------------------------
34double const ro_2_ = 84.54e-12 * unit_system::ERGS_PER_ANGSTROM2;
35double const ro1_ro2_ = -1.01e-12 * unit_system::ERGS_PER_ANGSTROM2;
36double const ro_theta_ = 2.288e-12 * unit_system::ERGS_PER_ANGSTROM2;
37double const theta_2_ = 7.607e-12 * unit_system::ERGS_PER_ANGSTROM2;
38
39// --------------------------------------------- Cubic
40// ----------------------------------------
41double const ro_3_ = -10.168e-12 * unit_system::ERGS_PER_ANGSTROM2;
42double const ro_ro1_ro2_ = 0.201e-12 * unit_system::ERGS_PER_ANGSTROM2;
43double const ro_2_theta_ = 4.308e-12 * unit_system::ERGS_PER_ANGSTROM2;
44double const ro1_ro2_theta_ = -4.020e-12 * unit_system::ERGS_PER_ANGSTROM2;
45double const ro_theta_2_ = -1.175e-12 * unit_system::ERGS_PER_ANGSTROM2;
46double const thetat_3_ = -1.595e-12 * unit_system::ERGS_PER_ANGSTROM2;
47
48// --------------------------------------------- Quartic
49// ----------------------------------------
50double const ro_4_ = -10.684e-12 * unit_system::ERGS_PER_ANGSTROM2;
51double const ro1_ro2_ro_2_ = -6.162e-12 * unit_system::ERGS_PER_ANGSTROM2;
52double const ro1_2_ro2_2_ = 2.717e-12 * unit_system::ERGS_PER_ANGSTROM2;
53
54double const ro_3_theta_ = 6.328e-12 * unit_system::ERGS_PER_ANGSTROM2;
55double const ro_ro1_ro2_theta_ = -4.020e-12 * unit_system::ERGS_PER_ANGSTROM2;
56
57double const ro_2_theta_2_ = -4.70e-12 * unit_system::ERGS_PER_ANGSTROM2;
58double const ro1_ro2_theta_2_ = 3.05e-12 * unit_system::ERGS_PER_ANGSTROM2;
59
60// ro_theta_3 = 0
61double const theta_4_ = -0.0318e-12 * unit_system::ERGS_PER_ANGSTROM2;
62#else
63double const re_ = 0.9572; // ANGSTROM
64double const thetae_ = 1.82421813418447321; // RADIANS
65
66double const re2_ = 0.91623184; // ANGSTROM^2
67
68// ------------------------ Quadratic ---------------------------
69double const ro_2_ = 52.7657211406036524; // eV / Angstrom^2
70double const ro1_ro2_ = -0.630392457440379528; // eV / Angstrom^2
71double const ro_theta_ = 1.42805736893424595; // eV / Angstrom^2
72double const theta_2_ = 4.74791626113759158; // eV / Angstrom^2
73
74// --------------------------------------------- Cubic
75// ----------------------------------------
76double const ro_3_ = -6.3463668388651282; // eV / Angstrom^2
77double const ro_ro1_ro2_ = 0.125454340540115145; // eV / Angstrom^2
78double const ro_2_theta_ = 2.68884228381500501; // eV / Angstrom^2
79double const ro1_ro2_theta_ = -2.509086810802303; // eV / Angstrom^2
80double const ro_theta_2_ = -0.733377363853906838; // eV / Angstrom^2
81double const thetat_3_ = -0.995520761997431114; // eV / Angstrom^2
82
83// --------------------------------------------- Quartic
84// ----------------------------------------
85double const ro_4_ = -6.668428728012886; // eV / Angstrom^2
86double const ro1_ro2_ro_2_ = -3.84601814133427622; // eV / Angstrom^2
87double const ro1_2_ro2_2_ = 1.69581812560941692; // eV / Angstrom^2
88double const ro_3_theta_ = 3.94962719869576429; // eV / Angstrom^2
89double const ro_ro1_ro2_theta_ = -2.509086810802303; // eV / Angstrom^2
90double const ro_2_theta_2_ = -2.93350945541562735; // eV / Angstrom^2
91double const ro1_ro2_theta_2_ = 1.90366039128035425; // eV / Angstrom^2
92double const theta_4_ = -0.0198480001451525438; // eV / Angstrom^2
93#endif
94} // namespace
95
106
108double const Ccl::re_ = re_;
110double const Ccl::thetae_ = thetae_;
111
125void Ccl::computeHH_O_(const int nAtoms, const double R[], double F[],
126 double &U, const double b[]) {
127 computeHH_O_(nAtoms, R, F, U, b, 0);
128}
129
143void Ccl::computeHH_O_(const int nAtoms, const double R[], double F[],
144 double &U, const double b[], const bool fixed[]) {
145 int const nMolecules = nAtoms / 3;
146 const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(R);
147 const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&R[3]);
148 const double (*const ro)[3] =
149 reinterpret_cast<const double (*)[3]>(&R[nMolecules * 6]);
150 double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(F);
151 double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&F[3]);
152 double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&F[nMolecules * 6]);
153 bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
154 bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
155 bool const(*const xo)[1] =
156 reinterpret_cast<bool const(*)[1]>(&fixed[nMolecules * 2]);
157 computeTemplate(nMolecules, rh1, rh2, ro, fh1, fh2, fo, U, b, xh1, xh2, xo);
158}
159
161char const *Ccl::getName() const { return "Ccl"; }
162
167Ccl::Ccl(double cutoff, double switchingWidth)
168 : PotentialBase(cutoff, switchingWidth) {}
169
177
192
194void Ccl::initialiseRho(Vector3 const &v, Rho &ro) {
195 ro._1 = (v._1 - re_) / v._1;
196 ro._2 = ro._1 * ro._1;
197 ro._3 = ro._2 * ro._1;
198 double const d = re_ / v._2 / v._1;
199 for (int i = 0; i < 3; ++i)
200 ro.n[i] = v.v[i] * d;
201}
202
204void Ccl::initialiseDtheta(Vector3 const &v1, Vector3 const &v2, Dtheta &dth,
205 double const thetaEquilibrium) {
206 double const cos_theta = dotProduct(v1.v, v2.v) / v1._1 / v2._1;
207 double const theta = std::acos(cos_theta);
208 dth._1 = theta - thetaEquilibrium;
209 dth._2 = dth._1 * dth._1;
210 dth._3 = dth._2 * dth._1;
211 double const d_theta = -1.0 / std::sqrt(1.0 - cos_theta * cos_theta);
212 for (int k = 0; k < 3; ++k) {
213 dth.n1[k] = 0.0;
214 dth.n2[k] = 0.0;
215 for (int j = 0; j < 3; ++j) {
216 dth.n1[k] -= v2.v[j] / v2._1 / v1._1 * v1.v[j] * v1.v[k] / v1._2;
217 dth.n2[k] -= v1.v[j] / v1._1 / v2._1 * v2.v[j] * v2.v[k] / v2._2;
218 };
219 dth.n1[k] += v2.v[k] / v2._1 / v1._1;
220 dth.n2[k] += v1.v[k] / v1._1 / v2._1;
221 dth.n1[k] *= d_theta;
222 dth.n2[k] *= d_theta;
223 };
224}
225
233void Ccl::intramolecular(double const rh1[], double const rh2[],
234 double const ro[], double fh1[], double fh2[],
235 double fo[], double &energy) {
236 Vector3 v1, v2;
237 distance(rh1, ro, v1);
238 distance(rh2, ro, v2);
239 // ------prepare _
240 Rho ro1, ro2;
241 initialiseRho(v1, ro1);
242 initialiseRho(v2, ro2);
243
244 // ------- prepare Delta theta
245 Dtheta dth;
246 initialiseDtheta(v1, v2, dth, thetae_);
247 intramolecular(ro1, ro2, dth, energy, fh1, fh2, fo);
248}
249
255void Ccl::intramolecular(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
256 double &energy, double fh1[], double fh2[],
257 double fo[]) {
258 double d1 = 0, d2 = 0,
259 d3 = 0; // d1=d(E)/d(rho1), d2=d(E)/d(rho2), d(E)/d(theta)
260
261 // -------- quadratic interaction
262 ro_2(ro1, ro2, energy, d1, d2);
263 ro1_ro2(ro1, ro2, energy, d1, d2);
264 ro_theta(ro1, ro2, dth, energy, d1, d2, d3);
265 theta_2(dth, energy, d3);
266 //*/
267
268 // -------- Cubic interaction
269 ro_3(ro1, ro2, energy, d1, d2);
270 ro_ro1_ro2(ro1, ro2, energy, d1, d2);
271 ro_2_theta(ro1, ro2, dth, energy, d1, d2, d3);
272 ro1_ro2_theta(ro1, ro2, dth, energy, d1, d2, d3);
273 ro_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
274 theta_3(dth, energy, d3);
275 //*/
276
277 // --------------------------------------------- Quartic
278 // -------------------------------------
279 ro_4(ro1, ro2, energy, d1, d2);
280 ro1_ro2_ro_2(ro1, ro2, energy, d1, d2);
281 ro1_2_ro2_2(ro1, ro2, energy, d1, d2);
282
283 ro_3_theta(ro1, ro2, dth, energy, d1, d2, d3);
284 ro_ro1_ro2_theta(ro1, ro2, dth, energy, d1, d2, d3);
285
286 ro_2_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
287 ro1_ro2_theta_2(ro1, ro2, dth, energy, d1, d2, d3);
288
289 // ro_theta_3 = 0
290 theta_4(dth, energy, d3);
291 //*/
292 for (int i = 0; i < 3; ++i) {
293 fh1[i] -= d1 * ro1.n[i] + d3 * dth.n1[i];
294 fh2[i] -= d2 * ro2.n[i] + d3 * dth.n2[i];
295 fo[i] += d1 * ro1.n[i] + d2 * ro2.n[i] + d3 * dth.n1[i] + d3 * dth.n2[i];
296 };
297}
298
317template <int H, int O>
318void Ccl::computeTemplate(
319 const int nMolecules, const double (*const rh1)[H * 3],
320 const double (*const rh2)[H * 3], const double (*const ro)[O * 3],
321 double (*const fh1)[H * 3], double (*const fh2)[H * 3],
322 double (*const fo)[O * 3], double &energy, double const b[],
323 bool const (*const xh1)[H], bool const (*const xh2)[H],
324 bool const (*const xo)[O]) {
325 for (int i = 0; i < nMolecules; ++i) {
326 for (int a = 0; a < 3; a++) {
327 fh1[i][a] = 0.0;
328 fh2[i][a] = 0.0;
329 fo[i][a] = 0.0;
330 };
331 };
332 energy = 0.0;
333 setPeriodicity(b);
334
335 for (int i = nMolecules - 1; i >= 0; --i) {
336 bool const areFixed = xh1[i][0] and xh2[i][0] and xo[i][0];
337 if (not areFixed) {
338 intramolecular(rh1[i], rh2[i], ro[i], fh1[i], fh2[i], fo[i], energy);
339 };
340 };
341 assert(not std::isnan(energy) and not std::isinf(energy));
342}
343
344// --------------------------------------------- Square
345// ----------------------------------------
346
356void Ccl::ro_2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
357 double &d2) {
358 energy += ro_2_ * re2_ * (ro1._2 + ro2._2) / 2.0;
359 d1 += ro_2_ * re2_ * ro1._1;
360 d2 += ro_2_ * re2_ * ro2._1;
361}
362
363void Ccl::ro1_ro2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
364 double &d2) {
365 energy += ro1_ro2_ * re2_ * ro1._1 * ro2._1;
366 d1 += ro1_ro2_ * re2_ * ro2._1;
367 d2 += ro1_ro2_ * re2_ * ro1._1;
368}
369
370void Ccl::ro_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
371 double &energy, double &d1, double &d2, double &d3) {
372 energy += ro_theta_ * re2_ * (ro1._1 + ro2._1) * dth._1;
373 d1 += ro_theta_ * re2_ * dth._1; // dE/drho1
374 d2 += ro_theta_ * re2_ * dth._1; // dE/drho2
375 d3 += ro_theta_ * re2_ * (ro1._1 + ro2._1); // dE/d(theta)
376}
377
387void Ccl::theta_2(Dtheta const &dth, double &energy, double &d3) {
388 energy += theta_2_ * re2_ * dth._2 / 2.0;
389 d3 += theta_2_ * re2_ * dth._1; // dE/d(theta)= - moment
390}
391
392// --------------------------------------------- Cubic
393// ----------------------------------------
394
395void Ccl::ro_3(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
396 double &d2) {
397 energy += ro_3_ * re2_ * (ro1._3 + ro2._3);
398 d1 += 3.0 * ro_3_ * re2_ * ro1._2;
399 d2 += 3.0 * ro_3_ * re2_ * ro2._2;
400}
401
402void Ccl::ro_ro1_ro2(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
403 double &d2) {
404 energy += ro_ro1_ro2_ * re2_ * (ro1._1 + ro2._1) * ro1._1 * ro2._1;
405 d1 += ro_ro1_ro2_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro2._2);
406 d2 += ro_ro1_ro2_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro1._2);
407}
408
409void Ccl::ro_2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
410 double &energy, double &d1, double &d2, double &d3) {
411 energy += ro_2_theta_ * re2_ * (ro1._2 + ro2._2) * dth._1;
412 d1 += ro_2_theta_ * re2_ * 2.0 * ro1._1 * dth._1; // dE/drho1
413 d2 += ro_2_theta_ * re2_ * 2.0 * ro2._1 * dth._1; // dE/drho2
414 d3 += ro_2_theta_ * re2_ * (ro1._2 + ro2._2); // dE/d(theta)
415}
416
417void Ccl::ro1_ro2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
418 double &energy, double &d1, double &d2, double &d3) {
419 energy += ro1_ro2_theta_ * re2_ * ro1._1 * ro2._1 * dth._1;
420 d1 += ro1_ro2_theta_ * re2_ * ro2._1 * dth._1; // dE/drho1
421 d2 += ro1_ro2_theta_ * re2_ * ro1._1 * dth._1; // dE/drho2
422 d3 += ro1_ro2_theta_ * re2_ * ro1._1 * ro2._1; // dE/d(theta)
423}
424
425void Ccl::ro_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
426 double &energy, double &d1, double &d2, double &d3) {
427 energy += ro_theta_2_ * re2_ * (ro1._1 + ro2._1) * dth._2;
428 d1 += ro_theta_2_ * re2_ * dth._2; // dE/drho1
429 d2 += ro_theta_2_ * re2_ * dth._2; // dE/drho2
430 d3 += ro_theta_2_ * re2_ * (ro1._1 + ro2._1) * 2.0 * dth._1; // dE/d(theta)
431}
432
433void Ccl::theta_3(Dtheta const &dth, double &energy, double &d3) {
434 energy += thetat_3_ * re2_ * dth._3;
435 d3 += thetat_3_ * re2_ * 3.0 * dth._2; // dE/d(theta)
436}
437
438// --------------------------------------------- Quartic
439// ----------------------------------------
440
441void Ccl::ro_4(Rho const &ro1, Rho const &ro2, double &energy, double &d1,
442 double &d2) {
443 energy += ro_4_ * re2_ * (ro1._2 * ro1._2 + ro2._2 * ro2._2);
444 d1 += ro_4_ * re2_ * 4.0 * ro1._3; // dE/drho1
445 d2 += ro_4_ * re2_ * 4.0 * ro2._3; // dE/drho2
446}
447
448void Ccl::ro1_ro2_ro_2(Rho const &ro1, Rho const &ro2, double &energy,
449 double &d1, double &d2) {
450 energy += ro1_ro2_ro_2_ * re2_ * ro1._1 * ro2._1 * (ro1._2 + ro2._2);
451 d1 += ro1_ro2_ro_2_ * re2_ * (3.0 * ro2._1 * ro1._2 + ro2._3); // dE/drho1
452 d2 += ro1_ro2_ro_2_ * re2_ * (3.0 * ro1._1 * ro2._2 + ro1._3); // dE/drho2
453}
454
455void Ccl::ro1_2_ro2_2(Rho const &ro1, Rho const &ro2, double &energy,
456 double &d1, double &d2) {
457 energy += ro1_2_ro2_2_ * re2_ * ro1._2 * ro2._2;
458 d1 += ro1_2_ro2_2_ * re2_ * 2.0 * ro1._1 * ro2._2; // dE/drho1
459 d2 += ro1_2_ro2_2_ * re2_ * 2.0 * ro1._2 * ro2._1; // dE/drho2
460}
461
462void Ccl::ro_3_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
463 double &energy, double &d1, double &d2, double &d3) {
464 energy += ro_3_theta_ * re2_ * (ro1._3 + ro2._3) * dth._1;
465 d1 += ro_3_theta_ * re2_ * 3.0 * ro1._2 * dth._1; // dE/drho1
466 d2 += ro_3_theta_ * re2_ * 3.0 * ro2._2 * dth._1; // dE/drho2
467 d3 += ro_3_theta_ * re2_ * (ro1._3 + ro2._3); // dE/d(theta)
468}
469
470void Ccl::ro_ro1_ro2_theta(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
471 double &energy, double &d1, double &d2, double &d3) {
472 energy +=
473 ro_ro1_ro2_theta_ * re2_ * (ro1._1 + ro2._1) * ro1._1 * ro2._1 * dth._1;
474 d1 += ro_ro1_ro2_theta_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro2._2) *
475 dth._1; // dE/drho1
476 d2 += ro_ro1_ro2_theta_ * re2_ * (2.0 * ro1._1 * ro2._1 + ro1._2) *
477 dth._1; // dE/drho2
478 d3 += ro_ro1_ro2_theta_ * re2_ * (ro1._1 + ro2._1) * ro1._1 *
479 ro2._1; // dE/d(theta)
480}
481
482void Ccl::ro_2_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
483 double &energy, double &d1, double &d2, double &d3) {
484 energy += ro_2_theta_2_ * re2_ * (ro1._2 + ro2._2) * dth._2;
485 d1 += ro_2_theta_2_ * re2_ * 2.0 * ro1._1 * dth._2; // dE/drho1
486 d2 += ro_2_theta_2_ * re2_ * 2.0 * ro2._1 * dth._2; // dE/drho2
487 d3 += ro_2_theta_2_ * re2_ * 2.0 * (ro1._2 + ro2._2) * dth._1; // dE/d(theta)
488}
489
490void Ccl::ro1_ro2_theta_2(Rho const &ro1, Rho const &ro2, Dtheta const &dth,
491 double &energy, double &d1, double &d2, double &d3) {
492 energy += ro1_ro2_theta_2_ * re2_ * ro1._1 * ro2._1 * dth._2;
493 d1 += ro1_ro2_theta_2_ * re2_ * ro2._1 * dth._2; // dE/drho1
494 d2 += ro1_ro2_theta_2_ * re2_ * ro1._1 * dth._2; // dE/drho2
495 d3 += ro1_ro2_theta_2_ * re2_ * ro1._1 * ro2._1 * 2.0 * dth._1; // dE/d(theta)
496}
497
498// rCcl::o_theta_3 = 0
499void Ccl::theta_4(Dtheta const &dth, double &energy, double &d3) {
500 energy += theta_4_ * re2_ * dth._2 * dth._2;
501 d3 += theta_4_ * re2_ * 4.0 * dth._3; // dE/d(theta)
502}
503} // namespace forcefields
ccl.hpp
Potential CCL Table II for intramolecular interactions in water.
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:470
forcefields::Ccl::ro_3
void ro_3(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:395
forcefields::Ccl::initialiseRho
void initialiseRho(Vector3 const &v, Rho &r)
Initialise Rho.
Definition ccl.cpp:194
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:462
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:482
forcefields::Ccl::theta_4
void theta_4(Dtheta const &s3, double &energy, double &d3)
Definition ccl.cpp:499
forcefields::Ccl::thetae_
static double const thetae_
Distance OH at equilibrium.
Definition ccl.hpp:57
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:233
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:370
forcefields::Ccl::initialiseDtheta
void initialiseDtheta(Vector3 const &v1, Vector3 const &v2, Dtheta &dth, double const thetaEquilibrium)
Initialise Dtheta.
Definition ccl.cpp:204
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:417
forcefields::Ccl::ro_ro1_ro2
void ro_ro1_ro2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:402
forcefields::Ccl::getName
char const * getName() const
Name of the potential.
Definition ccl.cpp:161
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:318
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:356
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:448
forcefields::Ccl::re_
static double const re_
Distance OH at equilibrium.
Definition ccl.hpp:56
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:409
forcefields::Ccl::ro1_ro2
void ro1_ro2(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:363
forcefields::Ccl::theta_2
void theta_2(Dtheta const &s3, double &energy, double &d3)
Energy and derivatives for Term 4.
Definition ccl.cpp:387
forcefields::Ccl::theta_3
void theta_3(Dtheta const &s3, double &energy, double &d3)
Definition ccl.cpp:433
forcefields::Ccl::ro_4
void ro_4(Rho const &s1, Rho const &s2, double &energy, double &d1, double &d2)
Definition ccl.cpp:441
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:490
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:455
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:125
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:425
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::PotentialBase
PotentialBase()
Non bond interaction cutoff.
Definition potential_base.cpp:40
forcefields::PotentialBase::dotProduct
static double dotProduct(double const v[], double const w[])
Dot product.
Definition potential_base.hpp:216
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