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 org.apache.commons.math3.util.FastMath.acos;
41 import static org.apache.commons.math3.util.FastMath.max;
42 import static org.apache.commons.math3.util.FastMath.min;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44 import static org.apache.commons.math3.util.FastMath.toDegrees;
45
46 import ffx.numerics.atomic.AtomicDoubleArray3D;
47 import ffx.potential.parameters.ForceField;
48 import ffx.potential.parameters.StretchBendType;
49
50 import java.io.Serial;
51 import java.util.logging.Logger;
52
53
54
55
56
57
58
59 public class StretchBend extends BondedTerm implements Comparable<BondedTerm> {
60
61 @Serial
62 private static final long serialVersionUID = 1L;
63
64 private static final Logger logger = Logger.getLogger(StretchBend.class.getName());
65
66 public double angleEq;
67
68 public double bond0Eq;
69
70 public double bond1Eq;
71
72 protected final Angle angle;
73
74 public double force0, force1;
75
76 private StretchBendType stretchBendType = null;
77
78 private double rigidScale = 1.0;
79
80
81
82
83
84
85 public StretchBend(Angle a) {
86 super();
87 angle = a;
88 atoms = a.atoms;
89 bonds = a.bonds;
90 angleEq = angle.angleType.angle[angle.nh];
91 bond0Eq = bonds[0].bondType.distance;
92 bond1Eq = bonds[1].bondType.distance;
93 setID_Key(false);
94 }
95
96
97
98
99
100
101
102
103 static StretchBend stretchBendFactory(Angle angle, ForceField forceField) {
104 StretchBendType stretchBendType = forceField.getStretchBendType(angle.getAngleType().getKey());
105 if (stretchBendType == null) {
106 return null;
107 }
108 StretchBend stretchBend = new StretchBend(angle);
109 stretchBend.setStretchBendType(stretchBendType);
110 return stretchBend;
111 }
112
113
114 @Override
115 public int compareTo(BondedTerm sb) {
116 if (!sb.getClass().isInstance(this)) {
117 return super.compareTo(sb);
118 }
119 return angle.compareTo(((StretchBend) sb).angle);
120 }
121
122
123
124
125
126
127 @Override
128 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
129 energy = 0.0;
130 value = 0.0;
131
132 if (!getUse()) {
133 return energy;
134 }
135 var atomA = atoms[0];
136 var atomB = atoms[1];
137 var atomC = atoms[2];
138 var ia = atomA.getIndex() - 1;
139 var ib = atomB.getIndex() - 1;
140 var ic = atomC.getIndex() - 1;
141 var va = atomA.getXYZ();
142 var vb = atomB.getXYZ();
143 var vc = atomC.getXYZ();
144 var vab = va.sub(vb);
145 var vcb = vc.sub(vb);
146 var rab2 = vab.length2();
147 var rcb2 = vcb.length2();
148 if (rab2 != 0.0 && rcb2 != 0.0) {
149 var rab = sqrt(rab2);
150 var rcb = sqrt(rcb2);
151 var vp = vcb.X(vab);
152 var rp = max(vp.length(), 0.000001);
153 value = toDegrees(acos(min(1.0, max(-1.0, vab.dot(vcb) / (rab * rcb)))));
154 var e0 = rab - bond0Eq;
155 var e1 = rcb - bond1Eq;
156 var dt = value - angleEq;
157 var dr = force0 * e0 + force1 * e1;
158 var prefactor = rigidScale;
159 energy = prefactor * dr * dt;
160 if (gradient) {
161 var vdta = vab.X(vp).scaleI(-prefactor * dr * toDegrees(1.0 / (rab2 * rp)));
162 var vdtc = vcb.X(vp).scaleI(prefactor * dr * toDegrees(1.0 / (rcb2 * rp)));
163 var ga = vdta.addI(vab.scaleI(prefactor * force0 * dt / rab));
164 var gc = vdtc.addI(vcb.scaleI(prefactor * force1 * dt / rcb));
165 grad.add(threadID, ia, ga);
166 grad.sub(threadID, ib, ga.add(gc));
167 grad.add(threadID, ic, gc);
168 }
169 }
170 return energy;
171 }
172
173
174 public void log() {
175 logger.info(String.format(" %s %6d-%s %6d-%s %6d-%s" + "%7.4f %10.4f", "Stretch-Bend",
176 atoms[0].getIndex(), atoms[0].getAtomType().name, atoms[1].getIndex(),
177 atoms[1].getAtomType().name, atoms[2].getIndex(), atoms[2].getAtomType().name, value,
178 energy));
179 }
180
181
182
183
184
185
186 public void setRigidScale(double rigidScale) {
187 this.rigidScale = rigidScale;
188 }
189
190
191
192
193
194
195 public void setStretchBendType(StretchBendType stretchBendType) {
196 this.stretchBendType = stretchBendType;
197
198 if (atoms[0].getAtomType().atomClass == stretchBendType.atomClasses[0]) {
199 force0 = stretchBendType.strbndunit * stretchBendType.forceConstants[0];
200 force1 = stretchBendType.strbndunit * stretchBendType.forceConstants[1];
201 } else {
202 force0 = stretchBendType.strbndunit * stretchBendType.forceConstants[1];
203 force1 = stretchBendType.strbndunit * stretchBendType.forceConstants[0];
204 }
205 atoms = angle.atoms;
206 bonds = angle.bonds;
207 angleEq = angle.angleType.angle[angle.nh];
208 bond0Eq = bonds[0].bondType.distance;
209 bond1Eq = bonds[1].bondType.distance;
210 }
211
212
213
214
215
216
217 @Override
218 public String toString() {
219 return String.format("%s (%7.2f,%7.2f,%7.1f,%7.2f)", id, bonds[0].value, bonds[1].value,
220 angle.value, energy);
221 }
222 }