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 );
      32         999 :   keys.add("atoms","GROUP","the atoms for which you would like to calculate the adjacency matrix");
      33         999 :   keys.add("atoms","GROUPA","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB");
      34         999 :   keys.add("atoms","GROUPB","when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA");
      35         999 :   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         999 :   keys.reserve("atoms","GROUPC","a group of atoms that must be summed over when calculating each element of the adjacency matrix");
      37         999 :   keys.addFlag("COMPONENTS",false,"also calculate the components of the vectors connecting the atoms in the contact matrix");
      38         999 :   keys.addFlag("NOPBC",false,"don't use pbc");
      39         999 :   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         999 :   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","matrix","a matrix containing the weights for the bonds between each pair of atoms");
      42        1998 :   keys.addOutputComponent("x","COMPONENTS","matrix","the projection of the bond on the x axis");
      43        1998 :   keys.addOutputComponent("y","COMPONENTS","matrix","the projection of the bond on the y axis");
      44        1998 :   keys.addOutputComponent("z","COMPONENTS","matrix","the projection of the bond on the z axis");
      45        1998 :   keys.setValueDescription("matrix","a matrix containing the weights for the bonds between each pair of atoms");
      46         999 :   keys.addDOI("10.1021/acs.jctc.6b01073");
      47         999 : }
      48             : 
      49         326 : AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
      50             :   Action(ao),
      51             :   ActionWithMatrix(ao),
      52         326 :   read_one_group(false),
      53         326 :   neighbour_list_updated(false),
      54         326 :   linkcells(comm),
      55         326 :   threecells(comm),
      56         326 :   maxcol(0),
      57         652 :   natoms_per_list(0) {
      58         326 :   std::vector<unsigned> shape(2);
      59             :   std::vector<AtomNumber> t;
      60         652 :   parseAtomList("GROUP", t );
      61         326 :   if( t.size()==0 ) {
      62         272 :     parseAtomList("ATOMS", t);
      63         136 :     if( t.size()>0 ) {
      64           0 :       warning("using depracated syntax for contact matrix.  You are strongly recommended to use GROUP instead of ATOMS");
      65             :     }
      66             :   }
      67             : 
      68         326 :   if( t.size()==0 ) {
      69             :     std::vector<AtomNumber> ta;
      70         272 :     parseAtomList("GROUPA",ta);
      71         136 :     if( ta.size()==0 && getName()=="HBOND_MATRIX") {
      72           0 :       parseAtomList("DONORS",ta);
      73             :     }
      74             :     std::vector<AtomNumber> tb;
      75         272 :     parseAtomList("GROUPB",tb);
      76         136 :     if( tb.size()==0 && getName()=="HBOND_MATRIX") {
      77           0 :       parseAtomList("ACCEPTORS",tb);
      78             :     }
      79         136 :     if( ta.size()==0 || tb.size()==0 ) {
      80           0 :       error("no atoms have been specified in input");
      81             :     }
      82             : 
      83             :     // Create list of tasks
      84         136 :     log.printf("  atoms in GROUPA ");
      85        2270 :     for(unsigned i=0; i<ta.size(); ++i) {
      86        2134 :       log.printf("%d ", ta[i].serial());
      87        2134 :       t.push_back(ta[i]);
      88             :     }
      89         136 :     log.printf("\n");
      90         136 :     log.printf("  atoms in GROUPB ");
      91         136 :     ablocks.resize( tb.size() );
      92             :     unsigned n=0;
      93       13747 :     for(unsigned i=0; i<tb.size(); ++i) {
      94       13611 :       log.printf("%d ", tb[i].serial());
      95       13611 :       ablocks[i]=ta.size()+n;
      96       13611 :       t.push_back(tb[i]);
      97       13611 :       n++;
      98             :     }
      99         136 :     log.printf("\n");
     100         136 :     shape[0]=ta.size();
     101         136 :     shape[1]=tb.size();
     102             :   } else {
     103             :     // Create list of tasks
     104         190 :     log.printf("  atoms in GROUP ");
     105         190 :     ablocks.resize( t.size() );
     106      117830 :     for(unsigned i=0; i<t.size(); ++i) {
     107      117640 :       log.printf("%d ", t[i].serial());
     108      117640 :       ablocks[i]=i;
     109             :     }
     110         190 :     log.printf("\n");
     111         190 :     shape[0]=shape[1]=t.size();
     112         190 :     read_one_group=true;
     113             :   }
     114         326 :   if( keywords.exists("GROUPC") ) {
     115             :     std::vector<AtomNumber> tc;
     116          12 :     parseAtomList("GROUPC",tc);
     117           6 :     if( tc.size()==0 ) {
     118           0 :       error("no atoms in GROUPC specified");
     119             :     }
     120           6 :     log.printf("  atoms in GROUPC ");
     121           6 :     setupThirdAtomBlock( tc, t );
     122         320 :   } else if( keywords.exists("BRIDGING_ATOMS") ) {
     123             :     std::vector<AtomNumber> tc;
     124          16 :     parseAtomList("BRIDGING_ATOMS",tc);
     125           8 :     if( tc.size()==0 ) {
     126           0 :       error("no BRIDGING_ATOMS specified");
     127             :     }
     128           8 :     log.printf("  bridging atoms are ");
     129           8 :     setupThirdAtomBlock( tc, t );
     130         312 :   } else if( keywords.exists("HYDROGENS") ) {
     131             :     std::vector<AtomNumber> tc;
     132           4 :     parseAtomList("HYDROGENS",tc);
     133           2 :     if( tc.size()==0 ) {
     134           0 :       error("no HYDROGEN atoms specified");
     135             :     }
     136           2 :     log.printf("  hydrogen atoms are ");
     137           2 :     setupThirdAtomBlock( tc, t );
     138         310 :   } else if( keywords.exists("BACKGROUND_ATOMS") ) {
     139             :     std::vector<AtomNumber> tc;
     140          16 :     parseAtomList("BACKGROUND_ATOMS",tc);
     141           8 :     if( tc.size()==0 ) {
     142           0 :       error("no ATOMS atoms specified");
     143             :     }
     144           8 :     log.printf("  atoms for background density are ");
     145           8 :     setupThirdAtomBlock( tc, t );
     146             :   }
     147             :   // Request the atoms from the ActionAtomistic
     148         326 :   requestAtoms( t );
     149         326 :   parseFlag("COMPONENTS",components);
     150         326 :   parseFlag("NOPBC",nopbc);
     151         326 :   if( !components ) {
     152         243 :     addValue( shape );
     153         243 :     setNotPeriodic();
     154             :   } else {
     155          83 :     addComponent( "w", shape );
     156         166 :     componentIsNotPeriodic("w");
     157             :   }
     158         326 :   getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
     159             :   // Stuff for neighbor list
     160         326 :   parse("NL_CUTOFF",nl_cut);
     161         326 :   nl_cut2=nl_cut*nl_cut;
     162         326 :   parse("NL_STRIDE",nl_stride);
     163         326 :   if( nl_cut==0 && nl_stride>1 ) {
     164           0 :     error("NL_CUTOFF must be set if NL_STRIDE is set greater than 1");
     165             :   }
     166         326 :   if( nl_cut>0 ) {
     167           0 :     log.printf("  using neighbor list with cutoff %f.  List is updated every %u steps.\n",nl_cut,nl_stride);
     168             :   }
     169             : 
     170         326 :   if( components ) {
     171          83 :     addComponent( "x", shape );
     172          83 :     componentIsNotPeriodic("x");
     173          83 :     addComponent( "y", shape );
     174          83 :     componentIsNotPeriodic("y");
     175          83 :     addComponent( "z", shape );
     176         166 :     componentIsNotPeriodic("z");
     177             :   }
     178         652 :   log<<"  Bibliography "<<plumed.cite("10.1021/acs.jctc.6b01073")<<"\n";
     179         326 : }
     180             : 
     181       11081 : unsigned AdjacencyMatrixBase::getNumberOfDerivatives() {
     182       11081 :   return 3*getNumberOfAtoms() + 9;
     183             : }
     184             : 
     185          24 : void AdjacencyMatrixBase::setupThirdAtomBlock( const std::vector<AtomNumber>& tc, std::vector<AtomNumber>& t ) {
     186          24 :   threeblocks.resize( tc.size() );
     187          24 :   unsigned base=t.size();
     188       53168 :   for(unsigned i=0; i<tc.size(); ++i) {
     189       53144 :     log.printf("%d ", tc[i].serial());
     190       53144 :     t.push_back(tc[i]);
     191       53144 :     threeblocks[i]=base+i;
     192             :   }
     193          24 :   log.printf("\n");
     194          24 : }
     195             : 
     196         326 : void AdjacencyMatrixBase::setLinkCellCutoff( const bool& symmetric, const double& lcut, double tcut ) {
     197         326 :   if( read_one_group && symmetric ) {
     198         184 :     getPntrToComponent(0)->setSymmetric( true );
     199             :   }
     200         326 :   if( nl_cut>0 && lcut>nl_cut ) {
     201           0 :     error("D_MAX for switching functions should be shorter than neighbor list cutoff");
     202             :   }
     203             : 
     204         326 :   if( tcut<0 ) {
     205         318 :     tcut=lcut;
     206             :   }
     207         326 :   if( nl_cut>0 ) {
     208           0 :     linkcells.setCutoff( nl_cut );
     209             :   } else {
     210         326 :     linkcells.setCutoff( lcut );
     211             :   }
     212         326 :   if( linkcells.getCutoff()<std::numeric_limits<double>::max() ) {
     213         187 :     log.printf("  set link cell cutoff to %f \n", linkcells.getCutoff() );
     214             :   }
     215         326 :   threecells.setCutoff( tcut );
     216         326 : }
     217             : 
     218       11700 : void AdjacencyMatrixBase::prepare() {
     219       11700 :   ActionWithVector::prepare();
     220       11700 :   neighbour_list_updated=false;
     221       11700 : }
     222             : 
     223       11750 : void AdjacencyMatrixBase::updateNeighbourList() {
     224       11750 :   neighbour_list_updated=true;
     225             :   // Build link cells here so that this is done in stream if it needed in stream
     226       11750 :   if( getStep()%nl_stride==0 ) {
     227             :     // Build the link cells
     228       11750 :     std::vector<Vector> ltmp_pos( ablocks.size() );
     229     1194628 :     for(unsigned i=0; i<ablocks.size(); ++i) {
     230     1182878 :       ltmp_pos[i]=ActionAtomistic::getPosition( ablocks[i] );
     231             :     }
     232       11750 :     linkcells.buildCellLists( ltmp_pos, ablocks, getPbc() );
     233             :     // 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.
     234             :     // In this way we ensure that the neighbour list doesn't get too big.  Also this number should always be large enough
     235       11750 :     natoms_per_list = 27*linkcells.getMaxInCell();
     236       11750 :     nlist.resize( getConstPntrToComponent(0)->getShape()[0]*( 2 + natoms_per_list ) );
     237             :     // Set the number of neighbors to zero for all ranks
     238       11750 :     nlist.assign(nlist.size(),0);
     239             :     // Now get stuff to do parallel implementation
     240       11750 :     unsigned stride=comm.Get_size();
     241       11750 :     unsigned rank=comm.Get_rank();
     242       11750 :     if( runInSerial() ) {
     243             :       stride=1;
     244             :       rank=0;
     245             :     }
     246       11750 :     unsigned nt=OpenMP::getNumThreads();
     247       11750 :     if( nt*stride*10>getConstPntrToComponent(0)->getShape()[0] ) {
     248       11217 :       nt=getConstPntrToComponent(0)->getShape()[0]/stride/10;
     249             :     }
     250             :     if( nt==0 ) {
     251             :       nt=1;
     252             :     }
     253             :     // Create a vector from the input set of tasks
     254       11750 :     std::vector<unsigned> & pTaskList( getListOfActiveTasks(this) );
     255             : 
     256       11750 :     #pragma omp parallel num_threads(nt)
     257             :     {
     258             :       // Get the number of tasks we have to deal with
     259             :       unsigned ntasks=getConstPntrToComponent(0)->getShape()[0];
     260             :       if( nl_stride==1 ) {
     261             :         ntasks=pTaskList.size();
     262             :       }
     263             :       // Build a tempory nlist so we can do omp parallelism
     264             :       std::vector<unsigned> omp_nlist;
     265             :       if( nt>1 ) {
     266             :         omp_nlist.resize( nlist.size(), 0 );
     267             :       }
     268             :       // Now run over all atoms and construct the link cells
     269             :       std::vector<Vector> t_atoms( 1+ablocks.size() );
     270             :       std::vector<unsigned> indices( 1+ablocks.size() ), cells_required( linkcells.getNumberOfCells() );
     271             :       #pragma omp for nowait
     272             :       for(unsigned i=rank; i<ntasks; i+=stride) {
     273             :         // Retrieve cells required from link cells - for matrix blocks
     274             :         unsigned ncells_required=0;
     275             :         linkcells.addRequiredCells( linkcells.findMyCell( ActionAtomistic::getPosition(pTaskList[i]) ), ncells_required, cells_required );
     276             :         // Now get the indices of the atoms in the link cells positions
     277             :         unsigned natoms=1;
     278             :         indices[0]=pTaskList[i];
     279             :         linkcells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
     280             :         if( nl_stride==1 ) {
     281             :           if( nt>1 ) {
     282             :             omp_nlist[indices[0]]=0;
     283             :           } else {
     284             :             nlist[indices[0]] = 0;
     285             :           }
     286             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + indices[0]*(1+natoms_per_list);
     287             :           for(unsigned j=0; j<natoms; ++j) {
     288             :             if( nt>1 ) {
     289             :               omp_nlist[ lstart + omp_nlist[indices[0]] ] = indices[j];
     290             :               omp_nlist[indices[0]]++;
     291             :             } else {
     292             :               nlist[ lstart + nlist[indices[0]] ] = indices[j];
     293             :               nlist[indices[0]]++;
     294             :             }
     295             :           }
     296             :         } else {
     297             :           // Get the positions of all the atoms in the link cells relative to the central atom
     298             :           for(unsigned j=0; j<natoms; ++j) {
     299             :             t_atoms[j] = ActionAtomistic::getPosition(indices[j]) - ActionAtomistic::getPosition(indices[0]);
     300             :           }
     301             :           if( !nopbc ) {
     302             :             pbcApply( t_atoms, natoms );
     303             :           }
     304             :           // Now construct the neighbor list
     305             :           if( nt>1 ) {
     306             :             omp_nlist[indices[0]] = 0;
     307             :           } else {
     308             :             nlist[indices[0]] = 0;
     309             :           }
     310             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + indices[0]*(1+natoms_per_list);
     311             :           for(unsigned j=0; j<natoms; ++j) {
     312             :             double d2;
     313             :             if ( (d2=t_atoms[j][0]*t_atoms[j][0])<nl_cut2 &&
     314             :                  (d2+=t_atoms[j][1]*t_atoms[j][1])<nl_cut2 &&
     315             :                  (d2+=t_atoms[j][2]*t_atoms[j][2])<nl_cut2 ) {
     316             :               if( nt>1 ) {
     317             :                 omp_nlist[ lstart + omp_nlist[indices[0]] ] = indices[j];
     318             :                 omp_nlist[indices[0]]++;
     319             :               } else {
     320             :                 nlist[ lstart + nlist[indices[0]] ] = indices[j];
     321             :                 nlist[indices[0]]++;
     322             :               }
     323             :             }
     324             :           }
     325             : 
     326             :         }
     327             :       }
     328             :       // Gather OMP stuff
     329             :       #pragma omp critical
     330             :       if(nt>1) {
     331             :         for(unsigned i=0; i<ntasks; ++i) {
     332             :           nlist[pTaskList[i]]+=omp_nlist[pTaskList[i]];
     333             :         }
     334             :         for(unsigned i=0; i<ntasks; ++i) {
     335             :           unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + pTaskList[i]*(1+natoms_per_list);
     336             :           for(unsigned j=0; j<omp_nlist[pTaskList[i]]; ++j) {
     337             :             nlist[ lstart + j ] += omp_nlist[ lstart + j ];
     338             :           }
     339             :         }
     340             :       }
     341             :     }
     342             :     // MPI gather
     343       11750 :     if( !runInSerial() ) {
     344       11750 :       comm.Sum( nlist );
     345             :     }
     346             :   }
     347       11750 :   if( threeblocks.size()>0 ) {
     348          89 :     std::vector<Vector> ltmp_pos2( threeblocks.size() );
     349       66184 :     for(unsigned i=0; i<threeblocks.size(); ++i) {
     350       66095 :       ltmp_pos2[i]=ActionAtomistic::getPosition( threeblocks[i] );
     351             :     }
     352          89 :     threecells.buildCellLists( ltmp_pos2, threeblocks, getPbc() );
     353             :   }
     354             :   // And finally work out maximum number of columns to use
     355       11750 :   maxcol = nlist[0];
     356     1140150 :   for(unsigned i=1; i<getConstPntrToComponent(0)->getShape()[0]; ++i) {
     357     1128400 :     if( nlist[i]>maxcol ) {
     358        1136 :       maxcol = nlist[i];
     359             :     }
     360             :   }
     361       11750 : }
     362             : 
     363         216 : void AdjacencyMatrixBase::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
     364         216 :   if( action==this ) {
     365         162 :     return;
     366             :   }
     367             :   // Update the neighbour list
     368          54 :   updateNeighbourList();
     369             : 
     370          54 :   unsigned nactive = atasks.size();
     371          54 :   std::vector<unsigned> indlist( 1 + ablocks.size() + threeblocks.size() );
     372        3844 :   for(unsigned i=0; i<nactive; ++i) {
     373        3790 :     unsigned num = retrieveNeighbours( atasks[i], indlist );
     374      269574 :     for(unsigned j=0; j<num; ++j) {
     375             :       bool found=false;
     376    67587417 :       for(unsigned k=0; k<atasks.size(); ++k ) {
     377    67573541 :         if( indlist[j]==atasks[k] ) {
     378             :           found=true;
     379             :           break;
     380             :         }
     381             :       }
     382      265784 :       if( !found ) {
     383       13876 :         atasks.push_back( indlist[j] );
     384             :       }
     385             :     }
     386             :   }
     387             : }
     388             : 
     389      409666 : unsigned AdjacencyMatrixBase::retrieveNeighbours( const unsigned& current, std::vector<unsigned> & indices ) const {
     390      409666 :   unsigned natoms=nlist[current];
     391      409666 :   indices[0]=current;
     392      409666 :   unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + current*(1+natoms_per_list);
     393             :   plumed_dbg_assert( nlist[lstart]==current );
     394    21916923 :   for(unsigned i=1; i<nlist[current]; ++i) {
     395    21507257 :     indices[i] = nlist[ lstart + i ];
     396             :   }
     397      409666 :   return natoms;
     398             : }
     399             : 
     400      405876 : void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<unsigned> & indices, MultiValue& myvals ) const {
     401             :   // Now retrieve bookeeping arrays
     402      405876 :   if( indices.size()!=(1+ablocks.size()+threeblocks.size()) ) {
     403       66071 :     indices.resize( 1+ablocks.size()+threeblocks.size() );
     404             :   }
     405             : 
     406             :   // Now get the positions
     407      405876 :   unsigned natoms=retrieveNeighbours( current, indices );
     408             :   unsigned ntwo_atoms=natoms;
     409             :   myvals.setSplitIndex( ntwo_atoms );
     410             : 
     411             :   // Now retrieve everything for the third atoms
     412      405876 :   if( threeblocks.size()>0 ) {
     413         920 :     unsigned ncells_required=0;
     414         920 :     std::vector<unsigned> cells_required( threecells.getNumberOfCells() );
     415         920 :     threecells.addRequiredCells( threecells.findMyCell( ActionAtomistic::getPosition(current) ), ncells_required, cells_required );
     416         920 :     threecells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
     417             :   }
     418      405876 :   myvals.setNumberOfIndices( natoms );
     419             : 
     420             : // Apply periodic boundary conditions to atom positions
     421             :   std::vector<std::vector<Vector> > & t_atoms( myvals.getFirstAtomDerivativeVector() );
     422      405876 :   if( t_atoms.size()!=1 ) {
     423       23826 :     t_atoms.resize(1);
     424             :   }
     425      405876 :   if( t_atoms[0].size()<getNumberOfAtoms() ) {
     426       22798 :     t_atoms[0].resize( getNumberOfAtoms() );
     427             :   }
     428    26226247 :   for(unsigned i=0; i<natoms; ++i) {
     429    25820371 :     t_atoms[0][i] = ActionAtomistic::getPosition(indices[i]) - ActionAtomistic::getPosition(current);
     430             :   }
     431      405876 :   if( !nopbc ) {
     432      405876 :     pbcApply( t_atoms[0], natoms );
     433             :   }
     434             :   // And collect atom position data
     435             :   std::vector<Vector> & atoms( myvals.getAtomVector() );
     436      405876 :   if( atoms.size()<getNumberOfAtoms() ) {
     437       23830 :     atoms.resize( getNumberOfAtoms() );
     438             :   }
     439    26226247 :   for(unsigned i=0; i<natoms; ++i) {
     440    25820371 :     atoms[ indices[i] ] = t_atoms[0][i];
     441             :   }
     442      405876 : }
     443             : 
     444    21245263 : void AdjacencyMatrixBase::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     445    21245263 :   Vector zero;
     446    21245263 :   zero.zero();
     447             :   plumed_dbg_assert( index2<myvals.getAtomVector().size() );
     448    21245263 :   double weight = calculateWeight( zero, myvals.getAtomVector()[index2], myvals.getNumberOfIndices()-myvals.getSplitIndex(), myvals );
     449    21245263 :   unsigned w_ind = getConstPntrToComponent(0)->getPositionInStream();
     450    21245263 :   myvals.setValue( w_ind, weight );
     451    21245263 :   if( fabs(weight)<epsilon ) {
     452             :     myvals.setValue( w_ind, 0 );
     453     7114953 :     return;
     454             :   }
     455             : 
     456    14130310 :   if( !doNotCalculateDerivatives() ) {
     457             :     // Update dynamic list indices for central atom
     458     6759016 :     myvals.updateIndex( w_ind, 3*index1+0 );
     459     6759016 :     myvals.updateIndex( w_ind, 3*index1+1 );
     460     6759016 :     myvals.updateIndex( w_ind, 3*index1+2 );
     461             :     // Update dynamic list indices for atom forming this bond
     462     6759016 :     myvals.updateIndex( w_ind, 3*index2+0 );
     463     6759016 :     myvals.updateIndex( w_ind, 3*index2+1 );
     464     6759016 :     myvals.updateIndex( w_ind, 3*index2+2 );
     465             :     // Now look after all the atoms in the third block
     466             :     std::vector<unsigned> & indices( myvals.getIndices() );
     467     6872649 :     for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
     468      113633 :       myvals.updateIndex( w_ind, 3*indices[i]+0 );
     469      113633 :       myvals.updateIndex( w_ind, 3*indices[i]+1 );
     470      113633 :       myvals.updateIndex( w_ind, 3*indices[i]+2 );
     471             :     }
     472             :     // Update dynamic list indices for virial
     473     6759016 :     unsigned base = 3*getNumberOfAtoms();
     474    67590160 :     for(unsigned j=0; j<9; ++j) {
     475    60831144 :       myvals.updateIndex( w_ind, base+j );
     476             :     }
     477             :     // And the indices for the derivatives of the row of the matrix
     478     6759016 :     unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     479             :     std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     480     6759016 :     matrix_indices[nmat_ind+0]=3*index2+0;
     481     6759016 :     matrix_indices[nmat_ind+1]=3*index2+1;
     482     6759016 :     matrix_indices[nmat_ind+2]=3*index2+2;
     483     6759016 :     myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
     484             :   }
     485             : 
     486             :   // Calculate the components if we need them
     487    14130310 :   if( components ) {
     488     2415123 :     unsigned x_index = getConstPntrToComponent(1)->getPositionInStream();
     489     2415123 :     unsigned y_index = getConstPntrToComponent(2)->getPositionInStream();
     490     2415123 :     unsigned z_index = getConstPntrToComponent(3)->getPositionInStream();
     491     2415123 :     Vector atom = myvals.getAtomVector()[index2];
     492     2415123 :     myvals.setValue( x_index, atom[0] );
     493     2415123 :     myvals.setValue( y_index, atom[1] );
     494     2415123 :     myvals.setValue( z_index, atom[2] );
     495     2415123 :     if( !doNotCalculateDerivatives() ) {
     496     1625864 :       myvals.addDerivative( x_index, 3*index1+0, -1 );
     497     1625864 :       myvals.addDerivative( x_index, 3*index2+0, +1 );
     498     1625864 :       myvals.addDerivative( x_index, 3*index1+1, 0 );
     499     1625864 :       myvals.addDerivative( x_index, 3*index2+1, 0 );
     500     1625864 :       myvals.addDerivative( x_index, 3*index1+2, 0 );
     501     1625864 :       myvals.addDerivative( x_index, 3*index2+2, 0 );
     502     1625864 :       myvals.addDerivative( y_index, 3*index1+0, 0 );
     503     1625864 :       myvals.addDerivative( y_index, 3*index2+0, 0 );
     504     1625864 :       myvals.addDerivative( y_index, 3*index1+1, -1 );
     505     1625864 :       myvals.addDerivative( y_index, 3*index2+1, +1 );
     506     1625864 :       myvals.addDerivative( y_index, 3*index1+2, 0 );
     507     1625864 :       myvals.addDerivative( y_index, 3*index2+2, 0 );
     508     1625864 :       myvals.addDerivative( z_index, 3*index1+0, 0 );
     509     1625864 :       myvals.addDerivative( z_index, 3*index2+0, 0 );
     510     1625864 :       myvals.addDerivative( z_index, 3*index1+1, 0 );
     511     1625864 :       myvals.addDerivative( z_index, 3*index2+1, 0 );
     512     1625864 :       myvals.addDerivative( z_index, 3*index1+2, -1 );
     513     1625864 :       myvals.addDerivative( z_index, 3*index2+2, +1 );
     514     6503456 :       for(unsigned k=0; k<3; ++k) {
     515             :         // Update dynamic lists for central atom
     516     4877592 :         myvals.updateIndex( x_index, 3*index1+k );
     517     4877592 :         myvals.updateIndex( y_index, 3*index1+k );
     518     4877592 :         myvals.updateIndex( z_index, 3*index1+k );
     519             :         // Update dynamic lists for bonded atom
     520     4877592 :         myvals.updateIndex( x_index, 3*index2+k );
     521     4877592 :         myvals.updateIndex( y_index, 3*index2+k );
     522     4877592 :         myvals.updateIndex( z_index, 3*index2+k );
     523             :       }
     524             :       // Add derivatives of virial
     525     1625864 :       unsigned base = 3*getNumberOfAtoms();
     526             :       // Virial for x
     527     1625864 :       myvals.addDerivative( x_index, base+0, -atom[0] );
     528     1625864 :       myvals.addDerivative( x_index, base+3, -atom[1] );
     529     1625864 :       myvals.addDerivative( x_index, base+6, -atom[2] );
     530     1625864 :       myvals.addDerivative( x_index, base+1, 0 );
     531     1625864 :       myvals.addDerivative( x_index, base+4, 0 );
     532     1625864 :       myvals.addDerivative( x_index, base+7, 0 );
     533     1625864 :       myvals.addDerivative( x_index, base+2, 0 );
     534     1625864 :       myvals.addDerivative( x_index, base+5, 0 );
     535     1625864 :       myvals.addDerivative( x_index, base+8, 0 );
     536             :       // Virial for y
     537     1625864 :       myvals.addDerivative( y_index, base+0, 0 );
     538     1625864 :       myvals.addDerivative( y_index, base+3, 0 );
     539     1625864 :       myvals.addDerivative( y_index, base+6, 0 );
     540     1625864 :       myvals.addDerivative( y_index, base+1, -atom[0] );
     541     1625864 :       myvals.addDerivative( y_index, base+4, -atom[1] );
     542     1625864 :       myvals.addDerivative( y_index, base+7, -atom[2] );
     543     1625864 :       myvals.addDerivative( y_index, base+2, 0 );
     544     1625864 :       myvals.addDerivative( y_index, base+5, 0 );
     545     1625864 :       myvals.addDerivative( y_index, base+8, 0 );
     546             :       // Virial for z
     547     1625864 :       myvals.addDerivative( z_index, base+0, 0 );
     548     1625864 :       myvals.addDerivative( z_index, base+3, 0 );
     549     1625864 :       myvals.addDerivative( z_index, base+6, 0 );
     550     1625864 :       myvals.addDerivative( z_index, base+1, 0 );
     551     1625864 :       myvals.addDerivative( z_index, base+4, 0 );
     552     1625864 :       myvals.addDerivative( z_index, base+7, 0 );
     553     1625864 :       myvals.addDerivative( z_index, base+2, -atom[0] );
     554     1625864 :       myvals.addDerivative( z_index, base+5, -atom[1] );
     555     1625864 :       myvals.addDerivative( z_index, base+8, -atom[2] );
     556    16258640 :       for(unsigned k=0; k<9; ++k) {
     557    14632776 :         myvals.updateIndex( x_index, base+k );
     558    14632776 :         myvals.updateIndex( y_index, base+k );
     559    14632776 :         myvals.updateIndex( z_index, base+k );
     560             :       }
     561     6503456 :       for(unsigned k=1; k<4; ++k) {
     562     4877592 :         unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     563             :         std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     564     4877592 :         matrix_indices[nmat_ind+0]=3*index2+0;
     565     4877592 :         matrix_indices[nmat_ind+1]=3*index2+1;
     566     4877592 :         matrix_indices[nmat_ind+2]=3*index2+2;
     567     4877592 :         myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 );
     568             :       }
     569             :     }
     570             :   }
     571             : }
     572             : 
     573      405876 : void AdjacencyMatrixBase::runEndOfRowJobs( const unsigned& ind, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     574      405876 :   if( doNotCalculateDerivatives() ) {
     575             :     return;
     576             :   }
     577             : 
     578      737807 :   for(int k=0; k<getNumberOfComponents(); ++k) {
     579      411298 :     unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     580             :     std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     581      411298 :     plumed_assert( nmat_ind<matrix_indices.size() );
     582      411298 :     matrix_indices[nmat_ind+0]=3*ind+0;
     583      411298 :     matrix_indices[nmat_ind+1]=3*ind+1;
     584      411298 :     matrix_indices[nmat_ind+2]=3*ind+2;
     585      411298 :     nmat_ind+=3;
     586      423041 :     for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
     587       11743 :       matrix_indices[nmat_ind+0]=3*indices[i]+0;
     588       11743 :       matrix_indices[nmat_ind+1]=3*indices[i]+1;
     589       11743 :       matrix_indices[nmat_ind+2]=3*indices[i]+2;
     590       11743 :       nmat_ind+=3;
     591             :     }
     592      411298 :     unsigned virbase = 3*getNumberOfAtoms();
     593     4112980 :     for(unsigned i=0; i<9; ++i) {
     594     3701682 :       matrix_indices[nmat_ind+i]=virbase+i;
     595             :     }
     596      411298 :     nmat_ind+=9;
     597             :     plumed_dbg_massert( nmat_ind<=3*getNumberOfAtoms() + 9, "found too many derivatives in " + getLabel() );
     598             :     myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     599             :   }
     600             : }
     601             : 
     602             : }
     603             : }