fix general samples

Author: pmblanco

pyMBE/lib/handy_functions.py

@@ -260,7 +260,6 @@ def define_peptide_AA_residues(sequence,model, pmb):
 
         - Residue names are constructed as `"AA-<residue>"`, e.g., `"AA-A"`, `"AA-L"`.
     """
-    defined_residues = []
     for residue_name in sequence:
         if model == '1beadAA':
             central_bead = residue_name
@@ -273,11 +272,15 @@ def define_peptide_AA_residues(sequence,model, pmb):
                 central_bead = 'CA'              
                 side_chains = [residue_name]
         residue_name='AA-'+residue_name
-        if residue_name not in defined_residues:   
+        if "residue" in pmb.db._templates:
+            if residue_name not in pmb.db._templates["residue"]:   
+                pmb.define_residue(name = residue_name, 
+                                    central_bead = central_bead,
+                                    side_chains = side_chains)
+        else:
             pmb.define_residue(name = residue_name, 
-                                central_bead = central_bead,
-                                side_chains = side_chains)              
-            defined_residues.append(residue_name)
+                                    central_bead = central_bead,
+                                    side_chains = side_chains)
 
 def do_reaction(algorithm, steps):
     """

samples/branched_polyampholyte.py

@@ -162,11 +162,10 @@
                                           anion_name=anion_name,
                                           c_salt=c_salt)
 
-#List of ionisable groups
-# collect ionisable particles*
-acidbase_particles = ["A","B"]
+# count acid/base particles
+pka_set = pmb.get_pka_set()
 acid_base_ids = []
-for name in acidbase_particles:
+for name in pka_set.keys():
     acid_base_ids+=pmb.db.find_instance_ids_by_name(pmb_type="particle",
                                                     name=name)        
 total_ionisable_groups = len(acid_base_ids)

samples/peptide_cpH.py

@@ -29,10 +29,7 @@
 pmb = pyMBE.pymbe_library(seed=42)
 
 # Load some functions from the handy_scripts library for convenience
-from pyMBE.lib.handy_functions import setup_electrostatic_interactions
-from pyMBE.lib.handy_functions import relax_espresso_system
-from pyMBE.lib.handy_functions import setup_langevin_dynamics
-from pyMBE.lib.handy_functions import do_reaction
+from pyMBE.lib.handy_functions import setup_electrostatic_interactions, relax_espresso_system, setup_langevin_dynamics, do_reaction, define_peptide_AA_residues
 from pyMBE.lib.analysis import built_output_name
 
 parser = argparse.ArgumentParser(description='Sample script to run the pre-made peptide models with pyMBE')
@@ -97,7 +94,13 @@
 path_to_interactions=pmb.root / "parameters" / "peptides" / "Lunkad2021"
 path_to_pka=pmb.root / "parameters" / "pka_sets" / "Hass2015.json"
 pmb.load_database(folder=path_to_interactions) 
-pmb.load_pka_set (path_to_pka)
+pmb.load_pka_set(path_to_pka)
+
+# Define acid/base particle states
+pka_set = pmb.get_pka_set()
+for particle_name in pka_set.keys():
+    pmb.define_monoprototic_particle_states(particle_name=particle_name,
+                                            acidity=pka_set[particle_name]["acidity"])
 
 generic_bond_length=0.4 * pmb.units.nm
 generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')
@@ -110,8 +113,11 @@
                         bond_parameters = HARMONIC_parameters)
 
 
-# Defines the peptide in the pyMBE data frame
+# Defines the peptide in the pyMBE database
 peptide_name = 'generic_peptide'
+define_peptide_AA_residues(sequence=sequence,
+                           model="2beadAA",
+                           pmb=pmb)
 pmb.define_peptide (name=peptide_name, 
                     sequence=sequence, 
                     model=model)
@@ -125,6 +131,8 @@
                     sigma=0.35*pmb.units.nm,  
                     epsilon=1*pmb.units('reduced_energy'))
 
+
+
 # Create an instance of an espresso system
 espresso_system=espressomd.System (box_l = [L.to('reduced_length').magnitude]*3)
 espresso_system.time_step=dt
@@ -152,20 +160,28 @@
     vtf.writevsf(espresso_system, coordinates)
     vtf.writevcf(espresso_system, coordinates)
 
