Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 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/KernelFunctions.h"
25 : #include "tools/Torsion.h"
26 : #include "tools/Matrix.h"
27 :
28 : //+PLUMEDOC MATRIX SMAC_MATRIX
29 : /*
30 : Adjacency matrix in which two molecules are adjacent if they are within a certain cutoff and if the angle between them is within certain ranges.
31 :
32 : In this case the elements of the adjacency matrix are calculated using:
33 :
34 : \f[
35 : A_{ij} = \sigma(r_{ij}) \sum_n K_n(\theta_{ij})
36 : \f]
37 :
38 : In this expression \f$r_{ij}\f$ is the distance between molecule \f$i\f$ and molecule \f$j\f$ and \f$\sigma(r_{ij}\f$ is a
39 : \ref switchingfunction that acts on this distance. The $K_n functions are \ref kernelfunctions that take the torsion angle, \f$\theta_{ij}\f$, between the
40 : internal orientation vectors for molecules \f$i\f$ and \f$j\f$ as input. These kernel functions should be set so that they are
41 : equal to one when the relative orientation of the molecules are as they are in the solid and equal to zero otherwise.
42 : As the above matrix element is a product of functions it is only equal to one when the centers of mass of molecules \f$i\f$ and\f$j\f$
43 : are with a certain distance of each other and when the molecules are aligned in some desirable way.
44 :
45 : \par Examples
46 :
47 : In the following example an adjacency matrix is constructed in which the \f$(i,j)\f$ element is equal to one if
48 : molecules \f$i\f$ and \f$j\f$ are within 6 angstroms of each other and if the torsional angle between the orientations
49 : of these molecules is close to 0 or \f$\pi\f$. The various connected components of this matrix are determined using the
50 : \ref DFSCLUSTERING algorithm and then the size of the largest cluster of connects molecules is output to a colvar file
51 :
52 : \plumedfile
53 : UNITS LENGTH=A
54 :
55 : MOLECULES ...
56 : MOL1=1,2,1
57 : MOL2=5,6,5
58 : MOL3=9,10,9
59 : MOL4=13,14,13
60 : MOL5=17,18,17
61 : LABEL=m1
62 : ... MOLECULES
63 :
64 : SMAC_MATRIX ...
65 : ATOMS=m1 SWITCH={RATIONAL D_0=5.99 R_0=0.1 D_MAX=6.0}
66 : KERNEL1={TRIANGULAR CENTER=0 SIGMA=1.0} KERNEL2={TRIANGULAR CENTER=pi SIGMA=0.6}
67 : LABEL=smac
68 : ... SMAC_MATRIX
69 :
70 : dfs1: DFSCLUSTERING MATRIX=smac
71 : cc2: CLUSTER_NATOMS CLUSTERS=dfs1 CLUSTER=1
72 : PRINT ARG=cc2 FILE=colvar
73 : \endplumedfile
74 :
75 : */
76 : //+ENDPLUMEDOC
77 :
78 : namespace PLMD {
79 : namespace adjmat {
80 :
81 : class SMACMatrix : public AlignedMatrixBase {
82 : private:
83 : Matrix<std::vector<KernelFunctions> > kernels;
84 : public:
85 : ///
86 : static void registerKeywords( Keywords& keys );
87 : ///
88 : explicit SMACMatrix(const ActionOptions&);
89 : void readOrientationConnector( const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) override;
90 : double computeVectorFunction( const unsigned& iv, const unsigned& jv,
91 : const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
92 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const override;
93 : };
94 :
95 10421 : PLUMED_REGISTER_ACTION(SMACMatrix,"SMAC_MATRIX")
96 :
97 2 : void SMACMatrix::registerKeywords( Keywords& keys ) {
98 2 : AlignedMatrixBase::registerKeywords( keys );
99 4 : keys.add("numbered","KERNEL","The various kernels that are used to determine whether or not the molecules are aligned");
100 2 : }
101 :
102 1 : SMACMatrix::SMACMatrix( const ActionOptions& ao ):
103 : Action(ao),
104 1 : AlignedMatrixBase(ao)
105 : {
106 1 : unsigned nrows, ncols, ig; retrieveTypeDimensions( nrows, ncols, ig );
107 1 : kernels.resize( nrows, ncols ); parseConnectionDescriptions("KERNEL",true,0);
108 1 : }
109 :
110 1 : void SMACMatrix::readOrientationConnector( const unsigned& iv, const unsigned& jv, const std::vector<std::string>& desc ) {
111 3 : for(int i=0; i<desc.size(); i++) {
112 2 : KernelFunctions mykernel( desc[i] );
113 2 : kernels(iv,jv).push_back( mykernel );
114 2 : if( jv!=iv ) kernels(jv,iv).push_back( mykernel );
115 2 : }
116 1 : if( kernels(iv,jv).size()==0 ) error("no kernels defined");
117 1 : }
118 :
119 5806 : double SMACMatrix::computeVectorFunction( const unsigned& iv, const unsigned& jv,
120 : const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
121 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
122 :
123 5806 : unsigned nvectors = ( vec1.size() - 2 ) / 3; plumed_assert( (vec1.size()-2)%3==0 );
124 5806 : std::vector<Vector> dv1(nvectors), dv2(nvectors), tdconn(nvectors); Torsion t; std::vector<Vector> v1(nvectors), v2(nvectors);
125 : std::vector<std::unique_ptr<Value>> pos;
126 11612 : for(unsigned i=0; i<nvectors; ++i) { pos.emplace_back( Tools::make_unique<Value>() ); pos[i]->setDomain( "-pi", "pi" ); }
127 :
128 11612 : for(unsigned j=0; j<nvectors; ++j) {
129 23224 : for(unsigned k=0; k<3; ++k) {
130 17418 : v1[j][k]=vec1[2+3*j+k]; v2[j][k]=vec2[2+3*j+k];
131 : }
132 5806 : double angle = t.compute( v1[j], conn, v2[j], dv1[j], tdconn[j], dv2[j] );
133 : pos[j]->set( angle );
134 : }
135 :
136 5806 : double ans=0; std::vector<double> deriv( nvectors ), df( nvectors, 0 );
137 :
138 5806 : auto pos_ptr=Tools::unique2raw(pos);
139 :
140 17418 : for(unsigned i=0; i<kernels(iv,jv).size(); ++i) {
141 11612 : ans += kernels(iv,jv)[i].evaluate( pos_ptr, deriv );
142 23224 : for(unsigned j=0; j<nvectors; ++j) df[j] += deriv[j];
143 : }
144 11612 : dconn.zero(); for(unsigned j=0; j<nvectors; ++j) dconn += df[j]*tdconn[j];
145 11612 : for(unsigned j=0; j<nvectors; ++j) {
146 23224 : for(unsigned k=0; k<3; ++k) { dvec1[2+3*j+k]=df[j]*dv1[j][k]; dvec2[2+3*j+k]=df[j]*dv2[j][k]; }
147 : }
148 5806 : return ans;
149 5806 : }
150 :
151 : }
152 : }
153 :
154 :
155 :
|
