diff --git a/Makefile b/Makefile index 6a33fcc..53a9b18 100644 --- a/Makefile +++ b/Makefile @@ -10,6 +10,7 @@ docs: tests: python3 testsuite/lj_tests.py python3 testsuite/set_particle_acidity_test.py + python3 testsuite/bond_tests.py python3 testsuite/generate_perpendicular_vectors_test.py python3 testsuite/create_molecule_position_test.py python3 testsuite/read-write-df_test.py diff --git a/README.md b/README.md index 40b49df..5ec731f 100644 --- a/README.md +++ b/README.md @@ -150,7 +150,7 @@ To make sure your code is valid, please run the testsuite before submitting your ```sh source pymbe/bin/activate -make testsuite +make tests deactivate ``` diff --git a/parameters/peptides/Avogadro_parametrization.txt b/parameters/peptides/Avogadro_parametrization.txt index 23f37d6..3d86852 100644 --- a/parameters/peptides/Avogadro_parametrization.txt +++ b/parameters/peptides/Avogadro_parametrization.txt @@ -8,12 +8,12 @@ {"object_type":"particle", "name": "n", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "H", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "K", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} -{"object_type":"bond", "name1": "CA", "name2": "CA", "bond_type": "harmonic", "r_0": {"value":0.382, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "A", "bond_type": "harmonic", "r_0": {"value":0.153, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "c", "bond_type": "harmonic", "r_0": {"value":0.246, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "D", "bond_type": "harmonic", "r_0": {"value":0.329, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "E", "bond_type": "harmonic", "r_0": {"value":0.435, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "Y", "bond_type": "harmonic", "r_0": {"value":0.648, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "n", "bond_type": "harmonic", "r_0": {"value":0.146, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.452, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.558, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","CA"]] , "bond_parameters" : {"r_0": {"value":0.382, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","D"]] , "bond_parameters" : {"r_0": {"value":0.329, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","E"]] , "bond_parameters" : {"r_0": {"value":0.435, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","H"]] , "bond_parameters" : {"r_0": {"value":0.452, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","Y"]] , "bond_parameters" : {"r_0": {"value":0.648, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","K"]] , "bond_parameters" : {"r_0": {"value":0.558, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","A"]] , "bond_parameters" : {"r_0": {"value":0.153, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","c"]] , "bond_parameters" : {"r_0": {"value":0.246, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","n"]] , "bond_parameters" : {"r_0": {"value":0.146, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} \ No newline at end of file diff --git a/parameters/peptides/Blanco2020.txt b/parameters/peptides/Blanco2020.txt index 379fd42..9987956 100644 --- a/parameters/peptides/Blanco2020.txt +++ b/parameters/peptides/Blanco2020.txt @@ -11,25 +11,4 @@ {"object_type":"particle", "name": "G", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "F", "q":0, "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "c", "sigma": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} -{"object_type":"bond", "name1": "n", "name2": "D", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "S", "name2": "D", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "S", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "A", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "A", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "E", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "E", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "K", "name2": "R", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "K", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "R", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "G", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "G", "name2": "Y", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "Y", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "R", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "K", "name2": "F", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "S", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "F", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "H", "name2": "R", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "R", "name2": "G", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "Y", "name2": "G", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} -{"object_type":"bond", "name1": "Y", "name2": "c", "bond_type": "harmonic", "r_0": {"value":0.4, "units":"nm"}, "k": {"value": 0.41, "units":"N / m"}} +{"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"]]} \ No newline at end of file diff --git a/parameters/peptides/Lunkad2021.txt b/parameters/peptides/Lunkad2021.txt index b302743..93edd5d 100644 --- a/parameters/peptides/Lunkad2021.txt +++ b/parameters/peptides/Lunkad2021.txt @@ -5,9 +5,9 @@ {"object_type":"particle", "name": "H", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "Y", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} {"object_type":"particle", "name": "K", "sigma": {"value":0.35, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}} -{"object_type":"bond", "name1": "CA", "name2": "CA", "bond_type": "harmonic", "r_0": {"value":0.382, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "D", "bond_type": "harmonic", "r_0": {"value":0.329, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "E", "bond_type": "harmonic", "r_0": {"value":0.435, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "H", "bond_type": "harmonic", "r_0": {"value":0.452, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "Y", "bond_type": "harmonic", "r_0": {"value":0.648, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} -{"object_type":"bond", "name1": "CA", "name2": "K", "bond_type": "harmonic", "r_0": {"value":0.558, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","CA"]] , "bond_parameters" : {"r_0": {"value":0.382, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","D"]] , "bond_parameters" : {"r_0": {"value":0.329, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","E"]] , "bond_parameters" : {"r_0": {"value":0.435, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","H"]] , "bond_parameters" : {"r_0": {"value":0.452, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","Y"]] , "bond_parameters" : {"r_0": {"value":0.648, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} +{"object_type":"bond", "bond_type": "harmonic", "particle_pairs": [["CA","K"]] , "bond_parameters" : {"r_0": {"value":0.558, "units":"nm"}, "k": {"value": 400, "units":"reduced_energy / nm**2"}}} \ No newline at end of file diff --git a/pyMBE.py b/pyMBE.py index 91d7be9..116e75d 100644 --- a/pyMBE.py +++ b/pyMBE.py @@ -678,6 +678,67 @@ def create_added_salt (self, espresso_system, cation_name, anion_name, c_salt, v print(f"\n Added salt concentration of {c_salt_calculated.to('mol/L')} given by {N_cation} cations and {N_anion} anions") return c_salt_calculated + def create_bond_in_espresso(self, bond_type, bond_parameters): + ''' + Creates either a harmonic or a FENE bond in ESPREesSo + + Args: + bond_type (`str`) : label to identify the potential to model the bond. + bond_parameters (`dict`) : parameters of the potential of the bond. + + Note: + Currently, only HARMONIC and FENE bonds are supported. + + For a HARMONIC bond the dictionary must contain: + + - k (`obj`) : Magnitude of the bond. It should have units of energy/length**2 + using the `pmb.units` UnitRegistry. + - r_0 (`obj`) : Equilibrium bond length. It should have units of length using + the `pmb.units` UnitRegistry. + + For a FENE bond the dictionay must additionally contain: + + - d_r_max (`obj`): Maximal stretching length for FENE. It should have + units of length using the `pmb.units` UnitRegistry. Default 'None'. + + Returns: + bond_object (`obj`): a harmonic or a FENE bond object in ESPREesSo + ''' + from espressomd import interactions + + valid_bond_types = ["harmonic", "FENE"] + + if 'k' in bond_parameters: + bond_magnitude = bond_parameters['k'].to('reduced_energy / reduced_length**2') + else: + raise ValueError("Magnitud of the potential (k) is missing") + + if bond_type == 'harmonic': + if 'r_0' in bond_parameters: + bond_length = bond_parameters['r_0'].to('reduced_length') + else: + raise ValueError("Equilibrium bond length (r_0) is missing") + bond_object = interactions.HarmonicBond(k = bond_magnitude.magnitude, + r_0 = bond_length.magnitude) + elif bond_type == 'FENE': + if 'r_0' in bond_parameters: + bond_length = bond_parameters['r_0'].to('reduced_length').magnitude + else: + print("WARNING: No value provided for r_0. Defaulting to r_0 = 0") + bond_length=0 + if 'd_r_max' in bond_parameters: + max_bond_stret = bond_parameters['d_r_max'].to('reduced_length') + else: + raise ValueError("Maximal stretching length (d_r_max) is missing") + bond_object = interactions.FeneBond(r_0 = bond_length, + k = bond_magnitude.magnitude, + d_r_max = max_bond_stret.magnitude) + else: + raise ValueError(f"Bond type {bond_type} currently not implemented in pyMBE, valid types are {valid_bond_types}") + + return bond_object + + def create_counterions(self, object_name, cation_name, anion_name, espresso_system,verbose=True): """ Creates particles of `cation_name` and `anion_name` in `espresso_system` to counter the net charge of `pmb_object`. @@ -1183,80 +1244,94 @@ def define_AA_residues(self, sequence, model): side_chains = side_chains) residue_list.append('AA-'+residue_name) return residue_list - - def define_bond(self, bond_object, bond_type, particle_name1, particle_name2): - """ - Defines a pmb object of type `bond` in `pymbe.df` + + def define_bond(self, bond_type, bond_parameters, particle_pairs): + ''' + Defines a pmb object of type `bond` in `pymbe.df`. + Args: - bond_object (`cls`): instance of a bond object from espressomd library - bond_type (`str`): label identifying the used bonded potential - particle_name1 (`str`): `name` of the first `particle` to be bonded. - particle_name2 (`str`): `name` of the second `particle` to be bonded. - + bond_type (`str`) : label to identify the potential to model the bond. + bond_parameters (`dict`) : parameters of the potential of the bond. + particle_pairs (`lst`) : list of the `names` of the `particles` to be bonded. + Note: - - Currently, only harmonic and FENE bonds are supported. - """ - - valid_bond_types=["harmonic", "FENE"] - if bond_type not in valid_bond_types: - raise ValueError(f"Bond type {bond_type} currently not implemented in pyMBE, valid types are {valid_bond_types}") + Currently, only HARMONIC and FENE bonds are supported. - lj_parameters=self.get_lj_parameters(particle_name1=particle_name1, - particle_name2=particle_name2, - combining_rule='Lorentz-Berthelot') + For a HARMONIC bond the dictionary must contain the following parameters: - cutoff = 2**(1.0/6.0) * lj_parameters["sigma"] - l0 = self.calculate_initial_bond_length(bond_object=bond_object, - bond_type=bond_type, - epsilon=lj_parameters["epsilon"], - sigma=lj_parameters["sigma"], - cutoff=cutoff) + - k (`obj`) : Magnitude of the bond. It should have units of energy/length**2 + using the `pmb.units` UnitRegistry. + - r_0 (`obj`) : Equilibrium bond length. It should have units of length using + the `pmb.units` UnitRegistry. + + For a FENE bond the dictionay must contain the same parameters as for a HARMONIC bond and: + + - d_r_max (`obj`): Maximal stretching length for FENE. It should have + units of length using the `pmb.units` UnitRegistry. Default 'None'. + ''' + + bond_object=self.create_bond_in_espresso(bond_type, bond_parameters) + + for particle_name1, particle_name2 in particle_pairs: + + lj_parameters=self.get_lj_parameters(particle_name1 = particle_name1, + particle_name2 = particle_name2, + combining_rule = 'Lorentz-Berthelot') + + cutoff = 2**(1.0/6.0) * lj_parameters["sigma"] + + l0 = self.calculate_initial_bond_length(bond_object = bond_object, + bond_type = bond_type, + epsilon = lj_parameters["epsilon"], + sigma = lj_parameters["sigma"], + cutoff = cutoff) + index = len(self.df) + for label in [particle_name1+'-'+particle_name2,particle_name2+'-'+particle_name1]: + self.check_if_name_is_defined_in_df(name=label, pmb_type_to_be_defined="bond") + self.df.at [index,'name']= particle_name1+'-'+particle_name2 + self.df.at [index,'bond_object'] = bond_object + self.df.at [index,'l0'] = l0 + self.add_value_to_df(index=index, + key=('pmb_type',''), + new_value='bond') + self.add_value_to_df(index=index, + key=('parameters_of_the_potential',''), + new_value=self.json.dumps(bond_object.get_params())) - index = len(self.df) - for label in [particle_name1+'-'+particle_name2,particle_name2+'-'+particle_name1]: - self.check_if_name_is_defined_in_df(name=label, pmb_type_to_be_defined="bond") - self.df.at [index,'name']= particle_name1+'-'+particle_name2 - self.df.at [index,'bond_object'] = bond_object - self.df.at [index,'l0'] = l0 - self.add_value_to_df(index=index, - key=('pmb_type',''), - new_value='bond') - self.add_value_to_df(index=index, - key=('parameters_of_the_potential',''), - new_value=self.json.dumps(bond_object.get_params())) self.df['parameters_of_the_potential'] = self.df['parameters_of_the_potential'].apply (lambda x: eval(str(x)) if self.pd.notnull(x) else x) return - def define_default_bond(self, bond_object, bond_type, epsilon=None, sigma=None, cutoff=None): + + def define_default_bond(self, bond_type, bond_parameters, epsilon=None, sigma=None, cutoff=None): """ Asigns `bond` in `pmb.df` as the default bond. The LJ parameters can be optionally provided to calculate the initial bond length Args: - bond (`obj`): instance of a bond object from the espressomd library. - bond_type (`str`): label identifying the used bonded potential - sigma (`float`, optional): LJ sigma of the interaction between the particles + bond_type (`str`) : label to identify the potential to model the bond. + bond_parameters (`dict`) : parameters of the potential of the bond. + sigma (`float`, optional) : LJ sigma of the interaction between the particles epsilon (`float`, optional): LJ epsilon for the interaction between the particles - cutoff (`float`, optional): cutoff-radius of the LJ interaction + cutoff (`float`, optional) : cutoff-radius of the LJ interaction Note: - Currently, only harmonic and FENE bonds are supported. """ - valid_bond_types=["harmonic", "FENE"] - if bond_type not in valid_bond_types: - raise ValueError(f"Bond type {bond_type} currently not implemented in pyMBE, valid types are {valid_bond_types}") + + bond_object=self.create_bond_in_espresso(bond_type, bond_parameters) + if epsilon is None: epsilon=1*self.units('reduced_energy') if sigma is None: sigma=1*self.units('reduced_length') if cutoff is None: cutoff=2**(1.0/6.0)*self.units('reduced_length') - l0 = self.calculate_initial_bond_length(bond_object=bond_object, - bond_type=bond_type, - epsilon=epsilon, - sigma=sigma, - cutoff=cutoff) + l0 = self.calculate_initial_bond_length(bond_object = bond_object, + bond_type = bond_type, + epsilon = epsilon, + sigma = sigma, + cutoff = cutoff) if self.check_if_name_is_defined_in_df(name='default',pmb_type_to_be_defined='bond'): return @@ -1264,16 +1339,15 @@ def define_default_bond(self, bond_object, bond_type, epsilon=None, sigma=None, index = max(self.df.index)+1 else: index = 0 - self.df.at [index,'name'] = 'default' + self.df.at [index,'name'] = 'default' self.df.at [index,'bond_object'] = bond_object - self.df.at [index,'l0'] = l0 - self.add_value_to_df(index=index, - key=('pmb_type',''), - new_value='bond') - self.add_value_to_df(index=index, - key=('parameters_of_the_potential',''), - new_value=self.json.dumps(bond_object.get_params())) - + self.df.at [index,'l0'] = l0 + self.add_value_to_df(index = index, + key = ('pmb_type',''), + new_value = 'bond') + self.add_value_to_df(index = index, + key = ('parameters_of_the_potential',''), + new_value = self.json.dumps(bond_object.get_params())) self.df['parameters_of_the_potential'] = self.df['parameters_of_the_potential'].apply (lambda x: eval(str(x)) if self.pd.notnull(x) else x) return @@ -1974,21 +2048,27 @@ def load_interaction_parameters(self, filename, verbose=False): sequence=param_dict.pop('sequence'), model=param_dict.pop('model')) elif object_type == 'bond': - name1 = param_dict.pop('name1') - name2 = param_dict.pop('name2') + particle_pairs = param_dict.pop('particle_pairs') + bond_parameters = param_dict.pop('bond_parameters') bond_type = param_dict.pop('bond_type') if bond_type == 'harmonic': - k = self.create_variable_with_units(variable=param_dict.pop('k')) - r_0 = self.create_variable_with_units(variable=param_dict.pop('r_0')) - bond = interactions.HarmonicBond(k=k.to('reduced_energy / reduced_length**2').magnitude, r_0=r_0.to('reduced_length').magnitude) + k = self.create_variable_with_units(variable=bond_parameters.pop('k')) + r_0 = self.create_variable_with_units(variable=bond_parameters.pop('r_0')) + bond = {'r_0' : r_0, + 'k' : k, + } + elif bond_type == 'FENE': - k = self.create_variable_with_units(variable=param_dict.pop('k')) - r_0 = self.create_variable_with_units(variable=param_dict.pop('r_0')) - d_r_max = self.create_variable_with_units(variable=param_dict.pop('d_r_max')) - bond = interactions.FeneBond(k=k.to('reduced_energy / reduced_length**2').magnitude, r_0=r_0.to('reduced_length').magnitude, d_r_max=d_r_max.to('reduced_length').magnitude) + k = self.create_variable_with_units(variable=bond_parameters.pop('k')) + r_0 = self.create_variable_with_units(variable=bond_parameters.pop('r_0')) + d_r_max = self.create_variable_with_units(variable=bond_parameters.pop('d_r_max')) + bond = {'r_0' : r_0, + 'k' : k, + 'd_r_max': d_r_max, + } else: raise ValueError("Current implementation of pyMBE only supports harmonic and FENE bonds") - self.define_bond(bond_object=bond, bond_type=bond_type, particle_name1=name1, particle_name2=name2) + self.define_bond(bond_type=bond_type, bond_parameters=bond, particle_pairs=particle_pairs) else: raise ValueError(object_type+' is not a known pmb object type') if verbose: diff --git a/samples/Beyer2024/weak_polyelectrolyte_dialysis.py b/samples/Beyer2024/weak_polyelectrolyte_dialysis.py index 3e6840c..8dfe5a3 100644 --- a/samples/Beyer2024/weak_polyelectrolyte_dialysis.py +++ b/samples/Beyer2024/weak_polyelectrolyte_dialysis.py @@ -93,12 +93,16 @@ pmb.define_residue(name='rA', central_bead="A", side_chains=[]) pmb.define_molecule(name=polyacid_name, residue_list=['rA']*Chain_length) + +bond_type = 'FENE' fene_spring_constant = 30 * pmb.units('reduced_energy / reduced_length**2') fene_r_max = 1.5 * pmb.units('reduced_length') -fene_bond = interactions.FeneBond(k=fene_spring_constant.to('reduced_energy / reduced_length**2').magnitude, - d_r_max=fene_r_max.to('reduced_length').magnitude) -pmb.define_bond(bond_object = fene_bond, particle_name1='A', particle_name2='A', bond_type="FENE") +fene_bond = {'k' : fene_spring_constant, + 'd_r_max': fene_r_max, + } + +pmb.define_bond(bond_type = bond_type, bond_parameters = fene_bond, particle_pairs = [['A','A']]) # Parameters of the small ions proton_name = 'Hplus' diff --git a/samples/branched_polyampholyte.py b/samples/branched_polyampholyte.py index e1c0b43..06735eb 100644 --- a/samples/branched_polyampholyte.py +++ b/samples/branched_polyampholyte.py @@ -97,15 +97,16 @@ residue_list = 5*["Res_1"] + 5*["Res_2"]) # Define bonds -generic_bond_length = 0.4 * pmb.units.nm +bond_type = 'harmonic' +generic_bond_lenght=0.4 * pmb.units.nm generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond = interactions.HarmonicBond( - k=generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude, - r_0=generic_bond_length.to('reduced_length').magnitude) +harmonic_bond = {'r_0' : generic_bond_lenght, + 'k' : generic_harmonic_constant, + } -pmb.define_default_bond(bond_object = generic_bond, - bond_type="harmonic") + +pmb.define_default_bond(bond_type = bond_type, bond_parameters = harmonic_bond) # Solution parameters cation_name = 'Na' diff --git a/samples/peptide.py b/samples/peptide.py index 4bda9d4..df703be 100644 --- a/samples/peptide.py +++ b/samples/peptide.py @@ -79,10 +79,14 @@ generic_bond_lenght=0.4 * pmb.units.nm generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond = interactions.HarmonicBond(k=generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude, - r_0=generic_bond_lenght.to('reduced_length').magnitude) -pmb.define_default_bond(bond_object = generic_bond, bond_type="harmonic") +HARMONIC_parameters = {'r_0' : generic_bond_lenght, + 'k' : generic_harmonic_constant, + } + +pmb.define_default_bond(bond_type = 'harmonic', + bond_parameters = HARMONIC_parameters) + # Defines the peptine in the pyMBE data frame peptide_name = 'generic_peptide' diff --git a/samples/peptide_mixture_grxmc_ideal.py b/samples/peptide_mixture_grxmc_ideal.py index dd1b6c4..76e6a25 100644 --- a/samples/peptide_mixture_grxmc_ideal.py +++ b/samples/peptide_mixture_grxmc_ideal.py @@ -91,12 +91,19 @@ warnings.simplefilter('error') pmb.load_pka_set(path_to_pka) + +# Defines the bonds + +bond_type = 'harmonic' generic_bond_lenght=0.4 * pmb.units.nm generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond = interactions.HarmonicBond(k=generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude, - r_0=generic_bond_lenght.to('reduced_length').magnitude) -pmb.define_default_bond(bond_object = generic_bond, bond_type="harmonic") +harmonic_bond = {'r_0' : generic_bond_lenght, + 'k' : generic_harmonic_constant, + } + + +pmb.define_default_bond(bond_type = bond_type, bond_parameters = harmonic_bond) # Defines the peptides in the pyMBE data frame peptide1 = 'generic_peptide1' diff --git a/testsuite/bond_tests.py b/testsuite/bond_tests.py new file mode 100644 index 0000000..95016c6 --- /dev/null +++ b/testsuite/bond_tests.py @@ -0,0 +1,353 @@ +# Import pyMBE and other libraries +import pyMBE +import numpy as np + +# Create an instance of pyMBE library +pmb = pyMBE.pymbe_library() + +def check_bond_setup(bond_object, input_parameters, bond_type): + """ + Checks that pyMBE sets up a bond object correctly. + + Args: + bond_object(`espressomd.interactions`): instance of a espresso bond object. + input_parameters(`dict`): dictionary with the parameters for the bond potential. + bond_type(`str`): label identifying the bond type + """ + # Check that pyMBE stores the correct type of bond + if bond_type.lower() not in str(bond_object).lower(): + raise Exception("*** Test failed: pyMBE does not store the correct type of bond ***") + # Check that pyMBE defines the bond with the correct parameters + bond_params = bond_object.get_params() + reduced_units = {'r_0' : 'reduced_length', + 'k' : 'reduced_energy / reduced_length**2', + 'd_r_max': 'reduced_length'} + for key in input_parameters.keys(): + np.testing.assert_equal(actual=bond_params[key], + desired=input_parameters[key].m_as(reduced_units[key]), + verbose=True) + return 0 + +#### DEFINE_BOND TESTS ### + +print(f"*** Unit test: check that define_bond sets up a harmonic bond correctly ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +pmb.define_bond(bond_type = bond_type, + bond_parameters = bond, + particle_pairs = [['A', 'A']]) + +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'].values[0], + input_parameters=bond, + bond_type=bond_type) + +print("*** Unit test passed ***") + +# Clean pmb.df +pmb.setup_df() + +print(f"*** Unit test: check that define_bond sets up a FENE bond correctly ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'FENE' + +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2'), + 'd_r_max': 0.8 * pmb.units.nm} + +pmb.define_bond(bond_type = bond_type, + bond_parameters = bond, + particle_pairs = [['A', 'A']]) + + +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'].values[0], + input_parameters=bond, + bond_type=bond_type) + +print("*** Unit test passed ***") + +# Clean pmb.df + +pmb.setup_df() + +print(f"*** Unit test: check that define_bond sets up a harmonic and a FENE bonds correctly in the same script ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) +pmb.define_particle(name='B', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type_1 = 'harmonic' +bond_1 = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +pmb.define_bond(bond_type = bond_type_1, + bond_parameters = bond_1, + particle_pairs = [['A', 'A']]) + +bond_type_2 = 'FENE' + +bond_2 = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2'), + 'd_r_max': 0.8 * pmb.units.nm} + +pmb.define_bond(bond_type = bond_type_2, + bond_parameters = bond_2, + particle_pairs = [['B', 'B']]) + +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'][2], + input_parameters=bond_1, + bond_type=bond_type_1) +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'][3], + input_parameters=bond_2, + bond_type=bond_type_2) + +print("*** Unit test passed ***") + +# Clean pmb.df + +pmb.setup_df() + +print(f"*** Unit test: check that define_bond raises a ValueError if the provided bond_type differs from the two currently implemented (harmonic and FENE) ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'Quartic' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + "particle_pairs" : [['A', 'A']] + } + +np.testing.assert_raises(ValueError, pmb.define_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_bond raises a ValueError if the provided dictionary does not have the Equilibrium length (r_0) for setting up a harmonic bond ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + "particle_pairs" : [['A', 'A']] + } + +np.testing.assert_raises(ValueError, pmb.define_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_bond raises a ValueError if the provided dictionary does not have the Magnitud of the potential (k) for setting up a harmonic bond ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'r_0' : 0.4*pmb.units.nm} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + "particle_pairs" : [['A', 'A']] + } + +np.testing.assert_raises(ValueError, pmb.define_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_bond raises a ValueError if the provided dictionary does not have the Magnitud of the potential (k) for setting up a FENE bond ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'FENE' +bond = {'r_0' : 0.4*pmb.units.nm, + 'd_r_max': 0.8 * pmb.units.nm} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + "particle_pairs" : [['A', 'A']] + } + +np.testing.assert_raises(ValueError, pmb.define_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_bond raises a ValueError if the provided dictionary does not have the Maximal stretching length (d_r_max) for setting up a FENE bond ***") + +# Define dummy particle + +pmb.define_particle(name='A', q=0, sigma=0.4*pmb.units.nm, epsilon=1*pmb.units('reduced_energy')) + +# Define dummy bond + +bond_type = 'FENE' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + "particle_pairs" : [['A', 'A']] + } + +np.testing.assert_raises(ValueError, pmb.define_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +#### DEFINE_DEFAULT_BOND TESTS ### + +print(f"*** Unit test: check that define_default_bond sets up a harmonic bond correctly ***") + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +pmb.define_default_bond(bond_type = bond_type, + bond_parameters = bond) + +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'].values[0], + input_parameters=bond, + bond_type=bond_type) + +print("*** Unit test passed ***") + +# Clean pmb.df + +pmb.setup_df() + +print(f"*** Unit test: check that define_dafult_bond sets up a FENE bond correctly ***") + +# Define dummy bond + +bond_type = 'FENE' + +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2'), + 'd_r_max': 0.8 * pmb.units.nm} + +pmb.define_default_bond(bond_type = bond_type, + bond_parameters = bond) + +check_bond_setup(bond_object=pmb.filter_df(pmb_type = 'bond')['bond_object'].values[0], + input_parameters=bond, + bond_type=bond_type) + +print("*** Unit test passed ***") + +# Clean pmb.df + +pmb.setup_df() + +print(f"*** Unit test: check that define_default_bond raises a ValueError if the provided bond_type differs from the two currently implemented (harmonic and FENE) ***") + +# Define dummy bond + +bond_type = 'Quartic' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + } + +np.testing.assert_raises(ValueError, pmb.define_default_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_default_bond raises a ValueError if the provided dictionary does not have the Equilibrium length (r_0) for setting up a harmonic bond ***") + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + } + +np.testing.assert_raises(ValueError, pmb.define_default_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_default_bond raises a ValueError if the provided dictionary does not have the Magnitud of the potential (k) for setting up a harmonic bond ***") + +# Define dummy bond + +bond_type = 'harmonic' +bond = {'r_0' : 0.4*pmb.units.nm} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + } + +np.testing.assert_raises(ValueError, pmb.define_default_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_default_bond raises a ValueError if the provided dictionary does not have the Magnitud of the potential (k) for setting up a FENE bond ***") + +# Define dummy bond + +bond_type = 'FENE' +bond = {'r_0' : 0.4*pmb.units.nm, + 'd_r_max': 0.8 * pmb.units.nm} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + } + +np.testing.assert_raises(ValueError, pmb.define_default_bond, **input_parameters) + +print(f"*** Unit test passed ***") + +print(f"*** Unit test: check that define_default_bond raises a ValueError if the provided dictionary does not have the Maximal stretching length (d_r_max) for setting up a FENE bond ***") + +# Define dummy bond + +bond_type = 'FENE' +bond = {'r_0' : 0.4*pmb.units.nm, + 'k' : 400 * pmb.units('reduced_energy / reduced_length**2')} + +input_parameters={"bond_type": bond_type, + "bond_parameters" : bond, + } + +np.testing.assert_raises(ValueError, pmb.define_default_bond, **input_parameters) + +print(f"*** Unit test passed ***") + diff --git a/testsuite/create_molecule_position_test.py b/testsuite/create_molecule_position_test.py index 7fffa5c..ff34bc7 100644 --- a/testsuite/create_molecule_position_test.py +++ b/testsuite/create_molecule_position_test.py @@ -31,14 +31,15 @@ side_chains = ['side_mon', 'side_mon'] ) - - +bond_type = 'harmonic' generic_bond_lenght=0.4 * pmb.units.nm generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond = interactions.HarmonicBond(k=generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude, - r_0=generic_bond_lenght.to('reduced_length').magnitude) -pmb.define_default_bond(bond_object = generic_bond, bond_type="harmonic") +harmonic_bond = {'r_0' : generic_bond_lenght, + 'k' : generic_harmonic_constant, + } + +pmb.define_default_bond(bond_type = bond_type, bond_parameters = harmonic_bond) # Defines the peptine in the pyMBE data frame molecule_name = 'generic_molecule' @@ -81,4 +82,4 @@ np.testing.assert_almost_equal(pos_list, central_bead_pos) -print(f"*** Unit test passed ***") \ No newline at end of file +print(f"*** Unit test passed ***") diff --git a/testsuite/read-write-df_test.py b/testsuite/read-write-df_test.py index 567437e..df22337 100644 --- a/testsuite/read-write-df_test.py +++ b/testsuite/read-write-df_test.py @@ -55,12 +55,16 @@ "Res_2"]) +bond_type = 'harmonic' generic_bond_lenght=0.4 * pmb.units.nm generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2') -generic_bond = interactions.HarmonicBond(k=generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude, - r_0=generic_bond_lenght.to('reduced_length').magnitude) -pmb.define_default_bond(bond_object = generic_bond, bond_type="harmonic") +harmonic_bond = {'r_0' : generic_bond_lenght, + 'k' : generic_harmonic_constant, + } + + +pmb.define_default_bond(bond_type = bond_type, bond_parameters = harmonic_bond) # Solution parameters cation_name = 'Na' diff --git a/tutorials/pyMBE_tutorial.ipynb b/tutorials/pyMBE_tutorial.ipynb index 816b748..ca0efd0 100644 --- a/tutorials/pyMBE_tutorial.ipynb +++ b/tutorials/pyMBE_tutorial.ipynb @@ -904,741 +904,87 @@ }, { "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [], - "source": [ - "generic_bond_lenght = 0.5*pmb.units.nm\n", - "generic_harmonic_constant = 400*pmb.units('reduced_energy / nm**2') \n", - "\n", - "generic_bond = interactions.HarmonicBond(k = generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude,\n", - " r_0 = generic_bond_lenght.to('reduced_length').magnitude)\n", - "\n", - "# backbone-backbone bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead,\n", - " particle_name2 = PDha_backbone_bead)\n", - "\n", - "# backbone-side_chain_1 (COOH) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead,\n", - " particle_name2 = PDha_carboxyl_bead)\n", - "\n", - "# backbone-side_chain_2 (NH3) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead,\n", - " particle_name2 = PDha_amine_bead)\n", - "\n", - "pmb.add_bonds_to_espresso(espresso_system = espresso_system)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "NOTE: Currently, only harmonic and FENE bonds are supported." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, one can use the residues to define the polymer sequence given by the argument `residue_list`. One needs to add one residue in `residue_list` per each residue in the polymer chain. For instance a decamer should be created as follows:" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "PDha_polymer = 'PDha'\n", - "N_monomers = 10\n", - "\n", - "pmb.define_molecule(name = PDha_polymer,\n", - " residue_list = [PDha_residue]*N_monomers)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "After defining the polymer, we are ready to create one PdHa polymer in the center of the simulation box." - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "N_polymers = 1\n", - "\n", - "pmb.create_pmb_object(name = PDha_polymer, \n", - " number_of_objects = N_polymers,\n", - " espresso_system = espresso_system, \n", - " position = [[Box_L.to('reduced_length').magnitude/2]*3]) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can always track our molecules..." - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
namepmb_typeparticle_idresidue_idmolecule_idaciditysigmacutoffoffsetepsilonstate_one
labeles_typecharge
0BB-PDhaparticle000inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
1COOH-PDhaparticle100inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
2NH3-PDhaparticle200inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
9BB-PDhaparticle310inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
10COOH-PDhaparticle410inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
12NH3-PDhaparticle510inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
15BB-PDhaparticle620inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
16COOH-PDhaparticle720inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
18NH3-PDhaparticle820inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
21BB-PDhaparticle930inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
22COOH-PDhaparticle1030inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
24NH3-PDhaparticle1130inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
28BB-PDhaparticle1240inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
29COOH-PDhaparticle1340inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
31NH3-PDhaparticle1440inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
35BB-PDhaparticle1550inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
36COOH-PDhaparticle1650inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
38NH3-PDhaparticle1750inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
42BB-PDhaparticle1860inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
43COOH-PDhaparticle1960inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
45NH3-PDhaparticle2060inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
49BB-PDhaparticle2170inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
50COOH-PDhaparticle2270inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
52NH3-PDhaparticle2370inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
56BB-PDhaparticle2480inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
57COOH-PDhaparticle2580inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
59NH3-PDhaparticle2680inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
63BB-PDhaparticle2790inert0.4 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyBB-PDha00
64COOH-PDhaparticle2890inert0.5 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyCOOH-PDha10
66NH3-PDhaparticle2990inert0.3 nanometer1.122462048309373 reduced_length0 reduced_length1 reduced_energyNH3-PDha20
\n", - "
" - ], - "text/plain": [ - " name pmb_type particle_id residue_id molecule_id acidity \\\n", - " \n", - "0 BB-PDha particle 0 0 0 inert \n", - "1 COOH-PDha particle 1 0 0 inert \n", - "2 NH3-PDha particle 2 0 0 inert \n", - "9 BB-PDha particle 3 1 0 inert \n", - "10 COOH-PDha particle 4 1 0 inert \n", - "12 NH3-PDha particle 5 1 0 inert \n", - "15 BB-PDha particle 6 2 0 inert \n", - "16 COOH-PDha particle 7 2 0 inert \n", - "18 NH3-PDha particle 8 2 0 inert \n", - "21 BB-PDha particle 9 3 0 inert \n", - "22 COOH-PDha particle 10 3 0 inert \n", - "24 NH3-PDha particle 11 3 0 inert \n", - "28 BB-PDha particle 12 4 0 inert \n", - "29 COOH-PDha particle 13 4 0 inert \n", - "31 NH3-PDha particle 14 4 0 inert \n", - "35 BB-PDha particle 15 5 0 inert \n", - "36 COOH-PDha particle 16 5 0 inert \n", - "38 NH3-PDha particle 17 5 0 inert \n", - "42 BB-PDha particle 18 6 0 inert \n", - "43 COOH-PDha particle 19 6 0 inert \n", - "45 NH3-PDha particle 20 6 0 inert \n", - "49 BB-PDha particle 21 7 0 inert \n", - "50 COOH-PDha particle 22 7 0 inert \n", - "52 NH3-PDha particle 23 7 0 inert \n", - "56 BB-PDha particle 24 8 0 inert \n", - "57 COOH-PDha particle 25 8 0 inert \n", - "59 NH3-PDha particle 26 8 0 inert \n", - "63 BB-PDha particle 27 9 0 inert \n", - "64 COOH-PDha particle 28 9 0 inert \n", - "66 NH3-PDha particle 29 9 0 inert \n", - "\n", - " sigma cutoff offset \\\n", - " \n", - "0 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "1 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "2 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "9 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "10 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "12 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "15 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "16 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "18 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "21 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "22 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "24 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "28 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "29 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "31 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "35 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "36 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "38 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "42 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "43 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "45 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "49 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "50 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "52 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "56 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "57 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "59 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "63 0.4 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "64 0.5 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "66 0.3 nanometer 1.122462048309373 reduced_length 0 reduced_length \n", - "\n", - " epsilon state_one \n", - " label es_type charge \n", - "0 1 reduced_energy BB-PDha 0 0 \n", - "1 1 reduced_energy COOH-PDha 1 0 \n", - "2 1 reduced_energy NH3-PDha 2 0 \n", - "9 1 reduced_energy BB-PDha 0 0 \n", - "10 1 reduced_energy COOH-PDha 1 0 \n", - "12 1 reduced_energy NH3-PDha 2 0 \n", - "15 1 reduced_energy BB-PDha 0 0 \n", - "16 1 reduced_energy COOH-PDha 1 0 \n", - "18 1 reduced_energy NH3-PDha 2 0 \n", - "21 1 reduced_energy BB-PDha 0 0 \n", - "22 1 reduced_energy COOH-PDha 1 0 \n", - "24 1 reduced_energy NH3-PDha 2 0 \n", - "28 1 reduced_energy BB-PDha 0 0 \n", - "29 1 reduced_energy COOH-PDha 1 0 \n", - "31 1 reduced_energy NH3-PDha 2 0 \n", - "35 1 reduced_energy BB-PDha 0 0 \n", - "36 1 reduced_energy COOH-PDha 1 0 \n", - "38 1 reduced_energy NH3-PDha 2 0 \n", - "42 1 reduced_energy BB-PDha 0 0 \n", - "43 1 reduced_energy COOH-PDha 1 0 \n", - "45 1 reduced_energy NH3-PDha 2 0 \n", - "49 1 reduced_energy BB-PDha 0 0 \n", - "50 1 reduced_energy COOH-PDha 1 0 \n", - "52 1 reduced_energy NH3-PDha 2 0 \n", - "56 1 reduced_energy BB-PDha 0 0 \n", - "57 1 reduced_energy COOH-PDha 1 0 \n", - "59 1 reduced_energy NH3-PDha 2 0 \n", - "63 1 reduced_energy BB-PDha 0 0 \n", - "64 1 reduced_energy COOH-PDha 1 0 \n", - "66 1 reduced_energy NH3-PDha 2 0 " - ] - }, - "execution_count": 18, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "bond_type = 'harmonic'\n", + "generic_bond_lenght=0.4 * pmb.units.nm\n", + "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n", + "\n", + "harmonic_bond = {'r_0' : generic_bond_lenght,\n", + " 'k' : generic_harmonic_constant\n", + " }\n", + "\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PDha_backbone_bead, PDha_backbone_bead],\n", + " [PDha_backbone_bead, PDha_carboxyl_bead],\n", + " [PDha_backbone_bead, PDha_amine_bead]])\n", + "\n", + "pmb.add_bonds_to_espresso(espresso_system = espresso_system)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "NOTE: Currently, only harmonic and FENE bonds are supported." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, one can use the residues to define the polymer sequence given by the argument `residue_list`. One needs to add one residue in `residue_list` per each residue in the polymer chain. For instance a decamer should be created as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "PDha_polymer = 'PDha'\n", + "N_monomers = 10\n", + "\n", + "pmb.define_molecule(name = PDha_polymer,\n", + " residue_list = [PDha_residue]*N_monomers)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After defining the polymer, we are ready to create one PdHa polymer in the center of the simulation box." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N_polymers = 1\n", + "\n", + "pmb.create_pmb_object(name = PDha_polymer, \n", + " number_of_objects = N_polymers,\n", + " espresso_system = espresso_system, \n", + " position = [[Box_L.to('reduced_length').magnitude/2]*3]) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can always track our particles..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "pmb.filter_df(pmb_type = 'particle')" ] @@ -1652,7 +998,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1672,45 +1018,9 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - "
\n", - "
" - ], - "text/plain": [ - "Empty DataFrame\n", - "Columns: []\n", - "Index: []" - ] - }, - "execution_count": 20, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "pmb.destroy_pmb_object_in_system(name = PDha_polymer, \n", " espresso_system = espresso_system)\n", @@ -1742,7 +1052,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1782,7 +1092,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1802,7 +1112,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1821,43 +1131,24 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "generic_bond_lenght = 0.5*pmb.units.nm\n", - "generic_harmonic_constant = 400*pmb.units('reduced_energy / nm**2') \n", - "\n", - "generic_bond = interactions.HarmonicBond(k = generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude,\n", - " r_0 = generic_bond_lenght.to('reduced_length').magnitude)\n", - "\n", - "# backbone-backbone bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_backbone_bead, \n", - " particle_name2 = PDAGA_backbone_bead)\n", - "\n", - "# backbone-cyclic amine bond (side -chain)\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_backbone_bead, \n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", - "\n", - "# cyclic amine-alpha carboxyl bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_alpha_carboxyl_bead,\n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", + "bond_type = 'harmonic'\n", + "generic_bond_lenght=0.4 * pmb.units.nm\n", + "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n", "\n", - "# cyclic amine-beta carboxyl bond\n", + "harmonic_bond = {'r_0' : generic_bond_lenght,\n", + " 'k' : generic_harmonic_constant,\n", + " }\n", "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_beta_carboxyl_bead,\n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PDAGA_backbone_bead, PDAGA_backbone_bead],\n", + " [PDAGA_backbone_bead, PDAGA_cyclic_amine_bead],\n", + " [PDAGA_alpha_carboxyl_bead, PDAGA_cyclic_amine_bead],\n", + " [PDAGA_beta_carboxyl_bead, PDAGA_cyclic_amine_bead]])\n", "\n", "pmb.add_bonds_to_espresso(espresso_system = espresso_system)" ] @@ -1871,7 +1162,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1891,7 +1182,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1912,7 +1203,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -1932,45 +1223,9 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - "
\n", - "
" - ], - "text/plain": [ - "Empty DataFrame\n", - "Columns: []\n", - "Index: []" - ] - }, - "execution_count": 28, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "pmb.destroy_pmb_object_in_system(name = PDAGA_polymer, \n", " espresso_system = espresso_system)\n", @@ -2002,7 +1257,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -2061,15 +1316,17 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "generic_bond_lenght = 0.5*pmb.units.nm\n", - "generic_harmonic_constant = 400*pmb.units('reduced_energy / nm**2') \n", + "bond_type = 'harmonic'\n", + "generic_bond_lenght=0.4 * pmb.units.nm\n", + "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n", "\n", - "generic_bond = interactions.HarmonicBond(k = generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude,\n", - " r_0 = generic_bond_lenght.to('reduced_length').magnitude)\n", + "harmonic_bond = {'r_0' : generic_bond_lenght,\n", + " 'k' : generic_harmonic_constant,\n", + " }\n", "\n", "######################################################################\n", "################################# PDha ###############################\n", @@ -2080,26 +1337,12 @@ " central_bead = PDha_backbone_bead ,\n", " side_chains = [PDha_carboxyl_bead,PDha_amine_bead])\n", "\n", - "# backbone-backbone bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead, \n", - " particle_name2 = PDha_backbone_bead)\n", "\n", - "# backbone-side_chain_1 (COOH) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead, \n", - " particle_name2 = PDha_carboxyl_bead)\n", - "\n", - "# backbone-side_chain_2 (NH3) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead, \n", - " particle_name2 = PDha_amine_bead)\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PDha_backbone_bead, PDha_backbone_bead],\n", + " [PDha_backbone_bead, PDha_carboxyl_bead],\n", + " [PDha_backbone_bead, PDha_amine_bead]])\n", "\n", "######################################################################\n", "################################ PDAGA ###############################\n", @@ -2115,44 +1358,20 @@ " central_bead = PDAGA_cyclic_amine_bead,\n", " side_chains = [PDAGA_alpha_carboxyl_bead,PDAGA_beta_carboxyl_bead])\n", "\n", - "# backbone-backbone bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_backbone_bead, \n", - " particle_name2 = PDAGA_backbone_bead)\n", - "\n", - "# backbone-cyclic amine bond (side -chain)\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_backbone_bead, \n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", - "\n", - "# cyclic amine-alpha carboxyl bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_alpha_carboxyl_bead, \n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", - "\n", - "# cyclic amine-beta carboxyl bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDAGA_beta_carboxyl_bead, \n", - " particle_name2 = PDAGA_cyclic_amine_bead)\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PDAGA_backbone_bead, PDAGA_backbone_bead],\n", + " [PDAGA_backbone_bead, PDAGA_cyclic_amine_bead],\n", + " [PDAGA_alpha_carboxyl_bead, PDAGA_cyclic_amine_bead],\n", + " [PDAGA_beta_carboxyl_bead, PDAGA_cyclic_amine_bead]])\n", "\n", "######################################################################\n", "############################# PDha - PDAGA ###########################\n", "######################################################################\n", "\n", - "# backbone PDha - backbone PDAGA bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PDha_backbone_bead, \n", - " particle_name2 = PDAGA_backbone_bead)\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PDha_backbone_bead, PDAGA_backbone_bead]])\n", "\n", "pmb.add_bonds_to_espresso(espresso_system = espresso_system)" ] @@ -2166,7 +1385,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -2187,7 +1406,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -2208,7 +1427,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -2228,45 +1447,9 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - "
\n", - "
" - ], - "text/plain": [ - "Empty DataFrame\n", - "Columns: []\n", - "Index: []" - ] - }, - "execution_count": 34, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "pmb.destroy_pmb_object_in_system(name = diblock_polymer, \n", " espresso_system = espresso_system)\n", @@ -2424,7 +1607,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -2442,22 +1625,9 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "FileNotFoundError", - "evalue": "[Errno 2] No such file or directory: '/home/pablo/Escritorio/pyMBE/reference_parameters/interaction_parameters/Lunkad2021.txt'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[36], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mpmb\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mload_interaction_parameters\u001b[49m\u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpmb\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_resource\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mreference_parameters/interaction_parameters/Lunkad2021.txt\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 2\u001b[0m pmb\u001b[38;5;241m.\u001b[39madd_bonds_to_espresso (espresso_system \u001b[38;5;241m=\u001b[39m espresso_system)\n", - "File \u001b[0;32m~/Escritorio/pyMBE/pyMBE.py:1836\u001b[0m, in \u001b[0;36mpymbe_library.load_interaction_parameters\u001b[0;34m(self, filename, verbose)\u001b[0m\n\u001b[1;32m 1834\u001b[0m without_units \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mq\u001b[39m\u001b[38;5;124m'\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mes_type\u001b[39m\u001b[38;5;124m'\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124macidity\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m 1835\u001b[0m with_units \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msigma\u001b[39m\u001b[38;5;124m'\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mepsilon\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[0;32m-> 1836\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m f:\n\u001b[1;32m 1837\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m line \u001b[38;5;129;01min\u001b[39;00m f:\n\u001b[1;32m 1838\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m line[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m#\u001b[39m\u001b[38;5;124m'\u001b[39m:\n", - "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/home/pablo/Escritorio/pyMBE/reference_parameters/interaction_parameters/Lunkad2021.txt'" - ] - } - ], + "outputs": [], "source": [ "pmb.load_interaction_parameters (filename = pmb.get_resource('parameters/peptides/Lunkad2021.txt'))\n", "pmb.add_bonds_to_espresso (espresso_system = espresso_system)" @@ -2578,7 +1748,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.8.10" }, "vscode": { "interpreter": { diff --git a/tutorials/solution_tutorial.ipynb b/tutorials/solution_tutorial.ipynb index ee846aa..06e1664 100644 --- a/tutorials/solution_tutorial.ipynb +++ b/tutorials/solution_tutorial.ipynb @@ -28,9 +28,24 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Current set of reduced units:\n", + "1.5 nanometer = 1 reduced_length\n", + "4.11640356238e-21 joule = 1 reduced_energy\n", + "Temperature: 298.15 kelvin\n", + "8.01088317e-19 coulomb = 1 reduced_charge\n", + "\n", + "The side of the simulation box is 7.5 nanometer = 4.999999999999999 reduced_length\n" + ] + } + ], "source": [ "# Import pyMBE and ESPResSo\n", "import pyMBE\n", @@ -77,7 +92,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -112,7 +127,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 2. Define the different residues for the PAA and PVAm." + "#### 2. Define the different residues and the bonds for the PAA and PVAm." ] }, { @@ -121,56 +136,28 @@ "metadata": {}, "outputs": [], "source": [ - "generic_bond_lenght = 0.5*pmb.units.nm\n", - "generic_harmonic_constant = 400*pmb.units('reduced_energy / nm**2') \n", - "\n", - "generic_bond = interactions.HarmonicBond(k = generic_harmonic_constant.to('reduced_energy / reduced_length**2').magnitude,\n", - " r_0 = generic_bond_lenght.to('reduced_length').magnitude)\n", - "\n", - "######################################################################\n", - "################################## PAA ###############################\n", - "######################################################################\n", - "\n", "PAA_residue = 'PAA_mon'\n", "pmb.define_residue(name = PAA_residue, \n", " central_bead = PAA_backbone_bead ,\n", " side_chains = [PAA_carboxyl_bead])\n", "\n", - "# backbone-side_chain (COOH) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PAA_backbone_bead, \n", - " particle_name2 = PAA_carboxyl_bead)\n", - "\n", - "\n", - "######################################################################\n", - "################################# PVam ###############################\n", - "######################################################################\n", - "\n", "PVam_residue = 'PVam_mon'\n", "pmb.define_residue( name = PVam_residue,\n", " central_bead = PVam_backbone_bead,\n", " side_chains = [PVam_amine_bead])\n", "\n", - "# backbone-side_chain (amine) bond\n", - "\n", - "pmb.define_bond(bond_object = generic_bond,\n", - " bond_type = 'harmonic',\n", - " particle_name1 = PVam_backbone_bead, \n", - " particle_name2 = PVam_amine_bead)\n", - "\n", - "\n", - "######################################################################\n", - "############################### PAA - PVam ###########################\n", - "######################################################################\n", + "bond_type = 'harmonic'\n", + "generic_bond_lenght=0.4 * pmb.units.nm\n", + "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n", "\n", - "# backbone PAA - backbone PVam bond\n", + "harmonic_bond = {'r_0' : generic_bond_lenght,\n", + " 'k' : generic_harmonic_constant}\n", "\n", - "pmb.define_bond(bond_object = generic_bond, \n", - " bond_type = 'harmonic',\n", - " particle_name1 = PAA_backbone_bead, \n", - " particle_name2 = PVam_backbone_bead)\n", + "pmb.define_bond(bond_type = bond_type,\n", + " bond_parameters = harmonic_bond,\n", + " particle_pairs = [[PAA_backbone_bead, PAA_carboxyl_bead],\n", + " [PVam_backbone_bead, PVam_amine_bead],\n", + " [PAA_backbone_bead, PVam_backbone_bead]])\n", "\n", "pmb.add_bonds_to_espresso(espresso_system = espresso_system)" ] @@ -272,7 +259,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.8.10" }, "vscode": { "interpreter": {