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