-#List of ionisable 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_ionisable_groups = basic_groups + acidic_groups
-total_ionisable_groups = len(list_ionisable_groups)
+# count acid/base particles
+pka_set = pmb.get_pka_set()
+acid_base_ids = []
+list_ionisable_groups = []
+for name in pka_set.keys():
+    part_ids = pmb.db.find_instance_ids_by_name(pmb_type="particle",
+                                                name=name)
+    if part_ids:
+        acid_base_ids+=part_ids
+        list_ionisable_groups+=[name]  
+total_ionisable_groups = len(acid_base_ids)
 
 if verbose:
     print(f"The box length of your system is {L.to('reduced_length')} {L.to('nm')}")
     print(f"The peptide concentration in your system is {calculated_peptide_concentration.to('mol/L')} with {N_peptide_chains} peptides")
     print(f"The ionisable groups in your peptide are {list_ionisable_groups}")
 
-cpH, labels = pmb.setup_cpH(counter_ion=cation_name, constant_pH=pH_value)
+cpH = pmb.setup_cpH(counter_ion=cation_name, 
+                    constant_pH=pH_value)
 if verbose:
-    print(f"The acid-base reaction has been successfully setup for {labels}")
+    print("The acid-base reaction has been successfully set up for:")
+    print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
 type_map =pmb.get_type_map()

samples/peptide_mixture_grxmc_ideal.py

@@ -24,7 +24,7 @@
 from espressomd.io.writer import vtf
 import pyMBE
 from pyMBE.lib.analysis import built_output_name
-from pyMBE.lib.handy_functions import do_reaction
+from pyMBE.lib.handy_functions import do_reaction, define_peptide_AA_residues
 
 # Create an instance of pyMBE library
 pmb = pyMBE.pymbe_library(seed=42)
@@ -95,11 +95,27 @@
 # Note that this parameterization only includes some of the natural aminoacids
 # For the other aminoacids the user needs to use  a parametrization including all the aminoacids in the peptide sequence
 path_to_pka=pmb.root / "parameters" / "pka_sets" / "Hass2015.json"
-path_to_interactions=pmb.root / "parameters" / "peptides" / "Lunkad2021.json"
-
-pmb.load_interaction_parameters(filename=path_to_interactions) 
+path_to_interactions=pmb.root / "parameters" / "peptides" / "Lunkad2021"
+
+pmb.load_database(folder=path_to_interactions) 
+# define templates for the c and n ends
+pmb.define_particle(name="n",
+                    sigma=1*pmb.units.reduced_length,
+                    epsilon=1*pmb.units.reduced_energy,
+                    acidity="basic")
+pmb.define_particle(name="c",
+                    sigma=1*pmb.units.reduced_length,
+                    epsilon=1*pmb.units.reduced_energy,
+                    acidity="acidic")
 pmb.load_pka_set(path_to_pka)
 
+# Define acid/base particle states
+pka_set = pmb.get_pka_set()
+for particle_name in pka_set.keys():
+    if particle_name not in ["c", "n"]: # Avoid redefing the ends
+        pmb.define_monoprototic_particle_states(particle_name=particle_name,
+                                                acidity=pka_set[particle_name]["acidity"])
+    
 # Defines the bonds
 bond_type = 'harmonic'
 generic_bond_length=0.4 * pmb.units.nm
@@ -120,6 +136,12 @@
 pmb.define_peptide (name=peptide2, 
                     sequence=sequence2, 
                     model=model)
+define_peptide_AA_residues(sequence=sequence1,
+                           model=model,
+                           pmb=pmb)
+define_peptide_AA_residues(sequence=sequence2,
+                           model=model,
+                           pmb=pmb)
 
 # Solution parameters
 c_salt=5e-3 * pmb.units.mol/ pmb.units.L
@@ -212,31 +234,35 @@
     vtf.writevsf(espresso_system, coordinates)
     vtf.writevcf(espresso_system, coordinates)
 
-#List of ionisable 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_ionisable_groups = basic_groups + acidic_groups
-total_ionisable_groups = len (list_ionisable_groups)
+# count acid/base particles
+pka_set = pmb.get_pka_set()
+acid_base_ids = []
+for name in pka_set.keys():
+    acid_base_ids+=pmb.db.find_instance_ids_by_name(pmb_type="particle",
+                                                    name=name)        
+total_ionisable_groups = len(acid_base_ids)
+
 # Get peptide net charge
 if verbose:
     print("The box length of your system is", L.to('reduced_length'), L.to('nm'))
 
 if args.mode == 'standard':
