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