Add test for cpH and G-RxMC (ideal). (#33)

* Add test of G-RxMC (ideal).

* Adjust G-RxMC test

* Add test on cpH

* add test mode to speed up CI tests, get rid of numpy warnings, clean up sample scripts, change order of the CI tests

* Remove plot option from samples

* CI check

---------

Co-authored-by: Pablo M. Blanco <pblanco@dama.icp.uni-stuttgart.de>

Makefile

@@ -8,15 +8,14 @@ docs:
 	pdoc ./pyMBE.py -o ./documentation --docformat google
 
 tests:
-	python3 testsuite/henderson_hasselbalch_tests.py
 	python3 testsuite/lj_tests.py
 	python3 testsuite/generate_perpendicular_vectors_test.py
 	python3 testsuite/read-write-df_test.py
+	python3 testsuite/henderson_hasselbalch_tests.py
+	python3 testsuite/cph_ideal_tests.py
+	python3 testsuite/grxmc_ideal_tests.py
 	python3 testsuite/peptide_tests.py
 
-sample:
-	python3 sample_scripts/peptide_simulation_example.py
-
 visual:
 	python3 handy_scripts/vmd-traj.py
 	vmd -e visualization.tcl

pyMBE.py

@@ -486,7 +486,10 @@ def check_if_df_cell_has_a_value(self, index,key):
             `bool`: `True` if the cell has a value, `False` otherwise.
         """
         idx = self.pd.IndexSlice
-        return not self.pd.isna(self.df.loc[index, idx[key]])
+        import warnings
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            return not self.pd.isna(self.df.loc[index, idx[key]])
 
     def check_if_name_is_defined_in_df(self, name, pmb_type_to_be_defined):
         """

samples/branched_polyampholyte.py

@@ -1,20 +1,30 @@
-# Load espresso, sugar and other necessary libraries
+# Load espresso, pyMBE and other necessary libraries
 import sys
 import os 
 import inspect
 from matplotlib.style import use
 import espressomd
 import numpy as np
 import matplotlib.pyplot as plt
-from tqdm import tqdm
+import argparse
+import pandas as pd
 from espressomd.io.writer import vtf
 from espressomd import interactions
-
+import pyMBE
 
 # Create an instance of pyMBE library
-import pyMBE
+
 pmb = pyMBE.pymbe_library()
 
+# Command line arguments
+parser = argparse.ArgumentParser(description='Script that runs a Monte Carlo simulation of an ideal branched polyampholyte using pyMBE and ESPResSo.')
+parser.add_argument('--test', 
+                    default=False, 
+                    action='store_true',
+                    help='to run a short simulation for testing the script')
+args = parser.parse_args()
+
+
 # Load some functions from the handy_scripts library for convinience
 from lib.handy_functions import setup_langevin_dynamics
 from lib.analysis import block_analyze
@@ -26,14 +36,24 @@
 # Simulation parameters
 pmb.set_reduced_units(unit_length=0.4*pmb.units.nm)
 pH_range = np.linspace(2, 12, num=20)
-Samples_per_pH = 100
-MD_steps_per_sample = 0
+Samples_per_pH = 1000
+MD_steps_per_sample = 1000
 steps_eq = int(Samples_per_pH/3)
 N_samples_print = 10  # Write the trajectory every 100 samples
 probability_reaction =0.5 
 SEED = 100
 dt = 0.001
 solvent_permitivity = 78.3
+N_polyampholyte_chains = 5
+polyampholyte_concentration = 5.56e-4 *pmb.units.mol/pmb.units.L
+volume = N_polyampholyte_chains/(pmb.N_A*polyampholyte_concentration)
+
+if args.test:
+    Samples_per_pH = 100
+    probability_reaction = 1 
+    N_samples_print = 1000
+    N_polyampholyte_chains = 1 
+    pH_range = np.linspace(2, 12, num=10)
 
 # Define different particles
 # Inert particle 
@@ -96,9 +116,7 @@
 pmb.define_particle(name=anion_name,  q=-1, sigma=0.35*pmb.units.nm,  epsilon=1*pmb.units('reduced_energy'))
 
 # System parameters
-N_polyampholyte_chains = 5
-polyampholyte_concentration = 5.56e-4 *pmb.units.mol/pmb.units.L
-volume = N_polyampholyte_chains/(pmb.N_A*polyampholyte_concentration)
+
 L = volume ** (1./3.) # Side of the simulation box
 calculated_polyampholyte_concentration = N_polyampholyte_chains/(volume*pmb.N_A)
 
@@ -121,8 +139,8 @@
 total_ionisible_groups = len (list_ionisible_groups)
 
 print("The box length of your system is", L.to('reduced_length'), L.to('nm'))
-print('The polyampholyte concentration in your system is ', calculated_polyampholyte_concentration.to('mol/L') , 'with', N_polyampholyte_chains, 'peptides')
-print('The ionisable groups in your peptide are ', list_ionisible_groups)
+print('The polyampholyte concentration in your system is ', calculated_polyampholyte_concentration.to('mol/L') , 'with', N_polyampholyte_chains, 'molecules')
+print('The ionisable groups in your polyampholyte are ', list_ionisible_groups)
 
 RE, sucessfull_reactions_labels = pmb.setup_cpH(counter_ion=cation_name, constant_pH=2, SEED=SEED)
 print('The acid-base reaction has been sucessfully setup for ', sucessfull_reactions_labels)
@@ -148,17 +166,24 @@
 
 N_frame=0
 Z_pH=[] # List of the average global charge at each pH
+err_Z_pH=[] # List of the error of the global charge at each pH
+
 
 #Save the pyMBE dataframe in a CSV file
 pmb.write_pmb_df (filename='df.csv')
 
 # Main loop for performing simulations at different pH-values
+labels_obs=["time","charge"]
 
-for index in tqdm(range(len(pH_range))):
+for index in range(len(pH_range)):
 
     pH_value=pH_range[index]
-    # Sample list inicialization
-    Z_sim=[]
+
+    time_series={}
+
+    for label in labels_obs:
+        time_series[label]=[]
+
     RE.constant_pH = pH_value
 
     # Inner loop for sampling each pH value
@@ -174,36 +199,54 @@
             # Get polyampholyte net charge
             charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system, 
                     molecule_name="polyampholyte")      
