Line data    Source code

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