Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "MetainferenceBase.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Pbc.h"
25 : #include "tools/Torsion.h"
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace isdb {
31 :
32 : //+PLUMEDOC ISDB_COLVAR JCOUPLING
33 : /*
34 : Calculates \f$^3J\f$ coupling constants for a dihedral angle.
35 :
36 : The J-coupling between two atoms is given by the Karplus relation:
37 :
38 : \f[
39 : ^3J(\theta)=A\cos^2(\theta+\Delta\theta)+B\cos(\theta+\Delta\theta)+C
40 : \f]
41 :
42 : where \f$A\f$, \f$B\f$ and \f$C\f$ are the Karplus parameters and \f$\Delta\theta\f$ is an additional constant
43 : added on to the dihedral angle \f$\theta\f$. The Karplus parameters are determined empirically and are dependent
44 : on the type of J-coupling.
45 :
46 : This collective variable computes the J-couplings for a set of atoms defining a dihedral angle. You can specify
47 : the atoms involved using the \ref MOLINFO notation. You can also specify the experimental couplings using the
48 : ADDCOUPLINGS flag and COUPLING keywords. These will be included in the output. You must choose the type of
49 : coupling using the type keyword, you can also supply custom Karplus parameters using TYPE=CUSTOM and the A, B, C
50 : and SHIFT keywords. You will need to make sure you are using the correct dihedral angle:
51 :
52 : - Ha-N: \f$\psi\f$
53 : - Ha-HN: \f$\phi\f$
54 : - N-C\f$\gamma\f$: \f$\chi_1\f$
55 : - CO-C\f$\gamma\f$: \f$\chi_1\f$
56 :
57 : J-couplings can be used to calculate a Metainference score using the internal keyword DOSCORE and all the options
58 : of \ref METAINFERENCE .
59 :
60 : \par Examples
61 :
62 : In the following example we calculate the Ha-N J-coupling from a set of atoms involved in
63 : dihedral \f$\psi\f$ angles in the peptide backbone. We also add the experimental datapoints and compute
64 : the correlation and other measures and finally print the results.
65 :
66 : \plumedfile
67 :
68 : MOLINFO MOLTYPE=protein STRUCTURE=peptide.pdb
69 : WHOLEMOLECULES ENTITY0=1-111
70 :
71 : JCOUPLING ...
72 : ADDCOUPLINGS
73 : TYPE=HAN
74 : ATOMS1=@psi-2 COUPLING1=-0.49
75 : ATOMS2=@psi-4 COUPLING2=-0.54
76 : ATOMS3=@psi-5 COUPLING3=-0.53
77 : ATOMS4=@psi-7 COUPLING4=-0.39
78 : ATOMS5=@psi-8 COUPLING5=-0.39
79 : LABEL=jhan
80 : ... JCOUPLING
81 :
82 : jhanst: STATS ARG=(jhan\.j_.*) PARARG=(jhan\.exp_.*)
83 :
84 : PRINT ARG=jhanst.*,jhan.* FILE=COLVAR STRIDE=100
85 : \endplumedfile
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 12 : class JCoupling :
91 : public MetainferenceBase
92 : {
93 : private:
94 : bool pbc;
95 : enum { HAN, HAHN, CCG, NCG, CUSTOM };
96 : unsigned ncoupl_;
97 : double ka_;
98 : double kb_;
99 : double kc_;
100 : double kshift_;
101 :
102 : public:
103 : static void registerKeywords(Keywords& keys);
104 : explicit JCoupling(const ActionOptions&);
105 : void calculate();
106 : void update();
107 : };
108 :
109 6458 : PLUMED_REGISTER_ACTION(JCoupling, "JCOUPLING")
110 :
111 7 : void JCoupling::registerKeywords(Keywords& keys) {
112 7 : componentsAreNotOptional(keys);
113 7 : useCustomisableComponents(keys);
114 7 : MetainferenceBase::registerKeywords(keys);
115 21 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
116 28 : keys.add("numbered", "ATOMS", "the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling. "
117 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one J-coupling will be "
118 : "calculated for each ATOMS keyword you specify.");
119 21 : keys.reset_style("ATOMS", "atoms");
120 21 : keys.addFlag("ADDCOUPLINGS", false, "Set this flag if you want to have fixed components with the experimental values.");
121 28 : keys.add("compulsory", "TYPE", "Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM)");
122 28 : keys.add("optional", "A", "Karplus parameter A");
123 28 : keys.add("optional", "B", "Karplus parameter B");
124 28 : keys.add("optional", "C", "Karplus parameter C");
125 28 : keys.add("optional", "SHIFT", "Angle shift in radians");
126 28 : keys.add("numbered", "COUPLING", "Add an experimental value for each coupling");
127 28 : keys.addOutputComponent("j", "default", "the calculated J-coupling");
128 28 : keys.addOutputComponent("exp", "ADDCOUPLINGS", "the experimental J-coupling");
129 7 : }
130 :
131 6 : JCoupling::JCoupling(const ActionOptions&ao):
132 : PLUMED_METAINF_INIT(ao),
133 6 : pbc(true)
134 : {
135 6 : bool nopbc = !pbc;
136 12 : parseFlag("NOPBC", nopbc);
137 6 : pbc =! nopbc;
138 :
139 : // Read in the atoms
140 : vector<AtomNumber> t, atoms;
141 28 : for (int i = 1; ; ++i) {
142 68 : parseAtomList("ATOMS", i, t );
143 34 : if (t.empty()) {
144 : break;
145 : }
146 :
147 28 : if (t.size() != 4) {
148 : std::string ss;
149 0 : Tools::convert(i, ss);
150 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
151 : }
152 :
153 : // This makes the distance calculation easier later on (see Torsion implementation)
154 28 : atoms.push_back(t[0]);
155 28 : atoms.push_back(t[1]);
156 28 : atoms.push_back(t[1]);
157 28 : atoms.push_back(t[2]);
158 28 : atoms.push_back(t[2]);
159 28 : atoms.push_back(t[3]);
160 28 : t.resize(0);
161 28 : }
162 :
163 : // We now have 6 atoms per datapoint
164 6 : ncoupl_ = atoms.size()/6;
165 :
166 : // Parse J-Coupling type, this will determine the Karplus parameters
167 : unsigned jtype_ = CUSTOM;
168 : string string_type;
169 12 : parse("TYPE", string_type);
170 6 : if(string_type == "HAN") {
171 : jtype_ = HAN;
172 5 : } else if(string_type == "HAHN") {
173 : jtype_ = HAHN;
174 2 : } else if(string_type == "CCG") {
175 : jtype_ = CCG;
176 1 : } else if(string_type == "NCG") {
177 : jtype_ = NCG;
178 0 : } else if(string_type == "CUSTOM") {
179 : jtype_ = CUSTOM;
180 : } else {
181 0 : error("Unknown J-coupling type!");
182 : }
183 :
184 : // Optionally add an experimental value (like with RDCs)
185 : vector<double> coupl;
186 6 : bool addcoupling = false;
187 12 : parseFlag("ADDCOUPLINGS", addcoupling);
188 12 : if (addcoupling||getDoScore()) {
189 1 : coupl.resize(ncoupl_);
190 : unsigned ntarget = 0;
191 8 : for (unsigned i = 0; i < ncoupl_; ++i) {
192 21 : if (!parseNumbered("COUPLING", i+1, coupl[i])) {
193 : break;
194 : }
195 : ntarget++;
196 : }
197 1 : if (ntarget != ncoupl_) {
198 0 : error("found wrong number of COUPLING values");
199 : }
200 : }
201 :
202 : // For custom types we allow use of custom Karplus parameters
203 6 : if (jtype_ == CUSTOM) {
204 0 : parse("A", ka_);
205 0 : parse("B", kb_);
206 0 : parse("C", kc_);
207 0 : parse("SHIFT", kshift_);
208 : }
209 :
210 6 : log << " Bibliography ";
211 :
212 : // Set Karplus parameters
213 6 : switch (jtype_) {
214 1 : case HAN:
215 1 : ka_ = -0.88;
216 1 : kb_ = -0.61;
217 1 : kc_ = -0.27;
218 1 : kshift_ = pi / 3.0;
219 1 : log.printf("J-coupling type is HAN, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
220 3 : log << plumed.cite("Wang A C, Bax A, J. Am. Chem. Soc. 117, 1810 (1995)");
221 1 : break;
222 3 : case HAHN:
223 3 : ka_ = 7.09;
224 3 : kb_ = -1.42;
225 3 : kc_ = 1.55;
226 3 : kshift_ = -pi / 3.0;
227 3 : log.printf("J-coupling type is HAHN, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
228 9 : log << plumed.cite("Hu J-S, Bax A, J. Am. Chem. Soc. 119, 6360 (1997)");
229 3 : break;
230 1 : case CCG:
231 1 : ka_ = 2.31;
232 1 : kb_ = -0.87;
233 1 : kc_ = 0.55;
234 1 : kshift_ = (2.0 * pi) / 3.0;
235 1 : log.printf("J-coupling type is CCG, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
236 3 : log << plumed.cite("Perez C, Löhr F, Rüterjans H, Schmidt J, J. Am. Chem. Soc. 123, 7081 (2001)");
237 1 : break;
238 1 : case NCG:
239 1 : ka_ = 1.29;
240 1 : kb_ = -0.49;
241 1 : kc_ = 0.37;
242 1 : kshift_ = 0.0;
243 1 : log.printf("J-coupling type is NCG, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
244 3 : log << plumed.cite("Perez C, Löhr F, Rüterjans H, Schmidt J, J. Am. Chem. Soc. 123, 7081 (2001)");
245 1 : break;
246 0 : case CUSTOM:
247 0 : log.printf("J-coupling type is custom, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
248 : break;
249 : }
250 18 : log<<plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
251 6 : log<<"\n";
252 :
253 34 : for (unsigned i = 0; i < ncoupl_; ++i) {
254 56 : log.printf(" The %uth J-Coupling is calculated from atoms : %d %d %d %d.",
255 112 : i+1, atoms[6*i].serial(), atoms[6*i+1].serial(), atoms[6*i+3].serial(), atoms[6*i+5].serial());
256 28 : if (addcoupling) {
257 0 : log.printf(" Experimental J-Coupling is %f.", coupl[i]);
258 : }
259 28 : log.printf("\n");
260 : }
261 :
262 6 : if (pbc) {
263 0 : log.printf(" using periodic boundary conditions\n");
264 : } else {
265 6 : log.printf(" without periodic boundary conditions\n");
266 : }
267 :
268 6 : if(!getDoScore()) {
269 47 : for (unsigned i = 0; i < ncoupl_; i++) {
270 21 : std::string num; Tools::convert(i, num);
271 42 : addComponentWithDerivatives("j_" + num);
272 42 : componentIsNotPeriodic("j_" + num);
273 : }
274 : } else {
275 15 : for (unsigned i = 0; i < ncoupl_; i++) {
276 7 : std::string num; Tools::convert(i, num);
277 14 : addComponent("j_" + num);
278 14 : componentIsNotPeriodic("j_" + num);
279 : }
280 : }
281 :
282 12 : if (addcoupling||getDoScore()) {
283 15 : for (unsigned i = 0; i < ncoupl_; i++) {
284 7 : std::string num; Tools::convert(i, num);
285 14 : addComponent("exp_" + num);
286 14 : componentIsNotPeriodic("exp_" + num);
287 14 : Value* comp = getPntrToComponent("exp_" + num);
288 14 : comp->set(coupl[i]);
289 : }
290 : }
291 :
292 6 : requestAtoms(atoms);
293 6 : if(getDoScore()) {
294 1 : setParameters(coupl);
295 1 : Initialise(ncoupl_);
296 : }
297 6 : setDerivatives();
298 6 : checkRead();
299 6 : }
300 :
301 16 : void JCoupling::calculate()
302 : {
303 16 : if (pbc) makeWhole();
304 16 : vector<Vector> deriv(ncoupl_*6);
305 16 : vector<double> j(ncoupl_,0.);
306 :
307 48 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
308 : {
309 : #pragma omp for
310 : // Loop through atoms, with steps of 6 atoms (one iteration per datapoint)
311 32 : for (unsigned r=0; r<ncoupl_; r++) {
312 : // Index is the datapoint index
313 98 : unsigned a0 = 6*r;
314 :
315 : // 6 atoms -> 3 vectors
316 294 : Vector d0 = delta(getPosition(a0+1), getPosition(a0));
317 294 : Vector d1 = delta(getPosition(a0+3), getPosition(a0+2));
318 294 : Vector d2 = delta(getPosition(a0+5), getPosition(a0+4));
319 :
320 : // Calculate dihedral with 3 vectors, get the derivatives
321 98 : Vector dd0, dd1, dd2;
322 : PLMD::Torsion t;
323 98 : double torsion = t.compute(d0, d1, d2, dd0, dd1, dd2);
324 :
325 : // Calculate the Karplus relation and its derivative
326 98 : double theta = torsion + kshift_;
327 98 : double cos_theta = cos(theta);
328 98 : double sin_theta = sin(theta);
329 196 : j[r] = ka_*cos_theta*cos_theta + kb_*cos_theta + kc_;
330 98 : double derj = -2.*ka_*sin_theta*cos_theta - kb_*sin_theta;
331 :
332 98 : dd0 *= derj;
333 98 : dd1 *= derj;
334 98 : dd2 *= derj;
335 :
336 182 : if(getDoScore()) setCalcData(r, j[r]);
337 196 : deriv[a0] = dd0;
338 196 : deriv[a0+1] = -dd0;
339 196 : deriv[a0+2] = dd1;
340 196 : deriv[a0+3] = -dd1;
341 196 : deriv[a0+4] = dd2;
342 196 : deriv[a0+5] = -dd2;
343 : }
344 : }
345 :
346 16 : if(getDoScore()) {
347 : /* Metainference */
348 6 : double score = getScore();
349 : setScore(score);
350 :
351 : /* calculate final derivatives */
352 6 : Tensor virial;
353 12 : Value* val=getPntrToComponent("score");
354 90 : for (unsigned r=0; r<ncoupl_; r++) {
355 42 : const unsigned a0 = 6*r;
356 84 : setAtomsDerivatives(val, a0, deriv[a0]*getMetaDer(r));
357 84 : setAtomsDerivatives(val, a0+1, deriv[a0+1]*getMetaDer(r));
358 84 : setAtomsDerivatives(val, a0+2, deriv[a0+2]*getMetaDer(r));
359 84 : setAtomsDerivatives(val, a0+3, deriv[a0+3]*getMetaDer(r));
360 84 : setAtomsDerivatives(val, a0+4, deriv[a0+4]*getMetaDer(r));
361 84 : setAtomsDerivatives(val, a0+5, deriv[a0+5]*getMetaDer(r));
362 84 : virial-=Tensor(getPosition(a0), deriv[a0]*getMetaDer(r));
363 84 : virial-=Tensor(getPosition(a0+1), deriv[a0+1]*getMetaDer(r));
364 84 : virial-=Tensor(getPosition(a0+2), deriv[a0+2]*getMetaDer(r));
365 84 : virial-=Tensor(getPosition(a0+3), deriv[a0+3]*getMetaDer(r));
366 84 : virial-=Tensor(getPosition(a0+4), deriv[a0+4]*getMetaDer(r));
367 84 : virial-=Tensor(getPosition(a0+5), deriv[a0+5]*getMetaDer(r));
368 : }
369 6 : setBoxDerivatives(val, virial);
370 : } else {
371 122 : for (unsigned r=0; r<ncoupl_; r++) {
372 56 : const unsigned a0 = 6*r;
373 56 : string num; Tools::convert(r,num);
374 112 : Value* val=getPntrToComponent("j_"+num);
375 112 : val->set(j[r]);
376 112 : setAtomsDerivatives(val, a0, deriv[a0]);
377 112 : setAtomsDerivatives(val, a0+1, deriv[a0+1]);
378 112 : setAtomsDerivatives(val, a0+2, deriv[a0+2]);
379 112 : setAtomsDerivatives(val, a0+3, deriv[a0+3]);
380 112 : setAtomsDerivatives(val, a0+4, deriv[a0+4]);
381 112 : setAtomsDerivatives(val, a0+5, deriv[a0+5]);
382 56 : Tensor virial;
383 112 : virial-=Tensor(getPosition(a0), deriv[a0]);
384 112 : virial-=Tensor(getPosition(a0+1), deriv[a0+1]);
385 112 : virial-=Tensor(getPosition(a0+2), deriv[a0+2]);
386 112 : virial-=Tensor(getPosition(a0+3), deriv[a0+3]);
387 112 : virial-=Tensor(getPosition(a0+4), deriv[a0+4]);
388 112 : virial-=Tensor(getPosition(a0+5), deriv[a0+5]);
389 56 : setBoxDerivatives(val, virial);
390 : }
391 : }
392 16 : }
393 :
394 16 : void JCoupling::update() {
395 : // write status file
396 32 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
397 16 : }
398 :
399 : }
400 4839 : }
|
