Loading...
Searching...
No Matches
eonc::pairhess Namespace Reference

Functions

void addPairBlock (Eigen::MatrixXd &P, int i, int j, const Eigen::Vector3d &rij, double k_par, double k_perp)
void ljDerivs (double r, double eps, double sigma, double &vp, double &vpp)
void morseDerivs (double r, double D, double a, double re, double &vp, double &vpp)
bool isModelHess (const std::string &kind)
void addScalarInternal (Eigen::MatrixXd &P, const int *atoms, const Eigen::Vector3d *g, int n_atoms, double k)
 Rank-1 update (k\,(\nabla q)(\nabla q)^\top) for a scalar internal.
bool angleGrads (const Eigen::Vector3d &ri, const Eigen::Vector3d &rj, const Eigen::Vector3d &rk, Eigen::Vector3d g[3])
 Wilson (B) rows for the valence angle (i{-}j{-}k) ((j) central).
bool torsionGrads (const Eigen::Vector3d &ri, const Eigen::Vector3d &rj, const Eigen::Vector3d &rk, const Eigen::Vector3d &rl, Eigen::Vector3d g[4])
 Wilson (B) rows for the torsion (i{-}j{-}k{-}l).
double lindhRho (double r, double r_nn, double alpha)
double swartScreen (double r, double r_nn)
void addModelHess (Eigen::MatrixXd &P, const Eigen::VectorXd &pos, const std::string &kind, double A, double mu, double r_nn, double rcut, ObjectiveFunction &objf)
Eigen::MatrixXd build (const Eigen::VectorXd &pos, const std::string &kind, PotType pot, double A, double mu, double rcut_in, ObjectiveFunction &objf)
 Analytic pair or model Hessian.

Function Documentation

◆ addModelHess()

void eonc::pairhess::addModelHess ( Eigen::MatrixXd & P,
const Eigen::VectorXd & pos,
const std::string & kind,
double A,
double mu,
double r_nn,
double rcut,
ObjectiveFunction & objf )
inline

Definition at line 150 of file PairHessian.h.