-            Z_sim.append(charge_dict["mean"])
+            if args.test:
+                time_series["time"].append(step)
+            else:
+                time_series["time"].append(espresso_system.time)
+            time_series["charge"].append(charge_dict["mean"])
 
         if (step % N_samples_print == 0) :
             N_frame+=1
             with open('frames/trajectory'+str(N_frame)+'.vtf', mode='w+t') as coordinates:
                 vtf.writevsf(espresso_system, coordinates)
                 vtf.writevcf(espresso_system, coordinates)
 
-    Z_pH.append(Z_sim)
+    # Estimate the statistical error and the autocorrelation time of the data
+    processed_data = block_analyze(full_data=pd.DataFrame(time_series, columns=labels_obs))
+    Z_pH.append(processed_data["mean", "charge"])
+    err_Z_pH.append(processed_data["err_mean", "charge"])
     print("pH = {:6.4g} done".format(pH_value))
 
-# Estimate the statistical error and the autocorrelation time of the data
-
-print("Net charge analysis")
-av_net_charge, err_net_charge, tau_net_charge, block_size_net_charge = block_analyze(input_data=Z_pH)
-
-# Calculate the ideal titration curve of the polyampholyte with Henderson-Hasselbach equation (manually)
-pH_range_HH = np.linspace(2, 12, num=1000)
-Z_HH_manually = [10 * (1/(1+10**(pH_value-9)) - 1/(1+10**(4-pH_value))) for pH_value in pH_range_HH]
-
-# Calculate the ideal titration curve of the polyampholyte with Henderson-Hasselbach equation (pyMBE)
-Z_HH = pmb.calculate_HH(molecule_name="polyampholyte", 
-                        pH_list=pH_range_HH)
-
-fig, ax = plt.subplots(figsize=(10, 7))
-plt.errorbar(pH_range, av_net_charge, yerr=err_net_charge, fmt = 'o', capsize=3, label='Simulation')
-plt.plot(pH_range_HH, Z_HH_manually, label="Henderson-Hasselbalch (manually)")
-ax.plot(pH_range_HH, Z_HH, "-k", label='Henderson-Hasselbach (pyMBE)', linestyle="dashed")
-plt.legend()
-plt.xlabel('pH')
-plt.ylabel('Charge of the polyampholyte / e')
 
-plt.show()
+if args.test:
+     # Calculate the ideal titration curve of the polyampholyte with Henderson-Hasselbach equation (pyMBE)
+    Z_HH = pmb.calculate_HH(molecule_name="polyampholyte", 
+                            pH_list=pH_range)
+
+    # Write out the data
+    data = {}
+    data["Z_sim"] = np.asarray(Z_pH)
+    data["Z_HH"] = np.asarray(Z_HH)
+    data = pd.DataFrame.from_dict(data) 
+
+    data_path = pmb.get_resource(path="samples")
+    data.to_csv(f"{data_path}/data_polyampholyte_cph.csv", index=False)
+
+else:
+    # Calculate the ideal titration curve of the polyampholyte with Henderson-Hasselbach equation (manually)
+    pH_range_HH = np.linspace(2, 12, num=1000)
+    Z_HH_manually = [10 * (1/(1+10**(pH_value-9)) - 1/(1+10**(4-pH_value))) for pH_value in pH_range_HH]
+
+    # Calculate the ideal titration curve of the polyampholyte with Henderson-Hasselbach equation (pyMBE)
+    Z_HH = pmb.calculate_HH(molecule_name="polyampholyte", 
+                            pH_list=pH_range_HH)
+
+    fig, ax = plt.subplots(figsize=(10, 7))
+    plt.errorbar(pH_range, Z_pH, yerr=err_Z_pH, fmt = 'o', capsize=3, label='Simulation')
+    plt.plot(pH_range_HH, Z_HH_manually, label="Henderson-Hasselbalch (manually)")
+    ax.plot(pH_range_HH, Z_HH, "-k", label='Henderson-Hasselbach (pyMBE)')
+    plt.legend()
+    plt.xlabel('pH')
+    plt.ylabel('Charge of the polyampholyte / e')
+
+    plt.show()

samples/peptide.py

@@ -1,4 +1,4 @@
-# Load espresso, sugar and other necessary libraries
+# Load espresso, pyMBE and other necessary libraries
 import sys
 import os 
 import inspect
@@ -10,9 +10,9 @@
 from espressomd.io.writer import vtf
 from espressomd import interactions
 from espressomd import electrostatics
+import pyMBE
 
 # Create an instance of pyMBE library
-import pyMBE
 pmb = pyMBE.pymbe_library()
 
 # Load some functions from the handy_scripts library for convinience