Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "RMSDVector.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "tools/PDB.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace colvar {
      29             : 
      30             : PLUMED_REGISTER_ACTION(RMSDVector,"RMSD_VECTOR")
      31             : 
      32         132 : void RMSDVector::registerKeywords(Keywords& keys) {
      33         132 :   ActionWithVector::registerKeywords(keys);
      34         132 :   keys.setDisplayName("RMSD");
      35         264 :   keys.addInputKeyword("compulsory","ARG","vector/matrix","the labels of two actions that you are calculating the RMSD between");
      36         132 :   keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed.  Should be OPTIMAL or SIMPLE.");
      37         132 :   keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
      38         132 :   keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
      39         132 :   keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
      40         132 :   keys.addFlag("UNORMALIZED",false,"by default the mean sequare deviation or root mean square deviation is calculated.  If this option is given no averaging is done");
      41         132 :   keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
      42         264 :   keys.addOutputComponent("disp","DISPLACEMENT","vector/matrix","the vector of displacements for the atoms");
      43         264 :   keys.addOutputComponent("dist","DISPLACEMENT","scalar/vector","the RMSD distance the atoms have moved");
      44         132 :   PTM::registerKeywords( keys );
      45         264 :   keys.setValueDescription("scalar/vector","a vector containing the RMSD between the instantaneous structure and each of the reference structures that were input");
      46         132 : }
      47             : 
      48          65 : RMSDVector::RMSDVector(const ActionOptions&ao):
      49             :   Action(ao),
      50             :   ActionWithVector(ao),
      51          65 :   firststep(true),
      52             :   input(ParallelActionsInput::create(getPbc())),
      53          65 :   taskmanager(this) {
      54          65 :   if( getPntrToArgument(0)->getRank()!=1 ) {
      55           0 :     error("first argument should be vector");
      56             :   }
      57          65 :   if( getPntrToArgument(1)->getRank()<1 ) {
      58           0 :     error("second argument should be matrix or a vector");
      59             :   }
      60          65 :   if( getPntrToArgument(1)->getRank()==1 ) {
      61          16 :     if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getNumberOfValues() ) {
      62           0 :       error("mismatch between sizes of input vectors");
      63             :     }
      64          49 :   } else if( getPntrToArgument(1)->getRank()==2 ) {
      65          49 :     if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getShape()[1] ) {
      66           0 :       error("mismatch between sizes of input vectors");
      67             :     }
      68             :   }
      69          65 :   if( getPntrToArgument(0)->getNumberOfValues()%3!=0 ) {
      70           0 :     error("number of components in input arguments should be multiple of three");
      71             :   }
      72             : 
      73             :   RMSDVectorData myinput;
      74          65 :   unsigned natoms = getPntrToArgument(0)->getNumberOfValues() / 3;
      75          65 :   myinput.type.assign("SIMPLE");
      76          65 :   parse("TYPE",myinput.type);
      77          65 :   parseFlag("SQUARED",myinput.squared);
      78          65 :   myinput.align.resize( natoms );
      79          65 :   parseVector("ALIGN",myinput.align);
      80          65 :   myinput.displace.resize( natoms );
      81          65 :   parseVector("DISPLACE",myinput.displace);
      82          65 :   bool unorm=false;
      83          65 :   parseFlag("UNORMALIZED",unorm);
      84          65 :   norm_weights=!unorm;
      85             :   double wa=0, wd=0;
      86          65 :   myinput.sqrtdisplace.resize( myinput.displace.size() );
      87         853 :   for(unsigned i=0; i<myinput.align.size(); ++i) {
      88         788 :     wa+=myinput.align[i];
      89         788 :     wd+=myinput.displace[i];
      90             :   }
      91             : 
      92          65 :   if( wa>epsilon ) {
      93          65 :     double iwa = 1. / wa;
      94         853 :     for(unsigned i=0; i<myinput.align.size(); ++i) {
      95         788 :       myinput.align[i] *= iwa;
      96             :     }
      97             :   } else {
      98           0 :     double iwa = 1. / natoms;
      99           0 :     for(unsigned i=0; i<myinput.align.size(); ++i) {
     100           0 :       myinput.align[i] = iwa;
     101             :     }
     102             :   }
     103          65 :   if( wd>epsilon ) {
     104          65 :     if( !norm_weights ) {
     105             :       wd = 1;
     106             :     }
     107          65 :     double iwd = 1. / wd;
     108         853 :     for(unsigned i=0; i<myinput.align.size(); ++i) {
     109         788 :       myinput.displace[i] *= iwd;
     110             :     }
     111             :   } else {
     112           0 :     double iwd = 1. / natoms;
     113           0 :     for(unsigned i=0; i<myinput.align.size(); ++i) {
     114           0 :       myinput.displace[i] = iwd;
     115             :     }
     116             :   }
     117         853 :   for(unsigned i=0; i<myinput.align.size(); ++i) {
     118         788 :     myinput.sqrtdisplace[i] = sqrt(myinput.displace[i]);
     119             :   }
     120             : 
     121          65 :   parseFlag("DISPLACEMENT",myinput.displacement);
     122          65 :   if( myinput.displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) {
     123          70 :     addComponentWithDerivatives("dist");
     124          70 :     componentIsNotPeriodic("dist");
     125          35 :     std::vector<std::size_t> shape( 1, getPntrToArgument(0)->getNumberOfValues() );
     126          35 :     addComponent( "disp", shape );
     127          35 :     componentIsNotPeriodic("disp");
     128          65 :   } else if( myinput.displacement ) {
     129           8 :     std::vector<std::size_t> shape( 1, getPntrToArgument(1)->getShape()[0] );
     130           8 :     addComponent( "dist", shape );
     131           8 :     componentIsNotPeriodic("dist");
     132           8 :     shape.resize(2);
     133           8 :     shape[0] = getPntrToArgument(1)->getShape()[0];
     134           8 :     shape[1] = getPntrToArgument(0)->getNumberOfValues();
     135           8 :     addComponent( "disp", shape );
     136           8 :     getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
     137           8 :     componentIsNotPeriodic("disp");
     138          30 :   } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) {
     139          13 :     addValue();
     140          13 :     setNotPeriodic();
     141             :   } else {
     142           9 :     std::vector<std::size_t> shape( 1, getPntrToArgument(1)->getShape()[0] );
     143           9 :     addValue( shape );
     144           9 :     setNotPeriodic();
     145           9 :   }
     146          65 :   if( getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getNumberOfValues()==0 ) {
     147          16 :     myinput.myrmsd.resize(1);
     148             :   } else {
     149          49 :     myinput.myrmsd.resize( getPntrToArgument(1)->getShape()[0] );
     150             :   }
     151             : 
     152          65 :   if( getPntrToArgument(1)->getRank()==1 )
     153          16 :     log.printf("  calculating RMSD distance between %d atoms. Distance between the avectors of atoms in %s and %s\n",
     154             :                natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
     155             :   else
     156          49 :     log.printf("  calculating RMSD distance between %zu sets of %d atoms. Distance between vector %s of atoms and matrix of configurations in %s\n",
     157             :                getPntrToArgument(1)->getShape()[0], natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
     158          65 :   log.printf("  method for alignment : %s \n",myinput.type.c_str() );
     159          65 :   if(myinput.squared) {
     160          34 :     log.printf("  chosen to use SQUARED option for MSD instead of RMSD\n");
     161             :   } else {
     162          31 :     log.printf("  using periodic boundary conditions\n");
     163             :   }
     164             :   // Setup the task manager
     165          65 :   if( getPntrToComponent(0)->getRank()>0 ) {
     166          17 :     taskmanager.setupParallelTaskManager( 3*natoms, 3*natoms );
     167             :   }
     168          65 :   taskmanager.setActionInput( myinput );
     169             :   // Setup the Parallel action input object
     170          65 :   input.noderiv = false;
     171          65 :   input.ncomponents = getNumberOfComponents();
     172             :   unsigned total_vals = 0;
     173         173 :   for(unsigned i=0; i<input.ncomponents; ++i) {
     174         108 :     total_vals += getPntrToComponent(i)->getNumberOfStoredValues();
     175             :   }
     176          65 :   input.nscalars = 1;
     177          65 :   if( myinput.displacement ) {
     178          43 :     input.nscalars = 1 + 3*natoms;
     179             :   }
     180          65 :   force_stash.resize( total_vals );
     181          65 :   ArgumentsBookkeeping abk;
     182          65 :   copyArgumentBookeeping( abk );
     183          65 :   input.setupArguments( abk );
     184          65 : }
     185             : 
     186       10100 : unsigned RMSDVector::getNumberOfDerivatives() {
     187       10100 :   return getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
     188             : }
     189             : 
     190      284646 : int RMSDVector::checkTaskIsActive( const unsigned& itask ) const {
     191      284646 :   return 1;
     192             : }
     193             : 
     194       25461 : void RMSDVector::setReferenceConfigurations() {
     195       25461 :   unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
     196             :   Vector center;
     197       25461 :   RMSDVectorData& myinput=taskmanager.getActionInput();
     198       25461 :   std::vector<Vector> pos( natoms );
     199       53027 :   for(unsigned jconf=0; jconf<myinput.myrmsd.size(); ++jconf) {
     200             :     center.zero();
     201      391094 :     for(unsigned i=0; i<pos.size(); ++i) {
     202     1454112 :       for(unsigned j=0; j<3; ++j) {
     203     1090584 :         pos[i][j] = getPntrToArgument(1)->get( (3*jconf+j)*pos.size() + i );
     204             :       }
     205      363528 :       center+=pos[i]*myinput.align[i];
     206             :     }
     207      391094 :     for(unsigned i=0; i<pos.size(); ++i) {
     208             :       pos[i] -= center;
     209             :     }
     210       27566 :     myinput.myrmsd[jconf].clear();
     211       27566 :     myinput.myrmsd[jconf].set(myinput.align,myinput.displace,pos,myinput.type,true,norm_weights);
     212             :   }
     213       25461 : }
     214             : 
     215             : // calculator
     216       34659 : void RMSDVector::calculate() {
     217       34659 :   if( firststep || !getPntrToArgument(1)->isConstant() ) {
     218       25357 :     setReferenceConfigurations();
     219       25357 :     firststep=false;
     220             :   }
     221       34659 :   input.noderiv = false;
     222       34659 :   getArgumentInputData( input_buffer );
     223       34659 :   input.dataSize = input_buffer.size();
     224       34659 :   input.inputdata = input_buffer.data();
     225             : 
     226       34659 :   if( getPntrToComponent(0)->getRank()==0 ) {
     227       26539 :     std::vector<double> buffer;
     228       26539 :     std::vector<double> vals( input.nscalars );
     229       26539 :     std::vector<double> deriv( input_buffer.size(), 0 );
     230       26539 :     auto output = ParallelActionsOutput::create( input.nscalars, vals.data(), input_buffer.size(), deriv.data(), 0, buffer.data() );
     231       26539 :     performTask( 0, taskmanager.getActionInput(), input, output );
     232             : 
     233       26539 :     getPntrToComponent(0)->set( vals[0] );
     234       26539 :     if( getNumberOfComponents()==2 ) {
     235       26518 :       Value* mydisp = getPntrToComponent(1);
     236     1070782 :       for(unsigned i=0; i<mydisp->getNumberOfStoredValues(); i++) {
     237     1044264 :         mydisp->set( i, vals[i+1] );
     238             :       }
     239             :     }
     240       26539 :     if( doNotCalculateDerivatives() ) {
     241             :       return;
     242             :     }
     243             : 
     244         612 :     Value* myval = getPntrToComponent(0);
     245       41772 :     for(unsigned i=0; i<deriv.size(); i++) {
     246       41160 :       myval->setDerivative( i, deriv[i] );
     247             :     }
     248             :   } else {
     249        8120 :     taskmanager.runAllTasks();
     250             :   }
     251             : }
     252             : 
     253       34639 : void RMSDVector::apply() {
     254       34639 :   if( doNotCalculateDerivatives() ) {
     255             :     return;
     256             :   }
     257             : 
     258        5530 :   if( getPntrToComponent(0)->getRank()==0 ) {
     259         612 :     Value* myval = getPntrToComponent(0);
     260         612 :     std::vector<double> deriv( input_buffer.size(), 0 );
     261       41772 :     for(unsigned i=0; i<deriv.size(); ++i) {
     262       41160 :       deriv[i] = myval->getDerivative( i );
     263             :     }
     264             : 
     265         612 :     std::vector<double> f( input.nscalars );
     266         612 :     const Value* disp = getConstPntrToComponent(1);
     267         612 :     f[0] = getConstPntrToComponent(0)->getForce(0);
     268       21192 :     for(unsigned i=0; i<disp->getNumberOfStoredValues(); ++i) {
     269       20580 :       f[i+1] = disp->getForce( i );
     270             :     }
     271             : 
     272         612 :     std::vector<double> newForces( getNumberOfDerivatives(), 0 );
     273         612 :     gatherForces( 0, taskmanager.getActionInput(), input, View<const double>( f.data(), input.nscalars ), deriv, newForces );
     274             : 
     275         612 :     unsigned ss=0;
     276         612 :     addForcesOnArguments( 0, newForces, ss  );
     277             :   } else {
     278        4918 :     ActionWithVector::apply();
     279             :   }
     280             : }
     281             : 
     282      429178 : void RMSDVector::getPositionsFromInputData( const ParallelActionsInput& input, std::vector<Vector>& pos ) {
     283      429178 :   View2D<const double> argpos( input.inputdata, 3, pos.size() );
     284             :   // At some stage it would be good to get rid of this and to be able to pass the 2D view directly to the methods in RMSD
     285     6010934 :   for(unsigned i=0; i<pos.size(); ++i) {
     286     5581756 :     pos[i][0] = argpos[0][i];
     287     5581756 :     pos[i][1] = argpos[1][i];
     288     5581756 :     pos[i][2] = argpos[2][i];
     289             :   }
     290      429178 : }
     291             : 
     292      212226 : bool RMSDVector::derivativesComputedInGatherForces( const RMSDVectorData& actiondata ) {
     293      212226 :   return actiondata.displacement && actiondata.type!="SIMPLE" && actiondata.myrmsd.size()>1;
     294             : }
     295             : 
     296      405655 : void RMSDVector::performTask( std::size_t task_index,
     297             :                               const RMSDVectorData& actiondata,
     298             :                               ParallelActionsInput& input,
     299             :                               ParallelActionsOutput& output ) {
     300             :   std::size_t natoms = actiondata.align.size();
     301      405655 :   std::vector<Vector> der(natoms), pos(natoms);
     302      405655 :   getPositionsFromInputData( input, pos );
     303      405655 :   if( actiondata.displacement && actiondata.type=="SIMPLE" ) {
     304       47569 :     std::vector<Vector> direction( natoms );
     305       47569 :     output.values[0] = actiondata.myrmsd[task_index].simpleAlignment( actiondata.align, actiondata.displace, View(pos.data(),pos.size()), actiondata.myrmsd[task_index].getReference(), der, direction, actiondata.squared );
     306      664946 :     for(unsigned i=0; i<direction.size(); ++i) {
     307      617377 :       output.values[1+i] = direction[i][0];
     308      617377 :       output.values[1+natoms+i] = direction[i][1];
     309      617377 :       output.values[1+2*natoms+i] = direction[i][2];
     310             :     }
     311      405655 :   } else if( actiondata.displacement ) {
     312      146217 :     if( !input.noderiv && derivativesComputedInGatherForces( actiondata ) ) {
     313             :       // This is activated on the backwards loop to ensure that we don't calculate the RMSD twice during that backwards loop
     314             :       // The derivatives of the distance are then taken from the alignment that is recomputed in gatherForces
     315       22944 :       return;
     316             :     }
     317      123273 :     Tensor rot;
     318      123273 :     std::vector<Vector> direction( natoms );
     319      123273 :     Matrix<std::vector<Vector> > DRotDPos(3,3);
     320      123273 :     std::vector<Vector> centeredpos( natoms ), centeredreference( natoms );
     321      123273 :     output.values[0] = actiondata.myrmsd[task_index].calc_PCAelements( pos, der, rot, DRotDPos, direction, centeredpos, centeredreference, actiondata.squared );
     322      123273 :     std::vector<Vector> ref( actiondata.myrmsd[task_index].getReference() );
     323     1730184 :     for(unsigned i=0; i<direction.size(); ++i) {
     324     1606911 :       output.values[1+i] = actiondata.sqrtdisplace[i]*( direction[i][0] - ref[i][0] );
     325     1606911 :       output.values[1+natoms+i] = actiondata.sqrtdisplace[i]*( direction[i][1] - ref[i][1] );
     326     1606911 :       output.values[1+2*natoms+i] = actiondata.sqrtdisplace[i]*( direction[i][2] - ref[i][2] );
     327             :     }
     328      246546 :   } else {
     329      211869 :     output.values[0] = actiondata.myrmsd[task_index].calculate( pos, der, actiondata.squared );
     330             :   }
     331      382711 :   if( !input.noderiv ) {
     332             :     // This also could be removed if someone was willing to do some additional work
     333     2018305 :     for(unsigned i=0; i<natoms; ++i) {
     334     1874376 :       output.derivatives[i] = der[i][0];
     335     1874376 :       output.derivatives[natoms+i] = der[i][1];
     336     1874376 :       output.derivatives[2*natoms+i] = der[i][2];
     337             :     }
     338             :   }
     339      405655 : }
     340             : 
     341        8120 : void RMSDVector::transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<double>& stash ) {
     342        8120 :   if( getNumberOfComponents()==1 ) {
     343        4368 :     ActionWithVector::transferStashToValues( partialTaskList, stash );
     344        4368 :     return;
     345             :   }
     346        3752 :   Value* dist = getPntrToComponent(0);
     347        3752 :   Value* disp = getPntrToComponent(1);
     348        3752 :   std::size_t ss = disp->getShape()[1];
     349      125144 :   for(unsigned i=0; i<partialTaskList.size(); ++i) {
     350      121392 :     dist->set(i,stash[(1+ss)*partialTaskList[i]]);
     351     4855644 :     for(unsigned j=0; j<ss; ++j) {
     352     4734252 :       disp->set( ss*partialTaskList[i] + j, stash[(1+ss)*partialTaskList[i]+1+j] );
     353             :     }
     354             :   }
     355             : }
     356             : 
     357        4918 : void RMSDVector::transferForcesToStash( const std::vector<unsigned>& partialTaskList, std::vector<double>& stash ) const {
     358        4918 :   if( getNumberOfComponents()==1 ) {
     359        3822 :     ActionWithVector::transferForcesToStash( partialTaskList, stash );
     360        3822 :     return;
     361             :   }
     362        1096 :   const Value* dist = getConstPntrToComponent(0);
     363        1096 :   const Value* disp = getConstPntrToComponent(1);
     364        1096 :   std::size_t ss = disp->getShape()[1];
     365       46972 :   for(unsigned i=0; i<partialTaskList.size(); ++i) {
     366       45876 :     stash[(1+ss)*partialTaskList[i]] = dist->getForce(i);
     367     1835004 :     for(unsigned j=0; j<ss; ++j) {
     368     1789128 :       stash[(1+ss)*partialTaskList[i]+1+j] = disp->getForce( ss*partialTaskList[i] + j );
     369             :     }
     370             :   }
     371             : }
     372             : 
     373      140946 : void RMSDVector::gatherForces( std::size_t task_index,
     374             :                                const RMSDVectorData& actiondata,
     375             :                                const ParallelActionsInput& input,
     376             :                                View<const double> f,
     377             :                                const std::vector<double>& deriv,
     378             :                                std::vector<double>& outforces ) {
     379             :   std::size_t natoms = actiondata.align.size();
     380      140946 :   if( actiondata.displacement && actiondata.type=="SIMPLE" ) {
     381             :     Vector comforce;
     382             :     comforce.zero();
     383      321235 :     for(unsigned i=0; i<natoms; ++i) {
     384      298270 :       comforce[0] += actiondata.align[i]*f[1+i];
     385      298270 :       comforce[1] += actiondata.align[i]*f[1+natoms+i];
     386      298270 :       comforce[2] += actiondata.align[i]*f[1+2*natoms+i];
     387             :     }
     388      321235 :     for(unsigned i=0; i<natoms; ++i) {
     389      298270 :       outforces[i] += f[1+i] - comforce[0];
     390      298270 :       outforces[natoms+i] += f[1+natoms+i] - comforce[1];
     391      298270 :       outforces[2*natoms+i] += f[1+2*natoms+i] - comforce[2];
     392             :     }
     393      117981 :   } else if( actiondata.displacement ) {
     394       23523 :     Tensor rot;
     395       23523 :     std::vector<Vector> der(natoms), pos(natoms);
     396       23523 :     getPositionsFromInputData( input, pos );
     397       23523 :     std::vector<Vector> direction( natoms );
     398       23523 :     Matrix<std::vector<Vector> > DRotDPos(3,3);
     399       23523 :     std::vector<Vector> centeredpos( natoms ), centeredreference( natoms );
     400       23523 :     /*double rmsd = */ actiondata.myrmsd[task_index].calc_PCAelements( pos,
     401             :         der,
     402             :         rot,
     403             :         DRotDPos,
     404             :         direction,
     405             :         centeredpos,
     406             :         centeredreference,
     407       23523 :         actiondata.squared );
     408       23523 :     Tensor trot=rot.transpose();
     409       23523 :     double prefactor = 1 / static_cast<double>( natoms );
     410             :     Vector v1;
     411             :     v1.zero();
     412      328489 :     for(unsigned n=0; n<natoms; n++) {
     413      304966 :       v1+=prefactor*matmul(trot, Vector(f[1+n],f[1+natoms+n],f[1+2*natoms+n]) );
     414             :     }
     415      328489 :     for(unsigned n=0; n<natoms; n++) {
     416      304966 :       Vector ff(f[1+n],f[1+natoms+n],f[1+2*natoms+n]);
     417      304966 :       Vector oforce = actiondata.sqrtdisplace[n]*( matmul(trot,ff) - v1 );
     418      304966 :       outforces[n] += oforce[0];
     419      304966 :       outforces[natoms + n] += oforce[1];
     420      304966 :       outforces[2*natoms + n] += oforce[2];
     421             :     }
     422       94092 :     for(unsigned a=0; a<3; a++) {
     423      282276 :       for(unsigned b=0; b<3; b++) {
     424             :         double tmp1=0.;
     425     2956401 :         for(unsigned m=0; m<natoms; m++) {
     426     2744694 :           tmp1+=centeredpos[m][b]*f[1+a*natoms + m];
     427             :         }
     428     2956401 :         for(unsigned i=0; i<natoms; i++) {
     429     2744694 :           outforces[i] += actiondata.sqrtdisplace[i]*tmp1*DRotDPos[a][b][i][0];
     430     2744694 :           outforces[natoms + i] += actiondata.sqrtdisplace[i]*tmp1*DRotDPos[a][b][i][1];
     431     2744694 :           outforces[2*natoms + i] += actiondata.sqrtdisplace[i]*tmp1*DRotDPos[a][b][i][2];
     432             :         }
     433             :       }
     434             :     }
     435       23523 :     if( derivativesComputedInGatherForces( actiondata ) ) {
     436             :       // Forces on the distance, using the derivatives that were calculated by the alignment above
     437      321204 :       for(unsigned i=0; i<natoms; ++i) {
     438      298260 :         outforces[i] += f[0]*der[i][0];
     439      298260 :         outforces[natoms + i] += f[0]*der[i][1];
     440      298260 :         outforces[2*natoms + i] += f[0]*der[i][2];
     441             :       }
     442             :     }
     443       47046 :   }
     444      140946 :   if( derivativesComputedInGatherForces( actiondata ) ) {
     445             :     // The derivatives of the distance were not calculated in performTask (deriv is zero here)
     446             :     // so the forces on the distance are applied in the branch above
     447             :     return;
     448             :   }
     449      118002 :   double ff = f[0];
     450   193082802 :   for(unsigned j=0; j<deriv.size(); ++j ) {
     451   192964800 :     outforces[j] += ff*deriv[j];
     452             :   }
     453             : }
     454             : 
     455        4918 : void RMSDVector::applyNonZeroRankForces( std::vector<double>& outforces ) {
     456             :   // Get the list of active tasks
     457        4918 :   std::vector<unsigned> & partialTaskList( getListOfActiveTasks() );
     458        4918 :   unsigned nactive_tasks=partialTaskList.size();
     459             :   // Clear force buffer
     460        4918 :   outforces.assign( outforces.size(), 0.0 );
     461             :   // Make sure that forces are calculated
     462        4918 :   input.noderiv = false;
     463             :   // Get the forces
     464        4918 :   transferForcesToStash( partialTaskList, force_stash );
     465             :   // Get the MPI details
     466        4918 :   unsigned stride=comm.Get_size();
     467        4918 :   unsigned rank=comm.Get_rank();
     468        4918 :   if( taskmanager.runInSerial() ) {
     469             :     stride=1;
     470             :     rank=0;
     471             :   }
     472             : 
     473             :   // Get number of threads for OpenMP
     474        4918 :   unsigned nt=OpenMP::getNumThreads();
     475        4918 :   if( nt*stride*10>nactive_tasks ) {
     476         550 :     nt=nactive_tasks/stride/10;
     477             :   }
     478             :   if( nt==0 ) {
     479             :     nt=1;
     480             :   }
     481             : 
     482        4918 :   #pragma omp parallel num_threads(nt)
     483             :   {
     484             :     unsigned nderivatives_per_component = getNumberOfDerivatives();
     485             :     std::vector<double> omp_forces( nderivatives_per_component, 0.0 );
     486             :     std::vector<double> buffer( 0 );
     487             :     std::vector<double> fake_vals( input.nscalars );
     488             :     std::vector<double> derivatives( nderivatives_per_component );
     489             :     #pragma omp for nowait
     490             :     for(unsigned i=rank; i<nactive_tasks; i+=stride) {
     491             :       std::size_t task_index = partialTaskList[i];
     492             :       auto myout = ParallelActionsOutput::create( input.nscalars,
     493             :                    fake_vals.data(),
     494             :                    derivatives.size(),
     495             :                    derivatives.data(),
     496             :                    0,
     497             :                    buffer.data() );
     498             :       // Calculate the stuff in the loop for this action
     499             :       performTask( task_index, taskmanager.getActionInput(), input, myout );
     500             : 
     501             :       // Gather the forces from the values
     502             :       gatherForces( task_index,
     503             :                     taskmanager.getActionInput(),
     504             :                     input,
     505             :                     View<const double>( force_stash.data()+input.nscalars*task_index, input.nscalars ),
     506             :                     derivatives,
     507             :                     omp_forces );
     508             :     }
     509             : 
     510             :     #pragma omp critical
     511             :     for(unsigned i=0; i<outforces.size(); ++i) {
     512             :       outforces[i] += omp_forces[i];
     513             :     }
     514             :   }
     515             :   // MPI Gather everything (this must be extended to the gpu thing, after makning it mpi-aware)
     516        4918 :   if( !taskmanager.runInSerial() ) {
     517        4918 :     comm.Sum( outforces );
     518             :   }
     519        4918 : }
     520             : 
     521             : }
     522             : }
     523             : 
     524             : 
     525             :