Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "AlignedMatrixBase.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Matrix.h"
25 :
26 : //+PLUMEDOC MATRIX ALIGNED_MATRIX
27 : /*
28 : Adjacency matrix in which two molecule are adjacent if they are within a certain cutoff and if they have the same orientation.
29 :
30 : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
31 : so called adjacency matrix. An adjacency matrix is an \f$N \times N\f$ matrix in which the \f$i\f$th, \f$j\f$th element tells you whether
32 : or not the \f$i\f$th and \f$j\f$th atoms/molecules from a set of \f$N\f$ atoms/molecules are adjacent or not. These matrices can then be further
33 : analysed using a number of other algorithms as is detailed in \cite tribello-clustering.
34 :
35 : For this action the elements of the adjacency matrix are calculated using:
36 :
37 : \f[
38 : a_{ij} = \sigma_1( |\mathbf{r}_{ij}| ) \sigma_2( \mathbf{v}_i . \mathbf{v}_j )
39 : \f]
40 :
41 : This form of adjacency matrix can only be used if the input species are objects that lie at a point in space and that have an orientation, \f$\mathbf{v}\f$.
42 : These orientations might represent the
43 : orientation of a molecule, which could be calculated using \ref MOLECULES or \ref PLANES, or it might be the complex vectors calculated using the
44 : Steinhardt parameters \ref Q3, \ref Q4 or \ref Q6. In the expression above \f$\mathbf{r}_{ij}\f$ is the vector connecting the points in space
45 : where objects \f$i\f$ and \f$j\f$ find themselves and \f$\sigma_1\f$ is a \ref switchingfunction that acts upon the magnitude of this vector.
46 : \f$\sigma_2\f$ is a second \ref switchingfunction that acts on the dot product of the directors of the vectors that define the orientations of
47 : objects \f$i\f$ and \f$j\f$.
48 :
49 : \par Examples
50 :
51 : The example input below is necessarily but gives you an idea of what can be achieved using this action.
52 : The orientations and positions of four molecules are defined using the \ref MOLECULES action as the position of the
53 : centeres of mass of the two atoms specified and the direction of the vector connecting the two atoms that were specified.
54 : A \f$4 \times 4\f$ matrix is then computed using the formula above. The \f$ij\f$-element of this matrix tells us whether
55 : or not atoms \f$i\f$ and \f$j\f$ are within 0.1 nm of each other and whether or not the dot-product of their orientation vectors
56 : is greater than 0.5. The sum of the rows of this matrix are then computed. The sums of the \f$i\f$th row of this matrix tells us how
57 : many of the molecules that are within the first coordination sphere of molecule \f$i\f$ have an orientation that is similar to that of
58 : molecule \f$i\f$. We thus calculate the number of these "coordination numbers" that are greater than 1.0 and output this quantity to a file.
59 :
60 : \verbatim
61 : m1: MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8
62 : mat: ALIGNED_MATRIX ATOMS=m1 SWITCH={RATIONAL R_0=0.1} ORIENTATION_SWITCH={RATIONAL R_0=0.1 D_MAX=0.5}
63 : rr: ROWSUMS MATRIX=mat MORE_THAN={RATIONAL D_0=1.0 R_0=0.1}
64 : PRINT ARG=rr.* FILE=colvar
65 : \endverbatim
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 : namespace PLMD {
71 : namespace adjmat {
72 :
73 0 : class ContactAlignedMatrix : public AlignedMatrixBase {
74 : private:
75 : Matrix<SwitchingFunction> sf;
76 : public:
77 : ///
78 : static void registerKeywords( Keywords& keys );
79 : ///
80 : explicit ContactAlignedMatrix(const ActionOptions&);
81 : void readOrientationConnector( const unsigned& i, const unsigned& j, const std::vector<std::string>& desc );
82 : double computeVectorFunction( const unsigned& iv, const unsigned& jv,
83 : const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
84 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const ;
85 : };
86 :
87 2523 : PLUMED_REGISTER_ACTION(ContactAlignedMatrix,"ALIGNED_MATRIX")
88 :
89 1 : void ContactAlignedMatrix::registerKeywords( Keywords& keys ) {
90 1 : AlignedMatrixBase::registerKeywords( keys );
91 1 : keys.add("numbered","ORIENTATION_SWITCH","A switching function that transforms the dot product of the input vectors.");
92 1 : }
93 :
94 0 : ContactAlignedMatrix::ContactAlignedMatrix( const ActionOptions& ao ):
95 : Action(ao),
96 0 : AlignedMatrixBase(ao)
97 : {
98 0 : unsigned nrows, ncols, ig; retrieveTypeDimensions( nrows, ncols, ig );
99 0 : sf.resize( nrows, ncols );
100 0 : parseConnectionDescriptions("ORIENTATION_SWITCH",false,0);
101 0 : }
102 :
103 0 : void ContactAlignedMatrix::readOrientationConnector( const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
104 0 : plumed_assert( desc.size()==1 ); std::string errors; sf(j,i).set(desc[0],errors);
105 0 : if( j!=i ) sf(i,j).set(desc[0],errors);
106 0 : log.printf(" vectors in %u th and %u th groups must have a dot product that is greater than %s \n",i+1,j+1,(sf(i,j).description()).c_str() );
107 0 : }
108 :
109 0 : double ContactAlignedMatrix::computeVectorFunction( const unsigned& iv, const unsigned& jv,
110 : const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
111 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
112 0 : double dot_df, dot=0; dconn.zero();
113 0 : for(unsigned k=2; k<vec1.size(); ++k) dot+=vec1[k]*vec2[k];
114 0 : double f_dot = sf(iv,jv).calculate( dot, dot_df );
115 0 : for(unsigned k=2; k<vec1.size(); ++k) { dvec1[k]=dot_df*vec2[k]; dvec2[k]=dot_df*vec1[k]; }
116 0 : return f_dot;
117 : }
118 :
119 : }
120 2523 : }
121 :
122 :
123 :
|
