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