Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 "tools/SwitchingFunction.h"
24 : #include "core/ActionRegister.h"
25 :
26 : #include <string>
27 : #include <cmath>
28 :
29 : namespace PLMD {
30 : namespace adjmat {
31 :
32 : //+PLUMEDOC MCOLVAR BRIDGE_MATRIX
33 : /*
34 : Calculate the number of atoms that bridge two parts of a structure
35 :
36 : This adjacency matrix is used to implement the [BRIDGE](BRIDGE.md) shortcut. The action outputs a adjacency matrix
37 : in the same way as [CONTACT_MATRIX](CONTACT_MATRIX.md). However, the $j,k$ element of the adjacency matrix is calculated
38 : using:
39 :
40 : $$
41 : M_{jk} = \sum_i s_A(r_{ij})s_B(r_{ik})
42 : $$
43 :
44 : In this expression, the sum runs over all the atoms that were specified using the `BRIDGING_ATOMS` keyword, $s_A$ and
45 : $s_B$ are switching functions, and $r_{ij}$ and $r_{ik}$ are the distances between atom $i$ and $j$ and between atoms
46 : $i$ and $k$. Less formally, this formula ensures that $j,k$ element of the output matrix is one if there is a bridging
47 : atom between atom $j$ and $k$.
48 :
49 : In the following example input atoms 100-200 can serve as bridging atoms between the atoms in GROUPA and GROUPB and the
50 : two switching functions $s_A$ and $s_B$ in the formula above are identical.
51 :
52 : ```plumed
53 : w1: BRIDGE_MATRIX ...
54 : BRIDGING_ATOMS=100-200
55 : GROUPA=1-10 GROUPB=11-20
56 : SWITCH={RATIONAL R_0=0.2}
57 : ...
58 : ```
59 :
60 : If you use a single GROUP keyword as in the input below as a single SWITCH keyword the output matrix is symmetric.
61 :
62 : ```plumed
63 : w2: BRIDGE_MATRIX ...
64 : BRIDGING_ATOMS=100-200 GROUP=1-10
65 : SWITCH={RATIONAL R_0=0.2}
66 : ...
67 : ```
68 :
69 : However, if the two switching functions are not identical, as in the following example, then the output matrix is __not__ symmetric
70 : even if GROUP is used rather than GROUPA/GROUPB.
71 :
72 : ```plumed
73 : w2: BRIDGE_MATRIX ...
74 : BRIDGING_ATOMS=100-200 GROUP=1-10
75 : SWITCHA={RATIONAL R_0=0.2}
76 : SWITCHB={RATIONAL R_0=0.4}
77 : ...
78 : ```
79 :
80 : Notice that in all the inputs above the $r_{ij}$ and $r_{ik}$ values that enter the formula above are calculated in a way that takes the
81 : periodic boundary conditions into account. If you want to ignore the periodic boundary conditions you can use the NOPBC flag as shown below.
82 :
83 : ```plumed
84 : w2: BRIDGE_MATRIX ...
85 : BRIDGING_ATOMS=100-200 GROUP=1-10
86 : SWITCH={RATIONAL R_0=0.2}
87 : NOPBC
88 : ...
89 : ```
90 :
91 : ## COMPONENTS flag
92 :
93 : If you add the flag COMPONENTS to the input as shown below:
94 :
95 : ```plumed
96 : c4: BRIDGE_MATRIX ...
97 : BRIDGING_ATOMS=100-200 GROUP=1-10
98 : SWITCH={RATIONAL R_0=0.2}
99 : COMPONENTS
100 : ...
101 : ```
102 :
103 : then four matrices with the labels `c4.w`, `c4.x`, `c4.y` and `c4.z` are output by the action. The matrix with the label `c4.w` is the adjacency matrix
104 : that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices `c4.x`, `c4.y` and `c4.z` contain the $x$, $y$ and $z$
105 : components of the vector connecting atoms $j$ and $k$. Importantly, however, the components of these vectors are only stored in `c4.x`, `c4.y` and `c4.z`
106 : if the elements of `c4.w` are non-zero. Using the COMPONENTS flag in this way ensures that you can use BRIDGE_MATRIX in tandem with many of the functionalities
107 : that are part of the [symfunc module](module_symfunc.md).
108 :
109 : ## The MASK keyword
110 :
111 : You use the MASK keyword with BRIDGE_MATRIX in the same way that is used in [CONTACT_MATRIX](CONTACT_MATRIX.md). This keyword thus expects a vector in input,
112 : which tells BRIDGE_MATRIX that it is safe to not calculate certain rows of the output matrix. An example where this keyword is used is shown below:
113 :
114 : ```plumed
115 : # The atoms that are of interest
116 : ow: GROUP ATOMS=1-1650
117 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
118 : center: FIXEDATOM AT=2.5,2.5,2.5
119 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
120 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
121 : # Calculates cooordination numbers
122 : cmap: BRIDGE_MATRIX ...
123 : GROUP=ow BRIDGING_ATOMS=1650-3000
124 : SWITCH={GAUSSIAN D_0=0.32 R_0=0.01 D_MAX=0.34} MASK=sphere
125 : ...
126 : ones: ONES SIZE=1650
127 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
128 : # Multiply coordination numbers by sphere vector
129 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
130 : # Sum of coordination numbers for atoms that are in the sphere of interest
131 : numer: SUM ARG=prod PERIODIC=NO
132 : # Number of atoms that are in sphere of interest
133 : denom: SUM ARG=sphere PERIODIC=NO
134 : # Average coordination number for atoms in sphere of interest
135 : av: CUSTOM ARG=prod,sphere FUNC=x/y PERIODIC=NO
136 : # And print out final CV to a file
137 : PRINT ARG=av FILE=colvar STRIDE=1
138 : ```
139 :
140 : This input calculates the average for a type of coordination number for those atoms that are within a spherical region that is centered on the point
141 : $(2.5,2.5,2.5)$. The coordination number for the atoms that is evaluated here counts an atom $k$ in the coordination sphere of atom $j$ if there is a
142 : bridging atom within 0.34 nm of atoms $j$ and $k$.
143 :
144 : */
145 : //+ENDPLUMEDOC
146 :
147 : class BridgeMatrix {
148 : public:
149 : SwitchingFunction sf1;
150 : SwitchingFunction sf2;
151 : static void registerKeywords( Keywords& keys );
152 : void parseInput( AdjacencyMatrixBase<BridgeMatrix>* action );
153 : static void calculateWeight( const BridgeMatrix& data,
154 : const AdjacencyMatrixInput& input,
155 : MatrixOutput output );
156 : };
157 :
158 : typedef AdjacencyMatrixBase<BridgeMatrix> bmap;
159 : PLUMED_REGISTER_ACTION(bmap,"BRIDGE_MATRIX")
160 :
161 18 : void BridgeMatrix::registerKeywords( Keywords& keys ) {
162 18 : keys.add("atoms","BRIDGING_ATOMS","The list of atoms that can form the bridge between the two interesting parts "
163 : "of the structure.");
164 18 : keys.add("optional","SWITCH","The parameters of the two switching functions in the above formula");
165 36 : keys.linkActionInDocs("SWITCH","LESS_THAN");
166 18 : keys.add("optional","SWITCHA","The switching function on the distance between bridging atoms and the atoms in "
167 : "group A");
168 36 : keys.linkActionInDocs("SWITCHA","LESS_THAN");
169 18 : keys.add("optional","SWITCHB","The switching function on the distance between the bridging atoms and the atoms in "
170 : "group B");
171 36 : keys.linkActionInDocs("SWITCHB","LESS_THAN");
172 18 : }
173 :
174 8 : void BridgeMatrix::parseInput( AdjacencyMatrixBase<BridgeMatrix>* action ) {
175 : bool oneswitch;
176 : std::string errors;
177 : std::string sf1input;
178 16 : action->parse("SWITCH",sf1input);
179 8 : if( sf1input.length()>0 ) {
180 8 : sf1.set(sf1input,errors);
181 8 : oneswitch=true;
182 8 : if( errors.length()!=0 ) {
183 0 : action->error("problem reading SWITCH keyword : " + errors );
184 : }
185 8 : sf2.set(sf1input,errors);
186 8 : if( errors.length()!=0 ) {
187 0 : action->error("problem reading SWITCH keyword : " + errors );
188 : }
189 : } else {
190 0 : action->parse("SWITCHA",sf1input);
191 0 : if(sf1input.length()>0) {
192 0 : sf1.set(sf1input,errors);
193 0 : oneswitch=false;
194 0 : if( errors.length()!=0 ) {
195 0 : action->error("problem reading SWITCHA keyword : " + errors );
196 : }
197 0 : std::string sf2input=sf1input;
198 0 : action->parse("SWITCHB",sf2input);
199 0 : if(sf2input.length()==0) {
200 0 : action->error("found SWITCHA keyword without SWITCHB");
201 : }
202 0 : sf2.set(sf2input,errors);
203 0 : if( errors.length()!=0 ) {
204 0 : action->error("problem reading SWITCHB keyword : " + errors );
205 : }
206 : } else {
207 0 : action->error("missing definition of switching functions");
208 : }
209 : }
210 8 : action->log.printf(" distance between bridging atoms and atoms in GROUPA must be less than %s\n",sf1.description().c_str());
211 8 : action->log.printf(" distance between bridging atoms and atoms in GROUPB must be less than %s\n",sf2.description().c_str());
212 :
213 : // Setup link cells
214 8 : action->setLinkCellCutoff( oneswitch, sf1.get_dmax() + sf2.get_dmax() );
215 8 : }
216 :
217 1037 : void BridgeMatrix::calculateWeight( const BridgeMatrix& data,
218 : const AdjacencyMatrixInput& input,
219 : MatrixOutput output ) {
220 1037 : output.val[0] = 0;
221 1037 : if( input.pos.modulo2()<epsilon ) {
222 : return;
223 : }
224 26859 : for(unsigned i=0; i<input.natoms; ++i) {
225 : Vector dij(input.extra_positions[i][0],
226 : input.extra_positions[i][1],
227 25822 : input.extra_positions[i][2]);
228 : double dijm = dij.modulo2();
229 25822 : double dw1, w1=data.sf1.calculateSqr( dijm, dw1 );
230 25822 : if( dijm<epsilon ) {
231 : w1=0.0;
232 0 : dw1=0.0;
233 : }
234 25822 : Vector dik=input.pbc->distance( dij, input.pos );
235 : double dikm=dik.modulo2();
236 25822 : double dw2, w2=data.sf2.calculateSqr( dikm, dw2 );
237 25822 : if( dikm<epsilon ) {
238 : w2=0.0;
239 0 : dw2=0.0;
240 : }
241 :
242 25822 : output.val[0] += w1*w2;
243 : // And finish the calculation
244 25822 : output.deriv[0] += -w2*dw1*dij[0];
245 25822 : output.deriv[1] += -w2*dw1*dij[1];
246 25822 : output.deriv[2] += -w2*dw1*dij[2];
247 25822 : output.deriv[3] += w1*dw2*dik[0];
248 25822 : output.deriv[4] += w1*dw2*dik[1];
249 25822 : output.deriv[5] += w1*dw2*dik[2];
250 25822 : output.deriv[6+i*3+0] = -w1*dw2*dik[0] + w2*dw1*dij[0];
251 25822 : output.deriv[6+i*3+1] = -w1*dw2*dik[1] + w2*dw1*dij[1];
252 25822 : output.deriv[6+i*3+2] = -w1*dw2*dik[2] + w2*dw1*dij[2];
253 51644 : Tensor vir = w1*(-dw2)*Tensor(dik,dik)+w2*(-dw1)*Tensor(dij,dij);
254 25822 : output.deriv[6 + 3*input.natoms + 0] += vir[0][0];
255 25822 : output.deriv[6 + 3*input.natoms + 1] += vir[0][1];
256 25822 : output.deriv[6 + 3*input.natoms + 2] += vir[0][2];
257 25822 : output.deriv[6 + 3*input.natoms + 3] += vir[1][0];
258 25822 : output.deriv[6 + 3*input.natoms + 4] += vir[1][1];
259 25822 : output.deriv[6 + 3*input.natoms + 5] += vir[1][2];
260 25822 : output.deriv[6 + 3*input.natoms + 6] += vir[2][0];
261 25822 : output.deriv[6 + 3*input.natoms + 7] += vir[2][1];
262 25822 : output.deriv[6 + 3*input.natoms + 8] += vir[2][2];
263 : }
264 : }
265 :
266 : }
267 : }
|
