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