Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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 "Colvar.h"
23 : #include "ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Torsion.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace colvar {
34 :
35 : //+PLUMEDOC COLVAR JCOUPLING
36 : /*
37 : Calculates \f$^3J\f$ coupling constants for a dihedral angle.
38 :
39 : The J-coupling between two atoms is given by the Karplus relation:
40 :
41 : \f[
42 : ^3J(\theta)=A\cos^2(\theta+\Delta\theta)+B\cos(\theta+\Delta\theta)+C
43 : \f]
44 :
45 : where \f$A\f$, \f$B\f$ and \f$C\f$ are the Karplus parameters and \f$\Delta\theta\f$ is an additional constant
46 : added on to the dihedral angle \f$\theta\f$. The Karplus parameters are determined empirically and are dependent
47 : on the type of J-coupling.
48 :
49 : This collective variable computes the J-couplings for a set of atoms defining a dihedral angle. You can specify
50 : the atoms involved using the \ref MOLINFO notation. You can also specify the experimental couplings using the
51 : ADDCOUPLINGS flag and COUPLING keywords. These will be included in the output. You must choose the type of
52 : coupling using the type keyword, you can also supply custom Karplus parameters using TYPE=CUSTOM and the A, B, C
53 : and SHIFT keywords. You will need to make sure you are using the correct dihedral angle:
54 :
55 : - Ha-N: \f$\psi\f$
56 : - Ha-HN: \f$\phi\f$
57 : - N-C\f$\gamma\f$: \f$\chi_1\f$
58 : - CO-C\f$\gamma\f$: \f$\chi_1\f$
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 : \verbatim
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 :
86 : ENDPLUMED
87 :
88 : \endverbatim
89 : (See also \ref PRINT, \ref STATS)
90 :
91 : */
92 : //+ENDPLUMEDOC
93 :
94 16 : class JCoupling : public Colvar {
95 : private:
96 : bool pbc;
97 : vector<double> coupl;
98 : unsigned ndata;
99 : unsigned jtype_;
100 : enum { HAN, HAHN, CCG, NCG, CUSTOM };
101 : double ka_;
102 : double kb_;
103 : double kc_;
104 : double kshift_;
105 :
106 : public:
107 : explicit JCoupling(const ActionOptions&);
108 : virtual void calculate();
109 : static void registerKeywords(Keywords& keys);
110 : };
111 :
112 2531 : PLUMED_REGISTER_ACTION(JCoupling, "JCOUPLING")
113 :
114 9 : void JCoupling::registerKeywords(Keywords& keys) {
115 9 : Colvar::registerKeywords(keys);
116 9 : componentsAreNotOptional(keys);
117 9 : useCustomisableComponents(keys);
118 : keys.add("numbered", "ATOMS", "the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling. "
119 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one J-coupling will be "
120 9 : "calculated for each ATOMS keyword you specify.");
121 9 : keys.reset_style("ATOMS", "atoms");
122 9 : keys.addFlag("ADDCOUPLINGS", false, "Set this flag if you want to have fixed components with the experimental values.");
123 9 : keys.add("compulsory", "TYPE", "Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM)");
124 9 : keys.add("optional", "A", "Karplus parameter A");
125 9 : keys.add("optional", "B", "Karplus parameter B");
126 9 : keys.add("optional", "C", "Karplus parameter C");
127 9 : keys.add("optional", "SHIFT", "Angle shift in radians");
128 9 : keys.add("numbered", "COUPLING", "Add an experimental value for each coupling");
129 9 : keys.addOutputComponent("j", "default", "the calculated J-coupling");
130 9 : keys.addOutputComponent("exp", "ADDCOUPLINGS", "the experimental J-coupling");
131 9 : }
132 :
133 8 : JCoupling::JCoupling(const ActionOptions&ao):
134 : PLUMED_COLVAR_INIT(ao),
135 : pbc(true),
136 8 : ndata(0)
137 : {
138 8 : bool nopbc = !pbc;
139 8 : parseFlag("NOPBC", nopbc);
140 8 : pbc =! nopbc;
141 :
142 : // Read in the atoms
143 16 : vector<AtomNumber> t, atoms;
144 36 : for (int i = 1; ; ++i) {
145 36 : parseAtomList("ATOMS", i, t );
146 36 : if (t.empty()) {
147 8 : break;
148 : }
149 :
150 28 : if (t.size() != 4) {
151 0 : std::string ss;
152 0 : Tools::convert(i, ss);
153 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
154 : }
155 :
156 : // This makes the distance calculation easier later on (see Torsion implementation)
157 28 : atoms.push_back(t[0]);
158 28 : atoms.push_back(t[1]);
159 28 : atoms.push_back(t[1]);
160 28 : atoms.push_back(t[2]);
161 28 : atoms.push_back(t[2]);
162 28 : atoms.push_back(t[3]);
163 28 : t.resize(0);
164 28 : }
165 :
166 : // We now have 6 atoms per datapoint
167 8 : ndata = atoms.size()/6;
168 :
169 : // Parse J-Coupling type, this will determine the Karplus parameters
170 16 : string string_type;
171 8 : parse("TYPE", string_type);
172 8 : if(string_type == "HAN") {
173 2 : jtype_ = HAN;
174 6 : } else if(string_type == "HAHN") {
175 2 : jtype_ = HAHN;
176 4 : } else if(string_type == "CCG") {
177 2 : jtype_ = CCG;
178 2 : } else if(string_type == "NCG") {
179 2 : jtype_ = NCG;
180 0 : } else if(string_type == "CUSTOM") {
181 0 : jtype_ = CUSTOM;
182 : } else {
183 0 : error("Unknown J-coupling type!");
184 : }
185 :
186 : // Optionally add an experimental value (like with RDCs)
187 8 : bool addcoupling = false;
188 8 : parseFlag("ADDCOUPLINGS", addcoupling);
189 8 : if (addcoupling) {
190 0 : coupl.resize(ndata);
191 0 : unsigned ntarget = 0;
192 0 : for (unsigned i = 0; i < ndata; ++i) {
193 0 : if (!parseNumbered("COUPLING", i+1, coupl[i])) {
194 0 : break;
195 : }
196 0 : ntarget++;
197 : }
198 0 : if (ntarget != ndata) {
199 0 : error("found wrong number of COUPLING values");
200 : }
201 : }
202 :
203 : // For custom types we allow use of custom Karplus parameters
204 8 : if (jtype_ == CUSTOM) {
205 0 : parse("A", ka_);
206 0 : parse("B", kb_);
207 0 : parse("C", kc_);
208 0 : parse("SHIFT", kshift_);
209 : }
210 :
211 8 : checkRead();
212 :
213 : // Set Karplus parameters
214 8 : switch (jtype_) {
215 : case HAN:
216 2 : ka_ = -0.88;
217 2 : kb_ = -0.61;
218 2 : kc_ = -0.27;
219 2 : kshift_ = pi / 3.0;
220 2 : log.printf("J-coupling type is HAN, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
221 2 : log << " Bibliography "
222 4 : << plumed.cite("Wang A C, Bax A, J. Am. Chem. Soc. 117, 1810 (1995)") << "\n";
223 2 : break;
224 : case HAHN:
225 2 : ka_ = 7.09;
226 2 : kb_ = -1.42;
227 2 : kc_ = 1.55;
228 2 : kshift_ = -pi / 3.0;
229 2 : log.printf("J-coupling type is HAHN, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
230 2 : log << " Bibliography "
231 4 : << plumed.cite("Hu J-S, Bax A, J. Am. Chem. Soc. 119, 6360 (1997)") << "\n";
232 2 : break;
233 : case CCG:
234 2 : ka_ = 2.31;
235 2 : kb_ = -0.87;
236 2 : kc_ = 0.55;
237 2 : kshift_ = (2.0 * pi) / 3.0;
238 2 : log.printf("J-coupling type is CCG, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
239 2 : log << " Bibliography "
240 4 : << plumed.cite("Perez C, Löhr F, Rüterjans H, Schmidt J, J. Am. Chem. Soc. 123, 7081 (2001)") << "\n";
241 2 : break;
242 : case NCG:
243 2 : ka_ = 1.29;
244 2 : kb_ = -0.49;
245 2 : kc_ = 0.37;
246 2 : kshift_ = 0.0;
247 2 : log.printf("J-coupling type is NCG, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
248 2 : log << " Bibliography "
249 4 : << plumed.cite("Perez C, Löhr F, Rüterjans H, Schmidt J, J. Am. Chem. Soc. 123, 7081 (2001)") << "\n";
250 2 : break;
251 : case CUSTOM:
252 0 : log.printf("J-coupling type is custom, with A: %f, B: %f, C: %f, angle shift: %f\n", ka_, kb_, kc_, kshift_);
253 0 : break;
254 : }
255 :
256 36 : for (unsigned i = 0; i < ndata; ++i) {
257 : log.printf(" The %uth J-Coupling is calculated from atoms : %d %d %d %d.",
258 28 : i+1, atoms[2*i].serial(), atoms[2*i+1].serial(), atoms[2*i+2].serial(), atoms[2*i+3].serial());
259 28 : if (addcoupling) {
260 0 : log.printf(" Experimental J-Coupling is %f.", coupl[i]);
261 : }
262 28 : log.printf("\n");
263 : }
264 :
265 8 : if (pbc) {
266 8 : log.printf(" using periodic boundary conditions\n");
267 : } else {
268 0 : log.printf(" without periodic boundary conditions\n");
269 : }
270 :
271 36 : for (unsigned i = 0; i < ndata; i++) {
272 28 : std::string num; Tools::convert(i, num);
273 28 : addComponentWithDerivatives("j_" + num);
274 28 : componentIsNotPeriodic("j_" + num);
275 28 : }
276 :
277 8 : if (addcoupling) {
278 0 : for (unsigned i = 0; i < ndata; i++) {
279 0 : std::string num; Tools::convert(i, num);
280 0 : addComponent("exp_" + num);
281 0 : componentIsNotPeriodic("exp_" + num);
282 0 : Value* comp = getPntrToComponent("exp_" + num);
283 0 : comp->set(coupl[i]);
284 0 : }
285 : }
286 :
287 16 : requestAtoms(atoms);
288 8 : }
289 :
290 1776 : void JCoupling::calculate() {
291 1776 : if (pbc) {
292 1776 : makeWhole();
293 : }
294 :
295 : // Loop through atoms, with steps of 6 atoms (one iteration per datapoint)
296 10908 : for (unsigned r = 0; r < (ndata * 6); r += 6) {
297 : // Index is the datapoint index
298 9132 : const unsigned index = r / 6;
299 :
300 : // 6 atoms -> 3 vectors
301 9132 : Vector d0, d1, d2;
302 9132 : d0 = delta(getPosition(r + 1), getPosition(r + 0));
303 9132 : d1 = delta(getPosition(r + 3), getPosition(r + 2));
304 9132 : d2 = delta(getPosition(r + 5), getPosition(r + 4));
305 :
306 : // Calculate dihedral with 3 vectors, get the derivatives
307 9132 : Vector dd0, dd1, dd2;
308 : PLMD::Torsion t;
309 9132 : const double torsion = t.compute(d0, d1, d2, dd0, dd1, dd2);
310 :
311 : // Calculate the Karplus relation and its derivative
312 9132 : const double theta = torsion + kshift_;
313 9132 : const double cos_theta = cos(theta);
314 9132 : const double j = ka_ * cos_theta * cos_theta + kb_ * cos_theta + kc_;
315 9132 : const double derj = sin(theta) * (-1.0 * (2.0 * ka_ * cos_theta + kb_));
316 :
317 9132 : Value* val = getPntrToComponent(index);
318 9132 : val->set(j);
319 :
320 9132 : setAtomsDerivatives(val, r + 0, derj * dd0);
321 9132 : setAtomsDerivatives(val, r + 1, derj * -dd0);
322 9132 : setAtomsDerivatives(val, r + 2, derj * dd1);
323 9132 : setAtomsDerivatives(val, r + 3, derj * -dd1);
324 9132 : setAtomsDerivatives(val, r + 4, derj * dd2);
325 9132 : setAtomsDerivatives(val, r + 5, derj * -dd2);
326 9132 : setBoxDerivativesNoPbc(val);
327 : }
328 1776 : }
329 : }
330 2523 : }
331 :
332 :
333 :
|
