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 "ActionAtomistic.h"
      23             : #include "PlumedMain.h"
      24             : #include "ActionSet.h"
      25             : #include "GenericMolInfo.h"
      26             : #include "PbcAction.h"
      27             : #include <vector>
      28             : #include <string>
      29             : #include "ActionWithValue.h"
      30             : #include "ActionShortcut.h"
      31             : #include "Group.h"
      32             : #include "ActionWithVirtualAtom.h"
      33             : #include "tools/Exception.h"
      34             : #include "tools/Pbc.h"
      35             : #include "tools/PDB.h"
      36             : 
      37             : namespace PLMD {
      38             : 
      39       18306 : ActionAtomistic::~ActionAtomistic() {
      40             : // empty destructor to delete unique_ptr
      41       18306 : }
      42             : 
      43       18306 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
      44             :   Action(ao),
      45       18306 :   unique_local_needs_update(true),
      46       18306 :   boxValue(NULL),
      47       18306 :   lockRequestAtoms(false),
      48       18306 :   donotretrieve(false),
      49       18306 :   donotforce(false),
      50       18306 :   massesWereSet(false),
      51       18306 :   chargesWereSet(false)
      52             : {
      53       18306 :   ActionWithValue* bv = plumed.getActionSet().selectWithLabel<ActionWithValue*>("Box");
      54       18306 :   if( bv ) boxValue=bv->copyOutput(0);
      55             :   // We now get all the information about atoms that are lying about
      56       18306 :   getAtomValuesFromPlumedObject( plumed, xpos, ypos, zpos, masv, chargev );
      57       18306 :   if( xpos.size()!=ypos.size() || xpos.size()!=zpos.size() || xpos.size()!=masv.size() || xpos.size()!=chargev.size() )
      58           0 :     error("mismatch between value arrays");
      59       18306 : }
      60             : 
      61       18314 : void ActionAtomistic::getAtomValuesFromPlumedObject( const PlumedMain& plumed, std::vector<Value*>& xpos, std::vector<Value*>& ypos, std::vector<Value*>& zpos, std::vector<Value*>& masv, std::vector<Value*>& chargev ) {
      62       18314 :   std::vector<ActionShortcut*> shortcuts = plumed.getActionSet().select<ActionShortcut*>(); bool foundpdb=false;
      63     6737322 :   for(const auto & ss : shortcuts ) {
      64     6719008 :     if( ss->getName()=="READMASSCHARGE" ) {
      65             :       foundpdb=true;
      66           2 :       ActionWithValue* mv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_mass");
      67           2 :       plumed_assert( mv ); masv.push_back( mv->copyOutput(0) );
      68           2 :       ActionWithValue* qv = plumed.getActionSet().selectWithLabel<ActionWithValue*>( ss->getShortcutLabel() + "_charges");
      69           2 :       plumed_assert( qv ); chargev.push_back( qv->copyOutput(0) );
      70             :     }
      71             :   }
      72       18314 :   std::vector<ActionWithValue*> vatoms = plumed.getActionSet().select<ActionWithValue*>();
      73     7539490 :   for(const auto & vv : vatoms ) {
      74     7521176 :     plumed_assert(vv); // needed for following calls, see #1046
      75     7521176 :     ActionToPutData* ap = vv->castToActionToPutData();
      76     7521176 :     if( ap ) {
      77      253584 :       if( ap->getRole()=="x" ) xpos.push_back( ap->copyOutput(0) );
      78      253584 :       if( ap->getRole()=="y" ) ypos.push_back( ap->copyOutput(0) );
      79      253584 :       if( ap->getRole()=="z" ) zpos.push_back( ap->copyOutput(0) );
      80      380348 :       if( !foundpdb && ap->getRole()=="m" ) masv.push_back( ap->copyOutput(0) );
      81      380348 :       if( !foundpdb && ap->getRole()=="q" ) chargev.push_back( ap->copyOutput(0) );
      82             :     }
      83     7521176 :     ActionWithVirtualAtom* av = vv->castToActionWithVirtualAtom();
      84     7521176 :     if( av || vv->getName()=="ARGS2VATOM" ) {
      85     6711834 :       xpos.push_back( vv->copyOutput( vv->getLabel() + ".x") );
      86     6711834 :       ypos.push_back( vv->copyOutput( vv->getLabel() + ".y") );
      87     6711834 :       zpos.push_back( vv->copyOutput( vv->getLabel() + ".z") );
      88     6711834 :       masv.push_back( vv->copyOutput( vv->getLabel() + ".mass") );
      89     6711834 :       chargev.push_back( vv->copyOutput( vv->getLabel() + ".charge") );
      90             :     }
      91             :   }
      92       18314 : }
      93             : 
      94       30874 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
      95             :   (void) keys; // avoid warning
      96       30874 : }
      97             : 
      98             : 
      99       14077 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
     100       14077 :   plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
     101       14077 :   int nat=a.size();
     102       14077 :   indexes=a;
     103       14077 :   positions.resize(nat);
     104       14077 :   forces.resize(nat);
     105       14077 :   masses.resize(nat);
     106       14077 :   charges.resize(nat);
     107       14077 :   atom_value_ind.resize( a.size() );
     108       14077 :   int n=getTotAtoms();
     109       14077 :   if(clearDep) clearDependencies();
     110       14606 :   unique.clear(); std::vector<bool> requirements( xpos.size(), false );
     111       14077 :   if( boxValue ) addDependency( boxValue->getPntrToAction() );
     112     1337488 :   for(unsigned i=0; i<indexes.size(); i++) {
     113     1323414 :     if(indexes[i].index()>=n) { std::string num; Tools::convert( indexes[i].serial(),num ); error("atom " + num + " out of range"); }
     114     1323411 :     atom_value_ind[i] = getValueIndices( indexes[i] ); requirements[atom_value_ind[i].first] = true;
     115     1323411 :     if( atom_value_ind[i].first==0 ) unique.push_back(indexes[i]);
     116       18024 :     else if( atom_value_ind[i].second>0 ) error("action atomistic is not set up to deal with multiple vectors in position input");
     117             :   }
     118             : 
     119       14076 :   atom_value_ind_grouped.clear();
     120             : 
     121       14076 :   if(atom_value_ind.size()>0) {
     122       14035 :     auto nn = atom_value_ind[0].first;
     123       14035 :     auto kk = atom_value_ind[0].second;
     124       14035 :     atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
     125       14035 :     atom_value_ind_grouped.back().second.push_back(kk);
     126             :     auto prev_nn=nn;
     127     1323404 :     for(unsigned i=1; i<atom_value_ind.size(); i++) {
     128     1309369 :       auto nn = atom_value_ind[i].first;
     129     1309369 :       auto kk = atom_value_ind[i].second;
     130     1329052 :       if(nn!=prev_nn) atom_value_ind_grouped.push_back(std::pair<std::size_t,std::vector<std::size_t>>(nn, {}));
     131     1309369 :       atom_value_ind_grouped.back().second.push_back(kk);
     132             :       prev_nn=nn;
     133             :     }
     134             :   }
     135             : 
     136             :   // Add the dependencies to the actions that we require
     137       14076 :   Tools::removeDuplicates(unique); value_depends.resize(0);
     138     6661371 :   for(unsigned i=0; i<requirements.size(); ++i ) {
     139     6647295 :     if( !requirements[i] ) continue;
     140       30783 :     value_depends.push_back( i );
     141       30782 :     addDependency( xpos[i]->getPntrToAction() );
     142       30782 :     addDependency( ypos[i]->getPntrToAction() );
     143       30782 :     addDependency( zpos[i]->getPntrToAction() );
     144       30782 :     addDependency( masv[i]->getPntrToAction() );
     145       30782 :     addDependency( chargev[i]->getPntrToAction() );
     146             :   }
     147       14076 :   unique_local_needs_update=true;
     148       14076 : }
     149             : 
     150      405876 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
     151      405876 :   pbc.apply(dlist, max_index);
     152      405876 : }
     153             : 
     154         156 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
     155         156 :   calculateAtomicNumericalDerivatives( a, 0 );
     156         156 : }
     157             : 
     158           0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
     159           0 :   pbc.setBox( newbox );
     160           0 : }
     161             : 
     162         240 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
     163         240 :   if(!a) {
     164         240 :     a=castToActionWithValue();
     165         240 :     plumed_massert(a,"only Actions with a value can be differentiated");
     166             :   }
     167             : 
     168         240 :   const size_t nval=a->getNumberOfComponents();
     169         240 :   const size_t natoms=getNumberOfAtoms();
     170         240 :   std::vector<Vector> value(nval*natoms);
     171         240 :   std::vector<Tensor> valuebox(nval);
     172         240 :   std::vector<Vector> savedPositions(natoms);
     173             :   const double delta=std::sqrt(epsilon);
     174             : 
     175       82068 :   for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
     176       61371 :       savedPositions[i][k]=positions[i][k];
     177       61371 :       positions[i][k]=positions[i][k]+delta;
     178       61371 :       a->calculate();
     179       61371 :       positions[i][k]=savedPositions[i][k];
     180      188208 :       for(int j=0; j<nval; j++) {
     181      126837 :         value[j*natoms+i][k]=a->getOutputQuantity(j);
     182             :       }
     183             :     }
     184         240 :   Tensor box(pbc.getBox());
     185        3120 :   for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
     186        2160 :       double arg0=box(i,k);
     187      186273 :       for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
     188        2160 :       box(i,k)=box(i,k)+delta;
     189        2160 :       pbc.setBox(box);
     190      186273 :       for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
     191        2160 :       a->calculate();
     192        2160 :       box(i,k)=arg0;
     193        2160 :       pbc.setBox(box);
     194      186273 :       for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
     195       14886 :       for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
     196             :     }
     197             : 
     198         240 :   a->calculate();
     199         240 :   a->clearDerivatives();
     200        1654 :   for(int j=0; j<nval; j++) {
     201        1414 :     Value* v=a->copyOutput(j);
     202        1414 :     double ref=v->get();
     203        1414 :     if(v->hasDerivatives()) {
     204       89779 :       for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
     205       66975 :           double d=(value[j*natoms+i][k]-ref)/delta;
     206       66975 :           v->addDerivative(startnum+3*i+k,d);
     207             :         }
     208         479 :       Tensor virial;
     209        6227 :       for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
     210             : // BE CAREFUL WITH NON ORTHOROMBIC CELL
     211         479 :       virial=-matmul(box.transpose(),virial);
     212        6227 :       for(int i=0; i<3; i++) for(int k=0; k<3; k++) v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
     213             :     }
     214             :   }
     215         240 : }
     216             : 
     217      105422 : bool ActionAtomistic::actionHasForces() {
     218      105422 :   ActionWithValue* av = castToActionWithValue();
     219      105422 :   if( av ) return !av->doNotCalculateDerivatives();
     220           0 :   if( indexes.size()>0 ) plumed_merror("you have to overwrite the function actionHasForce to tell plumed if you method applies forces");
     221             :   return true;
     222             : }
     223             : 
     224       13922 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
     225       13922 :   parseAtomList(key,-1,t);
     226       13921 : }
     227             : 
     228       51540 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
     229       51540 :   plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
     230             :   std::vector<std::string> strings;
     231       51540 :   if( num<0 ) {
     232       17481 :     parseVector(key,strings);
     233       17480 :     if(strings.empty()) return;
     234             :   } else {
     235       34059 :     if ( !parseNumberedVector(key,num,strings) ) return;
     236             :   }
     237       47363 :   t.resize(0); interpretAtomList( strings, xpos, this, t );
     238       51540 : }
     239             : 
     240          28 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
     241          28 :   interpretAtomList( strings, xpos, this, t );
     242          28 : }
     243             : 
     244       47955 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t) {
     245       47955 :   Tools::interpretRanges(strings);
     246     1159863 :   for(unsigned i=0; i<strings.size(); ++i) {
     247             :     AtomNumber atom;
     248     1111908 :     bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
     249     1111908 :     if(ok) t.push_back(atom);
     250             : // here we check if this is a special symbol for MOLINFO
     251     1111908 :     if( !ok && strings[i].compare(0,1,"@")==0 ) {
     252         705 :       std::string symbol=strings[i].substr(1);
     253         705 :       if(symbol=="allatoms") {
     254          24 :         auto n=0; for(unsigned i=0; i<xpos.size(); ++i) n += xpos[i]->getNumberOfValues();
     255         365 :         t.reserve(n); for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
     256             :         ok=true;
     257         695 :       } else if(symbol=="mdatoms") {
     258          14 :         const auto n=xpos[0]->getNumberOfValues();
     259          14 :         t.reserve(t.size()+n);
     260        1166 :         for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
     261             :         ok=true;
     262        1362 :       } else if(Tools::startWith(symbol,"ndx:")) {
     263          40 :         auto words=Tools::getWords(symbol.substr(4));
     264             :         std::string ndxfile,ndxgroup;
     265          20 :         if(words.size()==1) {
     266             :           ndxfile=words[0];
     267          18 :         } else if(words.size()==2) {
     268             :           ndxfile=words[0];
     269             :           ndxgroup=words[1];
     270           0 :         } else plumed_error()<<"Cannot intepret selection "<<symbol;
     271             : 
     272          38 :         if(ndxgroup.size()>0) action->log<<"  importing group '"+ndxgroup+"'";
     273           2 :         else                  action->log<<"  importing first group";
     274          20 :         action->log<<" from index file "<<ndxfile<<"\n";
     275             : 
     276          20 :         IFile ifile;
     277          20 :         ifile.open(ndxfile);
     278             :         std::string line;
     279             :         std::string groupname;
     280             :         bool firstgroup=true;
     281             :         bool groupfound=false;
     282        3849 :         while(ifile.getline(line)) {
     283        3829 :           std::vector<std::string> words=Tools::getWords(line);
     284        7765 :           if(words.size()>=3 && words[0]=="[" && words[2]=="]") {
     285         168 :             if(groupname.length()>0) firstgroup=false;
     286             :             groupname=words[1];
     287         168 :             if(groupname==ndxgroup || ndxgroup.length()==0) groupfound=true;
     288        3661 :           } else if(groupname==ndxgroup || (firstgroup && ndxgroup.length()==0)) {
     289        6186 :             for(unsigned i=0; i<words.size(); i++) {
     290        5782 :               AtomNumber at; Tools::convert(words[i],at);
     291        5782 :               t.push_back(at);
     292             :             }
     293             :           }
     294        3829 :         }
     295          20 :         if(!groupfound) plumed_error()<<"group has not been found in index file";
     296             :         ok=true;
     297          40 :       } else {
     298         661 :         auto* moldat=action->plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
     299         661 :         if( moldat ) {
     300         661 :           std::vector<AtomNumber> atom_list; moldat->interpretSymbol( symbol, atom_list );
     301         661 :           if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
     302           0 :           else { action->error(strings[i] + " is not a label plumed knows"); }
     303             :         } else {
     304           0 :           action->error("atoms specified using @ symbol but no MOLINFO was available");
     305             :         }
     306             :       }
     307             :     }
     308             : // here we check if the atom name is the name of a group
     309     1111203 :     if(!ok) {
     310       24226 :       Group* mygrp=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]);
     311       24226 :       if(mygrp) {
     312         445 :         std::vector<std::string> grp_str( mygrp->getGroupAtoms() );
     313         445 :         interpretAtomList( grp_str, xpos, action, t ); ok=true;
     314         445 :       } else {
     315       23781 :         Group* mygrp2=action->plumed.getActionSet().selectWithLabel<Group*>(strings[i]+"_grp");
     316       23781 :         if(mygrp2) {
     317         111 :           std::vector<std::string> grp_str( mygrp2->getGroupAtoms() );
     318         111 :           interpretAtomList( grp_str, xpos, action, t ); ok=true;
     319         111 :         }
     320             :       }
     321             :     }
     322             : // here we check if the atom name is the name of an added virtual atom
     323     1111352 :     if(!ok) {
     324             :       unsigned ind = 0;
     325    13269158 :       for(unsigned j=0; j<xpos.size(); ++j) {
     326    13269158 :         if( xpos[j]->getPntrToAction()->getLabel()==strings[i] ) {
     327       23670 :           t.push_back( AtomNumber::index(ind) ); ok=true; break;
     328             :         }
     329    13245488 :         ind = ind + xpos[j]->getNumberOfValues();
     330             :       }
     331             :     }
     332     1088238 :     if(!ok) action->error("it was not possible to interpret atom name " + strings[i]);
     333             :     // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
     334             :   }
     335       47955 : }
     336             : 
     337     1761306 : std::pair<std::size_t, std::size_t> ActionAtomistic::getValueIndices( const AtomNumber& i ) const {
     338     1761306 :   std::size_t valno=0, k = i.index();
     339    15486520 :   for(unsigned j=0; j<xpos.size(); ++j) {
     340    15486520 :     if( k<xpos[j]->getNumberOfValues() ) { valno=j; break; }
     341    13725214 :     k = k - xpos[j]->getNumberOfValues();
     342             :   }
     343     1761306 :   return std::pair<std::size_t, std::size_t>( valno, k );
     344             : }
     345             : 
     346      542621 : void ActionAtomistic::retrieveAtoms( const bool& force ) {
     347      542621 :   if( boxValue ) {
     348      526106 :     auto* ptr=boxValue->getPntrToAction();
     349      526106 :     plumed_assert(ptr); // needed for following calls, see #1046
     350      526106 :     PbcAction* pbca = ptr->castToPbcAction();
     351      526106 :     plumed_assert( pbca ); pbc=pbca->pbc;
     352             :   }
     353      542621 :   if( donotretrieve || indexes.size()==0 ) return;
     354      225004 :   auto * mtr=masv[0]->getPntrToAction();
     355      225004 :   plumed_assert(mtr); // needed for following calls, see #1046
     356      225004 :   ActionToPutData* mv = mtr->castToActionToPutData();
     357      225004 :   if(mv) massesWereSet=mv->hasBeenSet();
     358           2 :   else if( (masv[0]->getPntrToAction())->getName()=="CONSTANT" ) massesWereSet=true; // Read masses from PDB file
     359      225004 :   auto * ptr=chargev[0]->getPntrToAction();
     360      225004 :   plumed_assert(ptr); // needed for following calls, see #1046
     361      225004 :   ActionToPutData* cv = ptr->castToActionToPutData();
     362      225004 :   if(cv) chargesWereSet=cv->hasBeenSet();
     363           2 :   else if( (chargev[0]->getPntrToAction())->getName()=="CONSTANT" ) chargesWereSet=true; // Read masses from PDB file
     364             :   unsigned j = 0;
     365             : 
     366             : // for(const auto & a : atom_value_ind) {
     367             : //   std::size_t nn = a.first, kk = a.second;
     368             : //   positions[j][0] = xpos[nn]->data[kk];
     369             : //   positions[j][1] = ypos[nn]->data[kk];
     370             : //   positions[j][2] = zpos[nn]->data[kk];
     371             : //   charges[j] = chargev[nn]->data[kk];
     372             : //   masses[j] = masv[nn]->data[kk];
     373             : //   j++;
     374             : // }
     375             : 
     376      840649 :   for(const auto & a : atom_value_ind_grouped) {
     377      615645 :     const auto nn=a.first;
     378      615645 :     auto & xp=xpos[nn]->data;
     379      615645 :     auto & yp=ypos[nn]->data;
     380      615645 :     auto & zp=zpos[nn]->data;
     381      615645 :     auto & ch=chargev[nn]->data;
     382      615645 :     auto & ma=masv[nn]->data;
     383     6608291 :     for(const auto & kk : a.second) {
     384     5992646 :       positions[j][0] = xp[kk];
     385     5992646 :       positions[j][1] = yp[kk];
     386     5992646 :       positions[j][2] = zp[kk];
     387     5992646 :       charges[j] = ch[kk];
     388     5992646 :       masses[j] = ma[kk];
     389     5992646 :       j++;
     390             :     }
     391             :   }
     392             : 
     393             : }
     394             : 
     395      247404 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned& ind) {
     396      247404 :   if( donotforce || (indexes.size()==0 && getName()!="FIXEDATOM") ) return;
     397      486697 :   for(unsigned i=0; i<value_depends.size(); ++i) {
     398      343121 :     xpos[value_depends[i]]->hasForce = true;
     399      343121 :     ypos[value_depends[i]]->hasForce = true;
     400      343121 :     zpos[value_depends[i]]->hasForce = true;
     401             :   }
     402             : 
     403             : // for(const auto & a : atom_value_ind) {
     404             : //   plumed_dbg_massert( ind<forcesToApply.size(), "problem setting forces in " + getLabel() );
     405             : //   std::size_t nn = a.first, kk = a.second;
     406             : //   xpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
     407             : //   ypos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
     408             : //   zpos[nn]->inputForce[kk] += forcesToApply[ind]; ind++;
     409             : // }
     410             : 
     411      797085 :   for(const auto & a : atom_value_ind_grouped) {
     412      653509 :     const auto nn=a.first;
     413             :     plumed_dbg_assert(ind<forcesToApply.size()) << "problem setting forces in " << getLabel();
     414      653509 :     auto & xp=xpos[nn]->inputForce;
     415      653509 :     auto & yp=ypos[nn]->inputForce;
     416      653509 :     auto & zp=zpos[nn]->inputForce;
     417     3980595 :     for(const auto & kk : a.second) {
     418     3327086 :       xp[kk] += forcesToApply[ind]; ind++;
     419     3327086 :       yp[kk] += forcesToApply[ind]; ind++;
     420     3327086 :       zp[kk] += forcesToApply[ind]; ind++;
     421             :     }
     422             :   }
     423             : 
     424      143576 :   setForcesOnCell( forcesToApply, ind );
     425             : }
     426             : 
     427      143664 : void ActionAtomistic::setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind) {
     428      143664 :   setForcesOnCell(forcesToApply.data(),forcesToApply.size(),ind);
     429      143664 : }
     430             : 
     431      166565 : void ActionAtomistic::setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind) {
     432     1665650 :   for(unsigned i=0; i<9; ++i) {
     433             :     plumed_dbg_massert( ind<size, "problem setting forces in " + getLabel() );
     434     1499085 :     boxValue->addForce( i, forcesToApply[ind] ); ind++;
     435             :   }
     436      166565 : }
     437             : 
     438           6 : Tensor ActionAtomistic::getVirial() const {
     439          78 :   Tensor vir; for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) vir[i][j] = boxValue->getForce(3*i+j);
     440           6 :   return vir;
     441             : }
     442             : 
     443           0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
     444             : 
     445           0 :   for(unsigned j=0; j<indexes.size(); j++) {
     446           0 :     if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
     447           0 :     if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
     448           0 :     positions[j]=pdb.getPositions()[indexes[j].index()];
     449             :   }
     450           0 :   for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
     451           0 :   for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
     452           0 : }
     453             : 
     454       14340 : unsigned ActionAtomistic::getTotAtoms()const {
     455             :   unsigned natoms = 0;
     456     6661971 :   for(unsigned i=0; i<xpos.size(); ++i ) natoms += xpos[i]->getNumberOfValues();
     457       14340 :   return natoms;
     458             : }
     459             : 
     460      149867 : void ActionAtomistic::makeWhole() {
     461     1375530 :   for(unsigned j=0; j<positions.size()-1; ++j) {
     462             :     const Vector & first (positions[j]);
     463     1225663 :     Vector & second (positions[j+1]);
     464     1225663 :     second=first+pbcDistance(first,second);
     465             :   }
     466      149867 : }
     467             : 
     468        8208 : void ActionAtomistic::getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const {
     469        8208 :   std::size_t nn = atom_value_ind[ind].first;
     470        8208 :   if( nn==0 ) { gradients[indexes[ind]] += deriv; return; }
     471          39 :   xpos[nn]->passGradients( deriv[0], gradients );
     472          39 :   ypos[nn]->passGradients( deriv[1], gradients );
     473          39 :   zpos[nn]->passGradients( deriv[2], gradients );
     474             : }
     475             : 
     476      256390 : void ActionAtomistic::updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l ) {
     477      256390 :   if( useunique ) { unique_local=unique; return; }
     478             :   // Update unique local if it needs an update
     479        9154 :   unique_local_needs_update=false; unique_local.clear();
     480       75612 :   for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
     481       66458 :     if(g2l[pp->index()]>=0) unique_local.push_back(*pp); // already sorted
     482             :   }
     483             : }
     484             : 
     485             : }