Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2017 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 "core/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "multicolvar/MultiColvarShortcuts.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "core/ActionSet.h"
      27             : #include "core/ActionWithValue.h"
      28             : 
      29             : namespace PLMD {
      30             : namespace symfunc {
      31             : 
      32             : //+PLUMEDOC MCOLVARF LOCAL_Q1
      33             : /*
      34             : Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
      35             : 
      36             : \par Examples
      37             : 
      38             : */
      39             : //+ENDPLUMEDOC
      40             : 
      41             : //+PLUMEDOC MCOLVARF LOCAL_Q3
      42             : /*
      43             : Calculate the local degree of order around an atoms by taking the average dot product between the q_3 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
      44             : 
      45             : The \ref Q3 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
      46             : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
      47             : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
      48             : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
      49             : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
      50             : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
      51             : because the number of atoms is relatively small.
      52             : 
      53             : When the average \ref Q3 parameter is used to bias the dynamics a problems
      54             : can occur. These averaged coordinates cannot distinguish between the correct,
      55             : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
      56             : themselves into their solid-like configuration simultaneously. This second type
      57             : of pathway would be impossible in reality because there is a large entropic
      58             : barrier that prevents concerted processes like this from happening.  However,
      59             : in the finite sized systems that are commonly simulated this barrier is reduced
      60             : substantially. As a result in simulations where average Steinhardt parameters
      61             : are biased there are often quite dramatic system size effects
      62             : 
      63             : If one wants to simulate nucleation using some form on
      64             : biased dynamics what is really required is an order parameter that measures:
      65             : 
      66             : - Whether or not the coordination spheres around atoms are ordered
      67             : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
      68             : 
      69             : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
      70             : LOCAL_Q3 is another variable that can be used in these sorts of calculations. The LOCAL_Q3 parameter for a particular atom is a number that measures the extent to
      71             : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom.  It does this by calculating the following
      72             : quantity for each of the atoms in the system:
      73             : 
      74             : \f[
      75             :  s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-3}^3 q_{3m}^{*}(i)q_{3m}(j) }{ \sum_j \sigma( r_{ij} ) }
      76             : \f]
      77             : 
      78             : where \f$q_{3m}(i)\f$ and \f$q_{3m}(j)\f$ are the 3rd order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes complex
      79             : conjugation.  The function
      80             : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$.  The parameters of this function should be set
      81             : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.  The sum in the numerator
      82             : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
      83             : adjacent atoms is correlated.
      84             : 
      85             : \par Examples
      86             : 
      87             : The following command calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
      88             : quantity to a file called colvar.
      89             : 
      90             : \plumedfile
      91             : Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
      92             : LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq3
      93             : PRINT ARG=lq3.mean FILE=colvar
      94             : \endplumedfile
      95             : 
      96             : The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
      97             : 
      98             : \plumedfile
      99             : Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
     100             : LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq3
     101             : PRINT ARG=lq3.* FILE=colvar
     102             : \endplumedfile
     103             : 
     104             : The following calculates the LOCAL_Q3 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
     105             : are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
     106             : 
     107             : \plumedfile
     108             : Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3a
     109             : Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3b
     110             : 
     111             : LOCAL_Q3 SPECIES=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w3
     112             : PRINT ARG=w3.* FILE=colvar
     113             : \endplumedfile
     114             : 
     115             : */
     116             : //+ENDPLUMEDOC
     117             : 
     118             : //+PLUMEDOC MCOLVARF LOCAL_Q4
     119             : /*
     120             : Calculate the local degree of order around an atoms by taking the average dot product between the q_4 vector on the central atom and the q_4 vector on the atoms in the first coordination sphere.
     121             : 
     122             : The \ref Q4 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
     123             : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
     124             : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
     125             : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
     126             : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
     127             : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
     128             : because the number of atoms is relatively small.
     129             : 
     130             : When the average \ref Q4 parameter is used to bias the dynamics a problems
     131             : can occur. These averaged coordinates cannot distinguish between the correct,
     132             : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
     133             : themselves into their solid-like configuration simultaneously. This second type
     134             : of pathway would be impossible in reality because there is a large entropic
     135             : barrier that prevents concerted processes like this from happening.  However,
     136             : in the finite sized systems that are commonly simulated this barrier is reduced
     137             : substantially. As a result in simulations where average Steinhardt parameters
     138             : are biased there are often quite dramatic system size effects
     139             : 
     140             : If one wants to simulate nucleation using some form on
     141             : biased dynamics what is really required is an order parameter that measures:
     142             : 
     143             : - Whether or not the coordination spheres around atoms are ordered
     144             : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
     145             : 
     146             : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
     147             : LOCAL_Q4 is another variable that can be used in these sorts of calculations. The LOCAL_Q4 parameter for a particular atom is a number that measures the extent to
     148             : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom.  It does this by calculating the following
     149             : quantity for each of the atoms in the system:
     150             : 
     151             : \f[
     152             :  s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-4}^4 q_{4m}^{*}(i)q_{4m}(j) }{ \sum_j \sigma( r_{ij} ) }
     153             : \f]
     154             : 
     155             : where \f$q_{4m}(i)\f$ and \f$q_{4m}(j)\f$ are the fourth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
     156             : complex conjugation.  The function
     157             : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$.  The parameters of this function should be set
     158             : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.  The sum in the numerator
     159             : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
     160             : adjacent atoms is correlated.
     161             : 
     162             : \par Examples
     163             : 
     164             : The following command calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
     165             : quantity to a file called colvar.
     166             : 
     167             : \plumedfile
     168             : Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
     169             : LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq4
     170             : PRINT ARG=lq4.mean FILE=colvar
     171             : \endplumedfile
     172             : 
     173             : The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
     174             : 
     175             : \plumedfile
     176             : Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
     177             : LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq4
     178             : PRINT ARG=lq4.* FILE=colvar
     179             : \endplumedfile
     180             : 
     181             : The following calculates the LOCAL_Q4 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
     182             : are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
     183             : 
     184             : \plumedfile
     185             : Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4a
     186             : Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4b
     187             : 
     188             : LOCAL_Q4 SPECIES=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w4
     189             : PRINT ARG=w4.* FILE=colvar
     190             : \endplumedfile
     191             : 
     192             : */
     193             : //+ENDPLUMEDOC
     194             : 
     195             : //+PLUMEDOC MCOLVARF LOCAL_Q6
     196             : /*
     197             : Calculate the local degree of order around an atoms by taking the average dot product between the q_6 vector on the central atom and the q_6 vector on the atoms in the first coordination sphere.
     198             : 
     199             : The \ref Q6 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
     200             : around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
     201             : measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
     202             : examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
     203             : atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
     204             : change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
     205             : because the number of atoms is relatively small.
     206             : 
     207             : When the average \ref Q6 parameter is used to bias the dynamics a problems
     208             : can occur. These averaged coordinates cannot distinguish between the correct,
     209             : single-nucleus pathway and a concerted pathway in which all the atoms rearrange
     210             : themselves into their solid-like configuration simultaneously. This second type
     211             : of pathway would be impossible in reality because there is a large entropic
     212             : barrier that prevents concerted processes like this from happening.  However,
     213             : in the finite sized systems that are commonly simulated this barrier is reduced
     214             : substantially. As a result in simulations where average Steinhardt parameters
     215             : are biased there are often quite dramatic system size effects
     216             : 
     217             : If one wants to simulate nucleation using some form on
     218             : biased dynamics what is really required is an order parameter that measures:
     219             : 
     220             : - Whether or not the coordination spheres around atoms are ordered
     221             : - Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
     222             : 
     223             : \ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
     224             : LOCAL_Q6 is another variable that can be used in these sorts of calculations. The LOCAL_Q6 parameter for a particular atom is a number that measures the extent to
     225             : which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom.  It does this by calculating the following
     226             : quantity for each of the atoms in the system:
     227             : 
     228             : \f[
     229             :  s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
     230             : \f]
     231             : 
     232             : where \f$q_{6m}(i)\f$ and \f$q_{6m}(j)\f$ are the sixth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
     233             : complex conjugation.  The function
     234             : \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$.  The parameters of this function should be set
     235             : so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.  The sum in the numerator
     236             : of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
     237             : adjacent atoms is correlated.
     238             : 
     239             : \par Examples
     240             : 
     241             : The following command calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
     242             : quantity to a file called colvar.
     243             : 
     244             : \plumedfile
     245             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
     246             : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq6
     247             : PRINT ARG=lq6.mean FILE=colvar
     248             : \endplumedfile
     249             : 
     250             : The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
     251             : 
     252             : \plumedfile
     253             : Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
     254             : LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq6
     255             : PRINT ARG=lq6.* FILE=colvar
     256             : \endplumedfile
     257             : 
     258             : The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
     259             : are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
     260             : 
     261             : \plumedfile
     262             : Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6a
     263             : Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6b
     264             : 
     265             : LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w6
     266             : PRINT ARG=w6.* FILE=colvar
     267             : \endplumedfile
     268             : 
     269             : */
     270             : //+ENDPLUMEDOC
     271             : 
     272             : class LocalSteinhardt : public ActionShortcut {
     273             : private:
     274             :   std::string getSymbol( const int& m ) const ;
     275             :   std::string getArgsForStack( const int& l, const std::string& lab ) const;
     276             : public:
     277             :   static void registerKeywords( Keywords& keys );
     278             :   explicit LocalSteinhardt(const ActionOptions&);
     279             : };
     280             : 
     281             : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q1")
     282             : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q3")
     283             : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q4")
     284             : PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q6")
     285             : 
     286          20 : void LocalSteinhardt::registerKeywords( Keywords& keys ) {
     287          20 :   ActionShortcut::registerKeywords( keys );
     288          40 :   keys.add("optional","SPECIES","");
     289          40 :   keys.add("optional","SPECIESA","");
     290          40 :   keys.add("optional","SPECIESB","");
     291          40 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
     292             :            "The following provides information on the \\ref switchingfunction that are available. "
     293             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
     294          40 :   keys.addFlag("LOWMEM",false,"this flag does nothing and is present only to ensure back-compatibility");
     295          40 :   keys.setValueDescription("vector","the values of the local steinhardt parameters for the input atoms");
     296          20 :   multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
     297          60 :   keys.needsAction("CONTACT_MATRIX"); keys.needsAction("MATRIX_PRODUCT"); keys.needsAction("GROUP");
     298          60 :   keys.needsAction("ONES"); keys.needsAction("OUTER_PRODUCT"); keys.needsAction("VSTACK");
     299          60 :   keys.needsAction("CONCATENATE"); keys.needsAction("CUSTOM"); keys.needsAction("TRANSPOSE");
     300          20 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
     301          20 : }
     302             : 
     303          98 : std::string LocalSteinhardt::getSymbol( const int& m ) const {
     304          98 :   if( m<0 ) {
     305          40 :     std::string num; Tools::convert( -1*m, num );
     306          40 :     return "n" + num;
     307          58 :   } else if( m>0 ) {
     308          49 :     std::string num; Tools::convert( m, num );
     309          49 :     return "p" + num;
     310             :   }
     311           9 :   return "0";
     312             : }
     313             : 
     314           9 : std::string LocalSteinhardt::getArgsForStack( const int& l, const std::string& sp_lab ) const {
     315           9 :   std::string numstr; Tools::convert( l, numstr );
     316          18 :   std::string data_mat = " ARG=" + sp_lab + "_sp.rm-n" + numstr + "," + sp_lab + "_sp.im-n" + numstr;
     317         107 :   for(int i=-l+1; i<=l; ++i) {
     318          98 :     numstr = getSymbol( i );
     319         196 :     data_mat += "," + sp_lab + "_sp.rm-" + numstr + "," + sp_lab + "_sp.im-" + numstr;
     320             :   }
     321           9 :   return data_mat;
     322             : }
     323             : 
     324           7 : LocalSteinhardt::LocalSteinhardt(const ActionOptions& ao):
     325             :   Action(ao),
     326           7 :   ActionShortcut(ao)
     327             : {
     328           7 :   bool lowmem; parseFlag("LOWMEM",lowmem);
     329           7 :   if( lowmem ) warning("LOWMEM flag is deprecated and is no longer required for this action");
     330             :   // Get the Q value
     331          14 :   int l; Tools::convert( getName().substr(7), l);
     332             :   // Create a vector filled with ones
     333           7 :   std::string twolplusone; Tools::convert( 2*(2*l+1), twolplusone );
     334          14 :   readInputLine( getShortcutLabel() + "_uvec: ONES SIZE=" + twolplusone );
     335             :   // Read in species keyword
     336          14 :   std::string sp_str; parse("SPECIES",sp_str);
     337          14 :   std::string spa_str; parse("SPECIESA",spa_str);
     338           7 :   if( sp_str.length()>0 ) {
     339             :     // Create a group with these atoms
     340          12 :     readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
     341           6 :     std::vector<std::string> sp_lab = Tools::getWords(sp_str, "\t\n ,");
     342             :     // This creates the stash to hold all the vectors
     343           6 :     if( sp_lab.size()==1 ) {
     344             :       // The lengths of all the vectors in a vector
     345          12 :       readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + sp_lab[0] + "_norm," + getShortcutLabel() + "_uvec");
     346             :       // The unormalised vectors
     347          12 :       readInputLine( getShortcutLabel() + "_uvecs: VSTACK" + getArgsForStack( l, sp_lab[0] ) );
     348             :     } else {
     349           0 :       std::string len_vec = getShortcutLabel() + "_mags: CONCATENATE ARG=" + sp_lab[0] + "_norm";
     350           0 :       for(unsigned i=1; i<sp_lab.size(); ++i) len_vec += "," + sp_lab[i] + "_norm";
     351             :       // This is the vector that contains all the magnitudes
     352           0 :       readInputLine( len_vec );
     353           0 :       std::string concat_str = getShortcutLabel() + "_uvecs: CONCATENATE";
     354           0 :       for(unsigned i=0; i<sp_lab.size(); ++i) {
     355           0 :         std::string snum; Tools::convert( i+1, snum );
     356           0 :         concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecs" + snum;
     357           0 :         readInputLine( getShortcutLabel() + "_uvecs" + snum + ": VSTACK" + getArgsForStack( l, sp_lab[i] ) );
     358             :       }
     359             :       // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
     360           0 :       readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_mags," + getShortcutLabel() + "_uvec");
     361             :       // The unormalised vectors
     362           0 :       readInputLine( concat_str );
     363             :     }
     364             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     365          12 :     readInputLine( getShortcutLabel() + "_vecs: CUSTOM ARG=" + getShortcutLabel() + "_uvecs," + getShortcutLabel() + "_nmat FUNC=x/y PERIODIC=NO");
     366             :     // And transpose the matrix
     367          12 :     readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
     368          18 :     std::string sw_str; parse("SWITCH",sw_str); readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_str + " SWITCH={" + sw_str + "}");
     369             :     // And the matrix of dot products
     370          12 :     readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT" );
     371           7 :   } else if( spa_str.length()>0 ) {
     372             :     // Create a group with these atoms
     373           2 :     readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + spa_str );
     374           2 :     std::string spb_str; parse("SPECIESB",spb_str);
     375           1 :     if( spb_str.length()==0 ) plumed_merror("need both SPECIESA and SPECIESB in input");
     376           1 :     std::vector<std::string> sp_laba = Tools::getWords(spa_str, "\t\n ,");
     377           1 :     std::vector<std::string> sp_labb = Tools::getWords(spb_str, "\t\n ,");
     378           1 :     if( sp_laba.size()==1 ) {
     379             :       // The matrix that is used for normalising
     380           2 :       readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" +  sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
     381             :       // The unormalised vectors
     382           2 :       readInputLine( getShortcutLabel() + "_uvecsA: VSTACK" + getArgsForStack( l, sp_laba[0] ) );
     383             :     } else {
     384           0 :       std::string len_vec = getShortcutLabel() + "_magsA: CONCATENATE ARG=" + sp_laba[0] + "_norm";
     385           0 :       for(unsigned i=1; i<sp_laba.size(); ++i) len_vec += "," + sp_laba[i] + "_norm";
     386             :       //  This is the vector that contains all the magnitudes
     387           0 :       readInputLine( len_vec );
     388           0 :       std::string concat_str = getShortcutLabel() + "_uvecsA: CONCATENATE";
     389           0 :       for(unsigned i=0; i<sp_laba.size(); ++i) {
     390           0 :         std::string snum; Tools::convert( i+1, snum );
     391           0 :         concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecsA" + snum;
     392           0 :         readInputLine( getShortcutLabel() + "_uvecsA" + snum + ": VSTACK" + getArgsForStack( l, sp_laba[i] ) );
     393             :       }
     394             :       // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
     395           0 :       readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_magsA," + getShortcutLabel() + "_uvec");
     396             :       // The unormalised vector
     397           0 :       readInputLine( concat_str );
     398             :     }
     399             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     400           2 :     readInputLine( getShortcutLabel() + "_vecsA: CUSTOM ARG=" + getShortcutLabel() + "_uvecsA," + getShortcutLabel() + "_nmatA FUNC=x/y PERIODIC=NO");
     401             :     // Now do second matrix
     402           1 :     if( sp_labb.size()==1 ) {
     403           0 :       readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" +  sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
     404           0 :       readInputLine( getShortcutLabel() + "_uvecsBT: VSTACK" + getArgsForStack( l, sp_labb[0] ) );
     405           0 :       readInputLine( getShortcutLabel() + "_uvecsB: TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT");
     406             :     } else {
     407           2 :       std::string len_vec = getShortcutLabel() + "_magsB: CONCATENATE ARG=" +  sp_labb[0] + "_norm";
     408           2 :       for(unsigned i=1; i<sp_labb.size(); ++i) len_vec += "," + sp_labb[i] + "_norm";
     409             :       //  This is the vector that contains all the magnitudes
     410           1 :       readInputLine( len_vec );
     411           1 :       std::string concat_str = getShortcutLabel() + "_uvecsB: CONCATENATE";
     412           3 :       for(unsigned i=0; i<sp_labb.size(); ++i) {
     413           2 :         std::string snum; Tools::convert( i+1, snum );
     414           4 :         concat_str += " MATRIX1" + snum + "=" + getShortcutLabel() + "_uvecsB" + snum;
     415           4 :         readInputLine( getShortcutLabel() + "_uvecsBT" + snum + ": VSTACK" + getArgsForStack( l, sp_labb[i] ) );
     416           4 :         readInputLine( getShortcutLabel() + "_uvecsB" + snum + ": TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT" + snum );
     417             :       }
     418             :       // And the normalising matrix
     419           2 :       readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_uvec," + getShortcutLabel() + "_magsB");
     420             :       // The unormalised vectors
     421           1 :       readInputLine( concat_str );
     422             :     }
     423             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     424           2 :     readInputLine( getShortcutLabel() + "_vecsB: CUSTOM ARG=" + getShortcutLabel() + "_uvecsB," + getShortcutLabel() + "_nmatB FUNC=x/y PERIODIC=NO");
     425           3 :     std::string sw_str; parse("SWITCH",sw_str); readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + spa_str + " GROUPB=" + spb_str + " SWITCH={" + sw_str + "}");
     426           2 :     readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecsA," + getShortcutLabel() + "_vecsB" );
     427           1 :   }
     428             : 
     429             :   // Now create the product matrix
     430          14 :   readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_dpmat FUNC=x*y PERIODIC=NO");
     431             :   // Now the sum of coordination numbers times the switching functions
     432           7 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
     433           7 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     434           7 :   std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     435          14 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     436          14 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() +"_prod," + getShortcutLabel() +"_ones");
     437             :   // And just the sum of the coordination numbers
     438          14 :   readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() +"_ones");
     439             :   // And matheval to get the final quantity
     440          14 :   readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
     441             :   // And this expands everything
     442          14 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_av", "", this );
     443           7 : }
     444             : 
     445             : }
     446             : }