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.nonbonded;
39
40 import ffx.crystal.Crystal;
41 import ffx.crystal.SpaceGroup;
42 import ffx.crystal.SymOp;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.LambdaInterface;
45 import ffx.potential.parameters.ForceField;
46
47 import java.util.logging.Logger;
48
49 import static ffx.numerics.math.DoubleMath.length2;
50 import static java.lang.String.format;
51 import static java.util.Arrays.fill;
52 import static org.apache.commons.math3.util.FastMath.pow;
53
54
55
56
57
58
59
60
61 public class NCSRestraint implements LambdaInterface {
62
63 private static final Logger logger = Logger.getLogger(NCSRestraint.class.getName());
64 private final Atom[] atoms;
65 private final double forceConstant = 10.0;
66 private final double[][] transOp = new double[3][3];
67 private final double[] a1 = new double[3];
68 private final double[] a2 = new double[3];
69 private final double[] dx = new double[3];
70 private final double lambdaExp = 2.0;
71 private final double[] lambdaGradient;
72 private Crystal ncsCrystal;
73 private SpaceGroup spaceGroup;
74 private int nAtoms;
75 private int nSymm;
76 private double lambda = 1.0;
77 private double lambdaPow = pow(lambda, lambdaExp);
78 private double dLambdaPow = lambdaExp * pow(lambda, lambdaExp - 1.0);
79 private double d2LambdaPow = lambdaExp * (lambdaExp - 1.0) * pow(lambda, lambdaExp - 2.0);
80 private double dEdL = 0.0;
81 private double d2EdL2 = 0.0;
82 private boolean lambdaTerm;
83
84
85
86
87
88
89
90
91
92 public NCSRestraint(Atom[] atoms, ForceField forceField, Crystal crystal) {
93 this.ncsCrystal = crystal;
94 this.atoms = atoms;
95 nAtoms = atoms.length;
96 spaceGroup = this.ncsCrystal.spaceGroup;
97 nSymm = spaceGroup.getNumberOfSymOps();
98 assert (nAtoms % nSymm == 0);
99 lambdaTerm = forceField.getBoolean("LAMBDATERM", false);
100 if (lambdaTerm) {
101 lambdaGradient = new double[nAtoms * 3];
102 } else {
103 lambdaGradient = null;
104 this.lambda = 1.0;
105 lambdaPow = 1.0;
106 dLambdaPow = 0.0;
107 d2LambdaPow = 0.0;
108 }
109
110 logger.info(format("\n NCS Restraint%s", crystal));
111 }
112
113
114 @Override
115 public double getLambda() {
116 return lambda;
117 }
118
119
120 @Override
121 public void setLambda(double lambda) {
122 if (lambdaTerm) {
123 this.lambda = lambda;
124
125 double lambdaWindow = 1.0;
126
127 if (this.lambda < lambdaWindow) {
128 double dldgl = 1.0 / lambdaWindow;
129 double l = dldgl * this.lambda;
130 double l2 = l * l;
131 double l3 = l2 * l;
132 double l4 = l2 * l2;
133 double l5 = l4 * l;
134 double c3 = 10.0;
135 double c4 = -15.0;
136 double c5 = 6.0;
137 double threeC3 = 3.0 * c3;
138 double sixC3 = 6.0 * c3;
139 double fourC4 = 4.0 * c4;
140 double twelveC4 = 12.0 * c4;
141 double fiveC5 = 5.0 * c5;
142 double twentyC5 = 20.0 * c5;
143 lambdaPow = c3 * l3 + c4 * l4 + c5 * l5;
144 dLambdaPow = (threeC3 * l2 + fourC4 * l3 + fiveC5 * l4) * dldgl;
145 d2LambdaPow = (sixC3 * l + twelveC4 * l2 + twentyC5 * l3) * dldgl * dldgl;
146 } else {
147 lambdaPow = 1.0;
148 dLambdaPow = 0.0;
149 d2LambdaPow = 0.0;
150 }
151 } else {
152 this.lambda = 1.0;
153 lambdaPow = 1.0;
154 dLambdaPow = 0.0;
155 d2LambdaPow = 0.0;
156 }
157 }
158
159
160 @Override
161 public double getd2EdL2() {
162 if (lambdaTerm) {
163 return d2EdL2;
164 } else {
165 return 0.0;
166 }
167 }
168
169
170 @Override
171 public double getdEdL() {
172 if (lambdaTerm) {
173 return dEdL;
174 } else {
175 return 0.0;
176 }
177 }
178
179
180 @Override
181 public void getdEdXdL(double[] gradient) {
182 if (lambdaTerm) {
183 int n3 = nAtoms * 3;
184 for (int i = 0; i < n3; i++) {
185 gradient[i] += lambdaGradient[i];
186 }
187 }
188 }
189
190
191
192
193
194
195
196
197 public double residual(boolean gradient, boolean print) {
198
199
200 if (nAtoms % nSymm != 0) {
201 return 0;
202 }
203
204 if (lambdaTerm) {
205 dEdL = 0.0;
206 d2EdL2 = 0.0;
207 fill(lambdaGradient, 0.0);
208 }
209
210 double residual = 0.0;
211 int nAsymmAtoms = nAtoms / nSymm;
212 double fx2 = forceConstant * 2.0;
213 for (int i = 1; i < nSymm; i++) {
214 SymOp symOp = spaceGroup.getSymOp(i);
215 ncsCrystal.getTransformationOperator(symOp, transOp);
216 int offset = nAsymmAtoms * i;
217 for (int j = 0; j < nAsymmAtoms; j++) {
218
219 Atom atom1 = atoms[j];
220 atom1.getXYZ(a1);
221
222 if (atom1.isHydrogen()) {
223 continue;
224 }
225
226
227 ncsCrystal.applySymOp(a1, a1, symOp);
228
229
230 Atom atom2 = atoms[offset + j];
231 atom2.getXYZ(a2);
232
233 dx[0] = a1[0] - a2[0];
234 dx[1] = a1[1] - a2[1];
235 dx[2] = a1[2] - a2[2];
236
237 double r2 = length2(dx);
238
239
240 residual += r2;
241 if (gradient || lambdaTerm) {
242 final double dedx = dx[0] * fx2;
243 final double dedy = dx[1] * fx2;
244 final double dedz = dx[2] * fx2;
245 atom2.addToXYZGradient(-lambdaPow * dedx, -lambdaPow * dedy, -lambdaPow * dedz);
246
247
248 final double dedx1 = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
249 final double dedy1 = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
250 final double dedz1 = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
251 atom1.addToXYZGradient(lambdaPow * dedx1, lambdaPow * dedy1, lambdaPow * dedz1);
252 if (lambdaTerm) {
253 int j3 = j * 3;
254 lambdaGradient[j3] += dLambdaPow * dedx1;
255 lambdaGradient[j3 + 1] += dLambdaPow * dedy1;
256 lambdaGradient[j3 + 2] += dLambdaPow * dedz1;
257 int oj3 = (offset + j) * 3;
258 lambdaGradient[oj3] -= dLambdaPow * dedx;
259 lambdaGradient[oj3 + 1] -= dLambdaPow * dedy;
260 lambdaGradient[oj3 + 2] -= dLambdaPow * dedz;
261 }
262 }
263 }
264 }
265 if (lambdaTerm) {
266 dEdL = dLambdaPow * forceConstant * residual;
267 d2EdL2 = d2LambdaPow * forceConstant * residual;
268 }
269 return forceConstant * residual * lambdaPow;
270 }
271 }