solve inconsistencies between charge and charge number and other vari…

…ous small changes (#72)

* Make names of RE objects more descriptive

* Change SEED to seed

* Make the acitivity coefficient mandatory

* Fixed some typos

* Change charge to charge number

* Continue changing charge to charge number

* Add test for get_charge_number_map_tests

* Add units to calculate_net_charge

* Add generate_coordinates_tests

* Remove unused import

* Increase code coverage

Makefile

@@ -53,6 +53,8 @@ unit_tests:
 	${PYTHON} testsuite/setup_salt_ions_unit_tests.py
 	${PYTHON} testsuite/globular_protein_unit_tests.py
 	${PYTHON} testsuite/analysis_tests.py
+	${PYTHON} testsuite/charge_number_map_tests.py
+	${PYTHON} testsuite/generate_coordinates_tests.py
 
 functional_tests:
 	${PYTHON} testsuite/cph_ideal_tests.py

maintainer/standarize_data.py

@@ -23,7 +23,7 @@
 import argparse 
 
 # Create an instance of pyMBE library
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 # Expected inputs
 supported_filenames=["data_landsgesell.csv",

parameters/peptides/Blanco2021.json

@@ -9,14 +9,14 @@
 		"D": {"object_type":"particle", "name": "D", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"E": {"object_type":"particle", "name": "E", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"n": {"object_type":"particle", "name": "n", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
-		"S": {"object_type":"particle", "name": "S", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
+		"S": {"object_type":"particle", "name": "S", "z":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"H": {"object_type":"particle", "name": "H", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
-		"A": {"object_type":"particle", "name": "A", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
+		"A": {"object_type":"particle", "name": "A", "z":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"K": {"object_type":"particle", "name": "K", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"Y": {"object_type":"particle", "name": "Y", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"R": {"object_type":"particle", "name": "R", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
-		"G": {"object_type":"particle", "name": "G", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
-		"F": {"object_type":"particle", "name": "F", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
+		"G": {"object_type":"particle", "name": "G", "z":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
+		"F": {"object_type":"particle", "name": "F", "z":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"c": {"object_type":"particle", "name": "c", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"bond": {"object_type":"bond",  "bond_type": "harmonic", "bond_parameters" : {"r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N/m"}}, "particle_pairs": [["n","D"],["S","D"],["S","H"],["H","A"],["A","K"],["E","H"],["E","K"],["K","R"],["K","H"],["R","H"],["H","H"],["H","G"],["G","Y"],["Y","K"],["R","K"],["K","F"],["H","S"],["H","F"],["H","R"],["R","G"],["Y","G"],["Y","c"]]}
 	}

parameters/peptides/Lunkad2021.json

@@ -6,7 +6,7 @@
 		"citekey": "lunkad2021a"
 	},
 	"data": {
-		"CA": {"object_type":"particle", "name": "CA", "q":0, "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
+		"CA": {"object_type":"particle", "name": "CA", "z":0, "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"D": {"object_type":"particle", "name": "D",  "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"E": {"object_type":"particle", "name": "E",  "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},
 		"H": {"object_type":"particle", "name": "H",  "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}},

samples/Beyer2024/create_paper_data.py

@@ -26,7 +26,7 @@
 import subprocess
 
 # Create an instance of pyMBE library
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 valid_fig_labels=["7a", "7b", "7c", "8a", "8b", "9"]
 valid_modes=["short-run","long-run", "test"]

samples/Beyer2024/globular_protein.py

@@ -25,7 +25,7 @@
 from espressomd.io.writer import vtf
 # Create an instance of pyMBE library
 import pyMBE
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 #Import functions from handy_functions script 
 from lib.handy_functions import setup_electrostatic_interactions
@@ -143,13 +143,13 @@
 anion_name = 'Cl'
 
 pmb.define_particle(name = cation_name, 
-                    q = 1, 
+                    z = 1, 
                     sigma=sigma, 
                     epsilon=epsilon,
                     offset=ion_size-sigma)
 
 pmb.define_particle(name = anion_name,  
-                    q =-1, 
+                    z =-1, 
                     sigma=sigma, 
                     epsilon=epsilon,
                     offset=ion_size-sigma)
@@ -197,18 +197,18 @@
     print ('The total amount of ionizable groups are:',total_ionisible_groups)
 
 #Setup of the reactions in espresso 
-RE, sucessfull_reactions_labels = pmb.setup_cpH(counter_ion=cation_name, 
-                                                constant_pH= pH_value)
+cpH, labels = pmb.setup_cpH(counter_ion=cation_name, 
+                            constant_pH= pH_value)
 if verbose:
-    print('The acid-base reaction has been sucessfully setup for ', sucessfull_reactions_labels)
+    print('The acid-base reaction has been sucessfully setup for ', labels)
 
 type_map = pmb.get_type_map()
 types = list (type_map.values())
 espresso_system.setup_type_map( type_list = types)
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
 non_interacting_type = max(type_map.values())+1
-RE.set_non_interacting_type (type=non_interacting_type)
+cpH.set_non_interacting_type (type=non_interacting_type)
 if verbose:
     print('The non interacting type is set to ', non_interacting_type)
 
@@ -251,7 +251,8 @@
     time_series[label]=[]
 
 charge_dict=pmb.calculate_net_charge (espresso_system=espresso_system, 
-                                            molecule_name=protein_name)
+                                            molecule_name=protein_name,
+                                            dimensionless=True)
 
 net_charge_residues = charge_dict ['residues']
 net_charge_amino_save = {}
@@ -267,9 +268,10 @@
 
 for step in tqdm(range(N_samples),disable=not verbose):      
     espresso_system.integrator.run (steps = integ_steps)
-    RE.reaction( reaction_steps = total_ionisible_groups)
+    cpH.reaction(reaction_steps = total_ionisible_groups)
     charge_dict=pmb.calculate_net_charge (espresso_system=espresso_system, 
-                                            molecule_name=protein_name)
+                                            molecule_name=protein_name,
+                                            dimensionless=True)
     charge_residues = charge_dict['residues']
     charge_residues_per_type={}
 

samples/Beyer2024/peptide.py

@@ -29,7 +29,7 @@
 from lib import handy_functions as hf
 
 # Create an instance of pyMBE library
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 valid_modes=["short-run","long-run", "test"]
 parser = argparse.ArgumentParser(description='Script to run the peptide test cases for pyMBE')
@@ -129,13 +129,13 @@
 c_salt=5e-3 * pmb.units.mol/ pmb.units.L
 
 pmb.define_particle(name=cation_name,
-                    q=1,
+                    z=1,
                     sigma=sigma,
                     epsilon=1*pmb.units('reduced_energy'),
                     offset=offset_cation)
 
 pmb.define_particle(name=anion_name,
-                    q=-1,
+                    z=-1,
                     sigma=sigma,
                     epsilon=1*pmb.units('reduced_energy'),
                     offset=offset_anion)
@@ -168,20 +168,20 @@
                      c_salt=c_salt,
                     verbose=verbose)
 
-RE, successful_reactions_labels = pmb.setup_cpH(counter_ion=cation_name,
+cpH, labels = pmb.setup_cpH(counter_ion=cation_name,
                                                 constant_pH=pH)
 
 if verbose:
     print(f"The box length of your system is {L.to('reduced_length')} = {L.to('nm')}")
-    print(f"The acid-base reaction has been successfully setup for {successful_reactions_labels}")
+    print(f"The acid-base reaction has been successfully setup for {labels}")
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
 type_map =pmb.get_type_map()
 espresso_system.setup_type_map(type_list = list(type_map.values()))
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
 non_interacting_type = max(type_map.values())+1
-RE.set_non_interacting_type (type=non_interacting_type)
+cpH.set_non_interacting_type (type=non_interacting_type)
 if verbose:
     print(f"The non-interacting type is set to {non_interacting_type}")
 
@@ -218,10 +218,11 @@
     # Run LD
     espresso_system.integrator.run(steps=MD_steps_per_sample)
     # Run MC
-    RE.reaction(reaction_steps=len(sequence))
+    cpH.reaction(reaction_steps=len(sequence))
     # Sample observables
     charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system,
-                                            molecule_name=sequence)
+                                            molecule_name=sequence,
+                                            dimensionless=True)
 
     Rg = espresso_system.analysis.calc_rg(chain_start=0,
                                         number_of_chains=N_peptide_chains,

samples/Beyer2024/weak_polyelectrolyte_dialysis.py

@@ -34,7 +34,7 @@
 from lib import analysis
 
 # Create an instance of pyMBE library
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 # Load some functions from the handy_scripts library for convinience
 from lib.handy_functions import setup_electrostatic_interactions
@@ -125,10 +125,10 @@
 sodium_name = 'Na'
 chloride_name = 'Cl'
 
-pmb.define_particle(name=proton_name, q=1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
-pmb.define_particle(name=hydroxide_name,  q=-1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
-pmb.define_particle(name=sodium_name, q=1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
-pmb.define_particle(name=chloride_name,  q=-1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=proton_name, z=1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=hydroxide_name,  z=-1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=sodium_name, z=1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=chloride_name, z=-1, sigma=1*pmb.units('reduced_length'), epsilon=1*pmb.units('reduced_energy'))
 
 # System parameters (some are read from the command line)
 c_mon_sys = args.c_mon_sys * pmb.units.mol/ pmb.units.L
@@ -169,9 +169,9 @@
 activity_coefficient_monovalent_pair = lambda x: np.exp(excess_chemical_potential_monovalent_pair_interpolated(x.to('1/(reduced_length**3 * N_A)').magnitude))
 if verbose:
     print("Setting up reactions...")
-RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=pH_res, c_salt_res=c_salt_res, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, activity_coefficient=activity_coefficient_monovalent_pair, pka_set=pka_set)
+grxmc, labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=pH_res, c_salt_res=c_salt_res, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, activity_coefficient=activity_coefficient_monovalent_pair, pka_set=pka_set)
 if verbose:
-    print('The acid-base reaction has been sucessfully set up for ', sucessful_reactions_labels)
+    print('The acid-base reaction has been sucessfully set up for ', labels)
 
 # Setup espresso to track the ionization of the acid groups
 type_map = pmb.get_type_map()
@@ -180,7 +180,7 @@
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
 non_interacting_type = max(type_map.values())+1
-RE.set_non_interacting_type (type=non_interacting_type)
+grxmc.set_non_interacting_type (type=non_interacting_type)
 
 #Set up the interactions
 pmb.setup_lj_interactions(espresso_system=espresso_system)
@@ -196,7 +196,7 @@
     print("Running warmup without electrostatics")
 for i in tqdm(range(100),disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    RE.reaction(reaction_steps=1000)
+    grxmc.reaction(reaction_steps=1000)
 
 setup_electrostatic_interactions(units=pmb.units,
                                 espresso_system=espresso_system,
@@ -218,7 +218,7 @@
     N_warmup_loops = 100
 for i in tqdm(range(N_warmup_loops),disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    RE.reaction(reaction_steps=100)
+    grxmc.reaction(reaction_steps=100)
 
 
 # Main loop
@@ -236,13 +236,13 @@
     N_production_loops = 100
 for i in tqdm(range(N_production_loops),disable=not verbose):
     espresso_system.integrator.run(steps=1000)
-    RE.reaction(reaction_steps=100)
+    grxmc.reaction(reaction_steps=100)
 
     # Measure time
     time_series["time"].append(espresso_system.time)
 
     # Measure degree of ionization
-    charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system, molecule_name=polyacid_name)      
+    charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system, molecule_name=polyacid_name, dimensionless=True)
     time_series["alpha"].append(np.abs(charge_dict["mean"])/Chain_length)
 
 data_path = args.output

samples/branched_polyampholyte.py

@@ -28,7 +28,7 @@
 
 # Create an instance of pyMBE library
 
-pmb = pyMBE.pymbe_library(SEED=42)
+pmb = pyMBE.pymbe_library(seed=42)
 
 # Command line arguments
 parser = argparse.ArgumentParser(description='Script that runs a Monte Carlo simulation of an ideal branched polyampholyte using pyMBE and ESPResSo.')
@@ -74,7 +74,7 @@
 pmb.define_particle(
     name = "I",
     acidity = "inert",
-    q = 0,
+    z = 0,
     sigma = 1*pmb.units('reduced_length'),
     epsilon = 1*pmb.units('reduced_energy'))
 
@@ -127,8 +127,8 @@
 anion_name = 'Cl'
 c_salt=5e-3 * pmb.units.mol/ pmb.units.L
 
-pmb.define_particle(name=cation_name, q=1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
-pmb.define_particle(name=anion_name,  q=-1, sigma=0.35*pmb.units.nm,  epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=cation_name, z=1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle(name=anion_name, z=-1, sigma=0.35*pmb.units.nm,  epsilon=1*pmb.units('reduced_energy'))
 
 # System parameters
 
@@ -157,8 +157,8 @@
 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)
-print('The acid-base reaction has been sucessfully setup for ', sucessfull_reactions_labels)
+cpH, labels = pmb.setup_cpH(counter_ion=cation_name, constant_pH=2)
+print('The acid-base reaction has been sucessfully setup for ', labels)
 
 # Setup espresso to track the ionization of the acid/basic groups 
 type_map = pmb.get_type_map()
@@ -167,7 +167,7 @@
 
 # Setup the non-interacting type for speeding up the sampling of the reactions
 non_interacting_type = max(type_map.values())+1
-RE.set_non_interacting_type (type=non_interacting_type)
+cpH.set_non_interacting_type (type=non_interacting_type)
 print('The non interacting type is set to ', non_interacting_type)
 
 #Setup Langevin
@@ -197,7 +197,7 @@
     for label in labels_obs:
         time_series[label]=[]
 
-    RE.constant_pH = pH_value
+    cpH.constant_pH = pH_value
 
     # Inner loop for sampling each pH value
 
@@ -206,12 +206,13 @@
         if np.random.random() > probability_reaction:
             espresso_system.integrator.run(steps=MD_steps_per_sample)        
         else:
-            RE.reaction( reaction_steps = total_ionisible_groups)
+            cpH.reaction( reaction_steps = total_ionisible_groups)
 
         if step > steps_eq:
             # Get polyampholyte net charge
             charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system, 
-                    molecule_name="polyampholyte")      
+                    molecule_name="polyampholyte",
+                    dimensionless=True)
             if args.test:
                 time_series["time"].append(step)
             else: