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.utils;
39
40 import ffx.numerics.Potential;
41 import ffx.potential.cli.GradientOptions;
42
43 import java.util.ArrayList;
44 import java.util.List;
45 import java.util.logging.Logger;
46
47 import static ffx.utilities.StringUtils.parseAtomRanges;
48 import static java.lang.String.format;
49 import static org.apache.commons.math3.util.FastMath.sqrt;
50
51
52
53
54
55
56
57 public class GradientUtils {
58
59 private static final Logger logger = Logger.getLogger(GradientUtils.class.getName());
60
61 private final Potential potential;
62
63
64
65
66
67
68 public GradientUtils(Potential potential) {
69 this.potential = potential;
70 }
71
72
73
74
75
76
77
78 public int testGradient(GradientOptions gradientOptions) {
79 int nFailures = 0;
80
81
82 double step = gradientOptions.getDx();
83 logger.info(" Finite-difference step size:\t" + step);
84
85
86 boolean print = gradientOptions.getVerbose();
87 logger.info(" Verbose printing:\t\t" + print);
88
89
90 List<Integer> degreesOfFreedomToTest;
91 int nAtoms = potential.getNumberOfVariables() / 3;
92 String gradientAtoms = gradientOptions.getGradientAtoms();
93 if (gradientAtoms.equalsIgnoreCase("NONE")) {
94 logger.info(" The gradient of no atoms was evaluated.");
95 return nFailures;
96 } else if (gradientAtoms.equalsIgnoreCase("ALL")) {
97 logger.info(" Checking gradient for all degrees of freedom.\n");
98 degreesOfFreedomToTest = new ArrayList<>();
99 for (int i = 0; i < nAtoms; i++) {
100 degreesOfFreedomToTest.add(i);
101 }
102 } else {
103 degreesOfFreedomToTest = parseAtomRanges(" Gradient atoms", gradientAtoms, nAtoms);
104 logger.info(" Checking gradient for degrees of freedom in the range: " + gradientAtoms + "\n");
105 }
106
107
108 int n = potential.getNumberOfVariables();
109 double[] x = new double[n];
110 double[] g = new double[n];
111 potential.getCoordinates(x);
112 potential.energyAndGradient(x, g);
113
114
115 double expGrad = 1000.0;
116 double gradientTolerance = gradientOptions.getTolerance();
117 double width = 2.0 * step;
118 double avLen = 0.0;
119 double avGrad = 0.0;
120 double expGrad2 = expGrad * expGrad;
121
122 int nTested = 0;
123 int index = 0;
124 double[] numeric = new double[3];
125 for (int k = 0; k < nAtoms; k++) {
126 int i0 = index++;
127 int i1 = index++;
128 int i2 = index++;
129 if (!degreesOfFreedomToTest.contains(k)) {
130 continue;
131 }
132 nTested++;
133
134
135 double orig = x[i0];
136 x[i0] = x[i0] + step;
137 double e = potential.energy(x);
138 x[i0] = orig - step;
139 e -= potential.energy(x);
140 x[i0] = orig;
141 numeric[0] = e / width;
142
143
144 orig = x[i1];
145 x[i1] = x[i1] + step;
146 e = potential.energy(x);
147 x[i1] = orig - step;
148 e -= potential.energy(x);
149 x[i1] = orig;
150 numeric[1] = e / width;
151
152
153 orig = x[i2];
154 x[i2] = x[i2] + step;
155 e = potential.energy(x);
156 x[i2] = orig - step;
157 e -= potential.energy(x);
158 x[i2] = orig;
159 numeric[2] = e / width;
160
161 double dx = g[i0] - numeric[0];
162 double dy = g[i1] - numeric[1];
163 double dz = g[i2] - numeric[2];
164 double len = dx * dx + dy * dy + dz * dz;
165 avLen += len;
166 len = sqrt(len);
167
168 double grad2 = g[i0] * g[i0] + g[i1] * g[i1] + g[i2] * g[i2];
169 avGrad += grad2;
170
171 if (len > gradientTolerance) {
172 logger.info(format(" Degree of Freedom %d\n Failed: %10.6f\n", k + 1, len) +
173 format(" Analytic: (%12.4f, %12.4f, %12.4f)\n", g[i0], g[i1], g[i2]) +
174 format(" Numeric: (%12.4f, %12.4f, %12.4f)\n", numeric[0], numeric[1], numeric[2]));
175 ++nFailures;
176 } else {
177 logger.info(format(" Degree of Freedom %d\n Passed: %10.6f\n", k + 1, len) +
178 format(" Analytic: (%12.4f, %12.4f, %12.4f)\n", g[i0], g[i1], g[i2]) +
179 format(" Numeric: (%12.4f, %12.4f, %12.4f)", numeric[0], numeric[1], numeric[2]));
180 }
181
182 if (grad2 > expGrad2) {
183 logger.info(format(" Degree of Freedom %d has an unusually large gradient: %10.6f", k + 1, Math.sqrt(grad2)));
184 }
185 logger.info("\n");
186 }
187
188 avLen = avLen / nTested;
189 avLen = sqrt(avLen);
190 if (avLen > gradientTolerance) {
191 logger.info(format(" Test failure: RMSD from analytic solution is %10.6f > %10.6f", avLen, gradientTolerance));
192 } else {
193 logger.info(format(" Test success: RMSD from analytic solution is %10.6f < %10.6f", avLen, gradientTolerance));
194 }
195 logger.info(format(" Number of atoms failing analytic test: %d", nFailures));
196
197 avGrad = avGrad / nTested;
198 avGrad = sqrt(avGrad);
199 if (avGrad > expGrad) {
200 logger.info(format(" Unusually large RMS gradient: %10.6f > %10.6f", avGrad, expGrad));
201 } else {
202 logger.info(format(" RMS gradient: %10.6f", avGrad));
203 }
204
205 return nFailures;
206 }
207 }