Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 :
26 : #ifdef __PLUMED_HAS_GSL
27 : #include <gsl/gsl_vector.h>
28 : #include <gsl/gsl_matrix.h>
29 : #include <gsl/gsl_linalg.h>
30 : #include <gsl/gsl_blas.h>
31 : #endif
32 :
33 : using namespace std;
34 :
35 : namespace PLMD {
36 : namespace isdb {
37 :
38 : //+PLUMEDOC ISDB_COLVAR RDC
39 : /*
40 : Calculates the (Residual) Dipolar Coupling between two atoms.
41 :
42 : The Dipolar Coupling between two nuclei depends on the \f$\theta\f$ angle between
43 : the inter-nuclear vector and the external magnetic field.
44 :
45 : \f[
46 : D=D_{max}0.5(3\cos^2(\theta)-1)
47 : \f]
48 :
49 : where
50 :
51 : \f[
52 : D_{max}=-\mu_0\gamma_1\gamma_2h/(8\pi^3r^3)
53 : \f]
54 :
55 : that is the maximal value of the dipolar coupling for the two nuclear spins with gyromagnetic ratio \f$\gamma\f$.
56 : \f$\mu\f$ is the magnetic constant and h is the Planck constant.
57 :
58 : Common Gyromagnetic Ratios (C.G.S)
59 : - H(1) 26.7513
60 : - C(13) 6.7261
61 : - N(15) -2.7116
62 : and their products (this is what is given in input using the keyword GYROM)
63 : - N-H -72.5388
64 : - C-H 179.9319
65 : - C-N -18.2385
66 : - C-C 45.2404
67 :
68 : In isotropic media DCs average to zero because of the rotational
69 : averaging, but when the rotational symmetry is broken, either through the introduction of an alignment medium or for molecules
70 : with highly anisotropic paramagnetic susceptibility, then the average of the DCs is not zero and it is possible to measure a Residual Dipolar Coupling (RDCs).
71 :
72 : This collective variable calculates the Dipolar Coupling for a set of couple of atoms using
73 : the above definition.
74 :
75 : In a standard MD simulation the average over time of the DC should then be zero. If one wants to model the meaning of a set of measured RDCs it is possible to try to solve the following problem: "what is the distribution of structures and orientations that reproduce the measured RDCs".
76 :
77 : This collective variable can then be use to break the rotational symmetry of a simulation by imposing that the average of the DCs over the conformational ensemble must be equal to the measured RDCs \cite Camilloni:2015ka . Since measured RDCs are also a function of the fraction of aligned molecules in the sample it is better to compare them modulo a constant or looking at the correlation.
78 :
79 : Alternatively if the molecule is rigid it is possible to use the experimental data to calculate the alignment tensor and the use that to back calculate the RDCs, this is what is usually call the Single Value Decomposition approach. In this case the code rely on the
80 : a set of function from the GNU Scientific Library (GSL). (With SVD forces are not currently implemented).
81 :
82 : Replica-Averaged simulations can be perfomed using RDCs, \ref ENSEMBLE, \ref STATS and \ref RESTRAINT .
83 : \ref METAINFERENCE can be activated using DOSCORE and the other relevant keywords.
84 :
85 : Additional material and examples can be also found in the tutorial \ref belfast-9
86 :
87 : \par Examples
88 : In the following example five N-H RDCs are defined and averaged over multiple replicas,
89 : their correlation is then calculated with respect to a set of experimental data and restrained.
90 : In addition, and only for analysis purposes, the same RDCs each single conformation are calculated
91 : using a Single Value Decomposition algorithm, then averaged and again compared with the experimenta data.
92 :
93 : \plumedfile
94 : RDC ...
95 : GYROM=-72.5388
96 : SCALE=0.001
97 : ATOMS1=20,21
98 : ATOMS2=37,38
99 : ATOMS3=56,57
100 : ATOMS4=76,77
101 : ATOMS5=92,93
102 : LABEL=nh
103 : ... RDC
104 :
105 : erdc: ENSEMBLE ARG=nh.*
106 :
107 : st: STATS ARG=erdc.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
108 :
109 : rdce: RESTRAINT ARG=st.corr KAPPA=0. SLOPE=-25000.0 AT=1.
110 :
111 : RDC ...
112 : GYROM=-72.5388
113 : SVD
114 : ATOMS1=20,21 COUPLING1=8.17
115 : ATOMS2=37,38 COUPLING2=-8.271
116 : ATOMS3=56,57 COUPLING3=-10.489
117 : ATOMS4=76,77 COUPLING4=-9.871
118 : ATOMS5=92,93 COUPLING5=-9.152
119 : LABEL=svd
120 : ... RDC
121 :
122 : esvd: ENSEMBLE ARG=svd.*
123 :
124 : st_svd: STATS ARG=esvd.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
125 :
126 :
127 : PRINT ARG=st.corr,st_svd.corr,rdce.bias FILE=colvar
128 : \endplumedfile
129 :
130 : */
131 : //+ENDPLUMEDOC
132 :
133 : //+PLUMEDOC ISDB_COLVAR PCS
134 : /*
135 : Calculates the Pseudocontact shift of a nucleus determined by the presence of a metal ion susceptible to anisotropic magnetization.
136 :
137 : The PCS of an atomic nucleus depends on the \f$\theta\f$ angle between the vector from the spin-label to the nucleus
138 : and the external magnetic field and the module of the vector itself \cite Camilloni:2015jf . While in principle the averaging
139 : resulting from the tumbling should remove the pseudocontact shift, in presence of the NMR magnetic field the magnatically anisotropic molecule bound to system will break the rotational symmetry does resulting in measurable PCSs and RDCs.
140 :
141 : PCSs can also be calculated using a Single Value Decomposition approach, in this case the code rely on the
142 : a set of function from the GNU Scientific Library (GSL). (With SVD forces are not currently implemented).
143 :
144 : Replica-Averaged simulations can be perfomed using PCSs, \ref ENSEMBLE, \ref STATS and \ref RESTRAINT .
145 : Metainference simulations can be performed with this CV and \ref METAINFERENCE .
146 :
147 : \par Examples
148 :
149 : In the following example five PCSs are defined and their correlation with
150 : respect to a set of experimental data is calculated and restrained. In addition,
151 : and only for analysis purposes, the same PCSs are calculated using a Single Value
152 : Decomposition algorithm.
153 :
154 : \plumedfile
155 : PCS ...
156 : ATOMS1=20,21
157 : ATOMS2=20,38
158 : ATOMS3=20,57
159 : ATOMS4=20,77
160 : ATOMS5=20,93
161 : LABEL=nh
162 : ... PCS
163 :
164 : enh: ENSEMBLE ARG=nh.*
165 :
166 : st: STATS ARG=enh.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
167 :
168 : pcse: RESTRAINT ARG=st.corr KAPPA=0. SLOPE=-25000.0 AT=1.
169 :
170 : PRINT ARG=st.corr,pcse.bias FILE=colvar
171 : \endplumedfile
172 :
173 : */
174 : //+ENDPLUMEDOC
175 :
176 58 : class RDC :
177 : public MetainferenceBase
178 : {
179 : private:
180 : double Const;
181 : double mu_s;
182 : double scale;
183 : vector<double> coupl;
184 : bool svd;
185 : bool pbc;
186 :
187 : void do_svd();
188 : public:
189 : explicit RDC(const ActionOptions&);
190 : static void registerKeywords( Keywords& keys );
191 : virtual void calculate();
192 : void update();
193 : };
194 :
195 4850 : PLUMED_REGISTER_ACTION(RDC,"RDC")
196 4821 : PLUMED_REGISTER_ACTION(RDC,"PCS")
197 :
198 31 : void RDC::registerKeywords( Keywords& keys ) {
199 31 : componentsAreNotOptional(keys);
200 31 : useCustomisableComponents(keys);
201 31 : MetainferenceBase::registerKeywords(keys);
202 31 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
203 : keys.add("numbered","ATOMS","the couple of atoms involved in each of the bonds for which you wish to calculate the RDC. "
204 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dipolar coupling will be "
205 31 : "calculated for each ATOMS keyword you specify.");
206 31 : keys.reset_style("ATOMS","atoms");
207 31 : keys.add("compulsory","GYROM","1.","Add the product of the gyromagnetic constants for the bond. ");
208 31 : keys.add("compulsory","SCALE","1.","Add the scaling factor to take into account concentration and other effects. ");
209 31 : keys.addFlag("SVD",false,"Set to TRUE if you want to backcalculate using Single Value Decomposition (need GSL at compilation time).");
210 31 : keys.addFlag("ADDCOUPLINGS",false,"Set to TRUE if you want to have fixed components with the experimetnal values.");
211 31 : keys.add("numbered","COUPLING","Add an experimental value for each coupling (needed by SVD and usefull for \ref STATS).");
212 31 : keys.addOutputComponent("rdc","default","the calculated # RDC");
213 31 : keys.addOutputComponent("exp","SVD/ADDCOUPLINGS","the experimental # RDC");
214 31 : }
215 :
216 29 : RDC::RDC(const ActionOptions&ao):
217 : PLUMED_METAINF_INIT(ao),
218 : Const(1.),
219 : mu_s(1.),
220 : scale(1.),
221 29 : pbc(true)
222 : {
223 29 : bool nopbc=!pbc;
224 29 : parseFlag("NOPBC",nopbc);
225 29 : pbc=!nopbc;
226 :
227 29 : const double RDCConst = 0.3356806;
228 29 : const double PCSConst = 1.0;
229 :
230 29 : if( getName().find("RDC")!=std::string::npos) { Const *= RDCConst; }
231 0 : else if( getName().find("PCS")!=std::string::npos) { Const *= PCSConst; }
232 :
233 : // Read in the atoms
234 58 : vector<AtomNumber> t, atoms;
235 146 : for(int i=1;; ++i ) {
236 146 : parseAtomList("ATOMS", i, t );
237 146 : if( t.empty() ) break;
238 117 : if( t.size()!=2 ) {
239 0 : std::string ss; Tools::convert(i,ss);
240 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
241 : }
242 117 : atoms.push_back(t[0]);
243 117 : atoms.push_back(t[1]);
244 117 : t.resize(0);
245 117 : }
246 :
247 29 : const unsigned ndata = atoms.size()/2;
248 :
249 : // Read in GYROMAGNETIC constants
250 29 : parse("GYROM", mu_s);
251 29 : if(mu_s==0.) error("GYROM cannot be 0");
252 :
253 : // Read in SCALING factors
254 29 : parse("SCALE", scale);
255 29 : if(scale==0.) error("SCALE cannot be 0");
256 :
257 29 : svd=false;
258 29 : parseFlag("SVD",svd);
259 : #ifndef __PLUMED_HAS_GSL
260 : if(svd) error("You CANNOT use SVD without GSL. Recompile PLUMED with GSL!\n");
261 : #endif
262 29 : if(svd&&getDoScore()) error("It is not possible to use SVD and METAINFERENCE together");
263 :
264 29 : bool addexp=false;
265 29 : parseFlag("ADDCOUPLINGS",addexp);
266 29 : if(getDoScore()||svd) addexp=true;
267 :
268 29 : if(addexp) {
269 13 : coupl.resize( ndata );
270 13 : unsigned ntarget=0;
271 66 : for(unsigned i=0; i<ndata; ++i) {
272 53 : if( !parseNumbered( "COUPLING", i+1, coupl[i] ) ) break;
273 53 : ntarget++;
274 : }
275 13 : if( ntarget!=ndata ) error("found wrong number of COUPLING values");
276 : }
277 :
278 : // Ouput details of all contacts
279 29 : log.printf(" Gyromagnetic moment is %f. Scaling factor is %f.",mu_s,scale);
280 146 : for(unsigned i=0; i<ndata; ++i) {
281 117 : log.printf(" The %uth Bond Dipolar Coupling is calculated from atoms : %d %d.", i+1, atoms[2*i].serial(), atoms[2*i+1].serial());
282 117 : if(addexp) log.printf(" Experimental coupling is %f.", coupl[i]);
283 117 : log.printf("\n");
284 : }
285 :
286 29 : log<<" Bibliography "
287 87 : <<plumed.cite("Camilloni C, Vendruscolo M, J. Phys. Chem. B 119, 653 (2015)")
288 87 : <<plumed.cite("Camilloni C, Vendruscolo M, Biochemistry 54, 7470 (2015)");
289 29 : log<< plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
290 :
291 :
292 29 : if(!getDoScore()&&!svd) {
293 80 : for(unsigned i=0; i<ndata; i++) {
294 64 : std::string num; Tools::convert(i,num);
295 64 : addComponentWithDerivatives("rdc_"+num);
296 64 : componentIsNotPeriodic("rdc_"+num);
297 64 : }
298 16 : if(addexp) {
299 0 : for(unsigned i=0; i<ndata; i++) {
300 0 : std::string num; Tools::convert(i,num);
301 0 : addComponent("exp_"+num);
302 0 : componentIsNotPeriodic("exp_"+num);
303 0 : Value* comp=getPntrToComponent("exp_"+num);
304 0 : comp->set(coupl[i]);
305 0 : }
306 : }
307 : } else {
308 66 : for(unsigned i=0; i<ndata; i++) {
309 53 : std::string num; Tools::convert(i,num);
310 53 : addComponentWithDerivatives("rdc_"+num);
311 53 : componentIsNotPeriodic("rdc_"+num);
312 53 : }
313 66 : for(unsigned i=0; i<ndata; i++) {
314 53 : std::string num; Tools::convert(i,num);
315 53 : addComponent("exp_"+num);
316 53 : componentIsNotPeriodic("exp_"+num);
317 53 : Value* comp=getPntrToComponent("exp_"+num);
318 53 : comp->set(coupl[i]);
319 53 : }
320 : }
321 :
322 29 : if(svd) {
323 1 : addComponent("Sxx"); componentIsNotPeriodic("Sxx");
324 1 : addComponent("Syy"); componentIsNotPeriodic("Syy");
325 1 : addComponent("Szz"); componentIsNotPeriodic("Szz");
326 1 : addComponent("Sxy"); componentIsNotPeriodic("Sxy");
327 1 : addComponent("Sxz"); componentIsNotPeriodic("Sxz");
328 1 : addComponent("Syz"); componentIsNotPeriodic("Syz");
329 : }
330 :
331 29 : requestAtoms(atoms);
332 29 : if(getDoScore()) {
333 12 : setParameters(coupl);
334 12 : Initialise(coupl.size());
335 : }
336 29 : setDerivatives();
337 58 : checkRead();
338 29 : }
339 :
340 5 : void RDC::do_svd()
341 : {
342 : #ifdef __PLUMED_HAS_GSL
343 : gsl_vector *rdc_vec, *S, *Stmp, *work, *bc;
344 : gsl_matrix *coef_mat, *A, *V;
345 5 : rdc_vec = gsl_vector_alloc(coupl.size());
346 5 : bc = gsl_vector_alloc(coupl.size());
347 5 : Stmp = gsl_vector_alloc(5);
348 5 : S = gsl_vector_alloc(5);
349 5 : work = gsl_vector_alloc(5);
350 5 : coef_mat = gsl_matrix_alloc(coupl.size(),5);
351 5 : A = gsl_matrix_alloc(coupl.size(),5);
352 5 : V = gsl_matrix_alloc(5,5);
353 5 : gsl_matrix_set_zero(coef_mat);
354 5 : gsl_vector_set_zero(bc);
355 5 : unsigned index=0;
356 5 : vector<double> dmax(coupl.size());
357 30 : for(unsigned r=0; r<getNumberOfAtoms(); r+=2) {
358 25 : Vector distance;
359 25 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
360 0 : else distance = delta(getPosition(r),getPosition(r+1));
361 25 : double d = distance.modulo();
362 25 : double d2 = d*d;
363 25 : double d3 = d2*d;
364 25 : double id3 = 1./d3;
365 25 : double max = -Const*mu_s*scale;
366 25 : dmax[index] = id3*max;
367 25 : double mu_x = distance[0]/d;
368 25 : double mu_y = distance[1]/d;
369 25 : double mu_z = distance[2]/d;
370 25 : gsl_vector_set(rdc_vec,index,coupl[index]/dmax[index]);
371 25 : gsl_matrix_set(coef_mat,index,0,gsl_matrix_get(coef_mat,index,0)+(mu_x*mu_x-mu_z*mu_z));
372 25 : gsl_matrix_set(coef_mat,index,1,gsl_matrix_get(coef_mat,index,1)+(mu_y*mu_y-mu_z*mu_z));
373 25 : gsl_matrix_set(coef_mat,index,2,gsl_matrix_get(coef_mat,index,2)+(2.0*mu_x*mu_y));
374 25 : gsl_matrix_set(coef_mat,index,3,gsl_matrix_get(coef_mat,index,3)+(2.0*mu_x*mu_z));
375 25 : gsl_matrix_set(coef_mat,index,4,gsl_matrix_get(coef_mat,index,4)+(2.0*mu_y*mu_z));
376 25 : index++;
377 : }
378 5 : gsl_matrix_memcpy(A,coef_mat);
379 5 : gsl_linalg_SV_decomp(A, V, Stmp, work);
380 5 : gsl_linalg_SV_solve(A, V, Stmp, rdc_vec, S);
381 : /* tensor */
382 : Value* tensor;
383 5 : tensor=getPntrToComponent("Sxx");
384 5 : double Sxx = gsl_vector_get(S,0);
385 5 : tensor->set(Sxx);
386 5 : tensor=getPntrToComponent("Syy");
387 5 : double Syy = gsl_vector_get(S,1);
388 5 : tensor->set(Syy);
389 5 : tensor=getPntrToComponent("Szz");
390 5 : double Szz = -Sxx-Syy;
391 5 : tensor->set(Szz);
392 5 : tensor=getPntrToComponent("Sxy");
393 5 : double Sxy = gsl_vector_get(S,2);
394 5 : tensor->set(Sxy);
395 5 : tensor=getPntrToComponent("Sxz");
396 5 : double Sxz = gsl_vector_get(S,3);
397 5 : tensor->set(Sxz);
398 5 : tensor=getPntrToComponent("Syz");
399 5 : double Syz = gsl_vector_get(S,4);
400 5 : tensor->set(Syz);
401 :
402 5 : gsl_blas_dgemv(CblasNoTrans, 1.0, coef_mat, S, 0., bc);
403 30 : for(index=0; index<coupl.size(); index++) {
404 25 : double rdc = gsl_vector_get(bc,index)*dmax[index];
405 25 : Value* val=getPntrToComponent(index);
406 25 : val->set(rdc);
407 : }
408 5 : gsl_matrix_free(coef_mat);
409 5 : gsl_matrix_free(A);
410 5 : gsl_matrix_free(V);
411 5 : gsl_vector_free(rdc_vec);
412 5 : gsl_vector_free(bc);
413 5 : gsl_vector_free(Stmp);
414 5 : gsl_vector_free(S);
415 5 : gsl_vector_free(work);
416 : #endif
417 5 : }
418 :
419 2172 : void RDC::calculate()
420 : {
421 2172 : if(svd) {
422 5 : do_svd();
423 2177 : return;
424 : }
425 :
426 2167 : const double max = -Const*scale*mu_s;
427 2167 : const unsigned N=getNumberOfAtoms();
428 2167 : vector<Vector> dRDC(N/2, Vector{0.,0.,0.});
429 :
430 : /* RDC Calculations and forces */
431 6483 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
432 : {
433 4316 : #pragma omp for
434 : for(unsigned r=0; r<N; r+=2)
435 : {
436 8663 : const unsigned index=r/2;
437 8663 : Vector distance;
438 17288 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
439 0 : else distance = delta(getPosition(r),getPosition(r+1));
440 8640 : const double d2 = distance.modulo2();
441 8650 : const double ind = 1./sqrt(d2);
442 8650 : const double ind2 = 1./d2;
443 8650 : const double ind3 = ind2*ind;
444 8650 : const double x2 = distance[0]*distance[0]*ind2;
445 8667 : const double y2 = distance[1]*distance[1]*ind2;
446 8666 : const double z2 = distance[2]*distance[2]*ind2;
447 8663 : const double dmax = ind3*max;
448 8663 : const double ddmax = dmax*ind2;
449 :
450 8663 : const double rdc = 0.5*dmax*(3.*z2-1.);
451 8663 : const double prod_xy = (x2+y2-4.*z2);
452 8663 : const double prod_z = (3.*x2 + 3.*y2 - 2.*z2);
453 :
454 8663 : dRDC[index] = -1.5*ddmax*distance;
455 8653 : dRDC[index][0] *= prod_xy;
456 8656 : dRDC[index][1] *= prod_xy;
457 8659 : dRDC[index][2] *= prod_z;
458 :
459 8662 : string num; Tools::convert(index,num);
460 8666 : Value* val=getPntrToComponent("rdc_"+num);
461 8668 : val->set(rdc);
462 8663 : if(!getDoScore()) {
463 1355 : setBoxDerivatives(val, Tensor(distance,dRDC[index]));
464 1365 : setAtomsDerivatives(val, r, dRDC[index]);
465 1369 : setAtomsDerivatives(val, r+1, -dRDC[index]);
466 8658 : } else setCalcData(index, rdc);
467 : }
468 : }
469 :
470 2167 : if(getDoScore()) {
471 : /* Metainference */
472 1824 : Tensor dervir;
473 1824 : double score = getScore();
474 1824 : setScore(score);
475 :
476 : /* calculate final derivatives */
477 1824 : Value* val=getPntrToComponent("score");
478 9120 : for(unsigned r=0; r<N; r+=2)
479 : {
480 7296 : const unsigned index=r/2;
481 7296 : Vector distance;
482 7296 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
483 0 : else distance = delta(getPosition(r),getPosition(r+1));
484 7296 : const Vector der = dRDC[index]*getMetaDer(index);
485 7296 : dervir += Tensor(distance, der);
486 7296 : setAtomsDerivatives(val, r, der);
487 7296 : setAtomsDerivatives(val, r+1, -der);
488 : }
489 1824 : setBoxDerivatives(val, dervir);
490 2167 : }
491 : }
492 :
493 327 : void RDC::update() {
494 : // write status file
495 327 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
496 327 : }
497 :
498 : }
499 4821 : }
|
