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