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_LCPO
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 LCPO algorithm is used for the calculation (please, read and cite the first paper cited below). 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 the 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 : \verbatim
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 : \endverbatim
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)
85 : is freely available [here](https://github.com/andrea-arsiccio/DeltaG-calculation). The script automatically outputs a DeltaG.dat file compatible with this SASA module.
86 :
87 : 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.
88 :
89 : 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.
90 :
91 : In case you want to recover the old behavior you should use the NOPBC flag.
92 : In that case you need to take care that atoms are in the correct periodic image.
93 :
94 : The SASA may also be computed using the SASA_HASEL collective variable, which makes use of the algorithm described in the third paper cited below. SASA_HASEL is less accurate then SASA_LCPO, but the computation is faster.
95 :
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_LCPO 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_LCPO 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 blow, and take into account the effect of temperature using approach 2 as described in the paper.
118 :
119 : ```plumed
120 : SASA_LCPO 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_LCPO : 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 > > LCPOparam;
146 : unsigned natoms;
147 : vector < vector < double > > MaxSurf;
148 : vector < vector < double > > DeltaG;
149 : vector < vector < int > > Nlist;
150 : public:
151 : static void registerKeywords(Keywords& keys);
152 : explicit SASA_LCPO(const ActionOptions&);
153 : void readPDB();
154 : map<string, vector<double> > setupLCPOparam();
155 : void readLCPOparam();
156 : void calcNlist();
157 : map<string, vector<double> > setupMaxSurfMap();
158 : void readMaxSurf();
159 : void readDeltaG();
160 : void computeDeltaG();
161 : virtual void calculate();
162 : };
163 :
164 : PLUMED_REGISTER_ACTION(SASA_LCPO,"SASA_LCPO")
165 :
166 6 : void SASA_LCPO::registerKeywords(Keywords& keys) {
167 6 : Colvar::registerKeywords(keys);
168 6 : keys.add("atoms","ATOMS","the group of atoms that you are calculating the SASA for");
169 6 : keys.add("compulsory","TYPE","TOTAL","The type of calculation you want to perform. Can be TOTAL or TRANSFER");
170 6 : keys.add("compulsory", "NL_STRIDE", "The frequency with which the neighbor list is updated.");
171 6 : keys.add("optional","DELTAGFILE","a file containing the free energy values for backbone and sidechains. 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");
172 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");
173 12 : keys.setValueDescription("scalar","the solvent accessible surface area (SASA) of the molecule");
174 6 : keys.addDOI("10.1002/(SICI)1096-987X(19990130)20:2<217::AID-JCC4>3.0.CO;2-A");
175 6 : keys.addDOI("10.1021/acs.jpcb.1c01694");
176 6 : keys.addDOI("10.1021/acs.jpcb.2c00889");
177 6 : keys.addDOI("10.1021/acs.jpcb.1c04398");
178 6 : keys.addDOI("10.1016/0898-5529(88)90015-2");
179 6 : }
180 :
181 :
182 2 : SASA_LCPO::SASA_LCPO(const ActionOptions&ao):
183 : PLUMED_COLVAR_INIT(ao),
184 2 : nopbc(false),
185 4 : DeltaGValues("absent"),
186 2 : stride(10),
187 2 : nl_update(0),
188 2 : firstStepFlag(0),
189 2 : Ti(0) {
190 2 : rs = 0.14;
191 2 : parse("DELTAGFILE",DeltaGValues);
192 2 : parse("APPROACH", approach);
193 4 : parseAtomList("ATOMS",atoms);
194 2 : if(atoms.size()==0) {
195 0 : error("no atoms specified");
196 : }
197 : std::string Type;
198 2 : parse("TYPE",Type);
199 2 : parse("NL_STRIDE", stride);
200 2 : parseFlag("NOPBC",nopbc);
201 2 : checkRead();
202 :
203 2 : if(Type=="TOTAL") {
204 1 : sasa_type=TOTAL;
205 1 : } else if(Type=="TRANSFER") {
206 1 : sasa_type=TRANSFER;
207 : } else {
208 0 : error("Unknown SASA type");
209 : }
210 :
211 2 : switch(sasa_type) {
212 1 : case TOTAL:
213 1 : log.printf(" TOTAL SASA;");
214 : break;
215 1 : case TRANSFER:
216 1 : log.printf(" TRANSFER MODEL;");
217 : break;
218 : }
219 :
220 2 : log.printf(" atoms involved : ");
221 242 : for(unsigned i=0; i<atoms.size(); ++i) {
222 240 : if(i%25==0) {
223 10 : log<<"\n";
224 : }
225 240 : log.printf("%d ",atoms[i].serial());
226 : }
227 2 : log.printf("\n");
228 :
229 2 : if(nopbc) {
230 0 : log<<" PBC will be ignored\n";
231 : } else {
232 2 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
233 : }
234 :
235 :
236 2 : addValueWithDerivatives();
237 2 : setNotPeriodic();
238 2 : requestAtoms(atoms);
239 :
240 2 : natoms = getNumberOfAtoms();
241 2 : AtomResidueName.resize(2);
242 2 : LCPOparam.resize(natoms);
243 2 : MaxSurf.resize(natoms);
244 2 : DeltaG.resize(natoms+1);
245 2 : Nlist.resize(natoms);
246 :
247 :
248 2 : }
249 :
250 :
251 : //splits strings into tokens. Used to read into LCPO parameters file and into reference pdb file
252 : template <class Container>
253 0 : void split(const std::string& str, Container& cont) {
254 0 : std::istringstream iss(str);
255 0 : std::copy(std::istream_iterator<std::string>(iss),
256 0 : std::istream_iterator<std::string>(),
257 : std::back_inserter(cont));
258 0 : }
259 :
260 :
261 : //reads input PDB file provided with MOLINFO, and assigns atom and residue names to each atom involved in the CV
262 :
263 2 : void SASA_LCPO::readPDB() {
264 2 : auto* moldat = plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
265 2 : if( ! moldat ) {
266 0 : error("Unable to find MOLINFO in input");
267 : }
268 2 : AtomResidueName[0].clear();
269 2 : AtomResidueName[1].clear();
270 :
271 242 : for(unsigned i=0; i<natoms; i++) {
272 240 : string Aname = moldat->getAtomName(atoms[i]);
273 240 : string Rname = moldat->getResidueName(atoms[i]);
274 240 : AtomResidueName[0].push_back (Aname);
275 240 : AtomResidueName[1].push_back (Rname);
276 : }
277 :
278 2 : }
279 :
280 : //LCPO parameters database
281 2 : map<string, vector<double> > SASA_LCPO::setupLCPOparam() {
282 : map<string, vector<double> > lcpomap;
283 2 : lcpomap["ALA_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
284 2 : lcpomap["ALA_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
285 2 : lcpomap["ALA_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
286 2 : lcpomap["ALA_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
287 2 : lcpomap["ALA_CB"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
288 2 : lcpomap["ASP_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
289 2 : lcpomap["ASP_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
290 2 : lcpomap["ASP_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
291 2 : lcpomap["ASP_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
292 2 : lcpomap["ASP_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
293 2 : lcpomap["ASP_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
294 2 : lcpomap["ASP_OD1"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
295 2 : lcpomap["ASP_OD2"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
296 2 : lcpomap["ASN_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
297 2 : lcpomap["ASN_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
298 2 : lcpomap["ASN_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
299 2 : lcpomap["ASN_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
300 2 : lcpomap["ASN_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
301 2 : lcpomap["ASN_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
302 2 : lcpomap["ASN_OD1"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
303 2 : lcpomap["ASN_ND2"] = { 1.65, 0.73511, -0.22116, -0.00089148, 0.0002523};
304 2 : lcpomap["ARG_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
305 2 : lcpomap["ARG_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
306 2 : lcpomap["ARG_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
307 2 : lcpomap["ARG_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
308 2 : lcpomap["ARG_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
309 2 : lcpomap["ARG_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
310 2 : lcpomap["ARG_CD"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
311 2 : lcpomap["ARG_NE"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
312 2 : lcpomap["ARG_NH1"] = { 1.65, 0.73511, -0.22116, -0.00089148, 0.0002523};
313 2 : lcpomap["ARG_NH2"] = { 1.65, 0.73511, -0.22116, -0.00089148, 0.0002523};
314 2 : lcpomap["ARG_CZ"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
315 2 : lcpomap["CYS_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
316 2 : lcpomap["CYS_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
317 2 : lcpomap["CYS_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
318 2 : lcpomap["CYS_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
319 2 : lcpomap["CYS_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
320 2 : lcpomap["CYS_SG"] = { 1.9, 0.54581, -0.19477, -0.0012873, 0.00029247};
321 2 : lcpomap["GLU_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
322 2 : lcpomap["GLU_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
323 2 : lcpomap["GLU_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
324 2 : lcpomap["GLU_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
325 2 : lcpomap["GLU_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
326 2 : lcpomap["GLU_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
327 2 : lcpomap["GLU_CD"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
328 2 : lcpomap["GLU_OE1"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
329 2 : lcpomap["GLU_OE2"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
330 2 : lcpomap["GLN_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
331 2 : lcpomap["GLN_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
332 2 : lcpomap["GLN_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
333 2 : lcpomap["GLN_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
334 2 : lcpomap["GLN_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
335 2 : lcpomap["GLN_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
336 2 : lcpomap["GLN_CD"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
337 2 : lcpomap["GLN_OE1"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
338 2 : lcpomap["GLN_NE2"] = { 1.65, 0.73511, -0.22116, -0.00089148, 0.0002523};
339 2 : lcpomap["GLY_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
340 2 : lcpomap["GLY_CA"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
341 2 : lcpomap["GLY_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
342 2 : lcpomap["GLY_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
343 2 : lcpomap["HIS_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
344 2 : lcpomap["HIS_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
345 2 : lcpomap["HIS_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
346 2 : lcpomap["HIS_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
347 2 : lcpomap["HIS_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
348 2 : lcpomap["HIS_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
349 2 : lcpomap["HIS_ND1"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
350 2 : lcpomap["HIS_CE1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
351 2 : lcpomap["HIS_NE2"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
352 2 : lcpomap["HIS_CD2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
353 2 : lcpomap["ILE_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
354 2 : lcpomap["ILE_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
355 2 : lcpomap["ILE_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
356 2 : lcpomap["ILE_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
357 2 : lcpomap["ILE_CB"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
358 2 : lcpomap["ILE_CG2"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
359 2 : lcpomap["ILE_CG1"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
360 2 : lcpomap["ILE_CD1"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
361 2 : lcpomap["LEU_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
362 2 : lcpomap["LEU_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
363 2 : lcpomap["LEU_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
364 2 : lcpomap["LEU_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
365 2 : lcpomap["LEU_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
366 2 : lcpomap["LEU_CG"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
367 2 : lcpomap["LEU_CD1"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
368 2 : lcpomap["LEU_CD2"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
369 2 : lcpomap["LYS_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
370 2 : lcpomap["LYS_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
371 2 : lcpomap["LYS_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
372 2 : lcpomap["LYS_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
373 2 : lcpomap["LYS_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
374 2 : lcpomap["LYS_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
375 2 : lcpomap["LYS_CD"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
376 2 : lcpomap["LYS_CE"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
377 2 : lcpomap["LYS_NZ"] = { 1.65, 0.73511, -0.22116, -0.00089148, 0.0002523};
378 2 : lcpomap["MET_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
379 2 : lcpomap["MET_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
380 2 : lcpomap["MET_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
381 2 : lcpomap["MET_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
382 2 : lcpomap["MET_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
383 2 : lcpomap["MET_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
384 2 : lcpomap["MET_SD"] = { 1.9, 0.54581, -0.19477, -0.0012873, 0.00029247};
385 2 : lcpomap["MET_CE"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
386 2 : lcpomap["PHE_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
387 2 : lcpomap["PHE_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
388 2 : lcpomap["PHE_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
389 2 : lcpomap["PHE_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
390 2 : lcpomap["PHE_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
391 2 : lcpomap["PHE_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
392 2 : lcpomap["PHE_CD1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
393 2 : lcpomap["PHE_CE1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
394 2 : lcpomap["PHE_CZ"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
395 2 : lcpomap["PHE_CE2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
396 2 : lcpomap["PHE_CD2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
397 2 : lcpomap["PRO_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
398 2 : lcpomap["PRO_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
399 2 : lcpomap["PRO_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
400 2 : lcpomap["PRO_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
401 2 : lcpomap["PRO_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
402 2 : lcpomap["PRO_CG"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
403 2 : lcpomap["PRO_CD"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
404 2 : lcpomap["SER_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
405 2 : lcpomap["SER_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
406 2 : lcpomap["SER_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
407 2 : lcpomap["SER_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
408 2 : lcpomap["SER_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
409 2 : lcpomap["SER_OG"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
410 2 : lcpomap["THR_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
411 2 : lcpomap["THR_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
412 2 : lcpomap["THR_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
413 2 : lcpomap["THR_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
414 2 : lcpomap["THR_CB"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
415 2 : lcpomap["THR_CG2"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
416 2 : lcpomap["THR_OG1"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
417 2 : lcpomap["TRP_N"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
418 2 : lcpomap["TRP_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
419 2 : lcpomap["TRP_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
420 2 : lcpomap["TRP_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
421 2 : lcpomap["TRP_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
422 2 : lcpomap["TRP_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
423 2 : lcpomap["TRP_CD1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
424 2 : lcpomap["TRP_NE1"] = { 1.65, 0.41102, -0.12254, -7.5448e-05, 0.00011804};
425 2 : lcpomap["TRP_CE2"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
426 2 : lcpomap["TRP_CZ2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
427 2 : lcpomap["TRP_CH2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
428 2 : lcpomap["TRP_CZ3"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
429 2 : lcpomap["TRP_CE3"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
430 2 : lcpomap["TRP_CD2"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
431 2 : lcpomap["TYR_N"] = { 1.65, 0.062577, -0.017874, -8.312e-05, 1.9849e-05};
432 2 : lcpomap["TYR_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
433 2 : lcpomap["TYR_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
434 2 : lcpomap["TYR_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
435 2 : lcpomap["TYR_CB"] = { 1.7, 0.56482, -0.19608, -0.0010219, 0.0002658};
436 2 : lcpomap["TYR_CG"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
437 2 : lcpomap["TYR_CD1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
438 2 : lcpomap["TYR_CE1"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
439 2 : lcpomap["TYR_CZ"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
440 2 : lcpomap["TYR_OH"] = { 1.6, 0.77914, -0.25262, -0.0016056, 0.00035071};
441 2 : lcpomap["TYR_CE2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
442 2 : lcpomap["TYR_CD2"] = { 1.7, 0.51245, -0.15966, -0.00019781, 0.00016392};
443 2 : lcpomap["VAL_N"] = { 1.65, 0.062577, -0.017874, -8.312e-05, 1.9849e-05};
444 2 : lcpomap["VAL_CA"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
445 2 : lcpomap["VAL_C"] = { 1.7, 0.070344, -0.019015, -2.2009e-05, 1.6875e-05};
446 2 : lcpomap["VAL_O"] = { 1.6, 0.68563, -0.1868, -0.00135573, 0.00023743};
447 2 : lcpomap["VAL_CB"] = { 1.7, 0.23348, -0.072627, -0.00020079, 7.967e-05};
448 2 : lcpomap["VAL_CG1"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
449 2 : lcpomap["VAL_CG2"] = { 1.7, 0.77887, -0.28063, -0.0012968, 0.00039328};
450 2 : return lcpomap;
451 : }
452 :
453 : //assigns LCPO parameters to each atom reading from database
454 2 : void SASA_LCPO::readLCPOparam() {
455 :
456 242 : for(unsigned i=0; i<natoms; i++) {
457 240 : LCPOparam[i].clear();
458 : }
459 :
460 : map<string, vector<double> > lcpomap;
461 2 : lcpomap = setupLCPOparam();
462 2 : vector<double> LCPOparamVector;
463 : string identifier;
464 242 : for(unsigned i=0; i<natoms; i++) {
465 240 : if ((AtomResidueName[0][i][0]=='O') || (AtomResidueName[0][i][0]=='N') || (AtomResidueName[0][i][0]=='C') || (AtomResidueName[0][i][0]=='S')) {
466 120 : identifier = AtomResidueName[1][i]+"_"+AtomResidueName[0][i];
467 120 : LCPOparamVector = lcpomap.at(identifier);
468 120 : LCPOparam[i].push_back (LCPOparamVector[0]+rs*10);
469 120 : LCPOparam[i].push_back (LCPOparamVector[1]);
470 120 : LCPOparam[i].push_back (LCPOparamVector[2]);
471 120 : LCPOparam[i].push_back (LCPOparamVector[3]);
472 120 : LCPOparam[i].push_back (LCPOparamVector[4]);
473 : }
474 : }
475 :
476 :
477 242 : for(unsigned i=0; i<natoms; i++) {
478 240 : if (LCPOparam[i].size()==0 ) {
479 120 : if ((AtomResidueName[0][i][0]=='O') || (AtomResidueName[0][i][0]=='N') || (AtomResidueName[0][i][0]=='C') || (AtomResidueName[0][i][0]=='S')) {
480 0 : cout << "Could not find LCPO parameters for atom " << AtomResidueName[0][i] << " of residue " << AtomResidueName[1][i] << endl;
481 0 : error ("missing LCPO parameters\n");
482 : }
483 : }
484 : }
485 :
486 2 : if (AtomResidueName[0][0] == "N") {
487 2 : LCPOparam[0][1] = 7.3511e-01;
488 2 : LCPOparam[0][2] = -2.2116e-01;
489 2 : LCPOparam[0][3] = -8.9148e-04;
490 2 : LCPOparam[0][4] = 2.5230e-04;
491 : }
492 :
493 2 : if (AtomResidueName[0][natoms-1] == "O") {
494 2 : LCPOparam[natoms-1][1] = 8.8857e-01;
495 2 : LCPOparam[natoms-1][2] = -3.3421e-01;
496 2 : LCPOparam[natoms-1][3] = -1.8683e-03;
497 2 : LCPOparam[natoms-1][4] = 4.9372e-04;
498 : }
499 4 : }
500 :
501 :
502 : //Max Surf values, used only if TYPE=TRANSFER
503 1 : map<string, vector<double> > SASA_LCPO::setupMaxSurfMap() {
504 : // Max Surface Area for backbone and sidechain, in nm2
505 : map<string, vector<double> > valuemap;
506 1 : valuemap ["ALA"]= {0.671446,0.420263,};
507 1 : valuemap ["ARG"]= {0.578582,1.95454,};
508 1 : valuemap ["ASN"]= {0.559411,1.07102,};
509 1 : valuemap ["ASP"]= {0.558363,1.03971,};
510 1 : valuemap ["CYS"]= {0.609271,0.657612,};
511 1 : valuemap ["GLN"]= {0.565651,1.433031,};
512 1 : valuemap ["GLU"]= {0.572399,1.41848,};
513 1 : valuemap ["GLY"]= {0.861439,0,};
514 1 : valuemap ["HIS"]= {0.559929,1.143827,};
515 1 : valuemap ["ILE"]= {0.553491,1.25334,};
516 1 : valuemap ["LEU"]= {0.570103,1.260459,};
517 1 : valuemap ["LYS"]= {0.580304,1.687487,};
518 1 : valuemap ["MET"]= {0.5856,1.498954,};
519 1 : valuemap ["PHE"]= {0.54155,1.394861,};
520 1 : valuemap ["PRO"]= {0.456048,0.849461,};
521 1 : valuemap ["SER"]= {0.59074,0.714538,};
522 1 : valuemap ["THR"]= {0.559116,0.951644,};
523 1 : valuemap ["TRP"]= {0.55536,1.59324,};
524 1 : valuemap ["TYR"]= {0.451171,1.566918,};
525 1 : valuemap ["VAL"]= {0.454809,0.928685,};
526 1 : return valuemap;
527 : }
528 :
529 :
530 :
531 : //reads maximum surface values per residue type and assigns values to each atom, only used if sasa_type = TRANSFER
532 :
533 1 : void SASA_LCPO::readMaxSurf() {
534 : map<string, vector<double> > valuemap;
535 1 : valuemap = setupMaxSurfMap();
536 1 : vector<double> MaxSurfVector;
537 :
538 121 : for(unsigned i=0; i<natoms; i++) {
539 120 : MaxSurf[i].clear();
540 120 : MaxSurfVector = valuemap.at(AtomResidueName[1][i]);
541 120 : MaxSurf[i].push_back (MaxSurfVector[0]*100);
542 120 : MaxSurf[i].push_back (MaxSurfVector[1]*100);
543 : }
544 2 : }
545 :
546 : //reads file with free energy values for sidechains and for the backbone, and assigns values to each atom. Only used if sasa_type = TRANSFER
547 :
548 1 : void SASA_LCPO::readDeltaG() {
549 :
550 121 : for(unsigned i=0; i<natoms; i++) {
551 120 : DeltaG[i].clear();
552 : }
553 :
554 : string DeltaGline;
555 1 : fstream DeltaGFile;
556 1 : DeltaGFile.open(DeltaGValues);
557 1 : if (DeltaGFile) {
558 : int backboneflag = 0;
559 23 : while(getline(DeltaGFile, DeltaGline)) {
560 22 : if(!DeltaGline.empty()) {
561 21 : std::vector<std::string> DeltaGtoken;
562 21 : split(DeltaGline, DeltaGtoken);
563 2541 : for(unsigned i=0; i<natoms; i++) {
564 2520 : if (DeltaGtoken[0].compare(AtomResidueName[1][i])==0 ) {
565 120 : DeltaG[i].push_back (std::atof(DeltaGtoken[1].c_str()));
566 : }
567 : }
568 21 : if (DeltaGtoken[0].compare("BACKBONE")==0 ) {
569 : backboneflag = 1;
570 1 : DeltaG[natoms].push_back (std::atof(DeltaGtoken[1].c_str()));
571 : }
572 21 : }
573 : }
574 1 : if ( backboneflag == 0) {
575 0 : error("Cannot find backbone value in Delta G parameters file\n");
576 : }
577 : } else {
578 0 : error("Cannot open DeltaG file");
579 : }
580 :
581 121 : for(unsigned i=0; i<natoms; i++) {
582 120 : if (DeltaG[i].size()==0 ) {
583 0 : cout << "Delta G value for residue " << AtomResidueName[1][i] << " not found " << endl;
584 0 : error ("missing Delta G parameter\n");
585 : }
586 : }
587 :
588 2 : }
589 :
590 : //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.
591 :
592 0 : void SASA_LCPO::computeDeltaG() {
593 :
594 0 : for(unsigned i=0; i<natoms; i++) {
595 0 : DeltaG[i].clear();
596 : }
597 :
598 : double T;
599 0 : T = getkBT()/getKBoltzmann();
600 :
601 0 : if (T != Ti) {
602 0 : for(unsigned i=0; i<natoms; i++) {
603 0 : if (approach==2) {
604 0 : if (AtomResidueName[1][i].compare("ALA")==0) {
605 0 : DeltaG[i].push_back (-2.995/1000*std::pow(T,2)+1.808*T-272.895);
606 : }
607 0 : if (AtomResidueName[1][i].compare("ARG")==0) {
608 0 : DeltaG[i].push_back (-3.182/1000*std::pow(T,2)+1.894*T-282.032);
609 : }
610 0 : if (AtomResidueName[1][i].compare("ASN")==0) {
611 0 : DeltaG[i].push_back (-1.047/1000*std::pow(T,2)+0.6068*T-87.846);
612 : }
613 0 : if (AtomResidueName[1][i].compare("ASP")==0) {
614 0 : DeltaG[i].push_back (-0.1794/1000*std::pow(T,2)+0.1091*T-16.526);
615 : }
616 0 : if (AtomResidueName[1][i].compare("CYS")==0) {
617 0 : DeltaG[i].push_back (-3.09/1000*std::pow(T,2)+1.835*T-272.26);
618 : }
619 0 : if (AtomResidueName[1][i].compare("GLN")==0) {
620 0 : DeltaG[i].push_back (-2.23/1000*std::pow(T,2)+1.335*T-199.707);
621 : }
622 0 : if (AtomResidueName[1][i].compare("GLU")==0) {
623 0 : DeltaG[i].push_back (-1.511/1000*std::pow(T,2)+0.8904*T-131.168);
624 : }
625 0 : if (AtomResidueName[1][i].compare("GLY")==0) {
626 0 : DeltaG[i].push_back (0);
627 : }
628 0 : if (AtomResidueName[1][i].compare("HIS")==0) {
629 0 : DeltaG[i].push_back (-3.482/1000*std::pow(T,2)+2.084*T-311.694);
630 : }
631 0 : if (AtomResidueName[1][i].compare("ILE")==0) {
632 0 : DeltaG[i].push_back (-6.364/1000*std::pow(T,2)+3.8*T-567.444);
633 : }
634 0 : if (AtomResidueName[1][i].compare("LEU")==0) {
635 0 : DeltaG[i].push_back (-7.466/1000*std::pow(T,2)+4.417*T-653.394);
636 : }
637 0 : if (AtomResidueName[1][i].compare("LYS")==0) {
638 0 : DeltaG[i].push_back (-2.091/1000*std::pow(T,2)+1.2458*T-185.549);
639 : }
640 0 : if (AtomResidueName[1][i].compare("MET")==0) {
641 0 : DeltaG[i].push_back (-3.807/1000*std::pow(T,2)+2.272*T-339.007);
642 : }
643 0 : if (AtomResidueName[1][i].compare("PHE")==0) {
644 0 : DeltaG[i].push_back (-7.828/1000*std::pow(T,2)+4.644*T-688.874);
645 : }
646 0 : if (AtomResidueName[1][i].compare("PRO")==0) {
647 0 : DeltaG[i].push_back (-2.347/1000*std::pow(T,2)+1.411*T-212.059);
648 : }
649 0 : if (AtomResidueName[1][i].compare("SER")==0) {
650 0 : DeltaG[i].push_back (1.813/1000*std::pow(T,2)-1.05*T+151.957);
651 : }
652 0 : if (AtomResidueName[1][i].compare("THR")==0) {
653 0 : DeltaG[i].push_back (-2.64/1000*std::pow(T,2)+1.591*T-239.516);
654 : }
655 0 : if (AtomResidueName[1][i].compare("TRP")==0) {
656 0 : DeltaG[i].push_back (-11.66/1000*std::pow(T,2)+6.916*T-1025.293);
657 : }
658 0 : if (AtomResidueName[1][i].compare("TYR")==0) {
659 0 : DeltaG[i].push_back (-7.513/1000*std::pow(T,2)+4.478*T-667.261);
660 : }
661 0 : if (AtomResidueName[1][i].compare("VAL")==0) {
662 0 : DeltaG[i].push_back (-4.902/1000*std::pow(T,2)+2.921*T-435.309);
663 : }
664 0 : DeltaG[natoms].push_back (-0.6962/1000*std::pow(T,2)+0.4426*T-70.015);
665 : }
666 0 : if (approach==3) {
667 0 : if (AtomResidueName[1][i].compare("ALA")==0) {
668 0 : DeltaG[i].push_back (-2.995/1000*std::pow(T,2)+1.808*T-272.105);
669 : }
670 0 : if (AtomResidueName[1][i].compare("ARG")==0) {
671 0 : DeltaG[i].push_back (-3.182/1000*std::pow(T,2)+1.894*T-284.086);
672 : }
673 0 : if (AtomResidueName[1][i].compare("ASN")==0) {
674 0 : DeltaG[i].push_back (-1.047/1000*std::pow(T,2)+0.6068*T-90.597);
675 : }
676 0 : if (AtomResidueName[1][i].compare("ASP")==0) {
677 0 : DeltaG[i].push_back (-0.1794/1000*std::pow(T,2)+0.1091*T-19.143);
678 : }
679 0 : if (AtomResidueName[1][i].compare("CYS")==0) {
680 0 : DeltaG[i].push_back (-3.09/1000*std::pow(T,2)+1.835*T-268.527);
681 : }
682 0 : if (AtomResidueName[1][i].compare("GLN")==0) {
683 0 : DeltaG[i].push_back (-2.23/1000*std::pow(T,2)+1.335*T-201.559);
684 : }
685 0 : if (AtomResidueName[1][i].compare("GLU")==0) {
686 0 : DeltaG[i].push_back (-1.511/1000*std::pow(T,2)+0.8904*T-133.543);
687 : }
688 0 : if (AtomResidueName[1][i].compare("GLY")==0) {
689 0 : DeltaG[i].push_back (0);
690 : }
691 0 : if (AtomResidueName[1][i].compare("HIS")==0) {
692 0 : DeltaG[i].push_back (-3.482/1000*std::pow(T,2)+2.084*T-315.398);
693 : }
694 0 : if (AtomResidueName[1][i].compare("ILE")==0) {
695 0 : DeltaG[i].push_back (-6.364/1000*std::pow(T,2)+3.8*T-564.825);
696 : }
697 0 : if (AtomResidueName[1][i].compare("LEU")==0) {
698 0 : DeltaG[i].push_back (-7.466/1000*std::pow(T,2)+4.417*T-651.483);
699 : }
700 0 : if (AtomResidueName[1][i].compare("LYS")==0) {
701 0 : DeltaG[i].push_back (-2.091/1000*std::pow(T,2)+1.2458*T-187.485);
702 : }
703 0 : if (AtomResidueName[1][i].compare("MET")==0) {
704 0 : DeltaG[i].push_back (-3.807/1000*std::pow(T,2)+2.272*T-339.242);
705 : }
706 0 : if (AtomResidueName[1][i].compare("PHE")==0) {
707 0 : DeltaG[i].push_back (-7.828/1000*std::pow(T,2)+4.644*T-687.134);
708 : }
709 0 : if (AtomResidueName[1][i].compare("PRO")==0) {
710 0 : DeltaG[i].push_back (-2.347/1000*std::pow(T,2)+1.411*T-214.211);
711 : }
712 0 : if (AtomResidueName[1][i].compare("SER")==0) {
713 0 : DeltaG[i].push_back (1.813/1000*std::pow(T,2)-1.05*T+150.289);
714 : }
715 0 : if (AtomResidueName[1][i].compare("THR")==0) {
716 0 : DeltaG[i].push_back (-2.64/1000*std::pow(T,2)+1.591*T-240.267);
717 : }
718 0 : if (AtomResidueName[1][i].compare("TRP")==0) {
719 0 : DeltaG[i].push_back (-11.66/1000*std::pow(T,2)+6.916*T-1024.284);
720 : }
721 0 : if (AtomResidueName[1][i].compare("TYR")==0) {
722 0 : DeltaG[i].push_back (-7.513/1000*std::pow(T,2)+4.478*T-666.484);
723 : }
724 0 : if (AtomResidueName[1][i].compare("VAL")==0) {
725 0 : DeltaG[i].push_back (-4.902/1000*std::pow(T,2)+2.921*T-433.013);
726 : }
727 0 : DeltaG[natoms].push_back (-0.6927/1000*std::pow(T,2)+0.4404*T-68.724);
728 : }
729 : }
730 : }
731 :
732 0 : Ti = T;
733 :
734 0 : if (firstStepFlag ==0) {
735 0 : if (approach!=2 && approach!=3) {
736 0 : cout << "Unknown approach " << approach << endl;
737 : }
738 0 : for(unsigned i=0; i<natoms; i++) {
739 0 : if (DeltaG[i].size()==0 ) {
740 0 : cout << "Delta G value for residue " << AtomResidueName[1][i] << " not found " << endl;
741 0 : error ("missing Delta G parameter\n");
742 : }
743 : }
744 : }
745 0 : }
746 :
747 :
748 : //calculates neighbor list
749 24 : void SASA_LCPO::calcNlist() {
750 24 : if(!nopbc) {
751 24 : makeWhole();
752 : }
753 :
754 2904 : for(unsigned i = 0; i < natoms; i++) {
755 2880 : Nlist[i].clear();
756 : }
757 :
758 2904 : for(unsigned i = 0; i < natoms; i++) {
759 2880 : if (LCPOparam[i].size()>0) {
760 87264 : for (unsigned j = 0; j < i; j++) {
761 85824 : if (LCPOparam[j].size()>0) {
762 42480 : const Vector Delta_ij_vec = delta( getPosition(i), getPosition(j) );
763 42480 : double Delta_ij_mod = Delta_ij_vec.modulo()*10;
764 42480 : double overlapD = LCPOparam[i][0]+LCPOparam[j][0];
765 42480 : if (Delta_ij_mod < overlapD) {
766 19030 : Nlist.at(i).push_back (j);
767 19030 : Nlist.at(j).push_back (i);
768 : }
769 : }
770 : }
771 : }
772 : }
773 24 : }
774 :
775 :
776 : //calculates SASA according to LCPO algorithm
777 24 : void SASA_LCPO::calculate() {
778 24 : if(!nopbc) {
779 24 : makeWhole();
780 : }
781 :
782 24 : if(getExchangeStep()) {
783 0 : nl_update = 0;
784 : }
785 24 : if (firstStepFlag ==0) {
786 2 : readPDB();
787 2 : readLCPOparam();
788 : }
789 24 : if (nl_update == 0) {
790 24 : calcNlist();
791 : }
792 :
793 :
794 :
795 : double S1, Aij, Ajk, Aijk, Aijt, Ajkt, Aikt;
796 : double dAdd;
797 24 : double sasa = 0;
798 24 : vector<Vector> derivatives( natoms );
799 24 : Tensor virial;
800 24 : vector <double> dAijdc_2t(3);
801 24 : vector <double> dSASA_2_neigh_dc(3);
802 24 : vector <double> ddij_di(3);
803 24 : vector <double> ddik_di(3);
804 :
805 24 : if( sasa_type==TOTAL ) {
806 1452 : for(unsigned i = 0; i < natoms; i++) {
807 1440 : derivatives[i][0] = 0.;
808 1440 : derivatives[i][1] = 0.;
809 1440 : derivatives[i][2] = 0.;
810 1440 : if ( LCPOparam[i].size()>1) {
811 720 : if (LCPOparam[i][1]>0.0) {
812 : Aij = 0.0;
813 : Aijk = 0.0;
814 : Ajk = 0.0;
815 720 : double ri = LCPOparam[i][0];
816 720 : S1 = 4*M_PI*ri*ri;
817 720 : vector <double> dAijdc_2(3, 0);
818 720 : vector <double> dAijdc_4(3, 0);
819 :
820 :
821 19750 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
822 19030 : const Vector d_ij_vec = delta( getPosition(i), getPosition(Nlist[i][j]) );
823 19030 : double d_ij = d_ij_vec.modulo()*10;
824 :
825 19030 : double rj = LCPOparam[Nlist[i][j]][0];
826 19030 : Aijt = (2*M_PI*ri*(ri-d_ij/2-((ri*ri-rj*rj)/(2*d_ij))));
827 19030 : double sji = (2*M_PI*rj*(rj-d_ij/2+((ri*ri-rj*rj)/(2*d_ij))));
828 :
829 19030 : dAdd = M_PI*rj*(-(ri*ri-rj*rj)/(d_ij*d_ij)-1);
830 :
831 19030 : ddij_di[0] = -10*(getPosition(Nlist[i][j])[0]-getPosition(i)[0])/d_ij;
832 19030 : ddij_di[1] = -10*(getPosition(Nlist[i][j])[1]-getPosition(i)[1])/d_ij;
833 19030 : ddij_di[2] = -10*(getPosition(Nlist[i][j])[2]-getPosition(i)[2])/d_ij;
834 :
835 : Ajkt = 0.0;
836 : Aikt = 0.0;
837 :
838 19030 : vector <double> dSASA_3_neigh_dc(3, 0.0);
839 19030 : vector <double> dSASA_4_neigh_dc(3, 0.0);
840 19030 : vector <double> dSASA_3_neigh_dc2(3, 0.0);
841 19030 : vector <double> dSASA_4_neigh_dc2(3, 0.0);
842 :
843 19030 : dSASA_2_neigh_dc[0] = dAdd * ddij_di[0];
844 19030 : dSASA_2_neigh_dc[1] = dAdd * ddij_di[1];
845 19030 : dSASA_2_neigh_dc[2] = dAdd * ddij_di[2];
846 :
847 19030 : dAdd = M_PI*ri*((ri*ri-rj*rj)/(d_ij*d_ij)-1);
848 :
849 :
850 19030 : dAijdc_2t[0] = dAdd * ddij_di[0];
851 19030 : dAijdc_2t[1] = dAdd * ddij_di[1];
852 19030 : dAijdc_2t[2] = dAdd * ddij_di[2];
853 :
854 560038 : for (unsigned k = 0; k < Nlist[Nlist[i][j]].size(); k++) {
855 541008 : if (std::find (Nlist[i].begin(), Nlist[i].end(), Nlist[Nlist[i][j]][k]) != Nlist[i].end()) {
856 370872 : const Vector d_jk_vec = delta( getPosition(Nlist[i][j]), getPosition(Nlist[Nlist[i][j]][k]) );
857 : const Vector d_ik_vec = delta( getPosition(i), getPosition(Nlist[Nlist[i][j]][k]) );
858 :
859 370872 : double d_jk = d_jk_vec.modulo()*10;
860 370872 : double d_ik = d_ik_vec.modulo()*10;
861 :
862 370872 : double rk = LCPOparam[Nlist[Nlist[i][j]][k]][0];
863 370872 : double sjk = (2*M_PI*rj*(rj-d_jk/2-((rj*rj-rk*rk)/(2*d_jk))));
864 370872 : Ajkt += sjk;
865 370872 : Aikt += (2*M_PI*ri*(ri-d_ik/2-((ri*ri-rk*rk)/(2*d_ik))));
866 :
867 370872 : dAdd = M_PI*ri*((ri*ri-rk*rk)/(d_ik*d_ik)-1);
868 :
869 370872 : ddik_di[0] = -10*(getPosition(Nlist[Nlist[i][j]][k])[0]-getPosition(i)[0])/d_ik;
870 370872 : ddik_di[1] = -10*(getPosition(Nlist[Nlist[i][j]][k])[1]-getPosition(i)[1])/d_ik;
871 370872 : ddik_di[2] = -10*(getPosition(Nlist[Nlist[i][j]][k])[2]-getPosition(i)[2])/d_ik;
872 :
873 :
874 370872 : dSASA_3_neigh_dc[0] += dAdd*ddik_di[0];
875 370872 : dSASA_3_neigh_dc[1] += dAdd*ddik_di[1];
876 370872 : dSASA_3_neigh_dc[2] += dAdd*ddik_di[2];
877 :
878 370872 : dAdd = M_PI*rk*(-(ri*ri-rk*rk)/(d_ik*d_ik)-1);
879 :
880 370872 : dSASA_3_neigh_dc2[0] += dAdd*ddik_di[0];
881 370872 : dSASA_3_neigh_dc2[1] += dAdd*ddik_di[1];
882 370872 : dSASA_3_neigh_dc2[2] += dAdd*ddik_di[2];
883 :
884 370872 : dSASA_4_neigh_dc2[0] += sjk*dAdd*ddik_di[0];
885 370872 : dSASA_4_neigh_dc2[1] += sjk*dAdd*ddik_di[1];
886 370872 : dSASA_4_neigh_dc2[2] += sjk*dAdd*ddik_di[2];
887 :
888 : }
889 : }
890 19030 : dSASA_4_neigh_dc[0] = sji*dSASA_3_neigh_dc[0] + dSASA_4_neigh_dc2[0];
891 19030 : dSASA_4_neigh_dc[1] = sji*dSASA_3_neigh_dc[1] + dSASA_4_neigh_dc2[1];
892 19030 : dSASA_4_neigh_dc[2] = sji*dSASA_3_neigh_dc[2] + dSASA_4_neigh_dc2[2];
893 :
894 19030 : dSASA_3_neigh_dc[0] += dSASA_3_neigh_dc2[0];
895 19030 : dSASA_3_neigh_dc[1] += dSASA_3_neigh_dc2[1];
896 19030 : dSASA_3_neigh_dc[2] += dSASA_3_neigh_dc2[2];
897 :
898 19030 : dSASA_4_neigh_dc[0] += dSASA_2_neigh_dc[0] * Aikt;
899 19030 : dSASA_4_neigh_dc[1] += dSASA_2_neigh_dc[1] * Aikt;
900 19030 : dSASA_4_neigh_dc[2] += dSASA_2_neigh_dc[2] * Aikt;
901 :
902 :
903 19030 : derivatives[i][0] += (dSASA_2_neigh_dc[0]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[0]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[0]*LCPOparam[Nlist[i][j]][4])/10;
904 19030 : derivatives[i][1] += (dSASA_2_neigh_dc[1]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[1]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[1]*LCPOparam[Nlist[i][j]][4])/10;
905 19030 : derivatives[i][2] += (dSASA_2_neigh_dc[2]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[2]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[2]*LCPOparam[Nlist[i][j]][4])/10;
906 :
907 :
908 19030 : Aijk += (Aijt * Ajkt);
909 19030 : Aij += Aijt;
910 19030 : Ajk += Ajkt;
911 :
912 19030 : dAijdc_2[0] += dAijdc_2t[0];
913 19030 : dAijdc_2[1] += dAijdc_2t[1];
914 19030 : dAijdc_2[2] += dAijdc_2t[2];
915 :
916 :
917 19030 : dAijdc_4[0] += Ajkt*dAijdc_2t[0];
918 19030 : dAijdc_4[1] += Ajkt*dAijdc_2t[1];
919 19030 : dAijdc_4[2] += Ajkt*dAijdc_2t[2];
920 :
921 :
922 19030 : }
923 720 : double sasai = (LCPOparam[i][1]*S1+LCPOparam[i][2]*Aij+LCPOparam[i][3]*Ajk+LCPOparam[i][4]*Aijk);
924 720 : if (sasai > 0 ) {
925 715 : sasa += sasai/100;
926 : }
927 720 : derivatives[i][0] += (dAijdc_2[0]*LCPOparam[i][2]+dAijdc_4[0]*LCPOparam[i][4])/10;
928 720 : derivatives[i][1] += (dAijdc_2[1]*LCPOparam[i][2]+dAijdc_4[1]*LCPOparam[i][4])/10;
929 720 : derivatives[i][2] += (dAijdc_2[2]*LCPOparam[i][2]+dAijdc_4[2]*LCPOparam[i][4])/10;
930 720 : }
931 : }
932 2880 : virial -= Tensor(getPosition(i),derivatives[i]);
933 : }
934 : }
935 :
936 :
937 :
938 24 : if( sasa_type==TRANSFER ) {
939 :
940 12 : if (firstStepFlag ==0) {
941 1 : readMaxSurf();
942 : }
943 :
944 12 : if (firstStepFlag ==0 && DeltaGValues != "absent") {
945 1 : readDeltaG();
946 : }
947 :
948 12 : if (DeltaGValues == "absent") {
949 0 : computeDeltaG();
950 : }
951 :
952 :
953 1452 : for(unsigned i = 0; i < natoms; i++) {
954 :
955 :
956 :
957 1440 : derivatives[i][0] = 0.;
958 1440 : derivatives[i][1] = 0.;
959 1440 : derivatives[i][2] = 0.;
960 :
961 1440 : if ( LCPOparam[i].size()>1) {
962 720 : if (LCPOparam[i][1]>0.0) {
963 : Aij = 0.0;
964 : Aijk = 0.0;
965 : Ajk = 0.0;
966 720 : double ri = LCPOparam[i][0];
967 720 : S1 = 4*M_PI*ri*ri;
968 720 : vector <double> dAijdc_2(3, 0);
969 720 : vector <double> dAijdc_4(3, 0);
970 :
971 :
972 19750 : for (unsigned j = 0; j < Nlist[i].size(); j++) {
973 19030 : const Vector d_ij_vec = delta( getPosition(i), getPosition(Nlist[i][j]) );
974 19030 : double d_ij = d_ij_vec.modulo()*10;
975 :
976 19030 : double rj = LCPOparam[Nlist[i][j]][0];
977 19030 : Aijt = (2*M_PI*ri*(ri-d_ij/2-((ri*ri-rj*rj)/(2*d_ij))));
978 19030 : double sji = (2*M_PI*rj*(rj-d_ij/2+((ri*ri-rj*rj)/(2*d_ij))));
979 :
980 19030 : dAdd = M_PI*rj*(-(ri*ri-rj*rj)/(d_ij*d_ij)-1);
981 19030 : ddij_di[0] = -10*(getPosition(Nlist[i][j])[0]-getPosition(i)[0])/d_ij;
982 19030 : ddij_di[1] = -10*(getPosition(Nlist[i][j])[1]-getPosition(i)[1])/d_ij;
983 19030 : ddij_di[2] = -10*(getPosition(Nlist[i][j])[2]-getPosition(i)[2])/d_ij;
984 :
985 : Ajkt = 0.0;
986 : Aikt = 0.0;
987 :
988 19030 : vector <double> dSASA_3_neigh_dc(3, 0.0);
989 19030 : vector <double> dSASA_4_neigh_dc(3, 0.0);
990 19030 : vector <double> dSASA_3_neigh_dc2(3, 0.0);
991 19030 : vector <double> dSASA_4_neigh_dc2(3, 0.0);
992 :
993 19030 : dSASA_2_neigh_dc[0] = dAdd * ddij_di[0];
994 19030 : dSASA_2_neigh_dc[1] = dAdd * ddij_di[1];
995 19030 : dSASA_2_neigh_dc[2] = dAdd * ddij_di[2];
996 :
997 19030 : dAdd = M_PI*ri*((ri*ri-rj*rj)/(d_ij*d_ij)-1);
998 :
999 19030 : dAijdc_2t[0] = dAdd * ddij_di[0];
1000 19030 : dAijdc_2t[1] = dAdd * ddij_di[1];
1001 19030 : dAijdc_2t[2] = dAdd * ddij_di[2];
1002 :
1003 560038 : for (unsigned k = 0; k < Nlist[Nlist[i][j]].size(); k++) {
1004 541008 : if (std::find (Nlist[i].begin(), Nlist[i].end(), Nlist[Nlist[i][j]][k]) != Nlist[i].end()) {
1005 370872 : const Vector d_jk_vec = delta( getPosition(Nlist[i][j]), getPosition(Nlist[Nlist[i][j]][k]) );
1006 : const Vector d_ik_vec = delta( getPosition(i), getPosition(Nlist[Nlist[i][j]][k]) );
1007 :
1008 370872 : double d_jk = d_jk_vec.modulo()*10;
1009 370872 : double d_ik = d_ik_vec.modulo()*10;
1010 :
1011 370872 : double rk = LCPOparam[Nlist[Nlist[i][j]][k]][0];
1012 370872 : double sjk = (2*M_PI*rj*(rj-d_jk/2-((rj*rj-rk*rk)/(2*d_jk))));
1013 370872 : Ajkt += sjk;
1014 370872 : Aikt += (2*M_PI*ri*(ri-d_ik/2-((ri*ri-rk*rk)/(2*d_ik))));
1015 :
1016 370872 : dAdd = M_PI*ri*((ri*ri-rk*rk)/(d_ik*d_ik)-1);
1017 :
1018 370872 : ddik_di[0] = -10*(getPosition(Nlist[Nlist[i][j]][k])[0]-getPosition(i)[0])/d_ik;
1019 370872 : ddik_di[1] = -10*(getPosition(Nlist[Nlist[i][j]][k])[1]-getPosition(i)[1])/d_ik;
1020 370872 : ddik_di[2] = -10*(getPosition(Nlist[Nlist[i][j]][k])[2]-getPosition(i)[2])/d_ik;
1021 :
1022 :
1023 370872 : dSASA_3_neigh_dc[0] += dAdd*ddik_di[0];
1024 370872 : dSASA_3_neigh_dc[1] += dAdd*ddik_di[1];
1025 370872 : dSASA_3_neigh_dc[2] += dAdd*ddik_di[2];
1026 :
1027 370872 : dAdd = M_PI*rk*(-(ri*ri-rk*rk)/(d_ik*d_ik)-1);
1028 :
1029 370872 : dSASA_3_neigh_dc2[0] += dAdd*ddik_di[0];
1030 370872 : dSASA_3_neigh_dc2[1] += dAdd*ddik_di[1];
1031 370872 : dSASA_3_neigh_dc2[2] += dAdd*ddik_di[2];
1032 :
1033 370872 : dSASA_4_neigh_dc2[0] += sjk*dAdd*ddik_di[0];
1034 370872 : dSASA_4_neigh_dc2[1] += sjk*dAdd*ddik_di[1];
1035 370872 : dSASA_4_neigh_dc2[2] += sjk*dAdd*ddik_di[2];
1036 :
1037 : }
1038 : }
1039 19030 : dSASA_4_neigh_dc[0] = sji*dSASA_3_neigh_dc[0] + dSASA_4_neigh_dc2[0];
1040 19030 : dSASA_4_neigh_dc[1] = sji*dSASA_3_neigh_dc[1] + dSASA_4_neigh_dc2[1];
1041 19030 : dSASA_4_neigh_dc[2] = sji*dSASA_3_neigh_dc[2] + dSASA_4_neigh_dc2[2];
1042 :
1043 19030 : dSASA_3_neigh_dc[0] += dSASA_3_neigh_dc2[0];
1044 19030 : dSASA_3_neigh_dc[1] += dSASA_3_neigh_dc2[1];
1045 19030 : dSASA_3_neigh_dc[2] += dSASA_3_neigh_dc2[2];
1046 :
1047 19030 : dSASA_4_neigh_dc[0] += dSASA_2_neigh_dc[0] * Aikt;
1048 19030 : dSASA_4_neigh_dc[1] += dSASA_2_neigh_dc[1] * Aikt;
1049 19030 : dSASA_4_neigh_dc[2] += dSASA_2_neigh_dc[2] * Aikt;
1050 :
1051 19030 : if (AtomResidueName[0][Nlist[i][j]] == "N" || AtomResidueName[0][Nlist[i][j]] == "CA" || AtomResidueName[0][Nlist[i][j]] == "C" || AtomResidueName[0][Nlist[i][j]] == "O") {
1052 15720 : derivatives[i][0] += ((dSASA_2_neigh_dc[0]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[0]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[0]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][0]*DeltaG[natoms][0])*10;
1053 15720 : derivatives[i][1] += ((dSASA_2_neigh_dc[1]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[1]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[1]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][0]*DeltaG[natoms][0])*10;
1054 15720 : derivatives[i][2] += ((dSASA_2_neigh_dc[2]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[2]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[2]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][0]*DeltaG[natoms][0])*10;
1055 : }
1056 :
1057 19030 : if (AtomResidueName[0][Nlist[i][j]] != "N" && AtomResidueName[0][Nlist[i][j]] != "CA" && AtomResidueName[0][Nlist[i][j]] != "C" && AtomResidueName[0][Nlist[i][j]] != "O") {
1058 3310 : derivatives[i][0] += ((dSASA_2_neigh_dc[0]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[0]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[0]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][1]*DeltaG[Nlist[i][j]][0])*10;
1059 3310 : derivatives[i][1] += ((dSASA_2_neigh_dc[1]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[1]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[1]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][1]*DeltaG[Nlist[i][j]][0])*10;
1060 3310 : derivatives[i][2] += ((dSASA_2_neigh_dc[2]*LCPOparam[Nlist[i][j]][2] + dSASA_3_neigh_dc[2]*LCPOparam[Nlist[i][j]][3]+dSASA_4_neigh_dc[2]*LCPOparam[Nlist[i][j]][4])/MaxSurf[Nlist[i][j]][1]*DeltaG[Nlist[i][j]][0])*10;
1061 : }
1062 :
1063 19030 : Aijk += (Aijt * Ajkt);
1064 19030 : Aij += Aijt;
1065 19030 : Ajk += Ajkt;
1066 :
1067 19030 : dAijdc_2[0] += dAijdc_2t[0];
1068 19030 : dAijdc_2[1] += dAijdc_2t[1];
1069 19030 : dAijdc_2[2] += dAijdc_2t[2];
1070 :
1071 19030 : dAijdc_4[0] += Ajkt*dAijdc_2t[0];
1072 19030 : dAijdc_4[1] += Ajkt*dAijdc_2t[1];
1073 19030 : dAijdc_4[2] += Ajkt*dAijdc_2t[2];
1074 :
1075 19030 : }
1076 720 : double sasai = (LCPOparam[i][1]*S1+LCPOparam[i][2]*Aij+LCPOparam[i][3]*Ajk+LCPOparam[i][4]*Aijk);
1077 :
1078 720 : if (AtomResidueName[0][i] == "N" || AtomResidueName[0][i] == "CA" || AtomResidueName[0][i] == "C" || AtomResidueName[0][i] == "O") {
1079 576 : if (sasai > 0 ) {
1080 571 : sasa += (sasai/MaxSurf[i][0]*DeltaG[natoms][0]);
1081 : }
1082 576 : derivatives[i][0] += ((dAijdc_2[0]*LCPOparam[i][2]+dAijdc_4[0]*LCPOparam[i][4])/MaxSurf[i][0]*DeltaG[natoms][0])*10;
1083 576 : derivatives[i][1] += ((dAijdc_2[1]*LCPOparam[i][2]+dAijdc_4[1]*LCPOparam[i][4])/MaxSurf[i][0]*DeltaG[natoms][0])*10;
1084 576 : derivatives[i][2] += ((dAijdc_2[2]*LCPOparam[i][2]+dAijdc_4[2]*LCPOparam[i][4])/MaxSurf[i][0]*DeltaG[natoms][0])*10;
1085 : }
1086 :
1087 720 : if (AtomResidueName[0][i] != "N" && AtomResidueName[0][i] != "CA" && AtomResidueName[0][i] != "C" && AtomResidueName[0][i] != "O") {
1088 144 : if (sasai > 0. ) {
1089 144 : sasa += (sasai/MaxSurf[i][1]*DeltaG[i][0]);
1090 : }
1091 144 : derivatives[i][0] += ((dAijdc_2[0]*LCPOparam[i][2]+dAijdc_4[0]*LCPOparam[i][4])/MaxSurf[i][1]*DeltaG[i][0])*10;
1092 144 : derivatives[i][1] += ((dAijdc_2[1]*LCPOparam[i][2]+dAijdc_4[1]*LCPOparam[i][4])/MaxSurf[i][1]*DeltaG[i][0])*10;
1093 144 : derivatives[i][2] += ((dAijdc_2[2]*LCPOparam[i][2]+dAijdc_4[2]*LCPOparam[i][4])/MaxSurf[i][1]*DeltaG[i][0])*10;
1094 : }
1095 720 : }
1096 : }
1097 2880 : virial -= Tensor(getPosition(i),derivatives[i]);
1098 : }
1099 : }
1100 :
1101 :
1102 2904 : for(unsigned i=0; i<natoms; i++) {
1103 2880 : setAtomsDerivatives(i,derivatives[i]);
1104 : }
1105 24 : setBoxDerivatives(virial);
1106 24 : setValue(sasa);
1107 24 : firstStepFlag = 1;
1108 24 : ++nl_update;
1109 24 : if (nl_update == stride) {
1110 24 : nl_update = 0;
1111 : }
1112 48 : }
1113 :
1114 : }//namespace sasa
1115 : }//namespace PLMD
|
