n2p2 - A neural network potential package
Functions
nnp-dist.cpp File Reference
#include "ElementMap.h"
#include "Log.h"
#include "Structure.h"
#include "utility.h"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <limits>
#include <map>
#include <string>
#include <utility>
#include <vector>
Include dependency graph for nnp-dist.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 35 of file nnp-dist.cpp.

36{
37 if (argc < 4)
38 {
39 cout << "USAGE: " << argv[0] << " <rcut> <nbin> <adf> <elem1 "
40 << "<elem2 <elem3...>>>\n"
41 << " <rcut> .... Cutoff radius.\n"
42 << " <nbin> .... Number of histogram bins.\n"
43 << " <adf> ..... Calculate angle distribution (0/1).\n"
44 << " <elemN> ... Symbol for Nth element.\n"
45 << " Execute in directory with these NNP files present:\n"
46 << " - input.data (structure file)\n";
47 return 1;
48 }
49
50 ofstream logFile;
51 logFile.open("nnp-dist.log");
52 Log log;
53 log.registerStreamPointer(&logFile);
54
55 log << "\n";
56 log << "*** NNP-DIST ****************************"
57 "**************************************\n";
58 log << "\n";
59
60 size_t numElements = argc - 4;
61 log << strpr("Number of elements: %zu\n", numElements);
62 string elements;
63 elements += argv[4];
64 for (size_t i = 5; i < numElements + 4; ++i)
65 {
66 elements += " ";
67 elements += argv[i];
68 }
69 log << strpr("Element string : %s\n", elements.c_str());
70 double cutoffRadius = atof(argv[1]);
71 log << strpr("Cutoff radius : %f\n", cutoffRadius);
72 size_t numBins = (size_t)atoi(argv[2]);
73 log << strpr("Histogram bins : %zu\n", numBins);
74 bool calcAdf = (bool)atoi(argv[3]);
75 log << strpr("Calculate ADF : %d\n", calcAdf);
76
77 double minCellLength = numeric_limits<double>::max();
78 double maxCellLength = 0.0;
79 vector<size_t> periodicStructures;
80 vector<vector<double>> cellLengths(3);
81
82 size_t numRdf = numElements * (numElements + 1);
83 numRdf /= 2;
84 log << strpr("Number of RDFs : %zu\n", numRdf);
85 size_t numAdf = 0;
86 if (calcAdf)
87 {
88 numAdf = numElements * numElements * (numElements + 1);
89 numAdf /= 2;
90 }
91 log << strpr("Number of ADFs : %zu\n", numAdf);
92 log << "*****************************************"
93 "**************************************\n";
94
95 ElementMap elementMap;
96 elementMap.registerElements(elements);
97
98 double dr = cutoffRadius / numBins;
99 map<pair<size_t, size_t>, vector<double>* > rhist;
100 map<pair<size_t, size_t>, vector<double>* > rdf;
101 double da = 180.0 / numBins;
102 vector<map<pair<size_t, size_t>, vector<double>* > > ahist(numElements);
103 vector<map<pair<size_t, size_t>, vector<double>* > > adf(numElements);
104
105 for (size_t i = 0; i < numElements; ++i)
106 {
107 for (size_t j = i; j < numElements; ++j)
108 {
109 pair<size_t, size_t> e(i, j);
110 rhist[e] = new vector<double>(numBins, 0);
111 rdf [e] = new vector<double>(numBins, 0.0);
112 if (i != j)
113 {
114 rhist[make_pair(j, i)] = rhist[e];
115 rdf [make_pair(j, i)] = rdf [e];
116 }
117 }
118 }
119 if (calcAdf)
120 {
121 for (size_t i = 0; i < numElements; ++i)
122 {
123 for (size_t j = 0; j < numElements; ++j)
124 {
125 for (size_t k = j; k < numElements; ++k)
126 {
127 pair<size_t, size_t> e(j, k);
128 ahist.at(i)[e] = new vector<double>(numBins, 0);
129 adf.at(i) [e] = new vector<double>(numBins, 0.0);
130 if (j != k)
131 {
132 ahist.at(i)[make_pair(k, j)] = ahist.at(i)[e];
133 adf.at(i) [make_pair(k, j)] = adf.at(i) [e];
134 }
135 }
136 }
137 }
138 }
139
140 ifstream inputFile;
141 inputFile.open("input.data");
142 Structure structure;
143 structure.setElementMap(elementMap);
144
145 size_t countStructures = 0;
146 size_t countPeriodicStructures = 0;
147 vector<double> numberDensity(numElements, 0.0);
148 while (inputFile.peek() != EOF)
149 {
150 structure.readFromFile(inputFile);
151 structure.calculateNeighborList(cutoffRadius);
152 if (structure.isPeriodic)
153 {
154 periodicStructures.push_back(countStructures);
155 for (size_t i = 0; i < 3; ++i)
156 {
157 cellLengths.at(i).push_back(structure.box[i].norm());
158 }
159 }
160 for (vector<Atom>::const_iterator it = structure.atoms.begin();
161 it != structure.atoms.end(); ++it)
162 {
163 for (size_t j = 0; j < it->numNeighbors; ++j)
164 {
165 size_t const ei = it->element;
166 Atom::Neighbor const& nj = it->neighbors.at(j);
167 size_t const ej = nj.element;
168 double const rij = nj.d;
169 pair<size_t, size_t> e(ei, ej);
170 if (ei == ej)
171 {
172 rhist[e]->at((size_t)floor(rij / dr)) += 1.0;
173 }
174 else
175 {
176 rhist[e]->at((size_t)floor(rij / dr)) += 0.5;
177 }
178 if (calcAdf)
179 {
180 for (size_t k = j + 1; k < it->numNeighbors; ++k)
181 {
182 Atom::Neighbor const& nk = it->neighbors.at(k);
183 size_t const ek = nk.element;
184 e = make_pair(ej, ek);
185 double theta = nj.dr * nk.dr / (nj.d * nk.d);
186 // Use first bin for 0 degree angle and catch problems
187 // with rounding errors.
188 if (theta >= 1.0)
189 {
190 ahist.at(ei)[e]->at(0) += 1.0;
191 }
192 // Use last bin for 180 degree angle and catch problems
193 // with rounding errors.
194 else if (theta <= -1.0)
195 {
196 ahist.at(ei)[e]->at(numBins - 1) += 1.0;
197 }
198 else
199 {
200 theta = acos(theta) * 180.0 / M_PI;
201 ahist.at(ei)[e]
202 ->at((size_t)floor(theta / da)) += 1.0;
203 }
204 }
205 }
206 }
207 }
208 countStructures++;
209 if (structure.isPeriodic) countPeriodicStructures++;
210 double const volume = structure.volume;
211 for (size_t i = 0; i < numElements; ++i)
212 {
213 size_t const nni = structure.numAtomsPerElement.at(i);
214 if (structure.isPeriodic)
215 {
216 numberDensity.at(i) += nni / volume;
217 }
218 for (size_t j = i; j < numElements; ++j)
219 {
220 size_t const nnj = structure.numAtomsPerElement.at(j);
221 if (nni == 0 || nnj == 0) continue;
222 pair<size_t, size_t> e(i, j);
223 for (size_t n = 0; n < numBins; ++n)
224 {
225 double r = (n + 0.5) * dr;
226 if (structure.isPeriodic)
227 {
228 rhist[e]->at(n) *= volume / (nni * nnj);
229 }
230 rdf[e]->at(n) += rhist[e]->at(n)
231 / (4.0 * M_PI * r * r * dr);
232 }
233 rhist[e]->clear();
234 rhist[e]->resize(numBins, 0.0);
235 }
236 }
237 if (calcAdf)
238 {
239 for (size_t i = 0; i < numElements; ++i)
240 {
241 for (size_t j = 0; j < numElements; ++j)
242 {
243 for (size_t k = j; k < numElements; ++k)
244 {
245 pair<size_t, size_t> e(j, k);
246 double countAngles = 0.0;
247 for (size_t n = 0; n < numBins; ++n)
248 {
249 countAngles += ahist.at(i)[e]->at(n);
250 }
251 if (countAngles > 0)
252 {
253 for (size_t n = 0; n < numBins; ++n)
254 {
255 adf.at(i)[e]->at(n) += ahist.at(i)[e]->at(n)
256 / countAngles / da;
257 }
258 }
259 ahist.at(i)[e]->clear();
260 ahist.at(i)[e]->resize(numBins, 0.0);
261 }
262 }
263 }
264 }
265 log << strpr("Configuration %7zu: %7zu atoms\n",
266 countStructures,
267 structure.numAtoms);
268 structure.reset();
269 }
270 log << "*****************************************"
271 "**************************************\n";
272 log << strpr("Number of structures: %9zu\n", countStructures);
273 log << strpr("Number of periodic structures: %9zu\n",
274 countPeriodicStructures);
275 bool calcCn = false;
276 if (countStructures == countPeriodicStructures) calcCn = true;
277 if (calcCn)
278 {
279 vector<double>::const_iterator minnd = min_element(
280 numberDensity.begin(), numberDensity.end());
281 for (size_t i = 0; i < numElements; ++i)
282 {
283 log << strpr("Number density (ratio) of element %2s: %16.8E "
284 "(%.2f)\n", elementMap[i].c_str(),
285 numberDensity.at(i), numberDensity.at(i) / (*minnd));
286 }
287 }
288
289 ofstream outputFile;
290 if (countPeriodicStructures > 0)
291 {
292 outputFile.open("cell-lengths.out");
293
294 // File header.
295 vector<string> title;
296 vector<string> colName;
297 vector<string> colInfo;
298 vector<size_t> colSize;
299 title.push_back("Unit cell vector lengths for each periodic "
300 "structure.");
301 colSize.push_back(10);
302 colName.push_back("index");
303 colInfo.push_back("Index of periodic structure in data set "
304 "(starting with 1).");
305 colSize.push_back(16);
306 colName.push_back("norm(A)");
307 colInfo.push_back("Norm of 1st unit cell vector.");
308 colSize.push_back(16);
309 colName.push_back("norm(B)");
310 colInfo.push_back("Norm of 2nd unit cell vector.");
311 colSize.push_back(16);
312 colName.push_back("norm(C)");
313 colInfo.push_back("Norm of 3rd unit cell vector.");
314
315 appendLinesToFile(outputFile,
316 createFileHeader(title,
317 colSize,
318 colName,
319 colInfo));
320
321 for (size_t i = 0; i < periodicStructures.size(); ++i)
322 {
323 outputFile << strpr("%10zu %16.8E %16.8E %16.8E\n",
324 periodicStructures.at(i) + 1,
325 cellLengths.at(0).at(i),
326 cellLengths.at(1).at(i),
327 cellLengths.at(2).at(i));
328 for (size_t j = 0; j < 3; ++j)
329 {
330 minCellLength = min(minCellLength, cellLengths.at(j).at(i));
331 maxCellLength = max(maxCellLength, cellLengths.at(j).at(i));
332 }
333 }
334
335 outputFile.close();
336
337 log << strpr("Minimum unit cell vector length : %16.8E\n",
338 minCellLength);
339 log << strpr("Maximum unit cell vector length : %16.8E\n",
340 maxCellLength);
341 }
342 log << "*****************************************"
343 "**************************************\n";
344
345 for (size_t i = 0; i < numElements; ++i)
346 {
347 for (size_t j = i; j < numElements; ++j)
348 {
349 pair<size_t, size_t> e(i, j);
350 string fileName = strpr("rdf_%s_%s.out",
351 elementMap[i].c_str(),
352 elementMap[j].c_str());
353 log << strpr("Writing RDF for element combination (%2s/%2s) "
354 "to file %s.\n",
355 elementMap[i].c_str(),
356 elementMap[j].c_str(),
357 fileName.c_str());
358 outputFile.open(fileName.c_str());
359
360 // File header.
361 vector<string> title;
362 vector<string> colName;
363 vector<string> colInfo;
364 vector<size_t> colSize;
365 title.push_back(strpr("Radial distribution function for element "
366 "combination %2s-%2s.",
367 elementMap[i].c_str(),
368 elementMap[j].c_str()));
369 colSize.push_back(16);
370 colName.push_back("dist_bin_l");
371 colInfo.push_back("Distance, left bin limit.");
372 colSize.push_back(16);
373 colName.push_back("dist_bin_r");
374 colInfo.push_back("Distance, right bin limit.");
375 colSize.push_back(16);
376 colName.push_back("rdf");
377 colInfo.push_back("Radial distribution function, standard "
378 "normalization (-> 1 for r -> inf)");
379 colSize.push_back(16);
380 colName.push_back("rdf_max1");
381 colInfo.push_back("Radial distribution function, maximum "
382 "normalized to 1.");
383 if (calcCn)
384 {
385 colSize.push_back(16);
386 colName.push_back("cn");
387 colInfo.push_back("Coordination number.");
388 }
389 appendLinesToFile(outputFile,
390 createFileHeader(title,
391 colSize,
392 colName,
393 colInfo));
394
395 vector<double> cn(rdf[e]->size(), 0.0);
396 double integral = 0.0;
397 double const pre = 4.0 * M_PI * dr * 0.5
398 / countStructures / countStructures;
399 double maxRdf = *max_element(rdf[e]->begin(), rdf[e]->end());
400 for (size_t n = 0; n < numBins; ++n)
401 {
402 double const low = n * dr;
403 double const center = (n + 0.5) * dr;
404 double const high = (n + 1) * dr;
405 outputFile << strpr("%16.8E %16.8E %16.8E %16.8E",
406 low,
407 high,
408 rdf[e]->at(n) / countStructures,
409 rdf[e]->at(n) / maxRdf);
410 if (calcCn)
411 {
412 if (n > 0)
413 {
414 integral += pre * center * center * numberDensity.at(j)
415 * (rdf[e]->at(n-1) + rdf[e]->at(n));
416 }
417 outputFile << strpr(" %16.8E", integral);
418 }
419 outputFile << '\n';
420 }
421 outputFile.close();
422 delete rdf [e];
423 delete rhist[e];
424 }
425 }
426 if (calcAdf)
427 {
428 for (size_t i = 0; i < numElements; ++i)
429 {
430 for (size_t j = 0; j < numElements; ++j)
431 {
432 for (size_t k = j; k < numElements; ++k)
433 {
434 pair<size_t, size_t> e(j, k);
435 string fileName = strpr("adf_%s_%s_%s.out",
436 elementMap[i].c_str(),
437 elementMap[j].c_str(),
438 elementMap[k].c_str());
439 log << strpr("Writing ADF for element combination "
440 "(%2s/%2s/%2s) to file %s.\n",
441 elementMap[i].c_str(),
442 elementMap[j].c_str(),
443 elementMap[k].c_str(),
444 fileName.c_str());
445 outputFile.open(fileName.c_str());
446
447 // File header.
448 vector<string> title;
449 vector<string> colName;
450 vector<string> colInfo;
451 vector<size_t> colSize;
452 title.push_back(strpr("Angular distribution function for "
453 "element combination %2s-%2s-%2s.",
454 elementMap[i].c_str(),
455 elementMap[j].c_str(),
456 elementMap[k].c_str()));
457 colSize.push_back(16);
458 colName.push_back("angle_bin_l");
459 colInfo.push_back("Angle [degree], left bin limit.");
460 colSize.push_back(16);
461 colName.push_back("angle_bin_r");
462 colInfo.push_back("Angle [degree], right bin limit.");
463 colSize.push_back(16);
464 colName.push_back("adf");
465 colInfo.push_back("Angular distribution function, "
466 "probability normalization "
467 "(integral = 1)");
468 colSize.push_back(16);
469 colName.push_back("adf_max1");
470 colInfo.push_back("Angular distribution function, maximum "
471 "normalized to 1.");
472 appendLinesToFile(outputFile,
473 createFileHeader(title,
474 colSize,
475 colName,
476 colInfo));
477
478 double maxAdf = *max_element(adf.at(i)[e]->begin(),
479 adf.at(i)[e]->end());
480 for (size_t n = 0; n < numBins; ++n)
481 {
482 double const low = n * da;
483 double const high = (n + 1) * da;
484 outputFile << strpr("%16.8E %16.8E %16.8E %16.8E\n",
485 low,
486 high,
487 adf.at(i)[e]->at(n)
488 / countStructures,
489 adf.at(i)[e]->at(n) / maxAdf);
490 }
491 outputFile.close();
492 delete adf.at(i) [e];
493 delete ahist.at(i)[e];
494 }
495 }
496 }
497 }
498
499 log << "*****************************************"
500 "**************************************\n";
501 logFile.close();
502
503 return 0;
504}
nnp::ElementMap
Contains element map.
Definition: ElementMap.h:30
nnp::ElementMap::registerElements
std::size_t registerElements(std::string const &elementLine)
Extract all elements and store in element map.
Definition: ElementMap.cpp:36
nnp::Log
Logging class for library output.
Definition: Log.h:34
nnp::Log::registerStreamPointer
void registerStreamPointer(std::ofstream *const &streamPointer)
Register new C++ ofstream pointer.
Definition: Log.cpp:91
nnp::strpr
string strpr(const char *format,...)
String version of printf function.
Definition: utility.cpp:90
nnp::createFileHeader
vector< string > createFileHeader(vector< string > const &title, vector< size_t > const &colSize, vector< string > const &colName, vector< string > const &colInfo, char const &commentChar)
Definition: utility.cpp:104
nnp::appendLinesToFile
void appendLinesToFile(ofstream &file, vector< string > const lines)
Append multiple lines of strings to open file stream.
Definition: utility.cpp:218
logFile
ofstream logFile
Definition: nnp-cutoff.cpp:29
nnp::Atom::Neighbor
Struct to store information on neighbor atoms.
Definition: Atom.h:35
nnp::Atom::Neighbor::element
std::size_t element
Element index of neighbor atom.
Definition: Atom.h:41
nnp::Atom::Neighbor::d
double d
Distance to neighbor atom.
Definition: Atom.h:43
nnp::Atom::Neighbor::dr
Vec3D dr
Distance vector to neighbor atom.
Definition: Atom.h:45
nnp::Structure
Storage for one atomic configuration.
Definition: Structure.h:34
nnp::Structure::box
Vec3D box[3]
Simulation box vectors.
Definition: Structure.h:90
nnp::Structure::numAtomsPerElement
std::vector< std::size_t > numAtomsPerElement
Number of atoms of each element in this structure.
Definition: Structure.h:94
nnp::Structure::isPeriodic
bool isPeriodic
If periodic boundary conditions apply.
Definition: Structure.h:56
nnp::Structure::setElementMap
void setElementMap(ElementMap const &elementMap)
Set element map of structure.
Definition: Structure.cpp:60
nnp::Structure::readFromFile
void readFromFile(std::string const fileName="input.data")
Read configuration from file.
Definition: Structure.cpp:86
nnp::Structure::reset
void reset()
Reset everything but elementMap.
Definition: Structure.cpp:543
nnp::Structure::volume
double volume
Simulation box volume.
Definition: Structure.h:84
nnp::Structure::numAtoms
std::size_t numAtoms
Total number of atoms present in this structure.
Definition: Structure.h:68
nnp::Structure::calculateNeighborList
void calculateNeighborList(double cutoffRadius)
Calculate neighbor list for all atoms.
Definition: Structure.cpp:253
nnp::Structure::atoms
std::vector< Atom > atoms
Vector of all atoms in this structure.
Definition: Structure.h:96
nnp::Vec3D::norm
double norm() const
Calculate norm of vector.
Definition: Vec3D.h:256

References nnp::appendLinesToFile(), nnp::Structure::atoms, nnp::Structure::box, nnp::Structure::calculateNeighborList(), nnp::createFileHeader(), nnp::Atom::Neighbor::d, nnp::Atom::Neighbor::dr, nnp::Atom::Neighbor::element, nnp::Structure::isPeriodic, logFile, nnp::Vec3D::norm(), nnp::Structure::numAtoms, nnp::Structure::numAtomsPerElement, nnp::Structure::readFromFile(), nnp::ElementMap::registerElements(), nnp::Log::registerStreamPointer(), nnp::Structure::reset(), nnp::Structure::setElementMap(), nnp::strpr(), and nnp::Structure::volume.

Here is the call graph for this function: