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          20 :   multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
     296          60 :   keys.needsAction("CONTACT_MATRIX"); keys.needsAction("MATRIX_PRODUCT"); keys.needsAction("GROUP");
     297          60 :   keys.needsAction("ONES"); keys.needsAction("OUTER_PRODUCT"); keys.needsAction("VSTACK");
     298          60 :   keys.needsAction("CONCATENATE"); keys.needsAction("CUSTOM"); keys.needsAction("TRANSPOSE");
     299          20 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
     300          20 : }
     301             : 
     302          98 : std::string LocalSteinhardt::getSymbol( const int& m ) const {
     303          98 :   if( m<0 ) {
     304          40 :     std::string num; Tools::convert( -1*m, num );
     305          40 :     return "n" + num;
     306          58 :   } else if( m>0 ) {
     307          49 :     std::string num; Tools::convert( m, num );
     308          49 :     return "p" + num;
     309             :   }
     310           9 :   return "0";
     311             : }
     312             : 
     313           9 : std::string LocalSteinhardt::getArgsForStack( const int& l, const std::string& sp_lab ) const {
     314           9 :   std::string numstr; Tools::convert( l, numstr );
     315          18 :   std::string data_mat = " ARG=" + sp_lab + "_sp.rm-n" + numstr + "," + sp_lab + "_sp.im-n" + numstr;
     316         107 :   for(int i=-l+1; i<=l; ++i) {
     317          98 :     numstr = getSymbol( i );
     318         196 :     data_mat += "," + sp_lab + "_sp.rm-" + numstr + "," + sp_lab + "_sp.im-" + numstr;
     319             :   }
     320           9 :   return data_mat;
     321             : }
     322             : 
     323           7 : LocalSteinhardt::LocalSteinhardt(const ActionOptions& ao):
     324             :   Action(ao),
     325           7 :   ActionShortcut(ao)
     326             : {
     327           7 :   bool lowmem; parseFlag("LOWMEM",lowmem);
     328           7 :   if( lowmem ) warning("LOWMEM flag is deprecated and is no longer required for this action");
     329             :   // Get the Q value
     330          14 :   int l; Tools::convert( getName().substr(7), l);
     331             :   // Create a vector filled with ones
     332           7 :   std::string twolplusone; Tools::convert( 2*(2*l+1), twolplusone );
     333          14 :   readInputLine( getShortcutLabel() + "_uvec: ONES SIZE=" + twolplusone );
     334             :   // Read in species keyword
     335          14 :   std::string sp_str; parse("SPECIES",sp_str);
     336          14 :   std::string spa_str; parse("SPECIESA",spa_str);
     337           7 :   if( sp_str.length()>0 ) {
     338             :     // Create a group with these atoms
     339          12 :     readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
     340           6 :     std::vector<std::string> sp_lab = Tools::getWords(sp_str, "\t\n ,");
     341             :     // This creates the stash to hold all the vectors
     342           6 :     if( sp_lab.size()==1 ) {
     343             :       // The lengths of all the vectors in a vector
     344          12 :       readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + sp_lab[0] + "_norm," + getShortcutLabel() + "_uvec");
     345             :       // The unormalised vectors
     346          12 :       readInputLine( getShortcutLabel() + "_uvecs: VSTACK" + getArgsForStack( l, sp_lab[0] ) );
     347             :     } else {
     348           0 :       std::string len_vec = getShortcutLabel() + "_mags: CONCATENATE ARG=" + sp_lab[0] + "_norm";
     349           0 :       for(unsigned i=1; i<sp_lab.size(); ++i) len_vec += "," + sp_lab[i] + "_norm";
     350             :       // This is the vector that contains all the magnitudes
     351           0 :       readInputLine( len_vec );
     352           0 :       std::string concat_str = getShortcutLabel() + "_uvecs: CONCATENATE";
     353           0 :       for(unsigned i=0; i<sp_lab.size(); ++i) {
     354           0 :         std::string snum; Tools::convert( i+1, snum );
     355           0 :         concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecs" + snum;
     356           0 :         readInputLine( getShortcutLabel() + "_uvecs" + snum + ": VSTACK" + getArgsForStack( l, sp_lab[i] ) );
     357             :       }
     358             :       // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
     359           0 :       readInputLine( getShortcutLabel() + "_nmat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_mags," + getShortcutLabel() + "_uvec");
     360             :       // The unormalised vectors
     361           0 :       readInputLine( concat_str );
     362             :     }
     363             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     364          12 :     readInputLine( getShortcutLabel() + "_vecs: CUSTOM ARG=" + getShortcutLabel() + "_uvecs," + getShortcutLabel() + "_nmat FUNC=x/y PERIODIC=NO");
     365             :     // And transpose the matrix
     366          12 :     readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
     367          18 :     std::string sw_str; parse("SWITCH",sw_str); readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_str + " SWITCH={" + sw_str + "}");
     368             :     // And the matrix of dot products
     369          12 :     readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT" );
     370           7 :   } else if( spa_str.length()>0 ) {
     371             :     // Create a group with these atoms
     372           2 :     readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + spa_str );
     373           2 :     std::string spb_str; parse("SPECIESB",spb_str);
     374           1 :     if( spb_str.length()==0 ) plumed_merror("need both SPECIESA and SPECIESB in input");
     375           1 :     std::vector<std::string> sp_laba = Tools::getWords(spa_str, "\t\n ,");
     376           1 :     std::vector<std::string> sp_labb = Tools::getWords(spb_str, "\t\n ,");
     377           1 :     if( sp_laba.size()==1 ) {
     378             :       // The matrix that is used for normalising
     379           2 :       readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" +  sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
     380             :       // The unormalised vectors
     381           2 :       readInputLine( getShortcutLabel() + "_uvecsA: VSTACK" + getArgsForStack( l, sp_laba[0] ) );
     382             :     } else {
     383           0 :       std::string len_vec = getShortcutLabel() + "_magsA: CONCATENATE ARG=" + sp_laba[0] + "_norm";
     384           0 :       for(unsigned i=1; i<sp_laba.size(); ++i) len_vec += "," + sp_laba[i] + "_norm";
     385             :       //  This is the vector that contains all the magnitudes
     386           0 :       readInputLine( len_vec );
     387           0 :       std::string concat_str = getShortcutLabel() + "_uvecsA: CONCATENATE";
     388           0 :       for(unsigned i=0; i<sp_laba.size(); ++i) {
     389           0 :         std::string snum; Tools::convert( i+1, snum );
     390           0 :         concat_str += " MATRIX" + snum + "1=" + getShortcutLabel() + "_uvecsA" + snum;
     391           0 :         readInputLine( getShortcutLabel() + "_uvecsA" + snum + ": VSTACK" + getArgsForStack( l, sp_laba[i] ) );
     392             :       }
     393             :       // And the normalising matrix by taking the column vector of magnitudes and multiplying by the row vector of ones
     394           0 :       readInputLine( getShortcutLabel() + "_nmatA: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_magsA," + getShortcutLabel() + "_uvec");
     395             :       // The unormalised vector
     396           0 :       readInputLine( concat_str );
     397             :     }
     398             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     399           2 :     readInputLine( getShortcutLabel() + "_vecsA: CUSTOM ARG=" + getShortcutLabel() + "_uvecsA," + getShortcutLabel() + "_nmatA FUNC=x/y PERIODIC=NO");
     400             :     // Now do second matrix
     401           1 :     if( sp_labb.size()==1 ) {
     402           0 :       readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" +  sp_laba[0] + "_norm," + getShortcutLabel() + "_uvec");
     403           0 :       readInputLine( getShortcutLabel() + "_uvecsBT: VSTACK" + getArgsForStack( l, sp_labb[0] ) );
     404           0 :       readInputLine( getShortcutLabel() + "_uvecsB: TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT");
     405             :     } else {
     406           2 :       std::string len_vec = getShortcutLabel() + "_magsB: CONCATENATE ARG=" +  sp_labb[0] + "_norm";
     407           2 :       for(unsigned i=1; i<sp_labb.size(); ++i) len_vec += "," + sp_labb[i] + "_norm";
     408             :       //  This is the vector that contains all the magnitudes
     409           1 :       readInputLine( len_vec );
     410           1 :       std::string concat_str = getShortcutLabel() + "_uvecsB: CONCATENATE";
     411           3 :       for(unsigned i=0; i<sp_labb.size(); ++i) {
     412           2 :         std::string snum; Tools::convert( i+1, snum );
     413           4 :         concat_str += " MATRIX1" + snum + "=" + getShortcutLabel() + "_uvecsB" + snum;
     414           4 :         readInputLine( getShortcutLabel() + "_uvecsBT" + snum + ": VSTACK" + getArgsForStack( l, sp_labb[i] ) );
     415           4 :         readInputLine( getShortcutLabel() + "_uvecsB" + snum + ": TRANSPOSE ARG=" + getShortcutLabel() + "_uvecsBT" + snum );
     416             :       }
     417             :       // And the normalising matrix
     418           2 :       readInputLine( getShortcutLabel() + "_nmatB: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_uvec," + getShortcutLabel() + "_magsB");
     419             :       // The unormalised vectors
     420           1 :       readInputLine( concat_str );
     421             :     }
     422             :     // Now normalise all the vectors by doing Hadammard "product" with normalising matrix
     423           2 :     readInputLine( getShortcutLabel() + "_vecsB: CUSTOM ARG=" + getShortcutLabel() + "_uvecsB," + getShortcutLabel() + "_nmatB FUNC=x/y PERIODIC=NO");
     424           3 :     std::string sw_str; parse("SWITCH",sw_str); readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + spa_str + " GROUPB=" + spb_str + " SWITCH={" + sw_str + "}");
     425           2 :     readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_vecsA," + getShortcutLabel() + "_vecsB" );
     426           1 :   }
     427             : 
     428             :   // Now create the product matrix
     429          14 :   readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_dpmat FUNC=x*y PERIODIC=NO");
     430             :   // Now the sum of coordination numbers times the switching functions
     431           7 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
     432           7 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     433           7 :   std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     434          14 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     435          14 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() +"_prod," + getShortcutLabel() +"_ones");
     436             :   // And just the sum of the coordination numbers
     437          14 :   readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() +"_ones");
     438             :   // And matheval to get the final quantity
     439          14 :   readInputLine( getShortcutLabel() + "_av: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_denom FUNC=x/y PERIODIC=NO");
     440             :   // And this expands everything
     441          14 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_av", "", this );
     442           7 : }
     443             : 
     444             : }
     445             : }