merge conflict solved

.github/workflows/testsuite.yml

@@ -0,0 +1,32 @@
+name: run tests
+
+on:
+  push:
+  pull_request:
+
+permissions:
+  contents: read # to fetch code (actions/checkout)
+
+jobs:
+  ubuntu:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Setup EESSI
+        uses: eessi/github-action-eessi@v3
+        with:
+          eessi_stack_version: "2023.06"
+      - name: Checkout repository
+        uses: actions/checkout@main
+      - name: Install dependencies
+        run: |
+          module load ESPResSo/4.2.1-foss-2023a
+          python3 -m pip install --user -r requirements.txt
+      - name: Run testsuite
+        run: |
+          module load ESPResSo/4.2.1-foss-2023a
+          export OLD_PYTHONPATH="${PYTHONPATH}"
+          export PYTHONPATH="$(realpath .)${PYTHONPATH:+:$PYTHONPATH}"
+          sed -i "s/\${ESPResSo_build_path}\///" Makefile
+          make testsuite
+          export PYTHONPATH="${OLD_PYTHONPATH}"
+        shell: bash

CONTRIBUTING.md

@@ -1,16 +1,16 @@
 # Please contribute to pyMBE!
 
 ## Bug reporting
-If you discover any strange feature or unexpected behaviour of the library, please report it by opening an issue in this Gitlab repository or contacting Dr. Pablo M. Blanco (pablb@ntnu.no) or any other member of the pyMBE development team.
-Once a ticket is open in Gitlab, we will work to fix the issue as soon as possible.
+If you discover any strange feature or unexpected behaviour of the library, please report it by opening an issue in this GitHub repository or contacting Dr. Pablo M. Blanco (pablb@ntnu.no) or any other member of the pyMBE development team.
+Once a ticket is open in GitHub, we will work to fix the issue as soon as possible.
 
 ## New contributors
 New developers are welcomed to contribute to extend the functionalities of the pyMBE library. 
-To contribute to the pyMBE library, first one needs to be added as a member of this Gitlab repository.
+To contribute to the pyMBE library, first one needs to be added as a member of this GitHub repository.
 If you want to contribute to the development of the pyMBE library, please contact Dr. Pablo M. Blanco (pablb@ntnu.no)
 
 ## Rules to contribute to pyMBE
-All new features of pyMBE must be developed in a branch of this Gitlab repository and must be reported to a member of the pyMBE development team.
 Any new version of the code must reproduce all the data stored in `reference_data/`. 
 The scripts provided in `tests/` can be used to quickly set-up  simulations with pyMBE to reproduce such data sets.
 All new code will be reviewed by at least one member of the pyMBE development team before being merged into the stable version of the library to ensure that a functional version of the code is always available.
+Class methods are sorted in alphabetical order.

Makefile

@@ -3,13 +3,17 @@
 .PHONY: visual 
 .PHONY: clean
 .PHONY: tests
+.PHONY: testsuite
 .PHONY: docs
 
 ESPResSo_build_path=~/software/espresso_v4.2/build/
 
 docs:
 	pdoc ./pyMBE.py -o ./docs --docformat google 
 
+testsuite:
+	${ESPResSo_build_path}/pypresso testsuite/LYS_ASP_peptide.py
+
 sample:
 	${ESPResSo_build_path}/pypresso sample_scripts/peptide_simulation_example.py
 

README.md

