Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 :
26 : namespace PLMD {
27 : namespace isdb {
28 :
29 : //+PLUMEDOC ISDB_COLVAR RDC
30 : /*
31 : Calculates the (Residual) Dipolar Coupling between two atoms.
32 :
33 : The Dipolar Coupling between two nuclei depends on the $\theta$ angle between
34 : the inter-nuclear vector and the external magnetic field.
35 :
36 : $$
37 : D=D_{max}0.5(3\cos^2(\theta)-1)
38 : $$
39 :
40 : where
41 :
42 : $$
43 : D_{max}=-\mu_0\gamma_1\gamma_2h/(8\pi^3r^3)
44 : $$
45 :
46 : that is the maximal value of the dipolar coupling for the two nuclear spins with gyromagnetic ratio $\gamma$.
47 : $\mu$ is the magnetic constant and h is the Planck constant.
48 :
49 : Common Gyromagnetic Ratios (C.G.S)
50 : - H(1) 26.7513
51 : - C(13) 6.7261
52 : - N(15) -2.7116
53 : and their products (this is what is given in input using the keyword GYROM)
54 : - N-H -72.5388
55 : - C-H 179.9319
56 : - C-N -18.2385
57 : - C-C 45.2404
58 :
59 : In isotropic media DCs average to zero because of the rotational
60 : averaging, but when the rotational symmetry is broken, either through the introduction of an alignment medium or for molecules
61 : 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).
62 :
63 : This collective variable calculates the Dipolar Coupling for a set of couple of atoms using
64 : the above definition.
65 :
66 : 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".
67 :
68 : 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 as discussed in the first of the papers cited below. 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.
69 :
70 : Replica-Averaged simulations can be performed using RDCs, [ENSEMBLE](ENSEMBLE.md), [STATS](STATS.md) and [RESTRAINT](RESTRAINT.md).
71 : [METAINFERENCE](METAINFERENCE.md) can be activated using DOSCORE and the other relevant keywords.
72 :
73 : Additional material and examples can be also found in the tutorials
74 :
75 : ## Examples
76 :
77 : In the following example five N-H RDCs are defined and averaged over multiple replicas,
78 : their correlation is then calculated with respect to a set of experimental data and restrained.
79 : In addition, and only for analysis purposes, the same RDCs each single conformation are calculated
80 : using a Single Value Decomposition algorithm, then averaged and again compared with the experimental data.
81 :
82 : ```plumed
83 : #SETTINGS NREPLICAS=2
84 : RDC ...
85 : GYROM=-72.5388
86 : SCALE=0.001
87 : ATOMS1=20,21
88 : ATOMS2=37,38
89 : ATOMS3=56,57
90 : ATOMS4=76,77
91 : ATOMS5=92,93
92 : LABEL=nh
93 : ... RDC
94 :
95 : erdc: ENSEMBLE ARG=nh.*
96 :
97 : st: STATS ARG=erdc.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
98 :
99 : rdce: RESTRAINT ARG=st.corr KAPPA=0. SLOPE=-25000.0 AT=1.
100 :
101 : PRINT ARG=st.corr,rdce.bias FILE=colvar
102 : ```
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 : //+PLUMEDOC ISDB_COLVAR PCS
108 : /*
109 : Calculates the Pseudo-contact shift of a nucleus determined by the presence of a metal ion susceptible to anisotropic magnetization.
110 :
111 : The PCS of an atomic nucleus depends on the $\theta$ angle between the vector from the spin-label to the nucleus
112 : and the external magnetic field and the module of the vector itself the second of the two papers cited below. While in principle the averaging
113 : resulting from the tumbling should remove the pseudo-contact shift, in presence of the NMR magnetic field the magnetically anisotropic molecule bound to system will break the rotational symmetry does resulting in measurable values for the PCS and RDC.
114 :
115 : Replica-Averaged simulations can be performed using PCS values, [ENSEMBLE](ENSEMBLE.md), [STATS](STATS.md) and [RESTRAINT](RESTRAINT.md) .
116 : Metainference simulations can be performed with this CV and [METAINFERENCE](METAINFERENCE.md).
117 :
118 : ## Examples
119 :
120 : In the following example five PCS values are defined and their correlation with
121 : respect to a set of experimental data is calculated and restrained. In addition,
122 : and only for analysis purposes, the same PCS values are calculated using a Single Value
123 : Decomposition algorithm.
124 :
125 : ```plumed
126 : PCS ...
127 : ATOMS1=20,21
128 : ATOMS2=20,38
129 : ATOMS3=20,57
130 : ATOMS4=20,77
131 : ATOMS5=20,93
132 : LABEL=nh
133 : ... PCS
134 :
135 : enh: ENSEMBLE ARG=nh.*
136 :
137 : st: STATS ARG=enh.* PARAMETERS=8.17,-8.271,-10.489,-9.871,-9.152
138 :
139 : pcse: RESTRAINT ARG=st.corr KAPPA=0. SLOPE=-25000.0 AT=1.
140 :
141 : PRINT ARG=st.corr,pcse.bias FILE=colvar
142 : ```
143 :
144 : */
145 : //+ENDPLUMEDOC
146 :
147 : class RDC :
148 : public MetainferenceBase {
149 : private:
150 : // the components with the back-calculated data
151 : std::vector<Value*> dataValues_;
152 : double Const;
153 : double mu_s;
154 : double scale;
155 : std::vector<double> coupl;
156 : bool pbc;
157 :
158 : public:
159 : explicit RDC(const ActionOptions&);
160 : static void registerKeywords( Keywords& keys );
161 : void calculate() override;
162 : };
163 :
164 : PLUMED_REGISTER_ACTION(RDC,"RDC")
165 : PLUMED_REGISTER_ACTION(RDC,"PCS")
166 :
167 110 : void RDC::registerKeywords( Keywords& keys ) {
168 110 : MetainferenceBase::registerKeywords(keys);
169 110 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
170 110 : keys.add("numbered","ATOMS","the couple of atoms involved in each of the bonds for which you wish to calculate the RDC. "
171 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dipolar coupling will be "
172 : "calculated for each ATOMS keyword you specify.");
173 220 : keys.reset_style("ATOMS","atoms");
174 110 : keys.add("compulsory","GYROM","1.","Add the product of the gyromagnetic constants for the bond. ");
175 110 : keys.add("compulsory","SCALE","1.","Add the scaling factor to take into account concentration and other effects. ");
176 110 : keys.add("numbered","COUPLING","Add an experimental value for each coupling (useful for STATS).");
177 220 : keys.addOutputComponent("rdc","default","scalar","the calculated # RDC");
178 220 : keys.addOutputComponent("exp","COUPLING","scalar","the experimental # RDC");
179 110 : keys.addDOI("10.1021/jp5021824");
180 110 : keys.addDOI("10.1021/acs.biochem.5b01138");
181 110 : }
182 :
183 102 : RDC::RDC(const ActionOptions&ao):
184 : PLUMED_METAINF_INIT(ao),
185 102 : Const(1.),
186 102 : mu_s(1.),
187 102 : scale(1.),
188 102 : pbc(true) {
189 102 : bool nopbc=!pbc;
190 102 : parseFlag("NOPBC",nopbc);
191 102 : pbc=!nopbc;
192 :
193 : const double RDCConst = 0.3356806;
194 : const double PCSConst = 1.0;
195 :
196 102 : if( getName().find("RDC")!=std::string::npos) {
197 102 : Const *= RDCConst;
198 0 : } else if( getName().find("PCS")!=std::string::npos) {
199 : Const *= PCSConst;
200 : }
201 :
202 : // Read in the atoms
203 102 : std::vector<AtomNumber> t, atoms;
204 102 : for(int i=1;; ++i ) {
205 1784 : parseAtomList("ATOMS", i, t );
206 892 : if( t.empty() ) {
207 : break;
208 : }
209 790 : if( t.size()!=2 ) {
210 : std::string ss;
211 0 : Tools::convert(i,ss);
212 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
213 : }
214 790 : atoms.push_back(t[0]);
215 790 : atoms.push_back(t[1]);
216 790 : t.resize(0);
217 790 : }
218 :
219 102 : const unsigned ndata = atoms.size()/2;
220 :
221 : // Read in GYROMAGNETIC constants
222 102 : parse("GYROM", mu_s);
223 102 : if(mu_s==0.) {
224 0 : error("GYROM cannot be 0");
225 : }
226 :
227 : // Read in SCALING factors
228 102 : parse("SCALE", scale);
229 102 : if(scale==0.) {
230 0 : error("SCALE cannot be 0");
231 : }
232 :
233 : // Optionally add an experimental value
234 102 : coupl.resize( ndata );
235 : unsigned ntarget=0;
236 457 : for(unsigned i=0; i<ndata; ++i) {
237 832 : if( !parseNumbered( "COUPLING", i+1, coupl[i] ) ) {
238 : break;
239 : }
240 355 : ntarget++;
241 : }
242 : bool addexp=false;
243 102 : if(ntarget!=ndata && ntarget!=0) {
244 0 : error("found wrong number of COUPLING values");
245 : }
246 102 : if(ntarget==ndata) {
247 : addexp=true;
248 : }
249 102 : if(getDoScore()&&!addexp) {
250 0 : error("with DOSCORE you need to set the COUPLING values");
251 : }
252 :
253 : // Output details of all contacts
254 102 : log.printf(" Gyromagnetic moment is %f. Scaling factor is %f.",mu_s,scale);
255 892 : for(unsigned i=0; i<ndata; ++i) {
256 790 : log.printf(" The %uth Bond Dipolar Coupling is calculated from atoms : %d %d.", i+1, atoms[2*i].serial(), atoms[2*i+1].serial());
257 790 : if(addexp) {
258 355 : log.printf(" Experimental coupling is %f.", coupl[i]);
259 : }
260 790 : log.printf("\n");
261 : }
262 :
263 102 : log<<" Bibliography "
264 102 : <<plumed.cite("Camilloni C, Vendruscolo M, J. Phys. Chem. B 119, 653 (2015)")
265 306 : <<plumed.cite("Camilloni C, Vendruscolo M, Biochemistry 54, 7470 (2015)");
266 204 : log<< plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
267 :
268 :
269 102 : if(!getDoScore()) {
270 496 : for(unsigned i=0; i<ndata; i++) {
271 : std::string num;
272 435 : Tools::convert(i,num);
273 870 : addComponentWithDerivatives("rdc-"+num);
274 435 : componentIsNotPeriodic("rdc-"+num);
275 435 : dataValues_.push_back(getPntrToComponent(getNumberOfComponents()-1));
276 : }
277 61 : if(addexp) {
278 0 : for(unsigned i=0; i<ndata; i++) {
279 : std::string num;
280 0 : Tools::convert(i,num);
281 0 : addComponent("exp-"+num);
282 0 : componentIsNotPeriodic("exp-"+num);
283 0 : Value* comp=getPntrToComponent("exp-"+num);
284 0 : comp->set(coupl[i]);
285 : }
286 : }
287 : } else {
288 396 : for(unsigned i=0; i<ndata; i++) {
289 : std::string num;
290 355 : Tools::convert(i,num);
291 710 : addComponent("rdc-"+num);
292 355 : componentIsNotPeriodic("rdc-"+num);
293 355 : dataValues_.push_back(getPntrToComponent(getNumberOfComponents()-1));
294 : }
295 396 : for(unsigned i=0; i<ndata; i++) {
296 : std::string num;
297 355 : Tools::convert(i,num);
298 710 : addComponent("exp-"+num);
299 355 : componentIsNotPeriodic("exp-"+num);
300 355 : Value* comp=getPntrToComponent("exp-"+num);
301 355 : comp->set(coupl[i]);
302 : }
303 : }
304 :
305 102 : requestAtoms(atoms, false);
306 102 : if(getDoScore()) {
307 41 : setParameters(coupl);
308 41 : Initialise(coupl.size());
309 : }
310 102 : setDerivatives();
311 102 : checkRead();
312 102 : }
313 :
314 2347 : void RDC::calculate() {
315 2347 : const double max = -Const*scale*mu_s;
316 : const unsigned N=getNumberOfAtoms();
317 2347 : std::vector<Vector> dRDC(N/2, Vector{0.,0.,0.});
318 :
319 : /* RDC Calculations and forces */
320 2347 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
321 : {
322 : #pragma omp for
323 : for(unsigned r=0; r<N; r+=2) {
324 : const unsigned index=r/2;
325 : Vector distance;
326 : if(pbc) {
327 : distance = pbcDistance(getPosition(r),getPosition(r+1));
328 : } else {
329 : distance = delta(getPosition(r),getPosition(r+1));
330 : }
331 : const double d2 = distance.modulo2();
332 : const double ind = 1./std::sqrt(d2);
333 : const double ind2 = 1./d2;
334 : const double ind3 = ind2*ind;
335 : const double x2 = distance[0]*distance[0]*ind2;
336 : const double y2 = distance[1]*distance[1]*ind2;
337 : const double z2 = distance[2]*distance[2]*ind2;
338 : const double dmax = ind3*max;
339 : const double ddmax = dmax*ind2;
340 :
341 : const double rdc = 0.5*dmax*(3.*z2-1.);
342 : const double prod_xy = (x2+y2-4.*z2);
343 : const double prod_z = (3.*x2 + 3.*y2 - 2.*z2);
344 :
345 : dRDC[index] = -1.5*ddmax*distance;
346 : dRDC[index][0] *= prod_xy;
347 : dRDC[index][1] *= prod_xy;
348 : dRDC[index][2] *= prod_z;
349 :
350 : Value* val=dataValues_[index];
351 : val->set(rdc);
352 : if(!getDoScore()) {
353 : setBoxDerivatives(val, Tensor(distance,dRDC[index]));
354 : setAtomsDerivatives(val, r, dRDC[index]);
355 : setAtomsDerivatives(val, r+1, -dRDC[index]);
356 : } else {
357 : setCalcData(index, rdc);
358 : }
359 : }
360 : }
361 :
362 2347 : if(getDoScore()) {
363 : /* Metainference */
364 1896 : Tensor dervir;
365 1896 : double score = getScore();
366 1896 : setScore(score);
367 :
368 : /* calculate final derivatives */
369 1896 : Value* val=getPntrToComponent("score");
370 9762 : for(unsigned r=0; r<N; r+=2) {
371 7866 : const unsigned index=r/2;
372 7866 : Vector distance;
373 7866 : if(pbc) {
374 7848 : distance = pbcDistance(getPosition(r),getPosition(r+1));
375 : } else {
376 18 : distance = delta(getPosition(r),getPosition(r+1));
377 : }
378 : const Vector der = dRDC[index]*getMetaDer(index);
379 7866 : dervir += Tensor(distance, der);
380 7866 : setAtomsDerivatives(val, r, der);
381 7866 : setAtomsDerivatives(val, r+1, -der);
382 : }
383 1896 : setBoxDerivatives(val, dervir);
384 : }
385 10213 : }
386 :
387 : }
388 : }
|
