1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
#include "CbPotential.h"
#include <cmath>
namespace ATC
{
static const double EPS = 1.0e-8;
// Approximates the derivative of phi
double CbPotential::phi_r(const double &r) const
{
const double dr = r*EPS;
return (phi(r+dr)-phi(r)) / dr;
}
// Approximates the second derivative of phi
double CbPotential::phi_rr(const double &r) const
{
const double dr = r*EPS;
return (phi_r(r+dr)-phi_r(r)) / dr;
}
// Approximates the third derivative of phi
double CbPotential::phi_rrr(const double &r) const
{
const double dr = r*EPS;
return (phi_rr(r+dr)-phi_rr(r)) / dr;
}
// Approximates the derivative of rho
double CbPotential::rho_r(const double &r) const
{
const double dr = r*EPS;
return (rho(r+dr)-rho(r)) / dr;
}
// Approximates the second derivative of rho
double CbPotential::rho_rr(const double &r) const
{
const double dr = r*EPS;
return (rho_r(r+dr)-rho_r(r)) / dr;
}
// Approximates the third derivative of rho
double CbPotential::rho_rrr(const double &r) const
{
const double dr = r*EPS;
return (rho_rr(r+dr)-rho_rr(r)) / dr;
}
// Approximates the derivative of the embedding function
double CbPotential::F_p(const double &p) const
{
const double dp = p*EPS;
return (F(p+dp)-F(p)) / dp;
}
// Approximates the second derivative of the embedding function
double CbPotential::F_pp(const double &p) const
{
const double dp = p*EPS;
return (F_p(p+dp)-F_p(p)) / dp;
}
// Approximates the third derivative of the embedding function
double CbPotential::F_ppp(const double &p) const
{
const double dp = p*EPS;
return (F_pp(p+dp)-F_pp(p)) / dp;
}
// Approximates the derivative of phi3.
double CbPotential::phi3_q (const double &q) const
{
const double dq = q*EPS;
return (phi3(q+dq)-phi3(q)) / dq;
}
// Approximates the second derivative of phi3.
double CbPotential::phi3_qq(const double &q) const
{
const double dq = q*EPS;
return (phi3_q(q+dq)-phi3_q(q)) / dq;
}
// Compute bond angle jik from the squared length of vectors ij,ik,kj.
double calculate_theta(double ij2, double ik2, double jk2)
{
return acos( 0.5*(ik2+ij2-jk2)/sqrt(ij2*ik2) );
}
// Initializes atomic interactions for up to three different terms.
Interactions::Interactions(int a, int b, int c)
{
// bitwise OR combines the terms that are listed
const int abc = a|b|c;
pairwise = (abc&PAIRWISE)>0;
embedding = (abc&EAM)>0;
three_body = (abc&THREE_BDY)>0;
angle_bending = (abc&ANGLE_BND)>0;
}
}
|
