Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 :
23 : #include "Colvar.h"
24 : #include "ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 :
27 : #ifdef __PLUMED_HAS_GSL
28 : #include <gsl/gsl_vector.h>
29 : #include <gsl/gsl_matrix.h>
30 : #include <gsl/gsl_linalg.h>
31 : #include <gsl/gsl_blas.h>
32 : #endif
33 :
34 : using namespace std;
35 :
36 : namespace PLMD {
37 : namespace colvar {
38 :
39 : //+PLUMEDOC COLVAR RDC
40 : /*
41 : Calculates the (Residual) Dipolar Coupling between two atoms.
42 :
43 : The RDC between two atomic nuclei depends on the \f$\theta\f$ angle between
44 : the inter-nuclear vector and the external magnetic field. In isotropic media RDCs average to zero because of the orientational
45 : averaging, but when the rotational symmetry is broken, either through the introduction of an alignment medium or for molecules
46 : with highly anisotropic paramagnetic susceptibility, RDCs become measurable.
47 :
48 : \f[
49 : D=D_{max}0.5(3\cos^2(\theta)-1)
50 : \f]
51 :
52 : where
53 :
54 : \f[
55 : D_{max}=-\mu_0\gamma_1\gamma_2h/(8\pi^3r^3)
56 : \f]
57 :
58 : that is the maximal value of the dipolar coupling for the two nuclear spins with gyromagnetic ratio \f$\gamma\f$.
59 : \f$\mu\f$ is the magnetic constant and h is the Planck constant.
60 :
61 : Common Gyromagnetic Ratios (C.G.S)
62 : - H(1) 26.7513
63 : - C(13) 6.7261
64 : - N(15) -2.7116
65 : - NH -72.5388
66 : - CH 179.9319
67 : - CN -18.2385
68 : - CC 45.2404
69 :
70 : This collective variable calculates the Residual Dipolar Coupling for a set of couple of atoms using the above definition.
71 : From the calculated RDCs and a set of experimental values it calculates either their correlation or the squared quality factor
72 :
73 : \f[
74 : Q^2=\frac{\sum_i(D_i-D^{exp}_i)^2}{\sum_i(D^{exp}_i)^2}
75 : \f]
76 :
77 : RDCs report only on the fraction of molecules that is aligned, this means that comparing the RDCs from a single structure in
78 : a MD simulation to the experimental values is not particularly meaningfull, from this point of view it is better to compare
79 : their correlation. The fraction of aligned molecules can be obtained by maximising the correlation between the calculated and
80 : the experimental RDCs. This fraction can be used as a scaling factor in the calculation of the RDCs in order to compare their
81 : values. The averaging of the RDCs calculated with the above definition from a standard MD should result to 0 because of
82 : the rotational diffusion, but this variable can be used to break the rotational symmetry.
83 :
84 : RDCs can also be calculated using a Single Value Decomposition approach, in this case the code rely on the
85 : a set of function from the GNU Scientific Library (GSL). (With SVD forces are not currently implemented).
86 :
87 : Replica-Averaged restrained simulations can be performed with this CV and the function \ref ENSEMBLE.
88 :
89 : Additional material and examples can be also found in the tutorial \ref belfast-9
90 :
91 : \par Examples
92 :
93 : In the following example five N-H RDCs are defined and their correlation with
94 : respect to a set of experimental data is calculated and restrained. In addition,
95 : and only for analysis purposes, the same RDCs are calculated using a Single Value
96 : Decomposition algorithm.
97 :
98 : \verbatim
99 : RDC ...
100 : GYROM=-72.5388
101 : SCALE=1.0
102 : ATOMS1=20,21
103 : ATOMS2=37,38
104 : ATOMS3=56,57
105 : ATOMS4=76,77
106 : ATOMS5=92,93
107 : LABEL=nh
108 : ... RDC
109 :
110 : STATS ARG=nh.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
111 :
112 : rdce: RESTRAINT ARG=nh.corr KAPPA=0. SLOPE=-25000.0 AT=1.
113 :
114 : RDC ...
115 : GYROM=-72.5388
116 : SCALE=1.0
117 : SVD
118 : ATOMS1=20,21 COUPLING1=8.17
119 : ATOMS2=37,38 COUPLING2=-8.271
120 : ATOMS3=56,57 COUPLING3=-10.489
121 : ATOMS4=76,77 COUPLING4=-9.871
122 : ATOMS5=92,93 COUPLING5=-9.152
123 : LABEL=svd
124 : ... RDC
125 :
126 : PRINT ARG=nh.corr,rdce.bias FILE=colvar
127 : PRINT ARG=svd.* FILE=svd
128 : \endverbatim
129 : (See also \ref PRINT, \ref RESTRAINT)
130 :
131 : */
132 : //+ENDPLUMEDOC
133 :
134 26 : class RDC : public Colvar {
135 : private:
136 : const double Const;
137 : double mu_s;
138 : double scale;
139 : vector<double> coupl;
140 : bool svd;
141 : bool pbc;
142 : public:
143 : explicit RDC(const ActionOptions&);
144 : static void registerKeywords( Keywords& keys );
145 : virtual void calculate();
146 : };
147 :
148 2536 : PLUMED_REGISTER_ACTION(RDC,"RDC")
149 :
150 14 : void RDC::registerKeywords( Keywords& keys ) {
151 14 : Colvar::registerKeywords( keys );
152 14 : componentsAreNotOptional(keys);
153 14 : useCustomisableComponents(keys);
154 : keys.add("numbered","ATOMS","the couple of atoms involved in each of the bonds for which you wish to calculate the RDC. "
155 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dipolar coupling will be "
156 14 : "calculated for each ATOMS keyword you specify.");
157 14 : keys.reset_style("ATOMS","atoms");
158 14 : keys.add("compulsory","GYROM","Add the product of the gyromagnetic constants for the bond. ");
159 14 : keys.add("compulsory","SCALE","Add the scaling factor to take into account concentration and other effects. ");
160 14 : keys.addFlag("SVD",false,"Set to TRUE if you want to backcalculate using Single Value Decomposition (need GSL at compilation time).");
161 14 : keys.addFlag("ADDCOUPLINGS",false,"Set to TRUE if you want to have fixed components with the experimetnal values.");
162 14 : keys.add("numbered","COUPLING","Add an experimental value for each coupling (needed by SVD and usefull for \ref STATS).");
163 14 : keys.addOutputComponent("rdc","default","the calculated # RDC");
164 14 : keys.addOutputComponent("exp","SVD/ADDCOUPLINGS","the experimental # RDC");
165 14 : }
166 :
167 13 : RDC::RDC(const ActionOptions&ao):
168 : PLUMED_COLVAR_INIT(ao),
169 : Const(0.3356806),
170 : mu_s(0),
171 : scale(1),
172 13 : pbc(true)
173 : {
174 13 : bool nopbc=!pbc;
175 13 : parseFlag("NOPBC",nopbc);
176 13 : pbc=!nopbc;
177 :
178 : // Read in the atoms
179 26 : vector<AtomNumber> t, atoms;
180 66 : for(int i=1;; ++i ) {
181 66 : parseAtomList("ATOMS", i, t );
182 66 : if( t.empty() ) break;
183 53 : if( t.size()!=2 ) {
184 0 : std::string ss; Tools::convert(i,ss);
185 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
186 : }
187 53 : atoms.push_back(t[0]);
188 53 : atoms.push_back(t[1]);
189 53 : t.resize(0);
190 53 : }
191 :
192 13 : const unsigned ndata = atoms.size()/2;
193 :
194 : // Read in GYROMAGNETIC constants
195 13 : parse("GYROM", mu_s);
196 13 : if(mu_s==0.) error("GYROM must be set");
197 :
198 : // Read in SCALING factors
199 13 : parse("SCALE", scale);
200 13 : if(scale==0.) error("SCALE must be different from 0");
201 :
202 13 : svd=false;
203 13 : parseFlag("SVD",svd);
204 : #ifndef __PLUMED_HAS_GSL
205 : if(svd) error("You CANNOT use SVD without GSL. Recompile PLUMED with GSL!\n");
206 : #endif
207 :
208 13 : bool addcoupling=false;
209 13 : parseFlag("ADDCOUPLINGS",addcoupling);
210 :
211 13 : if(svd||addcoupling) {
212 1 : coupl.resize( ndata );
213 1 : unsigned ntarget=0;
214 6 : for(unsigned i=0; i<ndata; ++i) {
215 5 : if( !parseNumbered( "COUPLING", i+1, coupl[i] ) ) break;
216 5 : ntarget++;
217 : }
218 1 : if( ntarget!=ndata ) error("found wrong number of COUPLING values");
219 : }
220 :
221 : // Ouput details of all contacts
222 13 : log.printf(" Gyromagnetic moment is %f. Scaling factor is %f.",mu_s,scale);
223 66 : for(unsigned i=0; i<ndata; ++i) {
224 53 : log.printf(" The %uth Bond Dipolar Coupling is calculated from atoms : %d %d.", i+1, atoms[2*i].serial(), atoms[2*i+1].serial());
225 53 : if(svd||addcoupling) log.printf(" Experimental coupling is %f.", coupl[i]);
226 53 : log.printf("\n");
227 : }
228 :
229 13 : log<<" Bibliography "
230 39 : <<plumed.cite("Camilloni C, Vendruscolo M, J. Phys. Chem. B 119, 653 (2015)")
231 39 : <<plumed.cite("Camilloni C, Vendruscolo M, Biochemistry 54, 7470 (2015)") <<"\n";
232 :
233 13 : checkRead();
234 :
235 66 : for(unsigned i=0; i<ndata; i++) {
236 53 : std::string num; Tools::convert(i,num);
237 53 : if(!svd) {
238 48 : addComponentWithDerivatives("rdc_"+num);
239 48 : componentIsNotPeriodic("rdc_"+num);
240 : } else {
241 5 : addComponent("rdc_"+num);
242 5 : componentIsNotPeriodic("rdc_"+num);
243 : }
244 53 : }
245 :
246 13 : if(svd||addcoupling) {
247 6 : for(unsigned i=0; i<ndata; i++) {
248 5 : std::string num; Tools::convert(i,num);
249 5 : addComponent("exp_"+num);
250 5 : componentIsNotPeriodic("exp_"+num);
251 5 : Value* comp=getPntrToComponent("exp_"+num); comp->set(coupl[i]);
252 5 : }
253 : }
254 :
255 26 : requestAtoms(atoms);
256 13 : }
257 :
258 300 : void RDC::calculate()
259 : {
260 300 : if(!svd) {
261 295 : const double max = -Const*scale*mu_s;
262 295 : const unsigned N=getNumberOfAtoms();
263 : /* RDC Calculations and forces */
264 1475 : for(unsigned r=0; r<N; r+=2)
265 : {
266 1180 : const unsigned index=r/2;
267 1180 : Vector distance;
268 1180 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
269 0 : else distance = delta(getPosition(r),getPosition(r+1));
270 1180 : const double d = distance.modulo();
271 1180 : const double ind = 1./d;
272 1180 : const double id3 = ind*ind*ind;
273 1180 : const double dmax = id3*max;
274 1180 : const double cos_theta = distance[2]*ind;
275 :
276 1180 : const double rdc = 0.5*dmax*(3.*cos_theta*cos_theta-1.);
277 :
278 1180 : const double id7 = id3*id3*ind;
279 1180 : const double id9 = id7*ind*ind;
280 1180 : const double x2=distance[0]*distance[0];
281 1180 : const double y2=distance[1]*distance[1];
282 1180 : const double z2=distance[2]*distance[2];
283 1180 : const double prod = -max*id7*(1.5*x2 +1.5*y2 -6.*z2);
284 :
285 1180 : Vector dRDC;
286 1180 : dRDC[0] = prod*distance[0];
287 1180 : dRDC[1] = prod*distance[1];
288 1180 : dRDC[2] = -max*id9*distance[2]*(4.5*x2*x2 + 4.5*y2*y2 + 1.5*y2*z2 - 3.*z2*z2 + x2*(9.*y2 + 1.5*z2));
289 :
290 1180 : Value* val=getPntrToComponent(index);
291 1180 : val->set(rdc);
292 1180 : setBoxDerivatives(val, Tensor(distance,dRDC));
293 1180 : setAtomsDerivatives(val, r, dRDC);
294 1180 : setAtomsDerivatives(val, r+1, -dRDC);
295 : }
296 :
297 : } else {
298 :
299 : #ifdef __PLUMED_HAS_GSL
300 : gsl_vector *rdc_vec, *S, *Stmp, *work, *bc;
301 : gsl_matrix *coef_mat, *A, *V;
302 5 : rdc_vec = gsl_vector_alloc(coupl.size());
303 5 : bc = gsl_vector_alloc(coupl.size());
304 5 : Stmp = gsl_vector_alloc(5);
305 5 : S = gsl_vector_alloc(5);
306 5 : work = gsl_vector_alloc(5);
307 5 : coef_mat = gsl_matrix_alloc(coupl.size(),5);
308 5 : A = gsl_matrix_alloc(coupl.size(),5);
309 5 : V = gsl_matrix_alloc(5,5);
310 5 : gsl_matrix_set_zero(coef_mat);
311 5 : gsl_vector_set_zero(bc);
312 5 : unsigned index=0;
313 5 : vector<double> dmax(coupl.size());
314 30 : for(unsigned r=0; r<getNumberOfAtoms(); r+=2) {
315 25 : Vector distance;
316 25 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
317 0 : else distance = delta(getPosition(r),getPosition(r+1));
318 25 : double d = distance.modulo();
319 25 : double d2 = d*d;
320 25 : double d3 = d2*d;
321 25 : double id3 = 1./d3;
322 25 : double max = -Const*mu_s*scale;
323 25 : dmax[index] = id3*max;
324 25 : double mu_x = distance[0]/d;
325 25 : double mu_y = distance[1]/d;
326 25 : double mu_z = distance[2]/d;
327 25 : gsl_vector_set(rdc_vec,index,coupl[index]/dmax[index]);
328 25 : gsl_matrix_set(coef_mat,index,0,gsl_matrix_get(coef_mat,index,0)+(mu_x*mu_x-mu_z*mu_z));
329 25 : gsl_matrix_set(coef_mat,index,1,gsl_matrix_get(coef_mat,index,1)+(mu_y*mu_y-mu_z*mu_z));
330 25 : gsl_matrix_set(coef_mat,index,2,gsl_matrix_get(coef_mat,index,2)+(2.0*mu_x*mu_y));
331 25 : gsl_matrix_set(coef_mat,index,3,gsl_matrix_get(coef_mat,index,3)+(2.0*mu_x*mu_z));
332 25 : gsl_matrix_set(coef_mat,index,4,gsl_matrix_get(coef_mat,index,4)+(2.0*mu_y*mu_z));
333 25 : index++;
334 : }
335 5 : gsl_matrix_memcpy(A,coef_mat);
336 5 : gsl_linalg_SV_decomp(A, V, Stmp, work);
337 5 : gsl_linalg_SV_solve(A, V, Stmp, rdc_vec, S);
338 : /* tensor
339 : double Sxx = gsl_vector_get(S,0);
340 : double Syy = gsl_vector_get(S,1);
341 : double Szz = -Sxx-Syy;
342 : double Sxy = gsl_vector_get(S,2);
343 : double Sxz = gsl_vector_get(S,3);
344 : double Syz = gsl_vector_get(S,4);
345 : */
346 5 : gsl_blas_dgemv(CblasNoTrans, 1.0, coef_mat, S, 0., bc);
347 30 : for(index=0; index<coupl.size(); index++) {
348 25 : double rdc = gsl_vector_get(bc,index)*dmax[index];
349 25 : Value* val=getPntrToComponent(index);
350 25 : val->set(rdc);
351 : }
352 5 : gsl_matrix_free(coef_mat);
353 5 : gsl_matrix_free(A);
354 5 : gsl_matrix_free(V);
355 5 : gsl_vector_free(rdc_vec);
356 5 : gsl_vector_free(bc);
357 5 : gsl_vector_free(Stmp);
358 5 : gsl_vector_free(S);
359 5 : gsl_vector_free(work);
360 : #endif
361 :
362 : }
363 :
364 300 : }
365 :
366 : }
367 2523 : }
|
