31int main(
int argc,
char* argv[])
43 size_t thresholdEW = 0;
44 double thresholdEnergy = 0.0;
45 double thresholdForce = 0.0;
51 cout <<
"USAGE: " << argv[0] <<
" <mode> <t_ew> <t_en <t_f> "
52 "<elem1 <elem2 ...>>\n"
53 <<
" <mode> ... Either compare 2 NNPs (compare) or apply"
54 " threshold to existing comparison data set (select).\n"
55 <<
" If <mode> is \"select\":\n"
56 <<
" <t_ew> ... Extrapolation warning threshold.\n"
57 <<
" <t_en> ... Energy per atom threshold.\n"
58 <<
" <t_f> .... Force threshold.\n"
59 <<
" <elem1 <elem2 ...>> ... Element strings in data set "
61 <<
" Execute in directory with the data set file\n"
63 <<
" and 2 subdirectories\n"
66 <<
" each containing these NNP files:\n"
67 <<
" - input.nn (NNP settings)\n"
68 <<
" - scaling.data (symmetry function scaling data)\n"
69 <<
" - \"weights.%%03d.data\" (weights files)\n";
74 if (mode ==
"compare")
77 else if (mode ==
"select")
81 throw runtime_error(
"ERROR: Wrong number of arguments.\n");
83 thresholdEW = (size_t)atoi(argv[2]);
84 thresholdEnergy = atof(argv[3]);
85 thresholdForce = atof(argv[4]);
86 size_t numElements = argc - 5;
88 for (
size_t i = 6; i < numElements + 5; ++i)
96 throw runtime_error(
"ERROR: Unknown mode selected.\n");
99 MPI_Init(&argc, &argv);
100 MPI_Comm_size(MPI_COMM_WORLD, &numProcs);
101 MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
103 if ((mode ==
"compare") && (numProcs % 2 != 0))
105 throw runtime_error(
"ERROR: Please start with an even number of MPI"
108 else if ((mode ==
"select") && (numProcs != 1))
110 throw runtime_error(
"ERROR: Please start with a single MPI "
114 if (mode ==
"compare")
116 numWorkers = numProcs / 2;
118 myNNP = myRank / numWorkers;
119 MPI_Comm_split(MPI_COMM_WORLD, myNNP, myRank, &commNNP);
120 MPI_Comm_size(commNNP, &numProcsNNP);
121 MPI_Comm_rank(commNNP, &myRankNNP);
124 myData = myRankNNP % numWorkers;
125 MPI_Comm_split(MPI_COMM_WORLD, myData, myRank, &commData);
126 MPI_Comm_size(commData, &numProcsData);
127 MPI_Comm_rank(commData, &myRankData);
131 myLog.open(
strpr(
"nnp-comp2.log.%1d.%04d",
136 string dirNNP =
strpr(
"nnp-data-%1d", myNNP + 1);
139 dataset.
log <<
"*** 2-NNP COMPARISON INITIALIZATION *****"
140 "**************************************\n";
142 dataset.
log <<
strpr(
"Number of workers per NNP: %4d\n", numWorkers);
143 dataset.
log <<
strpr(
"NNP Id : %4d\n", myNNP);
144 dataset.
log <<
strpr(
"Data Id : %4d\n", myData);
145 dataset.
log <<
"My NNP directory : " + dirNNP +
"\n";
147 dataset.
log <<
"Global Communicator:\n";
148 dataset.
log <<
strpr(
"- numProcs: %4d\n", numProcs);
149 dataset.
log <<
strpr(
"- myRank : %4d\n", myRank);
151 dataset.
log <<
"NNP Communicator:\n";
152 dataset.
log <<
strpr(
"- numProcs: %4d\n", numProcsNNP);
153 dataset.
log <<
strpr(
"- myRank : %4d\n", myRankNNP);
155 dataset.
log <<
"Data Communicator:\n";
156 dataset.
log <<
strpr(
"- numProcs: %4d\n", numProcsData);
157 dataset.
log <<
strpr(
"- myRank : %4d\n", myRankData);
159 dataset.
log <<
"Starting NNP initialization...\n";
160 dataset.
log <<
"*****************************************"
161 "**************************************\n";
176 dataset.
log <<
"*** 2-NNP COMPARISON ********************"
177 "**************************************\n";
181 dataset.
log <<
"Calculating energies and forces for dataset.\n";
185 dataset.
log <<
"Calculating energies for dataset.\n";
188 for (vector<Structure>::iterator it = dataset.
structures.begin();
192#ifdef N2P2_NO_SF_GROUPS
205 it->clearNeighborList();
214 ofstream myComparisonData;
215 myComparisonData.open(
strpr(
"comp.data.%04d", myRankNNP).c_str());
217 myDiff.open(
strpr(
"diff.out.%04d", myRankNNP).c_str());
221 vector<string> title;
222 vector<string> colName;
223 vector<string> colInfo;
224 vector<size_t> colSize;
225 title.push_back(
"2-NNP EW, energy and force comparison.");
226 colSize.push_back(7);
227 colName.push_back(
"index");
228 colInfo.push_back(
"Structure index.");
229 colSize.push_back(10);
230 colName.push_back(
"nAtoms");
231 colInfo.push_back(
"Number of atoms in structure.");
232 colSize.push_back(10);
233 colName.push_back(
"EW_NNP1");
234 colInfo.push_back(
"Extrapolation warnings issued by NNP 1.");
235 colSize.push_back(10);
236 colName.push_back(
"EW_NNP2");
237 colInfo.push_back(
"Extrapolation warnings issued by NNP 2.");
238 colSize.push_back(16);
239 colName.push_back(
"E_NNP1");
240 colInfo.push_back(
"Energy prediction of NNP 1.");
241 colSize.push_back(16);
242 colName.push_back(
"E_NNP2");
243 colInfo.push_back(
"Energy prediction of NNP 2.");
244 colSize.push_back(16);
245 colName.push_back(
"FDiffMean");
246 colInfo.push_back(
"Mean absolute force difference.");
247 colSize.push_back(16);
248 colName.push_back(
"FDiffMax");
249 colInfo.push_back(
"Maximum absolute force difference.");
257 for (vector<Structure>::iterator it = dataset.
structures.begin();
260 vector<double> buffer;
261 size_t sizeBuffer = 2;
262 if (useForces) sizeBuffer += 3 * it->numAtoms;
263 buffer.resize(sizeBuffer, 0.0);
264 MPI_Recv(&(buffer.front()), buffer.size(), MPI_DOUBLE, 1, 0, commData, MPI_STATUS_IGNORE);
266 double eNNP2 = buffer.at(count++);
267 size_t ewNNP2 = (size_t)buffer.at(count++);
268 double meanForceDiff = 0.0;
269 double maxForceDiff = 0.0;
270 it->energyRef = it->energy - eNNP2;
273 for (vector<Atom>::iterator it2 = it->atoms.begin();
274 it2 != it->atoms.end(); ++it2)
276 for (
size_t i = 0; i < 3; ++i)
278 it2->fRef[i] = it2->f[i] - buffer.at(count++);
279 double const fAbsDiff = fabs(it2->fRef[i]);
280 meanForceDiff += fAbsDiff;
281 maxForceDiff = max(maxForceDiff, fAbsDiff);
284 meanForceDiff /= 3 * it->numAtoms;
287 myDiff <<
strpr(
"%7zu %10zu %10zu %10zu %16.8E %16.8E "
291 it->numElementsPresent,
297 it->numElementsPresent += ewNNP2;
298 it->comment =
strpr(
"EWSUM %zu ", it->numElementsPresent)
300 it->writeToFile(&myComparisonData);
302 myComparisonData.close();
306 else if (myRankData == 1)
308 for (vector<Structure>::const_iterator it = dataset.
structures.begin();
311 double ew = (double)it->numElementsPresent;
312 vector<double> buffer;
313 buffer.push_back(it->energy);
314 buffer.push_back(ew);
317 for (vector<Atom>::const_iterator it2 = it->atoms.begin();
318 it2 != it->atoms.end(); ++it2)
320 buffer.push_back(it2->f[0]);
321 buffer.push_back(it2->f[1]);
322 buffer.push_back(it2->f[2]);
325 MPI_Send(&(buffer.front()), buffer.size(), MPI_DOUBLE, 0, 0, commData);
329 MPI_Barrier(MPI_COMM_WORLD);
332 dataset.
log <<
"Difference structures written to \"comp.data\".\n";
334 dataset.
log <<
"Difference data collected in \"diff.out\".\n";
338 MPI_Comm_free(&commNNP);
339 MPI_Comm_free(&commData);
341 dataset.
log <<
"*****************************************"
342 "**************************************\n";
344 else if (mode ==
"select")
351 log <<
"*** 2-NNP COMPARISON SELECTION **********"
352 "**************************************\n";
354 log <<
strpr(
"Element string: %s\n", elements.c_str());
355 log <<
strpr(
"Threshold extrapolation warnings : %d\n",
357 log <<
strpr(
"Threshold energy difference per atom: %g\n",
359 log <<
strpr(
"Threshold force difference : %g\n",
361 log <<
"*****************************************"
362 "**************************************\n";
365 selectFile.open(
"select.out");
367 vector<string> title;
368 vector<string> colName;
369 vector<string> colInfo;
370 vector<size_t> colSize;
371 title.push_back(
"Selection results based on 2-NNP comparison.");
372 colSize.push_back(7);
373 colName.push_back(
"index");
374 colInfo.push_back(
"Structure index.");
375 colSize.push_back(5);
376 colName.push_back(
"fAny");
377 colInfo.push_back(
"Whether structure was selected.");
378 colSize.push_back(5);
379 colName.push_back(
"fEW");
380 colInfo.push_back(
"Structure selected due to extrapolation warning.");
381 colSize.push_back(5);
382 colName.push_back(
"fE");
383 colInfo.push_back(
"Structure selected due to energy difference.");
384 colSize.push_back(5);
385 colName.push_back(
"fF");
386 colInfo.push_back(
"Structure selected due to force difference.");
387 colSize.push_back(10);
388 colName.push_back(
"nEW");
389 colInfo.push_back(
"Number of extrapolation warnings (both NNPs "
391 colSize.push_back(10);
392 colName.push_back(
"nAtoms");
393 colInfo.push_back(
"Number of atoms in structure.");
394 colSize.push_back(16);
395 colName.push_back(
"Ediff");
396 colInfo.push_back(
"Energy difference per atom.");
397 colSize.push_back(10);
398 colName.push_back(
"nForce");
399 colInfo.push_back(
"Number of forces exceeding threshold.");
400 colSize.push_back(16);
401 colName.push_back(
"FDiffMean");
402 colInfo.push_back(
"Mean absolute force difference.");
403 colSize.push_back(16);
404 colName.push_back(
"FDiffMax");
405 colInfo.push_back(
"Maximum absolute force difference.");
415 ifstream comparisonData;
416 comparisonData.open(
"comp.data");
418 ofstream selectionData;
419 selectionData.open(
"comp-selection.data");
421 size_t countEWHit = 0;
422 size_t countEnergyHit = 0;
423 size_t countForceHit = 0;
424 size_t countStructures = 0;
425 size_t countSelected = 0;
426 while (comparisonData.peek() != EOF)
429 bool flagEnergy =
false;
430 bool flagForce =
false;
431 size_t countForce = 0;
432 double meanForceDiff = 0.0;
433 double maxForceDiff = 0.0;
435 structure.
index = countStructures++;
437 size_t countEW = (size_t)atoi(
split(s.
comment).at(1).c_str());
438 if (thresholdEW > 0 && countEW >= thresholdEW)
443 if (thresholdEnergy > 0 &&
449 if (thresholdForce > 0)
451 for (vector<Atom>::const_iterator it = s.
atoms.begin();
452 it != s.
atoms.end(); ++it)
454 for (
size_t i = 0; i < 3; ++i)
456 double const fAbsDiff = fabs(it->fRef[i]);
457 maxForceDiff = max(maxForceDiff, fAbsDiff);
458 meanForceDiff += fAbsDiff;
459 if (fAbsDiff > thresholdForce)
467 if (flagForce) countForceHit++;
469 if (flagEW || flagEnergy || flagForce)
473 log <<
strpr(
"Structure %7zu selected.\n", structure.
index);
475 selectFile <<
strpr(
"%7zu %5d %5d %5d %5d %10zu %10zu %16.8E "
476 "%10zu %16.8E %16.8E\n",
478 (
int)(flagEW || flagEnergy || flagForce),
492 log <<
"*****************************************"
493 "**************************************\n";
494 log <<
strpr(
"Total number of structures : %7zu\n", countStructures);
495 log <<
strpr(
"Number of selected structures: %7zu (%6.2f %)\n",
496 countSelected, countSelected * 100.0 / countStructures);
497 log <<
strpr(
"EW threshold exceeded : %7zu (%6.2f %)\n",
498 countEWHit, countEWHit * 100.0 / countStructures);
499 log <<
strpr(
"Energy threshold exceeded : %7zu (%6.2f %)\n",
500 countEnergyHit, countEnergyHit * 100.0 / countStructures);
501 log <<
strpr(
"Force threshold exceeded : %7zu (%6.2f %)\n",
502 countForceHit, countForceHit * 100.0 / countStructures);
503 log <<
"Selected structures written to \"comp-selection.data\".\n";
504 log <<
"*****************************************"
505 "**************************************\n";
Collect and process large data sets.
int distributeStructures(bool randomize, bool excludeRank0=false, std::string const &fileName="input.data")
Read data file and distribute structures among processors.
void setupMPI()
Initialize MPI with MPI_COMM_WORLD.
std::vector< Structure > structures
All structures in this dataset.
void combineFiles(std::string filePrefix) const
Combine individual MPI proc files to one.
void toNormalizedUnits()
Switch all structures to normalized units.
std::size_t registerElements(std::string const &elementLine)
Extract all elements and store in element map.
Logging class for library output.
void registerStreamPointer(std::ofstream *const &streamPointer)
Register new C++ ofstream pointer.
bool writeToStdout
Turn on/off output to stdout.
void addEnergyOffset(Structure &structure, bool ref=true)
Add atomic energy offsets to reference energy.
void initialize()
Write welcome message with version information.
void setupGeneric(std::string const &nnpDir="", bool skipNormalize=false, bool initialHardness=false)
Combine multiple setup routines and provide a basic NNP setup.
virtual void setupNeuralNetworkWeights(std::map< std::string, std::string > fileNameFormats=std::map< std::string, std::string >())
Set up neural network weights from files with given name format.
void calculateAtomicNeuralNetworks(Structure &structure, bool const derivatives, std::string id="")
Calculate atomic neural networks for all atoms in given structure.
double getMaxCutoffRadius() const
Getter for Mode::maxCutoffRadius.
bool useNormalization() const
Check if normalization is enabled.
void calculateEnergy(Structure &structure) const
Calculate potential energy for a given structure.
void calculateSymmetryFunctionGroups(Structure &structure, bool const derivatives)
Calculate all symmetry function groups for all atoms in given structure.
virtual void setupSymmetryFunctionScaling(std::string const &fileName="scaling.data")
Set up symmetry function scaling from file.
std::size_t getNumExtrapolationWarnings() const
Count total number of extrapolation warnings encountered for all elements and symmetry functions.
bool settingsKeywordExists(std::string const &keyword) const
Check if keyword was found in settings file.
void resetExtrapolationWarnings()
Erase all extrapolation warnings and reset counters.
void convertToPhysicalUnits(Structure &structure) const
Convert one structure to physical units.
void calculateSymmetryFunctions(Structure &structure, bool const derivatives)
Calculate all symmetry functions for all atoms in given structure.
void calculateForces(Structure &structure) const
Calculate forces for all atoms in given structure.
void setupSymmetryFunctionStatistics(bool collectStatistics, bool collectExtrapolationWarnings, bool writeExtrapolationWarnings, bool stopOnExtrapolationWarnings)
Set up symmetry function statistics collection.
void loadSettingsFile(std::string const &fileName="input.nn")
Open settings file and load all keywords into memory.
string strpr(const char *format,...)
String version of printf function.
vector< string > createFileHeader(vector< string > const &title, vector< size_t > const &colSize, vector< string > const &colName, vector< string > const &colInfo, char const &commentChar)
vector< string > split(string const &input, char delimiter)
Split string at each delimiter.
void appendLinesToFile(ofstream &file, vector< string > const lines)
Append multiple lines of strings to open file stream.
int main(int argc, char *argv[])
Storage for one atomic configuration.
void writeToFile(std::string const fileName="output.data", bool const ref=true, bool const append=false) const
Write configuration to file.
std::string comment
Structure comment.
void setElementMap(ElementMap const &elementMap)
Set element map of structure.
std::size_t index
Index number of this structure.
void readFromFile(std::string const fileName="input.data")
Read configuration from file.
double energyRef
Reference potential energy.
void reset()
Reset everything but elementMap.
std::size_t numAtoms
Total number of atoms present in this structure.
std::vector< Atom > atoms
Vector of all atoms in this structure.