Fix histatin-5 case, rename parameters folders

parameters/interaction_parameters/Blanco2020.txt → parameters/peptides/Blanco2020.txt

@@ -1,5 +1,6 @@
 # Parameters from Blanco et al. Soft Matter, 17(3), 655-669, 2021.
 {"object_type":"particle", "name": "D", "acidity": "acidic", "diameter": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}}
+{"object_type":"particle", "name": "E", "acidity": "acidic", "diameter": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}}
 {"object_type":"particle", "name": "n", "acidity": "basic", "diameter": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}}
 {"object_type":"particle", "name": "S", "q":0, "acidity": "inert", "diameter": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}}
 {"object_type":"particle", "name": "H", "acidity": "basic", "diameter": {"value":0.4, "units":"nm"}, "epsilon":{"value":1, "units":"reduced_energy"}}
@@ -15,7 +16,10 @@
 {"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"}}

samples/Beyer2024/create_paper_data.py

@@ -50,15 +50,15 @@
     elif fig_label == "6b":
         sequence="E"*5+"H"*5
     elif fig_label == "6c":
-        sequence="nDSHAKRHHGYKRKFHHSHRGYc"
+        sequence="nDSHAKRHHGYKRKFHEKHHSHRGYc"
     else:
         raise RuntimeError()
     pH_range = np.linspace(2, 12, num=21)
 
     for pH in pH_range:
         run_command=f"python3 {script_path} --sequence {sequence} --pH {pH} --mode {mode}"
         print(run_command)
-        #os.system(run_command)
+        os.system(run_command)
 
 # Read all files in the subdir
 data_files=[]
@@ -124,6 +124,10 @@
         pka_path=pmb.get_resource("parameters/pka_sets/CRC1991.txt")
     elif fig_label == "6c":
         pka_path=pmb.get_resource("parameters/pka_sets/Nozaki1967.txt")
+        # FIXME: this is only necessary due to an undesired feature in calculate_HH
+        # that forces to have all particles defined in pyMBE
+        par_path=pmb.get_resource("parameters/peptides/Blanco2020.txt")
+        pmb.load_interaction_parameters(par_path)
     else:
         raise RuntimeError()
     pmb.load_pka_set (filename=pka_path)
@@ -171,6 +175,7 @@
     # Plot data produced by pyMBE
 
     data=data.astype({'pH': 'float'}).sort_values(by="pH")
