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 "AdjacencyMatrixBase.h"
23 : #include "multicolvar/AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "tools/Matrix.h"
27 :
28 : //+PLUMEDOC MATRIX CONTACT_MATRIX
29 : /*
30 : Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
31 :
32 : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
33 : 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
34 : 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
35 : analysed using a number of other algorithms as is detailed in \cite tribello-clustering.
36 :
37 : For this action the elements of the contact matrix are calculated using:
38 :
39 : \f[
40 : a_{ij} = \sigma( |\mathbf{r}_{ij}| )
41 : \f]
42 :
43 : where \f$|\mathbf{r}_{ij}|\f$ is the magnitude of the vector connecting atoms \f$i\f$ and \f$j\f$ and where \f$\sigma\f$ is a \ref switchingfunction.
44 :
45 : \par Examples
46 :
47 : The input shown below calculates a \f$6 \times 6\f$ matrix whose elements are equal to one if atom \f$i\f$ and atom \f$j\f$ are within 0.3 nm
48 : of each other and which is zero otherwise. The columns in this matrix are then summed so as to give the coordination number for each atom.
49 : The final quantity output in the colvar file is thus the average coordination number.
50 :
51 : \verbatim
52 : aa: CONTACT_MATRIX ATOMS=1-6 SWITCH={EXP D_0=0.2 R_0=0.1 D_MAX=0.66}
53 : COLUMNSUMS MATRIX=mat MEAN LABEL=csums
54 : PRINT ARG=csums.* FILE=colvar
55 : \endverbatim
56 :
57 : */
58 : //+ENDPLUMEDOC
59 :
60 :
61 : namespace PLMD {
62 : namespace adjmat {
63 :
64 24 : class ContactMatrix : public AdjacencyMatrixBase {
65 : private:
66 : /// Number of types that are in rows
67 : unsigned ncol_t;
68 : /// switching function
69 : Matrix<SwitchingFunction> switchingFunction;
70 : public:
71 : /// Create manual
72 : static void registerKeywords( Keywords& keys );
73 : /// Constructor
74 : explicit ContactMatrix(const ActionOptions&);
75 : /// Create the ith, ith switching function
76 : void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc );
77 : /// This actually calculates the value of the contact function
78 : double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const ;
79 : /// This does nothing
80 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
81 : };
82 :
83 2535 : PLUMED_REGISTER_ACTION(ContactMatrix,"CONTACT_MATRIX")
84 :
85 13 : void ContactMatrix::registerKeywords( Keywords& keys ) {
86 13 : AdjacencyMatrixBase::registerKeywords( keys );
87 : keys.add("atoms","ATOMS","The list of atoms for which you would like to calculate the contact matrix. The atoms involved must be specified "
88 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
89 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
90 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
91 13 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
92 : keys.add("numbered","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
93 : "The following provides information on the \\ref switchingfunction that are available. "
94 13 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
95 : // I added these keywords so I can test the results I get for column and row sums against output from COORDINATIONNUMBERS
96 : /// These should never be used in production as I think they will be much slower than COORDINATIONNUMBERS
97 13 : keys.add("hidden","ATOMSA",""); keys.add("hidden","ATOMSB","");
98 13 : }
99 :
100 12 : ContactMatrix::ContactMatrix( const ActionOptions& ao ):
101 : Action(ao),
102 12 : AdjacencyMatrixBase(ao)
103 : {
104 : // Read in the atoms and setup the matrix
105 12 : readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
106 12 : unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ncol_t );
107 12 : switchingFunction.resize( nrows, ncols );
108 : // Read in the switching functions
109 12 : parseConnectionDescriptions("SWITCH",false,ncol_t);
110 :
111 : // Find the largest sf cutoff
112 12 : double sfmax=switchingFunction(0,0).get_dmax();
113 25 : for(unsigned i=0; i<switchingFunction.nrows(); ++i) {
114 28 : for(unsigned j=0; j<switchingFunction.ncols(); ++j) {
115 15 : double tsf=switchingFunction(i,j).get_dmax();
116 15 : if( tsf>sfmax ) sfmax=tsf;
117 : }
118 : }
119 : // And set the link cell cutoff
120 12 : setLinkCellCutoff( sfmax );
121 12 : }
122 :
123 14 : void ContactMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
124 14 : plumed_assert( id==0 && desc.size()==1 ); std::string errors; switchingFunction(j,i).set(desc[0],errors);
125 14 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
126 14 : if( j!=i) switchingFunction(i,j).set(desc[0],errors);
127 14 : log.printf(" %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
128 14 : }
129 :
130 55742 : double ContactMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
131 55742 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
132 55742 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) - ncol_t ).get_dmax2() ) return 1.0;
133 26858 : return 0.0;
134 : }
135 :
136 28884 : double ContactMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
137 28884 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
138 28887 : double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) - ncol_t ).calculate( distance.modulo(), dfunc );
139 :
140 28888 : if( !doNotCalculateDerivatives() ) {
141 28654 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
142 28653 : double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) - ncol_t ).calculate( distance.modulo(), dfunc );
143 28654 : addAtomDerivatives( 1, 0, (-dfunc)*distance, myatoms );
144 28654 : addAtomDerivatives( 1, 1, (+dfunc)*distance, myatoms );
145 28654 : myatoms.addBoxDerivatives( 1, (-dfunc)*Tensor(distance,distance) );
146 : }
147 28888 : return sw;
148 : }
149 :
150 : }
151 2523 : }
152 :
|