152 {
153 const int nat = static_cast<int>(pos.size()) / 3;
154 const double alpha = A > 0.0 ? A : 1.0;
155 const double scale = mu > 0.0 ? mu : 1.0;
156 const double r_bond = std::min(rcut, 1.35 * r_nn);
157 const double r_cov = 0.5 * r_nn;
158 auto mic = [&](int i, int j) {
159 Eigen::Vector3d dr = pos.segment<3>(3 * i) - pos.segment<3>(3 * j);
160 objf.minimumImage(dr);
161 return dr;
162 };
163 std::vector<std::vector<int>> nb(static_cast<size_t>(nat));
164 const double r_pair = (kind == "swart") ? rcut : r_bond;
165 for (int i = 0; i < nat; ++i) {
166 for (int j = i + 1; j < nat; ++j) {
167 const Eigen::Vector3d rij = mic(i, j);
168 const double r = rij.norm();
169 if (r >= r_pair || r < 1.0e-14) {
170 continue;
171 }
172 double k = 0.0;
173 if (kind == "schlegel") {
174 const double den = std::max(r - 0.45 * r_nn, 0.15);
175 k = scale * 1.734 / (den * den * den);
176 } else if (kind == "fischer") {
177 k = scale * 0.3601 * std::exp(-1.944 * (r - r_cov));
178 } else {
179 k = scale * lindhRho(r, r_nn, alpha);
180 if (kind == "swart") {
181 k *= swartScreen(r, r_nn);
182 }
183 }
184 addPairBlock(P, i, j, rij, k, 0.0);
185 if (r < r_bond) {
186 nb[static_cast<size_t>(i)].push_back(j);
187 nb[static_cast<size_t>(j)].push_back(i);
188 }
189 }
190 }
191 for (int j = 0; j < nat; ++j) {
192 const auto &nbr = nb[static_cast<size_t>(j)];
193 for (size_t a = 0; a < nbr.size(); ++a) {
194 for (size_t b = a + 1; b < nbr.size(); ++b) {
195 const int i = nbr[a];
196 const int k = nbr[b];
197 const Eigen::Vector3d ri = pos.segment<3>(3 * i);
198 const Eigen::Vector3d rj = pos.segment<3>(3 * j);
199 const Eigen::Vector3d rk = pos.segment<3>(3 * k);
200 Eigen::Vector3d g[3];
201 Eigen::Vector3d rji = ri - rj;
202 Eigen::Vector3d rjk = rk - rj;
203 objf.minimumImage(rji);
204 objf.minimumImage(rjk);
205 if (!angleGrads(rj + rji, rj, rj + rjk, g)) {
206 continue;
207 }
208 const double rij = rji.norm();
209 const double rkj = rjk.norm();
210 double kang = 0.16 * scale;
211 if (kind == "fischer") {
212 kang = scale * (0.089 + 0.11 *
213 std::pow(r_cov * r_cov, -0.42) *
214 std::exp(-0.44 * ((rij - r_cov) +
215 (rkj - r_cov))));
216 } else if (kind == "lindh_full" || kind == "swart") {
217 kang = 0.15 * scale * lindhRho(rij, r_nn, alpha) *
218 lindhRho(rkj, r_nn, alpha);
219 if (kind == "swart") {
220 kang *= swartScreen(rij, r_nn) * swartScreen(rkj, r_nn);
221 }
222 }
223 const int atoms[3] = {i, j, k};
224 addScalarInternal(P, atoms, g, 3, kang);
225 }
226 }
227 }
228 for (int j = 0; j < nat; ++j) {
229 const auto &nj = nb[static_cast<size_t>(j)];
230 for (int k : nj) {
231 if (k <= j) {
232 continue;
233 }
234 const auto &nk = nb[static_cast<size_t>(k)];
235 for (int i : nj) {
236 if (i == k) {
237 continue;
238 }
239 for (int l : nk) {
240 if (l == j || l == i) {
241 continue;
242 }
243 Eigen::Vector3d rijv = pos.segment<3>(3 * i) - pos.segment<3>(3 * j);
244 Eigen::Vector3d rkjv = pos.segment<3>(3 * k) - pos.segment<3>(3 * j);
245 Eigen::Vector3d rklv = pos.segment<3>(3 * k) - pos.segment<3>(3 * l);
246 objf.minimumImage(rijv);
247 objf.minimumImage(rkjv);
248 objf.minimumImage(rklv);
249 const Eigen::Vector3d rj = pos.segment<3>(3 * j);
250 Eigen::Vector3d g[4];
251 if (!torsionGrads(rj + rijv, rj, rj + rkjv, rj + rkjv - rklv, g)) {
252 continue;
253 }
254 const double rij = rijv.norm();
255 const double rjk = rkjv.norm();
256 const double rkl = rklv.norm();
257 double ktor = 0.01 * scale;
258 if (kind == "fischer") {
259 ktor = 0.0015 * scale;
260 } else if (kind == "lindh_full" || kind == "swart") {
261 ktor = 0.005 * scale * lindhRho(rij, r_nn, alpha) *
262 lindhRho(rjk, r_nn, alpha) * lindhRho(rkl, r_nn, alpha);
263 if (kind == "swart") {
264 ktor *= swartScreen(rij, r_nn) * swartScreen(rjk, r_nn) *
265 swartScreen(rkl, r_nn);
266 }
267 }
268 const int atoms[4] = {i, j, k, l};
269 addScalarInternal(P, atoms, g, 4, ktor);
270 }
271 }
272 }
273 }
274}
eonc::ObjectiveFunction::minimumImage
virtual void minimumImage(Eigen::Ref< Eigen::Vector3d >) const
Definition ObjectiveFunction.h:43
eonc::pairhess::angleGrads
bool angleGrads(const Eigen::Vector3d &ri, const Eigen::Vector3d &rj, const Eigen::Vector3d &rk, Eigen::Vector3d g[3])
Wilson (B) rows for the valence angle (i{-}j{-}k) ((j) central).
Definition PairHessian.h:95
eonc::pairhess::lindhRho
double lindhRho(double r, double r_nn, double alpha)
Definition PairHessian.h:142
eonc::pairhess::addPairBlock
void addPairBlock(Eigen::MatrixXd &P, int i, int j, const Eigen::Vector3d &rij, double k_par, double k_perp)
Definition PairHessian.h:27
eonc::pairhess::torsionGrads
bool torsionGrads(const Eigen::Vector3d &ri, const Eigen::Vector3d &rj, const Eigen::Vector3d &rk, const Eigen::Vector3d &rl, Eigen::Vector3d g[4])
Wilson (B) rows for the torsion (i{-}j{-}k{-}l).
Definition PairHessian.h:119
eonc::pairhess::addScalarInternal
void addScalarInternal(Eigen::MatrixXd &P, const int *atoms, const Eigen::Vector3d *g, int n_atoms, double k)
Rank-1 update (k\,(\nabla q)(\nabla q)^\top) for a scalar internal.
Definition PairHessian.h:77
eonc::pairhess::swartScreen
double swartScreen(double r, double r_nn)
Definition PairHessian.h:146

