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