@@ -9,6 +9,7 @@ pyMBE provides tools to facilitate building up molecules with complex architectu
 - [Pandas](https://pandas.pydata.org/) v1.5.3
 - [Pint-Pandas](https://pypi.org/project/Pint-Pandas/) v0.5
 - [Numpy](https://numpy.org/)
+- [SciPy](https://scipy.org/)
 - [pdoc](https://pdoc.dev/) (for building the docs)
 
 ## Branches
@@ -76,3 +77,16 @@ make tutorial
 
 provided that you have modified the `$ESPResSo_build_path` variable in `Makefile` to match the path where you have built ESPResSo.
 
+### Run the testsuite
+
+To make sure your code is valid, please run the testsuite before submitting your contribution:
+
+```sh
+PYTHONPATH=$(realpath .) make testsuite
+```
+
+When contributing new features, consider adding a unit test in the `testsuite/`
+folder and a corresponding line in the `testsuite` target of the Makefile.
+
+Every contribution is automatically tested in CI using EESSI (https://www.eessi.io)
+and the [EESSI GitHub Action](https://github.com/marketplace/actions/eessi).

requirements.txt

@@ -0,0 +1,8 @@
+numpy>=1.23
+pandas>=1.5.3
+pint>=0.20.01
+pint-pandas==0.5
+biopandas==0.5.1.dev0
+scipy
+matplotlib
+tqdm

sample_scripts/peptide_mixture_grxmc_ideal.py

@@ -37,7 +37,7 @@
 steps_eq = int(Samples_per_pH)
 N_samples_print = 1000  # Write the trajectory every 100 samples
 probability_reaction =1
-SEED = 100
+SEED = 42
 dt = 0.001
 solvent_permitivity = 78.3
 
@@ -137,7 +137,7 @@
 
 print("The box length of your system is", L.to('reduced_length'), L.to('nm'))
 
-RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=2, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, sodium_name=sodium_name, chloride_name=chloride_name, SEED=SEED)
+RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=2, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
 print('The acid-base reaction has been sucessfully setup for ', sucessful_reactions_labels)
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
@@ -184,7 +184,7 @@
     Z_sim=[]
     num_plus=[]
 
-    RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions (pH_res=pH_value, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, sodium_name=sodium_name, chloride_name=chloride_name, SEED=SEED)
+    RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions (pH_res=pH_value, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
 
     # Inner loop for sampling each pH value
 

sample_scripts/peptide_mixture_grxmc_unified_ideal.py

@@ -37,7 +37,7 @@
 steps_eq = int(Samples_per_pH)
 N_samples_print = 1000  # Write the trajectory every 100 samples
 probability_reaction =1
-SEED = 100
+SEED = 42
 dt = 0.001
 solvent_permitivity = 78.3
 

tests/simulation_script_grxmc.py

@@ -108,9 +108,9 @@
 # Set up the reactions
 ionic_strength, excess_chemical_potential_monovalent_pairs_in_bulk_data, bjerrums, excess_chemical_potential_monovalent_pairs_in_bulk_data_error =np.loadtxt("../../../../../../../reference_data/excess_chemical_potential.dat", unpack=True)
 excess_chemical_potential_monovalent_pair_interpolated = interpolate.interp1d(ionic_strength, excess_chemical_potential_monovalent_pairs_in_bulk_data)
-excess_chemical_potential_monovalent_pair = lambda x: excess_chemical_potential_monovalent_pair_interpolated(x.to('1/(reduced_length**3 * N_A)').magnitude) * pmb.units('reduced_energy')
+activity_coefficient_monovalent_pair = lambda x: np.exp(excess_chemical_potential_monovalent_pair_interpolated(x.to('1/(reduced_length**3 * N_A)').magnitude))
 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, sodium_name=sodium_name, chloride_name=chloride_name, SEED=REACTION_SEED, excess_chemical_potential_monovalent_pair=excess_chemical_potential_monovalent_pair, pka_set=pka_set)
+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, SEED=REACTION_SEED, activity_coefficient=activity_coefficient_monovalent_pair, pka_set=pka_set)
 print('The acid-base reaction has been sucessfully set up for ', sucessful_reactions_labels)
 
 # Setup espresso to track the ionization of the acid groups

testsuite/LYS_ASP_peptide.py

@@ -0,0 +1,199 @@
+# Load espresso, pyMBE and other necessary libraries
+import os
+import sys
+import inspect
+import numpy as np
+import pandas as pd
+from tqdm import tqdm
+import espressomd
+from espressomd import interactions
+from espressomd.io.writer import vtf
+from espressomd import electrostatics 
+
+# Create an instance of pyMBE library
+import pyMBE
+pmb = pyMBE.pymbe_library()
+
+# Load some functions from the handy_scripts library for convinience
+from handy_scripts.handy_functions import setup_electrostatic_interactions
+from handy_scripts.handy_functions import minimize_espresso_system_energy
+from handy_scripts.handy_functions import setup_langevin_dynamics
+from handy_scripts.handy_functions import block_analyze
+
+# The trajectories of the simulations will be stored using espresso built-up functions in separed files in the folder 'frames'
+if not os.path.exists('./frames'):
+    os.makedirs('./frames')
+
+# Simulation parameters
+pmb.set_reduced_units(unit_length=0.4*pmb.units.nm)
+pH_range = np.linspace(2, 12, num=21)
+Samples_per_pH = 36
+MD_steps_per_sample = 50
+steps_eq =int(Samples_per_pH/3)
+N_samples_print = 10 # Write the trajectory every 100 samples
+probability_reaction = 0.5 
+SEED = 1
+dt = 0.01
+solvent_permitivity = 78.3 
+
+L = 25.513*pmb.units.nm
+
+# Peptide parameters
+N_aminoacids = 5
+sequence = 'K'*N_aminoacids+'D'*N_aminoacids
+model = '2beadAA'  # Model with 2 beads per each aminoacid
+pep_concentration = 1e-4 *pmb.units.mol/pmb.units.L
+
+# Solution parameters
+cation_name = 'Na'
+anion_name = 'Cl'
+c_salt = 1e-2 * pmb.units.mol/ pmb.units.L
+
+# Define salt parameters
+
+pmb.define_particle( name=cation_name,  q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
+pmb.define_particle( name=anion_name,  q=-1, diameter=0.35*pmb.units.nm,  epsilon=1*pmb.units('reduced_energy'))
+
+# Load peptide parametrization from Lunkad, R. et al.  Molecular Systems Design & Engineering (2021), 6(2), 122-131.
+
+pmb.load_interaction_parameters (filename=pmb.get_resource('reference_parameters/interaction_parameters/Lunkad2021.txt'))
+pmb.load_pka_set (filename=pmb.get_resource('reference_parameters/pka_sets/CRC1991.txt'))
+
+# Define the peptide on the pyMBE dataframe 
+pmb.define_peptide( name=sequence, sequence=sequence, model=model)
+
+# System parameters
+volume = L**3
+N_peptide_chains = int ( volume * pmb.N_A * pep_concentration)
+L = volume ** (1./3.) # Side of the simulation box
+calculated_peptide_concentration = N_peptide_chains/(volume*pmb.N_A)
+
+# Create an instance of an espresso system
+espresso_system = espressomd.System(box_l=[L.to('reduced_length').magnitude]*3)
+
+# Add all bonds to espresso system
+pmb.add_bonds_to_espresso (espresso_system=espresso_system)
+
+# Create your molecules into the espresso system
+
+pmb.create_pmb_object (name=sequence, number_of_objects= N_peptide_chains,espresso_system=espresso_system, use_default_bond=True)
+
+
+# Create counterions for the peptide chains
+pmb.create_counterions (object_name=sequence,cation_name=cation_name,anion_name=anion_name,espresso_system=espresso_system) 
+c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system,cation_name=cation_name,anion_name=anion_name,c_salt=c_salt)
+
+
+#List of ionisible groups 
+basic_groups = pmb.df.loc[(~pmb.df['particle_id'].isna()) & (pmb.df['acidity']=='basic')].name.to_list()
+acidic_groups = pmb.df.loc[(~pmb.df['particle_id'].isna()) & (pmb.df['acidity']=='acidic')].name.to_list()
+list_ionisible_groups = basic_groups + acidic_groups
+total_ionisible_groups = len (list_ionisible_groups)
+
+print("The box length of your system is", L.to('reduced_length'), L.to('nm'))
+print('The peptide concentration in your system is ', calculated_peptide_concentration.to('mol/L') , 'with', N_peptide_chains, 'peptides')
+print('The ionisable groups in your peptide are ', list_ionisible_groups)
+
+# Setup the acid-base reactions of the peptide using the constant pH ensemble
+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)
+
+# Setup espresso to track the each type defined in type_map
+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)
+print('The non interacting type is set to ', non_interacting_type)
+
+# Setup the potential energy
+pmb.setup_lj_interactions(espresso_system=espresso_system)
+minimize_espresso_system_energy (espresso_system=espresso_system)
+setup_electrostatic_interactions(units=pmb.units,
+                                            espresso_system=espresso_system,
+                                            kT=pmb.kT)
+minimize_espresso_system_energy (espresso_system=espresso_system)
+
+
+setup_langevin_dynamics (espresso_system=espresso_system, 
+                                    kT = pmb.kT, 
+                                    SEED = SEED,
+                                    time_step=dt,
+                                    tune_skin=False)
+
+espresso_system.cell_system.skin=0.4
+
+# Save the initial state
+with open('frames/trajectory1.vtf', mode='w+t') as coordinates:
+    vtf.writevsf(espresso_system, coordinates)
+    vtf.writevcf(espresso_system, coordinates)
+
+N_frame=0
+Z_pH=[] # List of the average global charge at each pH
+Rg_pH=[] 
+
+particle_id_list = pmb.get_particle_id_map(object_name=sequence)["all"]
+first_peptide_id = min(particle_id_list)
+
+#Save the pyMBE dataframe in a CSV file
+pmb.df.to_csv('df.csv')
+
+for index in (pbar := tqdm(range(len(pH_range)))):
+    # Sample list inicialization
+    pH_value=pH_range[index]
+    Z_sim=[]
+    Rg_sim=[]
+    Z_groups_time_series=[]       
+    RE.constant_pH = pH_value
+    pbar.set_description(f"pH = {pH_value:2.1f}")
+
+    # Inner loop for sampling each pH value
+    for step in range(Samples_per_pH+steps_eq):
+
+        if np.random.random() > probability_reaction:
+            espresso_system.integrator.run(steps=MD_steps_per_sample)
+        else:
+            RE.reaction(reaction_steps=total_ionisible_groups)
+
+        if ( step > steps_eq):
+            # Get peptide net charge
+            charge_dict=pmb.calculate_net_charge(espresso_system=espresso_system, 
+                                                molecule_name=sequence)      
+            Z_sim.append(charge_dict["mean"])
+            # Get peptide radius of gyration
+            Rg = espresso_system.analysis.calc_rg(chain_start=first_peptide_id, number_of_chains=N_peptide_chains, chain_length=len(particle_id_list))
+            Rg_value = pmb.units.Quantity(Rg[0], 'reduced_length')
+            Rg_nm = Rg_value.to('nm').magnitude
+            Rg_sim.append(Rg_nm)
+
+        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(np.array(Z_sim))
+    Rg_pH.append(Rg_sim)
+
+print("Net charge analysis")
+av_charge, err_charge, tau_charge, block_size = block_analyze(input_data=pmb.np.array(Z_pH))
+
+print("Rg analysis")
+av_rg, err_rg, tau_rg, block_size = block_analyze(input_data=Rg_pH)
+
+# Calculate the ideal titration curve of the peptide with Henderson-Hasselbach equation
+Z_HH = pmb.calculate_HH(object_name=sequence,
+                         pH_list=pH_range)
+
+# Load the reference data 
+reference_file_Path = pmb.get_resource("reference_data/Lys-AspMSDE.csv")
+reference_data = pd.read_csv(reference_file_Path)
+
+Z_ref = N_aminoacids*-1*reference_data['aaa']+N_aminoacids*reference_data['aab']         
+Rg_ref = reference_data['arg']*0.37
+
+np.testing.assert_allclose(np.copy(av_charge), Z_ref.to_numpy(), atol=2.5, rtol=0.)