Skip to content

Commit

Permalink
Remove redundant sample script for G-RxMC
Browse files Browse the repository at this point in the history
  • Loading branch information
David Beyer committed Mar 12, 2024
1 parent cfc25cc commit 776b0ae
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 276 deletions.
117 changes: 82 additions & 35 deletions samples/peptide_mixture_grxmc_ideal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@
import espressomd
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import argparse
from tqdm import tqdm
from espressomd.io.writer import vtf
from espressomd import interactions
Expand All @@ -15,13 +17,25 @@
import pyMBE
pmb = pyMBE.pymbe_library()

# Command line arguments

valid_modes=["standard", "unified"]
parser = argparse.ArgumentParser(description='Script that runs a simulation of an ideal peptide mixture in the grand-reaction ensemble using pyMBE and ESPResSo.')
parser.add_argument('--mode',
type=str,
default= "standard",
help='set if the grand-reaction method is used with unified ions or not, valid modes are {valid_modes}')
args = parser.parse_args()



# The trajectories of the simulations will be stored using espresso built-up functions in separed files in the folder 'frames'
if not os.path.exists('./frames'):
os.makedirs('./frames')

#Import functions from handy_functions script
from handy_scripts.handy_functions import minimize_espresso_system_energy
from handy_scripts.handy_functions import block_analyze
from lib.handy_functions import minimize_espresso_system_energy
from lib.analysis import block_analyze

# Simulation parameters
pmb.set_reduced_units(unit_length=0.4*pmb.units.nm)
Expand All @@ -46,10 +60,10 @@
N_peptide2_chains = 10

# Load peptide parametrization from Lunkad, R. et al. Molecular Systems Design & Engineering (2021), 6(2), 122-131.
path_to_interactions=pmb.get_resource("reference_parameters/interaction_parameters/Lunkad2021.txt")
path_to_pka=pmb.get_resource("reference_parameters/pka_sets/Hass2015.txt")
pmb.load_interaction_parameters (filename=path_to_interactions)
pmb.load_pka_set (path_to_pka)
path_to_interactions=pmb.get_resource("parameters/peptides/Lunkad2021.txt")
path_to_pka=pmb.get_resource("parameters/pka_sets/Hass2015.txt")
pmb.load_interaction_parameters(filename=path_to_interactions)
pmb.load_pka_set(path_to_pka)

# Use a generic parametrization for the aminoacids not parametrized
not_parametrized_neutral_aminoacids = ['A','N','Q','G','I','L','M','F','P','O','S','U','T','W','V','J']
Expand Down Expand Up @@ -90,16 +104,25 @@
pmb.define_peptide (name=peptide2, sequence=sequence2, model=model)

# Solution parameters
proton_name = 'Hplus'
hydroxide_name = 'OHminus'
sodium_name = 'Na'
chloride_name = 'Cl'
c_salt=5e-3 * pmb.units.mol/ pmb.units.L

pmb.define_particle(name=proton_name, q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=hydroxide_name, q=-1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=sodium_name, q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=chloride_name, q=-1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
if args.mode == 'standard':
proton_name = 'Hplus'
hydroxide_name = 'OHminus'
sodium_name = 'Na'
chloride_name = 'Cl'

pmb.define_particle(name=proton_name, q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=hydroxide_name, q=-1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=sodium_name, q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=chloride_name, q=-1, diameter=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, q=1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))
pmb.define_particle(name=anion_name, q=-1, diameter=0.35*pmb.units.nm, epsilon=1*pmb.units('reduced_energy'))


# System parameters
Expand All @@ -116,10 +139,18 @@
# 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_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)
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

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

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:
vtf.writevsf(espresso_system, coordinates)
Expand All @@ -133,15 +164,17 @@

print("The box length of your system is", L.to('reduced_length'), L.to('nm'))

RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=2, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
if args.mode == 'standard':
RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=2, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
elif args.mode == 'unified':
RE, 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, SEED=SEED)
print('The acid-base reaction has been sucessfully setup for ', sucessful_reactions_labels)

# Setup espresso to track the ionization of the acid/basic groups in peptide
type_map =pmb.get_type_map()
types = list (type_map.values())
espresso_system.setup_type_map(type_list = types)


# Setup the non-interacting type for speeding up the sampling of the reactions
non_interacting_type = max(type_map.values())+1
RE.set_non_interacting_type (type=non_interacting_type)
Expand All @@ -164,23 +197,31 @@

N_frame=0
Z_pH=[] # List of the average global charge at each pH
err_Z_pH=[] # List of the error of the global charge at each pH
xi_plus=[] # List of the average partition coefficient of positive ions
err_xi_plus=[] # List of the error of the partition coefficient of positive ions

particle_id_list = pmb.df.loc[~pmb.df['molecule_id'].isna()].particle_id.dropna().to_list()

