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