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.numerics.math.DoubleMath.length;
41 import static ffx.numerics.math.DoubleMath.scale;
42 import static ffx.numerics.math.DoubleMath.sub;
43
44 import ffx.crystal.Crystal;
45 import ffx.numerics.atomic.AtomicDoubleArray3D;
46 import ffx.numerics.switching.ConstantSwitch;
47 import ffx.numerics.switching.UnivariateSwitchingFunction;
48 import ffx.potential.parameters.BondType;
49
50 import java.io.Serial;
51 import java.util.logging.Logger;
52
53
54
55
56
57
58
59 public class RestraintBond extends BondedTerm implements LambdaInterface {
60
61 @Serial
62 private static final long serialVersionUID = 1L;
63
64 public static final double DEFAULT_RB_LAM_START = 0.75;
65 public static final double DEFAULT_RB_LAM_END = 1.0;
66 private static final Logger logger = Logger.getLogger(RestraintBond.class.getName());
67 private final double restraintLambdaStart;
68 private final double restraintLambdaStop;
69 private final double restraintLambdaWindow;
70 private final double rlwInv;
71 private final UnivariateSwitchingFunction switchingFunction;
72 public BondType bondType = null;
73 private final boolean lambdaTerm;
74 private double lambda = 1.0;
75 private double switchVal = 1.0;
76 private double switchdUdL = 1.0;
77 private double switchd2UdL2 = 1.0;
78 private double dEdL = 0.0;
79 private double d2EdL2 = 0.0;
80 private final double[][] dEdXdL = new double[2][3];
81 private final Crystal crystal;
82
83
84
85
86
87
88
89
90
91
92
93
94 public RestraintBond(Atom a1, Atom a2, Crystal crystal, boolean lambdaTerm, double lamStart,
95 double lamEnd, UnivariateSwitchingFunction sf) {
96 restraintLambdaStart = lamStart;
97 restraintLambdaStop = lamEnd;
98 assert lamEnd > lamStart;
99 restraintLambdaWindow = lamEnd - lamStart;
100 rlwInv = 1.0 / restraintLambdaWindow;
101
102 atoms = new Atom[2];
103
104 this.crystal = crystal;
105
106 int i1 = a1.getIndex();
107 int i2 = a2.getIndex();
108 if (i1 < i2) {
109 atoms[0] = a1;
110 atoms[1] = a2;
111 } else {
112 atoms[0] = a2;
113 atoms[1] = a1;
114 }
115 setID_Key(false);
116 switchingFunction = (sf == null ? new ConstantSwitch() : sf);
117 this.lambdaTerm = lambdaTerm;
118 this.setLambda(1.0);
119 }
120
121
122
123
124
125
126 @Override
127 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
128 value = 0.0;
129 energy = 0.0;
130
131 if (!getUse()) {
132 return energy;
133 }
134 Atom atomA = atoms[0];
135 Atom atomB = atoms[1];
136 double[] a0 = new double[3];
137 double[] a1 = new double[3];
138
139 double[] v10 = new double[3];
140
141 double[] g0 = new double[3];
142 double[] g1 = new double[3];
143 atomA.getXYZ(a0);
144 atomB.getXYZ(a1);
145 sub(a0, a1, v10);
146 if (crystal != null) {
147 crystal.image(v10);
148 }
149
150 value = length(v10);
151 double dv = value - bondType.distance;
152 if (bondType.bondFunction.hasFlatBottom()) {
153 if (dv > 0) {
154 dv = Math.max(0, dv - bondType.flatBottomRadius);
155 } else if (dv < 0) {
156 dv = Math.min(0, dv + bondType.flatBottomRadius);
157 }
158 }
159
160 double dv2 = dv * dv;
161 double kx2 = bondType.bondUnit * bondType.forceConstant * dv2;
162 energy = switchVal * kx2;
163 dEdL = switchdUdL * kx2;
164 d2EdL2 = switchd2UdL2 * kx2;
165 double deddt = 2.0 * bondType.bondUnit * bondType.forceConstant * dv;
166 double de = 0.0;
167
168 if (value > 0.0) {
169 de = deddt / value;
170 }
171
172 scale(v10, switchVal * de, g0);
173 scale(v10, -switchVal * de, g1);
174 if (gradient) {
175 grad.add(threadID, atoms[0].getIndex() - 1, g0[0], g0[1], g0[2]);
176 grad.add(threadID, atoms[1].getIndex() - 1, g1[0], g1[1], g1[2]);
177 }
178
179
180 scale(v10, switchdUdL * de, g0);
181 scale(v10, -switchdUdL * de, g1);
182 dEdXdL[0][0] = g0[0];
183 dEdXdL[0][1] = g0[1];
184 dEdXdL[0][2] = g0[2];
185 dEdXdL[1][0] = g1[0];
186 dEdXdL[1][1] = g1[1];
187 dEdXdL[1][2] = g1[2];
188
189 value = dv;
190 return energy;
191 }
192
193
194
195
196
197
198
199
200 public Atom get1_2(Atom a) {
201 if (a == atoms[0]) {
202 return atoms[1];
203 }
204 if (a == atoms[1]) {
205 return atoms[0];
206 }
207 return null;
208 }
209
210
211
212
213
214
215 public BondType getBondType() {
216 return bondType;
217 }
218
219
220
221
222
223
224 public void setBondType(BondType bondType) {
225 this.bondType = bondType;
226 }
227
228
229 @Override
230 public double getLambda() {
231 return lambda;
232 }
233
234
235 @Override
236 public void setLambda(double lambda) {
237 this.lambda = lambda;
238 if (lambdaTerm) {
239 if (lambda < restraintLambdaStart) {
240 switchVal = 0.0;
241 switchdUdL = 0;
242 switchd2UdL2 = 0;
243 } else if (lambda > restraintLambdaStop) {
244 switchVal = 1.0;
245 switchdUdL = 0;
246 switchd2UdL2 = 0;
247 } else {
248 double restraintLambda = (lambda - restraintLambdaStart) / restraintLambdaWindow;
249 switchVal = switchingFunction.valueAt(restraintLambda);
250 switchdUdL = rlwInv * switchingFunction.firstDerivative(restraintLambda);
251 switchd2UdL2 = rlwInv * rlwInv * switchingFunction.secondDerivative(restraintLambda);
252 }
253 } else {
254 switchVal = 1.0;
255 switchdUdL = 0.0;
256 switchd2UdL2 = 0.0;
257 }
258 }
259
260
261 @Override
262 public double getd2EdL2() {
263 return d2EdL2;
264 }
265
266
267 @Override
268 public double getdEdL() {
269
270 return dEdL;
271 }
272
273
274 @Override
275 public void getdEdXdL(double[] gradient) {
276 int i1 = atoms[0].getIndex() - 1;
277 int index = i1 * 3;
278 gradient[index++] += dEdXdL[0][0];
279 gradient[index++] += dEdXdL[0][1];
280 gradient[index] += dEdXdL[0][2];
281 int i2 = atoms[1].getIndex() - 1;
282 index = i2 * 3;
283 gradient[index++] += dEdXdL[1][0];
284 gradient[index++] += dEdXdL[1][1];
285 gradient[index] += dEdXdL[1][2];
286 }
287
288 @Override
289 public boolean isLambdaScaled() {
290 return lambdaTerm;
291 }
292
293
294 public void log() {
295 logger.info(String.format(" %s %6d-%s %6d-%s %6.4f %6.4f %10.4f", "Restraint-Bond",
296 atoms[0].getIndex(), atoms[0].getAtomType().name, atoms[1].getIndex(),
297 atoms[1].getAtomType().name, bondType.distance, value, energy));
298 if (!(switchingFunction instanceof ConstantSwitch)) {
299 logger.info(String.format(" Switching function (lambda dependence): %s",
300 switchingFunction.toString()));
301 }
302 }
303
304 @Override
305 public String toString() {
306 StringBuilder sb = new StringBuilder(String.format(
307 " Distance restraint between atoms %s-%d %s-%d, "
308 + "current distance %10.4g, optimum %10.4g with a %10.4g Angstrom flat bottom, with force constant %10.4g.",
309 atoms[0], atoms[0].getIndex(), atoms[1], atoms[1].getIndex(), value, bondType.distance,
310 bondType.flatBottomRadius, bondType.forceConstant));
311 if (!(switchingFunction instanceof ConstantSwitch)) {
312 sb.append(String.format("\n Switching function (lambda dependence): %s",
313 switchingFunction.toString()));
314 }
315 return sb.toString();
316 }
317 }