Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "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 : \plumedfile
68 : MOLINFO STRUCTURE=rna.pdb MOLTYPE=rna
69 : PUCKERING ATOMS=@sugar-3 LABEL=puck
70 : PRINT ARG=puck.phs FILE=COLVAR
71 : \endplumedfile
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 106 : 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 6506 : PLUMED_REGISTER_ACTION(Puckering,"PUCKERING")
87 :
88 55 : void Puckering::registerKeywords(Keywords& keys) {
89 55 : Colvar::registerKeywords( keys );
90 110 : keys.remove("NOPBC");
91 220 : keys.add("atoms","ATOMS","the five or six atoms of the sugar ring in the proper order");
92 220 : keys.addOutputComponent("phs","default","Pseudorotation phase (5 membered rings)");
93 220 : keys.addOutputComponent("amp","default","Pseudorotation amplitude (5 membered rings)");
94 220 : keys.addOutputComponent("Zx","default","Pseudorotation x cartesian component (5 membered rings)");
95 220 : keys.addOutputComponent("Zy","default","Pseudorotation y cartesian component (5 membered rings)");
96 220 : keys.addOutputComponent("phi","default","Pseudorotation phase (6 membered rings)");
97 220 : keys.addOutputComponent("theta","default","Theta angle (6 membered rings)");
98 220 : keys.addOutputComponent("amplitude","default","Pseudorotation amplitude (6 membered rings)");
99 220 : keys.addOutputComponent("qx","default","Cartesian component x (6 membered rings)");
100 220 : keys.addOutputComponent("qy","default","Cartesian component y (6 membered rings)");
101 220 : keys.addOutputComponent("qz","default","Cartesian component z (6 membered rings)");
102 55 : }
103 :
104 54 : Puckering::Puckering(const ActionOptions&ao):
105 55 : PLUMED_COLVAR_INIT(ao)
106 : {
107 : vector<AtomNumber> atoms;
108 108 : parseAtomList("ATOMS",atoms);
109 55 : if(atoms.size()!=5 && atoms.size()!=6) error("only for 5 or 6-membered rings");
110 53 : checkRead();
111 :
112 53 : if(atoms.size()==5) {
113 104 : 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 2 : 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 53 : if(atoms.size()==5) {
119 260 : addComponentWithDerivatives("phs"); componentIsPeriodic("phs","-pi","pi");
120 156 : addComponentWithDerivatives("amp"); componentIsNotPeriodic("amp");
121 156 : addComponentWithDerivatives("Zx"); componentIsNotPeriodic("Zx");
122 156 : addComponentWithDerivatives("Zy"); componentIsNotPeriodic("Zy");
123 1 : } else if(atoms.size()==6) {
124 3 : addComponentWithDerivatives("qx"); componentIsNotPeriodic("qx");
125 3 : addComponentWithDerivatives("qy"); componentIsNotPeriodic("qy");
126 3 : addComponentWithDerivatives("qz"); componentIsNotPeriodic("qz");
127 5 : addComponentWithDerivatives("phi"); componentIsPeriodic("phi","0","2pi");
128 3 : addComponentWithDerivatives("theta"); componentIsNotPeriodic("theta");
129 3 : addComponentWithDerivatives("amplitude"); componentIsNotPeriodic("amplitude");
130 : }
131 :
132 53 : log<<" Bibliography ";
133 157 : if(atoms.size()==5) log<<plumed.cite("Huang, Giese, Lee, York, J. Chem. Theory Comput. 10, 1538 (2014)");
134 55 : if(atoms.size()==6) log<<plumed.cite("Cremer and Pople, J. Am. Chem. Soc. 97, 1354 (1975)");
135 53 : log<<"\n";
136 :
137 53 : requestAtoms(atoms);
138 53 : }
139 :
140 : // calculator
141 397 : void Puckering::calculate() {
142 397 : makeWhole();
143 397 : if(getNumberOfAtoms()==5) calculate5m();
144 103 : else calculate6m();
145 397 : }
146 :
147 294 : void Puckering::calculate5m() {
148 :
149 294 : Vector d0,d1,d2,d3,d4,d5;
150 :
151 294 : d0=delta(getPosition(2),getPosition(1));
152 294 : d1=delta(getPosition(3),getPosition(2));
153 294 : d2=delta(getPosition(4),getPosition(3));
154 294 : d3=delta(getPosition(4),getPosition(3));
155 294 : d4=delta(getPosition(0),getPosition(4));
156 294 : d5=delta(getPosition(1),getPosition(0));
157 :
158 294 : Vector dd0,dd1,dd2,dd3,dd4,dd5;
159 :
160 : PLMD::Torsion t;
161 :
162 294 : double v1=t.compute(d0,d1,d2,dd0,dd1,dd2);
163 294 : double v3=t.compute(d3,d4,d5,dd3,dd4,dd5);
164 :
165 294 : double Zx=(v1+v3)/(2.0*cos(4.0*pi/5.0));
166 294 : double Zy=(v1-v3)/(2.0*sin(4.0*pi/5.0));
167 294 : double phase=atan2(Zy,Zx);
168 294 : double amplitude=sqrt(Zx*Zx+Zy*Zy);
169 :
170 1764 : Vector dZx_dR[5];
171 1764 : Vector dZy_dR[5];
172 :
173 294 : dZx_dR[0]=(dd5-dd4);
174 294 : dZx_dR[1]=(dd0-dd5);
175 294 : dZx_dR[2]=(dd1-dd0);
176 294 : dZx_dR[3]=(dd2+dd3-dd1);
177 294 : dZx_dR[4]=(dd4-dd3-dd2);
178 :
179 294 : dZy_dR[0]=(dd4-dd5);
180 294 : dZy_dR[1]=(dd0+dd5);
181 294 : dZy_dR[2]=(dd1-dd0);
182 294 : dZy_dR[3]=(dd2-dd3-dd1);
183 294 : dZy_dR[4]=(dd3-dd4-dd2);
184 :
185 1764 : for(unsigned j=0; j<5; j++) dZx_dR[j]*=(1.0/(2.0*cos(4.0*pi/5.0)));
186 1764 : for(unsigned j=0; j<5; j++) dZy_dR[j]*=(1.0/(2.0*sin(4.0*pi/5.0)));
187 :
188 1764 : Vector dphase_dR[5];
189 1764 : for(unsigned j=0; j<5; j++) dphase_dR[j]=(1.0/(Zx*Zx+Zy*Zy))*(-Zy*dZx_dR[j] + Zx*dZy_dR[j]);
190 :
191 1764 : Vector damplitude_dR[5];
192 1764 : for(unsigned j=0; j<5; j++) damplitude_dR[j]=(1.0/amplitude)*(Zx*dZx_dR[j] + Zy*dZy_dR[j]);
193 :
194 588 : Value* vzx=getPntrToComponent("Zx");
195 : vzx->set(Zx);
196 294 : setAtomsDerivatives (vzx,0, dZx_dR[0]);
197 294 : setAtomsDerivatives (vzx,1, dZx_dR[1]);
198 294 : setAtomsDerivatives (vzx,2, dZx_dR[2]);
199 294 : setAtomsDerivatives (vzx,3, dZx_dR[3]);
200 294 : setAtomsDerivatives (vzx,4, dZx_dR[4]);
201 294 : setBoxDerivativesNoPbc(vzx);
202 :
203 588 : Value* vzy=getPntrToComponent("Zy");
204 : vzy->set(Zy);
205 294 : setAtomsDerivatives (vzy,0, dZy_dR[0]);
206 294 : setAtomsDerivatives (vzy,1, dZy_dR[1]);
207 294 : setAtomsDerivatives (vzy,2, dZy_dR[2]);
208 294 : setAtomsDerivatives (vzy,3, dZy_dR[3]);
209 294 : setAtomsDerivatives (vzy,4, dZy_dR[4]);
210 294 : setBoxDerivativesNoPbc(vzy);
211 :
212 :
213 588 : Value* vph=getPntrToComponent("phs");
214 : vph->set(phase);
215 294 : setAtomsDerivatives (vph,0, dphase_dR[0]);
216 294 : setAtomsDerivatives (vph,1, dphase_dR[1]);
217 294 : setAtomsDerivatives (vph,2, dphase_dR[2]);
218 294 : setAtomsDerivatives (vph,3, dphase_dR[3]);
219 294 : setAtomsDerivatives (vph,4, dphase_dR[4]);
220 294 : setBoxDerivativesNoPbc(vph);
221 :
222 588 : Value* vam=getPntrToComponent("amp");
223 : vam->set(amplitude);
224 294 : setAtomsDerivatives (vam,0, damplitude_dR[0]);
225 294 : setAtomsDerivatives (vam,1, damplitude_dR[1]);
226 294 : setAtomsDerivatives (vam,2, damplitude_dR[2]);
227 294 : setAtomsDerivatives (vam,3, damplitude_dR[3]);
228 294 : setAtomsDerivatives (vam,4, damplitude_dR[4]);
229 294 : setBoxDerivativesNoPbc(vam);
230 :
231 :
232 294 : }
233 :
234 103 : void Puckering::calculate6m() {
235 :
236 103 : vector<Vector> r(6);
237 1957 : for(unsigned i=0; i<6; i++) r[i]=getPosition(i);
238 :
239 103 : vector<Vector> R(6);
240 103 : Vector center;
241 1339 : for(unsigned j=0; j<6; j++) center+=r[j]/6.0;
242 1339 : for(unsigned j=0; j<6; j++) R[j]=(r[j]-center);
243 :
244 103 : Vector Rp,Rpp;
245 1339 : for(unsigned j=0; j<6; j++) Rp +=R[j]*sin(2.0/6.0*pi*j);
246 1339 : for(unsigned j=0; j<6; j++) Rpp+=R[j]*cos(2.0/6.0*pi*j);
247 :
248 103 : Vector n=crossProduct(Rp,Rpp);
249 103 : Vector nhat=n/modulo(n);
250 :
251 103 : Tensor dn_dRp=dcrossDv1(Rp,Rpp);
252 103 : Tensor dn_dRpp=dcrossDv2(Rp,Rpp);
253 :
254 103 : Tensor dnhat_dn= 1.0/modulo(n)*( Tensor::identity() - extProduct(nhat,nhat));
255 103 : Tensor dnhat_dRp=matmul(dnhat_dn,dn_dRp);
256 103 : Tensor dnhat_dRpp=matmul(dnhat_dn,dn_dRpp);
257 :
258 103 : vector<double> z(6);
259 1339 : for(unsigned j=0; j<6; j++) z[j]=dotProduct(R[j],nhat);
260 :
261 206 : vector<vector<Vector> > dz_dR(6);
262 1339 : for(unsigned j=0; j<6; j++) dz_dR[j].resize(6);
263 :
264 4429 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
265 4326 : if(i==j) dz_dR[i][j]+=nhat;
266 11124 : dz_dR[i][j]+=matmul(R[i],dnhat_dRp)*sin(2.0/6.0*pi*j);
267 11124 : dz_dR[i][j]+=matmul(R[i],dnhat_dRpp)*cos(2.0/6.0*pi*j);
268 : }
269 :
270 : double B=0.0;
271 1339 : for(unsigned j=0; j<6; j++) B+=z[j]*cos(4.0/6.0*pi*j);
272 :
273 103 : vector<Vector> dB_dR(6);
274 4429 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
275 11124 : dB_dR[i]+=dz_dR[j][i]*cos(4.0/6.0*pi*j);
276 : }
277 103 : Vector Bsum;
278 1339 : for(unsigned j=0; j<6; j++) Bsum+=dB_dR[j];
279 1339 : for(unsigned j=0; j<6; j++) dB_dR[j]-=Bsum/6.0;;
280 :
281 : double A=0.0;
282 1339 : for(unsigned j=0; j<6; j++) A+=z[j]*sin(4.0/6.0*pi*j);
283 :
284 103 : vector<Vector> dA_dR(6);
285 4429 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
286 11124 : dA_dR[i]+=dz_dR[j][i]*sin(4.0/6.0*pi*j);
287 : }
288 103 : Vector Asum;
289 1339 : for(unsigned j=0; j<6; j++) Asum+=dA_dR[j];
290 1339 : for(unsigned j=0; j<6; j++) dA_dR[j]-=Asum/6.0;;
291 :
292 : double C=0.0;
293 1957 : for(unsigned j=0; j<6; j++) C+=z[j]*Tools::fastpow(-1.0,(j));
294 :
295 103 : vector<Vector> dC_dR(6);
296 4429 : for(unsigned i=0; i<6; i++) for(unsigned j=0; j<6; j++) {
297 14832 : dC_dR[i]+=dz_dR[j][i]*Tools::fastpow(-1.0,(j));
298 : }
299 :
300 103 : Vector Csum;
301 1339 : for(unsigned j=0; j<6; j++) Csum+=dC_dR[j];
302 1339 : for(unsigned j=0; j<6; j++) dC_dR[j]-=Csum/6.0;;
303 :
304 :
305 : // qx
306 103 : double qx = A/sqrt(3);
307 :
308 : // qx derivaties
309 103 : vector<Vector> dqx_dR(6);
310 1339 : for(unsigned j=0; j<6; j++) {
311 1236 : dqx_dR[j]=1/sqrt(3) * dA_dR[j];
312 : }
313 :
314 206 : Value* vqx=getPntrToComponent("qx");
315 : vqx->set(qx);
316 206 : setAtomsDerivatives (vqx,0, dqx_dR[0] );
317 103 : setAtomsDerivatives (vqx,1, dqx_dR[1] );
318 103 : setAtomsDerivatives (vqx,2, dqx_dR[2] );
319 103 : setAtomsDerivatives (vqx,3, dqx_dR[3] );
320 103 : setAtomsDerivatives (vqx,4, dqx_dR[4] );
321 103 : setAtomsDerivatives (vqx,5, dqx_dR[5] );
322 103 : setBoxDerivativesNoPbc(vqx);
323 :
324 : // qy
325 103 : double qy = -B/sqrt(3);
326 :
327 : // qy derivatives
328 103 : vector<Vector> dqy_dR(6);
329 1339 : for(unsigned j=0; j<6; j++) {
330 1236 : dqy_dR[j]=-1/sqrt(3) * dB_dR[j];
331 : }
332 :
333 206 : Value* vqy=getPntrToComponent("qy");
334 : vqy->set(qy);
335 103 : setAtomsDerivatives (vqy,0, dqy_dR[0] );
336 103 : setAtomsDerivatives (vqy,1, dqy_dR[1] );
337 103 : setAtomsDerivatives (vqy,2, dqy_dR[2] );
338 103 : setAtomsDerivatives (vqy,3, dqy_dR[3] );
339 103 : setAtomsDerivatives (vqy,4, dqy_dR[4] );
340 103 : setAtomsDerivatives (vqy,5, dqy_dR[5] );
341 103 : setBoxDerivativesNoPbc(vqy);
342 :
343 : // qz
344 103 : double qz = C/sqrt(6);
345 :
346 : // qz derivatives
347 103 : vector<Vector> dqz_dR(6);
348 1339 : for(unsigned j=0; j<6; j++) {
349 1236 : dqz_dR[j]=1/sqrt(6) * dC_dR[j];
350 : }
351 :
352 206 : Value* vqz=getPntrToComponent("qz");
353 : vqz->set(qz);
354 103 : setAtomsDerivatives (vqz,0, dqz_dR[0] );
355 103 : setAtomsDerivatives (vqz,1, dqz_dR[1] );
356 103 : setAtomsDerivatives (vqz,2, dqz_dR[2] );
357 103 : setAtomsDerivatives (vqz,3, dqz_dR[3] );
358 103 : setAtomsDerivatives (vqz,4, dqz_dR[4] );
359 103 : setAtomsDerivatives (vqz,5, dqz_dR[5] );
360 103 : setBoxDerivativesNoPbc(vqz);
361 :
362 :
363 : // PHASE
364 103 : double phi=atan2(-A,B);
365 :
366 : // PHASE DERIVATIVES
367 103 : vector<Vector> dphi_dR(6);
368 1339 : for(unsigned j=0; j<6; j++) {
369 2472 : dphi_dR[j]=1.0/(A*A+B*B) * (-B*dA_dR[j] + A*dB_dR[j]);
370 : }
371 :
372 206 : Value* vphi=getPntrToComponent("phi");
373 : vphi->set(phi);
374 103 : setAtomsDerivatives (vphi,0, dphi_dR[0] );
375 103 : setAtomsDerivatives (vphi,1, dphi_dR[1] );
376 103 : setAtomsDerivatives (vphi,2, dphi_dR[2] );
377 103 : setAtomsDerivatives (vphi,3, dphi_dR[3] );
378 103 : setAtomsDerivatives (vphi,4, dphi_dR[4] );
379 103 : setAtomsDerivatives (vphi,5, dphi_dR[5] );
380 103 : setBoxDerivativesNoPbc(vphi);
381 :
382 : // AMPLITUDE
383 103 : double amplitude=sqrt((2*(A*A+B*B)+C*C)/6);
384 :
385 : // AMPLITUDE DERIVATIES
386 103 : vector<Vector> damplitude_dR(6);
387 1339 : for (unsigned j=0; j<6; j++) {
388 3090 : 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]);
389 : }
390 :
391 206 : Value* vamplitude=getPntrToComponent("amplitude");
392 : vamplitude->set(amplitude);
393 103 : setAtomsDerivatives (vamplitude,0, damplitude_dR[0] );
394 103 : setAtomsDerivatives (vamplitude,1, damplitude_dR[1] );
395 103 : setAtomsDerivatives (vamplitude,2, damplitude_dR[2] );
396 103 : setAtomsDerivatives (vamplitude,3, damplitude_dR[3] );
397 103 : setAtomsDerivatives (vamplitude,4, damplitude_dR[4] );
398 103 : setAtomsDerivatives (vamplitude,5, damplitude_dR[5] );
399 103 : setBoxDerivativesNoPbc(vamplitude);
400 :
401 : // THETA
402 103 : double theta=acos( C / sqrt(2.*(A*A+B*B) +C*C ) );
403 :
404 : // THETA DERIVATIVES
405 103 : vector<Vector> dtheta_dR(6);
406 1339 : for(unsigned j=0; j<6; j++) {
407 3090 : 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]);
408 : }
409 206 : Value* vtheta=getPntrToComponent("theta");
410 : vtheta->set(theta);
411 103 : setAtomsDerivatives (vtheta,0, dtheta_dR[0] );
412 103 : setAtomsDerivatives (vtheta,1, dtheta_dR[1] );
413 103 : setAtomsDerivatives (vtheta,2, dtheta_dR[2] );
414 103 : setAtomsDerivatives (vtheta,3, dtheta_dR[3] );
415 103 : setAtomsDerivatives (vtheta,4, dtheta_dR[4] );
416 103 : setAtomsDerivatives (vtheta,5, dtheta_dR[5] );
417 103 : setBoxDerivativesNoPbc(vtheta);
418 103 : }
419 :
420 :
421 : }
422 4839 : }
423 :
424 :
425 :
|
