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 package ffx.potential.bonded;
39
40 import ffx.numerics.atomic.AtomicDoubleArray3D;
41 import ffx.potential.parameters.TorsionType;
42
43 import java.io.Serial;
44 import java.util.Arrays;
45 import java.util.function.DoubleUnaryOperator;
46
47 import static java.lang.String.format;
48 import static org.apache.commons.math3.util.FastMath.acos;
49 import static org.apache.commons.math3.util.FastMath.sqrt;
50 import static org.apache.commons.math3.util.FastMath.toDegrees;
51
52
53
54
55
56
57
58 public class RestrainTorsion extends BondedTerm implements LambdaInterface {
59
60 @Serial
61 private static final long serialVersionUID = 1L;
62
63 private final Atom[] atoms;
64 public final TorsionType torsionType;
65 private final boolean lambdaTerm;
66 private final DoubleUnaryOperator lamMapper;
67 public final double units;
68
69 private double lambda = 1.0;
70 private double dEdL = 0.0;
71 private double d2EdL2 = 0.0;
72
73
74
75
76
77
78
79
80
81
82
83
84
85 public RestrainTorsion(Atom a1, Atom a2, Atom a3, Atom a4, TorsionType torsionType,
86 boolean lambdaEnabled, boolean reverseLambda, double units) {
87 atoms = new Atom[]{a1, a2, a3, a4};
88 this.torsionType = torsionType;
89 lambdaTerm = lambdaEnabled;
90 if (this.lambdaTerm) {
91 if (reverseLambda) {
92 lamMapper = (double l) -> 1.0 - l;
93 } else {
94 lamMapper = (double l) -> l;
95 }
96 } else {
97 lamMapper = (double l) -> 1.0;
98 }
99 this.units = units;
100 }
101
102 @Override
103 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
104 energy = 0.0;
105 value = 0.0;
106 dEdL = 0.0;
107
108 if (!getUse()) {
109 return energy;
110 }
111 var atomA = atoms[0];
112 var atomB = atoms[1];
113 var atomC = atoms[2];
114 var atomD = atoms[3];
115 var va = atomA.getXYZ();
116 var vb = atomB.getXYZ();
117 var vc = atomC.getXYZ();
118 var vd = atomD.getXYZ();
119 var vba = vb.sub(va);
120 var vcb = vc.sub(vb);
121 var vdc = vd.sub(vc);
122 var vt = vba.X(vcb);
123 var vu = vcb.X(vdc);
124 var rt2 = vt.length2();
125 var ru2 = vu.length2();
126 var rtru2 = rt2 * ru2;
127 if (rtru2 != 0.0) {
128 var rr = sqrt(rtru2);
129 var rcb = vcb.length();
130 var cosine = vt.dot(vu) / rr;
131 var sine = vcb.dot(vt.X(vu)) / (rcb * rr);
132 value = toDegrees(acos(cosine));
133 if (sine < 0.0) {
134 value = -value;
135 }
136 var amp = torsionType.amplitude;
137 var tsin = torsionType.sine;
138 var tcos = torsionType.cosine;
139 energy = amp[0] * (1.0 + cosine * tcos[0] + sine * tsin[0]);
140 var dedphi = amp[0] * (cosine * tsin[0] - sine * tcos[0]);
141 var cosprev = cosine;
142 var sinprev = sine;
143 var n = torsionType.terms;
144 for (int i = 1; i < n; i++) {
145 var cosn = cosine * cosprev - sine * sinprev;
146 var sinn = sine * cosprev + cosine * sinprev;
147 var phi = 1.0 + cosn * tcos[i] + sinn * tsin[i];
148 var dphi = (1.0 + i) * (cosn * tsin[i] - sinn * tcos[i]);
149 energy = energy + amp[i] * phi;
150 dedphi = dedphi + amp[i] * dphi;
151 cosprev = cosn;
152 sinprev = sinn;
153 }
154 energy = units * energy * lambda;
155 dEdL = units * energy;
156 if (gradient || lambdaTerm) {
157 dedphi = units * dedphi;
158 var vca = vc.sub(va);
159 var vdb = vd.sub(vb);
160 var dedt = vt.X(vcb).scaleI(dedphi / (rt2 * rcb));
161 var dedu = vu.X(vcb).scaleI(-dedphi / (ru2 * rcb));
162 var ga = dedt.X(vcb);
163 var gb = vca.X(dedt).addI(dedu.X(vdc));
164 var gc = dedt.X(vba).addI(vdb.X(dedu));
165 var gd = dedu.X(vcb);
166 int ia = atomA.getIndex() - 1;
167 int ib = atomB.getIndex() - 1;
168 int ic = atomC.getIndex() - 1;
169 int id = atomD.getIndex() - 1;
170 if (lambdaTerm) {
171 lambdaGrad.add(threadID, ia, ga);
172 lambdaGrad.add(threadID, ib, gb);
173 lambdaGrad.add(threadID, ic, gc);
174 lambdaGrad.add(threadID, id, gd);
175 }
176 if (gradient) {
177 grad.add(threadID, ia, ga.scaleI(lambda));
178 grad.add(threadID, ib, gb.scaleI(lambda));
179 grad.add(threadID, ic, gc.scaleI(lambda));
180 grad.add(threadID, id, gd.scaleI(lambda));
181 }
182 }
183 }
184
185 return energy;
186 }
187
188 @Override
189 public double getLambda() {
190 return lambda;
191 }
192
193 public Atom[] getAtoms() {
194 return Arrays.copyOf(atoms, atoms.length);
195 }
196
197 @Override
198 public Atom getAtom(int index) {
199 return atoms[index];
200 }
201
202 @Override
203 public boolean applyLambda() {
204 return lambdaTerm;
205 }
206
207 @Override
208 public void setLambda(double lambda) {
209 this.lambda = lamMapper.applyAsDouble(lambda);
210 }
211
212 @Override
213 public double getd2EdL2() {
214 return d2EdL2;
215 }
216
217 @Override
218 public double getdEdL() {
219 return dEdL;
220 }
221
222 public double mapLambda(double lambda) {
223 return lamMapper.applyAsDouble(lambda);
224 }
225
226 @Override
227 public void getdEdXdL(double[] gradient) {
228
229 }
230
231 @Override
232 public String toString() {
233 return format(" Restrain-Torsion %s, Angle: %.3f, Energy: %.3f", torsionType, value, energy);
234 }
235 }