Fix for non-default (not lambda 1.0->0.0) LRA calculations.

The physics behind calculating LRA contributions was simplified for the default 1.0->0.0 case, causing all calculations with other lambda values to be invalid. This does not include Group Contribution calculations (q_calc).
mpurg · Aug 15, 2017 · ca2125e · ca2125e
1 parent 30f1fc1
commit ca2125e
Show file tree

Hide file tree

Showing 13 changed files with 515 additions and 481 deletions.
diff --git a/README.md b/README.md
@@ -98,7 +98,7 @@ DOI of the latest release:
 [![DOI](https://zenodo.org/badge/80016679.svg)](https://zenodo.org/badge/latestdoi/80016679)  
 
 *Example:*  
-*...analysis was performed with qtools v0.5.9 (DOI: 10.5281/zenodo.842003).*
+*...analysis was performed with qtools v0.5.10 (DOI: 10.5281/zenodo.842003).*
 
 
 #### Author

diff --git a/packages/Qpyl/common.py b/packages/Qpyl/common.py
@@ -35,7 +35,7 @@
 
 logger = logging.getLogger(__name__)
 
-__version__ = "0.5.9"
+__version__ = "0.5.10"
 
 
 def get_version_full():

diff --git a/packages/Qpyl/core/qfep.py b/packages/Qpyl/core/qfep.py
@@ -187,34 +187,41 @@ def calc_lra(self, lambda_a, lambda_b):
         lra = DataContainer(["E_type", "(E2-E1)_10", "(E2-E1)_01",
                              "LRA", "REORG"])
 
-        e1_st1, e1_st2, e2_st1, e2_st2 = None, None, None, None
+        e1_a, e1_b, e2_a, e2_b = None, None, None, None
         # get the appropriate rows of energies
+        # note that these energies are not scaled by lambda
         # [4:] ignores 'file', 'state', 'points' and 'lambda'
         for row in self.data_state[0].get_rows():
             if abs(row[3] - lambda_a) < 1e-7:
-                e1_st1 = row[4:]
+                e1_a = row[4:]
             if abs(row[3] - lambda_b) < 1e-7:
-                e1_st2 = row[4:]
+                e1_b = row[4:]
         # lambda2 in data_state[1] is actually (1-lambda), correct for that
         for row in self.data_state[1].get_rows():
             if abs((1 - row[3]) - lambda_a) < 1e-7:
-                e2_st1 = row[4:]
+                e2_a = row[4:]
             if abs((1 - row[3]) - lambda_b) < 1e-7:
-                e2_st2 = row[4:]
+                e2_b = row[4:]
 
-        if not e1_st1:
+        if not e1_a:
             raise QFepOutputError("LRA: No energy values for lambda == '{}'"
                                   "".format(lambda_a))
-        if not e1_st2:
+        if not e1_b:
             raise QFepOutputError("LRA: No energy values for lambda == '{}'"
                                   "".format(lambda_b))
 
+        la, lb = lambda_a, lambda_b
+        # calculate total E=(l1*E1 + l2*E2) energies
+        e1_state1 = [la*e1a + (1-la)*e2a for e1a, e2a in zip(e1_a, e2_a)]
+        e1_state2 = [la*e1b + (1-la)*e2b for e1b, e2b in zip(e1_b, e2_b)]
 
+        e2_state1 = [lb*e1a + (1-lb)*e2a for e1a, e2a in zip(e1_a, e2_a)]
+        e2_state2 = [lb*e1b + (1-lb)*e2b for e1b, e2b in zip(e1_b, e2_b)]
 
         # (E2-E1)_10    (reactant state) = First row E2 - E1
         # (E2-E1)_01    (products state) = Last row E2 - E1
-        des_st1 = [e2 - e1 for e1, e2 in zip(e1_st1, e2_st1)]
-        des_st2 = [e2 - e1 for e1, e2 in zip(e1_st2, e2_st2)]
+        des_st1 = [e2 - e1 for e1, e2 in zip(e1_state1, e2_state1)]
+        des_st2 = [e2 - e1 for e1, e2 in zip(e1_state2, e2_state2)]
 
         # LRA=0.5*(<E2-E1>_10+<E2-E1>_01)
         # REO=0.5*(<E2-E1>_10-<E2-E1>_01)

diff --git a/tests/Qpyl/core/qfep_test.py b/tests/Qpyl/core/qfep_test.py
@@ -55,8 +55,12 @@ def test_dgfep1(self, qfo1):
         assert is_close(qfo1.part1.dg, -10.204)
 
     def test_lra1(self, qfo1):
-        # LRA of EQbond
-        # pen&paper:
+        # LRA of EQbond (1.0 -> 0.0)
+        # pen&paper
+        # E1_10 = 1.0*-65.90 + 0.0*-8.21
+        # E2_10 = 0.0*-65.90 + 1.0*-8.21
+        # E1_01 = 1.0*-6.81  + 0.0*-77.97
+        # E2_01 = 0.0*-6.81  + 1.0*-77.97
         # E = 0.5*( <E2-E1>_10 + <E2-E1>_01 )
         # E = 0.5*( (-8.21 --65.90) + (-77.97 --6.81) )
         # E = -6.7350
@@ -73,6 +77,29 @@ def test_lra2(self, qfo2):
                    "43.62 6.75 30.10 -0.61 -4.43 0.29 0.00"
         assert lras == expected
 
+    def test_lra3(self, qfo1):
+        # LRA of EQbond (0.9 -> 0.6)
+        # pen&paper
+        # E1_10 = 0.9*-65.52 + 0.1*-11.18 = -60.086
+        # E2_10 = 0.6*-65.52 + 0.4*-11.18 = -43.784
+        # E1_01 = 0.9*-56.26 + 0.1*-28.54 = -53.488
+        # E2_01 = 0.6*-56.26 + 0.4*-28.54 = -45.172
+        # E = 0.5*( <E2-E1>_10 + <E2-E1>_01 )
+        # E = 0.5*( (-43.784 --60.086) + (-45.172 --53.488) )
+        # E = 12.309
+        lras = qfo1.part0.calc_lra(0.9, 0.6)
+        print lras
+        lra_eqbond = lras.get_columns(["LRA"])[0][1]
+        assert is_close(lra_eqbond, 12.309)
+
+    def test_lra4(self, qfo2):
+        lras = qfo2.part0.calc_lra(0.9, 0.6)
+        print lras
+        lras = " ".join(["{:.2f}".format(x) for x in lras.get_columns()[3]])
+        expected = "187.28 89.52 51.98 0.27 -0.00 30.78 14.74 19.51 14.87 "\
+                   "10.84 -0.20 0.42 0.07 0.00"
+        assert lras == expected
+
     def test_dgs1(self, qfo2):
         assert is_close(qfo2.part3.dga, 23.98)
         assert is_close(qfo2.part3.dg0, -21.98)

diff --git a/tests/qscripts-cli/q_analysefeps/1-qcp_excl/output/qaf.pd.json b/tests/qscripts-cli/q_analysefeps/1-qcp_excl/output/qaf.pd.json
@@ -141,19 +141,19 @@
           NaN
         ], 
         "ydata": [
-          756.74, 
-          364.9, 
-          214.7, 
-          0.92, 
-          -0.02, 
-          117.46, 
-          58.77, 
-          70.65, 
-          59.52000000000001, 
-          38.10000000000001, 
-          -0.6600000000000001, 
-          8.71, 
-          -0.07999999999999996, 
+          378.37, 
+          182.45000000000002, 
+          107.35, 
+          0.4600000000000002, 
+          -0.010000000000000002, 
+          58.73000000000002, 
+          29.385000000000005, 
+          35.325, 
+          29.760000000000005, 
+          19.05000000000001, 
+          -0.33000000000000185, 
+          4.355, 
+          -0.040000000000000036, 
           0.0
         ], 
         "xdata": [
@@ -201,19 +201,19 @@
           NaN
         ], 
         "ydata": [
-          166.33999999999997, 
-          71.47, 
-          24.99000000000001, 
-          -0.10999999999999999, 
-          0.0, 
-          63.06999999999999, 
-          6.92, 
-          43.89, 
-          7.23, 
-          22.50999999999999, 
-          -0.3000000000000007, 
-          -3.3200000000000003, 
-          -0.010000000000000231, 
+          83.16999999999996, 
+          35.73499999999996, 
+          12.494999999999997, 
+          -0.05499999999999994, 
+          0.0, 
+          31.534999999999968, 
+          3.459999999999999, 
+          21.945000000000004, 
+          3.615000000000002, 
+          11.254999999999995, 
+          -0.15000000000000213, 
+          -1.6599999999999993, 
+          -0.004999999999999449, 
           0.0
         ], 
         "xdata": [
@@ -261,19 +261,19 @@
           NaN
         ], 
         "ydata": [
-          461.53999999999996, 
-          218.185, 
-          119.845, 
-          0.405, 
-          -0.01, 
-          90.26499999999999, 
-          32.845, 
-          57.27, 
-          33.37500000000001, 
-          30.305, 
-          -0.4800000000000004, 
-          2.6950000000000003, 
-          -0.045000000000000095, 
+          230.76999999999998, 
+          109.09249999999999, 
+          59.9225, 
+          0.20250000000000012, 
+          -0.005000000000000001, 
+          45.13249999999999, 
+          16.422500000000003, 
+          28.635000000000005, 
+          16.687500000000004, 
+          15.152500000000003, 
+          -0.240000000000002, 
+          1.3475000000000006, 
+          -0.022499999999999742, 
           0.0
         ], 
         "xdata": [
@@ -321,19 +321,19 @@
           NaN
         ], 
         "ydata": [
-          295.20000000000005, 
-          146.71499999999997, 
-          94.85499999999999, 
-          0.515, 
-          -0.01, 
-          27.195, 
-          25.925, 
-          13.380000000000003, 
-          26.145000000000003, 
-          7.795000000000009, 
-          -0.17999999999999972, 
-          6.015000000000001, 
-          -0.034999999999999865, 
+          147.60000000000002, 
+          73.35750000000003, 
+          47.427499999999995, 
+          0.25750000000000006, 
+          -0.005000000000000001, 
+          13.597500000000025, 
+          12.962500000000002, 
+          6.6899999999999995, 
+          13.072500000000002, 
+          3.897500000000008, 
+          -0.08999999999999986, 
+          3.0075, 
+          -0.017500000000000293, 
           0.0
         ], 
         "xdata": [

diff --git a/tests/qscripts-cli/q_analysefeps/1-qcp_excl/output/random.dir_name/qaf.QCP.json b/tests/qscripts-cli/q_analysefeps/1-qcp_excl/output/random.dir_name/qaf.QCP.json
@@ -135,19 +135,19 @@
           NaN
         ], 
         "ydata": [
-          750.99, 
-          360.09000000000003, 
-          214.06, 
-          0.6200000000000001, 
-          -0.02, 
-          117.47, 
-          58.77, 
-          70.65, 
-          59.53, 
-          38.10000000000001, 
-          -0.6700000000000017, 
-          8.71, 
-          -0.09999999999999998, 
+          375.49500000000006, 
+          180.04500000000004, 
+          107.03, 
+          0.31000000000000005, 
+          -0.010000000000000002, 
+          58.735000000000014, 
+          29.385, 
+          35.325, 
+          29.765000000000008, 
+          19.05000000000001, 
+          -0.3349999999999973, 
+          4.355, 
+          -0.050000000000000044, 
           0.0
         ], 
         "xdata": [
@@ -195,19 +195,19 @@
           NaN
         ], 
         "ydata": [
-          168.35999999999999, 
-          76.72, 
-          24.939999999999998, 
-          -0.030000000000000027, 
-          0.0, 
-          63.06999999999999, 
-          3.6500000000000004, 
-          43.870000000000005, 
-          3.9800000000000004, 
-          22.50999999999999, 
-          -0.3000000000000007, 
-          -3.3200000000000003, 
-          -0.020000000000000018, 
+          84.17999999999995, 
+          38.359999999999985, 
+          12.469999999999999, 
+          -0.014999999999999902, 
+          0.0, 
+          31.534999999999968, 
+          1.8249999999999993, 
+          21.935, 
+          1.9899999999999984, 
+          11.254999999999995, 
+          -0.15000000000000213, 
+          -1.6599999999999993, 
+          -0.009999999999999787, 
           0.0
         ], 
         "xdata": [
@@ -255,19 +255,19 @@
           NaN
         ], 
         "ydata": [
-          459.675, 
-          218.40500000000003, 
-          119.5, 
-          0.29500000000000004, 
-          -0.01, 
-          90.27, 
-          31.21, 
-          57.260000000000005, 
-          31.755000000000003, 
-          30.305, 
-          -0.4850000000000012, 
-          2.6950000000000003, 
-          -0.06, 
+          229.8375, 
+          109.20250000000001, 
+          59.75, 
+          0.14750000000000008, 
+          -0.005000000000000001, 
+          45.13499999999999, 
+          15.605, 
+          28.630000000000003, 
+          15.877500000000003, 
+          15.152500000000003, 
+          -0.24249999999999972, 
+          1.3475000000000006, 
+          -0.029999999999999916, 
           0.0
         ], 
         "xdata": [
@@ -315,19 +315,19 @@
           NaN
         ], 
         "ydata": [
-          291.315, 
-          141.685, 
-          94.56, 
-          0.32500000000000007, 
-          -0.01, 
-          27.200000000000003, 
-          27.560000000000002, 
-          13.39, 
-          27.775, 
-          7.795000000000009, 
-          -0.1850000000000005, 
-          6.015000000000001, 
-          -0.03999999999999998, 
+          145.65750000000006, 
+          70.84250000000003, 
+          47.28, 
+          0.16249999999999998, 
+          -0.005000000000000001, 
+          13.600000000000023, 
+          13.780000000000001, 
+          6.695000000000002, 
+          13.887500000000005, 
+          3.897500000000008, 
+          -0.09249999999999758, 
+          3.0075, 
+          -0.02000000000000013, 
           0.0
         ], 
         "xdata": [