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 13787 : PLUMED_REGISTER_ACTION(SMACMatrix,"SMAC_MATRIX")
96 :
97 5 : void SMACMatrix::registerKeywords( Keywords& keys ) {
98 5 : AlignedMatrixBase::registerKeywords( keys );
99 10 : keys.add("numbered","KERNEL","The various kernels that are used to determine whether or not the molecules are aligned");
100 5 : }
101 :
102 1 : SMACMatrix::SMACMatrix( const ActionOptions& ao ):
103 : Action(ao),
104 1 : AlignedMatrixBase(ao) {
105 : unsigned nrows, ncols, ig;
106 1 : retrieveTypeDimensions( nrows, ncols, ig );
107 1 : kernels.resize( nrows, ncols );
108 1 : parseConnectionDescriptions("KERNEL",true,0);
109 1 : }
110 :
111 1 : void SMACMatrix::readOrientationConnector( const unsigned& iv, const unsigned& jv, const std::vector<std::string>& desc ) {
112 3 : for(int i=0; i<desc.size(); i++) {
113 2 : KernelFunctions mykernel( desc[i] );
114 2 : kernels(iv,jv).push_back( mykernel );
115 2 : if( jv!=iv ) {
116 0 : kernels(jv,iv).push_back( mykernel );
117 : }
118 : }
119 1 : if( kernels(iv,jv).size()==0 ) {
120 0 : error("no kernels defined");
121 : }
122 1 : }
123 :
124 5806 : double SMACMatrix::computeVectorFunction( const unsigned& iv, const unsigned& jv,
125 : const Vector& conn, const std::vector<double>& vec1, const std::vector<double>& vec2,
126 : Vector& dconn, std::vector<double>& dvec1, std::vector<double>& dvec2 ) const {
127 :
128 5806 : unsigned nvectors = ( vec1.size() - 2 ) / 3;
129 5806 : plumed_assert( (vec1.size()-2)%3==0 );
130 5806 : std::vector<Vector> dv1(nvectors), dv2(nvectors), tdconn(nvectors);
131 : Torsion t;
132 5806 : std::vector<Vector> v1(nvectors), v2(nvectors);
133 : std::vector<std::unique_ptr<Value>> pos;
134 11612 : for(unsigned i=0; i<nvectors; ++i) {
135 5806 : pos.emplace_back( Tools::make_unique<Value>() );
136 11612 : pos[i]->setDomain( "-pi", "pi" );
137 : }
138 :
139 11612 : for(unsigned j=0; j<nvectors; ++j) {
140 23224 : for(unsigned k=0; k<3; ++k) {
141 17418 : v1[j][k]=vec1[2+3*j+k];
142 17418 : v2[j][k]=vec2[2+3*j+k];
143 : }
144 5806 : double angle = t.compute( v1[j], conn, v2[j], dv1[j], tdconn[j], dv2[j] );
145 : pos[j]->set( angle );
146 : }
147 :
148 : double ans=0;
149 5806 : std::vector<double> deriv( nvectors ), df( nvectors, 0 );
150 :
151 5806 : auto pos_ptr=Tools::unique2raw(pos);
152 :
153 17418 : for(unsigned i=0; i<kernels(iv,jv).size(); ++i) {
154 11612 : ans += kernels(iv,jv)[i].evaluate( pos_ptr, deriv );
155 23224 : for(unsigned j=0; j<nvectors; ++j) {
156 11612 : df[j] += deriv[j];
157 : }
158 : }
159 5806 : dconn.zero();
160 11612 : for(unsigned j=0; j<nvectors; ++j) {
161 5806 : dconn += df[j]*tdconn[j];
162 : }
163 11612 : for(unsigned j=0; j<nvectors; ++j) {
164 23224 : for(unsigned k=0; k<3; ++k) {
165 17418 : dvec1[2+3*j+k]=df[j]*dv1[j][k];
166 17418 : dvec2[2+3*j+k]=df[j]*dv2[j][k];
167 : }
168 : }
169 5806 : return ans;
170 5806 : }
171 :
172 : }
173 : }
174 :
175 :
176 :
|
