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