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