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.xray.commands.test;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.numerics.Potential;
42 import ffx.potential.MolecularAssembly;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.cli.AtomSelectionOptions;
45 import ffx.potential.cli.PotentialCommand;
46 import ffx.utilities.FFXBinding;
47 import ffx.xray.DiffractionData;
48 import ffx.xray.refine.RefinedParameter;
49 import ffx.xray.refine.RefinementMode;
50 import ffx.xray.refine.RefinementModel;
51 import ffx.xray.XRayEnergy;
52 import ffx.xray.cli.XrayOptions;
53 import org.apache.commons.configuration2.CompositeConfiguration;
54 import picocli.CommandLine.Command;
55 import picocli.CommandLine.Mixin;
56 import picocli.CommandLine.Option;
57 import picocli.CommandLine.Parameters;
58
59 import java.util.ArrayList;
60 import java.util.List;
61
62 import static ffx.utilities.StringUtils.parseAtomRanges;
63 import static java.lang.String.format;
64 import static org.apache.commons.math3.util.FastMath.sqrt;
65
66
67
68
69
70
71
72
73 @Command(description = " Test the potential energy gradient.", name = "test.Gradient")
74 public class Gradient extends AlgorithmsCommand {
75
76 @Mixin
77 AtomSelectionOptions atomSelectionOptions;
78
79 @Mixin
80 XrayOptions xrayOptions;
81
82
83
84
85 @Option(names = {"-d", "--dx"}, defaultValue = "1.0e-5", paramLabel = "1.0e-5",
86 description = "Finite-difference step size.")
87 public double dx = 1.0e-5;
88
89
90
91
92 @Option(names = {"--tol", "--tolerance"}, defaultValue = "1.0e-3",
93 paramLabel = "1.0e-3 kcal/mol/param", description = "Gradient error tolerance.")
94 public double tolerance = 1.0e-3;
95
96
97
98
99 @Option(names = {"--params", "--gradientParams"}, paramLabel = "ALL", defaultValue = "ALL",
100 description = "Ranges of degrees of freedom to test [ALL, NONE, Range(s): 1-3,6-N].")
101 public String gradientParams = "ALL";
102
103
104
105
106
107 @Parameters(arity = "1..*", paramLabel = "files", description = "PDB and Diffraction input files.")
108 private List<String> filenames;
109
110 private MolecularAssembly[] molecularAssemblies;
111 private DiffractionData diffractionData;
112 private XRayEnergy xrayEnergy;
113 public int nFailures = 0;
114
115
116
117
118 public Gradient() {
119 super();
120 }
121
122
123
124
125
126
127 public Gradient(FFXBinding binding) {
128 super(binding);
129 }
130
131
132
133
134
135
136 public Gradient(String[] args) {
137 super(args);
138 }
139
140
141
142
143 @Override
144 public Gradient run() {
145
146 if (!init()) {
147 return this;
148 }
149
150 xrayOptions.init();
151
152 String filename;
153 if (filenames != null && !filenames.isEmpty()) {
154
155 molecularAssemblies = algorithmFunctions.openAll(filenames.getFirst());
156 activeAssembly = molecularAssemblies[0];
157 } else if (activeAssembly == null) {
158 logger.info(helpString());
159 return this;
160 } else {
161 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
162 }
163
164
165 filename = activeAssembly.getFile().getAbsolutePath();
166
167
168 for (MolecularAssembly assembly : molecularAssemblies) {
169 atomSelectionOptions.setActiveAtoms(assembly);
170 }
171
172
173 CompositeConfiguration properties = activeAssembly.getProperties();
174 xrayOptions.setProperties(parseResult, properties);
175
176
177 diffractionData = xrayOptions.getDiffractionData(filenames, molecularAssemblies, properties);
178 diffractionData.scaleBulkFit();
179 diffractionData.printStats();
180 xrayEnergy = new XRayEnergy(diffractionData);
181
182 RefinementModel refinementModel = diffractionData.getRefinementModel();
183 RefinementMode refinementMode = refinementModel.getRefinementMode();
184
185 logger.info("\n Testing the gradient of " + filename);
186 logger.info(refinementMode.toString());
187 logger.info(refinementModel.toString());
188
189 int n = refinementModel.getNumParameters();
190 double[] x = new double[n];
191 double[] g = new double[n];
192
193 refinementModel.getParameters(x);
194 xrayEnergy.energyAndGradient(x, g);
195
196 List<RefinedParameter> refinedParameters = refinementModel.getRefinedParameters();
197 int numParameters = refinedParameters.size();
198
199 List<Integer> parametersToTest;
200 if (gradientParams.equalsIgnoreCase("NONE")) {
201 logger.info(" The gradient of no parameters will be evaluated.");
202 return this;
203 } else if (gradientParams.equalsIgnoreCase("ALL")) {
204 logger.info(" Checking the gradient for all parameters.\n");
205 parametersToTest = new ArrayList<>();
206 for (int i = 0; i < numParameters; i++) {
207 parametersToTest.add(i);
208 }
209 } else {
210 parametersToTest = parseAtomRanges(" Parameters", gradientParams, numParameters);
211 logger.info(" Checking the gradient for parameters in the range: " + gradientParams + "\n");
212 }
213
214 for (Integer parameter : parametersToTest) {
215 RefinedParameter refinedParameter = refinedParameters.get(parameter);
216 int index = refinedParameter.getIndex();
217 int numParams = refinedParameter.getNumberOfParameters();
218 double[] grad = new double[numParams];
219 double error = 0.0;
220
221 for (int i = 0; i < numParams; i++) {
222 double orig = x[index + i];
223 x[index + i] = orig + dx;
224 double eplus = xrayEnergy.energy(x);
225 x[index + i] = orig - dx;
226 double eminus = xrayEnergy.energy(x);
227 x[index + i] = orig;
228 grad[i] = (eplus - eminus) / (2.0 * dx);
229 double diff = grad[i] - g[index + i];
230 error += diff * diff;
231 }
232 error = sqrt(error);
233 Atom atom = refinedParameter.getAtom();
234 if (numParams == 1) {
235 logger.info(format("\n Refinement parameter %d with index %d: %s", parameter + 1, index, atom));
236 } else {
237 logger.info(format("\n Refinement parameter %d with indices %d - %s: %s",
238 parameter + 1, index, index + numParams - 1, atom));
239 }
240 logger.info(format(" %s", refinedParameter));
241 if (error > tolerance) {
242 logger.info(format(" Failed: %10.6f", error));
243 nFailures++;
244 } else {
245 logger.info(format(" Passed: %10.6f", error));
246 }
247 StringBuilder numeric = new StringBuilder(" Numeric: ");
248 StringBuilder analytic = new StringBuilder(" Analytic: ");
249 for (int i = 0; i < numParams; i++) {
250 numeric.append(format(" %10.6f", grad[i]));
251 analytic.append(format(" %10.6f", g[index + i]));
252 }
253 logger.info(numeric.toString());
254 logger.info(analytic.toString());
255 }
256
257 logger.info("\n %d / %d failures.".formatted(nFailures, numParameters));
258
259 return this;
260 }
261
262 @Override
263 public List<Potential> getPotentials() {
264 List<Potential> potentials = getPotentialsFromAssemblies(molecularAssemblies);
265 if (xrayEnergy != null) {
266 potentials.add(xrayEnergy);
267 }
268 return potentials;
269 }
270
271 }