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 : \par Examples
100 :
101 : The following input calculates the ENVIRONMENTSIMILARITY kernel for 250 atoms in the system
102 : using the BCC atomic environment as target, and then calculates and prints the average value
103 : for this quantity.
104 :
105 : \plumedfile
106 : ENVIRONMENTSIMILARITY SPECIES=1-250 SIGMA=0.05 LATTICE_CONSTANTS=0.423 CRYSTAL_STRUCTURE=BCC MEAN LABEL=es
107 :
108 : PRINT ARG=es.mean FILE=COLVAR
109 : \endplumedfile
110 :
111 : The next example compares the environments of the 96 selected atoms with a user specified reference
112 : environment. The reference environment is contained in the env1.pdb file. Once the kernel is computed
113 : the average and the number of atoms with a kernel larger than 0.5 are computed.
114 :
115 : \plumedfile
116 : ENVIRONMENTSIMILARITY ...
117 : SPECIES=1-288:3
118 : SIGMA=0.05
119 : CRYSTAL_STRUCTURE=CUSTOM
120 : REFERENCE=env1.pdb
121 : LABEL=es
122 : MEAN
123 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
124 : ... ENVIRONMENTSIMILARITY
125 :
126 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
127 : \endplumedfile
128 :
129 : The next example is similar to the one above but in this case 4 reference environments are specified.
130 : Each reference environment is given in a separate pdb file.
131 :
132 : \plumedfile
133 : ENVIRONMENTSIMILARITY ...
134 : SPECIES=1-288:3
135 : SIGMA=0.05
136 : CRYSTAL_STRUCTURE=CUSTOM
137 : REFERENCE_1=env1.pdb
138 : REFERENCE_2=env2.pdb
139 : REFERENCE_3=env3.pdb
140 : REFERENCE_4=env4.pdb
141 : LABEL=es
142 : MEAN
143 : MORE_THAN={RATIONAL R_0=0.5 NN=12 MM=24}
144 : ... ENVIRONMENTSIMILARITY
145 :
146 : PRINT ARG=es.mean,es.morethan FILE=COLVAR
147 : \endplumedfile
148 :
149 : */
150 : //+ENDPLUMEDOC
151 :
152 :
153 : class EnvironmentSimilarity : public multicolvar::MultiColvarBase {
154 : private:
155 : // All global variables end with underscore
156 : // square of cutoff, square of broadening parameter
157 : double rcut2_, sigmaSqr_;
158 : // lambda parameter for softmax function
159 : double lambda_;
160 : // Array of Vectors to store the reference environments, i.e. the templates
161 : std::vector<std::vector<Vector>> environments_;
162 : public:
163 : static void registerKeywords( Keywords& keys );
164 : explicit EnvironmentSimilarity(const ActionOptions&);
165 : // active methods:
166 : virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
167 : // Returns the number of coordinates of the field
168 15 : bool isPeriodic() { return false; }
169 : // Calculates maximum distance in an environment
170 : double maxDistance(std::vector<Vector> environment);
171 : // Parse everything connected to the definition of the reference environments
172 : // First argument is the array of Vectors that stores the reference environments
173 : // Second argument is the maximum distance in the ref environments and sets the
174 : // cutoff for the cell lists
175 : void parseReferenceEnvironments( std::vector<std::vector<Vector>>& environments, double& max_dist);
176 : };
177 :
178 10437 : PLUMED_REGISTER_ACTION(EnvironmentSimilarity,"ENVIRONMENTSIMILARITY")
179 :
180 10 : void EnvironmentSimilarity::registerKeywords( Keywords& keys ) {
181 10 : MultiColvarBase::registerKeywords( keys );
182 30 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
183 20 : keys.add("compulsory","SIGMA","0.1","Broadening parameter");
184 30 : keys.add("compulsory","CRYSTAL_STRUCTURE","FCC","Targeted crystal structure. Options are: "
185 : "SC: simple cubic, "
186 : "BCC: body center cubic, "
187 : "FCC: face centered cubic, "
188 : "HCP: hexagonal closed pack, "
189 : "DIAMOND: cubic diamond, "
190 : "CUSTOM: user defined "
191 : " ");
192 30 : keys.add("optional","LATTICE_CONSTANTS","Lattice constants. Two comma separated values for HCP, "
193 : "one value for all other CRYSTAL_STRUCTURES.");
194 20 : keys.add("compulsory","LAMBDA","100","Lambda parameter");
195 20 : keys.add("optional","REFERENCE","PDB file with relative distances from central atom."
196 : " Use this keyword if you are targeting a single reference environment.");
197 20 : keys.add("numbered","REFERENCE_","PDB files with relative distances from central atom."
198 : " Each file corresponds to one template."
199 : " Use these keywords if you are targeting more than one reference environment.");
200 : // Use actionWithDistributionKeywords
201 40 : keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MAX");
202 40 : keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
203 30 : keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
204 10 : }
205 :
206 9 : EnvironmentSimilarity::EnvironmentSimilarity(const ActionOptions&ao):
207 : Action(ao),
208 9 : MultiColvarBase(ao)
209 : {
210 9 : log.printf(" Please read and cite ");
211 18 : log << plumed.cite("Piaggi and Parrinello, J. Chem. Phys. 150 (24), 244119 (2019)");
212 9 : log.printf("\n");
213 :
214 : // Parse everything connected to the definition of the reference environments
215 : double max_dist_ref_vector;
216 9 : parseReferenceEnvironments(environments_, max_dist_ref_vector);
217 :
218 : double sigma;
219 9 : parse("SIGMA", sigma);
220 9 : log.printf(" representing local density as a sum of Gaussians with standard deviation %f\n",sigma);
221 9 : sigmaSqr_=sigma*sigma;
222 :
223 9 : lambda_=100;
224 18 : parse("LAMBDA", lambda_);
225 9 : if (environments_.size()>1) log.printf(" using a soft max function with lambda %f\n",lambda_);
226 :
227 : // Set the link cell cutoff
228 9 : double rcut = max_dist_ref_vector + 3*sigma;
229 9 : setLinkCellCutoff( rcut );
230 9 : rcut2_ = rcut * rcut;
231 :
232 : // And setup the ActionWithVessel
233 9 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms ); checkRead();
234 9 : }
235 :
236 30048 : double EnvironmentSimilarity::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
237 30048 : if (environments_.size()==1) {
238 : // One reference environment case
239 155464 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
240 : Vector& distance=myatoms.getPosition(i);
241 : double d2;
242 232452 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
243 77840 : (d2+=distance[1]*distance[1])<rcut2_ &&
244 194068 : (d2+=distance[2]*distance[2])<rcut2_ &&
245 : d2>epsilon ) {
246 : // Iterate over atoms in the reference environment
247 236856 : for(unsigned k=0; k<environments_[0].size(); ++k) {
248 220400 : Vector distanceFromRef=distance-environments_[0][k];
249 220400 : double value = std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[0].size() ;
250 : // CAREFUL! Off-diagonal virial is incorrect. Do not perform NPT simulations with flexible box angles.
251 220400 : accumulateSymmetryFunction( 1, i, value, -(value/(2*sigmaSqr_))*distanceFromRef, (value/(2*sigmaSqr_))*Tensor(distance,distanceFromRef), myatoms );
252 : }
253 : }
254 : }
255 : return myatoms.getValue(1);
256 : } else {
257 : // More than one reference environment case
258 29196 : std::vector<double> values(environments_.size()); //value for each template
259 : // First time calculate sums
260 2808756 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
261 : Vector& distance=myatoms.getPosition(i);
262 : double d2;
263 4534472 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
264 1754912 : (d2+=distance[1]*distance[1])<rcut2_ &&
265 3469608 : (d2+=distance[2]*distance[2])<rcut2_ &&
266 : d2>epsilon ) {
267 : // Iterate over templates
268 1216388 : for(unsigned j=0; j<environments_.size(); ++j) {
269 : // Iterate over atoms in the template
270 4893780 : for(unsigned k=0; k<environments_[j].size(); ++k) {
271 3921936 : Vector distanceFromRef=distance-environments_[j][k];
272 3921936 : values[j] += std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
273 : }
274 : }
275 : }
276 : }
277 : double sum=0;
278 145188 : for(unsigned j=0; j<environments_.size(); ++j) {
279 115992 : values[j] = std::exp(lambda_*values[j]);
280 115992 : sum += values[j];
281 : }
282 : // Second time find derivatives
283 2808756 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
284 : Vector& distance=myatoms.getPosition(i);
285 : double d2;
286 4534472 : if ( (d2=distance[0]*distance[0])<rcut2_ &&
287 1754912 : (d2+=distance[1]*distance[1])<rcut2_ &&
288 3469608 : (d2+=distance[2]*distance[2])<rcut2_ &&
289 : d2>epsilon ) {
290 : // Iterate over reference environment
291 1216388 : for(unsigned j=0; j<environments_.size(); ++j) {
292 : // Iterate over atoms in the reference environment
293 4893780 : for(unsigned k=0; k<environments_[j].size(); ++k) {
294 3921936 : Vector distanceFromRef=distance-environments_[j][k];
295 3921936 : double value = std::exp(-distanceFromRef.modulo2()/(4*sigmaSqr_) )/environments_[j].size() ;
296 3921936 : accumulateSymmetryFunction( 1, i, value, -(values[j]/sum)*(value/(2*sigmaSqr_))*distanceFromRef, (values[j]/sum)*(value/(2*sigmaSqr_))*Tensor(distance,distanceFromRef), myatoms );
297 : }
298 : }
299 : }
300 : }
301 29196 : if(sum==0) sum=std::numeric_limits<double>::min();
302 29196 : return std::log(sum)/lambda_;
303 : }
304 : }
305 :
306 13 : double EnvironmentSimilarity::maxDistance( std::vector<Vector> environment ) {
307 : double max_dist = 0.0;
308 65 : for(unsigned i=0; i<environment.size(); ++i) {
309 52 : double norm=environment[i].modulo();
310 52 : if (norm>max_dist) max_dist=norm;
311 : }
312 13 : return max_dist;
313 : }
314 :
315 9 : void EnvironmentSimilarity::parseReferenceEnvironments( std::vector<std::vector<Vector>>& environments, double& max_dist) {
316 : std::vector<double> lattice_constants;
317 18 : parseVector("LATTICE_CONSTANTS", lattice_constants);
318 : std::string crystal_structure;
319 18 : parse("CRYSTAL_STRUCTURE", crystal_structure);
320 : // find crystal structure
321 9 : if (crystal_structure == "FCC") {
322 1 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for FCC");
323 1 : environments.resize(1);
324 1 : environments[0].resize(12);
325 1 : environments[0][0] = Vector(+0.5,+0.5,+0.0)*lattice_constants[0];
326 1 : environments[0][1] = Vector(-0.5,-0.5,+0.0)*lattice_constants[0];
327 1 : environments[0][2] = Vector(+0.5,-0.5,+0.0)*lattice_constants[0];
328 1 : environments[0][3] = Vector(-0.5,+0.5,+0.0)*lattice_constants[0];
329 1 : environments[0][4] = Vector(+0.5,+0.0,+0.5)*lattice_constants[0];
330 1 : environments[0][5] = Vector(-0.5,+0.0,-0.5)*lattice_constants[0];
331 1 : environments[0][6] = Vector(-0.5,+0.0,+0.5)*lattice_constants[0];
332 1 : environments[0][7] = Vector(+0.5,+0.0,-0.5)*lattice_constants[0];
333 1 : environments[0][8] = Vector(+0.0,+0.5,+0.5)*lattice_constants[0];
334 1 : environments[0][9] = Vector(+0.0,-0.5,-0.5)*lattice_constants[0];
335 1 : environments[0][10] = Vector(+0.0,-0.5,+0.5)*lattice_constants[0];
336 1 : environments[0][11] = Vector(+0.0,+0.5,-0.5)*lattice_constants[0];
337 1 : max_dist = std::sqrt(2)*lattice_constants[0]/2.;
338 8 : } else if (crystal_structure == "SC") {
339 0 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for SC");
340 0 : environments.resize(1);
341 0 : environments[0].resize(6);
342 0 : environments[0][0] = Vector(+1.0,+0.0,+0.0)*lattice_constants[0];
343 0 : environments[0][1] = Vector(-1.0,+0.0,+0.0)*lattice_constants[0];
344 0 : environments[0][2] = Vector(+0.0,+1.0,+0.0)*lattice_constants[0];
345 0 : environments[0][3] = Vector(+0.0,-1.0,+0.0)*lattice_constants[0];
346 0 : environments[0][4] = Vector(+0.0,+0.0,+1.0)*lattice_constants[0];
347 0 : environments[0][5] = Vector(+0.0,+0.0,-1.0)*lattice_constants[0];
348 0 : max_dist = lattice_constants[0];
349 8 : } else if (crystal_structure == "BCC") {
350 2 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for BCC");
351 2 : environments.resize(1);
352 2 : environments[0].resize(14);
353 2 : environments[0][0] = Vector(+0.5,+0.5,+0.5)*lattice_constants[0];
354 2 : environments[0][1] = Vector(-0.5,-0.5,-0.5)*lattice_constants[0];
355 2 : environments[0][2] = Vector(-0.5,+0.5,+0.5)*lattice_constants[0];
356 2 : environments[0][3] = Vector(+0.5,-0.5,+0.5)*lattice_constants[0];
357 2 : environments[0][4] = Vector(+0.5,+0.5,-0.5)*lattice_constants[0];
358 2 : environments[0][5] = Vector(-0.5,-0.5,+0.5)*lattice_constants[0];
359 2 : environments[0][6] = Vector(+0.5,-0.5,-0.5)*lattice_constants[0];
360 2 : environments[0][7] = Vector(-0.5,+0.5,-0.5)*lattice_constants[0];
361 2 : environments[0][8] = Vector(+1.0,+0.0,+0.0)*lattice_constants[0];
362 2 : environments[0][9] = Vector(+0.0,+1.0,+0.0)*lattice_constants[0];
363 2 : environments[0][10] = Vector(+0.0,+0.0,+1.0)*lattice_constants[0];
364 2 : environments[0][11] = Vector(-1.0,+0.0,+0.0)*lattice_constants[0];
365 2 : environments[0][12] = Vector(+0.0,-1.0,+0.0)*lattice_constants[0];
366 2 : environments[0][13] = Vector(+0.0,+0.0,-1.0)*lattice_constants[0];
367 2 : max_dist = lattice_constants[0];
368 6 : } else if (crystal_structure == "HCP") {
369 1 : if (lattice_constants.size() != 2) error("Number of LATTICE_CONSTANTS arguments must be two for HCP");
370 1 : environments.resize(2);
371 1 : environments[0].resize(12);
372 1 : environments[1].resize(12);
373 : double sqrt3=std::sqrt(3);
374 1 : environments[0][0] = Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
375 1 : environments[0][1] = Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
376 1 : environments[0][2] = Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
377 1 : environments[0][3] = Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
378 1 : environments[0][4] = Vector(+1.0,+0.0,+0.0) *lattice_constants[0];
379 1 : environments[0][5] = Vector(-1.0,+0.0,+0.0) *lattice_constants[0];
380 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];
381 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];
382 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];
383 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];
384 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];
385 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];
386 1 : environments[1][0] = Vector(+0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
387 1 : environments[1][1] = Vector(-0.5,+sqrt3/2.0,+0.0)*lattice_constants[0];
388 1 : environments[1][2] = Vector(+0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
389 1 : environments[1][3] = Vector(-0.5,-sqrt3/2.0,+0.0)*lattice_constants[0];
390 1 : environments[1][4] = Vector(+1.0,+0.0,+0.0) *lattice_constants[0];
391 1 : environments[1][5] = Vector(-1.0,+0.0,+0.0) *lattice_constants[0];
392 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];
393 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];
394 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];
395 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];
396 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];
397 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];
398 1 : max_dist = lattice_constants[0];
399 5 : } else if (crystal_structure == "DIAMOND") {
400 1 : if (lattice_constants.size() != 1) error("Number of LATTICE_CONSTANTS arguments must be one for DIAMOND");
401 1 : environments.resize(2);
402 1 : environments[0].resize(4); environments[1].resize(4);
403 1 : environments[0][0] = Vector(+1.0,+1.0,+1.0)*lattice_constants[0]/4.0;
404 1 : environments[0][1] = Vector(-1.0,-1.0,+1.0)*lattice_constants[0]/4.0;
405 1 : environments[0][2] = Vector(+1.0,-1.0,-1.0)*lattice_constants[0]/4.0;
406 1 : environments[0][3] = Vector(-1.0,+1.0,-1.0)*lattice_constants[0]/4.0;
407 1 : environments[1][0] = Vector(+1.0,-1.0,+1.0)*lattice_constants[0]/4.0;
408 1 : environments[1][1] = Vector(-1.0,+1.0,+1.0)*lattice_constants[0]/4.0;
409 1 : environments[1][2] = Vector(+1.0,+1.0,-1.0)*lattice_constants[0]/4.0;
410 1 : environments[1][3] = Vector(-1.0,-1.0,-1.0)*lattice_constants[0]/4.0;
411 1 : max_dist = std::sqrt(3)*lattice_constants[0]/4.0;
412 4 : } else if (crystal_structure == "CUSTOM") {
413 : std::string reffile;
414 8 : parse("REFERENCE",reffile);
415 4 : if (!reffile.empty()) {
416 : // Case with one reference environment
417 1 : environments.resize(1);
418 1 : PDB pdb;
419 2 : if( !pdb.read(reffile,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) )
420 0 : error("missing input file " + reffile );
421 1 : unsigned natoms=pdb.getPositions().size(); environments[0].resize( natoms );
422 5 : for(unsigned i=0; i<natoms; ++i) environments[0][i]=pdb.getPositions()[i];
423 1 : max_dist=maxDistance(environments[0]);
424 1 : log.printf(" reading %d reference vectors from %s \n", natoms, reffile.c_str() );
425 1 : } else {
426 : // Case with several reference environments
427 3 : max_dist=0;
428 12 : for(unsigned int i=1;; i++) {
429 30 : if(!parseNumbered("REFERENCE_",i,reffile) ) {break;}
430 12 : PDB pdb;
431 24 : if( !pdb.read(reffile,plumed.getAtoms().usingNaturalUnits(),0.1/plumed.getAtoms().getUnits().getLength()) )
432 0 : error("missing input file " + reffile );
433 12 : unsigned natoms=pdb.getPositions().size(); std::vector<Vector> environment; environment.resize( natoms );
434 60 : for(unsigned i=0; i<natoms; ++i) environment[i]=pdb.getPositions()[i];
435 12 : environments.push_back(environment);
436 12 : double norm = maxDistance(environment);
437 12 : if (norm>max_dist) max_dist=norm;
438 12 : log.printf(" Reference environment %d : reading %d reference vectors from %s \n", i, natoms, reffile.c_str() );
439 12 : }
440 : }
441 4 : if (environments.size()==0) error("No environments have been found! Please specify a PDB file in the REFERENCE "
442 : "or in the REFERENCE_1, REFERENCE_2, etc keywords");
443 4 : log.printf(" Number of reference environments is %lu\n",environments.size() );
444 4 : log.printf(" Number of vectors per reference environment is %lu\n",environments[0].size() );
445 : } else {
446 0 : error("CRYSTAL_STRUCTURE=" + crystal_structure + " does not match any structures in the database");
447 : }
448 :
449 9 : log.printf(" targeting the %s crystal structure",crystal_structure.c_str());
450 9 : if (lattice_constants.size()>0) log.printf(" with lattice constants %f\n",lattice_constants[0]);
451 4 : else log.printf("\n");
452 :
453 9 : log.printf(" maximum distance in the reference environment is %f\n",max_dist);
454 9 : }
455 :
456 : }
457 : }
|
