From 04723421445b4c3bf1b23f13a5f3403fe3d62ab5 Mon Sep 17 00:00:00 2001 From: "Pablo M. Blanco" <75744061+pm-blanco@users.noreply.github.com> Date: Fri, 14 Jun 2024 16:44:46 +0200 Subject: [PATCH] change default of acidity to None (#69) * change default of acidity to None * remove call to warnings, improve docs * remove verbose from test * solve merge conflicts --- pyMBE.py | 50 ++++---- samples/Beyer2024/globular_protein.py | 4 +- samples/branched_polyampholyte.py | 1 - samples/peptide_mixture_grxmc_ideal.py | 130 ++++++++++++++++----- testsuite/create_molecule_position_test.py | 2 - testsuite/read-write-df_test.py | 3 +- testsuite/set_particle_acidity_test.py | 14 ++- 7 files changed, 141 insertions(+), 63 deletions(-) diff --git a/pyMBE.py b/pyMBE.py index 36cad88..9120892 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -32,12 +32,12 @@ class pymbe_library(): The library for the Molecular Builder for ESPResSo (pyMBE) Attributes: - N_A(`obj`): Avogadro number using the `pmb.units` UnitRegistry. - Kb(`obj`): Boltzmann constant using the `pmb.units` UnitRegistry. - e(`obj`): Elemental charge using the `pmb.units` UnitRegistry. - df(`obj`): PandasDataframe used to bookkeep all the information stored in pyMBE. Typically refered as `pmb.df`. - kT(`obj`): Thermal energy using the `pmb.units` UnitRegistry. It is used as the unit of reduced energy. - Kw(`obj`): Ionic product of water using the `pmb.units` UnitRegistry. Used in the setup of the G-RxMC method. + N_A(`pint.Quantity`): Avogadro number. + Kb(`pint.Quantity`): Boltzmann constant. + e(`pint.Quantity`): Elementary charge. + df(`Pandas.Dataframe`): Dataframe used to bookkeep all the information stored in pyMBE. Typically refered as `pmb.df`. + kT(`pint.Quantity`): Thermal energy. + Kw(`pint.Quantity`): Ionic product of water. Used in the setup of the G-RxMC method. """ units = pint.UnitRegistry() N_A=6.02214076e23 / units.mol @@ -1496,14 +1496,14 @@ def define_particle_entry_in_df(self,name): self.df.at [index,'pmb_type'] = 'particle' return index - def define_particle(self, name, z=0, acidity='inert', pka=None, sigma=None, epsilon=None, cutoff=None, offset=None,verbose=True,overwrite=False): + def define_particle(self, name, z=0, acidity=None, pka=None, sigma=None, epsilon=None, cutoff=None, offset=None,verbose=True,overwrite=False): """ Defines the properties of a particle object. Args: name(`str`): Unique label that identifies this particle type. z(`int`, optional): Permanent charge number of this particle type. Defaults to 0. - acidity(`str`, optional): Identifies whether if the particle is `acidic` or `basic`, used to setup constant pH simulations. Defaults to `inert`. + acidity(`str`, optional): Identifies whether if the particle is `acidic` or `basic`, used to setup constant pH simulations. Defaults to None. pka(`float`, optional): If `particle` is an acid or a base, it defines its pka-value. Defaults to None. sigma(`pint.Quantity`, optional): Sigma parameter used to set up Lennard-Jones interactions for this particle type. Defaults to None. cutoff(`pint.Quantity`, optional): Cutoff parameter used to set up Lennard-Jones interactions for this particle type. Defaults to None. @@ -2226,10 +2226,7 @@ def load_interaction_parameters(self, filename, verbose=False, overwrite=False): elif not_required_key in without_units: not_required_attributes[not_required_key]=param_dict.pop(not_required_key) else: - if not_required_key == 'acidity': - not_required_attributes[not_required_key] = 'inert' - else: - not_required_attributes[not_required_key]=None + not_required_attributes[not_required_key]=None self.define_particle(name=param_dict.pop('name'), z=not_required_attributes.pop('z'), sigma=not_required_attributes.pop('sigma'), @@ -2591,25 +2588,34 @@ def search_particles_in_residue(self, residue_name): return list_of_particles_in_residue - def set_particle_acidity(self, name, acidity='inert', default_charge_number=0, pka=None, verbose=True, overwrite=True): + def set_particle_acidity(self, name, acidity=None, default_charge_number=0, pka=None, verbose=True, overwrite=True): """ - Sets the particle acidity if it is acidic or basic, creates `state_one` and `state_two` with the protonated and - deprotonated states. In each state is set: `label`,`charge` and `es_type`. If it is inert, it will define only `state_one`. + Sets the particle acidity including the charges in each of its possible states. Args: name(`str`): Unique label that identifies the `particle`. - acidity(`str`): Identifies whether the particle is `acidic` or `basic`, used to setup constant pH simulations. Defaults to `inert`. + acidity(`str`): Identifies whether the particle is `acidic` or `basic`, used to setup constant pH simulations. Defaults to None. default_charge_number (`int`): Charge number of the particle. Defaults to 0. pka(`float`, optional): If `particle` is an acid or a base, it defines its pka-value. Defaults to None. verbose(`bool`, optional): Switch to activate/deactivate verbose. Defaults to True. overwrite(`bool`, optional): Switch to enable overwriting of already existing values in pmb.df. Defaults to False. + + Note: + - For particles with `acidity = acidic` or `acidity = basic`, `state_one` and `state_two` correspond to the protonated and + deprotonated states, respectively. + - For particles without an acidity `acidity = None`, only `state_one` is defined. + - Each state has the following properties as sub-indexes: `label`,`charge` and `es_type`. """ acidity_valid_keys = ['inert','acidic', 'basic'] - if acidity not in acidity_valid_keys: - raise ValueError(f"Acidity {acidity} provided for particle name {name} is not supproted. Valid keys are: {acidity_valid_keys}") - if acidity in ['acidic', 'basic'] and pka is None: - raise ValueError(f"pKa not provided for particle with name {name} with acidity {acidity}. pKa must be provided for acidic or basic particles.") - + if acidity is not None: + if acidity not in acidity_valid_keys: + raise ValueError(f"Acidity {acidity} provided for particle name {name} is not supproted. Valid keys are: {acidity_valid_keys}") + if acidity in ['acidic', 'basic'] and pka is None: + raise ValueError(f"pKa not provided for particle with name {name} with acidity {acidity}. pKa must be provided for acidic or basic particles.") + if acidity == "inert": + acidity = None + print("Deprecation warning: the keyword 'inert' for acidity has been replaced by setting acidity = None. For backwards compatibility, acidity has been set to None. Support for `acidity = 'inert'` may be deprecated in future releases of pyMBE") + self.define_particle_entry_in_df(name=name) for index in self.df[self.df['name']==name].index: @@ -2631,7 +2637,7 @@ def set_particle_acidity(self, name, acidity='inert', default_charge_number=0, p new_value=self.propose_unused_type(), verbose=verbose, overwrite=overwrite) - if self.df.loc [self.df['name'] == name].acidity.iloc[0] == 'inert': + if pd.isna(self.df.loc [self.df['name'] == name].acidity.iloc[0]): self.add_value_to_df(key=('state_one','z'), index=index, new_value=default_charge_number, diff --git a/samples/Beyer2024/globular_protein.py b/samples/Beyer2024/globular_protein.py index 38b9e76..6f170e0 100644 --- a/samples/Beyer2024/globular_protein.py +++ b/samples/Beyer2024/globular_protein.py @@ -258,7 +258,7 @@ net_charge_amino_save = {} AA_label_list=[] for amino in net_charge_residues.keys(): - amino_part_row=pmb.df[(pmb.df['residue_id']== amino) & (pmb.df['acidity'] != "inert")] + amino_part_row=pmb.df[(pmb.df['residue_id']== amino) & ((pmb.df['acidity'] == "acidic") | (pmb.df['acidity'] == "basic"))] if not amino_part_row.empty: label = f'charge_{amino_part_row["name"].values[0]}' if label not in AA_label_list: @@ -279,7 +279,7 @@ charge_residues_per_type[label]=[] for amino in charge_residues.keys(): - amino_part_row=pmb.df[(pmb.df['residue_id']== amino) & (pmb.df['acidity'] != "inert")] + amino_part_row=pmb.df[(pmb.df['residue_id']== amino) & ((pmb.df['acidity'] == "acidic") | (pmb.df['acidity'] == "basic"))] if not amino_part_row.empty: label = f'charge_{amino_part_row["name"].values[0]}' if label in AA_label_list: diff --git a/samples/branched_polyampholyte.py b/samples/branched_polyampholyte.py index 6bd5e3b..f515150 100644 --- a/samples/branched_polyampholyte.py +++ b/samples/branched_polyampholyte.py @@ -73,7 +73,6 @@ # Inert particle pmb.define_particle( name = "I", - acidity = "inert", z = 0, sigma = 1*pmb.units('reduced_length'), epsilon = 1*pmb.units('reduced_energy')) diff --git a/samples/peptide_mixture_grxmc_ideal.py b/samples/peptide_mixture_grxmc_ideal.py index 2df3701..1148f3a 100644 --- a/samples/peptide_mixture_grxmc_ideal.py +++ b/samples/peptide_mixture_grxmc_ideal.py @@ -29,10 +29,6 @@ # Create an instance of pyMBE library pmb = pyMBE.pymbe_library(seed=42) -import warnings - - - # Command line arguments valid_modes=["standard", "unified"] @@ -59,7 +55,8 @@ from lib.analysis import block_analyze # Simulation parameters -pmb.set_reduced_units(unit_length=0.4*pmb.units.nm, Kw=1e-14) +pmb.set_reduced_units(unit_length=0.4*pmb.units.nm, + Kw=1e-14) pH_range = np.linspace(2, 12, num=20) Samples_per_pH = 500 MD_steps_per_sample = 0 @@ -79,6 +76,7 @@ sequence2 = 'nEEHHc' pep2_concentration = 1e-2 *pmb.units.mol/pmb.units.L N_peptide2_chains = 10 +verbose=True if args.test: Samples_per_pH = 1000 @@ -87,6 +85,7 @@ N_peptide1_chains = 5 N_peptide2_chains = 5 pH_range = np.linspace(2, 12, num=10) + verbose = False # Load peptide parametrization from Lunkad, R. et al. Molecular Systems Design & Engineering (2021), 6(2), 122-131. @@ -96,9 +95,7 @@ path_to_interactions=pmb.get_resource("parameters/peptides/Lunkad2021.json") pmb.load_interaction_parameters(filename=path_to_interactions) -with warnings.catch_warnings(): - warnings.simplefilter('error') - pmb.load_pka_set(path_to_pka) +pmb.load_pka_set(path_to_pka) # Defines the bonds @@ -112,13 +109,18 @@ } -pmb.define_default_bond(bond_type = bond_type, bond_parameters = harmonic_bond) +pmb.define_default_bond(bond_type = bond_type, + bond_parameters = harmonic_bond) # Defines the peptides in the pyMBE data frame peptide1 = 'generic_peptide1' -pmb.define_peptide (name=peptide1, sequence=sequence1, model=model) +pmb.define_peptide (name=peptide1, + sequence=sequence1, + model=model) peptide2 = 'generic_peptide2' -pmb.define_peptide (name=peptide2, sequence=sequence2, model=model) +pmb.define_peptide (name=peptide2, + sequence=sequence2, + model=model) # Solution parameters c_salt=5e-3 * pmb.units.mol/ pmb.units.L @@ -129,17 +131,35 @@ sodium_name = 'Na' chloride_name = 'Cl' - pmb.define_particle(name=proton_name, z=1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) - pmb.define_particle(name=hydroxide_name, z=-1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) - pmb.define_particle(name=sodium_name, z=1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) - pmb.define_particle(name=chloride_name, z=-1, sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=proton_name, + z=1, + sigma=0.35*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=hydroxide_name, + z=-1, + sigma=0.35*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=sodium_name, + z=1, + sigma=0.35*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) + pmb.define_particle(name=chloride_name, + z=-1, + sigma=0.35*pmb.units.nm, + epsilon=1*pmb.units('reduced_energy')) elif args.mode == 'unified': cation_name = 'Na' anion_name = 'Cl' - 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')) + 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 @@ -154,19 +174,40 @@ pmb.add_bonds_to_espresso(espresso_system=espresso_system) # Create your molecules into the espresso system -pmb.create_pmb_object(name=peptide1, number_of_objects= N_peptide1_chains,espresso_system=espresso_system, use_default_bond=True) -pmb.create_pmb_object(name=peptide2, number_of_objects= N_peptide2_chains,espresso_system=espresso_system, use_default_bond=True) +pmb.create_pmb_object(name=peptide1, + number_of_objects= N_peptide1_chains, + espresso_system=espresso_system, + use_default_bond=True) +pmb.create_pmb_object(name=peptide2, + number_of_objects= N_peptide2_chains, + espresso_system=espresso_system, + use_default_bond=True) if args.mode == 'standard': - pmb.create_counterions(object_name=peptide1,cation_name=proton_name,anion_name=hydroxide_name,espresso_system=espresso_system) # Create counterions for the peptide chains - pmb.create_counterions(object_name=peptide2,cation_name=proton_name,anion_name=hydroxide_name,espresso_system=espresso_system) # Create counterions for the peptide chains + pmb.create_counterions(object_name=peptide1, + cation_name=proton_name, + anion_name=hydroxide_name, + espresso_system=espresso_system) # Create counterions for the peptide chains + pmb.create_counterions(object_name=peptide2, + cation_name=proton_name, + anion_name=hydroxide_name, + espresso_system=espresso_system) # Create counterions for the peptide chains c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system,cation_name=sodium_name,anion_name=chloride_name,c_salt=c_salt) elif args.mode == 'unified': - pmb.create_counterions(object_name=peptide1, cation_name=cation_name,anion_name=anion_name,espresso_system=espresso_system) # Create counterions for the peptide chains - pmb.create_counterions(object_name=peptide2, cation_name=cation_name,anion_name=anion_name,espresso_system=espresso_system) # Create counterions for the peptide chains + pmb.create_counterions(object_name=peptide1, + cation_name=cation_name, + anion_name=anion_name, + espresso_system=espresso_system) # Create counterions for the peptide chains + pmb.create_counterions(object_name=peptide2, + cation_name=cation_name, + anion_name=anion_name, + espresso_system=espresso_system) # Create counterions for the peptide chains - c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system,cation_name=cation_name,anion_name=anion_name,c_salt=c_salt) + c_salt_calculated = pmb.create_added_salt(espresso_system=espresso_system, + cation_name=cation_name, + anion_name=anion_name, + c_salt=c_salt) with open('frames/trajectory0.vtf', mode='w+t') as coordinates: @@ -182,10 +223,20 @@ print("The box length of your system is", L.to('reduced_length'), L.to('nm')) if args.mode == 'standard': - grxmc, 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, activity_coefficient=lambda x: 1.0) + grxmc, 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, + activity_coefficient=lambda x: 1.0) elif args.mode == 'unified': - grxmc, labels, ionic_strength_res = pmb.setup_grxmc_unified(pH_res=2, c_salt_res=c_salt, cation_name=cation_name, anion_name=anion_name, activity_coefficient=lambda x: 1.0) -print('The acid-base reaction has been sucessfully setup for ', labels) + grxmc, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_unified(pH_res=2, + c_salt_res=c_salt, + cation_name=cation_name, + anion_name=anion_name, + activity_coefficient=lambda x: 1.0) +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 type_map =pmb.get_type_map() @@ -231,9 +282,19 @@ time_series[label]=[] if args.mode == 'standard': - grxmc, 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, 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) elif args.mode == 'unified': - grxmc, 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, 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) # Inner loop for sampling each pH value @@ -268,7 +329,8 @@ vtf.writevcf(espresso_system, coordinates) # Estimate the statistical error and the autocorrelation time of the data - processed_data = block_analyze(full_data=pd.DataFrame(time_series, columns=labels_obs)) + processed_data = block_analyze(full_data=pd.DataFrame(time_series, columns=labels_obs), + verbose=verbose) Z_pH.append(processed_data["mean", "charge"]) err_Z_pH.append(processed_data["err_mean", "charge"]) @@ -281,7 +343,9 @@ if args.test: # Calculate the ideal titration curve of the peptide with Henderson-Hasselbach equation - HH_charge_dict = pmb.calculate_HH_Donnan(c_macro={peptide1: pep1_concentration, peptide2: pep2_concentration}, c_salt=c_salt, pH_list=pH_range) + HH_charge_dict = pmb.calculate_HH_Donnan(c_macro={peptide1: pep1_concentration, peptide2: pep2_concentration}, + c_salt=c_salt, + pH_list=pH_range) Z_HH_Donnan = HH_charge_dict["charges_dict"] xi = HH_charge_dict["partition_coefficients"] @@ -299,7 +363,9 @@ else: # Calculate the ideal titration curve of the peptide with Henderson-Hasselbach equation pH_range_HH = np.linspace(2, 12, num=100) - HH_charge_dict = pmb.calculate_HH_Donnan(c_macro={peptide1: pep1_concentration, peptide2: pep2_concentration}, c_salt=c_salt, pH_list=pH_range_HH) + HH_charge_dict = pmb.calculate_HH_Donnan(c_macro={peptide1: pep1_concentration, peptide2: pep2_concentration}, + c_salt=c_salt, + pH_list=pH_range_HH) Z_HH_Donnan = HH_charge_dict["charges_dict"] xi = HH_charge_dict["partition_coefficients"] diff --git a/testsuite/create_molecule_position_test.py b/testsuite/create_molecule_position_test.py index 07469e4..7c13b4a 100644 --- a/testsuite/create_molecule_position_test.py +++ b/testsuite/create_molecule_position_test.py @@ -34,11 +34,9 @@ pos_list = [[10,10,10], [20,20,20], [30,30,30]] pmb.define_particle(name='central_mon', - acidity='inert', sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) pmb.define_particle(name='side_mon', - acidity='inert', sigma=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) diff --git a/testsuite/read-write-df_test.py b/testsuite/read-write-df_test.py index 846f3eb..d493816 100644 --- a/testsuite/read-write-df_test.py +++ b/testsuite/read-write-df_test.py @@ -37,8 +37,7 @@ name = "I", sigma = 0.3*pmb.units.nm, epsilon = 1*pmb.units('reduced_energy'), - z = 0, - acidity = "inert") + z = 0) pmb.define_particle( name = "A", diff --git a/testsuite/set_particle_acidity_test.py b/testsuite/set_particle_acidity_test.py index 2eb9293..42d4785 100644 --- a/testsuite/set_particle_acidity_test.py +++ b/testsuite/set_particle_acidity_test.py @@ -33,7 +33,7 @@ def check_acid_base_setup(input_parameters,acidity_setup): """ pmb.define_particle(**input_parameters) - if input_parameters["acidity"] == "inert": + if input_parameters["acidity"] is None: input_parameters.pop("z") # Checks that the input parameters are stored properly for parameter_key in input_parameters.keys(): @@ -55,7 +55,7 @@ def check_acid_base_setup(input_parameters,acidity_setup): # Clean pmb.df pmb.setup_df() input_parameters={"name":"I", - "acidity": "inert", + "acidity": None, "pka": np.nan, "z":2} acidity_setup={"state_one":{"label":f"{input_parameters['name']}", @@ -66,6 +66,15 @@ def check_acid_base_setup(input_parameters,acidity_setup): check_acid_base_setup(input_parameters=input_parameters, acidity_setup=acidity_setup) +print("*** Unit test passed ***") +print("*** Unit test: check that a deprecation warning is raised if the keyword 'inert' is used for acidity ***") +# Clean pmb.df +pmb.setup_df() +input_parameters={"name":"I", + "acidity": "inert", + "pka": np.nan, + "z":2} +pmb.define_particle(**input_parameters) print("*** Unit test passed ***") print("*** Unit test: check that all acid/base input parameters in define_particle for an acid are correctly stored in pmb.df***") # Clean pmb.df @@ -95,6 +104,7 @@ def check_acid_base_setup(input_parameters,acidity_setup): check_acid_base_setup(input_parameters=input_parameters, acidity_setup=acidity_setup) print("*** Unit test passed ***") + print("*** Unit test: check that set_particle_acidity raises a ValueError if pKa is not provided and pKa is acidic or basic ***") input_parametersA={"name":"A", "acidity": "acidic" }