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