Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/ActionWithValue.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "core/ActionSet.h"
      27             : #include "multicolvar/MultiColvarShortcuts.h"
      28             : 
      29             : //+PLUMEDOC MCOLVAR SMAC
      30             : /*
      31             : Calculate the SMAC order parameter for a set of molecules
      32             : 
      33             : This shortcut action provides a quantity that can be thought of as a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) for molecules.
      34             : You can thus use it to detect whether the molecules in a sample are ordered in the way that they would expect to be ordered in the crystal. If you are using this CV you normally use
      35             : the [COM](COM.md) command to define centers of mass for your molecules and the [DISTANCE](DISTANCE.md) or [PLANE](PLANE.md) command to define each molecules orientation.  You can thus calculate
      36             : the relative orientation molecules $i$ and $j$ by calculating the torsional angle, $\theta_{ij}$, between the two molecular orientations around the vector $r_{ij}$ connecting the centers of
      37             : mass of the two molecules. The value of the molecular order parameter for molecule $i$ is thus calculated as:
      38             : 
      39             : $$
      40             : s_i = frac{1-\gamma(c_i)}{c_i} \sum_j \sigma(|r_{ij}|) \sum_k K\left( \frac{\theta_{ij} - \phi_k}{b} \right) \qquad \textrm{where} \qquad c_i = \sum_j \sigma(|r_{ij}|)
      41             : $$
      42             : 
      43             : In this expression $\sigma$ is a switching function that acts on the distance between the centers of mass of molecules $i$ and $j$. $c_i$ is thus the number of molecules that are within
      44             : a certain cutoff of molecule $i$ and $\gamma$ is another switching function that acts upon this quantity. This switching function ensures that the symmetry function is zero for atoms that are
      45             : regions where the density of molecules are low.  $K$ is then a kernel function with bandwidth $b$ that is centered at $\phi_k$.  The function above is thus only large if molecule $i$ is surrounded
      46             : by molecules whose relative orientations are as the user has requested by specifying $\phi_k$ parameters.
      47             : 
      48             : The following example illustrates how the SMAC order parameter in PLUMED is used:
      49             : 
      50             : ```plumed
      51             : m3: DISTANCES ...
      52             :    ATOMS1=9,10 LOCATION1=9
      53             :    ATOMS2=89,90 LOCATION2=89
      54             :    ATOMS3=473,474 LOCATION3=473
      55             :    ATOMS4=1161,1162 LOCATION4=1161
      56             :    ATOMS5=1521,1522 LOCATION5=1521
      57             :    ATOMS6=1593,1594 LOCATION6=1593
      58             :    ATOMS7=1601,1602 LOCATION7=1601
      59             :    ATOMS8=2201,2202 LOCATION8=2201
      60             :    COMPONENTS
      61             : ...
      62             : 
      63             : s2: SMAC SPECIES=m3 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
      64             : 
      65             : PRINT ARG=s2_morethan FILE=colvar
      66             : ```
      67             : 
      68             : Here the orientations of the molecules are specified by calculating the vectors connecting pairs of atoms in the molecules.  The LOCATION keywords in the distance command are used to specify the
      69             : positions of the molecules from which the $r_{ij}$ vectors in the above expression are calculated.  The sum over $k$ in the above expression has two terms corresponding to the two Gaussian kernels
      70             : that have been specified using KERNEL keywords.  The SWITCH keyword has been used to specify the parameters of the switching function $\sigma$, while the SWITCH_COORD keyword has been used to specify
      71             : the parameters of the switching function $\gamma$.
      72             : 
      73             : A vector of 8 values for $s_i$ is calculated.  The elements of this vector are transformed by a switching function that is one if the $s_i$ value is larger than 0.7.  The eight elements of the resulting vector
      74             : of transformed $s_i$ are then added together to determine the final quantity that is output to the colvar file.
      75             : 
      76             : Incidentally, the authors who designed the SMAC symmetry function have forgotten what the letters in this acronym stand for.
      77             : 
      78             : */
      79             : //+ENDPLUMEDOC
      80             : 
      81             : namespace PLMD {
      82             : namespace symfunc {
      83             : 
      84             : class SMAC : public ActionShortcut {
      85             : public:
      86             :   static void registerKeywords(Keywords& keys);
      87             :   explicit SMAC(const ActionOptions&);
      88             : };
      89             : 
      90             : PLUMED_REGISTER_ACTION(SMAC,"SMAC")
      91             : 
      92          13 : void SMAC::registerKeywords(Keywords& keys) {
      93          13 :   ActionShortcut::registerKeywords( keys );
      94          13 :   keys.add("optional","SPECIES","the label of the DISTANCES or PLANES action that computes the orientations of the molecules for which you would like to compute SMAC.  The coordination sphere contains the same list of molecules if you use this keyword.");
      95          13 :   keys.add("optional","SPECIESA","the label of the DISTANCES or PLANES action that computes the orientations of the molecules for which you would like to compute SMAC.  This keyword must be used with SPECIESB.");
      96          13 :   keys.add("optional","SPECIESB","the label of the DISTANCES or PLANES action that computes the orientations of the molecules that should be considered as part of the coordination sphere.  This keyword must be used with SPECIESA.");
      97          13 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
      98             :            "The following provides information on the \\ref switchingfunction that are available. "
      99             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
     100          26 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
     101          13 :   keys.add("numbered","KERNEL","The kernels used in the function of the angle");
     102          13 :   keys.add("optional","SWITCH_COORD","This keyword is used to define the coordination switching function.");
     103          26 :   keys.linkActionInDocs("SWITCH_COORD","LESS_THAN");
     104          26 :   keys.reset_style("KERNEL","optional");
     105          26 :   keys.setValueDescription("vector","the value of the smac parameter for each of the input molecules");
     106          13 :   multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
     107          13 :   keys.needsAction("VSTACK");
     108          13 :   keys.needsAction("TRANSPOSE");
     109          13 :   keys.needsAction("CONTACT_MATRIX");
     110          13 :   keys.needsAction("TORSIONS_MATRIX");
     111          13 :   keys.needsAction("COMBINE");
     112          13 :   keys.needsAction("CUSTOM");
     113          13 :   keys.needsAction("ONES");
     114          13 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
     115          13 :   keys.needsAction("MORE_THAN");
     116          13 :   keys.addDOI("10.1016/j.ces.2014.08.032");
     117          13 :   keys.addDOI("10.1021/acs.jctc.6b01073");
     118          13 : }
     119             : 
     120           5 : SMAC::SMAC(const ActionOptions& ao):
     121             :   Action(ao),
     122           5 :   ActionShortcut(ao) {
     123             :   // Create the matrices
     124             :   std::string sw_input;
     125          10 :   parse("SWITCH",sw_input);
     126             :   std::string sp_lab, sp_laba;
     127           5 :   parse("SPECIES",sp_lab);
     128          10 :   parse("SPECIESA",sp_laba);
     129           5 :   if( sp_lab.length()>0 ) {
     130           8 :     readInputLine( getShortcutLabel() + "_vecs: VSTACK ARG=" + sp_lab + ".x," + sp_lab + ".y," + sp_lab + ".z" );
     131           8 :     readInputLine( getShortcutLabel() + "_vecsT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecs" );
     132           8 :     readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_lab + " SWITCH={" + sw_input + "}");
     133           8 :     readInputLine( getShortcutLabel() + "_tpmat: TORSIONS_MATRIX ARG=" + getShortcutLabel() + "_vecs," + getShortcutLabel() + "_vecsT POSITIONS1=" + sp_lab + " POSITIONS2=" + sp_lab );
     134           1 :   } else if( sp_laba.length()>0 ) {
     135             :     std::string sp_labb;
     136           1 :     parse("SPECIESB",sp_labb);
     137           2 :     readInputLine( getShortcutLabel() + "_vecsa: VSTACK ARG=" + sp_laba + ".x," + sp_laba + ".y," + sp_laba + ".z" );
     138           2 :     readInputLine( getShortcutLabel() + "_vecsb: VSTACK ARG=" + sp_labb + ".x," + sp_labb + ".y," + sp_labb + ".z" );
     139           2 :     readInputLine( getShortcutLabel() + "_vecsbT: TRANSPOSE ARG=" + getShortcutLabel() + "_vecsb" );
     140           2 :     readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + sp_laba + " GROUPB=" + sp_labb + " SWITCH={" + sw_input + "}");
     141           2 :     readInputLine( getShortcutLabel() + "_tpmat: TORSIONS_MATRIX ARG=" + getShortcutLabel() + "_vecsa," + getShortcutLabel() + "_vecsbT POSITIONS1=" + sp_laba + " POSITIONS2=" + sp_labb );
     142             :   }
     143             :   // Now need the Gaussians
     144           5 :   std::string kmap_input= getShortcutLabel() + "_ksum: COMBINE PERIODIC=NO";
     145           5 :   for(unsigned i=1;; ++i) {
     146             :     std::string kstr_inpt, istr;
     147          15 :     Tools::convert( i, istr );
     148          30 :     if( !parseNumbered("KERNEL",i,kstr_inpt ) ) {
     149             :       break;
     150             :     }
     151          10 :     std::vector<std::string> words = Tools::getWords(kstr_inpt);
     152             :     std::string center, var;
     153          10 :     Tools::parse(words,"CENTER",center);
     154          10 :     Tools::parse(words,"SIGMA",var);
     155             :     double nsig;
     156          10 :     Tools::convert( var, nsig );
     157             :     std::string coeff;
     158          10 :     Tools::convert( 1/(nsig*nsig), coeff );
     159          20 :     readInputLine( getShortcutLabel() + "_kf" + istr + "_r2: COMBINE PERIODIC=NO ARG=" + getShortcutLabel() + "_tpmat COEFFICIENTS=" + coeff + " PARAMETERS=" + center + " POWERS=2");
     160          10 :     if( words[0]=="GAUSSIAN" ) {
     161          16 :       readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM PERIODIC=NO FUNC=exp(-x/2) ARG=" +  getShortcutLabel() + "_kf" + istr + "_r2" );
     162           2 :     } else if( words[0]=="TRIANGULAR" ) {
     163           4 :       readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM PERIODIC=NO FUNC=step(1-sqrt(x))*(1-sqrt(x)) ARG=" + getShortcutLabel() + "_kf" + istr + "_r2" );
     164             :     } else {
     165           0 :       readInputLine( getShortcutLabel() + "_kf" + istr + ": CUSTOM PERIODIC=NO FUNC=" + words[0] + " ARG=" + getShortcutLabel() + "_kf" + istr + "_r2" );
     166             :     }
     167          10 :     if( i==1 ) {
     168          10 :       kmap_input += " ARG=" + getShortcutLabel() + "_kf" + istr;
     169             :     } else {
     170          10 :       kmap_input += "," + getShortcutLabel() + "_kf" + istr;
     171             :     }
     172          20 :   }
     173           5 :   readInputLine( kmap_input );
     174             :   // Now create the product matrix
     175          10 :   readInputLine( getShortcutLabel() + "_prod: CUSTOM ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_ksum FUNC=x*y PERIODIC=NO");
     176             :   // Now the sum of coordination numbers times the switching functions
     177           5 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
     178           5 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     179             :   std::string size;
     180           5 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     181          10 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     182          10 :   readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_prod," + getShortcutLabel() + "_ones");
     183             :   // And just the sum of the coordination numbers
     184          10 :   readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap," + getShortcutLabel() + "_ones");
     185             :   // And the transformed switching functions
     186             :   std::string swcoord_str;
     187           5 :   parse("SWITCH_COORD",swcoord_str);
     188          10 :   readInputLine( getShortcutLabel() + "_mtdenom: MORE_THAN ARG=" + getShortcutLabel() + "_denom SWITCH={" + swcoord_str +"}");
     189             : // And matheval to get the final quantity
     190          10 :   readInputLine( getShortcutLabel() + "_smac: CUSTOM ARG=" + getShortcutLabel() + "," + getShortcutLabel() + "_mtdenom," + getShortcutLabel() + "_denom FUNC=(x*y)/z PERIODIC=NO");
     191             :   // And this expands everything
     192          10 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_smac", "", this );
     193           5 : }
     194             : 
     195             : }
     196             : }