◆ addPairBlock()

void eonc::pairhess::addPairBlock ( Eigen::MatrixXd & P,
int i,
int j,
const Eigen::Vector3d & rij,
double k_par,
double k_perp )
inline

Definition at line 27 of file PairHessian.h.

29 {
30 const double r = rij.norm();
31 if (r < 1.0e-14 || (!std::isfinite(k_par) && !std::isfinite(k_perp))) {
32 return;
33 }
34 if (k_par == 0.0 && k_perp == 0.0) {
35 return;
36 }
37 const Eigen::Vector3d u = rij / r;
38 const Eigen::Matrix3d H = k_perp * Eigen::Matrix3d::Identity() +
39 (k_par - k_perp) * (u * u.transpose());
40 for (int a = 0; a < 3; ++a) {
41 for (int b = 0; b < 3; ++b) {
42 const int ia = 3 * i + a;
43 const int ib = 3 * i + b;
44 const int ja = 3 * j + a;
45 const int jb = 3 * j + b;
46 P(ia, ib) += H(a, b);
47 P(ja, jb) += H(a, b);
48 P(ia, jb) -= H(a, b);
49 P(ja, ib) -= H(a, b);
50 }
51 }
52}

◆ addScalarInternal()

void eonc::pairhess::addScalarInternal ( Eigen::MatrixXd & P,
const int * atoms,
const Eigen::Vector3d * g,
int n_atoms,
double k )
inline

Rank-1 update (k\,(\nabla q)(\nabla q)^\top) for a scalar internal.

Definition at line 77 of file PairHessian.h.

79 {
80 if (!(k > 0.0) || !std::isfinite(k)) {
81 return;
82 }
83 for (int a = 0; a < n_atoms; ++a) {
84 for (int b = 0; b < n_atoms; ++b) {
85 for (int x = 0; x < 3; ++x) {
86 for (int y = 0; y < 3; ++y) {
87 P(3 * atoms[a] + x, 3 * atoms[b] + y) += k * g[a][x] * g[b][y];
88 }
89 }
90 }
91 }
92}

◆ angleGrads()

bool eonc::pairhess::angleGrads ( const Eigen::Vector3d & ri,
const Eigen::Vector3d & rj,
const Eigen::Vector3d & rk,
Eigen::Vector3d g[3] )
inline

Wilson (B) rows for the valence angle (i{-}j{-}k) ((j) central).

Definition at line 95 of file PairHessian.h.

96 {
97 Eigen::Vector3d u = ri - rj;
98 Eigen::Vector3d v = rk - rj;
99 const double pu = u.norm();
100 const double pv = v.norm();
101 if (pu < 1.0e-14 || pv < 1.0e-14) {
102 return false;
103 }
104 u /= pu;
105 v /= pv;
106 double c = u.dot(v);
107 c = std::max(-1.0 + 1.0e-14, std::min(1.0 - 1.0e-14, c));
108 const double s = std::sqrt(std::max(0.0, 1.0 - c * c));
109 if (s < 1.0e-8) {
110 return false;
111 }
112 g[0] = (c * u - v) / (pu * s);
113 g[2] = (c * v - u) / (pv * s);
114 g[1] = -g[0] - g[2];
115 return true;
116}