#Save the pyMBE dataframe in a CSV file
pmb.df.to_csv('df.csv')

# Main loop for performing simulations at different pH-values
labels_obs=["time","charge","num_plus"]

for index in tqdm(range(len(pH_range))):

pH_value=pH_range[index]
# Sample list inicialization
Z_sim=[]
num_plus=[]

RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions (pH_res=pH_value, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
time_series={}

for label in labels_obs:
time_series[label]=[]

if args.mode == 'standard':
RE, sucessful_reactions_labels, ionic_strength_res = pmb.setup_grxmc_reactions(pH_res=pH_value, c_salt_res=c_salt, proton_name=proton_name, hydroxide_name=hydroxide_name, salt_cation_name=sodium_name, salt_anion_name=chloride_name, SEED=SEED)
elif args.mode == 'unified':
RE, 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, SEED=SEED)

# Inner loop for sampling each pH value

Expand All @@ -197,26 +238,32 @@
for pid in particle_id_list:
z_one_object +=espresso_system.part.by_id(pid).q

Z_sim.append(np.mean((z_one_object)))
num_plus.append(espresso_system.number_of_particles(type=type_map["Na"])+espresso_system.number_of_particles(type=type_map["Hplus"]))
time_series["time"].append(espresso_system.time)
time_series["charge"].append(np.mean((z_one_object)))

if args.mode == 'standard':
time_series["num_plus"].append(espresso_system.number_of_particles(type=type_map["Na"])+espresso_system.number_of_particles(type=type_map["Hplus"]))
elif args.mode == 'unified':
time_series["num_plus"].append(espresso_system.number_of_particles(type=type_map["Na"]))

if (step % N_samples_print == 0) :
N_frame+=1
with open('frames/trajectory'+str(N_frame)+'.vtf', mode='w+t') as coordinates:
vtf.writevsf(espresso_system, coordinates)
vtf.writevcf(espresso_system, coordinates)

Z_pH.append(Z_sim)
concentration_plus = (num_plus/(pmb.N_A * L**3)).to('mol/L')
# Estimate the statistical error and the autocorrelation time of the data
print("Net charge analysis")
processed_data = block_analyze(full_data=pd.DataFrame(time_series, columns=labels_obs))

Z_pH.append(processed_data["mean", "charge"])
err_Z_pH.append(processed_data["err_mean", "charge"])
concentration_plus = (processed_data["mean", "num_plus"]/(pmb.N_A * L**3)).to('mol/L')
err_concentration_plus = (processed_data["err_mean", "num_plus"]/(pmb.N_A * L**3)).to('mol/L')
xi_plus.append(concentration_plus/ionic_strength_res)
err_xi_plus.append(err_concentration_plus/ionic_strength_res)
print("pH = {:6.4g} done".format(pH_value))

# Estimate the statistical error and the autocorrelation time of the data

print("Net charge analysis")
av_net_charge, err_net_charge, tau_net_charge, block_size_net_charge = block_analyze(input_data=Z_pH)
av_xi_plus, err_xi_plus, tau_xi_plus, block_size_xi_plus = block_analyze(input_data=xi_plus)

# 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(espresso_system=espresso_system, object_names=[peptide1, peptide2], c_salt=c_salt, pH_list=pH_range_HH)
Expand All @@ -225,7 +272,7 @@
xi = HH_charge_dict["partition_coefficients"]

fig, ax = plt.subplots(figsize=(10, 7))
plt.errorbar(pH_range, np.asarray(av_net_charge)/N_peptide1_chains, yerr=err_net_charge/N_peptide1_chains, fmt = 'o', capsize=3, label='Simulation')
plt.errorbar(pH_range, np.asarray(Z_pH)/N_peptide1_chains, yerr=np.asarray(err_Z_pH)/N_peptide1_chains, fmt = 'o', capsize=3, label='Simulation')
ax.plot(pH_range_HH, np.asarray(Z_HH_Donnan[peptide1])+np.asarray(Z_HH_Donnan[peptide2]), "--r", label='HH+Donnan')
plt.legend()
plt.xlabel('pH')
Expand All @@ -234,7 +281,7 @@
plt.close()

fig, ax = plt.subplots(figsize=(10, 7))
plt.errorbar(pH_range, av_xi_plus, yerr=err_xi_plus, fmt = 'o', capsize=3, label='Simulation')
plt.errorbar(pH_range, xi_plus, yerr=err_xi_plus, fmt = 'o', capsize=3, label='Simulation')
ax.plot(pH_range_HH, np.asarray(xi), "-k", label='HH+Donnan')
plt.legend()
plt.xlabel('pH')
Expand Down
Loading

0 comments on commit 776b0ae

Please sign in to comment.