Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 : #include "Colvar.h"
23 : #include "core/PlumedMain.h"
24 : #include "ActionRegister.h"
25 : #include "tools/Torsion.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 : #include <iostream>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace colvar {
35 :
36 : //+PLUMEDOC COLVAR PUCKERING
37 : /*
38 : Calculate sugar pseudorotation coordinates.
39 :
40 : This command can be used to calculate ring's pseudorotations in sugars (puckers). It works for both
41 : 5-membered and 6-membered rings. Notice that there are two different implementations depending if
42 : one passes 5 or 6 atoms in the ATOMS keyword.
43 :
44 : For 5-membered rings the implementation is the one discussed in \cite huang2014improvement .
45 : This implementation is simple and can be used in RNA to distinguish C2'-endo and C3'-endo conformations.
46 : Both the polar coordinates (phs and amp) and the cartesian coordinates (Zx and Zy) are provided.
47 : C2'-endo conformations have negative Zx, whereas C3'-endo conformations have positive Zy.
48 : Notation is consistent with \cite huang2014improvement .
49 : The five atoms should be provided as C4',O4',C1',C2',C3'.
50 : Notice that this is the same order that can be obtained using the \ref MOLINFO syntax (see example below).
51 :
52 : For 6-membered rings the implementation is the general Cremer-Pople one \cite cremer1975general
53 : as also discussed in \cite biarnes2007conformational .
54 : This implementation provides both a triplet with Cartesian components (qx, qy, and qz)
55 : and a triplet of polar components (amplitude, phi, and theta).
56 : Applications of this particular implementation are to be published (paper in preparation).
57 :
58 : \bug The 6-membered ring implementation distributed with this version of PLUMED leads to
59 : qx and qy values that have an opposite sign with respect to those originally defined in \cite cremer1975general.
60 : The bug will be fixed in a later version but is here kept to preserve reproducibility.
61 :
62 : Components of this action are:
63 :
64 : \par Examples
65 :
66 : This input tells plumed to print the puckering phase angle of the 3rd nucleotide of a RNA molecule on file COLVAR.
67 : \verbatim
68 : MOLINFO STRUCTURE=rna.pdb MOLTYPE=rna
69 : PUCKERING ATOMS=@sugar-3 LABEL=puck
70 : PRINT ARG=puck.phs FILE=COLVAR
71 : \endverbatim
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 30 : class Puckering : public Colvar {
77 :
78 : public:
79 : explicit Puckering(const ActionOptions&);
80 : virtual void calculate();
81 : static void registerKeywords(Keywords& keys);
82 : void calculate5m();
83 : void calculate6m();
84 : };
85 :
86 2538 : PLUMED_REGISTER_ACTION(Puckering,"PUCKERING")
87 :
88 16 : void Puckering::registerKeywords(Keywords& keys) {
89 16 : Colvar::registerKeywords( keys );
90 16 : keys.remove("NOPBC");
91 16 : keys.add("atoms","ATOMS","the five or six atoms of the sugar ring in the proper order");
92 16 : keys.addOutputComponent("phs","default","Pseudorotation phase (5 membered rings)");
93 16 : keys.addOutputComponent("amp","default","Pseudorotation amplitude (5 membered rings)");
94 16 : keys.addOutputComponent("Zx","default","Pseudorotation x cartesian component (5 membered rings)");
95 16 : keys.addOutputComponent("Zy","default","Pseudorotation y cartesian component (5 membered rings)");
96 16 : keys.addOutputComponent("phi","default","Pseudorotation phase (6 membered rings)");
97 16 : keys.addOutputComponent("theta","default","Theta angle (6 membered rings)");
98 16 : keys.addOutputComponent("amplitude","default","Pseudorotation amplitude (6 membered rings)");
99 16 : keys.addOutputComponent("qx","default","Cartesian component x (6 membered rings)");
100 16 : keys.addOutputComponent("qy","default","Cartesian component y (6 membered rings)");
101 16 : keys.addOutputComponent("qz","default","Cartesian component z (6 membered rings)");
102 16 : }
103 :
104 15 : Puckering::Puckering(const ActionOptions&ao):
105 15 : PLUMED_COLVAR_INIT(ao)
106 : {
107 15 : vector<AtomNumber> atoms;
108 15 : parseAtomList("ATOMS",atoms);
109 15 : if(atoms.size()!=5 && atoms.size()!=6) error("only for 5 or 6-membered rings");
110 15 : checkRead();
111 :
112 15 : if(atoms.size()==5) {
113 14 : log.printf(" between atoms %d %d %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial(),atoms[4].serial());
114 1 : } else if(atoms.size()==6) {
115 1 : log.printf(" between atoms %d %d %d %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial(),atoms[4].serial(),atoms[5].serial());
116 0 : } else error("ATOMS should specify 5 atoms");
117 :
118 15 : if(atoms.size()==5) {
119 14 : addComponentWithDerivatives("phs"); componentIsPeriodic("phs","-pi","pi");
120 14 : addComponentWithDerivatives("amp"); componentIsNotPeriodic("amp");
121 14 : addComponentWithDerivatives("Zx"); componentIsNotPeriodic("Zx");
122 14 : addComponentWithDerivatives("Zy"); componentIsNotPeriodic("Zy");
123 1 : } else if(atoms.size()==6) {
124 1 : addComponentWithDerivatives("qx"); componentIsNotPeriodic("qx");
125 1 : addComponentWithDerivatives("qy"); componentIsNotPeriodic("qy");
126 1 : addComponentWithDerivatives("qz"); componentIsNotPeriodic("qz");
127 1 : addComponentWithDerivatives("phi"); componentIsPeriodic("phi","0","2pi");
128 1 : addComponentWithDerivatives("theta"); componentIsNotPeriodic("theta");
129 1 : addComponentWithDerivatives("amplitude"); componentIsNotPeriodic("amplitude");
130 : }
131 :
132 15 : log<<" Bibliography ";
133 15 : if(atoms.size()==5) log<<plumed.cite("Huang, Giese, Lee, York, J. Chem. Theory Comput. 10, 1538 (2014)");
134 15 : if(atoms.size()==6) log<<plumed.cite("Cremer and Pople, J. Am. Chem. Soc. 97, 1354 (1975)");
135 15 : log<<"\n";
136 :
137 15 : requestAtoms(atoms);
138 15 : }
139 :
140 : // calculator
141 383 : void Puckering::calculate() {
142 383 : makeWhole();
143 383 : if(getNumberOfAtoms()==5) calculate5m();
144 103 : else calculate6m();
145 383 : }
146 :
147 280 : void Puckering::calculate5m() {
148 :
149 280 : Vector d0,d1,d2,d3,d4,d5;
150 280 : makeWhole();
151 :
152 280 : d0=delta(getPosition(2),getPosition(1));
153 280 : d1=delta(getPosition(3),getPosition(2));
154 280 : d2=delta(getPosition(4),getPosition(3));
155 280 : d3=delta(getPosition(4),getPosition(3));
156 280 : d4=delta(getPosition(0),getPosition(4));
157 280 : d5=delta(getPosition(1),getPosition(0));
158 :
159 280 : Vector r[5];
160 280 : r[0]=getPosition(0);
161 1400 : for(unsigned i=1; i<5; i++) {
162 1120 : r[i]=r[i-1]+pbcDistance(getPosition(i-1),getPosition(i));
163 : }
164 :
165 280 : Vector dd0,dd1,dd2,dd3,dd4,dd5;
166 :
167 : PLMD::Torsion t;
168 :
169 280 : double v1=t.compute(d0,d1,d2,dd0,dd1,dd2);
170 280 : double v3=t.compute(d3,d4,d5,dd3,dd4,dd5);
171 :
172 280 : double Zx=(v1+v3)/(2.0*cos(4.0*pi/5.0));
173 280 : double Zy=(v1-v3)/(2.0*sin(4.0*pi/5.0));
174 280 : double phase=atan2(Zy,Zx);
175 280 : double amplitude=sqrt(Zx*Zx+Zy*Zy);
176 :
177 280 : Vector dZx_dR[5];
178 280 : Vector dZy_dR[5];
179 :
180 280 : dZx_dR[0]=(dd5-dd4);
181 280 : dZx_dR[1]=(dd0-dd5);
182 280 : dZx_dR[2]=(dd1-dd0);
183 280 : dZx_dR[3]=(dd2+dd3-dd1);
184 280 : dZx_dR[4]=(dd4-dd3-dd2);
185 :
186 280 : dZy_dR[0]=(dd4-dd5);
187 280 : dZy_dR[1]=(dd0+dd5);
188 280 : dZy_dR[2]=(dd1-dd0);
189 280 : dZy_dR[3]=(dd2-dd3-dd1);
190 280 : dZy_dR[4]=(dd3-dd4-dd2);
191 :
192 280 : for(unsigned j=0; j<5; j++) dZx_dR[j]*=(1.0/(2.0*cos(4.0*pi/5.0)));
193 280 : for(unsigned j=0; j<5; j++) dZy_dR[j]*=(1.0/(2.0*sin(4.0*pi/5.0)));
194 :
195 280 : Vector dphase_dR[5];
196 280 : for(unsigned j=0; j<5; j++) dphase_dR[j]=(1.0/(Zx*Zx+Zy*Zy))*(-Zy*dZx_dR[j] + Zx*dZy_dR[j]);
197 :
198 280 : Vector damplitude_dR[5];
199 280 : for(unsigned j=0; j<5; j++) damplitude_dR[j]=(1.0/amplitude)*(Zx*dZx_dR[j] + Zy*dZy_dR[j]);
200 :
201 280 : Value* vzx=getPntrToComponent("Zx");
202 280 : vzx->set(Zx);
203 280 : setAtomsDerivatives (vzx,0, dZx_dR[0]);
204 280 : setAtomsDerivatives (vzx,1, dZx_dR[1]);
205 280 : setAtomsDerivatives (vzx,2, dZx_dR[2]);
206 280 : setAtomsDerivatives (vzx,3, dZx_dR[3]);
207 280 : setAtomsDerivatives (vzx,4, dZx_dR[4]);
208 280 : setBoxDerivativesNoPbc(vzx);
209 :
210 280 : Value* vzy=getPntrToComponent("Zy");
211 280 : vzy->set(Zy);
212 280 : setAtomsDerivatives (vzy,0, dZy_dR[0]);
213 280 : setAtomsDerivatives (vzy,1, dZy_dR[1]);
214 280 : setAtomsDerivatives (vzy,2, dZy_dR[2]);
215 280 : setAtomsDerivatives (vzy,3, dZy_dR[3]);
216 280 : setAtomsDerivatives (vzy,4, dZy_dR[4]);
217 280 : setBoxDerivativesNoPbc(vzy);
218 :
219 :
220 280 : Value* vph=getPntrToComponent("phs");
221 280 : vph->set(phase);
222 280 : setAtomsDerivatives (vph,0, dphase_dR[0]);
223 280 : setAtomsDerivatives (vph,1, dphase_dR[1]);
224 280 : setAtomsDerivatives (vph,2, dphase_dR[2]);
225 280 : setAtomsDerivatives (vph,3, dphase_dR[3]);
226 280 : setAtomsDerivatives (vph,4, dphase_dR[4]);
227 280 : setBoxDerivativesNoPbc(vph);
228 :
229 280 : Value* vam=getPntrToComponent("amp");
230 280 : vam->set(amplitude);
231 280 : setAtomsDerivatives (vam,0, damplitude_dR[0]);
232 280 : setAtomsDerivatives (vam,1, damplitude_dR[1]);
233 280 : setAtomsDerivatives (vam,2, damplitude_dR[2]);
234 280 : setAtomsDerivatives (vam,3, damplitude_dR[3]);
235 280 : setAtomsDerivatives (vam,4, damplitude_dR[4]);
236 280 : setBoxDerivativesNoPbc(vam);
237 :
238 :
239 280 : }
240 :
241 103 : void Puckering::calculate6m() {
242 :
243 103 : vector<Vector> r(6);
244 103 : r[0]=getPosition(0);
245 618 : for(unsigned i=1; i<6; i++) {
246 515 : r[i]=r[i-1]+pbcDistance(getPosition(i-1),getPosition(i));
247 : }
248 :
249 206 : vector<Vector> R(6);
250 103 : Vector center;
251 103 : for(unsigned j=0; j<6; j++) center+=r[j]/6.0;
252 103 : for(unsigned j=0; j<6; j++) R[j]=(r[j]-center);
253 :
254 103 : Vector Rp,Rpp;
255 103 : for(unsigned j=0; j<6; j++) Rp +=R[j]*sin(2.0/6.0*pi*j);
256 103 : for(unsigned j=0; j<6; j++) Rpp+=R[j]*cos(2.0/6.0*pi*j);
257 :
258 103 : Vector n=crossProduct(Rp,Rpp);
259 103 : Vector nhat=n/modulo(n);
260 :
261 103 : Tensor dn_dRp=dcrossDv1(Rp,Rpp);
262 103 : Tensor dn_dRpp=dcrossDv2(Rp,Rpp);
263 :
264 103 : Tensor dnhat_dn= 1.0/modulo(n)*( Tensor::identity() - extProduct(nhat,nhat));
265 103 : Tensor dnhat_dRp=matmul(dnhat_dn,dn_dRp);
266 103 : Tensor dnhat_dRpp=matmul(dnhat_dn,dn_dRpp);
267 :
268 206 : vector<double> z(6);
269 103 : for(unsigned j=0; j<6; j++) z[j]=dotProduct(R[j],nhat);
270 :
271 206 : vector<vector<Vector> > dz_dR(6);
272 103 : for(unsigned j=0; j<6; j++) dz_dR[j].resize(6);
273 :
274 3811 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
275 3708 : if(i==j) dz_dR[i][j]+=nhat;
276 3708 : dz_dR[i][j]+=matmul(R[i],dnhat_dRp)*sin(2.0/6.0*pi*j);
277 3708 : dz_dR[i][j]+=matmul(R[i],dnhat_dRpp)*cos(2.0/6.0*pi*j);
278 : }
279 :
280 103 : double B=0.0;
281 103 : for(unsigned j=0; j<6; j++) B+=z[j]*cos(4.0/6.0*pi*j);
282 :
283 206 : vector<Vector> dB_dR(6);
284 3811 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
285 3708 : dB_dR[i]+=dz_dR[j][i]*cos(4.0/6.0*pi*j);
286 : }
287 103 : Vector Bsum;
288 103 : for(unsigned j=0; j<6; j++) Bsum+=dB_dR[j];
289 103 : for(unsigned j=0; j<6; j++) dB_dR[j]-=Bsum/6.0;;
290 :
291 103 : double A=0.0;
292 103 : for(unsigned j=0; j<6; j++) A+=z[j]*sin(4.0/6.0*pi*j);
293 :
294 206 : vector<Vector> dA_dR(6);
295 3811 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
296 3708 : dA_dR[i]+=dz_dR[j][i]*sin(4.0/6.0*pi*j);
297 : }
298 103 : Vector Asum;
299 103 : for(unsigned j=0; j<6; j++) Asum+=dA_dR[j];
300 103 : for(unsigned j=0; j<6; j++) dA_dR[j]-=Asum/6.0;;
301 :
302 103 : double C=0.0;
303 103 : for(unsigned j=0; j<6; j++) C+=z[j]*Tools::fastpow(-1.0,(j));
304 :
305 206 : vector<Vector> dC_dR(6);
306 3811 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
307 3708 : dC_dR[i]+=dz_dR[j][i]*Tools::fastpow(-1.0,(j));
308 : }
309 :
310 103 : Vector Csum;
311 103 : for(unsigned j=0; j<6; j++) Csum+=dC_dR[j];
312 103 : for(unsigned j=0; j<6; j++) dC_dR[j]-=Csum/6.0;;
313 :
314 :
315 : // qx
316 103 : double qx = A/sqrt(3);
317 :
318 : // qx derivaties
319 206 : vector<Vector> dqx_dR(6);
320 721 : for(unsigned j=0; j<6; j++) {
321 618 : dqx_dR[j]=1/sqrt(3) * dA_dR[j];
322 : }
323 :
324 103 : Value* vqx=getPntrToComponent("qx");
325 103 : vqx->set(qx);
326 103 : setAtomsDerivatives (vqx,0, dqx_dR[0] );
327 103 : setAtomsDerivatives (vqx,1, dqx_dR[1] );
328 103 : setAtomsDerivatives (vqx,2, dqx_dR[2] );
329 103 : setAtomsDerivatives (vqx,3, dqx_dR[3] );
330 103 : setAtomsDerivatives (vqx,4, dqx_dR[4] );
331 103 : setAtomsDerivatives (vqx,5, dqx_dR[5] );
332 103 : setBoxDerivativesNoPbc(vqx);
333 :
334 : // qy
335 103 : double qy = -B/sqrt(3);
336 :
337 : // qy derivatives
338 206 : vector<Vector> dqy_dR(6);
339 721 : for(unsigned j=0; j<6; j++) {
340 618 : dqy_dR[j]=-1/sqrt(3) * dB_dR[j];
341 : }
342 :
343 103 : Value* vqy=getPntrToComponent("qy");
344 103 : vqy->set(qy);
345 103 : setAtomsDerivatives (vqy,0, dqy_dR[0] );
346 103 : setAtomsDerivatives (vqy,1, dqy_dR[1] );
347 103 : setAtomsDerivatives (vqy,2, dqy_dR[2] );
348 103 : setAtomsDerivatives (vqy,3, dqy_dR[3] );
349 103 : setAtomsDerivatives (vqy,4, dqy_dR[4] );
350 103 : setAtomsDerivatives (vqy,5, dqy_dR[5] );
351 103 : setBoxDerivativesNoPbc(vqy);
352 :
353 : // qz
354 103 : double qz = C/sqrt(6);
355 :
356 : // qz derivatives
357 206 : vector<Vector> dqz_dR(6);
358 721 : for(unsigned j=0; j<6; j++) {
359 618 : dqz_dR[j]=1/sqrt(6) * dC_dR[j];
360 : }
361 :
362 103 : Value* vqz=getPntrToComponent("qz");
363 103 : vqz->set(qz);
364 103 : setAtomsDerivatives (vqz,0, dqz_dR[0] );
365 103 : setAtomsDerivatives (vqz,1, dqz_dR[1] );
366 103 : setAtomsDerivatives (vqz,2, dqz_dR[2] );
367 103 : setAtomsDerivatives (vqz,3, dqz_dR[3] );
368 103 : setAtomsDerivatives (vqz,4, dqz_dR[4] );
369 103 : setAtomsDerivatives (vqz,5, dqz_dR[5] );
370 103 : setBoxDerivativesNoPbc(vqz);
371 :
372 :
373 : // PHASE
374 103 : double phi=atan2(-A,B);
375 :
376 : // PHASE DERIVATIVES
377 206 : vector<Vector> dphi_dR(6);
378 721 : for(unsigned j=0; j<6; j++) {
379 618 : dphi_dR[j]=1.0/(A*A+B*B) * (-B*dA_dR[j] + A*dB_dR[j]);
380 : }
381 :
382 103 : Value* vphi=getPntrToComponent("phi");
383 103 : vphi->set(phi);
384 103 : setAtomsDerivatives (vphi,0, dphi_dR[0] );
385 103 : setAtomsDerivatives (vphi,1, dphi_dR[1] );
386 103 : setAtomsDerivatives (vphi,2, dphi_dR[2] );
387 103 : setAtomsDerivatives (vphi,3, dphi_dR[3] );
388 103 : setAtomsDerivatives (vphi,4, dphi_dR[4] );
389 103 : setAtomsDerivatives (vphi,5, dphi_dR[5] );
390 103 : setBoxDerivativesNoPbc(vphi);
391 :
392 : // AMPLITUDE
393 103 : double amplitude=sqrt((2*(A*A+B*B)+C*C)/6);
394 :
395 : // AMPLITUDE DERIVATIES
396 206 : vector<Vector> damplitude_dR(6);
397 721 : for (unsigned j=0; j<6; j++) {
398 618 : damplitude_dR[j]=0.5*sqrt(2.0/6.0)/(sqrt(A*A+B*B+0.5*C*C)) * (2*A*dA_dR[j] + 2*B*dB_dR[j] + C*dC_dR[j]);
399 : }
400 :
401 103 : Value* vamplitude=getPntrToComponent("amplitude");
402 103 : vamplitude->set(amplitude);
403 103 : setAtomsDerivatives (vamplitude,0, damplitude_dR[0] );
404 103 : setAtomsDerivatives (vamplitude,1, damplitude_dR[1] );
405 103 : setAtomsDerivatives (vamplitude,2, damplitude_dR[2] );
406 103 : setAtomsDerivatives (vamplitude,3, damplitude_dR[3] );
407 103 : setAtomsDerivatives (vamplitude,4, damplitude_dR[4] );
408 103 : setAtomsDerivatives (vamplitude,5, damplitude_dR[5] );
409 103 : setBoxDerivativesNoPbc(vamplitude);
410 :
411 : // THETA
412 103 : double theta=acos( C / sqrt(2.*(A*A+B*B) +C*C ) );
413 :
414 : // THETA DERIVATIVES
415 206 : vector<Vector> dtheta_dR(6);
416 721 : for(unsigned j=0; j<6; j++) {
417 618 : dtheta_dR[j]=1.0/(3.0*sqrt(2)*amplitude*amplitude) * (C/(sqrt(A*A+B*B)) * (A*dA_dR[j] + B*dB_dR[j]) - sqrt(A*A+B*B)*dC_dR[j]);
418 : }
419 103 : Value* vtheta=getPntrToComponent("theta");
420 103 : vtheta->set(theta);
421 103 : setAtomsDerivatives (vtheta,0, dtheta_dR[0] );
422 103 : setAtomsDerivatives (vtheta,1, dtheta_dR[1] );
423 103 : setAtomsDerivatives (vtheta,2, dtheta_dR[2] );
424 103 : setAtomsDerivatives (vtheta,3, dtheta_dR[3] );
425 103 : setAtomsDerivatives (vtheta,4, dtheta_dR[4] );
426 103 : setAtomsDerivatives (vtheta,5, dtheta_dR[5] );
427 206 : setBoxDerivativesNoPbc(vtheta);
428 103 : }
429 :
430 :
431 : }
432 2523 : }
433 :
434 :
435 :
|
