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 "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/Group.h"
27 : #include "ContactMatrix.h"
28 :
29 : //+PLUMEDOC MATRIX CONTACT_MATRIX
30 : /*
31 : Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
32 :
33 : A useful tool for developing complex collective variables is the notion of the
34 : so called adjacency matrix. An adjacency matrix can be an $N \times N$ matrix in which the $i$th, $j$th element tells you whether
35 : or not the $i$th and $j$th atoms/molecules from a set of $N$ atoms/molecules are adjacent or not. Alternatively, you can calculate
36 : an adjacency matrix between a set of $N$ atoms and a second set of $M$ atoms. For this type of matrix the $i$th, $j$th element tells you
37 : whether the whether the $i$th atom in the first group and the $j$th atom in the second group are adjacent or not. The adjacency matrix in this
38 : case is thus $N \times M$.
39 :
40 : The simplest type of adjacency matrix is a contact matrix. The elements of a contact matrix are calculated using:
41 :
42 : $$
43 : a_{ij} = \sigma( r_{ij} )
44 : $$
45 :
46 : where $r_{ij}$ is the magnitude of the vector connecting atoms $i$ and $j$ and where $\sigma$ is a switching function. If you want to calculate a
47 : contact matrix for one group of atoms the input would look something like this:
48 :
49 : ```plumed
50 : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
51 : ```
52 :
53 : Alternatively, if you want to calculate the contact matrix between two groups of atoms you would use an input like following:
54 :
55 : ```plumed
56 : c2: CONTACT_MATRIX GROUPA=1-7 GROUPB=8-20 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
57 : ```
58 :
59 : Once you have calculated a contact matrix you can perform various matrix operations by using the tools in the matrixtools or clusters modules.
60 :
61 : ## Coordination numbers
62 :
63 : If a contact matrix, $\mathbf{C}$, is multiplied from the back by a vector of ones the $i$th element of the resulting matrix tells you the number of atoms that are
64 : within $r_c$ of atom $i$. In other words, the coordination numbers of the atoms can be calculated from the contact matrix by doing matrix vector multiplication.
65 :
66 : This realisation about the relationship between the contact map and the coordination number is heavily used in PLUMED. For example, to calculate
67 : and print the coordination numbers of the first 7 atoms in the system with themselves you would use an input something like this:
68 :
69 : ```plumed
70 : c1: CONTACT_MATRIX GROUP=1-7 D_0=0.0 R_0=2.6 NN=6 MM=12
71 : ones: ONES SIZE=7
72 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
73 : PRINT ARG=cc FILE=colvar
74 : ```
75 :
76 : This input transforms the distances using the RATIONAL switching function that is discussed in the documentation for [LESS_THAN](LESS_THAN.md)
77 : The parameters for this type of switching function can be specified using `D_0`, `R_0`, `NN` and `MM` keywords as indicated above. However, we would recommend
78 : using the following syntax instead as this allows you to use the full range of switching function types. Most importantly, if you use the following syntax you can
79 : set the D_MAX parameter, as shown below. As discussed in the optimization details section below, __if you want good computational performance on large systems it is essential to set the D_MAX parameter.__
80 :
81 :
82 : ```plumed
83 : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12 D_MAX=5.0}
84 : ones: ONES SIZE=7
85 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
86 : PRINT ARG=cc FILE=colvar
87 : ```
88 :
89 : Implmenting the coordination number as a product of a matrix and a vector is useful as there are many different ways to define whether two atoms/molecules and to construct a "contact" matrix based on
90 : the result. For example:
91 :
92 : * You could say that two molecules are connected if they are within a certain distance of each other and if they have the same orientation (see [TORSIONS_MATRIX](TORSIONS_MATRIX.md)).
93 : * You could say that two water molecules are connected if they are hydrogen bonded to each other (see [HBOND_MATRIX](HBOND_MATRIX.md)).
94 : * You could say that two atoms are connected if they are within a certain distance of each other and if they have similar values for a CV (see [OUTER_PRODUCT](OUTER_PRODUCT.md)).
95 :
96 : When the coordination numbers is implemented in the way described above (by doing the matrix-vector multiplication) you have the flexibility to define the contact matrix that
97 : is used in the multiplication in whatever way you choose. In other words, this implementation of the coordination number is much more flexible. For example, suppose you want
98 : to calculate the number of atoms that have a coordination is greater than 3.0. You can do this with PLUMED using the following input:
99 :
100 : ```plumed
101 : # Calculate the contact matrix for the first seven atoms in the system
102 : c1: CONTACT_MATRIX GROUP=1-7 SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
103 : # Calculate the coordination numbers for the first seven atoms in the system
104 : ones: ONES SIZE=7
105 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
106 : # Set the ith element of the vector mtc equal to one if the coordination number of atom i is greater than 3.
107 : mtc: MORE_THAN ARG=cc SWITCH={RATIONAL D_0=3 R_0=1}
108 : # Calculate the number of atoms with a coordination number greater than 3.
109 : s: SUM ARG=mtc PERIODIC=NO
110 : PRINT ARG=s FILE=colvar
111 : ```
112 :
113 : Alternatively, consider the CV that was used to study perovskite nucleation in [this paper](https://pubs.acs.org/doi/abs/10.1021/acs.chemmater.9b04259). The CV here measures the number of
114 : methylamonium molecules that are attached to at least 5 other methylamoniusm molecules, at least 7 lead atoms and at least 11 ionide ions. We can calculate something akin to this
115 : CV and apply a restraint on the resulting quantity by using the following input file:
116 :
117 : ```plumed
118 : # Lead ions
119 : Pb: GROUP ATOMS=1-64
120 : # Iodide atoms
121 : I: GROUP ATOMS=65-256
122 : # Methylamonium "atoms" -- in the real CV these are centers of mass rather than single atoms
123 : cn: GROUP ATOMS=257-320
124 :
125 : ones64: ONES SIZE=64
126 : # Contact matrix that determines if methylamonium molecules are within 8 A of each other
127 : cm_cncn: CONTACT_MATRIX GROUP=cn SWITCH={RATIONAL R_0=0.8}
128 : # Coordination number of methylamounium with methylamonium
129 : cc_cncn: MATRIX_VECTOR_PRODUCT ARG=cm_cncn,ones64
130 : # Vector with elements that are one if coordiantion of methylamonium with methylamonium >5
131 : mt_cncn: MORE_THAN ARG=cc_cncn SWITCH={RATIONAL R_0=5 NN=12 MM=24}
132 :
133 : # Contact matrix that determines if methylamoinium moleulcule and Pb atom are within 7.5 A of each other
134 : cm_cnpb: CONTACT_MATRIX GROUPA=cn GROUPB=Pb SWITCH={RATIONAL R_0=0.75}
135 : # Coordination number of methylamonium with Pb
136 : cc_cnpb: MATRIX_VECTOR_PRODUCT ARG=cm_cnpb,ones64
137 : # Vector with elements that are one if coordination of methylamounium with lead is >7
138 : mt_cnpb: MORE_THAN ARG=cc_cnpb SWITCH={RATIONAL R_0=7 NN=12 MM=24}
139 :
140 : ones192: ONES SIZE=192
141 : # Contact matrix that determines if methylamoinium moleulcule and I atom are within 6.5 A of each other
142 : cm_cnI: CONTACT_MATRIX GROUPA=cn GROUPB=I SWITCH={RATIONAL R_0=0.65}
143 : # Coordination number of methylamonium with I
144 : cc_cnI: MATRIX_VECTOR_PRODUCT ARG=cm_cnI,ones192
145 : # Vector with elements that are one if coordination of methylamounium with lead is >11
146 : mt_cnI: MORE_THAN ARG=cc_cnI SWITCH={RATIONAL R_0=11 NN=12 MM=24}
147 :
148 : # Element wise product of these three input vectors.
149 : # mm[i]==1 if coordination number of corrsponding methylamounium with methylamonium is >5
150 : # and if coordination of methylamounium with Pb is >7 and if coordination of methylamounium with I > 11
151 : mm: CUSTOM ARG=mt_cncn,mt_cnpb,mt_cnI FUNC=x*y*z PERIODIC=NO
152 :
153 : # Sum of coordination numbers and thus equal to number of methylamoniums with desired coordination numbers
154 : ff: SUM ARG=mm PERIODIC=NO
155 :
156 : rr: RESTRAINT ARG=ff AT=62 KAPPA=10
157 : ```
158 :
159 : ## COMPONENTS flag
160 :
161 : If you add the flag COMPONENTS to the input as shown below:
162 :
163 : ```plumed
164 : c4: CONTACT_MATRIX GROUP=1-7 COMPONENTS SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
165 : ```
166 :
167 : 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
168 : 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$
169 : components of the vector connecting atoms $i$ and $j$. Importantly, however, the components of these vectors are only stored in `c4.x`, `c4.y` and `c4.z`
170 : if the elements of `c4.w` are non-zero.
171 :
172 : ## NOPBC flag
173 :
174 : By default PLUMED calculates the distances that are transformed by the switching function in the CONTACT_MATRIX action in a way that accounts for the periodic boundary
175 : conditions. If for any reason you do not want PLUMED to account for the periodic boundary conditions you can use the NOPBC flag as shown below:
176 :
177 : ```plumed
178 : c: CONTACT_MATRIX GROUP=1-7 NOPBC SWITCH={RATIONAL R_0=2.6 NN=6 MM=12}
179 : ```
180 :
181 : If you also add the COMPONENTS flag in this input then the vector connecting each pairs of atoms will also be calculated without taking the periodic boundary conditions into account.
182 :
183 : ## Optimisation details
184 :
185 : Adjacency matrices are sparse. Each atom is only be connected to a small number of neighbours and the vast majority of the elements of the contact matrix are thus zero. To reduce
186 : the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact matrix only the elements of the matrix that are
187 : non-zero are stored. The same thing holds for the additional values that are created when you use the COMPONENTS flag. The components of the vectors connecting atoms are only stored
188 : when the elements of `c4.w` are non-zero.
189 :
190 : We can also use the sparsity of the adjacency matrix to make the time required to compute a contact matrix scale linearly rather than quadratically with the number of atoms. Element
191 : $i,j$ of the contact matrix is only non-zero if two atoms are within a cutoff, $r_c$. We can determine that many pairs of atoms are further appart than $r_c$ without computing the
192 : distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. __To turn on these features you need to set the `D_MAX` parameter in the
193 : switching functions.__ The value you pass to the `D_MAX` keyword is used as the cutoff in the link cell algorithm.
194 :
195 : In theory we could further optimize the implementation of the CONTACT_MATRIX action by exploiting neighbor lists. If we were to do this we would likely add two further keywords as shown
196 : below:
197 :
198 : ```plumed
199 : c4: CONTACT_MATRIX GROUP=1-10000 SWITCH={RATIONAL R_0=0.1 NN=6 MM=12 D_MAX=0.5} NL_CUTOFF=0.7 NL_STRIDE=5
200 : ```
201 :
202 : The `NL_CUTOFF` keyword would be used to specify the cutoff (in nm) to use when constructing neighbor lists. This value would need to be slightly larger than the D_MAX parameter for the switching function.
203 : The `NL_STRIDE` keyword would then be used to specify how frequently the neighbour list should be updated. Thus far we have not found it necessary to implement this algorithm. We have been happy with the
204 : performance even if we use the linked list algorithm to update the neighbors on every step. If you feel that you need this CV to perform better please get in touch as adding a neighbor list for this action
205 : should be relatively straightforward.
206 :
207 : ## The MASK keyword
208 :
209 : Suppose that you want to calculate the average coordination number for the atoms that are within a sphere in the center of your simulation box. You can do so by exploiting an input similar to the one shown
210 : below:
211 :
212 : ```plumed
213 : # The atoms that are of interest
214 : ow: GROUP ATOMS=1-16500
215 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
216 : center: FIXEDATOM AT=2.5,2.5,2.5
217 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
218 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
219 : # Calculates cooordination numbers
220 : cmap: CONTACT_MATRIX GROUP=ow SWITCH={GAUSSIAN D_0=0.32 R_0=0.01 D_MAX=0.34}
221 : ones: ONES SIZE=16500
222 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
223 : # Multiply coordination numbers by sphere vector
224 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
225 : # Sum of coordination numbers for atoms that are in the sphere of interest
226 : numer: SUM ARG=prod PERIODIC=NO
227 : # Number of atoms that are in sphere of interest
228 : denom: SUM ARG=sphere PERIODIC=NO
229 : # Average coordination number for atoms in sphere of interest
230 : av: CUSTOM ARG=prod,sphere FUNC=x/y PERIODIC=NO
231 : # And print out final CV to a file
232 : PRINT ARG=av FILE=colvar STRIDE=1
233 : ```
234 :
235 : This calculation is slow because you have to calculate the coordination numbers of all the atoms even though only a small subset of these quanitties are required to compute the average coordination number in the
236 : sphere. To avoid all these unecessary calculations you use the `MASK` keyword as shown below:
237 :
238 : ```plumed
239 : # The atoms that are of interest
240 : ow: GROUP ATOMS=1-16500
241 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
242 : center: FIXEDATOM AT=2.5,2.5,2.5
243 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
244 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
245 : # Calculates cooordination numbers
246 : cmap: CONTACT_MATRIX GROUP=ow SWITCH={GAUSSIAN D_0=0.32 R_0=0.01 D_MAX=0.34} MASK=sphere
247 : ones: ONES SIZE=16500
248 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
249 : # Multiply coordination numbers by sphere vector
250 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
251 : # Sum of coordination numbers for atoms that are in the sphere of interest
252 : numer: SUM ARG=prod PERIODIC=NO
253 : # Number of atoms that are in sphere of interest
254 : denom: SUM ARG=sphere PERIODIC=NO
255 : # Average coordination number for atoms in sphere of interest
256 : av: CUSTOM ARG=prod,sphere FUNC=x/y PERIODIC=NO
257 : # And print out final CV to a file
258 : PRINT ARG=av FILE=colvar STRIDE=1
259 : ```
260 :
261 : Adding the instruction `MASK=sphere` to the CONTACT_MATRIX line in this input tells PLUMED to only calculate the $i$th row in the adjacency matrix if the $i$th element of the vector `sphere` is non-zero.
262 : In other words, by adding this command we have ensured that we are not calculating coordination numbers for atoms that are not in the sphere that is of interest. In this way we can thus reduce the computational
263 : expense of the calculation enormously.
264 :
265 : Notice, that there are other places where we can use this same trick. For example, we could have used MASK as shown below for our calculation of the CV for perovskite nucleation as shown below:
266 :
267 : ```plumed
268 : # Lead ions
269 : Pb: GROUP ATOMS=1-64
270 : # Iodide atoms
271 : I: GROUP ATOMS=65-256
272 : # Methylamonium "atoms" -- in the real CV these are centers of mass rather than single atoms
273 : cn: GROUP ATOMS=257-320
274 :
275 : ones192: ONES SIZE=192
276 : # Contact matrix that determines if methylamoinium moleulcule and I atom are within 6.5 A of each other
277 : cm_cnI: CONTACT_MATRIX GROUPA=cn GROUPB=I SWITCH={RATIONAL R_0=0.65 D_MAX=1.0}
278 : # Coordination number of methylamonium with I
279 : cc_cnI: MATRIX_VECTOR_PRODUCT ARG=cm_cnI,ones192
280 : # Vector with elements that are one if coordination of methylamounium with lead is >11
281 : mt_cnI: MORE_THAN ARG=cc_cnI SWITCH={RATIONAL R_0=11 NN=12 MM=24}
282 :
283 : ones64: ONES SIZE=64
284 : # Contact matrix that determines if methylamoinium moleulcule and Pb atom are within 7.5 A of each other
285 : cm_cnpb: CONTACT_MATRIX GROUPA=cn GROUPB=Pb SWITCH={RATIONAL R_0=0.75 D_MAX=1.5} MASK=mt_cnI
286 : # Coordination number of methylamonium with Pb
287 : cc_cnpb: MATRIX_VECTOR_PRODUCT ARG=cm_cnpb,ones64
288 : # Vector with elements that are one if coordination of methylamounium with lead is >7
289 : mt_cnpb: MORE_THAN ARG=cc_cnpb SWITCH={RATIONAL R_0=7 NN=12 MM=24}
290 :
291 : # Contact matrix that determines if methylamonium molecules are within 8 A of each other
292 : cm_cncn: CONTACT_MATRIX GROUP=cn SWITCH={RATIONAL R_0=0.8 D_MAX=1.75} MASK=mt_cnI
293 : # Coordination number of methylamounium with methylamonium
294 : cc_cncn: MATRIX_VECTOR_PRODUCT ARG=cm_cncn,ones64
295 : # Vector with elements that are one if coordiantion of methylamonium with methylamonium >5
296 : mt_cncn: MORE_THAN ARG=cc_cncn SWITCH={RATIONAL R_0=5 NN=12 MM=24}
297 :
298 : # Element wise product of these three input vectors.
299 : # mm[i]==1 if coordination number of corrsponding methylamounium with methylamonium is >5
300 : # and if coordination of methylamounium with Pb is >7 and if coordination of methylamounium with I > 11
301 : mm: CUSTOM ARG=mt_cncn,mt_cnpb,mt_cnI FUNC=x*y*z PERIODIC=NO
302 :
303 : # Sum of coordination numbers and thus equal to number of methylamoniums with desired coordination numbers
304 : ff: SUM ARG=mm PERIODIC=NO
305 :
306 : rr: RESTRAINT ARG=ff AT=62 KAPPA=10
307 : ```
308 :
309 : This trick works here because when we find that there are methylamoinium moleulcules with fewer than 11 lead atoms in there first coordination sphere we know that there
310 : is no point calculating the second two coordination numbers.
311 :
312 : */
313 : //+ENDPLUMEDOC
314 :
315 : namespace PLMD {
316 : namespace adjmat {
317 :
318 : class ContactMatrixShortcut : public ActionShortcut {
319 : public:
320 : static void registerKeywords(Keywords& keys);
321 : explicit ContactMatrixShortcut(const ActionOptions&);
322 : };
323 :
324 : PLUMED_REGISTER_ACTION(ContactMatrixShortcut,"CONTACT_MATRIX")
325 :
326 346 : void ContactMatrixShortcut::registerKeywords(Keywords& keys) {
327 346 : AdjacencyMatrixBase<ContactMatrix>::registerKeywords( keys );
328 346 : keys.remove("GROUP");
329 692 : keys.remove("SWITCH");
330 346 : keys.add("numbered","GROUP","specifies the list of atoms that should be assumed indistinguishable");
331 346 : keys.add("numbered","SWITCH","the input for the switching function that acts upon the distance between each pair of atoms");
332 692 : keys.linkActionInDocs("SWITCH","LESS_THAN");
333 346 : keys.addActionNameSuffix("_PROPER");
334 346 : keys.needsAction("TRANSPOSE");
335 346 : keys.needsAction("CONCATENATE");
336 346 : }
337 :
338 217 : ContactMatrixShortcut::ContactMatrixShortcut(const ActionOptions& ao):
339 : Action(ao),
340 217 : ActionShortcut(ao) {
341 : std::vector<std::string> grp_str;
342 217 : std::string atomsstr="";
343 : std::vector<std::string> atomsvec;
344 434 : parseVector("ATOMS",atomsvec);
345 217 : if( atomsvec.size()>0 ) {
346 0 : for(unsigned i=0; i<atomsvec.size(); ++i) {
347 0 : Group* gg = plumed.getActionSet().selectWithLabel<Group*>( atomsvec[i] );
348 0 : if( gg ) {
349 0 : grp_str.push_back( atomsvec[i] );
350 : }
351 : }
352 0 : if( grp_str.size()!=atomsvec.size() ) {
353 0 : grp_str.resize(0);
354 0 : atomsstr = " ATOMS=" + atomsvec[0];
355 0 : for(unsigned i=1; i<atomsvec.size(); ++i) {
356 0 : atomsstr += "," + atomsvec[i];
357 : }
358 : }
359 : } else {
360 : std::string grp_inpt;
361 2 : for(unsigned i=1;; ++i) {
362 438 : if( !parseNumbered("GROUP",i,grp_inpt) ) {
363 : break;
364 : }
365 2 : grp_str.push_back( grp_inpt );
366 : }
367 : }
368 217 : if( grp_str.size()>9 ) {
369 0 : error("cannot handle more than 9 groups");
370 : }
371 217 : if( grp_str.size()==0 ) {
372 432 : readInputLine( getShortcutLabel() + ": CONTACT_MATRIX_PROPER " + atomsstr + " " + convertInputLineToString() );
373 : return;
374 : }
375 :
376 3 : for(unsigned i=0; i<grp_str.size(); ++i) {
377 : std::string sw_str, num;
378 2 : Tools::convert( i+1, num );
379 4 : parseNumbered("SWITCH", (i+1)*10 + 1 + i, sw_str );
380 2 : if( sw_str.length()==0 ) {
381 0 : error("missing SWITCH" + num + num + " keyword");
382 : }
383 4 : readInputLine( getShortcutLabel() + num + num + ": CONTACT_MATRIX_PROPER GROUP=" + grp_str[i] + " SWITCH={" + sw_str + "}" );
384 3 : for(unsigned j=0; j<i; ++j) {
385 : std::string sw_str2, jnum;
386 1 : Tools::convert( j+1, jnum );
387 2 : parseNumbered("SWITCH", (j+1)*10 + 1 + i, sw_str2);
388 1 : if( sw_str2.length()==0 ) {
389 0 : error("missing SWITCH" + jnum + num + " keyword");
390 : }
391 2 : readInputLine( getShortcutLabel() + jnum + num + ": CONTACT_MATRIX_PROPER GROUPA=" + grp_str[j] + " GROUPB=" + grp_str[i] + " SWITCH={" + sw_str2 +"}");
392 2 : readInputLine( getShortcutLabel() + num + jnum + ": TRANSPOSE ARG=" + getShortcutLabel() + jnum + num );
393 : }
394 : }
395 1 : std::string join_matrices = getShortcutLabel() + ": CONCATENATE";
396 3 : for(unsigned i=0; i<grp_str.size(); ++i) {
397 : std::string inum;
398 2 : Tools::convert(i+1,inum);
399 6 : for(unsigned j=0; j<grp_str.size(); ++j) {
400 : std::string jnum;
401 4 : Tools::convert(j+1,jnum);
402 4 : if( i>j ) {
403 2 : join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum + jnum;
404 : } else {
405 6 : join_matrices += " MATRIX" + inum + jnum + "=" + getShortcutLabel() + inum + jnum;
406 : }
407 : }
408 : }
409 1 : readInputLine( join_matrices );
410 434 : }
411 :
412 : }
413 : }
|
