View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * @author Timothy D. Fenn
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      // load the structure
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     // load any properties associated with it
107     CompositeConfiguration properties = molecularAssembly.getProperties();
108 
109     // read in Fo/sigFo/FreeR
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     // associate molecular assembly with the structure, set up forcefield
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     // initialize atomic form factors
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     // set up FFT and run it
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     // delta for numerical finite differences
217     double delta = 1e-4;
218 
219     // compute gradients
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 }