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/PlumedMain.h"
      25             : #include "core/ActionSet.h"
      26             : #include "core/ActionWithValue.h"
      27             : #include "multicolvar/MultiColvarShortcuts.h"
      28             : 
      29             : //+PLUMEDOC MCOLVAR ATOMIC_SMAC
      30             : /*
      31             : Calculate the atomic smac CV
      32             : 
      33             : This shortcut offers another example of a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html).
      34             : This action was inspired by [SMAC](SMAC.md) and the distribution of angles in the first coordination sphere around an atom to determine if the
      35             : environment around the atom is ordered. For atom $i$ the symmetry function is calculated using:
      36             : 
      37             : $$
      38             : s_i = [ 1 - \gamma(c_i) ] \frac{ \sum_j \sum{k \ne j} \sigma(r_{ij})\sigma(r_{ik}) \sum_n G\left( \frac{ \theta_{jik} - \phi_n}{b_n}\right) }{\sum{k \ne j} \sigma(r_{ij})\sigma(r_{ik})} \qquad \textrm{where} \qquad c_i = \sum_j \sigma(r_{ij})
      39             : $$
      40             : 
      41             : In this expression $r_{ij}$ is the distance between atom $i$ and atom $j$ and $\sigma$ is a switching function that acts upon this distance.   $c_i$ is thus the number of atoms that are within
      42             : a certain cutoff of atom $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
      43             : regions where the density is low.  $\theta_{jik}$ is the angle between the vector connecting atoms $i$ and $j$ and the vector connecting atoms $i$ and $k$.  This angle is the argument for the
      44             : set of Gaussian kernel functions, $G$, that are centered on $\phi_n$ and that have bandwidths of $b_n$.  The function above is thus determining if the angles between the bonds in the first coordination
      45             : sphere around atom $i$ are similar to the $\phi_n$ values that have been specified by the user or not.
      46             : 
      47             : The following example demonstrates how this symmetry function can be used in practise.
      48             : 
      49             : ```plumed
      50             : smac: ATOMIC_SMAC SPECIES=1-64 KERNEL1={GAUSSIAN CENTER=pi/2 SIGMA=1.0} KERNEL2={GAUSSIAN CENTER=pi/4 SIGMA=1.0} SWITCH_COORD={EXP R_0=4.0} SWITCH={RATIONAL R_0=2.0 D_0=2.0} SUM
      51             : PRINT ARG=smac.* FILE=colvar
      52             : ```
      53             : 
      54             : The input above would calculate 64 instances of $s_i$ using the formula above.  In each of these two Gaussian Kernels are used in the sum over $n$.  The parameters for the switching function
      55             : $\sigma$ are specified using the SWITCH keyword, while the parameters for $\gamma$ are specified using SWITCH_COORD.  As you can see if you expand the shortcut in the input above, the 64 values
      56             : for $s_i$ are stored in a vector.  All the elements of this vector are then added together to produce the single quantity that is output in the colvar file.
      57             : 
      58             : */
      59             : //+ENDPLUMEDOC
      60             : 
      61             : namespace PLMD {
      62             : namespace symfunc {
      63             : 
      64             : class AtomicSMAC : public ActionShortcut {
      65             : public:
      66             :   static void registerKeywords(Keywords& keys);
      67             :   explicit AtomicSMAC(const ActionOptions&);
      68             : };
      69             : 
      70             : PLUMED_REGISTER_ACTION(AtomicSMAC,"ATOMIC_SMAC")
      71             : 
      72           5 : void AtomicSMAC::registerKeywords(Keywords& keys) {
      73           5 :   ActionShortcut::registerKeywords( keys );
      74           5 :   keys.add("atoms-3","SPECIES","the list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments");
      75           5 :   keys.add("atoms-4","SPECIESA","the list of atoms for which the symmetry function is being calculated.  This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment.");
      76           5 :   keys.add("atoms-4","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated.  This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which the symmetry function is being calculated.");
      77           5 :   keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
      78           5 :   keys.add("numbered","KERNEL","The kernels used in the function of the angle");
      79           5 :   keys.add("optional","SWITCH_COORD","This keyword is used to define the coordination switching function.");
      80          10 :   keys.reset_style("KERNEL","optional");
      81           5 :   multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
      82           5 :   keys.needsAction("CONTACT_MATRIX");
      83           5 :   keys.needsAction("GSYMFUNC_THREEBODY");
      84           5 :   keys.needsAction("ONES");
      85           5 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
      86           5 : }
      87             : 
      88           1 : AtomicSMAC::AtomicSMAC(const ActionOptions& ao):
      89             :   Action(ao),
      90           1 :   ActionShortcut(ao) {
      91             :   // Create the matrices
      92             :   std::string sw_input;
      93           2 :   parse("SWITCH",sw_input);
      94             :   std::string sp_lab, sp_laba;
      95           1 :   parse("SPECIES",sp_lab);
      96           1 :   parse("SPECIESA",sp_laba);
      97           1 :   std::string cmap_input = getShortcutLabel() + "_cmap: CONTACT_MATRIX";
      98           1 :   if( sp_lab.length()>0 ) {
      99           2 :     readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUP=" + sp_lab + " COMPONENTS SWITCH={" + sw_input + "}");
     100           0 :   } else if( sp_laba.length()>0 ) {
     101             :     std::string sp_labb;
     102           0 :     parse("SPECIESB",sp_labb);
     103           0 :     readInputLine( getShortcutLabel() + "_cmap: CONTACT_MATRIX GROUPA=" + sp_laba + " GROUPB=" + sp_labb + " COMPONENTS SWITCH={" + sw_input + "}");
     104             :   }
     105             :   // Now need the Gaussians
     106             :   std::string mykernels;
     107           1 :   for(unsigned i=1;; ++i) {
     108             :     std::string kstr_inpt, istr, kern_str;
     109           3 :     Tools::convert( i, istr );
     110           6 :     if( !parseNumbered("KERNEL",i,kstr_inpt ) ) {
     111             :       break;
     112             :     }
     113           2 :     std::vector<std::string> words = Tools::getWords(kstr_inpt);
     114           2 :     if( words[0]=="GAUSSIAN" ) {
     115             :       kern_str="gaussian";
     116             :     } else {
     117           0 :       error("unknown kernel type");
     118             :     }
     119             :     std::string center, var;
     120           2 :     Tools::parse(words,"CENTER",center);
     121           4 :     Tools::parse(words,"SIGMA",var);
     122           2 :     if( mykernels.length()==0 ) {
     123           2 :       mykernels = "exp(-(ajik-" + center + ")^2/(2*" + var + "*" + var + "))";
     124             :     } else {
     125           2 :       mykernels = mykernels + "+exp(-(ajik-" + center + ")^2/(2*" + var + "*" + var + "))";
     126             :     }
     127           4 :   }
     128             :   // Hard coded switching function on minimum distance here -- should be improved
     129           3 :   readInputLine( getShortcutLabel() + "_ksum: GSYMFUNC_THREEBODY WEIGHT=" + getShortcutLabel() + "_cmap.w " +
     130           3 :                  "ARG=" + getShortcutLabel() + "_cmap.x," + getShortcutLabel() + "_cmap.y," + getShortcutLabel() + "_cmap.z"
     131           2 :                  " FUNCTION1={FUNC=" + mykernels + " LABEL=n} FUNCTION2={FUNC=1 LABEL=d}" );
     132             :   // And just the sum of the coordination numbers
     133           1 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_cmap");
     134           1 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     135             :   std::string size;
     136           1 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     137           2 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     138           2 :   readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cmap.w," + getShortcutLabel() + "_ones");
     139             :   // And the transformed switching functions
     140             :   std::string swcoord_str;
     141           1 :   parse("SWITCH_COORD",swcoord_str);
     142           2 :   readInputLine( getShortcutLabel() + "_mtdenom: MORE_THAN ARG=" + getShortcutLabel() + "_denom SWITCH={" + swcoord_str +"}");
     143             :   // And matheval to get the final quantity
     144           2 :   readInputLine( getShortcutLabel() + "_smac: CUSTOM ARG=" + getShortcutLabel() + "_ksum.n," + getShortcutLabel() + "_mtdenom," + getShortcutLabel() + "_ksum.d FUNC=x*y/z PERIODIC=NO");
     145             :   // And this expands everything
     146           2 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel() + "_smac", "", this );
     147           1 : }
     148             : 
     149             : }
     150             : }