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 ); keys.needsAction("GROUP");
223 84 : keys.needsAction("DISTANCE_MATRIX"); keys.needsAction("ONES"); keys.needsAction("CONSTANT");
224 84 : keys.needsAction("CUSTOM"); keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("COMBINE");
225 28 : }
226 :
227 10 : EnvironmentSimilarity::EnvironmentSimilarity(const ActionOptions&ao):
228 : Action(ao),
229 10 : ActionShortcut(ao)
230 : {
231 20 : std::string atomNamesFile; parse("ATOM_NAMES_FILE",atomNamesFile); PDB atomnamepdb;
232 10 : if( !atomNamesFile.empty() && !atomnamepdb.read(atomNamesFile,usingNaturalUnits(),0.1/getUnits().getLength()) ) error("missing input file " + atomNamesFile);
233 :
234 10 : double maxdist=0; std::vector<std::string> allspec(1);
235 20 : std::string crystal_structure; parse("CRYSTAL_STRUCTURE", crystal_structure);
236 : std::vector<std::vector<std::pair<unsigned,Vector> > > environments;
237 10 : if( crystal_structure=="CUSTOM" ) {
238 5 : if( !atomNamesFile.empty() ) {
239 2 : allspec[0]=atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[0]); unsigned natoms=atomnamepdb.getPositions().size();
240 385 : for(unsigned i=0; i<natoms; ++i) {
241 : bool found=false;
242 576 : for(unsigned j=0; j<allspec.size(); ++j) {
243 575 : if( allspec[j]==atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i] ) ) { found=true; break; }
244 : }
245 385 : if( !found ) allspec.push_back( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i]) );
246 : }
247 : }
248 10 : std::string reffile; parse("REFERENCE",reffile);
249 5 : if( reffile.length()>0 ) {
250 2 : PDB pdb; pdb.read(reffile,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
251 2 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
252 2 : log.printf(" reading %d reference vectors from %s \n", environments[0].size(), reffile.c_str() );
253 2 : } else {
254 12 : for(unsigned int i=1;; i++) {
255 30 : PDB pdb; if( !parseNumbered("REFERENCE_",i,reffile) ) break;
256 12 : if( !pdb.read(reffile,usingNaturalUnits(),0.1/getUnits().getLength()) ) error("missing input file " + reffile );
257 12 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
258 12 : log.printf(" Reference environment %d : reading %d reference vectors from %s \n", i, environments[i-1].size(), reffile.c_str() );
259 15 : }
260 : }
261 : } else {
262 10 : std::vector<double> lattice_constants; parseVector("LATTICE_CONSTANTS", lattice_constants);
263 5 : if (crystal_structure == "FCC") {
264 1 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for FCC");
265 1 : environments.resize(1); environments[0].resize(12);
266 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.0)*lattice_constants[0] );
267 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.0)*lattice_constants[0] );
268 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.0)*lattice_constants[0] );
269 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.0)*lattice_constants[0] );
270 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,+0.5)*lattice_constants[0] );
271 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,-0.5)*lattice_constants[0] );
272 1 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,+0.5)*lattice_constants[0] );
273 1 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,-0.5)*lattice_constants[0] );
274 1 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,+0.5)*lattice_constants[0] );
275 1 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,-0.5)*lattice_constants[0] );
276 1 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,+0.5)*lattice_constants[0] );
277 1 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,-0.5)*lattice_constants[0] );
278 1 : maxdist = std::sqrt(2)*lattice_constants[0]/2.;
279 4 : } else if (crystal_structure == "SC") {
280 0 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for SC");
281 0 : environments.resize(1); environments[0].resize(6);
282 0 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
283 0 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
284 0 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
285 0 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
286 0 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
287 0 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
288 0 : maxdist = lattice_constants[0];
289 4 : } else if( crystal_structure == "BCC") {
290 2 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for BCC");
291 2 : environments.resize(1); environments[0].resize(14);
292 2 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.5)*lattice_constants[0] );
293 2 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,-0.5)*lattice_constants[0] );
294 2 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.5)*lattice_constants[0] );
295 2 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.5)*lattice_constants[0] );
296 2 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,-0.5)*lattice_constants[0] );
297 2 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.5)*lattice_constants[0] );
298 2 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,-0.5)*lattice_constants[0] );
299 2 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,-0.5)*lattice_constants[0] );
300 2 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
301 2 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
302 2 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
303 2 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
304 2 : environments[0][12] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
305 2 : environments[0][13] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
306 2 : maxdist = lattice_constants[0];
307 2 : } else if (crystal_structure == "HCP") {
308 1 : if (lattice_constants.size() != 2) error("Number of LATTICE_CONSTANTS arguments must be two for HCP");
309 1 : environments.resize(2); environments[0].resize(12); environments[1].resize(12); double sqrt3=std::sqrt(3);
310 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
311 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
312 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
313 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
314 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
315 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
316 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] );
317 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] );
318 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] );
319 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] );
320 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] );
321 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] );
322 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
323 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
324 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
325 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
326 1 : environments[1][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
327 1 : environments[1][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
328 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] );
329 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] );
330 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] );
331 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] );
332 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] );
333 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] );
334 1 : maxdist = lattice_constants[0];
335 1 : } else if (crystal_structure == "DIAMOND") {
336 1 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for DIAMOND");
337 1 : environments.resize(2); environments[0].resize(4); environments[1].resize(4);
338 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
339 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
340 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
341 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
342 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
343 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
344 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
345 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
346 1 : maxdist = std::sqrt(3)*lattice_constants[0]/4.0;
347 0 : } else error( crystal_structure + " is not a valid input for keyword CRYSTAL_STRUCTURE");
348 : }
349 10 : std::string matlab = getShortcutLabel() + "_cmat";
350 40 : double cutoff, sig; parse("SIGMA",sig); parse("CUTOFF",cutoff); std::string lcutoff; parse("LCUTOFF",lcutoff);
351 10 : std::string sig2; Tools::convert( sig*sig, sig2 ); std::vector<std::vector<std::string> > funcstr(environments.size());
352 10 : std::string str_cutoff; Tools::convert( maxdist + cutoff*sig, str_cutoff );
353 10 : std::string str_natoms, xpos, ypos, zpos; Tools::convert( environments[0].size(), str_natoms );
354 31 : for(unsigned j=0; j<environments.size(); ++j) {
355 21 : funcstr[j].resize( allspec.size() );
356 46 : for(unsigned k=0; k<allspec.size(); ++k) {
357 177 : for(unsigned i=0; i<environments[j].size(); ++i) {
358 152 : if( environments[j][i].first!=k ) continue ;
359 136 : Tools::convert( environments[j][i].second[0], xpos ); Tools::convert( environments[j][i].second[1], ypos ); Tools::convert( environments[j][i].second[2], zpos );
360 157 : if( i==0 ) 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 + "))";
361 230 : else funcstr[j][k] += "+exp(-((x-" + xpos + ")^2+(y-" + ypos + ")^2+(z-" + zpos + ")^2)/(4*" + sig2 + "))";
362 : }
363 25 : if( funcstr[j][k].length()>0 ) funcstr[j][k] += ")"; else funcstr[j][k] ="FUNC=0";
364 : }
365 : }
366 :
367 : // Create the constact matrix
368 40 : std::string sp_str, specA, specB; parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB);
369 10 : if( sp_str.length()>0 ) {
370 18 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUP=" + sp_str + " CUTOFF=" + str_cutoff );
371 18 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
372 : } else {
373 1 : if( specA.length()==0 ) error("no atoms were specified use SPECIES or SPECIESA+SPECIESB");
374 1 : if( specB.length()==0 ) error("no atoms were specified for SPECIESB");
375 2 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUPA=" + specA + " GROUPB=" + specB + " CUTOFF=" + str_cutoff );
376 2 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + specA );
377 : }
378 :
379 : // Make a vector containing all ones
380 10 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
381 10 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
382 10 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size );
383 10 : if( allspec.size()==1 ) {
384 18 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
385 : } else {
386 1 : unsigned natoms=atomnamepdb.getPositions().size();
387 1 : unsigned firstneigh=0; if( sp_str.length()==0 ) firstneigh = (av->copyOutput(0))->getShape()[0];
388 3 : for(unsigned i=0; i<allspec.size(); ++i) {
389 2 : std::string onesstr="0"; if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[firstneigh])==allspec[i] ) onesstr = "1";
390 576 : for(unsigned j=firstneigh+1; j<natoms; ++j) {
391 574 : if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[j])==allspec[i] ) onesstr += ",1";
392 : else onesstr += ",0";
393 : }
394 4 : readInputLine( getShortcutLabel() + "_ones_" + allspec[i] + ": CONSTANT VALUES=" + onesstr );
395 : }
396 : }
397 :
398 15 : std::string envargstr,varstr, maxfuncstr, lambda; if( funcstr.size()>1 ) parse("LAMBDA",lambda);
399 : // And now do the funcstr bit
400 31 : for(unsigned j=0; j<funcstr.size(); ++j) {
401 21 : std::string jnum; Tools::convert( j+1, jnum );
402 21 : if(j==0) {
403 30 : varstr = "v" + jnum; maxfuncstr = "(1/" + lambda + ")*log(exp(" + lambda + "*v1)";
404 20 : envargstr = getShortcutLabel() + "_env" + jnum;
405 : } else {
406 33 : varstr += ",v" + jnum; maxfuncstr += "+exp(" + lambda + "*v" + jnum + ")";
407 22 : envargstr += "," + getShortcutLabel() + "_env" + jnum;
408 : }
409 : // And coordination numbers
410 21 : if( allspec.size()>1 ) {
411 : std::string argnames;
412 12 : for(unsigned i=0; i<allspec.size(); ++i) {
413 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] );
414 16 : readInputLine( getShortcutLabel() + "_" + allspec[i] + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_" + allspec[i] + "_matenv" + jnum + "," + getShortcutLabel() + "_ones_" + allspec[i] );
415 16 : if( i==0 ) argnames = getShortcutLabel() + "_" + allspec[i] + "_env" + jnum; else argnames += "," + getShortcutLabel() + "_" + allspec[i] + "_env" + jnum;
416 : }
417 4 : if( funcstr.size()==1) readInputLine( getShortcutLabel() + ": COMBINE PERIODIC=NO ARG=" + argnames );
418 8 : else readInputLine( getShortcutLabel() + "_env" + jnum + ": COMBINE PERIODIC=NO ARG=" + argnames );
419 : } else {
420 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] );
421 22 : if( funcstr.size()==1) readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
422 24 : else readInputLine( getShortcutLabel() + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
423 : }
424 : }
425 : // And get the maximum
426 15 : if( funcstr.size()>1 ) readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + envargstr + " PERIODIC=NO VAR=" + varstr + " FUNC=" + maxfuncstr + ")" );
427 : // Read in all the shortcut stuff
428 10 : std::map<std::string,std::string> keymap; multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
429 20 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
430 40 : }
431 :
432 14 : std::vector<std::pair<unsigned,Vector> > EnvironmentSimilarity::getReferenceEnvironment( const PDB& pdb, const std::vector<std::string>& anames, double& maxdist ) {
433 14 : unsigned natoms = pdb.getPositions().size(); std::vector<std::pair<unsigned,Vector> > env( natoms );
434 78 : for(unsigned i=0; i<natoms; ++i) {
435 : unsigned identity=0;
436 80 : for(unsigned j=1; j<anames.size(); ++j) {
437 16 : if( pdb.getAtomName(pdb.getAtomNumbers()[i])==anames[j] ) { identity=j; break; }
438 : }
439 64 : env[i] = std::pair<unsigned,Vector>( identity, pdb.getPositions()[i] );
440 64 : double dist = env[i].second.modulo();
441 64 : if( dist>maxdist ) maxdist = dist;
442 : }
443 14 : return env;
444 : }
445 :
446 : }
447 : }
|
