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