1 //******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2023.
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.numerics.estimator;
39
40 import ffx.utilities.Constants;
41 import ffx.utilities.FFXTest;
42 import org.junit.Test;
43
44 import static org.junit.Assert.assertEquals;
45
46 public class MBARHarmonicOscillatorsTest extends FFXTest {
47
48 /**
49 * Test the MBAR estimator numerics with harmonic oscillators. This test uses SCI & NR, so it
50 * isn't exactly deterministic. The error is set to 1.0E-1, so the first decimal is fine for this.
51 * Pymbar does the same thing.
52 */
53 @Test
54 public void testMBAROscillators() {
55 double[] O_k = {1, 2, 3, 4};
56 double[] K_k = {.5, 1.0, 1.5, 2};
57 int[] N_k = {10000, 10000, 10000, 10000}; // No support for different number of snapshots
58 double beta = 1.0;
59
60 // Create an instance of HarmonicOscillatorsTestCase
61 MultistateBennettAcceptanceRatio.HarmonicOscillatorsTestCase testCase = new MultistateBennettAcceptanceRatio.HarmonicOscillatorsTestCase(O_k, K_k, beta);
62
63 // Generate sample data
64 String setting = "u_kln";
65 Object[] sampleResult = testCase.sample(N_k, setting, (long) 0);
66 double[][][] u_kln = (double[][][]) sampleResult[1];
67 double[] temps = {1 / Constants.R};
68
69 MultistateBennettAcceptanceRatio mbar = new MultistateBennettAcceptanceRatio(O_k, u_kln, temps, 1.0E-7, MultistateBennettAcceptanceRatio.SeedType.ZEROS);
70 double[] mbarFEEstimates = mbar.getMBARFreeEnergies();
71 double[] mbarErrorEstimates = mbar.getBinUncertainties();
72 double[][] mbarDiffMatrix = mbar.getDiffMatrix();
73 double[] mbarFEExpected = new double[]{0.0, 0.3468272332334239, 0.554882810046907, 0.6909139007747198};
74 double[] mbarErrorExpected = new double[]{0.00647778279366289, 0.006176323555016366, 0.008170508071621832};
75 double[][] mbarDiffMatrixExpected = new double[][]{
76 {0.0, 0.00647778279366289, 0.010874539771152386, 0.01568591641036036},
77 {0.005859375, 0.0, 0.006176323555016366, 0.012881744099875898},
78 {0.010697706201272776, 0.00647778279366289, 0.0, 0.008170508071621832},
79 {0.015563845166512918, 0.013028968812623373, 0.008170508071621832, 0.0}
80 };
81
82 // Compare tolerance to pymbar --> I tested this with a for loop over 1000 times and it always passed
83 for (int i = 0; i < mbarFEExpected.length; i++) {
84 assertEquals(mbarFEExpected[i], mbarFEEstimates[i], 1.0E-1);
85 }
86 for (int i = 0; i < mbarErrorExpected.length; i++) {
87 assertEquals(mbarErrorExpected[i], mbarErrorEstimates[i], 1.0E-1);
88 }
89 for (int i = 0; i < mbarDiffMatrixExpected.length; i++) {
90 for (int j = 0; j < mbarDiffMatrixExpected[i].length; j++) {
91 assertEquals(mbarDiffMatrixExpected[i][j], mbarDiffMatrix[i][j], 1.0E-1);
92 }
93 }
94 }
95 }