+    data=data[data.sequence == f"{sequence}"]
     plt.errorbar(data["pH"], 
                    data["mean","charge"], 
                    yerr=data["err_mean","charge"], 

samples/Beyer2024/peptide.py

@@ -52,30 +52,29 @@
 
 # Simulation parameters
 if mode == "short-run":
-    Nsamples = 100
-    MD_steps_per_sample = 100
+    Nsamples = 1000
+    MD_steps_per_sample = 1000
 
 if mode == "long-run":
     Nsamples = 5000
-    MD_steps_per_sample = 1000
+    MD_steps_per_sample = 5000
 
 SEED = 100
 dt = 0.01
 solvent_permitivity = 78.3
 pep_concentration = 5.56e-4 *pmb.units.mol/pmb.units.L
 
 # Sanity check
-
 Lunkad_test_sequences=["E"*5+"H"*5,"K"*5+"D"*5]
-Blanco_test_sequence=["nDSHAKRHHGYKRKFHHSHRGYc"]
+Blanco_test_sequence=["nDSHAKRHHGYKRKFHEKHHSHRGYc"]
 
 valid_sequences=Lunkad_test_sequences+Blanco_test_sequence
 
 if sequence not in valid_sequences:
     raise ValueError(f"ERROR: the only valid peptide sequence for this test script are {valid_sequences}")
 
 if sequence in Lunkad_test_sequences:
-    pmb.load_interaction_parameters (filename='parameters/interaction_parameters/Lunkad2021.txt') 
+    pmb.load_interaction_parameters (filename='parameters/peptides/Lunkad2021.txt') 
     pmb.load_pka_set (filename='parameters/pka_sets/CRC1991.txt')
     model = '2beadAA'  # Model with 2 beads per each aminoacid
     N_peptide_chains = 4
@@ -85,7 +84,7 @@
     chain_length=len(sequence)*2
 
 elif sequence in Blanco_test_sequence:
-    pmb.load_interaction_parameters (pmb.get_resource(path='parameters/interaction_parameters/Blanco2020.txt')) 
+    pmb.load_interaction_parameters (pmb.get_resource(path='parameters/peptides/Blanco2020.txt')) 
     pmb.load_pka_set (pmb.get_resource(path='parameters/pka_sets/Nozaki1967.txt'))
     model = '1beadAA'
     N_peptide_chains = 1

testsuite/data/Blanco2020a.csv

@@ -0,0 +1,22 @@
+pH,charge,charge_error
+2.0,13.674,0.00692820323028
+2.5,12.891,0.00640312423743
+3.0,12.266,0.00469041575982
+3.5,11.788,0.00761577310586
+4.0,11.169,0.00721110255093
+4.5,10.167,0.00721110255093
+5.0,8.853,0.00714142842854
+5.5,7.447,0.00707106781187
+6.0,6.252,0.0636631761696
+6.5,5.416,0.0637024332345
+7.0,4.837,0.0637024332345
+7.5,4.278,0.0895321171424
+8.0,3.69,0.0895823643358
+8.5,3.112,0.089554452709
+9.0,2.501,0.089554452709
+9.5,1.729,0.0895377015564
+10.0,0.826,0.0895879456177
+10.5,-0.115,0.0634428877022
+11.0,-0.946,0.0634428877022
+11.5,-1.654,0.0634428877022
+12.0,-2.275,0.00721110255093

testsuite/data/Lunkad2021a.csv

@@ -0,0 +1,22 @@
+pH,charge,charge_error
+2.0,4.328030899613755,0.0538962451412445
+2.5,3.676594042574468,0.06369778710295225
+3.0,2.846670666616667,0.0691114121594015
+3.5,1.9688966387920153,0.06892040676649949
+4.0,1.1622067224159691,0.06312241346943875
+4.5,0.5701116236047055,0.05065347517351694
+5.0,0.22044849439381942,0.03501732616728689
+5.5,0.07649279384007635,0.021428102306625365
+6.0,0.024360945488182573,0.01246182548157354
+6.5,0.0072299096261296825,0.007105971896087083
+7.0,0.0009874876564035517,0.004790283141830334
+7.5,-0.003301208734892036,0.005269950644776093
+8.0,-0.012316096048799352,0.009043092692059515
+8.5,-0.04268946638167215,0.016169568639180957
+9.0,-0.12507093661329272,0.02716407000198067
+9.5,-0.3434894563817936,0.04259887933043685
+10.0,-0.8331858351770594,0.05840629764257237
+10.5,-1.612156098048775,0.06884781788167588
+11.0,-2.5490993862576716,0.07198577122196975
+11.5,-3.48038649516881,0.06908347786356787
+12.0,-4.236439544505694,0.05791528997068287

testsuite/data/Lunkad2021b.csv

@@ -0,0 +1,22 @@
+pH,charge,charge_error
+2.0,4.769976625292183,0.03589861247790809
+2.5,4.403126210922363,0.05183214139258131
+3.0,3.767890401369984,0.06473042658758646
+3.5,2.910399870001626,0.07297264097441067
+4.0,1.9509056136798293,0.076159435380892
+4.5,1.01215734803315,0.07623737867900791
+5.0,0.1300896238797007,0.07517492641612307
+5.5,-0.7458756765540429,0.07613015058623873
+6.0,-1.6825464681691473,0.07653887855278552
+6.5,-2.684907688653892,0.07382705658513389
+7.0,-3.5904288696391298,0.0687849849211043
+7.5,-4.312922338470769,0.05559908936108033
+8.0,-4.725088436394546,0.0387987480746436
+8.5,-4.9062336720791,0.02385030140762049
+9.0,-4.967755403057463,0.013989601304319406
+9.5,-4.989578880263996,0.008034723183410933
+10.0,-4.997462531718353,0.004017910847679971
+10.5,-4.998956263046711,0.0025471341039558855
+11.0,-4.999727503406207,0.001303259008174
+11.5,-4.999813752328096,0.0010779120839355
+12.0,-5.0,0.0