◆ build()

Eigen::MatrixXd eonc::pairhess::build ( const Eigen::VectorXd & pos,
const std::string & kind,
PotType pot,
double A,
double mu,
double rcut_in,
ObjectiveFunction & objf )
inline

Analytic pair or model Hessian.

kind is pair / pair_abs / pair_full / lindh / lindh_full / exp / c1 / fischer / schlegel / swart.

Definition at line 279 of file PairHessian.h.

280 {
281 const int n = static_cast<int>(pos.size());
282 const int nat = n / 3;
283 std::vector<double> nn(static_cast<size_t>(nat),
284 std::numeric_limits<double>::infinity());
285 auto mic = [&](int i, int j) {
286 Eigen::Vector3d dr = pos.segment<3>(3 * i) - pos.segment<3>(3 * j);
287 objf.minimumImage(dr);
288 return dr;
289 };
290 for (int i = 0; i < nat; ++i) {
291 for (int j = 0; j < nat; ++j) {
292 if (i == j) {
293 continue;
294 }
295 const double r = mic(i, j).norm();
296 if (r > 1.0e-12 && r < nn[static_cast<size_t>(i)]) {
297 nn[static_cast<size_t>(i)] = r;
298 }
299 }
300 }
301 double r_nn = 0.0;
302 for (double v : nn) {
303 if (std::isfinite(v)) {
304 r_nn = std::max(r_nn, v);
305 }
306 }
307 if (r_nn <= 0.0) {
308 r_nn = 1.0;
309 }
310 const double rcut = rcut_in > 0.0 ? rcut_in : 2.5 * r_nn;
311 Eigen::MatrixXd P = Eigen::MatrixXd::Zero(n, n);
312 if (isModelHess(kind)) {
313 addModelHess(P, pos, kind, A, mu, r_nn, rcut, objf);
314 const double shift = 1.0e-6 * std::max(1.0, P.diagonal().mean());
315 P.diagonal().array() += shift;
316 return P;
317 }
318 for (int i = 0; i < nat; ++i) {
319 for (int j = i + 1; j < nat; ++j) {
320 const Eigen::Vector3d rij = mic(i, j);
321 const double r = rij.norm();
322 if (r >= rcut || r < 1.0e-14) {
323 continue;
324 }
325 if (kind == "pair" || kind == "pair_abs" || kind == "pair_full") {
326 double vp = 0.0;
327 double vpp = 0.0;
328 if (pot == PotType::MORSE_PT) {
329 morseDerivs(r, 0.7102, 1.6047, 2.897, vp, vpp);
330 } else {
331 ljDerivs(r, 1.0, 1.0, vp, vpp);
332 }
333 double k_par = std::max(vpp, 0.0);
334 double k_perp = std::max(vp / r, 0.0);
335 if (kind == "pair_abs") {
336 k_par = std::abs(vpp);
337 k_perp = std::abs(vp / r);
338 } else if (kind == "pair_full") {
339 k_par = vpp;
340 k_perp = vp / r;
341 }
342 addPairBlock(P, i, j, rij, k_par, k_perp);
343 } else if (kind == "lindh") {
344 const double alpha = A > 0.0 ? A : 1.0;
345 const double k = mu * std::exp(alpha * (r_nn * r_nn - r * r));
346 addPairBlock(P, i, j, rij, k, 0.0);
347 } else if (kind == "c1") {
348 const double x = r / rcut;
349 const double w = mu * (1.0 - x) * (1.0 - x) * (1.0 + 2.0 * x);
350 addPairBlock(P, i, j, rij, w, w);
351 } else {
352 const double w = mu * std::exp(-A * (r / r_nn - 1.0));
353 addPairBlock(P, i, j, rij, w, w);
354 }
355 }
356 }
357 const double shift = 1.0e-6 * std::max(1.0, P.diagonal().mean());
358 P.diagonal().array() += shift;
359 return P;
360}
eonc::pairhess::addModelHess
void addModelHess(Eigen::MatrixXd &P, const Eigen::VectorXd &pos, const std::string &kind, double A, double mu, double r_nn, double rcut, ObjectiveFunction &objf)
Definition PairHessian.h:150
eonc::pairhess::ljDerivs
void ljDerivs(double r, double eps, double sigma, double &vp, double &vpp)
Definition PairHessian.h:54
eonc::pairhess::morseDerivs
void morseDerivs(double r, double D, double a, double re, double &vp, double &vpp)
Definition PairHessian.h:63
eonc::pairhess::isModelHess
bool isModelHess(const std::string &kind)
Definition PairHessian.h:71

