Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2018 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/ActionRegister.h"
      23             : #include "core/ActionShortcut.h"
      24             : #include "core/ActionWithValue.h"
      25             : #include "CoordinationNumbers.h"
      26             : #include "multicolvar/MultiColvarShortcuts.h"
      27             : #include "core/PlumedMain.h"
      28             : #include "core/ActionSet.h"
      29             : 
      30             : //+PLUMEDOC MCOLVAR TETRA_RADIAL
      31             : /*
      32             : Calculate the radial tetra CV
      33             : 
      34             : This shortcut calculates a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html). The particular function that is being
      35             : evaluated for the coordination sphere here is as follows:
      36             : 
      37             : $$
      38             : s_i = 1 - \frac{\sum_{j=1}^4 r_{ij}^2 - z_i\sum_{j=1}^4 r_{ij}}{12 z_i^2} \qquad \textrm{where} \qquad z_i = \frac{1}{4} \sum_{j=1}^4 r_{ij}
      39             : $$
      40             : 
      41             : In this expression the 4 atoms in the sums over $j$ are the four atoms that are nearest to atom $i$ and $r_{ij}$ is the distance between atoms $i$ and $j$.
      42             : The CV is large if the four atoms nearest atom $i$ are arranged on the vertices of a regular tetrahedron
      43             : and small otherwise.  The following example shows how you can use this action to measure the degree of tetrahedral order in a system.
      44             : 
      45             : ```plumed
      46             : # Calculate a vector that contains 64 values for the symmetry function.
      47             : acv: TETRA_RADIAL SPECIES=1-64
      48             : # Sum the elements of the vector and calculate the mean value on the atoms from this sum.
      49             : acv_sum: SUM ARG=acv PERIODIC=NO
      50             : acv_mean: CUSTOM ARG=acv FUNC=x/64 PERIODIC=NO
      51             : # Print out the positions of the 64 atoms for which the symmetry function was calculated
      52             : # to an xyz file along with the values of the symmetry function
      53             : DUMPATOMS ATOMS=1-64 ARG=acv FILE=mcolv.xyz
      54             : # Print out the average value of the symmetry function and the sum of all the symmetry functions
      55             : PRINT ARG=acv_sum,acv_mean FILE=colvar
      56             : ```
      57             : 
      58             : In the input above we have only one type of atom.  If you want to measure the degree of tetrahedral order amongst the bonds from atoms of SPECIESA to atoms of SPECIESB
      59             : you use an input like the one shown below:
      60             : 
      61             : ```plumed
      62             : acv: TETRA_RADIAL SPECIESA=1-64 SPECIESB=65-128
      63             : acv_sum: SUM ARG=acv PERIODIC=NO
      64             : acv_mean: CUSTOM ARG=acv FUNC=x/64 PERIODIC=NO
      65             : DUMPATOMS ATOMS=1-64 ARG=acv FILE=mcolv.xyz
      66             : PRINT ARG=acv_sum,acv_mean FILE=colvar
      67             : ```
      68             : 
      69             : By default the vectors connecting atoms that appear in the expression above are calculated in a way that takes periodic boundary conditions into account.  If you want
      70             : not to take the periodic boundary conditions into account for any reason you use the NOPBC flag as shown below:
      71             : 
      72             : ```plumed
      73             : acv: TETRA_RADIAL SPECIES=1-64 NOPBC
      74             : acv_sum: SUM ARG=acv PERIODIC=NO
      75             : acv_mean: CUSTOM ARG=acv FUNC=x/64 PERIODIC=NO
      76             : DUMPATOMS ATOMS=1-64 ARG=acv FILE=mcolv.xyz
      77             : PRINT ARG=acv_sum,acv_mean FILE=colvar
      78             : ```
      79             : 
      80             : ## Deprecated syntax
      81             : 
      82             : More information on the deprecated keywords that are given below is available in the documentation for the [DISTANCES](DISTANCES.md) command.
      83             : 
      84             : */
      85             : //+ENDPLUMEDOC
      86             : 
      87             : namespace PLMD {
      88             : namespace symfunc {
      89             : 
      90             : class RadialTetra : public ActionShortcut {
      91             : public:
      92             :   static void registerKeywords( Keywords& keys );
      93             :   explicit RadialTetra(const ActionOptions&ao);
      94             : };
      95             : 
      96             : PLUMED_REGISTER_ACTION(RadialTetra,"TETRA_RADIAL")
      97             : 
      98           3 : void RadialTetra::registerKeywords( Keywords& keys ) {
      99           3 :   CoordinationNumbers::shortcutKeywords( keys );
     100           6 :   keys.remove("MASK");
     101           3 :   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
     102           3 :   keys.add("compulsory","CUTOFF","-1","ignore distances that have a value larger than this cutoff");
     103           6 :   keys.setValueDescription("vector","the value of the radial tetrahedrality parameter for each of the input atoms");
     104           3 :   keys.remove("NN");
     105           3 :   keys.remove("MM");
     106           3 :   keys.remove("D_0");
     107           3 :   keys.remove("R_0");
     108           3 :   keys.remove("SWITCH");
     109           3 :   keys.needsAction("DISTANCE_MATRIX");
     110           3 :   keys.needsAction("NEIGHBORS");
     111           3 :   keys.needsAction("CUSTOM");
     112           3 :   keys.needsAction("ONES");
     113           3 :   keys.needsAction("MATRIX_VECTOR_PRODUCT");
     114           3 : }
     115             : 
     116           1 : RadialTetra::RadialTetra( const ActionOptions& ao):
     117             :   Action(ao),
     118           1 :   ActionShortcut(ao) {
     119             :   // Read species input and create the matrix
     120             :   bool nopbc;
     121           1 :   parseFlag("NOPBC",nopbc);
     122           1 :   std::string pbcstr="";
     123           1 :   if( nopbc ) {
     124             :     pbcstr = " NOPBC";
     125             :   }
     126             :   std::string sp_str, rcut;
     127           1 :   parse("SPECIES",sp_str);
     128           2 :   parse("CUTOFF",rcut);
     129           1 :   if( sp_str.length()>0 ) {
     130           2 :     readInputLine( getShortcutLabel() + "_mat: DISTANCE_MATRIX GROUP=" + sp_str + " CUTOFF=" + rcut + pbcstr );
     131             :   } else {
     132             :     std::string specA, specB;
     133           0 :     parse("SPECIESA",specA);
     134           0 :     parse("SPECIESB",specB);
     135           0 :     if( specA.length()==0 ) {
     136           0 :       error("missing input atoms");
     137             :     }
     138           0 :     if( specB.length()==0 ) {
     139           0 :       error("missing SPECIESB keyword");
     140             :     }
     141           0 :     readInputLine( getShortcutLabel() + "_mat: DISTANCE_MATRIX GROUPA=" + specA + " GROUPB=" + specB + " CUTOFF=" + rcut + pbcstr);
     142             :   }
     143             :   // Get the neighbors matrix
     144           2 :   readInputLine( getShortcutLabel() + "_neigh: NEIGHBORS ARG=" + getShortcutLabel() + "_mat NLOWEST=4");
     145             :   // Now get distance matrix that just contains four nearest distances
     146           2 :   readInputLine( getShortcutLabel() + "_near4: CUSTOM ARG=" + getShortcutLabel() + "_mat," + getShortcutLabel() + "_neigh MASK=" + getShortcutLabel() + "_neigh FUNC=x*y PERIODIC=NO");
     147             :   //Now compute sum of four nearest distances
     148           1 :   ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat");
     149           1 :   plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
     150             :   std::string size;
     151           1 :   Tools::convert( (av->copyOutput(0))->getShape()[1], size );
     152           2 :   readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
     153           2 :   readInputLine( getShortcutLabel() + "_sum4: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_near4," + getShortcutLabel() + "_ones");
     154             :   // Now compute squares of four nearest distance
     155           2 :   readInputLine( getShortcutLabel() + "_near4_2: CUSTOM ARG=" + getShortcutLabel() + "_near4 FUNC=x*x PERIODIC=NO");
     156             :   // Now compute sum of the squares of the four nearest distances
     157           2 :   readInputLine( getShortcutLabel() + "_sum4_2: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_near4_2," + getShortcutLabel() + "_ones");
     158             :   // Evaluate the average distance to the four nearest neighbors
     159           2 :   readInputLine( getShortcutLabel() + "_meanr: CUSTOM ARG=" + getShortcutLabel() + "_sum4 FUNC=0.25*x PERIODIC=NO");
     160             :   // Now evaluate the actual per atom CV
     161           2 :   readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_sum4," + getShortcutLabel() + "_sum4_2," + getShortcutLabel() + "_meanr " +
     162             :                  "FUNC=(1-(y-x*z)/(12*z*z)) PERIODIC=NO");
     163             :   // And get the things to do with the quantities we have computed
     164             :   std::map<std::string,std::string> keymap;
     165           1 :   multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
     166           2 :   multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
     167           1 : }
     168             : 
     169             : }
     170             : }