Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2021, Andrea Arsiccio
3 :
4 : This software is provided 'as-is', without any express or implied
5 : warranty. In no event will the authors be held liable for any damages
6 : arising from the use of this software.
7 :
8 : Permission is granted to anyone to use this software for any purpose,
9 : including commercial applications, and to alter it and redistribute it
10 : freely, subject to the following restrictions:
11 :
12 : 1. The origin of this software must not be misrepresented; you must not
13 : claim that you wrote the original software. If you use this software
14 : in a product, an acknowledgment in the product documentation would be
15 : appreciated but is not required.
16 : 2. Altered source versions must be plainly marked as such, and must not be
17 : misrepresented as being the original software.
18 : 3. This notice may not be removed or altered from any source distribution.
19 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
20 :
21 : #include "Sasa.h"
22 : #include "core/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/GenericMolInfo.h"
25 : #include "core/ActionSet.h"
26 : #include <cstdio>
27 : #include <iostream>
28 : #include <string>
29 : #include <cmath>
30 : #include <sstream>
31 : #include <algorithm>
32 : #include <iterator>
33 : #include <fstream>
34 :
35 :
36 : using namespace std;
37 :
38 : namespace PLMD {
39 : namespace sasa {
40 :
41 : //+PLUMEDOC COLVAR SASA_HASEL
42 : /*
43 : Calculates the solvent accessible surface area (SASA) of a protein molecule, or other properties related to it.
44 :
45 : The atoms for which the SASA is desired should be indicated with the keyword ATOMS, and a pdb file of the protein must be provided in input with the MOLINFO keyword. The algorithm described in the first paper cited below is used for the calculation. The radius of the solvent is assumed to be 0.14 nm, which is the radius of water molecules. Using the keyword NL_STRIDE it is also possible to specify the frequency with which the neighbor list for the calculation of SASA is updated (the default is every 10 steps).
46 :
47 : Different properties can be calculated and selected using the TYPE keyword:
48 :
49 : 1) the total SASA (TOTAL);
50 :
51 : 2) the free energy of transfer for the protein according to the transfer model (TRANSFER). This keyword can be used, for instance, to compute the transfer of a protein to different temperatures, as detailed in the second paper cited below, or to different pressures, as detailed in the fourth paper cited below, or to different osmolyte solutions, as detailed in the third paper cited below.
52 :
53 :
54 : When the TRANSFER keyword is used, a file with the free energy of transfer values for the sidechains and backbone atoms should be provided (using the keyword DELTAGFILE). Such file should have the following format:
55 :
56 : ````
57 : ----------------Sample DeltaG.dat file---------------------
58 : ALA 0.711019999999962
59 : ARG -2.24832799999996
60 : ASN -2.74838799999999
61 : ASP -2.5626376
62 : CYS 3.89864000000006
63 : GLN -1.76192
64 : GLU -2.38664400000002
65 : GLY 0
66 : HIS -3.58152799999999
67 : ILE 2.42634399999986
68 : LEU 1.77233599999988
69 : LYS -1.92576400000002
70 : MET -0.262827999999956
71 : PHE 1.62028800000007
72 : PRO -2.15598800000001
73 : SER -1.60934800000004
74 : THR -0.591559999999987
75 : TRP 1.22936000000027
76 : TYR 0.775547999999958
77 : VAL 2.12779200000011
78 : BACKBONE 1.00066920000002
79 : -----------------------------------------------------------
80 : ````
81 :
82 : where the second column is the free energy of transfer for each sidechain/backbone, in kJ/mol.
83 :
84 : A Python script for the computation of free energy of transfer values to describe the effect of osmolyte concentration, temperature and pressure (according to the second, third and fourth papers cited below) is freely available at
85 : [here](https://github.com/andrea-arsiccio/DeltaG-calculation). The script automatically outputs a DeltaG.dat file compatible with this SASA module.
86 :
87 :
88 : If the DELTAGFILE is not provided, the program computes the free energy of transfer values as if they had to take into account the effect of temperature according to approaches 2 or 3 in the paper the second paper cited below. Please read and cite this paper if using the transfer model for computing the effect of temperature in implicit solvent simulations. For this purpose, the keyword APPROACH should be added, and set to either 2 or 3.
89 :
90 : The SASA usually makes sense when atoms used for the calculation are all part of the same molecule. When running with periodic boundary conditions, the atoms should be in the proper periodic image. This is done automatically since PLUMED 2.2, by rebuilding the broken entities using a procedure that is equivalent to that done in [WHOLEMOLECULES](WHOLEMOLECULES.md). Notice that rebuilding is local to this action. This is different from [WHOLEMOLECULES](WHOLEMOLECULES.md) which actually modifies the coordinates stored in PLUMED.
91 :
92 : In case you want to recover the old behavior you should use the NOPBC flag.
93 : In that case you need to take care that atoms are in the correct periodic image.
94 :
95 : The SASA may also be computed using the SASA_LCPO collective variable, which makes use of the LCPO algorithm that is discussed in the last paper cited below. SASA_LCPO is more accurate then SASA_HASEL, but the computation is slower.
96 :
97 :
98 : ## Examples
99 :
100 : The following input tells plumed to print the total SASA for atoms 10 to 20 in a protein chain.
101 :
102 : ```plumed
103 : SASA_HASEL TYPE=TOTAL ATOMS=10-20 NL_STRIDE=10 LABEL=sasa
104 : PRINT ARG=sasa STRIDE=1 FILE=colvar
105 : ```
106 :
107 : The following input tells plumed to compute the transfer free energy for the protein chain containing atoms 10 to 20. Such transfer free energy is then used as a bias in the simulation (e.g., implicit solvent simulations). The free energy of transfer values are read from a file called DeltaG.dat.
108 :
109 : ```plumed
110 : SASA_HASEL TYPE=TRANSFER ATOMS=10-20 NL_STRIDE=10 DELTAGFILE=DeltaG.dat LABEL=sasa
111 :
112 : bias: BIASVALUE ARG=sasa
113 :
114 : PRINT ARG=sasa,bias.* STRIDE=1 FILE=colvar
115 : ```
116 :
117 : The following input tells plumed to compute the transfer free energy for the protein chain containing atoms 10 to 20. Such transfer free energy is then used as a bias in the simulation (e.g., implicit solvent simulations). The free energy of transfer values are computed according to the second paper cited below, and take into account the effect of temperature using approach 2 as described in the paper.
118 :
119 : ```plumed
120 : SASA_HASEL TYPE=TRANSFER ATOMS=10-20 NL_STRIDE=10 APPROACH=2 LABEL=sasa
121 :
122 : bias: BIASVALUE ARG=sasa
123 :
124 : PRINT ARG=sasa,bias.* STRIDE=1 FILE=colvar
125 : ```
126 :
127 : */
128 : //+ENDPLUMEDOC
129 :
130 : class SASA_HASEL : public Colvar {
131 : private:
132 : enum CV_TYPE {TOTAL, TRANSFER};
133 : int sasa_type;
134 : bool nopbc;
135 : double rs;
136 : string DeltaGValues;
137 : int approach;
138 : unsigned stride;
139 : unsigned nl_update;
140 : int firstStepFlag;
141 : double Ti;
142 : // cppcheck-suppress duplInheritedMember
143 : std::vector<AtomNumber> atoms;
144 : vector < vector < std::string > > AtomResidueName;
145 : vector < vector < double > > SASAparam;
146 : vector < vector < std::string > > CONNECTparam;
147 : unsigned natoms;
148 : vector < vector < double > > MaxSurf;
149 : vector < vector < double > > DeltaG;
150 : vector < vector < int > > Nlist;
151 : public:
152 : static void registerKeywords(Keywords& keys);
153 : explicit SASA_HASEL(const ActionOptions&);
154 : void readPDB();
155 : map<string, vector<std::string> > setupHASELparam();
156 : void readSASAparam();
157 : void calcNlist();
158 : map<string, vector<double> > setupMaxSurfMap();
159 : void readMaxSurf();
160 : void readDeltaG();
161 : void computeDeltaG();
162 : virtual void calculate();
163 : };
164 :
165 : PLUMED_REGISTER_ACTION(SASA_HASEL,"SASA_HASEL")
166 :
167 6 : void SASA_HASEL::registerKeywords(Keywords& keys) {
168 6 : Colvar::registerKeywords(keys);
169 6 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the SASA for");
170 6 : keys.add("compulsory","TYPE","TOTAL","The type of calculation you want to perform. Can be TOTAL or TRANSFER");
171 6 : keys.add("compulsory", "NL_STRIDE", "The frequency with which the neighbor list for the calculation of SASA is updated.");
172 6 : keys.add("optional","DELTAGFILE","a file containing the free energy of transfer values for backbone and sidechains atoms. Necessary only if TYPE = TRANSFER. A Python script for the computation of free energy of transfer values to describe the effect of osmolyte concentration, temperature and pressure is freely available at https://github.com/andrea-arsiccio/DeltaG-calculation. The script automatically outputs a DeltaG.dat file compatible with this SASA module. If TYPE = TRANSFER and no DELTAGFILE is provided, the free energy values are those describing the effect of temperature, and are computed using the temperature value passed by the MD engine");
173 6 : keys.add("optional","APPROACH","either approach 2 or 3. Necessary only if TYPE = TRANSFER and no DELTAGFILE is provided. If TYPE = TRANSFER and no DELTAGFILE is provided, the free energy values are those describing the effect of temperature, and the program must know if approach 2 or 3 (as described in Arsiccio and Shea, Protein Cold Denaturation in Implicit Solvent Simulations: A Transfer Free Energy Approach, J. Phys. Chem. B, 2021) needs to be used to compute them");
174 12 : keys.setValueDescription("scalar","the solvent accessible surface area (SASA) of the molecule");
175 6 : keys.addDOI("10.1016/0898-5529(88)90015-2");
176 6 : keys.addDOI("10.1021/acs.jpcb.1c01694");
177 6 : keys.addDOI("10.1021/acs.jpcb.2c00889");
178 6 : keys.addDOI("10.1021/acs.jpcb.1c04398");
179 6 : keys.addDOI("10.1002/(SICI)1096-987X(19990130)20:2<217::AID-JCC4>3.0.CO;2-A");
180 6 : }
181 :
182 :
183 2 : SASA_HASEL::SASA_HASEL(const ActionOptions&ao):
184 : PLUMED_COLVAR_INIT(ao),
185 2 : nopbc(false),
186 4 : DeltaGValues("absent"),
187 2 : stride(10),
188 2 : nl_update(0),
189 2 : firstStepFlag(0),
190 2 : Ti(0) {
191 2 : rs = 0.14;
192 2 : parse("DELTAGFILE",DeltaGValues);
193 2 : parse("APPROACH", approach);
194 4 : parseAtomList("ATOMS",atoms);
195 2 : if(atoms.size()==0) {
196 0 : error("no atoms specified");
197 : }
198 : std::string Type;
199 2 : parse("TYPE",Type);
200 2 : parse("NL_STRIDE", stride);
201 2 : parseFlag("NOPBC",nopbc);
202 2 : checkRead();
203 :
204 2 : if(Type=="TOTAL") {
205 1 : sasa_type=TOTAL;
206 1 : } else if(Type=="TRANSFER") {
207 1 : sasa_type=TRANSFER;
208 : } else {
209 0 : error("Unknown SASA type");
210 : }
211 :
212 2 : switch(sasa_type) {
213 1 : case TOTAL:
214 1 : log.printf(" TOTAL SASA;");
215 : break;
216 1 : case TRANSFER:
217 1 : log.printf(" TRANSFER MODEL;");
218 : break;
219 : }
220 :
221 2 : log.printf(" atoms involved : ");
222 242 : for(unsigned i=0; i<atoms.size(); ++i) {
223 240 : if(i%25==0) {
224 10 : log<<"\n";
225 : }
226 240 : log.printf("%d ",atoms[i].serial());
227 : }
228 2 : log.printf("\n");
229 :
230 2 : if(nopbc) {
231 0 : log<<" PBC will be ignored\n";
232 : } else {
233 2 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
234 : }
235 :
236 :
237 2 : addValueWithDerivatives();
238 2 : setNotPeriodic();
239 2 : requestAtoms(atoms);
240 :
241 2 : natoms = getNumberOfAtoms();
242 2 : AtomResidueName.resize(2);
243 2 : SASAparam.resize(natoms);
244 2 : CONNECTparam.resize(natoms);
245 2 : MaxSurf.resize(natoms);
246 2 : DeltaG.resize(natoms+1);
247 2 : Nlist.resize(natoms);
248 :
249 :
250 2 : }
251 :
252 :
253 : //splits strings into tokens. Used to read into SASA parameters file and into reference pdb file
254 : template <class Container>
255 42 : void split(const std::string& str, Container& cont) {
256 42 : std::istringstream iss(str);
257 84 : std::copy(std::istream_iterator<std::string>(iss),
258 42 : std::istream_iterator<std::string>(),
259 : std::back_inserter(cont));
260 42 : }
261 :
262 :
263 : //reads input PDB file provided with MOLINFO, and assigns atom and residue names to each atom involved in the CV
264 :
265 2 : void SASA_HASEL::readPDB() {
266 2 : auto* moldat = plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
267 2 : if( ! moldat ) {
268 0 : error("Unable to find MOLINFO in input");
269 : }
270 2 : AtomResidueName[0].clear();
271 2 : AtomResidueName[1].clear();
272 :
273 242 : for(unsigned i=0; i<natoms; i++) {
274 240 : string Aname = moldat->getAtomName(atoms[i]);
275 240 : string Rname = moldat->getResidueName(atoms[i]);
276 240 : AtomResidueName[0].push_back (Aname);
277 240 : AtomResidueName[1].push_back (Rname);
278 : }
279 :
280 2 : }
281 :
282 :
283 : //Hasel et al. parameters database
284 2 : map<string, vector<std::string> > SASA_HASEL::setupHASELparam() {
285 : map<string, vector<std::string> > haselmap;
286 16 : haselmap["ALA_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
287 16 : haselmap["ALA_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
288 16 : haselmap["ALA_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
289 16 : haselmap["ALA_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
290 16 : haselmap["ALA_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
291 16 : haselmap["ALA_CB"] = { "2.0", "0.88", "CA", "Z", "Z", "Z", };
292 16 : haselmap["ASP_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
293 16 : haselmap["ASP_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
294 16 : haselmap["ASP_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
295 16 : haselmap["ASP_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
296 16 : haselmap["ASP_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
297 16 : haselmap["ASP_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
298 16 : haselmap["ASP_CG"] = { "1.72", "1.554", "CB", "OD1", "OD2", "Z", };
299 16 : haselmap["ASP_OD1"] = { "1.5", "0.926", "CG", "Z", "Z", "Z", };
300 16 : haselmap["ASP_OD2"] = { "1.7", "0.922", "CG", "Z", "Z", "Z", };
301 16 : haselmap["ASN_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
302 16 : haselmap["ASN_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
303 16 : haselmap["ASN_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
304 16 : haselmap["ASN_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
305 16 : haselmap["ASN_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
306 16 : haselmap["ASN_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
307 16 : haselmap["ASN_CG"] = { "1.7", "2.149", "CB", "OD1", "ND2", "Z", };
308 16 : haselmap["ASN_OD1"] = { "1.5", "0.926", "CG", "Z", "Z", "Z", };
309 16 : haselmap["ASN_ND2"] = { "1.6", "1.215", "CG", "1HD2", "1HD2", "Z", };
310 16 : haselmap["ASN_1HD2"] = { "1.1", "1.128", "ND2", "Z", "Z", "Z", };
311 16 : haselmap["ASN_2HD2"] = { "1.1", "1.128", "ND2", "Z", "Z", "Z", };
312 16 : haselmap["ARG_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
313 16 : haselmap["ARG_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
314 16 : haselmap["ARG_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
315 16 : haselmap["ARG_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
316 16 : haselmap["ARG_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
317 16 : haselmap["ARG_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
318 16 : haselmap["ARG_CG"] = { "1.9", "1.045", "CB", "CD", "Z", "Z", };
319 16 : haselmap["ARG_CD"] = { "1.9", "1.045", "CG", "NE", "Z", "Z", };
320 16 : haselmap["ARG_NE"] = { "1.55", "1.028", "CD", "HE", "CZ", "Z", };
321 16 : haselmap["ARG_NH1"] = { "1.55", "1.028", "CZ", "1HH1", "2HH1", "Z", };
322 16 : haselmap["ARG_NH2"] = { "1.55", "1.028", "CZ", "1HH2", "2HH2", "Z", };
323 16 : haselmap["ARG_CZ"] = { "1.72", "1.554", "NE", "NH1", "NH2", "Z", };
324 16 : haselmap["ARG_HE"] = { "1.1", "1.128", "NE", "Z", "Z", "Z", };
325 16 : haselmap["ARG_1HH2"] = { "1.1", "1.128", "NH2", "Z", "Z", "Z", };
326 16 : haselmap["ARG_2HH2"] = { "1.1", "1.128", "NH2", "Z", "Z", "Z", };
327 16 : haselmap["ARG_2HH1"] = { "1.1", "1.128", "NH1", "Z", "Z", "Z", };
328 16 : haselmap["ARG_1HH1"] = { "1.1", "1.128", "NH1", "Z", "Z", "Z", };
329 16 : haselmap["CYS_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
330 16 : haselmap["CYS_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
331 16 : haselmap["CYS_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
332 16 : haselmap["CYS_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
333 16 : haselmap["CYS_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
334 16 : haselmap["CYS_CB"] = { "1.9", "1.045", "CA", "SG", "Z", "Z", };
335 16 : haselmap["CYS_SG"] = { "1.8", "1.121", "CB", "HG", "Z", "Z", };
336 16 : haselmap["CYS_HG"] = { "1.2", "0.928", "SG", "Z", "Z", "Z", };
337 16 : haselmap["GLU_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
338 16 : haselmap["GLU_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
339 16 : haselmap["GLU_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
340 16 : haselmap["GLU_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
341 16 : haselmap["GLU_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
342 16 : haselmap["GLU_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
343 16 : haselmap["GLU_CG"] = { "1.9", "1.045", "CB", "CD", "Z", "Z", };
344 16 : haselmap["GLU_CD"] = { "1.72", "1.554", "CG", "OE1", "OE2", "Z", };
345 16 : haselmap["GLU_OE1"] = { "1.5", "0.926", "CD", "Z", "Z", "Z", };
346 16 : haselmap["GLU_OE2"] = { "1.7", "0.922", "CD", "Z", "Z", "Z", };
347 16 : haselmap["GLN_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
348 16 : haselmap["GLN_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
349 16 : haselmap["GLN_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
350 16 : haselmap["GLN_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
351 16 : haselmap["GLN_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
352 16 : haselmap["GLN_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
353 16 : haselmap["GLN_CG"] = { "1.9", "1.045", "CB", "CD", "Z", "Z", };
354 16 : haselmap["GLN_CD"] = { "1.72", "1.554", "CG", "OE1", "NE2", "Z", };
355 16 : haselmap["GLN_OE1"] = { "1.5", "0.926", "CD", "Z", "Z", "Z", };
356 16 : haselmap["GLN_NE2"] = { "1.6", "1.215", "CD", "2HE2", "1HE2", "Z", };
357 16 : haselmap["GLN_2HE2"] = { "1.1", "1.128", "NE2", "Z", "Z", "Z", };
358 16 : haselmap["GLN_1HE2"] = { "1.1", "1.128", "NE2", "Z", "Z", "Z", };
359 16 : haselmap["GLY_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
360 16 : haselmap["GLY_CA"] = { "1.7", "2.149", "N", "C", "Z", "Z", };
361 16 : haselmap["GLY_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
362 16 : haselmap["GLY_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
363 16 : haselmap["GLY_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
364 16 : haselmap["HIS_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
365 16 : haselmap["HIS_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
366 16 : haselmap["HIS_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
367 16 : haselmap["HIS_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
368 16 : haselmap["HIS_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
369 16 : haselmap["HIS_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
370 16 : haselmap["HIS_CG"] = { "1.72", "1.554", "CB", "CD2", "ND1", "Z", };
371 16 : haselmap["HIS_ND1"] = { "1.55", "1.028", "CG", "CE1", "Z", "Z", };
372 16 : haselmap["HIS_CE1"] = { "1.8", "1.073", "ND1", "NE2", "Z", "Z", };
373 16 : haselmap["HIS_NE2"] = { "1.55", "1.413", "CE1", "2HE", "CD2", "Z", };
374 16 : haselmap["HIS_CD2"] = { "1.8", "1.073", "NE2", "CG", "Z", "Z", };
375 16 : haselmap["HIS_2HE"] = { "1.1", "1.128", "NE2", "Z", "Z", "Z", };
376 16 : haselmap["ILE_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
377 16 : haselmap["ILE_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
378 16 : haselmap["ILE_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
379 16 : haselmap["ILE_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
380 16 : haselmap["ILE_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
381 16 : haselmap["ILE_CB"] = { "1.8", "1.276", "CA", "CG2", "CG1", "Z", };
382 16 : haselmap["ILE_CG2"] = { "2.0", "0.88", "CB", "Z", "Z", "Z", };
383 16 : haselmap["ILE_CG1"] = { "1.9", "1.045", "CB", "CD1", "Z", "Z", };
384 16 : haselmap["ILE_CD1"] = { "2.0", "0.88", "CG1", "Z", "Z", "Z", };
385 16 : haselmap["LEU_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
386 16 : haselmap["LEU_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
387 16 : haselmap["LEU_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
388 16 : haselmap["LEU_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
389 16 : haselmap["LEU_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
390 16 : haselmap["LEU_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
391 16 : haselmap["LEU_CG"] = { "1.8", "1.276", "CB", "CD1", "CD2", "Z", };
392 16 : haselmap["LEU_CD1"] = { "2.0", "0.88", "CG", "Z", "Z", "Z", };
393 16 : haselmap["LEU_CD2"] = { "2.0", "0.88", "CG", "Z", "Z", "Z", };
394 16 : haselmap["LYS_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
395 16 : haselmap["LYS_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
396 16 : haselmap["LYS_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
397 16 : haselmap["LYS_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
398 16 : haselmap["LYS_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
399 16 : haselmap["LYS_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
400 16 : haselmap["LYS_CG"] = { "1.9", "1.045", "CB", "CD", "Z", "Z", };
401 16 : haselmap["LYS_CD"] = { "1.9", "1.045", "CG", "CE", "Z", "Z", };
402 16 : haselmap["LYS_CE"] = { "1.9", "1.045", "CD", "NZ", "Z", "Z", };
403 16 : haselmap["LYS_NZ"] = { "1.6", "1.215", "CE", "1HZ", "2HZ", "3HZ", };
404 16 : haselmap["LYS_1HZ"] = { "1.1", "1.128", "NZ", "Z", "Z", "Z", };
405 16 : haselmap["LYS_2HZ"] = { "1.1", "1.128", "NZ", "Z", "Z", "Z", };
406 16 : haselmap["LYS_3HZ"] = { "1.1", "1.128", "NZ", "Z", "Z", "Z", };
407 16 : haselmap["MET_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
408 16 : haselmap["MET_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
409 16 : haselmap["MET_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
410 16 : haselmap["MET_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
411 16 : haselmap["MET_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
412 16 : haselmap["MET_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
413 16 : haselmap["MET_CG"] = { "1.9", "1.045", "CB", "SD", "Z", "Z", };
414 16 : haselmap["MET_SD"] = { "1.8", "1.121", "CG", "CE", "Z", "Z", };
415 16 : haselmap["MET_CE"] = { "2.0", "0.88", "SD", "Z", "Z", "Z", };
416 16 : haselmap["PHE_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
417 16 : haselmap["PHE_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
418 16 : haselmap["PHE_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
419 16 : haselmap["PHE_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
420 16 : haselmap["PHE_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
421 16 : haselmap["PHE_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
422 16 : haselmap["PHE_CG"] = { "1.72", "1.554", "CB", "CD1", "CD2", "Z", };
423 16 : haselmap["PHE_CD1"] = { "1.8", "1.073", "CG", "CE1", "Z", "Z", };
424 16 : haselmap["PHE_CE1"] = { "1.8", "1.073", "CD1", "CZ", "Z", "Z", };
425 16 : haselmap["PHE_CZ"] = { "1.8", "1.073", "CE1", "CE2", "Z", "Z", };
426 16 : haselmap["PHE_CE2"] = { "1.8", "1.073", "CZ", "CD2", "Z", "Z", };
427 16 : haselmap["PHE_CD2"] = { "1.8", "1.073", "CE2", "CG", "Z", "Z", };
428 16 : haselmap["PRO_N"] = { "1.55", "1.028", "CD", "CA", "Z", "Z", };
429 16 : haselmap["PRO_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
430 16 : haselmap["PRO_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
431 16 : haselmap["PRO_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
432 16 : haselmap["PRO_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
433 16 : haselmap["PRO_CG"] = { "1.9", "1.045", "CB", "CD", "Z", "Z", };
434 16 : haselmap["PRO_CD"] = { "1.9", "1.045", "CG", "N", "Z", "Z", };
435 16 : haselmap["SER_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
436 16 : haselmap["SER_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
437 16 : haselmap["SER_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
438 16 : haselmap["SER_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
439 16 : haselmap["SER_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
440 16 : haselmap["SER_CB"] = { "1.9", "1.045", "CA", "OG", "Z", "Z", };
441 16 : haselmap["SER_OG"] = { "1.52", "1.08", "CB", "HG", "Z", "Z", };
442 16 : haselmap["SER_HG"] = { "1.0", "0.944", "OG", "Z", "Z", "Z", };
443 16 : haselmap["THR_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
444 16 : haselmap["THR_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
445 16 : haselmap["THR_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
446 16 : haselmap["THR_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
447 16 : haselmap["THR_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
448 16 : haselmap["THR_CB"] = { "1.8", "1.276", "CA", "CG2", "OG1", "Z", };
449 16 : haselmap["THR_CG2"] = { "2.0", "0.88", "CB", "Z", "Z", "Z", };
450 16 : haselmap["THR_OG1"] = { "1.52", "1.08", "1HG", "CB", "Z", "Z", };
451 16 : haselmap["THR_1HG"] = { "1.0", "0.944", "OG1", "Z", "Z", "Z", };
452 16 : haselmap["TRP_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
453 16 : haselmap["TRP_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
454 16 : haselmap["TRP_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
455 16 : haselmap["TRP_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
456 16 : haselmap["TRP_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
457 16 : haselmap["TRP_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
458 16 : haselmap["TRP_CG"] = { "1.72", "1.554", "CB", "CD2", "CD1", "Z", };
459 16 : haselmap["TRP_CD1"] = { "1.8", "1.073", "CG", "NE1", "Z", "Z", };
460 16 : haselmap["TRP_NE1"] = { "1.55", "1.413", "CD1", "CE2", "1HE", "Z", };
461 16 : haselmap["TRP_CE2"] = { "1.72", "1.554", "NE1", "CD2", "CZ2", "Z", };
462 16 : haselmap["TRP_CZ2"] = { "1.8", "1.073", "CE2", "CH2", "Z", "Z", };
463 16 : haselmap["TRP_CH2"] = { "1.8", "1.073", "CZ2", "CZ3", "Z", "Z", };
464 16 : haselmap["TRP_CZ3"] = { "1.8", "1.073", "CH2", "CE3", "Z", "Z", };
465 16 : haselmap["TRP_CE3"] = { "1.8", "1.073", "CZ3", "CD2", "Z", "Z", };
466 16 : haselmap["TRP_CD2"] = { "1.72", "1.554", "CE3", "CE2", "CG", "Z", };
467 16 : haselmap["TRP_1HE"] = { "1.1", "1.128", "NE1", "Z", "Z", "Z", };
468 16 : haselmap["TYR_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
469 16 : haselmap["TYR_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
470 16 : haselmap["TYR_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
471 16 : haselmap["TYR_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
472 16 : haselmap["TYR_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
473 16 : haselmap["TYR_CB"] = { "1.9", "1.045", "CA", "CG", "Z", "Z", };
474 16 : haselmap["TYR_CG"] = { "1.72", "1.554", "CB", "CD1", "CD2", "Z", };
475 16 : haselmap["TYR_CD1"] = { "1.8", "1.073", "CG", "CE1", "Z", "Z", };
476 16 : haselmap["TYR_CE1"] = { "1.8", "1.073", "CD1", "CZ", "Z", "Z", };
477 16 : haselmap["TYR_CZ"] = { "1.72", "1.554", "CE1", "OH", "CE2", "Z", };
478 16 : haselmap["TYR_OH"] = { "1.52", "1.08", "CZ", "HH", "Z", "Z", };
479 16 : haselmap["TYR_HH"] = { "1.0", "0.944", "OH", "Z", "Z", "Z", };
480 16 : haselmap["TYR_CE2"] = { "1.8", "1.073", "CZ", "CD2", "Z", "Z", };
481 16 : haselmap["TYR_CD2"] = { "1.8", "1.073", "CE2", "CG", "Z", "Z", };
482 16 : haselmap["VAL_N"] = { "1.55", "1.028", "H", "CA", "Z", "Z", };
483 16 : haselmap["VAL_CA"] = { "1.7", "2.149", "N", "C", "CB", "Z", };
484 16 : haselmap["VAL_C"] = { "1.72", "1.554", "CA", "O", "Z", "Z", };
485 16 : haselmap["VAL_O"] = { "1.5", "0.926", "C", "Z", "Z", "Z", };
486 16 : haselmap["VAL_H"] = { "1.1", "1.128", "N", "Z", "Z", "Z", };
487 16 : haselmap["VAL_CB"] = { "1.8", "1.276", "CA", "CG1", "CG2", "Z", };
488 16 : haselmap["VAL_CG1"] = { "2.0", "0.88", "CB", "Z", "Z", "Z", };
489 16 : haselmap["VAL_CG2"] = { "2.0", "0.88", "CB", "Z", "Z", "Z", };
490 2 : return haselmap;
491 410 : }
492 :
493 : //assigns SASA parameters to each atom reading from HASEL parameter database
494 2 : void SASA_HASEL::readSASAparam() {
495 :
496 242 : for(unsigned i=0; i<natoms; i++) {
497 240 : SASAparam[i].clear();
498 240 : CONNECTparam[i].clear();
499 : }
500 :
501 : map<string, vector<std::string> > haselmap;
502 2 : haselmap = setupHASELparam();
503 2 : vector<std::string> HASELparamVector;
504 : string identifier;
505 :
506 :
507 242 : for(unsigned i=0; i<natoms; i++) {
508 480 : identifier = AtomResidueName[1][i]+"_"+AtomResidueName[0][i];
509 240 : if (haselmap.find(identifier)!=haselmap.end()) {
510 144 : HASELparamVector = haselmap.at(identifier);
511 144 : SASAparam[i].push_back (std::atof(HASELparamVector[0].c_str())+rs*10);
512 144 : SASAparam[i].push_back (std::atof(HASELparamVector[1].c_str()));
513 288 : CONNECTparam[i].push_back (HASELparamVector[2].c_str());
514 288 : CONNECTparam[i].push_back (HASELparamVector[3].c_str());
515 288 : CONNECTparam[i].push_back (HASELparamVector[4].c_str());
516 288 : CONNECTparam[i].push_back (HASELparamVector[5].c_str());
517 : }
518 : }
519 :
520 :
521 242 : for(unsigned i=0; i<natoms; i++) {
522 240 : if (SASAparam[i].size()==0 ) {
523 96 : if ((AtomResidueName[0][i][0]=='O') || (AtomResidueName[0][i][0]=='N') || (AtomResidueName[0][i][0]=='C') || (AtomResidueName[0][i][0]=='S')) {
524 0 : cout << "Could not find SASA parameters for atom " << AtomResidueName[0][i] << " of residue " << AtomResidueName[1][i] << endl;
525 0 : error ("missing SASA parameters\n");
526 : }
527 : }
528 : }
529 :
530 :
531 4 : }
532 :
533 :
534 :
535 : //Max Surf values, used only if TYPE=TRANSFER
536 1 : map<string, vector<double> > SASA_HASEL::setupMaxSurfMap() {
537 : // Max Surface Area for backbone and sidechain, in nm2
538 : map<string, vector<double> > valuemap;
539 1 : valuemap ["ALA"]= {0.56425,0.584851,};
540 1 : valuemap ["ARG"]= {0.498656,1.808093,};
541 1 : valuemap ["ASN"]= {0.473409,0.818394,};
542 1 : valuemap ["ASP"]= {0.477057,0.977303,};
543 1 : valuemap ["CYS"]= {0.507512,0.791483,};
544 1 : valuemap ["GLN"]= {0.485859,1.281534,};
545 1 : valuemap ["GLU"]= {0.495054,1.464718,};
546 1 : valuemap ["GLY"]= {0.658632,0,};
547 1 : valuemap ["HIS"]= {0.48194,1.118851,};
548 1 : valuemap ["ILE"]= {0.461283,1.450569,};
549 1 : valuemap ["LEU"]= {0.476315,1.498843,};
550 1 : valuemap ["LYS"]= {0.493533,1.619731,};
551 1 : valuemap ["MET"]= {0.507019,1.631904,};
552 1 : valuemap ["PHE"]= {0.457462, 1.275125,};
553 1 : valuemap ["PRO"]= {0.315865,0.859456,};
554 1 : valuemap ["SER"]= {0.48636,0.627233,};
555 1 : valuemap ["THR"]= {0.45064,0.91088,};
556 1 : valuemap ["TRP"]= {0.45762,1.366369,};
557 1 : valuemap ["TYR"]= {0.461826,1.425822,};
558 1 : valuemap ["VAL"]= {0.477054,1.149101,};
559 1 : return valuemap;
560 : }
561 :
562 :
563 :
564 : //reads maximum surface values per residue type and assigns values to each atom, only used if sasa_type = TRANSFER
565 :
566 1 : void SASA_HASEL::readMaxSurf() {
567 : map<string, vector<double> > valuemap;
568 1 : valuemap = setupMaxSurfMap();
569 1 : vector<double> MaxSurfVector;
570 :
571 121 : for(unsigned i=0; i<natoms; i++) {
572 120 : MaxSurf[i].clear();
573 120 : MaxSurfVector = valuemap.at(AtomResidueName[1][i]);
574 120 : MaxSurf[i].push_back (MaxSurfVector[0]*100);
575 120 : MaxSurf[i].push_back (MaxSurfVector[1]*100);
576 : }
577 2 : }
578 :
579 : //reads file with free energy values for sidechains and for the backbone, and assigns values to each atom. Only used if sasa_type = TRANSFER
580 :
581 1 : void SASA_HASEL::readDeltaG() {
582 :
583 121 : for(unsigned i=0; i<natoms; i++) {
584 120 : DeltaG[i].clear();
585 : }
586 :
587 : string DeltaGline;
588 1 : fstream DeltaGFile;
589 1 : DeltaGFile.open(DeltaGValues);
590 1 : if (DeltaGFile) {
591 : int backboneflag = 0;
592 23 : while(getline(DeltaGFile, DeltaGline)) {
593 22 : if(!DeltaGline.empty()) {
594 21 : std::vector<std::string> DeltaGtoken;
595 21 : split(DeltaGline, DeltaGtoken);
596 2541 : for(unsigned i=0; i<natoms; i++) {
597 2520 : if (DeltaGtoken[0].compare(AtomResidueName[1][i])==0 ) {
598 120 : DeltaG[i].push_back (std::atof(DeltaGtoken[1].c_str()));
599 : }
600 : }
601 21 : if (DeltaGtoken[0].compare("BACKBONE")==0 ) {
602 : backboneflag = 1;
603 1 : DeltaG[natoms].push_back (std::atof(DeltaGtoken[1].c_str()));
604 : }
605 21 : }
606 : }
607 1 : if ( backboneflag == 0) {
608 0 : error("Cannot find backbone value in Delta G parameters file\n");
609 : }
610 : } else {
611 0 : error("Cannot open DeltaG file");
612 : }
613 :
614 121 : for(unsigned i=0; i<natoms; i++) {
615 120 : if (DeltaG[i].size()==0 ) {
616 0 : cout << "Delta G value for residue " << AtomResidueName[1][i] << " not found " << endl;
617 0 : error ("missing Delta G parameter\n");
618 : }
619 : }
620 :
621 2 : }
622 :
623 : //computes free energy values for the sidechains and for the backbone, and assigns values to each atom. Only used if sasa_type = TRANSFER, and if no DELTAGFILE is provided. In this case, the free energy values are those describing the effect of temperature, and the program must know if approach 2 or 3 (as described in Arsiccio and Shea, Protein Cold Denaturation in Implicit Solvent Simulations: A Transfer Free Energy Approach, JPCB, 2021) needs to be used to compute them.
624 :
625 0 : void SASA_HASEL::computeDeltaG() {
626 :
627 0 : for(unsigned i=0; i<natoms; i++) {
628 0 : DeltaG[i].clear();
629 : }
630 :
631 : double T;
632 0 : T = getkBT()/getKBoltzmann();
633 :
634 0 : if (T != Ti) {
635 0 : for(unsigned i=0; i<natoms; i++) {
636 0 : if (approach==2) {
637 0 : if (AtomResidueName[1][i].compare("ALA")==0) {
638 0 : DeltaG[i].push_back (-2.995/1000*std::pow(T,2)+1.808*T-272.895);
639 : }
640 0 : if (AtomResidueName[1][i].compare("ARG")==0) {
641 0 : DeltaG[i].push_back (-3.182/1000*std::pow(T,2)+1.894*T-282.032);
642 : }
643 0 : if (AtomResidueName[1][i].compare("ASN")==0) {
644 0 : DeltaG[i].push_back (-1.047/1000*std::pow(T,2)+0.6068*T-87.846);
645 : }
646 0 : if (AtomResidueName[1][i].compare("ASP")==0) {
647 0 : DeltaG[i].push_back (-0.1794/1000*std::pow(T,2)+0.1091*T-16.526);
648 : }
649 0 : if (AtomResidueName[1][i].compare("CYS")==0) {
650 0 : DeltaG[i].push_back (-3.09/1000*std::pow(T,2)+1.835*T-272.26);
651 : }
652 0 : if (AtomResidueName[1][i].compare("GLN")==0) {
653 0 : DeltaG[i].push_back (-2.23/1000*std::pow(T,2)+1.335*T-199.707);
654 : }
655 0 : if (AtomResidueName[1][i].compare("GLU")==0) {
656 0 : DeltaG[i].push_back (-1.511/1000*std::pow(T,2)+0.8904*T-131.168);
657 : }
658 0 : if (AtomResidueName[1][i].compare("GLY")==0) {
659 0 : DeltaG[i].push_back (0);
660 : }
661 0 : if (AtomResidueName[1][i].compare("HIS")==0) {
662 0 : DeltaG[i].push_back (-3.482/1000*std::pow(T,2)+2.084*T-311.694);
663 : }
664 0 : if (AtomResidueName[1][i].compare("ILE")==0) {
665 0 : DeltaG[i].push_back (-6.364/1000*std::pow(T,2)+3.8*T-567.444);
666 : }
667 0 : if (AtomResidueName[1][i].compare("LEU")==0) {
668 0 : DeltaG[i].push_back (-7.466/1000*std::pow(T,2)+4.417*T-653.394);
669 : }
670 0 : if (AtomResidueName[1][i].compare("LYS")==0) {
671 0 : DeltaG[i].push_back (-2.091/1000*std::pow(T,2)+1.2458*T-185.549);
672 : }
673 0 : if (AtomResidueName[1][i].compare("MET")==0) {
674 0 : DeltaG[i].push_back (-3.807/1000*std::pow(T,2)+2.272*T-339.007);
675 : }
676 0 : if (AtomResidueName[1][i].compare("PHE")==0) {
677 0 : DeltaG[i].push_back (-7.828/1000*std::pow(T,2)+4.644*T-688.874);
678 : }
679 0 : if (AtomResidueName[1][i].compare("PRO")==0) {
680 0 : DeltaG[i].push_back (-2.347/1000*std::pow(T,2)+1.411*T-212.059);
681 : }
682 0 : if (AtomResidueName[1][i].compare("SER")==0) {
683 0 : DeltaG[i].push_back (1.813/1000*std::pow(T,2)-1.05*T+151.957);
684 : }
685 0 : if (AtomResidueName[1][i].compare("THR")==0) {
686 0 : DeltaG[i].push_back (-2.64/1000*std::pow(T,2)+1.591*T-239.516);
687 : }
688 0 : if (AtomResidueName[1][i].compare("TRP")==0) {
689 0 : DeltaG[i].push_back (-11.66/1000*std::pow(T,2)+6.916*T-1025.293);
690 : }
691 0 : if (AtomResidueName[1][i].compare("TYR")==0) {
692 0 : DeltaG[i].push_back (-7.513/1000*std::pow(T,2)+4.478*T-667.261);
693 : }
694 0 : if (AtomResidueName[1][i].compare("VAL")==0) {
695 0 : DeltaG[i].push_back (-4.902/1000*std::pow(T,2)+2.921*T-435.309);
696 : }
697 0 : DeltaG[natoms].push_back (-0.6962/1000*std::pow(T,2)+0.4426*T-70.015);
698 : }
699 0 : if (approach==3) {
700 0 : if (AtomResidueName[1][i].compare("ALA")==0) {
701 0 : DeltaG[i].push_back (-2.995/1000*std::pow(T,2)+1.808*T-272.105);
702 : }
703 0 : if (AtomResidueName[1][i].compare("ARG")==0) {
704 0 : DeltaG[i].push_back (-3.182/1000*std::pow(T,2)+1.894*T-284.086);
705 : }
706 0 : if (AtomResidueName[1][i].compare("ASN")==0) {
707 0 : DeltaG[i].push_back (-1.047/1000*std::pow(T,2)+0.6068*T-90.597);
708 : }
709 0 : if (AtomResidueName[1][i].compare("ASP")==0) {
710 0 : DeltaG[i].push_back (-0.1794/1000*std::pow(T,2)+0.1091*T-19.143);
711 : }
712 0 : if (AtomResidueName[1][i].compare("CYS")==0) {
713 0 : DeltaG[i].push_back (-3.09/1000*std::pow(T,2)+1.835*T-268.527);
714 : }
715 0 : if (AtomResidueName[1][i].compare("GLN")==0) {
716 0 : DeltaG[i].push_back (-2.23/1000*std::pow(T,2)+1.335*T-201.559);
717 : }
718 0 : if (AtomResidueName[1][i].compare("GLU")==0) {
719 0 : DeltaG[i].push_back (-1.511/1000*std::pow(T,2)+0.8904*T-133.543);
720 : }
721 0 : if (AtomResidueName[1][i].compare("GLY")==0) {
722 0 : DeltaG[i].push_back (0);
723 : }
724 0 : if (AtomResidueName[1][i].compare("HIS")==0) {
725 0 : DeltaG[i].push_back (-3.482/1000*std::pow(T,2)+2.084*T-315.398);
726 : }
727 0 : if (AtomResidueName[1][i].compare("ILE")==0) {
728 0 : DeltaG[i].push_back (-6.364/1000*std::pow(T,2)+3.8*T-564.825);
729 : }
730 0 : if (AtomResidueName[1][i].compare("LEU")==0) {
731 0 : DeltaG[i].push_back (-7.466/1000*std::pow(T,2)+4.417*T-651.483);
732 : }
733 0 : if (AtomResidueName[1][i].compare("LYS")==0) {
734 0 : DeltaG[i].push_back (-2.091/1000*std::pow(T,2)+1.2458*T-187.485);
735 : }
736 0 : if (AtomResidueName[1][i].compare("MET")==0) {
737 0 : DeltaG[i].push_back (-3.807/1000*std::pow(T,2)+2.272*T-339.242);
738 : }
739 0 : if (AtomResidueName[1][i].compare("PHE")==0) {
740 0 : DeltaG[i].push_back (-7.828/1000*std::pow(T,2)+4.644*T-687.134);
741 : }
742 0 : if (AtomResidueName[1][i].compare("PRO")==0) {
743 0 : DeltaG[i].push_back (-2.347/1000*std::pow(T,2)+1.411*T-214.211);
744 : }
745 0 : if (AtomResidueName[1][i].compare("SER")==0) {
746 0 : DeltaG[i].push_back (1.813/1000*std::pow(T,2)-1.05*T+150.289);
747 : }
748 0 : if (AtomResidueName[1][i].compare("THR")==0) {
749 0 : DeltaG[i].push_back (-2.64/1000*std::pow(T,2)+1.591*T-240.267);
750 : }
751 0 : if (AtomResidueName[1][i].compare("TRP")==0) {
752 0 : DeltaG[i].push_back (-11.66/1000*std::pow(T,2)+6.916*T-1024.284);
753 : }
754 0 : if (AtomResidueName[1][i].compare("TYR")==0) {
755 0 : DeltaG[i].push_back (-7.513/1000*std::pow(T,2)+4.478*T-666.484);
756 : }
757 0 : if (AtomResidueName[1][i].compare("VAL")==0) {
758 0 : DeltaG[i].push_back (-4.902/1000*std::pow(T,2)+2.921*T-433.013);
759 : }
760 0 : DeltaG[natoms].push_back (-0.6927/1000*std::pow(T,2)+0.4404*T-68.724);
761 : }
762 : }
763 : }
764 :
765 0 : Ti = T;
766 :
767 0 : if (firstStepFlag ==0) {
768 0 : if (approach!=2 && approach!=3) {
769 0 : cout << "Unknown approach " << approach << endl;
770 : }
771 0 : for(unsigned i=0; i<natoms; i++) {
772 0 : if (DeltaG[i].size()==0 ) {
773 0 : cout << "Delta G value for residue " << AtomResidueName[1][i] << " not found " << endl;
774 0 : error ("missing Delta G parameter\n");
775 : }
776 : }
777 : }
778 0 : }
779 :
780 :
781 : //calculates neighbor list
782 24 : void SASA_HASEL::calcNlist() {
783 24 : if(!nopbc) {
784 24 : makeWhole();
785 : }
786 :
787 2904 : for(unsigned i = 0; i < natoms; i++) {
788 2880 : Nlist[i].clear();
789 : }
790 :
791 2904 : for(unsigned i = 0; i < natoms; i++) {
792 2880 : if (SASAparam[i].size()>0) {
793 103680 : for (unsigned j = 0; j < i; j++) {
794 101952 : if (SASAparam[j].size()>0) {
795 61344 : const Vector Delta_ij_vec = delta( getPosition(i), getPosition(j) );
796 61344 : double Delta_ij_mod = Delta_ij_vec.modulo()*10;
797 61344 : double overlapD = SASAparam[i][0]+SASAparam[j][0];
798 61344 : if (Delta_ij_mod < overlapD) {
799 26748 : Nlist.at(i).push_back (j);
800 26748 : Nlist.at(j).push_back (i);
801 : }
802 : }
803 : }
804 : }
805 : }
806 24 : }
807 :
808 :
809 : //calculates SASA according to Hasel et al., Tetrahedron Computer Methodology Vol. 1, No. 2, pp. 103-116, 1988
810 24 : void SASA_HASEL::calculate() {
811 24 : if(!nopbc) {
812 24 : makeWhole();
813 : }
814 :
815 24 : if(getExchangeStep()) {
816 0 : nl_update = 0;
817 : }
818 24 : if (firstStepFlag ==0) {
819 2 : readPDB();
820 2 : readSASAparam();
821 : }
822 24 : if (nl_update == 0) {
823 24 : calcNlist();
824 : }
825 :
826 :
827 24 : auto* moldat = plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
828 24 : if( ! moldat ) {
829 0 : error("Unable to find MOLINFO in input");
830 : }
831 : double Si, sasai, bij;
832 24 : double sasa = 0;
833 24 : vector<Vector> derivatives( natoms );
834 2904 : for(unsigned i = 0; i < natoms; i++) {
835 2880 : derivatives[i][0] = 0.;
836 2880 : derivatives[i][1] = 0.;
837 2880 : derivatives[i][2] = 0.;
838 : }
839 :
840 24 : Tensor virial;
841 24 : vector <double> ddij_di(3);
842 24 : vector <double> dbij_di(3);
843 24 : vector <double> dAijt_di(3);
844 :
845 24 : if( sasa_type==TOTAL ) {
846 1452 : for(unsigned i = 0; i < natoms; i++) {
847 1440 : if(SASAparam[i].size() > 0) {
848 864 : double ri = SASAparam[i][0];
849 864 : Si = 4*M_PI*ri*ri;
850 : sasai = 1.0;
851 :
852 864 : vector <vector <double> > derTerm( Nlist[i].size(), vector <double>(3));
853 :
854 864 : dAijt_di[0] = 0;
855 864 : dAijt_di[1] = 0;
856 864 : dAijt_di[2] = 0;
857 864 : int NumRes_i = moldat->getResidueNumber(atoms[i]);
858 :
859 27612 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
860 : double pij = 0.3516;
861 :
862 26748 : int NumRes_j = moldat->getResidueNumber(atoms[Nlist[i][j]]);
863 26748 : if (NumRes_i==NumRes_j) {
864 4320 : if (CONNECTparam[i][0].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][1].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][2].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][3].compare(AtomResidueName[0][Nlist[i][j]])==0) {
865 : pij = 0.8875;
866 : }
867 : }
868 26748 : if ( abs(NumRes_i-NumRes_j) == 1 ) {
869 8868 : if ((AtomResidueName[0][i] == "N" && AtomResidueName[0][Nlist[i][j]]== "CA") || (AtomResidueName[0][Nlist[i][j]] == "N" && AtomResidueName[0][i]== "CA")) {
870 : pij = 0.8875;
871 : }
872 : }
873 :
874 26748 : const Vector d_ij_vec = delta( getPosition(i), getPosition(Nlist[i][j]) );
875 26748 : double d_ij = d_ij_vec.modulo()*10;
876 :
877 26748 : double rj = SASAparam[Nlist[i][j]][0];
878 26748 : bij = M_PI*ri*(ri+rj-d_ij)*(1+(rj-ri)/d_ij); //Angstrom2
879 :
880 26748 : sasai = sasai*(1-SASAparam[i][1]*pij*bij/Si); //nondimensional
881 :
882 26748 : ddij_di[0] = -10*(getPosition(Nlist[i][j])[0]-getPosition(i)[0])/d_ij; //nondimensional
883 26748 : ddij_di[1] = -10*(getPosition(Nlist[i][j])[1]-getPosition(i)[1])/d_ij;
884 26748 : ddij_di[2] = -10*(getPosition(Nlist[i][j])[2]-getPosition(i)[2])/d_ij;
885 :
886 26748 : dbij_di[0] = -M_PI*ri*ddij_di[0]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij)); //Angstrom
887 26748 : dbij_di[1] = -M_PI*ri*ddij_di[1]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij));
888 26748 : dbij_di[2] = -M_PI*ri*ddij_di[2]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij));
889 :
890 26748 : dAijt_di[0] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[0]; //Angstrom-1
891 26748 : dAijt_di[1] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[1];
892 26748 : dAijt_di[2] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[2];
893 :
894 26748 : derTerm[j][0] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[0]; //Angstrom-1
895 26748 : derTerm[j][1] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[1];
896 26748 : derTerm[j][2] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[2];
897 :
898 : }
899 :
900 864 : sasa += Si*sasai/100; //nm2
901 :
902 864 : derivatives[i][0] += Si*sasai/10*dAijt_di[0]; //nm
903 864 : derivatives[i][1] += Si*sasai/10*dAijt_di[1];
904 864 : derivatives[i][2] += Si*sasai/10*dAijt_di[2];
905 :
906 27612 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
907 26748 : derivatives[Nlist[i][j]][0] += Si*sasai/10*derTerm[j][0]; //nm
908 26748 : derivatives[Nlist[i][j]][1] += Si*sasai/10*derTerm[j][1];
909 26748 : derivatives[Nlist[i][j]][2] += Si*sasai/10*derTerm[j][2];
910 : }
911 864 : }
912 : }
913 : }
914 :
915 :
916 24 : if( sasa_type==TRANSFER ) {
917 :
918 12 : if (firstStepFlag ==0) {
919 1 : readMaxSurf();
920 : }
921 :
922 12 : if (firstStepFlag ==0 && DeltaGValues != "absent") {
923 1 : readDeltaG();
924 : }
925 :
926 12 : if (DeltaGValues == "absent") {
927 0 : computeDeltaG();
928 : }
929 :
930 :
931 1452 : for(unsigned i = 0; i < natoms; i++) {
932 :
933 :
934 :
935 1440 : if(SASAparam[i].size() > 0) {
936 864 : double ri = SASAparam[i][0];
937 864 : Si = 4*M_PI*ri*ri;
938 : sasai = 1.0;
939 :
940 864 : vector <vector <double> > derTerm( Nlist[i].size(), vector <double>(3));
941 :
942 864 : dAijt_di[0] = 0;
943 864 : dAijt_di[1] = 0;
944 864 : dAijt_di[2] = 0;
945 864 : int NumRes_i = moldat->getResidueNumber(atoms[i]);
946 :
947 27612 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
948 : double pij = 0.3516;
949 :
950 26748 : int NumRes_j = moldat->getResidueNumber(atoms[Nlist[i][j]]);
951 26748 : if (NumRes_i==NumRes_j) {
952 4320 : if (CONNECTparam[i][0].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][1].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][2].compare(AtomResidueName[0][Nlist[i][j]])==0 || CONNECTparam[i][3].compare(AtomResidueName[0][Nlist[i][j]])==0) {
953 : pij = 0.8875;
954 : }
955 : }
956 26748 : if ( abs(NumRes_i-NumRes_j) == 1 ) {
957 8868 : if ((AtomResidueName[0][i] == "N" && AtomResidueName[0][Nlist[i][j]]== "CA") || (AtomResidueName[0][Nlist[i][j]] == "N" && AtomResidueName[0][i]== "CA")) {
958 : pij = 0.8875;
959 : }
960 : }
961 :
962 26748 : const Vector d_ij_vec = delta( getPosition(i), getPosition(Nlist[i][j]) );
963 26748 : double d_ij = d_ij_vec.modulo()*10;
964 :
965 26748 : double rj = SASAparam[Nlist[i][j]][0];
966 26748 : bij = M_PI*ri*(ri+rj-d_ij)*(1+(rj-ri)/d_ij); //Angstrom2
967 :
968 26748 : sasai = sasai*(1-SASAparam[i][1]*pij*bij/Si); //nondimensional
969 :
970 26748 : ddij_di[0] = -10*(getPosition(Nlist[i][j])[0]-getPosition(i)[0])/d_ij; //nondimensional
971 26748 : ddij_di[1] = -10*(getPosition(Nlist[i][j])[1]-getPosition(i)[1])/d_ij;
972 26748 : ddij_di[2] = -10*(getPosition(Nlist[i][j])[2]-getPosition(i)[2])/d_ij;
973 :
974 26748 : dbij_di[0] = -M_PI*ri*ddij_di[0]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij)); //Angstrom
975 26748 : dbij_di[1] = -M_PI*ri*ddij_di[1]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij));
976 26748 : dbij_di[2] = -M_PI*ri*ddij_di[2]*(1+(ri+rj)*(rj-ri)/(d_ij*d_ij));
977 :
978 26748 : dAijt_di[0] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[0]; //Angstrom-1
979 26748 : dAijt_di[1] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[1];
980 26748 : dAijt_di[2] += -1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[2];
981 :
982 26748 : derTerm[j][0] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[0]; //Angstrom-1
983 26748 : derTerm[j][1] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[1];
984 26748 : derTerm[j][2] = 1/(Si/(SASAparam[i][1]*pij)-bij)*dbij_di[2];
985 :
986 : }
987 :
988 864 : if (AtomResidueName[0][i] == "N" || AtomResidueName[0][i] == "CA" || AtomResidueName[0][i] == "C" || AtomResidueName[0][i] == "O" || AtomResidueName[0][i] == "H") {
989 :
990 720 : sasa += Si*sasai/MaxSurf[i][0]*DeltaG[natoms][0]; //kJ/mol
991 :
992 :
993 720 : derivatives[i][0] += Si*sasai*dAijt_di[0]/MaxSurf[i][0]*DeltaG[natoms][0]*10; //kJ/mol/nm
994 720 : derivatives[i][1] += Si*sasai*dAijt_di[1]/MaxSurf[i][0]*DeltaG[natoms][0]*10;
995 720 : derivatives[i][2] += Si*sasai*dAijt_di[2]/MaxSurf[i][0]*DeltaG[natoms][0]*10;
996 : }
997 :
998 864 : if (AtomResidueName[0][i] != "N" && AtomResidueName[0][i] != "CA" && AtomResidueName[0][i] != "C" && AtomResidueName[0][i] != "O" && AtomResidueName[0][i] != "H") {
999 144 : sasa += Si*sasai/MaxSurf[i][1]*DeltaG[i][0]; //kJ/mol
1000 :
1001 144 : derivatives[i][0] += Si*sasai*dAijt_di[0]/MaxSurf[i][1]*DeltaG[i][0]*10; //kJ/mol/nm
1002 144 : derivatives[i][1] += Si*sasai*dAijt_di[1]/MaxSurf[i][1]*DeltaG[i][0]*10;
1003 144 : derivatives[i][2] += Si*sasai*dAijt_di[2]/MaxSurf[i][1]*DeltaG[i][0]*10;
1004 : }
1005 :
1006 :
1007 27612 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
1008 26748 : if (AtomResidueName[0][i] == "N" || AtomResidueName[0][i] == "CA" || AtomResidueName[0][i] == "C" || AtomResidueName[0][i] == "O" || AtomResidueName[0][i] == "H") {
1009 22438 : derivatives[Nlist[i][j]][0] += Si*sasai*10*derTerm[j][0]/MaxSurf[i][0]*DeltaG[natoms][0]; //kJ/mol/nm
1010 22438 : derivatives[Nlist[i][j]][1] += Si*sasai*10*derTerm[j][1]/MaxSurf[i][0]*DeltaG[natoms][0];
1011 22438 : derivatives[Nlist[i][j]][2] += Si*sasai*10*derTerm[j][2]/MaxSurf[i][0]*DeltaG[natoms][0];
1012 : }
1013 :
1014 26748 : if (AtomResidueName[0][i] != "N" && AtomResidueName[0][i] != "CA" && AtomResidueName[0][i] != "C" && AtomResidueName[0][i] != "O" && AtomResidueName[0][i] != "H") {
1015 4310 : derivatives[Nlist[i][j]][0] += Si*sasai*10*derTerm[j][0]/MaxSurf[i][1]*DeltaG[i][0]; //kJ/mol/nm
1016 4310 : derivatives[Nlist[i][j]][1] += Si*sasai*10*derTerm[j][1]/MaxSurf[i][1]*DeltaG[i][0];
1017 4310 : derivatives[Nlist[i][j]][2] += Si*sasai*10*derTerm[j][2]/MaxSurf[i][1]*DeltaG[i][0];
1018 : }
1019 : }
1020 864 : }
1021 : }
1022 : }
1023 :
1024 :
1025 2904 : for(unsigned i=0; i<natoms; i++) {
1026 2880 : setAtomsDerivatives(i,derivatives[i]);
1027 5760 : virial -= Tensor(getPosition(i),derivatives[i]);
1028 : }
1029 :
1030 24 : setBoxDerivatives(virial);
1031 24 : setValue(sasa);
1032 24 : firstStepFlag = 1;
1033 24 : ++nl_update;
1034 24 : if (nl_update == stride) {
1035 24 : nl_update = 0;
1036 : }
1037 48 : }
1038 :
1039 : }//namespace sasa
1040 : }//namespace PLMD
|
