Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) envsim 2023-2024 The code team
3 : (see the PEOPLE-envsim file at the root of the distribution for a list of names)
4 :
5 : This file is part of envsim code module.
6 :
7 : The envsim code module is free software: you can redistribute it and/or modify
8 : it under the terms of the GNU Lesser General Public License as published by
9 : the Free Software Foundation, either version 3 of the License, or
10 : (at your option) any later version.
11 :
12 : The envsim code module is distributed in the hope that it will be useful,
13 : but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : GNU Lesser General Public License for more details.
16 :
17 : You should have received a copy of the GNU Lesser General Public License
18 : along with the envsim code module. If not, see <http://www.gnu.org/licenses/>.
19 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
20 : #include "core/ActionShortcut.h"
21 : #include "core/ActionRegister.h"
22 : #include "core/ActionWithValue.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "tools/PDB.h"
26 : #include "multicolvar/MultiColvarShortcuts.h"
27 : #include <string>
28 : #include <cmath>
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace envsim {
34 :
35 : //+PLUMEDOC MCOLVAR ENVIRONMENTSIMILARITY
36 : /*
37 : Measure how similar the environment around atoms is to that found in some reference crystal structure.
38 :
39 : This CV was introduced in this article \cite Piaggi-JCP-2019.
40 : The starting point for the definition of the CV is the local atomic density around an atom.
41 : We consider an environment \f$\chi\f$ around this atom and we define the density by
42 : \f[
43 : \rho_{\chi}(\mathbf{r})=\sum\limits_{i\in\chi} \exp\left(- \frac{|\mathbf{r}_i-\mathbf{r}|^2} {2\sigma^2} \right),
44 : \f]
45 : where \f$i\f$ runs over the neighbors in the environment \f$\chi\f$, \f$\sigma\f$ is a broadening parameter, and \f$\mathbf{r}_i\f$ are the
46 : coordinates of the neighbors relative to the central atom.
47 : We now define a reference environment or template \f$\chi_0\f$ that contains \f$n\f$ reference positions \f$\{\mathbf{r}^0_1,...,\mathbf{r}^0_n\}\f$
48 : that describe, for instance, the nearest neighbors in a given lattice.
49 : \f$\sigma\f$ is set using the SIGMA keyword and \f$\chi_0\f$ is chosen with the CRYSTAL_STRUCTURE keyword.
50 : If only the SPECIES keyword is given then the atoms defined there will be the central and neighboring atoms.
51 : If instead the SPECIESA and SPECIESB keywords are given then SPECIESA determines the central atoms and SPECIESB the neighbors.
52 :
53 : The environments \f$\chi\f$ and \f$\chi_0\f$ are compared using the kernel,
54 : \f[
55 : k_{\chi_0}(\chi)= \int d\mathbf{r} \rho_{\chi}(\mathbf{r}) \rho_{\chi_0}(\mathbf{r}) .
56 : \f]
57 : Combining the two equations above and performing the integration analytically we obtain,
58 : \f[
59 : k_{\chi_0}(\chi)= \sum\limits_{i\in\chi} \sum\limits_{j\in\chi_0} \pi^{3/2} \sigma^3 \exp\left(- \frac{|\mathbf{r}_i-\mathbf{r}^0_j|^2} {4\sigma^2} \right).
60 : \f]
61 : The kernel is finally normalized,
62 : \f[
63 : \tilde{k}_{\chi_0}(\chi) = \frac{1}{n} \sum\limits_{i\in\chi} \sum\limits_{j\in\chi_0} \exp\left( - \frac{|\mathbf{r}_i-\mathbf{r}^0_j|^2} {4\sigma^2} \right),
64 : \f]
65 : such that \f$\tilde{k}_{\chi_0}(\chi_0) = 1\f$.
66 : The above kernel is computed for each atom in the SPECIES or SPECIESA keywords.
67 : This quantity is a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
68 : the average value for the atoms in your system, the number of atoms that have an \f$\tilde{k}_{\chi_0}\f$ value that is more that some target and
69 : so on.
70 :
71 : The kernel can be generalized to crystal structures described as a lattice with a basis of more than one atom.
72 : In this case there is more than one type of environment.
73 : We consider the case of \f$M\f$ environments \f$X = \chi_1,\chi_2,...,\chi_M\f$ and we define the kernel through a best match strategy:
74 : \f[
75 : \tilde{k}_X(\chi)= \frac{1}{\lambda} \log \left ( \sum\limits_{l=1}^{M}\exp \left (\lambda \: \tilde{k}_{\chi_l}(\chi) \right ) \right ).
76 : \f]
77 : For a large enough \f$\lambda\f$ this expression will select the largest \f$\tilde{k}_{\chi_l}(\chi)\f$ with \f$\chi_l \in X\f$.
78 : This approach can be used, for instance, to target the hexagonal closed packed (HCP keyword) or the diamond structure (DIAMOND keyword).
79 :
80 : The CRYSTAL_STRUCTURE keyword can take the values SC (simple cubic), BCC (body centered cubic), FCC (face centered cubic),
81 : HCP (hexagonal closed pack), DIAMOND (cubic diamond), and CUSTOM (user defined).
82 : All options follow the same conventions as in the [lattice command](https://lammps.sandia.gov/doc/lattice.html) of [LAMMPS](https://lammps.sandia.gov/).
83 : If a CRYSTAL_STRUCTURE other than CUSTOM is used, then the lattice constants have to be specified using the keyword LATTICE_CONSTANTS.
84 : One value has to be specified for SC, BCC, FCC, and DIAMOND and two values have to be set for HCP (a and c lattice constants in that order).
85 :
86 : If the CUSTOM option is used then the reference environments have to be specified by the user.
87 : The reference environments are specified in pdb files containing the distance vectors from the central atom to the neighbors.
88 : Make sure your PDB file is correctly formatted as explained \ref pdbreader "in this page"
89 : If only one reference environment is specified then the filename should be given as argument of the keyword REFERENCE.
90 : If instead several reference environments are given, then they have to be provided in separate pdb files and given as arguments of the
91 : keywords REFERENCE_1, REFERENCE_2, etc.
92 : If you have a reference crystal structure configuration you can use the [Environment Finder](https://github.com/PabloPiaggi/EnvironmentFinder) app to determine the reference environments that you should use.
93 :
94 : If multiple chemical species are involved in the calculation, it is possible to provide the atom types (names) both for atoms in the reference environments and in the simulation box.
95 : This information is provided in pdb files using the atom name field.
96 : The comparison between environments is performed taking into account whether the atom names match.
97 :
98 : \par Examples
99 :
100 : The following input calculates the ENVIRONMENTSIMILARITY kernel for 250 atoms in the system
101 : using the BCC atomic environment as target, and then calculates and prints the average value
102 : for this quantity.
103 :
104 : \plumedfile
105 : ENVIRONMENTSIMILARITY SPECIES=1-250 SIGMA=0.05 LATTICE_CONSTANTS=0.423 CRYSTAL_STRUCTURE=BCC MEAN LABEL=es
106 :
107 : PRINT ARG=es.mean FILE=COLVAR
108 : \endplumedfile
109 :
110 : The next example compares the environments of the 96 selected atoms with a user specified reference
111 : environment. The reference environment is contained in the env1.pdb file. Once the kernel is computed
112 : the average and the number of atoms with a kernel larger than 0.5 are computed.
113 :
114 : \plumedfile
115 : ENVIRONMENTSIMILARITY ...
116 : SPECIES=1-288:3
117 : SIGMA=0.05
118 : CRYSTAL_STRUCTURE=CUSTOM
119 : REFERENCE=env1.pdb
120 : LABEL=es
121 : MEAN
122 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
123 : ... ENVIRONMENTSIMILARITY
124 :
125 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
126 : \endplumedfile
127 :
128 : The next example is similar to the one above but in this case 4 reference environments are specified.
129 : Each reference environment is given in a separate pdb file.
130 :
131 : \plumedfile
132 : ENVIRONMENTSIMILARITY ...
133 : SPECIES=1-288:3
134 : SIGMA=0.05
135 : CRYSTAL_STRUCTURE=CUSTOM
136 : REFERENCE_1=env1.pdb
137 : REFERENCE_2=env2.pdb
138 : REFERENCE_3=env3.pdb
139 : REFERENCE_4=env4.pdb
140 : LABEL=es
141 : MEAN
142 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
143 : ... ENVIRONMENTSIMILARITY
144 :
145 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
146 : \endplumedfile
147 :
148 : The following examples illustrates the use of pdb files to provide information about different chemical species:
149 : \plumedfile
150 : ENVIRONMENTSIMILARITY ...
151 : SPECIES=1-6
152 : SIGMA=0.05
153 : CRYSTAL_STRUCTURE=CUSTOM
154 : REFERENCE=env.pdb
155 : LABEL=es
156 : MEAN
157 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
158 : ATOM_NAMES_FILE=atom-names.pdb
159 : ... ENVIRONMENTSIMILARITY
160 : \endplumedfile
161 : Here the file env.pdb is:
162 : \verbatim
163 : ATOM 1 O MOL 1 -2.239 -1.296 -0.917 1.00 0.00 O
164 : ATOM 2 O MOL 1 0.000 0.000 2.751 1.00 0.00 O
165 : \endverbatim
166 : where atoms are of type O, and the atom-names.pdb file is:
167 : \verbatim
168 : ATOM 1 O X 1 0.000 2.593 4.126 0.00 0.00 O
169 : ATOM 2 H X 1 0.000 3.509 3.847 0.00 0.00 H
170 : ATOM 3 H X 1 0.000 2.635 5.083 0.00 0.00 H
171 : ATOM 4 O X 1 0.000 2.593 11.462 0.00 0.00 O
172 : ATOM 5 H X 1 0.000 3.509 11.183 0.00 0.00 H
173 : ATOM 6 H X 1 0.000 2.635 12.419 0.00 0.00 H
174 : \endverbatim
175 : where atoms are of type O and H.
176 : In this case, all atoms are used as centers, but only neighbors of type O are taken into account.
177 :
178 : */
179 : //+ENDPLUMEDOC
180 :
181 : class EnvironmentSimilarity : public ActionShortcut {
182 : private:
183 : std::vector<std::pair<unsigned,Vector> > getReferenceEnvironment( const PDB& pdb, const std::vector<std::string>& anames, double& maxdist );
184 : public:
185 : static void registerKeywords( Keywords& keys );
186 : explicit EnvironmentSimilarity(const ActionOptions&);
187 : };
188 :
189 : PLUMED_REGISTER_ACTION(EnvironmentSimilarity,"ENVIRONMENTSIMILARITY")
190 :
191 28 : void EnvironmentSimilarity::registerKeywords( Keywords& keys ) {
192 28 : ActionShortcut::registerKeywords( keys );
193 56 : keys.add("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
194 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
195 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
196 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
197 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
198 : "coordination number more than four for example");
199 56 : keys.add("atoms-4","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
200 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
201 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
202 : "using the label of another multicolvar");
203 56 : keys.add("atoms-4","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
204 : "the documentation for that keyword");
205 56 : keys.add("compulsory","CRYSTAL_STRUCTURE","FCC","Targeted crystal structure. Options are: "
206 : "SC: simple cubic, "
207 : "BCC: body center cubic, "
208 : "FCC: face centered cubic, "
209 : "HCP: hexagonal closed pack, "
210 : "DIAMOND: cubic diamond, "
211 : "CUSTOM: user defined "
212 : " ");
213 56 : keys.add("compulsory","LATTICE_CONSTANTS","Lattice constants. Two comma separated values for HCP, "
214 : "one value for all other CRYSTAL_STRUCTURES.");
215 56 : keys.add("compulsory","SIGMA","0.1","the width to use for the gaussian kernels");
216 56 : keys.add("compulsory","LCUTOFF","0.0001","any atoms separated by less than this tolerance should be ignored");
217 56 : keys.add("optional","REFERENCE","PDB files with relative distances from central atom. Use this keyword if you are targeting a single reference environment.");
218 56 : keys.add("numbered","REFERENCE_","PDB files with relative distances from central atom. Each file corresponds to one template. Use these keywords if you are targeting more than one reference environment.");
219 56 : keys.add("compulsory","LAMBDA","100","Lambda parameter. This is only used if you have more than one reference environment");
220 56 : keys.add("compulsory","CUTOFF","3","how many multiples of sigma would you like to consider beyond the maximum distance in the environment");
221 56 : keys.add("optional","ATOM_NAMES_FILE","PDB file with atom names for all atoms in SPECIES. Atoms in reference environments will be compared only if atom names match.");
222 28 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
223 28 : keys.needsAction("GROUP");
224 28 : keys.needsAction("DISTANCE_MATRIX");
225 28 : keys.needsAction("ONES");
226 28 : keys.needsAction("CONSTANT");
227 28 : keys.needsAction("CUSTOM");
228 28 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
229 28 : keys.needsAction("COMBINE");
230 28 : }
231 :
232 10 : EnvironmentSimilarity::EnvironmentSimilarity(const ActionOptions&ao):
233 : Action(ao),
234 10 : ActionShortcut(ao) {
235 : std::string atomNamesFile;
236 10 : parse("ATOM_NAMES_FILE",atomNamesFile);
237 10 : PDB atomnamepdb;
238 10 : if( !atomNamesFile.empty() && !atomnamepdb.read(atomNamesFile,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
239 0 : error("missing input file " + atomNamesFile);
240 : }
241 :
242 10 : double maxdist=0;
243 10 : std::vector<std::string> allspec(1);
244 : std::string crystal_structure;
245 20 : parse("CRYSTAL_STRUCTURE", crystal_structure);
246 : std::vector<std::vector<std::pair<unsigned,Vector> > > environments;
247 10 : if( crystal_structure=="CUSTOM" ) {
248 5 : if( !atomNamesFile.empty() ) {
249 1 : allspec[0]=atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[0]);
250 1 : unsigned natoms=atomnamepdb.getPositions().size();
251 385 : for(unsigned i=0; i<natoms; ++i) {
252 : bool found=false;
253 576 : for(unsigned j=0; j<allspec.size(); ++j) {
254 575 : if( allspec[j]==atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i] ) ) {
255 : found=true;
256 : break;
257 : }
258 : }
259 384 : if( !found ) {
260 2 : allspec.push_back( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i]) );
261 : }
262 : }
263 : }
264 : std::string reffile;
265 10 : parse("REFERENCE",reffile);
266 5 : if( reffile.length()>0 ) {
267 2 : PDB pdb;
268 2 : pdb.read(reffile,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
269 2 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
270 2 : log.printf(" reading %d reference vectors from %s \n", environments[0].size(), reffile.c_str() );
271 2 : } else {
272 12 : for(unsigned int i=1;; i++) {
273 15 : PDB pdb;
274 30 : if( !parseNumbered("REFERENCE_",i,reffile) ) {
275 : break;
276 : }
277 12 : if( !pdb.read(reffile,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
278 0 : error("missing input file " + reffile );
279 : }
280 12 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
281 12 : log.printf(" Reference environment %d : reading %d reference vectors from %s \n", i, environments[i-1].size(), reffile.c_str() );
282 15 : }
283 : }
284 : } else {
285 : std::vector<double> lattice_constants;
286 10 : parseVector("LATTICE_CONSTANTS", lattice_constants);
287 5 : if (crystal_structure == "FCC") {
288 1 : if (lattice_constants.size() != 1) {
289 0 : error("Number of LATTICE_CONSTANTS arguments must be one for FCC");
290 : }
291 1 : environments.resize(1);
292 1 : environments[0].resize(12);
293 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.0)*lattice_constants[0] );
294 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.0)*lattice_constants[0] );
295 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.0)*lattice_constants[0] );
296 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.0)*lattice_constants[0] );
297 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,+0.5)*lattice_constants[0] );
298 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,-0.5)*lattice_constants[0] );
299 1 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,+0.5)*lattice_constants[0] );
300 1 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,-0.5)*lattice_constants[0] );
301 1 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,+0.5)*lattice_constants[0] );
302 1 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,-0.5)*lattice_constants[0] );
303 1 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,+0.5)*lattice_constants[0] );
304 1 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,-0.5)*lattice_constants[0] );
305 1 : maxdist = std::sqrt(2)*lattice_constants[0]/2.;
306 4 : } else if (crystal_structure == "SC") {
307 0 : if (lattice_constants.size() != 1) {
308 0 : error("Number of LATTICE_CONSTANTS arguments must be one for SC");
309 : }
310 0 : environments.resize(1);
311 0 : environments[0].resize(6);
312 0 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
313 0 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
314 0 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
315 0 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
316 0 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
317 0 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
318 0 : maxdist = lattice_constants[0];
319 4 : } else if( crystal_structure == "BCC") {
320 2 : if (lattice_constants.size() != 1) {
321 0 : error("Number of LATTICE_CONSTANTS arguments must be one for BCC");
322 : }
323 2 : environments.resize(1);
324 2 : environments[0].resize(14);
325 2 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.5)*lattice_constants[0] );
326 2 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,-0.5)*lattice_constants[0] );
327 2 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.5)*lattice_constants[0] );
328 2 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.5)*lattice_constants[0] );
329 2 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,-0.5)*lattice_constants[0] );
330 2 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.5)*lattice_constants[0] );
331 2 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,-0.5)*lattice_constants[0] );
332 2 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,-0.5)*lattice_constants[0] );
333 2 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
334 2 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
335 2 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
336 2 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
337 2 : environments[0][12] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
338 2 : environments[0][13] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
339 2 : maxdist = lattice_constants[0];
340 2 : } else if (crystal_structure == "HCP") {
341 1 : if (lattice_constants.size() != 2) {
342 0 : error("Number of LATTICE_CONSTANTS arguments must be two for HCP");
343 : }
344 1 : environments.resize(2);
345 1 : environments[0].resize(12);
346 1 : environments[1].resize(12);
347 : double sqrt3=std::sqrt(3);
348 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
349 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
350 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
351 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
352 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
353 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
354 1 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
355 1 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
356 1 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
357 1 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
358 1 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
359 1 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
360 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
361 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
362 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
363 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
364 1 : environments[1][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
365 1 : environments[1][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
366 1 : environments[1][6] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
367 1 : environments[1][7] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
368 1 : environments[1][8] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1] );
369 1 : environments[1][9] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
370 1 : environments[1][10] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
371 1 : environments[1][11] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1] );
372 1 : maxdist = lattice_constants[0];
373 1 : } else if (crystal_structure == "DIAMOND") {
374 1 : if (lattice_constants.size() != 1) {
375 0 : error("Number of LATTICE_CONSTANTS arguments must be one for DIAMOND");
376 : }
377 1 : environments.resize(2);
378 1 : environments[0].resize(4);
379 1 : environments[1].resize(4);
380 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
381 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
382 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
383 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
384 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
385 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
386 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
387 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
388 1 : maxdist = std::sqrt(3)*lattice_constants[0]/4.0;
389 : } else {
390 0 : error( crystal_structure + " is not a valid input for keyword CRYSTAL_STRUCTURE");
391 : }
392 : }
393 10 : std::string matlab = getShortcutLabel() + "_cmat";
394 : double cutoff, sig;
395 10 : parse("SIGMA",sig);
396 20 : parse("CUTOFF",cutoff);
397 : std::string lcutoff;
398 20 : parse("LCUTOFF",lcutoff);
399 : std::string sig2;
400 10 : Tools::convert( sig*sig, sig2 );
401 10 : std::vector<std::vector<std::string> > funcstr(environments.size());
402 : std::string str_cutoff;
403 10 : Tools::convert( maxdist + cutoff*sig, str_cutoff );
404 : std::string str_natoms, xpos, ypos, zpos;
405 10 : Tools::convert( environments[0].size(), str_natoms );
406 31 : for(unsigned j=0; j<environments.size(); ++j) {
407 21 : funcstr[j].resize( allspec.size() );
408 46 : for(unsigned k=0; k<allspec.size(); ++k) {
409 177 : for(unsigned i=0; i<environments[j].size(); ++i) {
410 152 : if( environments[j][i].first!=k ) {
411 16 : continue ;
412 : }
413 136 : Tools::convert( environments[j][i].second[0], xpos );
414 136 : Tools::convert( environments[j][i].second[1], ypos );
415 136 : Tools::convert( environments[j][i].second[2], zpos );
416 136 : if( i==0 ) {
417 42 : funcstr[j][k] = "FUNC=(step(w-" + lcutoff + ")*step(" + str_cutoff + "-w)/" + str_natoms + ")*(exp(-((x-" + xpos + ")^2+(y-" + ypos + ")^2+(z-" + zpos + ")^2)/(4*" + sig2 + "))";
418 : } else {
419 230 : funcstr[j][k] += "+exp(-((x-" + xpos + ")^2+(y-" + ypos + ")^2+(z-" + zpos + ")^2)/(4*" + sig2 + "))";
420 : }
421 : }
422 25 : if( funcstr[j][k].length()>0 ) {
423 : funcstr[j][k] += ")";
424 : } else {
425 : funcstr[j][k] ="FUNC=0";
426 : }
427 : }
428 : }
429 :
430 : // Create the constact matrix
431 : std::string sp_str, specA, specB;
432 10 : parse("SPECIES",sp_str);
433 10 : parse("SPECIESA",specA);
434 20 : parse("SPECIESB",specB);
435 10 : if( sp_str.length()>0 ) {
436 18 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUP=" + sp_str + " CUTOFF=" + str_cutoff );
437 18 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
438 : } else {
439 1 : if( specA.length()==0 ) {
440 0 : error("no atoms were specified use SPECIES or SPECIESA+SPECIESB");
441 : }
442 1 : if( specB.length()==0 ) {
443 0 : error("no atoms were specified for SPECIESB");
444 : }
445 2 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUPA=" + specA + " GROUPB=" + specB + " CUTOFF=" + str_cutoff );
446 2 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + specA );
447 : }
448 :
449 : // Make a vector containing all ones
450 10 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
451 10 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
452 : std::string size;
453 10 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
454 10 : if( allspec.size()==1 ) {
455 18 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
456 : } else {
457 1 : unsigned natoms=atomnamepdb.getPositions().size();
458 : unsigned firstneigh=0;
459 1 : if( sp_str.length()==0 ) {
460 1 : firstneigh = (av->copyOutput(0))->getShape()[0];
461 : }
462 3 : for(unsigned i=0; i<allspec.size(); ++i) {
463 2 : std::string onesstr="0";
464 2 : if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[firstneigh])==allspec[i] ) {
465 : onesstr = "1";
466 : }
467 576 : for(unsigned j=firstneigh+1; j<natoms; ++j) {
468 574 : if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[j])==allspec[i] ) {
469 : onesstr += ",1";
470 : } else {
471 : onesstr += ",0";
472 : }
473 : }
474 4 : readInputLine( getShortcutLabel() + "_ones_" + allspec[i] + ": CONSTANT VALUES=" + onesstr );
475 : }
476 : }
477 :
478 : std::string envargstr,varstr, maxfuncstr, lambda;
479 10 : if( funcstr.size()>1 ) {
480 10 : parse("LAMBDA",lambda);
481 : }
482 : // And now do the funcstr bit
483 31 : for(unsigned j=0; j<funcstr.size(); ++j) {
484 : std::string jnum;
485 21 : Tools::convert( j+1, jnum );
486 21 : if(j==0) {
487 10 : varstr = "v" + jnum;
488 20 : maxfuncstr = "(1/" + lambda + ")*log(exp(" + lambda + "*v1)";
489 20 : envargstr = getShortcutLabel() + "_env" + jnum;
490 : } else {
491 11 : varstr += ",v" + jnum;
492 22 : maxfuncstr += "+exp(" + lambda + "*v" + jnum + ")";
493 22 : envargstr += "," + getShortcutLabel() + "_env" + jnum;
494 : }
495 : // And coordination numbers
496 21 : if( allspec.size()>1 ) {
497 : std::string argnames;
498 12 : for(unsigned i=0; i<allspec.size(); ++i) {
499 16 : readInputLine( getShortcutLabel() + "_" + allspec[i] + "_matenv" + jnum + ": CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w PERIODIC=NO " + funcstr[j][i] );
500 16 : readInputLine( getShortcutLabel() + "_" + allspec[i] + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_" + allspec[i] + "_matenv" + jnum + "," + getShortcutLabel() + "_ones_" + allspec[i] );
501 8 : if( i==0 ) {
502 8 : argnames = getShortcutLabel() + "_" + allspec[i] + "_env" + jnum;
503 : } else {
504 8 : argnames += "," + getShortcutLabel() + "_" + allspec[i] + "_env" + jnum;
505 : }
506 : }
507 4 : if( funcstr.size()==1) {
508 0 : readInputLine( getShortcutLabel() + ": COMBINE PERIODIC=NO ARG=" + argnames );
509 : } else {
510 8 : readInputLine( getShortcutLabel() + "_env" + jnum + ": COMBINE PERIODIC=NO ARG=" + argnames );
511 : }
512 : } else {
513 34 : readInputLine( getShortcutLabel() + "_matenv" + jnum + ": CUSTOM ARG=" + matlab + ".x," + matlab + ".y," + matlab + ".z," + matlab + ".w VAR=x,y,z,w PERIODIC=NO " + funcstr[j][0] );
514 17 : if( funcstr.size()==1) {
515 10 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
516 : } else {
517 24 : readInputLine( getShortcutLabel() + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
518 : }
519 : }
520 : }
521 : // And get the maximum
522 10 : if( funcstr.size()>1 ) {
523 10 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + envargstr + " PERIODIC=NO VAR=" + varstr + " FUNC=" + maxfuncstr + ")" );
524 : }
525 : // Read in all the shortcut stuff
526 : std::map<std::string,std::string> keymap;
527 10 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
528 20 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
529 40 : }
530 :
531 14 : std::vector<std::pair<unsigned,Vector> > EnvironmentSimilarity::getReferenceEnvironment( const PDB& pdb, const std::vector<std::string>& anames, double& maxdist ) {
532 14 : unsigned natoms = pdb.getPositions().size();
533 14 : std::vector<std::pair<unsigned,Vector> > env( natoms );
534 78 : for(unsigned i=0; i<natoms; ++i) {
535 : unsigned identity=0;
536 80 : for(unsigned j=1; j<anames.size(); ++j) {
537 16 : if( pdb.getAtomName(pdb.getAtomNumbers()[i])==anames[j] ) {
538 : identity=j;
539 : break;
540 : }
541 : }
542 64 : env[i] = std::pair<unsigned,Vector>( identity, pdb.getPositions()[i] );
543 64 : double dist = env[i].second.modulo();
544 64 : if( dist>maxdist ) {
545 13 : maxdist = dist;
546 : }
547 : }
548 14 : return env;
549 : }
550 :
551 : }
552 : }
|
