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 : The next example compares the environments of the 96 selected atoms with a user specified reference
124 : environment. The reference environment is contained in the env1.pdb file. Once the kernel is computed
125 : the average and the number of atoms with a kernel larger than 0.5 are computed.
126 :
127 : ```plumed
128 : #SETTINGS INPUTFILES=regtest/envsim/rt-env-sim-atom-names-match/env1.pdb
129 : es: ENVIRONMENTSIMILARITY ...
130 : SPECIES=1-288:3
131 : SIGMA=0.05
132 : CRYSTAL_STRUCTURE=CUSTOM
133 : REFERENCE=regtest/envsim/rt-env-sim-atom-names-match/env1.pdb
134 : MEAN
135 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
136 : ...
137 :
138 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
139 : ```
140 :
141 : The next example is similar to the one above but in this case 4 reference environments are specified.
142 : Each reference environment is given in a separate pdb file.
143 :
144 : ```plumed
145 : #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
146 : es: ENVIRONMENTSIMILARITY ...
147 : SPECIES=1-288:3
148 : SIGMA=0.05
149 : CRYSTAL_STRUCTURE=CUSTOM
150 : REFERENCE_1=regtest/envsim/rt-env-sim-atom-names-match/env1.pdb
151 : REFERENCE_2=regtest/envsim/rt-env-sim-atom-names-match/env2.pdb
152 : REFERENCE_3=regtest/envsim/rt-env-sim-atom-names-match/env3.pdb
153 : REFERENCE_4=regtest/envsim/rt-env-sim-atom-names-match/env4.pdb
154 : MEAN
155 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
156 : ...
157 :
158 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
159 : ```
160 :
161 : The following examples illustrates the use of pdb files to provide information about different chemical species:
162 :
163 :
164 : ```plumed
165 : #SETTINGS INPUTFILES=regtest/envsim/rt-env-sim-custom-1env/env1.pdb,regtest/envsim/rt-env-sim-atom-names-match/IceIh-atom-names.pdb
166 : es: ENVIRONMENTSIMILARITY ...
167 : SPECIES=1-384
168 : SIGMA=0.05
169 : CRYSTAL_STRUCTURE=CUSTOM
170 : REFERENCE=regtest/envsim/rt-env-sim-custom-1env/env1.pdb
171 : MEAN
172 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
173 : ATOM_NAMES_FILE=regtest/envsim/rt-env-sim-atom-names-match/IceIh-atom-names.pdb
174 : ...
175 : ```
176 :
177 : In this case, all atoms are used as centers, but only neighbors of type O are taken into account.
178 :
179 : */
180 : //+ENDPLUMEDOC
181 :
182 : class EnvironmentSimilarity : public ActionShortcut {
183 : private:
184 : std::vector<std::pair<unsigned,Vector> > getReferenceEnvironment( const PDB& pdb, const std::vector<std::string>& anames, double& maxdist );
185 : public:
186 : static void registerKeywords( Keywords& keys );
187 : explicit EnvironmentSimilarity(const ActionOptions&);
188 : };
189 :
190 : PLUMED_REGISTER_ACTION(EnvironmentSimilarity,"ENVIRONMENTSIMILARITY")
191 :
192 28 : void EnvironmentSimilarity::registerKeywords( Keywords& keys ) {
193 28 : ActionShortcut::registerKeywords( keys );
194 28 : keys.add("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
195 : "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
196 : "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
197 : "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
198 : "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
199 : "coordination number more than four for example");
200 28 : 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 "
201 : "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
202 : "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
203 : "using the label of another multicolvar");
204 28 : 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 "
205 : "the documentation for that keyword");
206 28 : keys.add("compulsory","CRYSTAL_STRUCTURE","FCC","Targeted crystal structure. Options are: "
207 : "SC: simple cubic, "
208 : "BCC: body center cubic, "
209 : "FCC: face centered cubic, "
210 : "HCP: hexagonal closed pack, "
211 : "DIAMOND: cubic diamond, "
212 : "CUSTOM: user defined "
213 : " ");
214 28 : keys.add("compulsory","LATTICE_CONSTANTS","Lattice constants. Two comma separated values for HCP, "
215 : "one value for all other CRYSTAL_STRUCTURES.");
216 28 : keys.add("compulsory","SIGMA","0.1","the width to use for the gaussian kernels");
217 28 : keys.add("compulsory","LCUTOFF","0.0001","any atoms separated by less than this tolerance should be ignored");
218 28 : keys.add("optional","REFERENCE","PDB files with relative distances from central atom. Use this keyword if you are targeting a single reference environment.");
219 28 : 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.");
220 28 : keys.add("compulsory","LAMBDA","100","Lambda parameter. This is only used if you have more than one reference environment");
221 28 : keys.add("compulsory","CUTOFF","3","how many multiples of sigma would you like to consider beyond the maximum distance in the environment");
222 28 : 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.");
223 56 : keys.setValueDescription("vector","the environmental similar parameter for each of the input atoms");
224 28 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
225 28 : keys.needsAction("GROUP");
226 28 : keys.needsAction("DISTANCE_MATRIX");
227 28 : keys.needsAction("ONES");
228 28 : keys.needsAction("CONSTANT");
229 28 : keys.needsAction("CUSTOM");
230 28 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
231 28 : keys.needsAction("COMBINE");
232 28 : keys.addDOI("10.1063/1.5102104");
233 28 : }
234 :
235 10 : EnvironmentSimilarity::EnvironmentSimilarity(const ActionOptions&ao):
236 : Action(ao),
237 10 : ActionShortcut(ao) {
238 : std::string atomNamesFile;
239 10 : parse("ATOM_NAMES_FILE",atomNamesFile);
240 10 : PDB atomnamepdb;
241 10 : if( !atomNamesFile.empty() && !atomnamepdb.read(atomNamesFile,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
242 0 : error("missing input file " + atomNamesFile);
243 : }
244 :
245 10 : double maxdist=0;
246 10 : std::vector<std::string> allspec(1);
247 : std::string crystal_structure;
248 20 : parse("CRYSTAL_STRUCTURE", crystal_structure);
249 : std::vector<std::vector<std::pair<unsigned,Vector> > > environments;
250 10 : if( crystal_structure=="CUSTOM" ) {
251 5 : if( !atomNamesFile.empty() ) {
252 1 : allspec[0]=atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[0]);
253 1 : unsigned natoms=atomnamepdb.getPositions().size();
254 385 : for(unsigned i=0; i<natoms; ++i) {
255 : bool found=false;
256 576 : for(unsigned j=0; j<allspec.size(); ++j) {
257 575 : if( allspec[j]==atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i] ) ) {
258 : found=true;
259 : break;
260 : }
261 : }
262 384 : if( !found ) {
263 2 : allspec.push_back( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[i]) );
264 : }
265 : }
266 : }
267 : std::string reffile;
268 10 : parse("REFERENCE",reffile);
269 5 : if( reffile.length()>0 ) {
270 2 : PDB pdb;
271 2 : pdb.read(reffile,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
272 2 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
273 2 : log.printf(" reading %d reference vectors from %s \n", environments[0].size(), reffile.c_str() );
274 2 : } else {
275 12 : for(unsigned int i=1;; i++) {
276 15 : PDB pdb;
277 30 : if( !parseNumbered("REFERENCE_",i,reffile) ) {
278 : break;
279 : }
280 12 : if( !pdb.read(reffile,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
281 0 : error("missing input file " + reffile );
282 : }
283 12 : environments.push_back( getReferenceEnvironment( pdb, allspec, maxdist ) );
284 12 : log.printf(" Reference environment %d : reading %d reference vectors from %s \n", i, environments[i-1].size(), reffile.c_str() );
285 15 : }
286 : }
287 : } else {
288 : std::vector<double> lattice_constants;
289 10 : parseVector("LATTICE_CONSTANTS", lattice_constants);
290 5 : if (crystal_structure == "FCC") {
291 1 : if (lattice_constants.size() != 1) {
292 0 : error("Number of LATTICE_CONSTANTS arguments must be one for FCC");
293 : }
294 1 : environments.resize(1);
295 1 : environments[0].resize(12);
296 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.0)*lattice_constants[0] );
297 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.0)*lattice_constants[0] );
298 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.0)*lattice_constants[0] );
299 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.0)*lattice_constants[0] );
300 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,+0.5)*lattice_constants[0] );
301 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,-0.5)*lattice_constants[0] );
302 1 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.0,+0.5)*lattice_constants[0] );
303 1 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.0,-0.5)*lattice_constants[0] );
304 1 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,+0.5)*lattice_constants[0] );
305 1 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,-0.5)*lattice_constants[0] );
306 1 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-0.5,+0.5)*lattice_constants[0] );
307 1 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.5,-0.5)*lattice_constants[0] );
308 1 : maxdist = std::sqrt(2)*lattice_constants[0]/2.;
309 4 : } else if (crystal_structure == "SC") {
310 0 : if (lattice_constants.size() != 1) {
311 0 : error("Number of LATTICE_CONSTANTS arguments must be one for SC");
312 : }
313 0 : environments.resize(1);
314 0 : environments[0].resize(6);
315 0 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
316 0 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
317 0 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
318 0 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
319 0 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
320 0 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
321 0 : maxdist = lattice_constants[0];
322 4 : } else if( crystal_structure == "BCC") {
323 2 : if (lattice_constants.size() != 1) {
324 0 : error("Number of LATTICE_CONSTANTS arguments must be one for BCC");
325 : }
326 2 : environments.resize(1);
327 2 : environments[0].resize(14);
328 2 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,+0.5)*lattice_constants[0] );
329 2 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,-0.5)*lattice_constants[0] );
330 2 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,+0.5)*lattice_constants[0] );
331 2 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,+0.5)*lattice_constants[0] );
332 2 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+0.5,-0.5)*lattice_constants[0] );
333 2 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-0.5,+0.5)*lattice_constants[0] );
334 2 : environments[0][6] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-0.5,-0.5)*lattice_constants[0] );
335 2 : environments[0][7] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+0.5,-0.5)*lattice_constants[0] );
336 2 : environments[0][8] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0)*lattice_constants[0] );
337 2 : environments[0][9] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+1.0,+0.0)*lattice_constants[0] );
338 2 : environments[0][10] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,+1.0)*lattice_constants[0] );
339 2 : environments[0][11] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0)*lattice_constants[0] );
340 2 : environments[0][12] = std::pair<unsigned,Vector>( 0, Vector(+0.0,-1.0,+0.0)*lattice_constants[0] );
341 2 : environments[0][13] = std::pair<unsigned,Vector>( 0, Vector(+0.0,+0.0,-1.0)*lattice_constants[0] );
342 2 : maxdist = lattice_constants[0];
343 2 : } else if (crystal_structure == "HCP") {
344 1 : if (lattice_constants.size() != 2) {
345 0 : error("Number of LATTICE_CONSTANTS arguments must be two for HCP");
346 : }
347 1 : environments.resize(2);
348 1 : environments[0].resize(12);
349 1 : environments[1].resize(12);
350 : double sqrt3=std::sqrt(3);
351 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
352 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
353 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
354 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
355 1 : environments[0][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
356 1 : environments[0][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
357 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] );
358 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] );
359 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] );
360 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] );
361 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] );
362 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] );
363 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
364 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0] );
365 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
366 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0] );
367 1 : environments[1][4] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+0.0,+0.0) *lattice_constants[0] );
368 1 : environments[1][5] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+0.0,+0.0) *lattice_constants[0] );
369 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] );
370 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] );
371 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] );
372 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] );
373 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] );
374 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] );
375 1 : maxdist = lattice_constants[0];
376 1 : } else if (crystal_structure == "DIAMOND") {
377 1 : if (lattice_constants.size() != 1) {
378 0 : error("Number of LATTICE_CONSTANTS arguments must be one for DIAMOND");
379 : }
380 1 : environments.resize(2);
381 1 : environments[0].resize(4);
382 1 : environments[1].resize(4);
383 1 : environments[0][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
384 1 : environments[0][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
385 1 : environments[0][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
386 1 : environments[0][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
387 1 : environments[1][0] = std::pair<unsigned,Vector>( 0, Vector(+1.0,-1.0,+1.0)*lattice_constants[0]/4.0 );
388 1 : environments[1][1] = std::pair<unsigned,Vector>( 0, Vector(-1.0,+1.0,+1.0)*lattice_constants[0]/4.0 );
389 1 : environments[1][2] = std::pair<unsigned,Vector>( 0, Vector(+1.0,+1.0,-1.0)*lattice_constants[0]/4.0 );
390 1 : environments[1][3] = std::pair<unsigned,Vector>( 0, Vector(-1.0,-1.0,-1.0)*lattice_constants[0]/4.0 );
391 1 : maxdist = std::sqrt(3)*lattice_constants[0]/4.0;
392 : } else {
393 0 : error( crystal_structure + " is not a valid input for keyword CRYSTAL_STRUCTURE");
394 : }
395 : }
396 10 : std::string matlab = getShortcutLabel() + "_cmat";
397 : double cutoff, sig;
398 10 : parse("SIGMA",sig);
399 20 : parse("CUTOFF",cutoff);
400 : std::string lcutoff;
401 20 : parse("LCUTOFF",lcutoff);
402 : std::string sig2;
403 10 : Tools::convert( sig*sig, sig2 );
404 10 : std::vector<std::vector<std::string> > funcstr(environments.size());
405 : std::string str_cutoff;
406 10 : Tools::convert( maxdist + cutoff*sig, str_cutoff );
407 : std::string str_natoms, xpos, ypos, zpos;
408 10 : Tools::convert( environments[0].size(), str_natoms );
409 31 : for(unsigned j=0; j<environments.size(); ++j) {
410 21 : funcstr[j].resize( allspec.size() );
411 46 : for(unsigned k=0; k<allspec.size(); ++k) {
412 177 : for(unsigned i=0; i<environments[j].size(); ++i) {
413 152 : if( environments[j][i].first!=k ) {
414 16 : continue ;
415 : }
416 136 : Tools::convert( environments[j][i].second[0], xpos );
417 136 : Tools::convert( environments[j][i].second[1], ypos );
418 136 : Tools::convert( environments[j][i].second[2], zpos );
419 136 : if( i==0 ) {
420 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 + "))";
421 : } else {
422 230 : funcstr[j][k] += "+exp(-((x-" + xpos + ")^2+(y-" + ypos + ")^2+(z-" + zpos + ")^2)/(4*" + sig2 + "))";
423 : }
424 : }
425 25 : if( funcstr[j][k].length()>0 ) {
426 : funcstr[j][k] += ")";
427 : } else {
428 : funcstr[j][k] ="FUNC=0";
429 : }
430 : }
431 : }
432 :
433 : // Create the constact matrix
434 : std::string sp_str, specA, specB;
435 10 : parse("SPECIES",sp_str);
436 10 : parse("SPECIESA",specA);
437 20 : parse("SPECIESB",specB);
438 10 : if( sp_str.length()>0 ) {
439 18 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUP=" + sp_str + " CUTOFF=" + str_cutoff );
440 18 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + sp_str );
441 : } else {
442 1 : if( specA.length()==0 ) {
443 0 : error("no atoms were specified use SPECIES or SPECIESA+SPECIESB");
444 : }
445 1 : if( specB.length()==0 ) {
446 0 : error("no atoms were specified for SPECIESB");
447 : }
448 2 : readInputLine( matlab + ": DISTANCE_MATRIX COMPONENTS GROUPA=" + specA + " GROUPB=" + specB + " CUTOFF=" + str_cutoff );
449 2 : readInputLine( getShortcutLabel() + "_grp: GROUP ATOMS=" + specA );
450 : }
451 :
452 : // Make a vector containing all ones
453 10 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( matlab );
454 10 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 );
455 : std::string size;
456 10 : Tools::convert( (av->copyOutput(0))->getShape()[1], size );
457 10 : if( allspec.size()==1 ) {
458 18 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size );
459 : } else {
460 1 : unsigned natoms=atomnamepdb.getPositions().size();
461 : unsigned firstneigh=0;
462 1 : if( sp_str.length()==0 ) {
463 1 : firstneigh = (av->copyOutput(0))->getShape()[0];
464 : }
465 3 : for(unsigned i=0; i<allspec.size(); ++i) {
466 2 : std::string onesstr="0";
467 2 : if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[firstneigh])==allspec[i] ) {
468 : onesstr = "1";
469 : }
470 576 : for(unsigned j=firstneigh+1; j<natoms; ++j) {
471 574 : if( atomnamepdb.getAtomName(atomnamepdb.getAtomNumbers()[j])==allspec[i] ) {
472 : onesstr += ",1";
473 : } else {
474 : onesstr += ",0";
475 : }
476 : }
477 4 : readInputLine( getShortcutLabel() + "_ones_" + allspec[i] + ": CONSTANT VALUES=" + onesstr );
478 : }
479 : }
480 :
481 : std::string envargstr,varstr, maxfuncstr, lambda;
482 10 : if( funcstr.size()>1 ) {
483 10 : parse("LAMBDA",lambda);
484 : }
485 : // And now do the funcstr bit
486 31 : for(unsigned j=0; j<funcstr.size(); ++j) {
487 : std::string jnum;
488 21 : Tools::convert( j+1, jnum );
489 21 : if(j==0) {
490 10 : varstr = "v" + jnum;
491 20 : maxfuncstr = "(1/" + lambda + ")*log(exp(" + lambda + "*v1)";
492 20 : envargstr = getShortcutLabel() + "_env" + jnum;
493 : } else {
494 11 : varstr += ",v" + jnum;
495 22 : maxfuncstr += "+exp(" + lambda + "*v" + jnum + ")";
496 22 : envargstr += "," + getShortcutLabel() + "_env" + jnum;
497 : }
498 : // And coordination numbers
499 21 : if( allspec.size()>1 ) {
500 : std::string argnames;
501 12 : for(unsigned i=0; i<allspec.size(); ++i) {
502 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] );
503 16 : readInputLine( getShortcutLabel() + "_" + allspec[i] + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_" + allspec[i] + "_matenv" + jnum + "," + getShortcutLabel() + "_ones_" + allspec[i] );
504 8 : if( i==0 ) {
505 8 : argnames = getShortcutLabel() + "_" + allspec[i] + "_env" + jnum;
506 : } else {
507 8 : argnames += "," + getShortcutLabel() + "_" + allspec[i] + "_env" + jnum;
508 : }
509 : }
510 4 : if( funcstr.size()==1) {
511 0 : readInputLine( getShortcutLabel() + ": COMBINE PERIODIC=NO ARG=" + argnames );
512 : } else {
513 8 : readInputLine( getShortcutLabel() + "_env" + jnum + ": COMBINE PERIODIC=NO ARG=" + argnames );
514 : }
515 : } else {
516 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] );
517 17 : if( funcstr.size()==1) {
518 10 : readInputLine( getShortcutLabel() + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
519 : } else {
520 24 : readInputLine( getShortcutLabel() + "_env" + jnum + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_matenv" + jnum + "," + getShortcutLabel() + "_ones");
521 : }
522 : }
523 : }
524 : // And get the maximum
525 10 : if( funcstr.size()>1 ) {
526 10 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + envargstr + " PERIODIC=NO VAR=" + varstr + " FUNC=" + maxfuncstr + ")" );
527 : }
528 : // Read in all the shortcut stuff
529 : std::map<std::string,std::string> keymap;
530 10 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
531 20 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", keymap, this );
532 40 : }
533 :
534 14 : std::vector<std::pair<unsigned,Vector> > EnvironmentSimilarity::getReferenceEnvironment( const PDB& pdb, const std::vector<std::string>& anames, double& maxdist ) {
535 14 : unsigned natoms = pdb.getPositions().size();
536 14 : std::vector<std::pair<unsigned,Vector> > env( natoms );
537 78 : for(unsigned i=0; i<natoms; ++i) {
538 : unsigned identity=0;
539 80 : for(unsigned j=1; j<anames.size(); ++j) {
540 16 : if( pdb.getAtomName(pdb.getAtomNumbers()[i])==anames[j] ) {
541 : identity=j;
542 : break;
543 : }
544 : }
545 64 : env[i] = std::pair<unsigned,Vector>( identity, pdb.getPositions()[i] );
546 64 : double dist = env[i].second.modulo();
547 64 : if( dist>maxdist ) {
548 13 : maxdist = dist;
549 : }
550 : }
551 14 : return env;
552 : }
553 :
554 : }
555 : }
|