-    grxmc, 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,
-                                                                                   activity_coefficient=lambda x: 1.0)
+    grxmc, 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,
+                                                        activity_coefficient=lambda x: 1.0)
 elif args.mode == 'unified':
-    grxmc, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_unified(pH_res=pH_value, 
-                                                                                 c_salt_res=c_salt, 
-                                                                                 cation_name=cation_name, 
-                                                                                 anion_name=anion_name,
-                                                                                 activity_coefficient=lambda x: 1.0)
+    grxmc, ionic_strength_res = pmb.setup_grxmc_unified(pH_res=pH_value, 
+                                                        c_salt_res=c_salt, 
+                                                        cation_name=cation_name, 
+                                                        anion_name=anion_name,
+                                                        activity_coefficient=lambda x: 1.0)
 if verbose:
-    print('The acid-base reaction has been sucessfully setup for ', sucessful_reactions_labels)
+    print("The acid-base reaction has been successfully set up for:")
+    print(pmb.get_reactions_df())
 
 # Setup espresso to track the ionization of the acid/basic groups in peptide
 type_map =pmb.get_type_map()

samples/plot_peptide_cpH.py

@@ -23,6 +23,7 @@
 import argparse
 import pathlib
 import pandas as pd
+import pyMBE.lib.handy_functions as hf
 # Create an instance of pyMBE library
 import pyMBE
 pmb = pyMBE.pymbe_library(seed=42)
@@ -52,15 +53,18 @@
 
 # Define peptide parameters
 sequence = args.sequence
-# Define the peptide in the pyMBE dataframe and load the pka set
+# Define the peptide in the pyMBE database and load the pka set
 # This is necessary to calculate the analytical solution from the Henderson-Hasselbach equation
 peptide = 'generic_peptide'
-pmb.define_peptide (name=peptide, 
-                    sequence=sequence,
-                   model="1beadAA") # not really relevant for plotting
+model="1beadAA" # not really relevant for plotting
+pmb.define_peptide(name=peptide, 
+                   sequence=sequence,
+                   model=model) 
 path_to_pka=pmb.root / "parameters" / "pka_sets" / "Hass2015.json"
 pmb.load_pka_set(path_to_pka)
-
+hf.define_peptide_AA_residues(sequence=sequence,
+                              model=model,
+                              pmb=pmb)
 # Calculate the ideal titration curve of the peptide with Henderson-Hasselbach equation
 if args.mode == "plot":
     pH_range_HH = np.linspace(2, 12, num=100)

samples/plot_peptide_mixture_grxmc_ideal.py

@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2024 pyMBE-dev team
+# Copyright (C) 2024-2026 pyMBE-dev team
 #
 # This file is part of pyMBE.
 #
@@ -22,6 +22,7 @@
 import pandas as pd
 import argparse
 import pathlib
+from pyMBE.lib.handy_functions import define_peptide_AA_residues
 
 parser = argparse.ArgumentParser(description='Plots the titration data from peptide.py and the corresponding analytical solution.')
 parser.add_argument('--sequence1',
@@ -65,6 +66,8 @@
 # Create an instance of pyMBE library
 pmb = pyMBE.pymbe_library(seed=42)
 c_salt=args.csalt * pmb.units.mol/ pmb.units.L
+model = "1beadAA"
+
 
 # Define peptide parameters
 sequence1 = args.sequence1
@@ -75,6 +78,13 @@
 
 # Define the peptides in the pyMBE data frame and load the pka set
 # This is necesary to calculate the analytical solution from the Henderson-Hasselbach equation
+
+define_peptide_AA_residues(sequence=sequence1,
+                           model=model,
+                           pmb=pmb)
+define_peptide_AA_residues(sequence=sequence2,
+                           model=model,
+                           pmb=pmb)
 peptide1 = 'generic_peptide1'
 pmb.define_peptide (name=peptide1, 
                     sequence=sequence1,
@@ -86,6 +96,7 @@
 path_to_pka=pmb.root / "parameters" / "pka_sets" / "Hass2015.json"
 pmb.load_pka_set(path_to_pka)
 
+
 # Calculate the ideal titration curve of the peptide with Henderson-Hasselbach equation
 if args.mode == "plot":
     pH_range_HH = np.linspace(2, 12, num=100)