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;
39
40 import static ffx.numerics.math.ScalarMath.b2u;
41 import static ffx.xray.CrystalReciprocalSpace.SolventModel.NONE;
42 import static java.lang.String.format;
43 import static org.junit.Assert.assertEquals;
44 import static org.junit.Assert.assertTrue;
45
46 import edu.rit.pj.ParallelTeam;
47 import ffx.algorithms.misc.AlgorithmsTest;
48 import ffx.crystal.Crystal;
49 import ffx.crystal.ReflectionList;
50 import ffx.crystal.Resolution;
51 import ffx.potential.ForceFieldEnergy;
52 import ffx.potential.MolecularAssembly;
53 import ffx.potential.bonded.Atom;
54 import ffx.potential.utils.PotentialsUtils;
55 import ffx.xray.CrystalReciprocalSpace.SolventModel;
56 import ffx.xray.RefinementMinimize.RefinementMode;
57 import ffx.xray.parsers.MTZFilter;
58
59 import java.io.File;
60 import java.util.Arrays;
61 import java.util.Collection;
62 import java.util.List;
63
64 import org.apache.commons.configuration2.CompositeConfiguration;
65 import org.junit.Test;
66 import org.junit.runner.RunWith;
67 import org.junit.runners.Parameterized;
68 import org.junit.runners.Parameterized.Parameters;
69
70
71
72
73 @RunWith(Parameterized.class)
74 public class FiniteDifferenceTest extends AlgorithmsTest {
75
76 private final boolean ciOnly;
77 private final Atom[] atomArray;
78 private final int[] atoms;
79 private final DiffractionRefinementData refinementData;
80 private final SigmaAMinimize sigmaAMinimize;
81
82 public FiniteDifferenceTest(boolean ciOnly, String info, SolventModel solventModel,
83 int[] atoms, String pdbName, String mtzName) {
84 this.ciOnly = ciOnly;
85 this.atoms = atoms;
86
87 if (!ffxCI && ciOnly) {
88 atomArray = null;
89 refinementData = null;
90 sigmaAMinimize = null;
91 return;
92 }
93
94
95 File structure = getResourceFile(pdbName);
96 File mtzFile = getResourceFile(mtzName);
97 PotentialsUtils potentialsUtils = new PotentialsUtils();
98 MolecularAssembly molecularAssembly = null;
99 try {
100 molecularAssembly = potentialsUtils.open(structure);
101 } catch (Exception ex) {
102 System.err.println(format(" Exception ex: %s", ex));
103 ex.printStackTrace();
104 }
105
106
107 CompositeConfiguration properties = molecularAssembly.getProperties();
108
109
110 MTZFilter mtzFilter = new MTZFilter();
111 ReflectionList reflectionList;
112 Crystal crystal = Crystal.checkProperties(properties);
113 Resolution resolution = Resolution.checkProperties(properties);
114 if (crystal == null || resolution == null) {
115 reflectionList = mtzFilter.getReflectionList(mtzFile);
116 } else {
117 reflectionList = new ReflectionList(crystal, resolution);
118 }
119
120 refinementData = new DiffractionRefinementData(properties, reflectionList);
121 assertTrue(info + " mtz file should be read in without errors",
122 mtzFilter.readFile(mtzFile, reflectionList, refinementData, properties));
123
124
125 molecularAssembly.finalize(true, molecularAssembly.getForceField());
126 ForceFieldEnergy energy = molecularAssembly.getPotentialEnergy();
127
128 List<Atom> atomList = molecularAssembly.getAtomList();
129 atomArray = atomList.toArray(new Atom[0]);
130 boolean use_3g = properties.getBoolean("use-3g", true);
131
132
133 for (Atom atom : atomArray) {
134 XRayFormFactor xrayFormFactor = new XRayFormFactor(atom, use_3g, 2.0);
135 atom.setFormFactorIndex(xrayFormFactor.ffIndex);
136 if (atom.getOccupancy() == 0.0) {
137 atom.setFormFactorWidth(1.0);
138 continue;
139 }
140 double aRad = 2.4;
141 double[] xyz = new double[3];
142 xyz[0] = atom.getX() + aRad;
143 xyz[1] = atom.getY();
144 xyz[2] = atom.getZ();
145 while (true) {
146 double rho = xrayFormFactor.rho(0.0, 1.0, xyz);
147 if (rho > 0.1) {
148 aRad += 0.5;
149 } else if (rho > 0.001) {
150 aRad += 0.1;
151 } else {
152 aRad += 0.75;
153 atom.setFormFactorWidth(aRad);
154 break;
155 }
156 xyz[0] = atom.getX() + aRad;
157 }
158 }
159
160 String gridString = properties.getString("grid-method", "SLICE").toUpperCase();
161 CrystalReciprocalSpace.GridMethod gridMethod;
162 try {
163 gridMethod = CrystalReciprocalSpace.GridMethod.valueOf(gridString);
164 } catch (Exception e) {
165 gridMethod = CrystalReciprocalSpace.GridMethod.SLICE;
166 }
167
168
169 ParallelTeam parallelTeam = new ParallelTeam();
170 CrystalReciprocalSpace crs = new CrystalReciprocalSpace(
171 reflectionList, atomArray, parallelTeam, parallelTeam, false, false, NONE, gridMethod);
172 crs.computeDensity(refinementData.fc);
173 refinementData.setCrystalReciprocalSpaceFc(crs);
174 crs = new CrystalReciprocalSpace(reflectionList, atomArray, parallelTeam, parallelTeam, true,
175 false, solventModel, gridMethod);
176 crs.computeDensity(refinementData.fs);
177 refinementData.setCrystalReciprocalSpaceFs(crs);
178
179 ScaleBulkMinimize scaleBulkMinimize = new ScaleBulkMinimize(reflectionList, refinementData, crs, parallelTeam);
180 scaleBulkMinimize.minimize(6, 1.0e-4);
181
182 sigmaAMinimize = new SigmaAMinimize(reflectionList, refinementData, parallelTeam);
183 sigmaAMinimize.minimize(7, 2.0e-2);
184
185 SplineMinimize splineMinimize = new SplineMinimize(reflectionList, refinementData, refinementData.spline, SplineEnergy.Type.FOFC);
186 splineMinimize.minimize(7, 1e-5);
187
188 CrystalStats crystalstats = new CrystalStats(reflectionList, refinementData);
189 crystalstats.printScaleStats();
190 crystalstats.printHKLStats();
191 crystalstats.printSNStats();
192 crystalstats.printRStats();
193 }
194
195 @Parameters
196 public static Collection<Object[]> data() {
197 return Arrays.asList(
198 new Object[][]{
199 {
200 true,
201 "ala met anisou",
202 NONE,
203 new int[]{91, 105, 119},
204 "alamet.pdb",
205 "alamet.mtz"
206 }
207 });
208 }
209
210 @Test
211 public void testFiniteDifferences() {
212 if (!ffxCI && ciOnly) {
213 return;
214 }
215
216
217 double delta = 1e-4;
218
219
220 refinementData.crystalReciprocalSpaceFc.computeAtomicGradients(refinementData.dFc, refinementData.freeR,
221 refinementData.rFreeFlag, RefinementMode.COORDINATES_AND_BFACTORS);
222
223 double mean = 0.0;
224 double nmean = 0.0;
225 double[] gxyz = new double[3];
226 for (int i = 0; i < atoms.length; i++) {
227 Atom atom = atomArray[atoms[i]];
228 int index = atom.getXyzIndex() - 1;
229 if (atom.getOccupancy() == 0.0) {
230 continue;
231 }
232
233 System.out.println(" " + i + ": " + atom);
234 atom.getXYZGradient(gxyz);
235 double bg = atom.getTempFactorGradient();
236
237 refinementData.crystalReciprocalSpaceFc.deltaX(index, delta);
238 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
239 double llk1 = sigmaAMinimize.calculateLikelihood();
240 refinementData.crystalReciprocalSpaceFc.deltaX(index, -delta);
241 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
242 double llk2 = sigmaAMinimize.calculateLikelihood();
243 double fd = (llk1 - llk2) / (2.0 * delta);
244 System.out.print(format(" X A: %16.8f FD: %16.8f Error: %16.8f\n", gxyz[0], fd, gxyz[0] - fd));
245 refinementData.crystalReciprocalSpaceFc.deltaX(index, 0.0);
246
247 nmean++;
248 mean += (gxyz[0] / fd - mean) / nmean;
249
250 refinementData.crystalReciprocalSpaceFc.deltaY(index, delta);
251 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
252 llk1 = sigmaAMinimize.calculateLikelihood();
253 refinementData.crystalReciprocalSpaceFc.deltaY(index, -delta);
254 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
255 llk2 = sigmaAMinimize.calculateLikelihood();
256 fd = (llk1 - llk2) / (2.0 * delta);
257 System.out.print(format(" Y A: %16.8f FD: %16.8f Error: %16.8f\n", gxyz[1], fd, gxyz[1] - fd));
258 refinementData.crystalReciprocalSpaceFc.deltaY(index, 0.0);
259
260 nmean++;
261 mean += (gxyz[1] / fd - mean) / nmean;
262
263 refinementData.crystalReciprocalSpaceFc.deltaZ(index, delta);
264 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
265 llk1 = sigmaAMinimize.calculateLikelihood();
266 refinementData.crystalReciprocalSpaceFc.deltaZ(index, -delta);
267 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
268 llk2 = sigmaAMinimize.calculateLikelihood();
269 fd = (llk1 - llk2) / (2.0 * delta);
270 System.out.print(format(" Z A: %16.8f FD: %16.8f Error: %16.8f\n", gxyz[2], fd, gxyz[2] - fd));
271 refinementData.crystalReciprocalSpaceFc.deltaZ(index, 0.0);
272
273 nmean++;
274 mean += (gxyz[2] / fd - mean) / nmean;
275
276 if (atom.getAnisou(null) == null) {
277 double b = atom.getTempFactor();
278 atom.setTempFactor(b + delta);
279 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
280 llk1 = sigmaAMinimize.calculateLikelihood();
281 atom.setTempFactor(b - delta);
282 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
283 llk2 = sigmaAMinimize.calculateLikelihood();
284 fd = (llk1 - llk2) / (2.0 * delta);
285 System.out.print(format(" B A: %16.8f FD: %16.8f Error: %16.8f\n", bg, fd, bg - fd));
286 atom.setTempFactor(b);
287 nmean++;
288 mean += (bg / fd - mean) / nmean;
289 } else {
290 double[] anisou = atom.getAnisou(null);
291 double[] anisouG = atom.getAnisouGradient(null);
292 for (int j = 0; j < 6; j++) {
293 double tmpu = anisou[j];
294 anisou[j] = tmpu + b2u(delta);
295 atom.setAnisou(anisou);
296 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
297 llk1 = sigmaAMinimize.calculateLikelihood();
298 anisou[j] = tmpu - b2u(delta);
299 atom.setAnisou(anisou);
300 refinementData.crystalReciprocalSpaceFc.computeDensity(refinementData.fc);
301 llk2 = sigmaAMinimize.calculateLikelihood();
302 fd = (llk1 - llk2) / (2.0 * b2u(delta));
303 atom.getAnisouGradient(anisouG);
304 System.out.print(format(" %d A: %16.8f FD: %16.8f Error: %16.8f\n", j, anisouG[j], fd, anisouG[j] - fd));
305 anisou[j] = tmpu;
306 atom.setAnisou(anisou);
307 nmean++;
308 mean += (anisouG[j] / fd - mean) / nmean;
309 }
310 }
311 }
312
313 System.out.println(" Mean df/fd ratio (" + nmean + "): " + mean);
314 assertEquals(" Gradient: finite-difference ratio should be approx. 1.0", 1.0, mean, 0.01);
315 }
316 }