◆ isModelHess()

bool eonc::pairhess::isModelHess ( const std::string & kind)
inline

Definition at line 71 of file PairHessian.h.

71 {
72 return kind == "lindh_full" || kind == "swart" || kind == "fischer" ||
73 kind == "schlegel";
74}

◆ lindhRho()

double eonc::pairhess::lindhRho ( double r,
double r_nn,
double alpha )
inline

Definition at line 142 of file PairHessian.h.

142 {
143 return std::exp(alpha * (r_nn * r_nn - r * r));
144}

◆ ljDerivs()

void eonc::pairhess::ljDerivs ( double r,
double eps,
double sigma,
double & vp,
double & vpp )
inline

Definition at line 54 of file PairHessian.h.

55 {
56 const double inv = sigma / r;
57 const double x6 = inv * inv * inv * inv * inv * inv;
58 const double x12 = x6 * x6;
59 vp = 4.0 * eps * (-12.0 * x12 + 6.0 * x6) / r;
60 vpp = 4.0 * eps * (156.0 * x12 - 42.0 * x6) / (r * r);
61}

◆ morseDerivs()

void eonc::pairhess::morseDerivs ( double r,
double D,
double a,
double re,
double & vp,
double & vpp )
inline

Definition at line 63 of file PairHessian.h.

64 {
65 const double e1 = std::exp(-a * (r - re));
66 const double e2 = e1 * e1;
67 vp = 2.0 * D * a * (e1 - e2);
68 vpp = 2.0 * D * a * a * (2.0 * e2 - e1);
69}

◆ swartScreen()

double eonc::pairhess::swartScreen ( double r,
double r_nn )
inline

Definition at line 146 of file PairHessian.h.

146 {
147 return 1.0 / (1.0 + std::exp(12.0 * (r / r_nn - 1.2)));
148}

◆ torsionGrads()

bool eonc::pairhess::torsionGrads ( const Eigen::Vector3d & ri,
const Eigen::Vector3d & rj,
const Eigen::Vector3d & rk,
const Eigen::Vector3d & rl,
Eigen::Vector3d g[4] )
inline

Wilson (B) rows for the torsion (i{-}j{-}k{-}l).

Definition at line 119 of file PairHessian.h.

121 {
122 const Eigen::Vector3d vij = ri - rj;
123 const Eigen::Vector3d vkj = rk - rj;
124 const Eigen::Vector3d vkl = rk - rl;
125 Eigen::Vector3d n1 = vij.cross(vkj);
126 Eigen::Vector3d n2 = vkl.cross(vkj);
127 const double n1n = n1.norm();
128 const double n2n = n2.norm();
129 const double b2 = vkj.norm();
130 if (n1n < 1.0e-14 || n2n < 1.0e-14 || b2 < 1.0e-14) {
131 return false;
132 }
133 g[0] = -(b2 / (n1n * n1n)) * n1;
134 g[3] = (b2 / (n2n * n2n)) * n2;
135 const double c1 = vij.dot(vkj) / (b2 * b2);
136 const double c2 = vkl.dot(vkj) / (b2 * b2);
137 g[1] = -g[0] + c1 * g[0] - c2 * g[3];
138 g[2] = -g[3] - c1 * g[0] + c2 * g[3];
139 return true;
140}