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 performed 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 experimental 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 Pseudo-contact 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 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.
140 :
141 : PCS values 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 performed using PCS values, \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 PCS values 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 PCS values 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 87 : 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 : #ifdef __PLUMED_HAS_GSL
188 : /// Auxiliary class to delete a gsl_vector.
189 : /// If used somewhere else we can move it.
190 : struct gsl_vector_deleter {
191 : void operator()(gsl_vector* p) {
192 25 : gsl_vector_free(p);
193 : }
194 : };
195 :
196 : /// unique_ptr to a gsl_vector.
197 : /// Gets deleted when going out of scope.
198 : typedef std::unique_ptr<gsl_vector,gsl_vector_deleter> gsl_vector_unique_ptr;
199 :
200 : /// Auxiliary class to delete a gsl_matrix.
201 : /// If used somewhere else we can move it.
202 : struct gsl_matrix_deleter {
203 : void operator()(gsl_matrix* p) {
204 15 : gsl_matrix_free(p);
205 : }
206 : };
207 :
208 : /// unique_ptr to a gsl_matrix.
209 : /// Gets deleted when going out of scope.
210 : typedef std::unique_ptr<gsl_matrix,gsl_matrix_deleter> gsl_matrix_unique_ptr;
211 : #endif
212 :
213 :
214 : void do_svd();
215 : public:
216 : explicit RDC(const ActionOptions&);
217 : static void registerKeywords( Keywords& keys );
218 : void calculate() override;
219 : void update() override;
220 : };
221 :
222 7890 : PLUMED_REGISTER_ACTION(RDC,"RDC")
223 7832 : PLUMED_REGISTER_ACTION(RDC,"PCS")
224 :
225 31 : void RDC::registerKeywords( Keywords& keys ) {
226 31 : componentsAreNotOptional(keys);
227 31 : MetainferenceBase::registerKeywords(keys);
228 93 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
229 : keys.add("numbered","ATOMS","the couple of atoms involved in each of the bonds for which you wish to calculate the RDC. "
230 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dipolar coupling will be "
231 124 : "calculated for each ATOMS keyword you specify.");
232 93 : keys.reset_style("ATOMS","atoms");
233 155 : keys.add("compulsory","GYROM","1.","Add the product of the gyromagnetic constants for the bond. ");
234 155 : keys.add("compulsory","SCALE","1.","Add the scaling factor to take into account concentration and other effects. ");
235 93 : keys.addFlag("SVD",false,"Set to TRUE if you want to back calculate using Single Value Decomposition (need GSL at compilation time).");
236 124 : keys.add("numbered","COUPLING","Add an experimental value for each coupling (needed by SVD and useful for \\ref STATS).");
237 124 : keys.addOutputComponent("rdc","default","the calculated # RDC");
238 124 : keys.addOutputComponent("exp","SVD/COUPLING","the experimental # RDC");
239 124 : keys.addOutputComponent("Sxx","SVD","Tensor component");
240 124 : keys.addOutputComponent("Syy","SVD","Tensor component");
241 124 : keys.addOutputComponent("Szz","SVD","Tensor component");
242 124 : keys.addOutputComponent("Sxy","SVD","Tensor component");
243 124 : keys.addOutputComponent("Sxz","SVD","Tensor component");
244 124 : keys.addOutputComponent("Syz","SVD","Tensor component");
245 31 : }
246 :
247 29 : RDC::RDC(const ActionOptions&ao):
248 : PLUMED_METAINF_INIT(ao),
249 : Const(1.),
250 : mu_s(1.),
251 : scale(1.),
252 58 : pbc(true)
253 : {
254 29 : bool nopbc=!pbc;
255 58 : parseFlag("NOPBC",nopbc);
256 29 : pbc=!nopbc;
257 :
258 : const double RDCConst = 0.3356806;
259 : const double PCSConst = 1.0;
260 :
261 29 : if( getName().find("RDC")!=std::string::npos) { Const *= RDCConst; }
262 0 : else if( getName().find("PCS")!=std::string::npos) { Const *= PCSConst; }
263 :
264 : // Read in the atoms
265 : vector<AtomNumber> t, atoms;
266 117 : for(int i=1;; ++i ) {
267 292 : parseAtomList("ATOMS", i, t );
268 146 : if( t.empty() ) break;
269 117 : if( t.size()!=2 ) {
270 0 : std::string ss; Tools::convert(i,ss);
271 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
272 : }
273 117 : atoms.push_back(t[0]);
274 117 : atoms.push_back(t[1]);
275 117 : t.resize(0);
276 117 : }
277 :
278 29 : const unsigned ndata = atoms.size()/2;
279 :
280 : // Read in GYROMAGNETIC constants
281 58 : parse("GYROM", mu_s);
282 29 : if(mu_s==0.) error("GYROM cannot be 0");
283 :
284 : // Read in SCALING factors
285 58 : parse("SCALE", scale);
286 29 : if(scale==0.) error("SCALE cannot be 0");
287 :
288 29 : svd=false;
289 58 : parseFlag("SVD",svd);
290 : #ifndef __PLUMED_HAS_GSL
291 : if(svd) error("You CANNOT use SVD without GSL. Recompile PLUMED with GSL!\n");
292 : #endif
293 30 : if(svd&&getDoScore()) error("It is not possible to use SVD and METAINFERENCE together");
294 :
295 : // Optionally add an experimental value
296 29 : coupl.resize( ndata );
297 : unsigned ntarget=0;
298 82 : for(unsigned i=0; i<ndata; ++i) {
299 207 : if( !parseNumbered( "COUPLING", i+1, coupl[i] ) ) break;
300 53 : ntarget++;
301 : }
302 : bool addexp=false;
303 29 : if(ntarget!=ndata && ntarget!=0) error("found wrong number of COUPLING values");
304 29 : if(ntarget==ndata) addexp=true;
305 29 : if(getDoScore()&&!addexp) error("with DOSCORE you need to set the COUPLING values");
306 29 : if(svd&&!addexp) error("with SVD you need to set the COUPLING values");
307 :
308 :
309 : // Ouput details of all contacts
310 29 : log.printf(" Gyromagnetic moment is %f. Scaling factor is %f.",mu_s,scale);
311 146 : for(unsigned i=0; i<ndata; ++i) {
312 351 : log.printf(" The %uth Bond Dipolar Coupling is calculated from atoms : %d %d.", i+1, atoms[2*i].serial(), atoms[2*i+1].serial());
313 170 : if(addexp) log.printf(" Experimental coupling is %f.", coupl[i]);
314 117 : log.printf("\n");
315 : }
316 :
317 29 : log<<" Bibliography "
318 116 : <<plumed.cite("Camilloni C, Vendruscolo M, J. Phys. Chem. B 119, 653 (2015)")
319 116 : <<plumed.cite("Camilloni C, Vendruscolo M, Biochemistry 54, 7470 (2015)");
320 87 : log<< plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
321 :
322 :
323 29 : if(!getDoScore()&&!svd) {
324 64 : for(unsigned i=0; i<ndata; i++) {
325 64 : std::string num; Tools::convert(i,num);
326 128 : addComponentWithDerivatives("rdc-"+num);
327 128 : componentIsNotPeriodic("rdc-"+num);
328 : }
329 16 : if(addexp) {
330 0 : for(unsigned i=0; i<ndata; i++) {
331 0 : std::string num; Tools::convert(i,num);
332 0 : addComponent("exp-"+num);
333 0 : componentIsNotPeriodic("exp-"+num);
334 0 : Value* comp=getPntrToComponent("exp-"+num);
335 0 : comp->set(coupl[i]);
336 : }
337 : }
338 : } else {
339 53 : for(unsigned i=0; i<ndata; i++) {
340 53 : std::string num; Tools::convert(i,num);
341 106 : addComponentWithDerivatives("rdc-"+num);
342 106 : componentIsNotPeriodic("rdc-"+num);
343 : }
344 53 : for(unsigned i=0; i<ndata; i++) {
345 53 : std::string num; Tools::convert(i,num);
346 106 : addComponent("exp-"+num);
347 106 : componentIsNotPeriodic("exp-"+num);
348 106 : Value* comp=getPntrToComponent("exp-"+num);
349 106 : comp->set(coupl[i]);
350 : }
351 : }
352 :
353 29 : if(svd) {
354 3 : addComponent("Sxx"); componentIsNotPeriodic("Sxx");
355 3 : addComponent("Syy"); componentIsNotPeriodic("Syy");
356 3 : addComponent("Szz"); componentIsNotPeriodic("Szz");
357 3 : addComponent("Sxy"); componentIsNotPeriodic("Sxy");
358 3 : addComponent("Sxz"); componentIsNotPeriodic("Sxz");
359 3 : addComponent("Syz"); componentIsNotPeriodic("Syz");
360 : }
361 :
362 29 : requestAtoms(atoms, false);
363 29 : if(getDoScore()) {
364 12 : setParameters(coupl);
365 12 : Initialise(coupl.size());
366 : }
367 29 : setDerivatives();
368 29 : checkRead();
369 29 : }
370 :
371 5 : void RDC::do_svd()
372 : {
373 : #ifdef __PLUMED_HAS_GSL
374 5 : gsl_vector_unique_ptr rdc_vec(gsl_vector_alloc(coupl.size())),
375 5 : S(gsl_vector_alloc(5)),
376 5 : Stmp(gsl_vector_alloc(5)),
377 5 : work(gsl_vector_alloc(5)),
378 5 : bc(gsl_vector_alloc(coupl.size()));
379 :
380 5 : gsl_matrix_unique_ptr coef_mat(gsl_matrix_alloc(coupl.size(),5)),
381 5 : A(gsl_matrix_alloc(coupl.size(),5)),
382 5 : V(gsl_matrix_alloc(5,5));
383 :
384 5 : gsl_matrix_set_zero(coef_mat.get());
385 5 : gsl_vector_set_zero(bc.get());
386 :
387 : unsigned index=0;
388 5 : vector<double> dmax(coupl.size());
389 60 : for(unsigned r=0; r<getNumberOfAtoms(); r+=2) {
390 25 : Vector distance;
391 75 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
392 0 : else distance = delta(getPosition(r),getPosition(r+1));
393 25 : double d = distance.modulo();
394 25 : double d2 = d*d;
395 25 : double d3 = d2*d;
396 25 : double id3 = 1./d3;
397 25 : double max = -Const*mu_s*scale;
398 50 : dmax[index] = id3*max;
399 25 : double mu_x = distance[0]/d;
400 25 : double mu_y = distance[1]/d;
401 25 : double mu_z = distance[2]/d;
402 50 : gsl_vector_set(rdc_vec.get(),index,coupl[index]/dmax[index]);
403 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));
404 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));
405 25 : gsl_matrix_set(coef_mat.get(),index,2,gsl_matrix_get(coef_mat.get(),index,2)+(2.0*mu_x*mu_y));
406 25 : gsl_matrix_set(coef_mat.get(),index,3,gsl_matrix_get(coef_mat.get(),index,3)+(2.0*mu_x*mu_z));
407 25 : gsl_matrix_set(coef_mat.get(),index,4,gsl_matrix_get(coef_mat.get(),index,4)+(2.0*mu_y*mu_z));
408 25 : index++;
409 : }
410 5 : gsl_matrix_memcpy(A.get(),coef_mat.get());
411 5 : gsl_linalg_SV_decomp(A.get(), V.get(), Stmp.get(), work.get());
412 5 : gsl_linalg_SV_solve(A.get(), V.get(), Stmp.get(), rdc_vec.get(), S.get());
413 : /* tensor */
414 : Value* tensor;
415 10 : tensor=getPntrToComponent("Sxx");
416 5 : double Sxx = gsl_vector_get(S.get(),0);
417 : tensor->set(Sxx);
418 10 : tensor=getPntrToComponent("Syy");
419 5 : double Syy = gsl_vector_get(S.get(),1);
420 : tensor->set(Syy);
421 10 : tensor=getPntrToComponent("Szz");
422 5 : double Szz = -Sxx-Syy;
423 : tensor->set(Szz);
424 10 : tensor=getPntrToComponent("Sxy");
425 5 : double Sxy = gsl_vector_get(S.get(),2);
426 : tensor->set(Sxy);
427 10 : tensor=getPntrToComponent("Sxz");
428 5 : double Sxz = gsl_vector_get(S.get(),3);
429 : tensor->set(Sxz);
430 10 : tensor=getPntrToComponent("Syz");
431 5 : double Syz = gsl_vector_get(S.get(),4);
432 : tensor->set(Syz);
433 :
434 5 : gsl_blas_dgemv(CblasNoTrans, 1.0, coef_mat.get(), S.get(), 0., bc.get());
435 55 : for(index=0; index<coupl.size(); index++) {
436 50 : double rdc = gsl_vector_get(bc.get(),index)*dmax[index];
437 25 : Value* val=getPntrToComponent(index);
438 : val->set(rdc);
439 : }
440 : #endif
441 5 : }
442 :
443 2172 : void RDC::calculate()
444 : {
445 2172 : if(svd) {
446 5 : do_svd();
447 2177 : return;
448 : }
449 :
450 2167 : const double max = -Const*scale*mu_s;
451 : const unsigned N=getNumberOfAtoms();
452 2167 : vector<Vector> dRDC(N/2, Vector{0.,0.,0.});
453 :
454 : /* RDC Calculations and forces */
455 6348 : #pragma omp parallel num_threads(OpenMP::getNumThreads())
456 : {
457 : #pragma omp for
458 : for(unsigned r=0; r<N; r+=2)
459 : {
460 8605 : const unsigned index=r/2;
461 8605 : Vector distance;
462 34340 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
463 0 : else distance = delta(getPosition(r),getPosition(r+1));
464 8575 : const double d2 = distance.modulo2();
465 8613 : const double ind = 1./sqrt(d2);
466 8613 : const double ind2 = 1./d2;
467 8613 : const double ind3 = ind2*ind;
468 8613 : const double x2 = distance[0]*distance[0]*ind2;
469 8635 : const double y2 = distance[1]*distance[1]*ind2;
470 8658 : const double z2 = distance[2]*distance[2]*ind2;
471 8657 : const double dmax = ind3*max;
472 8657 : const double ddmax = dmax*ind2;
473 :
474 8657 : const double rdc = 0.5*dmax*(3.*z2-1.);
475 8657 : const double prod_xy = (x2+y2-4.*z2);
476 8657 : const double prod_z = (3.*x2 + 3.*y2 - 2.*z2);
477 :
478 17314 : dRDC[index] = -1.5*ddmax*distance;
479 17246 : dRDC[index][0] *= prod_xy;
480 17264 : dRDC[index][1] *= prod_xy;
481 17290 : dRDC[index][2] *= prod_z;
482 :
483 8649 : string num; Tools::convert(index,num);
484 17255 : Value* val=getPntrToComponent("rdc-"+num);
485 : val->set(rdc);
486 8631 : if(!getDoScore()) {
487 2724 : setBoxDerivatives(val, Tensor(distance,dRDC[index]));
488 2738 : setAtomsDerivatives(val, r, dRDC[index]);
489 2726 : setAtomsDerivatives(val, r+1, -dRDC[index]);
490 : } else setCalcData(index, rdc);
491 : }
492 : }
493 :
494 2167 : if(getDoScore()) {
495 : /* Metainference */
496 1824 : Tensor dervir;
497 1824 : double score = getScore();
498 : setScore(score);
499 :
500 : /* calculate final derivatives */
501 3648 : Value* val=getPntrToComponent("score");
502 9120 : for(unsigned r=0; r<N; r+=2)
503 : {
504 7296 : const unsigned index=r/2;
505 7296 : Vector distance;
506 21888 : if(pbc) distance = pbcDistance(getPosition(r),getPosition(r+1));
507 0 : else distance = delta(getPosition(r),getPosition(r+1));
508 7296 : const Vector der = dRDC[index]*getMetaDer(index);
509 7296 : dervir += Tensor(distance, der);
510 7296 : setAtomsDerivatives(val, r, der);
511 7296 : setAtomsDerivatives(val, r+1, -der);
512 : }
513 1824 : setBoxDerivatives(val, dervir);
514 : }
515 : }
516 :
517 327 : void RDC::update() {
518 : // write status file
519 654 : if(getWstride()>0&& (getStep()%getWstride()==0 || getCPT()) ) writeStatus();
520 327 : }
521 :
522 : }
523 5874 : }
|
