Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : /* ----------------------------------------------------------------------
24 : Contributing author: Pablo Piaggi (Princeton University)
25 : ------------------------------------------------------------------------- */
26 :
27 : #include "multicolvar/MultiColvarBase.h"
28 : #include "multicolvar/AtomValuePack.h"
29 : #include "core/ActionRegister.h"
30 : #include "core/PlumedMain.h"
31 : #include "tools/PDB.h"
32 : #include <limits>
33 :
34 : using namespace std;
35 :
36 : namespace PLMD {
37 : namespace crystallization {
38 :
39 :
40 : //+PLUMEDOC MCOLVAR ENVIRONMENTSIMILARITY
41 : /*
42 : Measure how similar the environment around atoms is to that found in some reference crystal structure.
43 :
44 : This CV was introduced in this article \cite Piaggi-JCP-2019.
45 : The starting point for the definition of the CV is the local atomic density around an atom.
46 : We consider an environment \f$\chi\f$ around this atom and we define the density by
47 : \f[
48 : \rho_{\chi}(\mathbf{r})=\sum\limits_{i\in\chi} \exp\left(- \frac{|\mathbf{r}_i-\mathbf{r}|^2} {2\sigma^2} \right),
49 : \f]
50 : 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
51 : coordinates of the neighbors relative to the central atom.
52 : 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$
53 : that describe, for instance, the nearest neighbors in a given lattice.
54 : \f$\sigma\f$ is set using the SIGMA keyword and \f$\chi_0\f$ is chosen with the CRYSTAL_STRUCTURE keyword.
55 : If only the SPECIES keyword is given then the atoms defined there will be the central and neighboring atoms.
56 : If instead the SPECIESA and SPECIESB keywords are given then SPECIESA determines the central atoms and SPECIESB the neighbors.
57 :
58 : The environments \f$\chi\f$ and \f$\chi_0\f$ are compared using the kernel,
59 : \f[
60 : k_{\chi_0}(\chi)= \int d\mathbf{r} \rho_{\chi}(\mathbf{r}) \rho_{\chi_0}(\mathbf{r}) .
61 : \f]
62 : Combining the two equations above and performing the integration analytically we obtain,
63 : \f[
64 : 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).
65 : \f]
66 : The kernel is finally normalized,
67 : \f[
68 : \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),
69 : \f]
70 : such that \f$\tilde{k}_{\chi_0}(\chi_0) = 1\f$.
71 : The above kernel is computed for each atom in the SPECIES or SPECIESA keywords.
72 : This quantity is a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
73 : 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
74 : so on.
75 :
76 : The kernel can be generalized to crystal structures described as a lattice with a basis of more than one atom.
77 : In this case there is more than one type of environment.
78 : 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:
79 : \f[
80 : \tilde{k}_X(\chi)= \frac{1}{\lambda} \log \left ( \sum\limits_{l=1}^{M}\exp \left (\lambda \: \tilde{k}_{\chi_l}(\chi) \right ) \right ).
81 : \f]
82 : 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$.
83 : This approach can be used, for instance, to target the hexagonal closed packed (HCP keyword) or the diamond structure (DIAMOND keyword).
84 :
85 : The CRYSTAL_STRUCTURE keyword can take the values SC (simple cubic), BCC (body centered cubic), FCC (face centered cubic),
86 : HCP (hexagonal closed pack), DIAMOND (cubic diamond), and CUSTOM (user defined).
87 : All options follow the same conventions as in the [lattice command](https://lammps.sandia.gov/doc/lattice.html) of [LAMMPS](https://lammps.sandia.gov/).
88 : If a CRYSTAL_STRUCTURE other than CUSTOM is used, then the lattice constants have to be specified using the keyword LATTICE_CONSTANTS.
89 : 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).
90 :
91 : If the CUSTOM option is used then the reference environments have to be specified by the user.
92 : The reference environments are specified in pdb files containing the distance vectors from the central atom to the neighbors.
93 : Make sure your PDB file is correctly formatted as explained \ref pdbreader "in this page"
94 : If only one reference environment is specified then the filename should be given as argument of the keyword REFERENCE.
95 : If instead several reference environments are given, then they have to be provided in separate pdb files and given as arguments of the
96 : keywords REFERENCE_1, REFERENCE_2, etc.
97 : 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.
98 :
99 : 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.
100 : This information is provided in pdb files using the atom name field.
101 : The comparison between environments is performed taking into account whether the atom names match.
102 :
103 : \par Examples
104 :
105 : The following input calculates the ENVIRONMENTSIMILARITY kernel for 250 atoms in the system
106 : using the BCC atomic environment as target, and then calculates and prints the average value
107 : for this quantity.
108 :
109 : \plumedfile
110 : ENVIRONMENTSIMILARITY SPECIES=1-250 SIGMA=0.05 LATTICE_CONSTANTS=0.423 CRYSTAL_STRUCTURE=BCC MEAN LABEL=es
111 :
112 : PRINT ARG=es.mean FILE=COLVAR
113 : \endplumedfile
114 :
115 : The next example compares the environments of the 96 selected atoms with a user specified reference
116 : environment. The reference environment is contained in the env1.pdb file. Once the kernel is computed
117 : the average and the number of atoms with a kernel larger than 0.5 are computed.
118 :
119 : \plumedfile
120 : ENVIRONMENTSIMILARITY ...
121 : SPECIES=1-288:3
122 : SIGMA=0.05
123 : CRYSTAL_STRUCTURE=CUSTOM
124 : REFERENCE=env1.pdb
125 : LABEL=es
126 : MEAN
127 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
128 : ... ENVIRONMENTSIMILARITY
129 :
130 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
131 : \endplumedfile
132 :
133 : The next example is similar to the one above but in this case 4 reference environments are specified.
134 : Each reference environment is given in a separate pdb file.
135 :
136 : \plumedfile
137 : ENVIRONMENTSIMILARITY ...
138 : SPECIES=1-288:3
139 : SIGMA=0.05
140 : CRYSTAL_STRUCTURE=CUSTOM
141 : REFERENCE_1=env1.pdb
142 : REFERENCE_2=env2.pdb
143 : REFERENCE_3=env3.pdb
144 : REFERENCE_4=env4.pdb
145 : LABEL=es
146 : MEAN
147 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
148 : ... ENVIRONMENTSIMILARITY
149 :
150 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
151 : \endplumedfile
152 :
153 : The following examples illustrates the use of pdb files to provide information about different chemical species:
154 : \plumedfile
155 : ENVIRONMENTSIMILARITY ...
156 : SPECIES=1-6
157 : SIGMA=0.05
158 : CRYSTAL_STRUCTURE=CUSTOM
159 : REFERENCE=env.pdb
160 : LABEL=es
161 : MEAN
162 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
163 : ATOM_NAMES_FILE=atom-names.pdb
164 : ... ENVIRONMENTSIMILARITY
165 : \endplumedfile
166 : Here the file env.pdb is:
167 : \verbatim
168 : ATOM 1 O MOL 1 -2.239 -1.296 -0.917 1.00 0.00 O
169 : ATOM 2 O MOL 1 0.000 0.000 2.751 1.00 0.00 O
170 : \endverbatim
171 : where atoms are of type O, and the atom-names.pdb file is:
172 : \verbatim
173 : ATOM 1 O X 1 0.000 2.593 4.126 0.00 0.00 O
174 : ATOM 2 H X 1 0.000 3.509 3.847 0.00 0.00 H
175 : ATOM 3 H X 1 0.000 2.635 5.083 0.00 0.00 H
176 : ATOM 4 O X 1 0.000 2.593 11.462 0.00 0.00 O
177 : ATOM 5 H X 1 0.000 3.509 11.183 0.00 0.00 H
178 : ATOM 6 H X 1 0.000 2.635 12.419 0.00 0.00 H
179 : \endverbatim
180 : where atoms are of type O and H.
181 : In this case, all atoms are used as centers, but only neighbors of type O are taken into account.
182 :
183 : */
184 : //+ENDPLUMEDOC
185 :
186 :
187 : class EnvironmentSimilarity : public multicolvar::MultiColvarBase {
188 : private:
189 : // All global variables end with underscore
190 : // square of cutoff, square of broadening parameter
191 : double rcut2_, sigmaSqr_;
192 : // lambda parameter for softmax function
193 : double lambda_;
194 : // Array of Vectors to store the reference environments, i.e. the templates
195 : std::vector<std::vector<Vector>> environments_;
196 : // Array of strings to store the atom names of reference environments
197 : std::vector<std::vector<std::string>> environmentsAtomNames_;
198 : // Use of atom names reference file
199 : std::vector<std::string> atomNames_;
200 : public:
201 : static void registerKeywords( Keywords& keys );
202 : explicit EnvironmentSimilarity(const ActionOptions&);
203 : // active methods:
204 : virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
205 : // Returns the number of coordinates of the field
206 17 : bool isPeriodic() {
207 17 : return false;
208 : }
209 : // Calculates maximum distance in an environment
210 : double maxDistance(std::vector<Vector> environment);
211 : // Parse everything connected to the definition of the reference environments
212 : // First argument is the array of Vectors that stores the reference environments
213 : // Second argument is the array of strings that stores the atom names of reference environments
214 : // Third argument is the maximum distance in the ref environments and sets the
215 : // cutoff for the cell lists
216 : void parseReferenceEnvironments( std::vector<std::vector<Vector>>& environments, std::vector<std::vector<std::string>>& environmentsAtomNames, double& max_dist);
217 : };
218 :
219 13805 : PLUMED_REGISTER_ACTION(EnvironmentSimilarity,"ENVIRONMENTSIMILARITY")
220 :
221 14 : void EnvironmentSimilarity::registerKeywords( Keywords& keys ) {
222 14 : MultiColvarBase::registerKeywords( keys );
223 14 : keys.use("SPECIES");
224 14 : keys.use("SPECIESA");
225 14 : keys.use("SPECIESB");
226 28 : keys.add("compulsory","SIGMA","0.1","Broadening parameter");
227 28 : keys.add("compulsory","CRYSTAL_STRUCTURE","FCC","Targeted crystal structure. Options are: "
228 : "SC: simple cubic, "
229 : "BCC: body center cubic, "
230 : "FCC: face centered cubic, "
231 : "HCP: hexagonal closed pack, "
232 : "DIAMOND: cubic diamond, "
233 : "CUSTOM: user defined "
234 : " ");
235 28 : keys.add("optional","LATTICE_CONSTANTS","Lattice constants. Two comma separated values for HCP, "
236 : "one value for all other CRYSTAL_STRUCTURES.");
237 28 : keys.add("compulsory","LAMBDA","100","Lambda parameter");
238 28 : keys.add("optional","REFERENCE","PDB file with relative distances from central atom."
239 : " Use this keyword if you are targeting a single reference environment.");
240 28 : keys.add("numbered","REFERENCE_","PDB files with relative distances from central atom."
241 : " Each file corresponds to one template."
242 : " Use these keywords if you are targeting more than one reference environment.");
243 28 : keys.add("optional","ATOM_NAMES_FILE","PDB file with atom names for all atoms in SPECIES."
244 : " Atoms in reference environments will be compared only if atom names match.");
245 : // Use actionWithDistributionKeywords
246 14 : keys.use("MEAN");
247 14 : keys.use("MORE_THAN");
248 14 : keys.use("LESS_THAN");
249 14 : keys.use("MAX");
250 14 : keys.use("MIN");
251 14 : keys.use("BETWEEN");
252 14 : keys.use("HISTOGRAM");
253 14 : keys.use("MOMENTS");
254 14 : keys.use("ALT_MIN");
255 14 : keys.use("LOWEST");
256 14 : keys.use("HIGHEST");
257 14 : }
258 :
259 10 : EnvironmentSimilarity::EnvironmentSimilarity(const ActionOptions&ao):
260 : Action(ao),
261 10 : MultiColvarBase(ao) {
262 10 : log.printf(" Please read and cite ");
263 20 : log << plumed.cite("Piaggi and Parrinello, J. Chem. Phys. 150 (24), 244119 (2019)");
264 10 : log.printf("\n");
265 :
266 : // Parse everything connected to the definition of the reference environments
267 : double max_dist_ref_vector;
268 10 : parseReferenceEnvironments(environments_, environmentsAtomNames_, max_dist_ref_vector);
269 :
270 : double sigma;
271 10 : parse("SIGMA", sigma);
272 10 : log.printf(" representing local density as a sum of Gaussians with standard deviation %f\n",sigma);
273 10 : sigmaSqr_=sigma*sigma;
274 :
275 10 : lambda_=100;
276 20 : parse("LAMBDA", lambda_);
277 10 : if (environments_.size()>1) {
278 6 : log.printf(" using a soft max function with lambda %f\n",lambda_);
279 : }
280 :
281 : // Set the link cell cutoff
282 10 : double rcut = max_dist_ref_vector + 3*sigma;
283 10 : setLinkCellCutoff( rcut );
284 10 : rcut2_ = rcut * rcut;
285 :
286 : // And setup the ActionWithVessel
287 : std::vector<AtomNumber> all_atoms;
288 10 : setupMultiColvarBase( all_atoms );
289 10 : checkRead();
290 :
291 10 : plumed_massert(atomNames_.empty() || atomNames_.size()==getNumberOfAtoms(),"mismatch between atoms in SPECIES and ATOM_NAMES_FILE");
292 10 : }
293 :
294 30144 : double EnvironmentSimilarity::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
295 30144 : if (environments_.size()==1 && atomNames_.empty() ) {
296 : // One reference environment case - no atom names
297 155464 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
298 : const Vector& distance=myatoms.getPosition(i);
299 : double d2;
300 232452 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
301 77840 : (d2+=distance[1]*distance[1])<rcut2_ &&
302 194068 : (d2+=distance[2]*distance[2])<rcut2_ &&
303 : d2>epsilon ) {
304 : // Iterate over atoms in the reference environment
305 236856 : for(unsigned k=0; k<environments_[0].size(); ++k) {
306 220400 : Vector distanceFromRef=distance-environments_[0][k];
307 220400 : double value = std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[0].size() ;
308 220400 : accumulateSymmetryFunction( 1, i, value, -(value/(2*sigmaSqr_))*distanceFromRef, (value/(2*sigmaSqr_))*Tensor(distance,distanceFromRef), myatoms );
309 : }
310 : }
311 : }
312 : return myatoms.getValue(1);
313 29292 : } else if (atomNames_.empty()) {
314 : // More than one reference environment case - no atom names
315 29196 : std::vector<double> values(environments_.size()); //value for each template
316 : // First time calculate sums
317 2808756 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
318 : const Vector& distance=myatoms.getPosition(i);
319 : double d2;
320 4534472 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
321 1754912 : (d2+=distance[1]*distance[1])<rcut2_ &&
322 3469608 : (d2+=distance[2]*distance[2])<rcut2_ &&
323 : d2>epsilon ) {
324 : // Iterate over templates
325 1216388 : for(unsigned j=0; j<environments_.size(); ++j) {
326 : // Iterate over atoms in the template
327 4893780 : for(unsigned k=0; k<environments_[j].size(); ++k) {
328 3921936 : Vector distanceFromRef=distance-environments_[j][k];
329 3921936 : values[j] += std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
330 : }
331 : }
332 : }
333 : }
334 : double sum=0;
335 145188 : for(unsigned j=0; j<environments_.size(); ++j) {
336 115992 : values[j] = std::exp(lambda_*values[j]);
337 115992 : sum += values[j];
338 : }
339 : // Second time find derivatives
340 2808756 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
341 : const Vector& distance=myatoms.getPosition(i);
342 : double d2;
343 4534472 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
344 1754912 : (d2+=distance[1]*distance[1])<rcut2_ &&
345 3469608 : (d2+=distance[2]*distance[2])<rcut2_ &&
346 : d2>epsilon ) {
347 : // Iterate over reference environment
348 1216388 : for(unsigned j=0; j<environments_.size(); ++j) {
349 : // Iterate over atoms in the reference environment
350 4893780 : for(unsigned k=0; k<environments_[j].size(); ++k) {
351 3921936 : Vector distanceFromRef=distance-environments_[j][k];
352 3921936 : double value = std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
353 3921936 : accumulateSymmetryFunction( 1, i, value, -(values[j]/sum)*(value/(2*sigmaSqr_))*distanceFromRef, (values[j]/sum)*(value/(2*sigmaSqr_))*Tensor(distance,distanceFromRef), myatoms );
354 : }
355 : }
356 : }
357 : }
358 29196 : return std::log(sum)/lambda_;
359 : } else {
360 : // Reference environments with atom names
361 96 : std::vector<double> values(environments_.size()); //value for each template
362 : // First time calculate sums
363 27648 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
364 : const Vector& distance=myatoms.getPosition(i);
365 : double d2;
366 42816 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
367 15264 : (d2+=distance[1]*distance[1])<rcut2_ &&
368 33216 : (d2+=distance[2]*distance[2])<rcut2_ &&
369 : d2>epsilon ) {
370 : // Iterate over templates
371 9600 : for(unsigned j=0; j<environments_.size(); ++j) {
372 : // Iterate over atoms in the template
373 38400 : for(unsigned k=0; k<environments_[j].size(); ++k) {
374 30720 : if (atomNames_[myatoms.getIndex(i)]==environmentsAtomNames_[j][k]) {
375 6144 : Vector distanceFromRef=distance-environments_[j][k];
376 6144 : values[j] += std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
377 : }
378 : }
379 : }
380 : }
381 : }
382 : double sum=0;
383 480 : for(unsigned j=0; j<environments_.size(); ++j) {
384 384 : values[j] = std::exp(lambda_*values[j]);
385 384 : sum += values[j];
386 : }
387 : // Second time find derivatives
388 27648 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
389 : const Vector& distance=myatoms.getPosition(i);
390 : double d2;
391 42816 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
392 15264 : (d2+=distance[1]*distance[1])<rcut2_ &&
393 33216 : (d2+=distance[2]*distance[2])<rcut2_ &&
394 : d2>epsilon ) {
395 : // Iterate over reference environment
396 9600 : for(unsigned j=0; j<environments_.size(); ++j) {
397 : // Iterate over atoms in the reference environment
398 38400 : for(unsigned k=0; k<environments_[j].size(); ++k) {
399 30720 : if (atomNames_[myatoms.getIndex(i)]==environmentsAtomNames_[j][k]) {
400 6144 : Vector distanceFromRef=distance-environments_[j][k];
401 6144 : double value = std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
402 6144 : accumulateSymmetryFunction( 1, i, value, -(values[j]/sum)*(value/(2*sigmaSqr_))*distanceFromRef, (values[j]/sum)*(value/(2*sigmaSqr_))*Tensor(distance,distanceFromRef), myatoms );
403 : }
404 : }
405 : }
406 : }
407 : }
408 96 : return std::log(sum)/lambda_;
409 : }
410 : }
411 :
412 17 : double EnvironmentSimilarity::maxDistance( std::vector<Vector> environment ) {
413 : double max_dist = 0.0;
414 85 : for(unsigned i=0; i<environment.size(); ++i) {
415 68 : double norm=environment[i].modulo();
416 68 : if (norm>max_dist) {
417 : max_dist=norm;
418 : }
419 : }
420 17 : return max_dist;
421 : }
422 :
423 10 : void EnvironmentSimilarity::parseReferenceEnvironments( std::vector<std::vector<Vector>>& environments, std::vector<std::vector<std::string>>& environmentsAtomNames, double& max_dist) {
424 : std::vector<double> lattice_constants;
425 20 : parseVector("LATTICE_CONSTANTS", lattice_constants);
426 : std::string crystal_structure;
427 20 : parse("CRYSTAL_STRUCTURE", crystal_structure);
428 : // find crystal structure
429 10 : if (crystal_structure == "FCC") {
430 1 : if (lattice_constants.size() != 1) {
431 0 : error("Number of LATTICE_CONSTANTS arguments must be one for FCC");
432 : }
433 1 : environments.resize(1);
434 1 : environments[0].resize(12);
435 1 : environments[0][0] = Vector(+0.5,+0.5,+0.0)*lattice_constants[0];
436 1 : environments[0][1] = Vector(-0.5,-0.5,+0.0)*lattice_constants[0];
437 1 : environments[0][2] = Vector(+0.5,-0.5,+0.0)*lattice_constants[0];
438 1 : environments[0][3] = Vector(-0.5,+0.5,+0.0)*lattice_constants[0];
439 1 : environments[0][4] = Vector(+0.5,+0.0,+0.5)*lattice_constants[0];
440 1 : environments[0][5] = Vector(-0.5,+0.0,-0.5)*lattice_constants[0];
441 1 : environments[0][6] = Vector(-0.5,+0.0,+0.5)*lattice_constants[0];
442 1 : environments[0][7] = Vector(+0.5,+0.0,-0.5)*lattice_constants[0];
443 1 : environments[0][8] = Vector(+0.0,+0.5,+0.5)*lattice_constants[0];
444 1 : environments[0][9] = Vector(+0.0,-0.5,-0.5)*lattice_constants[0];
445 1 : environments[0][10] = Vector(+0.0,-0.5,+0.5)*lattice_constants[0];
446 1 : environments[0][11] = Vector(+0.0,+0.5,-0.5)*lattice_constants[0];
447 1 : max_dist = std::sqrt(2)*lattice_constants[0]/2.;
448 9 : } else if (crystal_structure == "SC") {
449 0 : if (lattice_constants.size() != 1) {
450 0 : error("Number of LATTICE_CONSTANTS arguments must be one for SC");
451 : }
452 0 : environments.resize(1);
453 0 : environments[0].resize(6);
454 0 : environments[0][0] = Vector(+1.0,+0.0,+0.0)*lattice_constants[0];
455 0 : environments[0][1] = Vector(-1.0,+0.0,+0.0)*lattice_constants[0];
456 0 : environments[0][2] = Vector(+0.0,+1.0,+0.0)*lattice_constants[0];
457 0 : environments[0][3] = Vector(+0.0,-1.0,+0.0)*lattice_constants[0];
458 0 : environments[0][4] = Vector(+0.0,+0.0,+1.0)*lattice_constants[0];
459 0 : environments[0][5] = Vector(+0.0,+0.0,-1.0)*lattice_constants[0];
460 0 : max_dist = lattice_constants[0];
461 9 : } else if (crystal_structure == "BCC") {
462 2 : if (lattice_constants.size() != 1) {
463 0 : error("Number of LATTICE_CONSTANTS arguments must be one for BCC");
464 : }
465 2 : environments.resize(1);
466 2 : environments[0].resize(14);
467 2 : environments[0][0] = Vector(+0.5,+0.5,+0.5)*lattice_constants[0];
468 2 : environments[0][1] = Vector(-0.5,-0.5,-0.5)*lattice_constants[0];
469 2 : environments[0][2] = Vector(-0.5,+0.5,+0.5)*lattice_constants[0];
470 2 : environments[0][3] = Vector(+0.5,-0.5,+0.5)*lattice_constants[0];
471 2 : environments[0][4] = Vector(+0.5,+0.5,-0.5)*lattice_constants[0];
472 2 : environments[0][5] = Vector(-0.5,-0.5,+0.5)*lattice_constants[0];
473 2 : environments[0][6] = Vector(+0.5,-0.5,-0.5)*lattice_constants[0];
474 2 : environments[0][7] = Vector(-0.5,+0.5,-0.5)*lattice_constants[0];
475 2 : environments[0][8] = Vector(+1.0,+0.0,+0.0)*lattice_constants[0];
476 2 : environments[0][9] = Vector(+0.0,+1.0,+0.0)*lattice_constants[0];
477 2 : environments[0][10] = Vector(+0.0,+0.0,+1.0)*lattice_constants[0];
478 2 : environments[0][11] = Vector(-1.0,+0.0,+0.0)*lattice_constants[0];
479 2 : environments[0][12] = Vector(+0.0,-1.0,+0.0)*lattice_constants[0];
480 2 : environments[0][13] = Vector(+0.0,+0.0,-1.0)*lattice_constants[0];
481 2 : max_dist = lattice_constants[0];
482 7 : } else if (crystal_structure == "HCP") {
483 1 : if (lattice_constants.size() != 2) {
484 0 : error("Number of LATTICE_CONSTANTS arguments must be two for HCP");
485 : }
486 1 : environments.resize(2);
487 1 : environments[0].resize(12);
488 1 : environments[1].resize(12);
489 : double sqrt3=std::sqrt(3);
490 1 : environments[0][0] = Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
491 1 : environments[0][1] = Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
492 1 : environments[0][2] = Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
493 1 : environments[0][3] = Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
494 1 : environments[0][4] = Vector(+1.0,+0.0,+0.0) *lattice_constants[0];
495 1 : environments[0][5] = Vector(-1.0,+0.0,+0.0) *lattice_constants[0];
496 1 : environments[0][6] = Vector(+0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
497 1 : environments[0][7] = Vector(-0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
498 1 : environments[0][8] = Vector(+0.0,-sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
499 1 : environments[0][9] = Vector(+0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
500 1 : environments[0][10] = Vector(-0.5,+sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
501 1 : environments[0][11] = Vector(+0.0,-sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
502 1 : environments[1][0] = Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
503 1 : environments[1][1] = Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
504 1 : environments[1][2] = Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
505 1 : environments[1][3] = Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
506 1 : environments[1][4] = Vector(+1.0,+0.0,+0.0) *lattice_constants[0];
507 1 : environments[1][5] = Vector(-1.0,+0.0,+0.0) *lattice_constants[0];
508 1 : environments[1][6] = Vector(+0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
509 1 : environments[1][7] = Vector(-0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
510 1 : environments[1][8] = Vector(+0.0,+sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,+0.5)*lattice_constants[1];
511 1 : environments[1][9] = Vector(+0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
512 1 : environments[1][10] = Vector(-0.5,-sqrt3/6.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
513 1 : environments[1][11] = Vector(+0.0,+sqrt3/3.0,+0.0)*lattice_constants[0] + Vector(+0.0,+0.0,-0.5)*lattice_constants[1];
514 1 : max_dist = lattice_constants[0];
515 6 : } else if (crystal_structure == "DIAMOND") {
516 1 : if (lattice_constants.size() != 1) {
517 0 : error("Number of LATTICE_CONSTANTS arguments must be one for DIAMOND");
518 : }
519 1 : environments.resize(2);
520 1 : environments[0].resize(4);
521 1 : environments[1].resize(4);
522 1 : environments[0][0] = Vector(+1.0,+1.0,+1.0)*lattice_constants[0]/4.0;
523 1 : environments[0][1] = Vector(-1.0,-1.0,+1.0)*lattice_constants[0]/4.0;
524 1 : environments[0][2] = Vector(+1.0,-1.0,-1.0)*lattice_constants[0]/4.0;
525 1 : environments[0][3] = Vector(-1.0,+1.0,-1.0)*lattice_constants[0]/4.0;
526 1 : environments[1][0] = Vector(+1.0,-1.0,+1.0)*lattice_constants[0]/4.0;
527 1 : environments[1][1] = Vector(-1.0,+1.0,+1.0)*lattice_constants[0]/4.0;
528 1 : environments[1][2] = Vector(+1.0,+1.0,-1.0)*lattice_constants[0]/4.0;
529 1 : environments[1][3] = Vector(-1.0,-1.0,-1.0)*lattice_constants[0]/4.0;
530 1 : max_dist = std::sqrt(3)*lattice_constants[0]/4.0;
531 5 : } else if (crystal_structure == "CUSTOM") {
532 : std::string reffile;
533 10 : parse("REFERENCE",reffile);
534 5 : if (!reffile.empty()) {
535 : // Case with one reference environment
536 1 : environments.resize(1);
537 1 : environmentsAtomNames.resize(1);
538 1 : PDB pdb;
539 2 : if( !pdb.read(reffile,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) ) {
540 0 : error("missing input file " + reffile );
541 : }
542 1 : unsigned natoms=pdb.getPositions().size();
543 1 : environments[0].resize( natoms );
544 1 : environmentsAtomNames[0].resize( natoms );
545 5 : for(unsigned i=0; i<natoms; ++i) {
546 4 : environments[0][i]=pdb.getPositions()[i];
547 : }
548 5 : for(unsigned i=0; i<natoms; ++i) {
549 8 : environmentsAtomNames[0][i]=pdb.getAtomName(pdb.getAtomNumbers()[i]);
550 : }
551 1 : max_dist=maxDistance(environments[0]);
552 1 : log.printf(" reading %d reference vectors from %s \n", natoms, reffile.c_str() );
553 1 : } else {
554 : // Case with several reference environments
555 4 : max_dist=0;
556 16 : for(unsigned int i=1;; i++) {
557 40 : if(!parseNumbered("REFERENCE_",i,reffile) ) {
558 : break;
559 : }
560 16 : PDB pdb;
561 32 : if( !pdb.read(reffile,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) ) {
562 0 : error("missing input file " + reffile );
563 : }
564 16 : unsigned natoms=pdb.getPositions().size();
565 : std::vector<Vector> environment;
566 : std::vector<std::string> environmentAtomNames;
567 16 : environment.resize( natoms );
568 16 : environmentAtomNames.resize( natoms);
569 80 : for(unsigned i=0; i<natoms; ++i) {
570 64 : environment[i]=pdb.getPositions()[i];
571 : }
572 80 : for(unsigned i=0; i<natoms; ++i) {
573 128 : environmentAtomNames[i]=pdb.getAtomName(pdb.getAtomNumbers()[i]);
574 : }
575 16 : environments.push_back(environment);
576 16 : environmentsAtomNames.push_back(environmentAtomNames);
577 16 : double norm = maxDistance(environment);
578 16 : if (norm>max_dist) {
579 4 : max_dist=norm;
580 : }
581 16 : log.printf(" Reference environment %d : reading %d reference vectors from %s \n", i, natoms, reffile.c_str() );
582 32 : }
583 : }
584 5 : if (environments.size()==0)
585 0 : error("No environments have been found! Please specify a PDB file in the REFERENCE "
586 : "or in the REFERENCE_1, REFERENCE_2, etc keywords");
587 5 : log.printf(" Number of reference environments is %lu\n",environments.size() );
588 5 : log.printf(" Number of vectors per reference environment is %lu\n",environments[0].size() );
589 : std::string atomNamesFile;
590 10 : parse("ATOM_NAMES_FILE",atomNamesFile);
591 5 : if (!atomNamesFile.empty()) {
592 1 : PDB pdb;
593 2 : if( !pdb.read(atomNamesFile,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) ) {
594 0 : error("missing input file " + atomNamesFile);
595 : }
596 1 : unsigned natoms=pdb.getPositions().size();
597 1 : atomNames_.resize( natoms );
598 385 : for(unsigned i=0; i<natoms; ++i) {
599 768 : atomNames_[i]=pdb.getAtomName(pdb.getAtomNumbers()[i]);
600 : }
601 1 : log.printf(" Read %d atoms from atom names file %s \n", natoms, atomNamesFile.c_str() );
602 1 : log.printf(" Atoms in environments will be considered only if atom names match.\n");
603 1 : } else {
604 4 : log.printf(" No atom names file has been given. Atoms in reference environments will be matched disregarding atom type.\n");
605 : }
606 : } else {
607 0 : error("CRYSTAL_STRUCTURE=" + crystal_structure + " does not match any structures in the database");
608 : }
609 :
610 10 : log.printf(" targeting the %s crystal structure",crystal_structure.c_str());
611 10 : if (lattice_constants.size()>0) {
612 5 : log.printf(" with lattice constants %f\n",lattice_constants[0]);
613 : } else {
614 5 : log.printf("\n");
615 : }
616 :
617 10 : log.printf(" maximum distance in the reference environment is %f\n",max_dist);
618 10 : }
619 :
620 : }
621 : }
|
