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 static ffx.potential.parameters.AngleType.AngleMode.IN_PLANE;
41 import static org.apache.commons.math3.util.FastMath.acos;
42 import static org.apache.commons.math3.util.FastMath.max;
43 import static org.apache.commons.math3.util.FastMath.min;
44 import static org.apache.commons.math3.util.FastMath.signum;
45 import static org.apache.commons.math3.util.FastMath.sqrt;
46 import static org.apache.commons.math3.util.FastMath.toDegrees;
47
48 import ffx.numerics.atomic.AtomicDoubleArray3D;
49 import ffx.potential.parameters.ForceField;
50 import ffx.potential.parameters.OutOfPlaneBendType;
51
52 import java.io.Serial;
53 import java.util.logging.Logger;
54
55
56
57
58
59
60
61 public class OutOfPlaneBend extends BondedTerm {
62
63 @Serial
64 private static final long serialVersionUID = 1L;
65
66 private static final Logger logger = Logger.getLogger(OutOfPlaneBend.class.getName());
67
68 public OutOfPlaneBendType outOfPlaneBendType = null;
69
70
71
72
73
74
75
76 public OutOfPlaneBend(Angle angle, Atom atom) {
77 super();
78 atoms = new Atom[4];
79 atoms[0] = angle.atoms[0];
80 atoms[1] = angle.atoms[1];
81 atoms[2] = angle.atoms[2];
82 atoms[3] = atom;
83 bonds = new Bond[3];
84 bonds[0] = angle.bonds[0];
85 bonds[1] = angle.bonds[1];
86 bonds[2] = atoms[1].getBond(atom);
87 setID_Key(false);
88 }
89
90
91
92
93
94
95 public Atom getFourthAtom() {
96 return atoms[3];
97 }
98
99
100
101
102
103
104 public Atom getTrigonalAtom() {
105 return atoms[1];
106 }
107
108
109
110
111
112
113 public Atom getFirstAngleAtom() {
114 return atoms[0];
115 }
116
117
118
119
120
121
122 public Atom getLastAngleAtom() {
123 return atoms[2];
124 }
125
126
127
128
129
130
131
132
133
134 public static OutOfPlaneBend outOfPlaneBendFactory(Angle angle, ForceField forceField) {
135 Atom centralAtom = angle.atoms[1];
136 if (centralAtom.isTrigonal()) {
137 Atom fourthAtom = angle.getFourthAtomOfTrigonalCenter();
138 Atom[] atoms = angle.atoms;
139
140 OutOfPlaneBendType outOfPlaneBendType = forceField.getOutOfPlaneBendType(
141 fourthAtom.getAtomType(), atoms[0].getAtomType(), atoms[1].getAtomType(),
142 atoms[2].getAtomType());
143
144 if (outOfPlaneBendType != null) {
145 if (angle.getAngleMode() == IN_PLANE) {
146 angle.setInPlaneAtom(fourthAtom);
147 }
148 OutOfPlaneBend outOfPlaneBend = new OutOfPlaneBend(angle, fourthAtom);
149 outOfPlaneBend.setOutOfPlaneBendType(outOfPlaneBendType);
150 return outOfPlaneBend;
151 }
152 }
153 return null;
154 }
155
156
157 @Override
158 public int compareTo(BondedTerm o) {
159 if (o == null) {
160 throw new NullPointerException();
161 }
162 if (o == this) {
163 return 0;
164 }
165 if (!o.getClass().isInstance(this)) {
166 return super.compareTo(o);
167 }
168 int this1 = atoms[1].getIndex();
169 int a1 = o.atoms[1].getIndex();
170 if (this1 < a1) {
171 return -1;
172 }
173 if (this1 > a1) {
174 return 1;
175 }
176 int this3 = atoms[3].getIndex();
177 int a3 = o.atoms[3].getIndex();
178 return Integer.compare(this3, a3);
179 }
180
181
182
183
184
185
186 @Override
187 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
188 energy = 0.0;
189 value = 0.0;
190
191 if (!getUse()) {
192 return energy;
193 }
194 var atomA = atoms[0];
195 var atomB = atoms[1];
196 var atomC = atoms[2];
197 var atomD = atoms[3];
198 var va = atomA.getXYZ();
199 var vb = atomB.getXYZ();
200 var vc = atomC.getXYZ();
201 var vd = atomD.getXYZ();
202 var vab = va.sub(vb);
203 var vcb = vc.sub(vb);
204 var vdb = vd.sub(vb);
205 var vad = va.sub(vd);
206 var vcd = vc.sub(vd);
207 var rdb2 = vdb.length2();
208 var rad2 = vad.length2();
209 var rcd2 = vcd.length2();
210 var vp = vcb.X(vdb);
211 var ee = vab.dot(vp);
212 var rac2 = vad.dot(vcd);
213 var cc = rad2 * rcd2 - rac2 * rac2;
214 if (rdb2 != 0.0 && cc != 0.0) {
215 var bkk2 = rdb2 - ee * ee / cc;
216 var cosine = min(1.0, max(-1.0, sqrt(bkk2 / rdb2)));
217 value = toDegrees(acos(cosine));
218 var dv = value;
219 var dv2 = dv * dv;
220 var dv3 = dv2 * dv;
221 var dv4 = dv2 * dv2;
222 energy = outOfPlaneBendType.opBendUnit * outOfPlaneBendType.forceConstant * dv2 * (1.0
223 + outOfPlaneBendType.cubic * dv + outOfPlaneBendType.quartic * dv2
224 + outOfPlaneBendType.pentic * dv3 + outOfPlaneBendType.sextic * dv4);
225 if (gradient) {
226 var deddt = outOfPlaneBendType.opBendUnit * outOfPlaneBendType.forceConstant * dv * toDegrees(
227 2.0 + 3.0 * outOfPlaneBendType.cubic * dv + 4.0 * outOfPlaneBendType.quartic * dv2
228 + 5.0 * outOfPlaneBendType.pentic * dv3 + 6.0 * outOfPlaneBendType.sextic * dv4);
229 var dedcos = 0.0;
230 if (ee != 0.0) {
231 dedcos = -deddt * signum(ee) / sqrt(cc * bkk2);
232 } else {
233 dedcos = -deddt / sqrt(cc * bkk2);
234 }
235 var term = ee / cc;
236
237
238 var svad = vad.scale(rcd2);
239 var svcd = vcd.scale(rac2);
240 var dcda = svad.sub(svcd).scaleI(term);
241 svad = vad.scale(rac2);
242 svcd = vcd.scale(rad2);
243 var dadc = svcd.sub(svad).scaleI(term);
244 var dcdd = dcda.add(dadc).scaleI(-1.0);
245 var deda = vdb.X(vcb);
246 var dedc = vab.X(vdb);
247 var dedd = vcb.X(vab);
248 dedd.addI(vdb.scaleI(ee / rdb2));
249
250
251 var ga = dcda.add(deda).scaleI(dedcos);
252 var gc = dadc.add(dedc).scaleI(dedcos);
253 var gd = dcdd.add(dedd).scaleI(dedcos);
254 var gb = ga.add(gc).addI(gd).scaleI(-1.0);
255 var ia = atomA.getIndex() - 1;
256 var ib = atomB.getIndex() - 1;
257 var ic = atomC.getIndex() - 1;
258 var id = atomD.getIndex() - 1;
259 grad.add(threadID, ia, ga);
260 grad.add(threadID, ib, gb);
261 grad.add(threadID, ic, gc);
262 grad.add(threadID, id, gd);
263 }
264 }
265 return energy;
266 }
267
268
269 public void log() {
270 logger.info(
271 String.format(" %s %6d-%s %6d-%s %6.4f %10.4f", "Out-of-Plane Bend", atoms[1].getIndex(),
272 atoms[1].getAtomType().name, atoms[3].getIndex(), atoms[3].getAtomType().name, value,
273 energy));
274 }
275
276
277
278
279
280
281 public void setOutOfPlaneBendType(OutOfPlaneBendType a) {
282 outOfPlaneBendType = a;
283 }
284
285
286
287
288
289
290 @Override
291 public String toString() {
292 return String.format("%s (%7.1f,%7.2f)", id, value, energy);
293 }
294 }