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 "AdjacencyMatrixBase.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "tools/Communicator.h"
      25             : #include "tools/OpenMP.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace adjmat {
      29             : 
      30         999 : void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
      31         999 :   ActionWithMatrix::registerKeywords( keys ); keys.remove("ARG");
      32        1998 :   keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
      33        1998 :   keys.add("atoms","GROUPA","");
      34        1998 :   keys.add("atoms","GROUPB","");
      35        1998 :   keys.add("atoms-2","ATOMS","the atoms for which you would like to calculate the adjacency matrix. This is a depracated syntax that is equivalent to GROUP.  You are strongly recommened to use GROUP instead of ATOMS.");
      36        1998 :   keys.reserve("atoms","GROUPC","");
      37        1998 :   keys.addFlag("COMPONENTS",false,"also calculate the components of the vector connecting the atoms in the contact matrix");
      38        1998 :   keys.addFlag("NOPBC",false,"don't use pbc");
      39        1998 :   keys.add("compulsory","NL_CUTOFF","0.0","The cutoff for the neighbor list.  A value of 0 means we are not using a neighbor list");
      40        1998 :   keys.add("compulsory","NL_STRIDE","1","The frequency with which we are updating the atoms in the neighbor list");
      41        1998 :   keys.addOutputComponent("w","COMPONENTS","a matrix containing the weights for the bonds between each pair of atoms");
      42        1998 :   keys.addOutputComponent("x","COMPONENTS","the projection of the bond on the x axis");
      43        1998 :   keys.addOutputComponent("y","COMPONENTS","the projection of the bond on the y axis");
      44        1998 :   keys.addOutputComponent("z","COMPONENTS","the projection of the bond on the z axis");
      45         999 :   keys.setValueDescription("a matrix containing the weights for the bonds between each pair of atoms");
      46         999 : }
      47             : 
      48         326 : AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
      49             :   Action(ao),
      50             :   ActionWithMatrix(ao),
      51         326 :   read_one_group(false),
      52         326 :   neighbour_list_updated(false),
      53         326 :   linkcells(comm),
      54         326 :   threecells(comm),
      55         326 :   maxcol(0),
      56         652 :   natoms_per_list(0)
      57             : {
      58         652 :   std::vector<unsigned> shape(2); std::vector<AtomNumber> t; parseAtomList("GROUP", t );
      59         326 :   if( t.size()==0 ) {
      60         272 :     parseAtomList("ATOMS", t);
      61         136 :     if( t.size()>0 ) warning("using depracated syntax for contact matrix.  You are strongly recommended to use GROUP instead of ATOMS");
      62             :   }
      63             : 
      64         326 :   if( t.size()==0 ) {
      65         272 :     std::vector<AtomNumber> ta; parseAtomList("GROUPA",ta);
      66         136 :     if( ta.size()==0 && getName()=="HBOND_MATRIX") parseAtomList("DONORS",ta);
      67         272 :     std::vector<AtomNumber> tb; parseAtomList("GROUPB",tb);
      68         136 :     if( tb.size()==0 && getName()=="HBOND_MATRIX") parseAtomList("ACCEPTORS",tb);
      69         136 :     if( ta.size()==0 || tb.size()==0 ) error("no atoms have been specified in input");
      70             : 
      71             :     // Create list of tasks
      72         136 :     log.printf("  atoms in GROUPA ");
      73        2270 :     for(unsigned i=0; i<ta.size(); ++i) { log.printf("%d ", ta[i].serial()); t.push_back(ta[i]); }
      74         136 :     log.printf("\n");
      75         136 :     log.printf("  atoms in GROUPB "); ablocks.resize( tb.size() ); unsigned n=0;
      76       13747 :     for(unsigned i=0; i<tb.size(); ++i) {
      77       13611 :       log.printf("%d ", tb[i].serial());
      78       13611 :       ablocks[i]=ta.size()+n; t.push_back(tb[i]); n++;
      79             :     }
      80         136 :     log.printf("\n"); shape[0]=ta.size(); shape[1]=tb.size();
      81             :   } else {
      82             :     // Create list of tasks
      83         190 :     log.printf("  atoms in GROUP "); ablocks.resize( t.size() );
      84      117830 :     for(unsigned i=0; i<t.size(); ++i) { log.printf("%d ", t[i].serial()); ablocks[i]=i; }
      85         190 :     log.printf("\n"); shape[0]=shape[1]=t.size(); read_one_group=true;
      86             :   }
      87         652 :   if( keywords.exists("GROUPC") ) {
      88          12 :     std::vector<AtomNumber> tc; parseAtomList("GROUPC",tc);
      89           6 :     if( tc.size()==0 ) error("no atoms in GROUPC specified");
      90           6 :     log.printf("  atoms in GROUPC "); setupThirdAtomBlock( tc, t );
      91         640 :   } else if( keywords.exists("BRIDGING_ATOMS") ) {
      92          16 :     std::vector<AtomNumber> tc; parseAtomList("BRIDGING_ATOMS",tc);
      93           8 :     if( tc.size()==0 ) error("no BRIDGING_ATOMS specified");
      94           8 :     log.printf("  bridging atoms are "); setupThirdAtomBlock( tc, t );
      95         624 :   } else if( keywords.exists("HYDROGENS") ) {
      96           4 :     std::vector<AtomNumber> tc; parseAtomList("HYDROGENS",tc);
      97           2 :     if( tc.size()==0 ) error("no HYDROGEN atoms specified");
      98           2 :     log.printf("  hydrogen atoms are "); setupThirdAtomBlock( tc, t );
      99         620 :   } else if( keywords.exists("BACKGROUND_ATOMS") ) {
     100          16 :     std::vector<AtomNumber> tc; parseAtomList("BACKGROUND_ATOMS",tc);
     101           8 :     if( tc.size()==0 ) error("no ATOMS atoms specified");
     102           8 :     log.printf("  atoms for background density are "); setupThirdAtomBlock( tc, t );
     103             :   }
     104             :   // Request the atoms from the ActionAtomistic
     105         652 :   requestAtoms( t ); parseFlag("COMPONENTS",components); parseFlag("NOPBC",nopbc);
     106         326 :   if( !components ) { addValue( shape ); setNotPeriodic(); }
     107         249 :   else { addComponent( "w", shape ); componentIsNotPeriodic("w"); }
     108         326 :   getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
     109             :   // Stuff for neighbor list
     110         652 :   parse("NL_CUTOFF",nl_cut); nl_cut2=nl_cut*nl_cut; parse("NL_STRIDE",nl_stride);
     111         326 :   if( nl_cut==0 && nl_stride>1 ) error("NL_CUTOFF must be set if NL_STRIDE is set greater than 1");
     112         326 :   if( nl_cut>0 ) log.printf("  using neighbor list with cutoff %f.  List is updated every %u steps.\n",nl_cut,nl_stride);
     113             : 
     114         326 :   if( components ) {
     115         166 :     addComponent( "x", shape ); componentIsNotPeriodic("x");
     116         166 :     addComponent( "y", shape ); componentIsNotPeriodic("y");
     117         249 :     addComponent( "z", shape ); componentIsNotPeriodic("z");
     118             :   }
     119         652 :   log<<"  Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
     120         326 : }
     121             : 
     122       11081 : unsigned AdjacencyMatrixBase::getNumberOfDerivatives() {
     123       11081 :   return 3*getNumberOfAtoms() + 9;
     124             : }
     125             : 
     126          24 : void AdjacencyMatrixBase::setupThirdAtomBlock( const std::vector<AtomNumber>& tc, std::vector<AtomNumber>& t ) {
     127          24 :   threeblocks.resize( tc.size() ); unsigned base=t.size();
     128       53168 :   for(unsigned i=0; i<tc.size(); ++i) { log.printf("%d ", tc[i].serial()); t.push_back(tc[i]); threeblocks[i]=base+i; }
     129          24 :   log.printf("\n");
     130          24 : }
     131             : 
     132         326 : void AdjacencyMatrixBase::setLinkCellCutoff( const bool& symmetric, const double& lcut, double tcut ) {
     133         326 :   if( read_one_group && symmetric ) getPntrToComponent(0)->setSymmetric( true );
     134         326 :   if( nl_cut>0 && lcut>nl_cut ) error("D_MAX for switching functions should be shorter than neighbor list cutoff");
     135             : 
     136         326 :   if( tcut<0 ) tcut=lcut;
     137         326 :   if( nl_cut>0 ) linkcells.setCutoff( nl_cut ); else linkcells.setCutoff( lcut );
     138         326 :   if( linkcells.getCutoff()<std::numeric_limits<double>::max() ) log.printf("  set link cell cutoff to %f \n", linkcells.getCutoff() );
     139         326 :   threecells.setCutoff( tcut );
     140         326 : }
     141             : 
     142       11700 : void AdjacencyMatrixBase::prepare() {
     143       11700 :   ActionWithVector::prepare(); neighbour_list_updated=false;
     144       11700 : }
     145             : 
     146       11750 : void AdjacencyMatrixBase::updateNeighbourList() {
     147       11750 :   neighbour_list_updated=true;
     148             :   // Build link cells here so that this is done in stream if it needed in stream
     149       11750 :   if( getStep()%nl_stride==0 ) {
     150             :     // Build the link cells
     151       11750 :     std::vector<Vector> ltmp_pos( ablocks.size() );
     152     1194628 :     for(unsigned i=0; i<ablocks.size(); ++i) ltmp_pos[i]=ActionAtomistic::getPosition( ablocks[i] );
     153       11750 :     linkcells.buildCellLists( ltmp_pos, ablocks, getPbc() );
     154             :     // This ensures the link cell does not get too big.  We find the cell that contains the maximum number of atoms and multiply this by 27.
     155             :     // In this way we ensure that the neighbour list doesn't get too big.  Also this number should always be large enough
     156       11750 :     natoms_per_list = 27*linkcells.getMaxInCell();
     157       11750 :     nlist.resize( getConstPntrToComponent(0)->getShape()[0]*( 2 + natoms_per_list ) );
     158             :     // Set the number of neighbors to zero for all ranks
     159       11750 :     nlist.assign(nlist.size(),0);
     160             :     // Now get stuff to do parallel implementation
     161       11750 :     unsigned stride=comm.Get_size(); unsigned rank=comm.Get_rank();
     162       11750 :     if( runInSerial() ) { stride=1; rank=0; }
     163       11750 :     unsigned nt=OpenMP::getNumThreads();
     164       11750 :     if( nt*stride*10>getConstPntrToComponent(0)->getShape()[0] ) nt=getConstPntrToComponent(0)->getShape()[0]/stride/10;
     165             :     if( nt==0 ) nt=1;
     166             :     // Create a vector from the input set of tasks
     167       11750 :     std::vector<unsigned> & pTaskList( getListOfActiveTasks(this) );
     168             : 
     169       11750 :     #pragma omp parallel num_threads(nt)
     170             :     {
     171             :       // Get the number of tasks we have to deal with
     172             :       unsigned ntasks=getConstPntrToComponent(0)->getShape()[0];
     173             :       if( nl_stride==1 ) ntasks=pTaskList.size();
     174             :       // Build a tempory nlist so we can do omp parallelism
     175             :       std::vector<unsigned> omp_nlist;
     176             :       if( nt>1 ) omp_nlist.resize( nlist.size(), 0 );
     177             :       // Now run over all atoms and construct the link cells
     178             :       std::vector<Vector> t_atoms( 1+ablocks.size() );
     179             :       std::vector<unsigned> indices( 1+ablocks.size() ), cells_required( linkcells.getNumberOfCells() );
     180             :       #pragma omp for nowait
     181             :       for(unsigned i=rank; i<ntasks; i+=stride) {
     182             :         // Retrieve cells required from link cells - for matrix blocks
     183             :         unsigned ncells_required=0;
     184             :         linkcells.addRequiredCells( linkcells.findMyCell( ActionAtomistic::getPosition(pTaskList[i]) ), ncells_required, cells_required );
     185             :         // Now get the indices of the atoms in the link cells positions
     186             :         unsigned natoms=1; indices[0]=pTaskList[i];
     187             :         linkcells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
     188             :         if( nl_stride==1 ) {
     189             :           if( nt>1 ) omp_nlist[indices[0]]=0; else nlist[indices[0]] = 0;
     190             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + indices[0]*(1+natoms_per_list);
     191             :           for(unsigned j=0; j<natoms; ++j) {
     192             :             if( nt>1 ) { omp_nlist[ lstart + omp_nlist[indices[0]] ] = indices[j]; omp_nlist[indices[0]]++; }
     193             :             else { nlist[ lstart + nlist[indices[0]] ] = indices[j]; nlist[indices[0]]++; }
     194             :           }
     195             :         } else {
     196             :           // Get the positions of all the atoms in the link cells relative to the central atom
     197             :           for(unsigned j=0; j<natoms; ++j) t_atoms[j] = ActionAtomistic::getPosition(indices[j]) - ActionAtomistic::getPosition(indices[0]);
     198             :           if( !nopbc ) pbcApply( t_atoms, natoms );
     199             :           // Now construct the neighbor list
     200             :           if( nt>1 ) omp_nlist[indices[0]] = 0; else nlist[indices[0]] = 0;
     201             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + indices[0]*(1+natoms_per_list);
     202             :           for(unsigned j=0; j<natoms; ++j) {
     203             :             double d2;
     204             :             if ( (d2=t_atoms[j][0]*t_atoms[j][0])<nl_cut2 &&
     205             :                  (d2+=t_atoms[j][1]*t_atoms[j][1])<nl_cut2 &&
     206             :                  (d2+=t_atoms[j][2]*t_atoms[j][2])<nl_cut2 ) {
     207             :               if( nt>1 ) { omp_nlist[ lstart + omp_nlist[indices[0]] ] = indices[j]; omp_nlist[indices[0]]++; }
     208             :               else { nlist[ lstart + nlist[indices[0]] ] = indices[j]; nlist[indices[0]]++; }
     209             :             }
     210             :           }
     211             : 
     212             :         }
     213             :       }
     214             :       // Gather OMP stuff
     215             :       #pragma omp critical
     216             :       if(nt>1) {
     217             :         for(unsigned i=0; i<ntasks; ++i) nlist[pTaskList[i]]+=omp_nlist[pTaskList[i]];
     218             :         for(unsigned i=0; i<ntasks; ++i) {
     219             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + pTaskList[i]*(1+natoms_per_list);
     220             :           for(unsigned j=0; j<omp_nlist[pTaskList[i]]; ++j) nlist[ lstart + j ] += omp_nlist[ lstart + j ];
     221             :         }
     222             :       }
     223             :     }
     224             :     // MPI gather
     225       11750 :     if( !runInSerial() ) comm.Sum( nlist );
     226             :   }
     227       11750 :   if( threeblocks.size()>0 ) {
     228          89 :     std::vector<Vector> ltmp_pos2( threeblocks.size() );
     229       66184 :     for(unsigned i=0; i<threeblocks.size(); ++i) {
     230       66095 :       ltmp_pos2[i]=ActionAtomistic::getPosition( threeblocks[i] );
     231             :     }
     232          89 :     threecells.buildCellLists( ltmp_pos2, threeblocks, getPbc() );
     233             :   }
     234             :   // And finally work out maximum number of columns to use
     235       11750 :   maxcol = nlist[0];
     236     1140150 :   for(unsigned i=1; i<getConstPntrToComponent(0)->getShape()[0]; ++i) {
     237     1128400 :     if( nlist[i]>maxcol ) maxcol = nlist[i];
     238             :   }
     239       11750 : }
     240             : 
     241         216 : void AdjacencyMatrixBase::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
     242         216 :   if( action==this ) return;
     243             :   // Update the neighbour list
     244          54 :   updateNeighbourList();
     245             : 
     246          54 :   unsigned nactive = atasks.size();
     247          54 :   std::vector<unsigned> indlist( 1 + ablocks.size() + threeblocks.size() );
     248        3844 :   for(unsigned i=0; i<nactive; ++i) {
     249        3790 :     unsigned num = retrieveNeighbours( atasks[i], indlist );
     250      269574 :     for(unsigned j=0; j<num; ++j) {
     251             :       bool found=false;
     252    67587417 :       for(unsigned k=0; k<atasks.size(); ++k ) {
     253    67573541 :         if( indlist[j]==atasks[k] ) { found=true; break; }
     254             :       }
     255      265784 :       if( !found ) atasks.push_back( indlist[j] );
     256             :     }
     257             :   }
     258             : }
     259             : 
     260      409666 : unsigned AdjacencyMatrixBase::retrieveNeighbours( const unsigned& current, std::vector<unsigned> & indices ) const {
     261      409666 :   unsigned natoms=nlist[current]; indices[0]=current;
     262      409666 :   unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + current*(1+natoms_per_list); plumed_dbg_assert( nlist[lstart]==current );
     263    21916923 :   for(unsigned i=1; i<nlist[current]; ++i) { indices[i] = nlist[ lstart + i ]; }
     264      409666 :   return natoms;
     265             : }
     266             : 
     267      405876 : void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<unsigned> & indices, MultiValue& myvals ) const {
     268             :   // Now retrieve bookeeping arrays
     269      405876 :   if( indices.size()!=(1+ablocks.size()+threeblocks.size()) ) indices.resize( 1+ablocks.size()+threeblocks.size() );
     270             : 
     271             :   // Now get the positions
     272      405876 :   unsigned natoms=retrieveNeighbours( current, indices );
     273             :   unsigned ntwo_atoms=natoms; myvals.setSplitIndex( ntwo_atoms );
     274             : 
     275             :   // Now retrieve everything for the third atoms
     276      405876 :   if( threeblocks.size()>0 ) {
     277         920 :     unsigned ncells_required=0; std::vector<unsigned> cells_required( threecells.getNumberOfCells() );
     278         920 :     threecells.addRequiredCells( threecells.findMyCell( ActionAtomistic::getPosition(current) ), ncells_required, cells_required );
     279         920 :     threecells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
     280             :   }
     281      405876 :   myvals.setNumberOfIndices( natoms );
     282             : 
     283             : // Apply periodic boundary conditions to atom positions
     284      405876 :   std::vector<std::vector<Vector> > & t_atoms( myvals.getFirstAtomDerivativeVector() ); if( t_atoms.size()!=1 ) t_atoms.resize(1);
     285      405876 :   if( t_atoms[0].size()<getNumberOfAtoms() ) t_atoms[0].resize( getNumberOfAtoms() );
     286    26226247 :   for(unsigned i=0; i<natoms; ++i) t_atoms[0][i] = ActionAtomistic::getPosition(indices[i]) - ActionAtomistic::getPosition(current);
     287      405876 :   if( !nopbc ) pbcApply( t_atoms[0], natoms );
     288             :   // And collect atom position data
     289             :   std::vector<Vector> & atoms( myvals.getAtomVector() );
     290      405876 :   if( atoms.size()<getNumberOfAtoms() ) atoms.resize( getNumberOfAtoms() );
     291    26226247 :   for(unsigned i=0; i<natoms; ++i) atoms[ indices[i] ] = t_atoms[0][i];
     292      405876 : }
     293             : 
     294    21245263 : void AdjacencyMatrixBase::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     295    21245263 :   Vector zero; zero.zero(); plumed_dbg_assert( index2<myvals.getAtomVector().size() );
     296    21245263 :   double weight = calculateWeight( zero, myvals.getAtomVector()[index2], myvals.getNumberOfIndices()-myvals.getSplitIndex(), myvals );
     297    21245263 :   unsigned w_ind = getConstPntrToComponent(0)->getPositionInStream(); myvals.setValue( w_ind, weight );
     298    21245263 :   if( fabs(weight)<epsilon ) { myvals.setValue( w_ind, 0 ); return; }
     299             : 
     300    14130310 :   if( !doNotCalculateDerivatives() ) {
     301             :     // Update dynamic list indices for central atom
     302     6759016 :     myvals.updateIndex( w_ind, 3*index1+0 ); myvals.updateIndex( w_ind, 3*index1+1 ); myvals.updateIndex( w_ind, 3*index1+2 );
     303             :     // Update dynamic list indices for atom forming this bond
     304     6759016 :     myvals.updateIndex( w_ind, 3*index2+0 ); myvals.updateIndex( w_ind, 3*index2+1 ); myvals.updateIndex( w_ind, 3*index2+2 );
     305             :     // Now look after all the atoms in the third block
     306             :     std::vector<unsigned> & indices( myvals.getIndices() );
     307     6872649 :     for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
     308      113633 :       myvals.updateIndex( w_ind, 3*indices[i]+0 ); myvals.updateIndex( w_ind, 3*indices[i]+1 ); myvals.updateIndex( w_ind, 3*indices[i]+2 );
     309             :     }
     310             :     // Update dynamic list indices for virial
     311    67590160 :     unsigned base = 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) myvals.updateIndex( w_ind, base+j );
     312             :     // And the indices for the derivatives of the row of the matrix
     313     6759016 :     unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     314             :     std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     315     6759016 :     matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2;
     316     6759016 :     myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
     317             :   }
     318             : 
     319             :   // Calculate the components if we need them
     320    14130310 :   if( components ) {
     321     2415123 :     unsigned x_index = getConstPntrToComponent(1)->getPositionInStream();
     322     2415123 :     unsigned y_index = getConstPntrToComponent(2)->getPositionInStream();
     323     2415123 :     unsigned z_index = getConstPntrToComponent(3)->getPositionInStream();
     324     2415123 :     Vector atom = myvals.getAtomVector()[index2];
     325     2415123 :     myvals.setValue( x_index, atom[0] ); myvals.setValue( y_index, atom[1] ); myvals.setValue( z_index, atom[2] );
     326     2415123 :     if( !doNotCalculateDerivatives() ) {
     327     1625864 :       myvals.addDerivative( x_index, 3*index1+0, -1 ); myvals.addDerivative( x_index, 3*index2+0, +1 );
     328     1625864 :       myvals.addDerivative( x_index, 3*index1+1, 0 ); myvals.addDerivative( x_index, 3*index2+1, 0 );
     329     1625864 :       myvals.addDerivative( x_index, 3*index1+2, 0 ); myvals.addDerivative( x_index, 3*index2+2, 0 );
     330     1625864 :       myvals.addDerivative( y_index, 3*index1+0, 0 ); myvals.addDerivative( y_index, 3*index2+0, 0 );
     331     1625864 :       myvals.addDerivative( y_index, 3*index1+1, -1 ); myvals.addDerivative( y_index, 3*index2+1, +1 );
     332     1625864 :       myvals.addDerivative( y_index, 3*index1+2, 0 ); myvals.addDerivative( y_index, 3*index2+2, 0 );
     333     1625864 :       myvals.addDerivative( z_index, 3*index1+0, 0 ); myvals.addDerivative( z_index, 3*index2+0, 0 );
     334     1625864 :       myvals.addDerivative( z_index, 3*index1+1, 0 ); myvals.addDerivative( z_index, 3*index2+1, 0 );
     335     1625864 :       myvals.addDerivative( z_index, 3*index1+2, -1 ); myvals.addDerivative( z_index, 3*index2+2, +1 );
     336     6503456 :       for(unsigned k=0; k<3; ++k) {
     337             :         // Update dynamic lists for central atom
     338     4877592 :         myvals.updateIndex( x_index, 3*index1+k ); myvals.updateIndex( y_index, 3*index1+k ); myvals.updateIndex( z_index, 3*index1+k );
     339             :         // Update dynamic lists for bonded atom
     340     4877592 :         myvals.updateIndex( x_index, 3*index2+k ); myvals.updateIndex( y_index, 3*index2+k ); myvals.updateIndex( z_index, 3*index2+k );
     341             :       }
     342             :       // Add derivatives of virial
     343     1625864 :       unsigned base = 3*getNumberOfAtoms();
     344             :       // Virial for x
     345     1625864 :       myvals.addDerivative( x_index, base+0, -atom[0] ); myvals.addDerivative( x_index, base+3, -atom[1] ); myvals.addDerivative( x_index, base+6, -atom[2] );
     346     1625864 :       myvals.addDerivative( x_index, base+1, 0 ); myvals.addDerivative( x_index, base+4, 0 ); myvals.addDerivative( x_index, base+7, 0 );
     347     1625864 :       myvals.addDerivative( x_index, base+2, 0 ); myvals.addDerivative( x_index, base+5, 0 ); myvals.addDerivative( x_index, base+8, 0 );
     348             :       // Virial for y
     349     1625864 :       myvals.addDerivative( y_index, base+0, 0 ); myvals.addDerivative( y_index, base+3, 0 ); myvals.addDerivative( y_index, base+6, 0 );
     350     1625864 :       myvals.addDerivative( y_index, base+1, -atom[0] ); myvals.addDerivative( y_index, base+4, -atom[1] ); myvals.addDerivative( y_index, base+7, -atom[2] );
     351     1625864 :       myvals.addDerivative( y_index, base+2, 0 ); myvals.addDerivative( y_index, base+5, 0 ); myvals.addDerivative( y_index, base+8, 0 );
     352             :       // Virial for z
     353     1625864 :       myvals.addDerivative( z_index, base+0, 0 ); myvals.addDerivative( z_index, base+3, 0 ); myvals.addDerivative( z_index, base+6, 0 );
     354     1625864 :       myvals.addDerivative( z_index, base+1, 0 ); myvals.addDerivative( z_index, base+4, 0 ); myvals.addDerivative( z_index, base+7, 0 );
     355     1625864 :       myvals.addDerivative( z_index, base+2, -atom[0] ); myvals.addDerivative( z_index, base+5, -atom[1] ); myvals.addDerivative( z_index, base+8, -atom[2] );
     356    16258640 :       for(unsigned k=0; k<9; ++k) { myvals.updateIndex( x_index, base+k ); myvals.updateIndex( y_index, base+k ); myvals.updateIndex( z_index, base+k ); }
     357     6503456 :       for(unsigned k=1; k<4; ++k) {
     358     4877592 :         unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     359             :         std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     360     4877592 :         matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2;
     361     4877592 :         myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
     362             :       }
     363             :     }
     364             :   }
     365             : }
     366             : 
     367      405876 : void AdjacencyMatrixBase::runEndOfRowJobs( const unsigned& ind, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     368      405876 :   if( doNotCalculateDerivatives() ) return;
     369             : 
     370      737807 :   for(int k=0; k<getNumberOfComponents(); ++k) {
     371      411298 :     unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     372             :     std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     373      411298 :     plumed_assert( nmat_ind<matrix_indices.size() );
     374      411298 :     matrix_indices[nmat_ind+0]=3*ind+0;
     375      411298 :     matrix_indices[nmat_ind+1]=3*ind+1;
     376      411298 :     matrix_indices[nmat_ind+2]=3*ind+2;
     377      411298 :     nmat_ind+=3;
     378      423041 :     for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
     379       11743 :       matrix_indices[nmat_ind+0]=3*indices[i]+0;
     380       11743 :       matrix_indices[nmat_ind+1]=3*indices[i]+1;
     381       11743 :       matrix_indices[nmat_ind+2]=3*indices[i]+2;
     382       11743 :       nmat_ind+=3;
     383             :     }
     384      411298 :     unsigned virbase = 3*getNumberOfAtoms();
     385     4112980 :     for(unsigned i=0; i<9; ++i) matrix_indices[nmat_ind+i]=virbase+i;
     386      411298 :     nmat_ind+=9; plumed_dbg_massert( nmat_ind<=3*getNumberOfAtoms() + 9, "found too many derivatives in " + getLabel() );
     387             :     myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     388             :   }
     389             : }
     390             : 
     391             : }
     392             : }