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 "AdjacencyMatrixBase.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/Matrix.h"
      25             : 
      26             : //+PLUMEDOC MATRIX DISTANCE_MATRIX
      27             : /*
      28             : Calculate a matrix of distances between atoms.
      29             : 
      30             : To calculate the matrix of distances between every distinct pair of atoms in a single group you can use the following command:
      31             : 
      32             : ```plumed
      33             : d1: DISTANCE_MATRIX GROUP=1-7
      34             : ```
      35             : 
      36             : If you would like to calculate the matrix of distances between the atoms in two different groups of atoms you can use the following command:
      37             : 
      38             : ```plumed
      39             : d2: DISTANCE_MATRIX GROUPA=1-7 GROUPB=8-20
      40             : ```
      41             : 
      42             : For both these inputs the distances between atoms are calculated in a way that takes the periodic boundary conditions into account. If you want to
      43             : ignore the periodic boundaries when calculating distances you use the NOPBC flag as shown below:
      44             : 
      45             : ```plumed
      46             : d3: DISTANCE_MATRIX GROUP=1-7 NOPBC
      47             : ```
      48             : 
      49             : Once you have calculated your distance matrix in this way you can do many of the operations that were discussed for [CONTACT_MATRIX](CONTACT_MATRIX.md) with the output.
      50             : For example, you can use the COMPONENTS flag to calcuate the $x$, $y$ and $z$ components of the vectors connecting the atoms in your two groups by using
      51             : an input like that shown below:
      52             : 
      53             : ```plumed
      54             : d1: DISTANCE_MATRIX GROUP=1-7 COMPONENTS
      55             : ```
      56             : 
      57             : ## Optimisation details
      58             : 
      59             : If for some reaon, you only want to calculate the distances if they are less than a certain CUTOFF you can add the cutoff keyword as follows:
      60             : 
      61             : ```plumed
      62             : d3: DISTANCE_MATRIX GROUP=1-7 CUTOFF=1.0
      63             : ```
      64             : 
      65             : This command will only store those distances that are less than 1 nm.  Notice, however, that the cutoff implemented here is __not__ a continuous function.
      66             : You should thus be careful when commands such as this one above to ensure that any quantities that are forced have continuous derivatives.  If you use
      67             : the CUTOFF keyword, however, many of the features that are used to optimise [CONTACT_MATRIX](CONTACT_MATRIX.md) are used for this action.
      68             : 
      69             : Also notice that you can use MASK to calculate a subset of the rows in the DISTANCE matrix as is done in the following example:
      70             : 
      71             : ```plumed
      72             : # The atoms that are of interest
      73             : ow: GROUP ATOMS=1-16500
      74             : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
      75             : center: FIXEDATOM AT=2.5,2.5,2.5
      76             : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
      77             : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
      78             : # The distance matrix
      79             : dmap: DISTANCE_MATRIX COMPONENTS GROUP=ow CUTOFF=1.0 MASK=sphere
      80             : # Find the four nearest neighbors
      81             : acv_neigh: NEIGHBORS ARG=dmap.w NLOWEST=4 MASK=sphere
      82             : # Compute a function for the atoms that are in the first coordination sphere
      83             : acv_g8: GSYMFUNC_THREEBODY ...
      84             :   WEIGHT=acv_neigh ARG=dmap.x,dmap.y,dmap.z
      85             :   FUNCTION1={FUNC=(cos(ajik)+1/3)^2 LABEL=g8}
      86             :   MASK=sphere
      87             : ...
      88             : # Now compute the value of the function above for those atoms that are in the
      89             : # sphere of interest
      90             : acv: CUSTOM ARG=acv_g8.g8,sphere FUNC=y*(1-(3*x/8)) PERIODIC=NO
      91             : # And now compute the final average
      92             : acv_sum: SUM ARG=acv PERIODIC=NO
      93             : acv_norm: SUM ARG=sphere PERIODIC=NO
      94             : mean: CUSTOM ARG=acv_sum,acv_norm FUNC=x/y PERIODIC=NO
      95             : PRINT ARG=mean FILE=colvar
      96             : ```
      97             : 
      98             : This input calculates the average value for measure of tetrahedral order that is introduced in the documentation for the [TETRA_ANGULAR](TETRA_ANGULAR.md) shortcut
      99             : for those atom that are within a sphere that is centered on the point $(2.5,2.5,2.5)$.
     100             : 
     101             : */
     102             : //+ENDPLUMEDOC
     103             : 
     104             : 
     105             : namespace PLMD {
     106             : namespace adjmat {
     107             : 
     108             : class DistanceMatrix {
     109             : public:
     110             :   double cutoff;
     111             :   static void registerKeywords( Keywords& keys );
     112             :   void parseInput( AdjacencyMatrixBase<DistanceMatrix>* action );
     113             :   static void calculateWeight( const DistanceMatrix& data,
     114             :                                const AdjacencyMatrixInput& input,
     115             :                                MatrixOutput output );
     116             : };
     117             : 
     118             : typedef AdjacencyMatrixBase<DistanceMatrix> dmap;
     119             : PLUMED_REGISTER_ACTION(dmap,"DISTANCE_MATRIX")
     120             : 
     121          57 : void DistanceMatrix::registerKeywords( Keywords& keys ) {
     122          57 :   keys.add("compulsory","CUTOFF","-1","ignore distances that have a value larger than this cutoff");
     123          57 : }
     124             : 
     125          33 : void DistanceMatrix::parseInput( AdjacencyMatrixBase<DistanceMatrix>* action ) {
     126             :   // And set the link cell cutoff
     127          33 :   action->log.printf("  weight is distance between atoms \n");
     128          33 :   action->parse("CUTOFF",cutoff);
     129          33 :   if( cutoff<0 ) {
     130           6 :     action->setLinkCellCutoff( true, std::numeric_limits<double>::max() );
     131             :   } else {
     132          27 :     action->log.printf("  ignoring distances that are larger than %f \n", cutoff);
     133          27 :     action->setLinkCellCutoff( true, cutoff );
     134             :   }
     135          33 : }
     136             : 
     137     7241483 : void DistanceMatrix::calculateWeight( const DistanceMatrix& data,
     138             :                                       const AdjacencyMatrixInput& input,
     139             :                                       MatrixOutput output ) {
     140     7241483 :   output.val[0] = input.pos.modulo();
     141     7241483 :   if( data.cutoff<0 || output.val[0]<data.cutoff ) {
     142     6023202 :     double invd = 1.0/output.val[0];
     143     6023202 :     Vector v = (-invd)*input.pos;
     144     6023202 :     output.deriv[0] = v[0];
     145     6023202 :     output.deriv[1] = v[1];
     146     6023202 :     output.deriv[2] = v[2];
     147     6023202 :     output.deriv[3] = -v[0];
     148     6023202 :     output.deriv[4] = -v[1];
     149     6023202 :     output.deriv[5] = -v[2];
     150             : 
     151     6023202 :     output.assignOuterProduct(6,v,input.pos);
     152             : 
     153             :   }
     154     7241483 : }
     155             : 
     156             : }
     157             : }
     158             :