18#include "Eigen/src/Core/Matrix.h"
30#include <gsl/gsl_rng.h>
40 parallelMode (PM_TRAIN_RK0 ),
41 jacobianMode (JM_SUM ),
42 updateStrategy (US_COMBINED ),
44 hasStructures (false ),
46 repeatedEnergyUpdates (false ),
48 writeTrainingLog (false ),
53 writeWeightsEvery (0 ),
54 writeWeightsAlways (0 ),
55 writeNeuronStatisticsEvery (0 ),
56 writeNeuronStatisticsAlways(0 ),
59 trainingLogFileName (
"train-log.out")
66 for (vector<Updater*>::iterator it =
updaters.begin();
85 log <<
"*** DEFINE TRAINING/TEST SETS ***********"
86 "**************************************\n";
89 vector<string> pCheck = {
"force",
"charge"};
93 check |= (find(pCheck.begin(), pCheck.end(), k) != pCheck.end());
97 double testSetFraction = atof(
settings[
"test_fraction"].c_str());
98 log <<
strpr(
"Desired test set ratio : %f\n", testSetFraction);
103 for (
size_t i = 0; i <
structures.size(); ++i)
109 double const r = gsl_rng_uniform(
rng);
116 k =
"energy";
if (
p.
exists(k))
p[k].numTestPatterns++;
117 k =
"force";
if (
p.
exists(k))
p[k].numTestPatterns += 3 * na;
120 k =
"charge";
if (
p.
exists(k))
p[k].numTestPatterns++;
124 k =
"charge";
if (
p.
exists(k))
p[k].numTestPatterns += na;
137 p[k].numTrainPatterns++;
139 p[k].updateCandidates.back().s = i;
149 p[k].updateCandidates.back().s = i;
151 for (
auto &ai : s.
atoms)
153 for (
size_t j = 0; j < 3; ++j)
158 p[k].updateCandidates.back().s = i;
160 p[k].updateCandidates.back().subCandidates
162 p[k].updateCandidates.back().subCandidates
163 .back().a = ai.index;
164 p[k].updateCandidates.back().subCandidates
174 p[k].numTrainPatterns++;
176 p[k].updateCandidates.back().s = i;
181 for (vector<Atom>::const_iterator it = s.
atoms.begin();
182 it != s.
atoms.end(); ++it)
185 p[k].updateCandidates.back().s = i;
186 p[k].updateCandidates.back().subCandidates
188 p[k].updateCandidates.back().subCandidates.back().a
197 log <<
strpr(
"WARNING: Structure %zu not assigned to either "
198 "training or test set.\n", s.
index);
205 log <<
strpr(
"WARNING: Process %d has no atoms of element "
214 MPI_Allreduce(MPI_IN_PLACE, &(
p[k].numTrainPatterns), 1,
MPI_SIZE_T, MPI_SUM,
comm);
215 MPI_Allreduce(MPI_IN_PLACE, &(
p[k].numTestPatterns) , 1,
MPI_SIZE_T, MPI_SUM,
comm);
216 double sum =
p[k].numTrainPatterns +
p[k].numTestPatterns;
217 log <<
"Training/test split of data set for property \"" + k +
"\":\n";
218 log <<
strpr(
"- Total patterns : %.0f\n", sum);
219 log <<
strpr(
"- Training patterns : %d\n",
p[k].numTrainPatterns);
220 log <<
strpr(
"- Test patterns : %d\n",
p[k].numTestPatterns);
221 log <<
strpr(
"- Test set fraction : %f\n",
p[k].numTestPatterns / sum);
224 log <<
"*****************************************"
225 "**************************************\n";
233 log <<
"*** WRITE TRAINING/TEST SETS ************"
234 "**************************************\n";
237 string fileName =
strpr(
"train.data.%04d",
myRank);
239 fileTrain.open(fileName.c_str());
240 if (!fileTrain.is_open())
242 runtime_error(
strpr(
"ERROR: Could not open file %s\n",
247 fileTest.open(fileName.c_str());
248 if (!fileTest.is_open())
250 runtime_error(
strpr(
"ERROR: Could not open file %s\n",
253 for (vector<Structure>::iterator it =
structures.begin();
261 it->writeToFile(&fileTrain);
265 it->writeToFile(&fileTest);
277 log <<
"Writing training/test set to files:\n";
278 log <<
" - train.data\n";
279 log <<
" - test.data\n";
280 fileName =
"train.data";
282 fileName =
"test.data";
286 log <<
"*****************************************"
287 "**************************************\n";
295 log <<
"*** WEIGHT INITIALIZATION ***************"
296 "**************************************\n";
302 throw runtime_error(
"ERROR: Nguyen Widrow and preconditioning weights"
303 " initialization are incompatible\n");
310 log <<
"Setting up " +
nns.at(
"elec").name +
" neural networks:\n";
316 log <<
"Setting up " + nn.
name +
" neural networks:\n";
319 log <<
"Reading old weights from files.\n";
324 log <<
"*****************************************"
325 "**************************************\n";
336 log <<
strpr(
"Combined updater for all elements selected: "
344 .getNumConnections();
355 log <<
strpr(
"Separate updaters for elements selected: "
363 .getNumConnections();
366 log <<
strpr(
"Fit parameters for element %2s: %zu\n",
373 throw runtime_error(
"ERROR: Unknown update strategy.\n");
387 pk.push_back(
"energy");
399 pk.push_back(
"charge");
404 pk.push_back(
"energy");
407 pk.push_back(
"force");
413 throw runtime_error(
strpr(
"ERROR: No or incorrect training stage "
414 "specified: %zu (must be 1 or 2).\n",
420 auto initP = [
this](
string key) {
p.emplace(piecewise_construct,
421 forward_as_tuple(key),
422 forward_as_tuple(key));
424 for (
auto k :
pk) initP(k);
432 log <<
"*** DATA SET NORMALIZATION **************"
433 "**************************************\n";
436 log <<
"Computing statistics from reference data and initial "
445 throw runtime_error(
"ERROR: Normalization of charges not yet "
456 fileEvsV.open(
strpr(
"evsv.dat.%04d",
myRank).c_str());
460 vector<string> title;
461 vector<string> colName;
462 vector<string> colInfo;
463 vector<size_t> colSize;
464 title.push_back(
"Energy vs. volume comparison.");
465 colSize.push_back(16);
466 colName.push_back(
"V_atom");
467 colInfo.push_back(
"Volume per atom.");
468 colSize.push_back(16);
469 colName.push_back(
"Eref_atom");
470 colInfo.push_back(
"Reference energy per atom.");
471 colSize.push_back(10);
472 colName.push_back(
"N");
473 colInfo.push_back(
"Number of atoms.");
474 colSize.push_back(16);
475 colName.push_back(
"V");
476 colInfo.push_back(
"Volume of structure.");
477 colSize.push_back(16);
478 colName.push_back(
"Eref");
479 colInfo.push_back(
"Reference energy of structure.");
480 colSize.push_back(16);
481 colName.push_back(
"Eref_offset");
482 colInfo.push_back(
"Reference energy of structure (including offset).");
487 size_t numAtomsTotal = 0;
489 double meanEnergyPerAtomRef = 0.0;
490 double meanEnergyPerAtomNnp = 0.0;
491 double sigmaEnergyPerAtomRef = 0.0;
492 double sigmaEnergyPerAtomNnp = 0.0;
493 double meanForceRef = 0.0;
494 double meanForceNnp = 0.0;
495 double sigmaForceRef = 0.0;
496 double sigmaForceNnp = 0.0;
500 fileEvsV <<
strpr(
"%16.8E %16.8E %10zu %16.8E %16.8E %16.8E\n",
501 s.volume / s.numAtoms,
502 s.energyRef / s.numAtoms,
509#ifdef N2P2_NO_SF_GROUPS
517 s.clearNeighborList();
520 numAtomsTotal += s.numAtoms;
521 meanEnergyPerAtomRef += s.energyRef / s.numAtoms;
522 meanEnergyPerAtomNnp += s.energy / s.numAtoms;
523 for (
auto& a : s.atoms)
525 meanForceRef += a.fRef[0] + a.fRef[1] + a.fRef[2];
526 meanForceNnp += a.f [0] + a.f [1] + a.f [2];
532 MPI_Barrier(MPI_COMM_WORLD);
533 log <<
"Writing energy/atom vs. volume/atom data to \"evsv.dat\".\n";
536 MPI_Allreduce(MPI_IN_PLACE, &numAtomsTotal , 1,
MPI_SIZE_T, MPI_SUM, MPI_COMM_WORLD);
537 MPI_Allreduce(MPI_IN_PLACE, &meanEnergyPerAtomRef, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
538 MPI_Allreduce(MPI_IN_PLACE, &meanEnergyPerAtomNnp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
539 MPI_Allreduce(MPI_IN_PLACE, &meanForceRef , 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
540 MPI_Allreduce(MPI_IN_PLACE, &meanForceNnp , 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
543 meanForceRef /= 3 * numAtomsTotal;
544 meanForceNnp /= 3 * numAtomsTotal;
547 double ediffRef = s.energyRef / s.numAtoms - meanEnergyPerAtomRef;
548 double ediffNnp = s.energy / s.numAtoms - meanEnergyPerAtomNnp;
549 sigmaEnergyPerAtomRef += ediffRef * ediffRef;
550 sigmaEnergyPerAtomNnp += ediffNnp * ediffNnp;
551 for (
auto const& a : s.atoms)
553 double fdiffRef = a.fRef[0] - meanForceRef;
554 double fdiffNnp = a.f [0] - meanForceNnp;
555 sigmaForceRef += fdiffRef * fdiffRef;
556 sigmaForceNnp += fdiffNnp * fdiffNnp;
557 fdiffRef = a.fRef[1] - meanForceRef;
558 fdiffNnp = a.f [1] - meanForceNnp;
559 sigmaForceRef += fdiffRef * fdiffRef;
560 sigmaForceNnp += fdiffNnp * fdiffNnp;
561 fdiffRef = a.fRef[2] - meanForceRef;
562 fdiffNnp = a.f [2] - meanForceNnp;
563 sigmaForceRef += fdiffRef * fdiffRef;
564 sigmaForceNnp += fdiffNnp * fdiffNnp;
567 MPI_Allreduce(MPI_IN_PLACE, &sigmaEnergyPerAtomRef, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
568 MPI_Allreduce(MPI_IN_PLACE, &sigmaEnergyPerAtomNnp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
569 MPI_Allreduce(MPI_IN_PLACE, &sigmaForceRef , 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
570 MPI_Allreduce(MPI_IN_PLACE, &sigmaForceNnp , 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
571 sigmaEnergyPerAtomRef = sqrt(sigmaEnergyPerAtomRef / (
numStructures - 1));
572 sigmaEnergyPerAtomNnp = sqrt(sigmaEnergyPerAtomNnp / (
numStructures - 1));
573 sigmaForceRef = sqrt(sigmaForceRef / (3 * numAtomsTotal - 1));
574 sigmaForceNnp = sqrt(sigmaForceNnp / (3 * numAtomsTotal - 1));
577 log <<
strpr(
"Total number of atoms : %zu\n", numAtomsTotal);
578 log <<
"----------------------------------\n";
579 log <<
"Reference data statistics:\n";
580 log <<
"----------------------------------\n";
581 log <<
strpr(
"Mean/sigma energy per atom : %16.8E / %16.8E\n",
582 meanEnergyPerAtomRef,
583 sigmaEnergyPerAtomRef);
584 log <<
strpr(
"Mean/sigma force : %16.8E / %16.8E\n",
587 log <<
"----------------------------------\n";
588 log <<
"Initial NNP prediction statistics:\n";
589 log <<
"----------------------------------\n";
590 log <<
strpr(
"Mean/sigma energy per atom : %16.8E / %16.8E\n",
591 meanEnergyPerAtomNnp,
592 sigmaEnergyPerAtomNnp);
593 log <<
strpr(
"Mean/sigma force : %16.8E / %16.8E\n",
596 log <<
"----------------------------------\n";
598 if (
settings[
"normalize_data_set"] ==
"stats-only")
600 log <<
"Data set statistics computation completed, now make up for \n";
601 log <<
"initially skipped normalization setup...\n";
605 else if (
settings[
"normalize_data_set"] ==
"ref")
607 log <<
"Normalization based on standard deviation of reference data "
610 log <<
" mean(e_ref) = 0, sigma(e_ref) = sigma(F_ref) = 1\n";
614 if (useForcesLocal)
convLength = sigmaForceRef / sigmaEnergyPerAtomRef;
618 else if (
settings[
"normalize_data_set"] ==
"force")
622 throw runtime_error(
"ERROR: Selected normalization mode only "
623 "possible when forces are available.\n");
625 log <<
"Normalization based on standard deviation of reference forces "
627 log <<
"initial prediction selected:\n";
629 log <<
" mean(e_ref) = 0, sigma(F_NNP) = sigma(F_ref) = 1\n";
638 throw runtime_error(
"ERROR: Unknown data set normalization mode.\n");
641 if (
settings[
"normalize_data_set"] !=
"stats-only")
643 log <<
"Final conversion data:\n";
647 log <<
"----------------------------------\n";
651 (
settings[
"normalize_data_set"] !=
"stats-only"))
654 log <<
"Writing backup of original settings file to "
655 "\"input.nn.bak\".\n";
656 ofstream fileSettings;
657 fileSettings.open(
"input.nn.bak");
659 fileSettings.close();
662 log <<
"Writing normalization data to settings file \"input.nn\".\n";
663 string n1 =
strpr(
"mean_energy %24.16E # nnp-train\n",
664 meanEnergyPerAtomRef);
665 string n2 =
strpr(
"conv_energy %24.16E # nnp-train\n",
667 string n3 =
strpr(
"conv_length %24.16E # nnp-train\n",
672 map<size_t, string> replace;
673 for (
size_t i = 0; i < lines.size(); ++i)
675 vector<string> sl =
split(lines.at(i));
678 if (sl.at(0) ==
"mean_energy") replace[i] = n1;
679 if (sl.at(0) ==
"conv_energy") replace[i] = n2;
680 if (sl.at(0) ==
"conv_length") replace[i] = n3;
683 if (!replace.empty())
685 log <<
"WARNING: Preexisting normalization data was found and "
686 "replaced in original \"input.nn\" file.\n";
689 fileSettings.open(
"input.nn");
692 fileSettings <<
"#########################################"
693 "######################################\n";
694 fileSettings <<
"# DATA SET NORMALIZATION\n";
695 fileSettings <<
"#########################################"
696 "######################################\n";
700 fileSettings <<
"#########################################"
701 "######################################\n";
702 fileSettings <<
"\n";
705 fileSettings.close();
713 log <<
"Silently repeating symmetry function setup...\n";
715 for (
auto& e :
elements) e.clearSymmetryFunctions();
717#ifndef N2P2_FULL_SFD_MEMORY
720#ifndef N2P2_NO_SF_CACHE
723#ifndef N2P2_NO_SF_GROUPS
731 log <<
"*****************************************"
732 "**************************************\n";
740 log <<
"*** SETUP: TRAINING *********************"
741 "**************************************\n";
748 if (
stage == 1)
log <<
"electrostatic NN fitting.\n";
749 else if (
stage == 2)
log <<
"short-range NN fitting.\n";
750 else throw runtime_error(
"\nERROR: Unknown training stage.\n");
760 log <<
"Forces will be used for training.\n";
767 log <<
"WARNING: Force weight not set, using default value.\n";
774 log <<
"Only energies will be used for training.\n";
777 log <<
"Training will act on \"" <<
nns.at(
nnId).name
778 <<
"\" neural networks.\n";
783 if (
settings[
"main_error_metric"] ==
"RMSEpa")
785 k =
"energy";
if (
p.
exists(k))
p[k].displayMetric =
"RMSEpa";
786 k =
"force";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
787 k =
"charge";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
789 else if (
settings[
"main_error_metric"] ==
"RMSE")
791 k =
"energy";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
792 k =
"force";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
793 k =
"charge";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
795 else if (
settings[
"main_error_metric"] ==
"MAEpa")
797 k =
"energy";
if (
p.
exists(k))
p[k].displayMetric =
"MAEpa";
798 k =
"force";
if (
p.
exists(k))
p[k].displayMetric =
"MAE";
799 k =
"charge";
if (
p.
exists(k))
p[k].displayMetric =
"MAE";
801 else if (
settings[
"main_error_metric"] ==
"MAE")
803 k =
"energy";
if (
p.
exists(k))
p[k].displayMetric =
"MAE";
804 k =
"force";
if (
p.
exists(k))
p[k].displayMetric =
"MAE";
805 k =
"charge";
if (
p.
exists(k))
p[k].displayMetric =
"MAE";
809 throw runtime_error(
"ERROR: Unknown error metric.\n");
815 k =
"energy";
if (
p.
exists(k))
p[k].displayMetric =
"RMSEpa";
816 k =
"force";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
817 k =
"charge";
if (
p.
exists(k))
p[k].displayMetric =
"RMSE";
823 log <<
strpr(
"Weight update via gradient descent selected: "
824 "updaterType::UT_GD (%d)\n",
829 log <<
strpr(
"Weight update via Kalman filter selected: "
830 "updaterType::UT_KF (%d)\n",
835 throw runtime_error(
"ERROR: LM algorithm not yet implemented.\n");
836 log <<
strpr(
"Weight update via Levenberg-Marquardt algorithm "
837 "selected: updaterType::UT_LM (%d)\n",
842 throw runtime_error(
"ERROR: Unknown updater type.\n");
854 log <<
strpr(
"Parallel training (rank 0 updates) selected: "
855 "ParallelMode::PM_TRAIN_RK0 (%d)\n",
860 log <<
strpr(
"Parallel training (all ranks update) selected: "
861 "ParallelMode::PM_TRAIN_ALL (%d)\n",
866 throw runtime_error(
"ERROR: Unknown parallelization mode.\n");
872 log <<
strpr(
"Gradient summation only selected: "
874 log <<
"No Jacobi matrix, gradients of all training candidates are "
875 "summed up instead.\n";
879 log <<
strpr(
"Per-task Jacobian selected: "
880 "JacobianMode::JM_TASK (%d)\n",
882 log <<
"One Jacobi matrix row per MPI task is stored, within each "
883 "task gradients are summed up.\n";
887 log <<
strpr(
"Full Jacobian selected: "
888 "JacobianMode::JM_FULL (%d)\n",
890 log <<
"Each update candidate generates one Jacobi matrix "
895 throw runtime_error(
"ERROR: Unknown Jacobian mode.\n");
900 throw runtime_error(
"ERROR: Gradient descent methods can only be "
901 "combined with Jacobian mode JM_SUM.\n");
906 throw runtime_error(
"ERROR: US_ELEMENT only implemented for "
922 log <<
"-----------------------------------------"
923 "--------------------------------------\n";
927 throw runtime_error(
"ERROR: Repeated energy updates are not correctly"
928 " implemented at the moment.\n");
936 log <<
"Symmetry function memory is cleared after each calculation.\n";
940 log <<
"Symmetry function memory is reused (HIGH MEMORY USAGE!).\n";
962 log <<
strpr(
"Training log with update information will be written to:"
967 vector<string> title;
968 vector<string> colName;
969 vector<string> colInfo;
970 vector<size_t> colSize;
971 title.push_back(
"Detailed information on each weight update.");
972 colSize.push_back(3);
973 colName.push_back(
"U");
974 colInfo.push_back(
"Update type (E = energy, F = force, Q = charge).");
975 colSize.push_back(5);
976 colName.push_back(
"epoch");
977 colInfo.push_back(
"Current training epoch.");
978 colSize.push_back(10);
979 colName.push_back(
"count");
980 colInfo.push_back(
"Update counter (Multiple lines with identical count"
981 " for multi-streaming!).");
982 colSize.push_back(5);
983 colName.push_back(
"proc");
984 colInfo.push_back(
"MPI process providing this update candidate.");
985 colSize.push_back(3);
986 colName.push_back(
"tl");
987 colInfo.push_back(
"Threshold loop counter.");
988 colSize.push_back(10);
989 colName.push_back(
"rmse_frac");
990 colInfo.push_back(
"Update candidates error divided by this "
992 colSize.push_back(10);
993 colName.push_back(
"s_ind_g");
994 colInfo.push_back(
"Global structure index.");
995 colSize.push_back(5);
996 colName.push_back(
"s_ind");
997 colInfo.push_back(
"Local structure index on this MPI process.");
998 colSize.push_back(5);
999 colName.push_back(
"a_ind");
1000 colInfo.push_back(
"Atom index.");
1001 colSize.push_back(2);
1002 colName.push_back(
"c");
1003 colInfo.push_back(
"Force component (0 = x, 1 = y, 2 = z).");
1009 log <<
"-----------------------------------------"
1010 "--------------------------------------\n";
1016 log <<
strpr(
"Energy to force ratio : "
1018 static_cast<double>(
1021 log <<
strpr(
"Energy to force percentages : "
1022 "%5.1f%% : %5.1f%%\n",
1030 double totalUpdates = 0.0;
1031 for (
auto k :
pk) totalUpdates +=
p[k].numUpdates;
1032 log <<
"-----------------------------------------"
1033 "--------------------------------------\n";
1037 log <<
"-----------------------------------------"
1038 "--------------------------------------\n";
1049 atoi(
settings[
"gradient_type"].c_str());
1055 atoi(
settings[
"kalman_type"].c_str());
1076 updaters.back()->resetTimingLoop();
1086 double const eta = atof(
settings[
"gradient_eta"].c_str());
1096 double const eta = atof(
settings[
"gradient_adam_eta"].c_str());
1097 double const beta1 = atof(
settings[
"gradient_adam_beta1"].c_str());
1098 double const beta2 = atof(
settings[
"gradient_adam_beta2"].c_str());
1099 double const eps = atof(
settings[
"gradient_adam_epsilon"].c_str());
1112 double const epsilon = atof(
settings[
"kalman_epsilon"].c_str());
1113 double const q0 = atof(
settings[
"kalman_q0" ].c_str());
1114 double const qtau = atof(
settings[
"kalman_qtau" ].c_str())
1116 log <<
"qtau is divided by number "
1117 "of projected updates per epoch.\n";
1118 double const qmin = atof(
settings[
"kalman_qmin" ].c_str());
1119 double const eta0 = atof(
settings[
"kalman_eta" ].c_str());
1120 double etatau = 1.0;
1121 double etamax = eta0;
1125 etatau = atof(
settings[
"kalman_etatau"].c_str())
1127 log <<
"etatau is divided by number "
1128 "of projected updates per epoch.\n";
1129 etamax = atof(
settings[
"kalman_etamax"].c_str());
1131 for (
size_t i = 0; i <
updaters.size(); ++i)
1145 double const epsilon = atof(
settings[
"kalman_epsilon"].c_str());
1146 double const q0 = atof(
settings[
"kalman_q0" ].c_str());
1147 double const qtau = atof(
settings[
"kalman_qtau" ].c_str())
1149 log <<
"qtau is divided by number "
1150 "of projected updates per epoch.\n";
1151 double const qmin = atof(
settings[
"kalman_qmin"].c_str());
1152 double const lambda =
1153 atof(
settings[
"kalman_lambda_short"].c_str());
1157 double const nu = atof(
settings[
"kalman_nue_short"].c_str());
1158 for (
size_t i = 0; i <
updaters.size(); ++i)
1171 log <<
"-----------------------------------------"
1172 "--------------------------------------\n";
1173 for (
size_t i = 0; i <
updaters.size(); ++i)
1177 log <<
strpr(
"Combined weight updater:\n");
1181 log <<
strpr(
"Weight updater for element %2s :\n",
1182 elements.at(i).getSymbol().c_str());
1184 log <<
"-----------------------------------------"
1185 "--------------------------------------\n";
1189 log <<
"Note: During training loop the actual observation\n";
1190 log <<
" size corresponds to error vector size:\n";
1193 log <<
strpr(
"sizeObservation = %zu (%s updates)\n",
1194 p[k].error.at(i).size(), k.c_str());
1197 log <<
"-----------------------------------------"
1198 "--------------------------------------\n";
1201 log <<
strpr(
"TIMING Finished setup: %.2f seconds.\n",
1202 sw[
"setup"].stop());
1203 log <<
"*****************************************"
1204 "**************************************\n";
1212 log <<
"*** SETUP WEIGHT DERIVATIVES CHECK ******"
1213 "**************************************\n";
1216 log <<
"Weight derivatives will be checked for these properties:\n";
1217 for (
auto k :
pk)
log <<
" - " +
p[k].plural +
"\n";
1241 log <<
"*****************************************"
1242 "**************************************\n";
1251 log <<
"*** CALCULATE NEIGHBOR LISTS ************"
1252 "**************************************\n";
1256 int num_threads = omp_get_max_threads();
1257 omp_set_num_threads(1);
1258 log <<
strpr(
"Temporarily disabling OpenMP parallelization: %d threads.\n",
1259 omp_get_max_threads());
1261 log <<
"Calculating neighbor lists for all structures.\n";
1267 log <<
strpr(
"Cutoff radius for neighbor lists: %f\n",
1268 maxCutoffRadiusPhys);
1269 double maxCutoffRadiusAllStructures = 0.0;
1270 for (vector<Structure>::iterator it =
structures.begin();
1276 it->calculateMaxCutoffRadiusOverall(
1282 if (it->maxCutoffRadiusOverall > maxCutoffRadiusAllStructures)
1283 maxCutoffRadiusAllStructures = it->maxCutoffRadiusOverall;
1292 log <<
strpr(
"Largest cutoff radius for neighbor lists: %f\n",
1293 maxCutoffRadiusAllStructures);
1295 omp_set_num_threads(num_threads);
1296 log <<
strpr(
"Restoring OpenMP parallelization: max. %d threads.\n",
1297 omp_get_max_threads());
1300 log <<
"-----------------------------------------"
1301 "--------------------------------------\n";
1302 log <<
strpr(
"TIMING Finished neighbor lists: %.2f seconds.\n",
1304 log <<
"*****************************************"
1305 "**************************************\n";
1311 map<
string, pair<string, string>>
const fileNames)
1314 int num_threads = omp_get_max_threads();
1315 omp_set_num_threads(1);
1317 vector<string> write;
1318 for (
auto i : fileNames)
1320 if (i.second.first.size() == 0 || i.second.second.size() == 0)
1322 throw runtime_error(
"ERROR: No filename provided for comparison "
1325 write.push_back(i.first);
1327 auto doWrite = [&write](
string key){
1328 return find(write.begin(),
1330 key) != write.end();
1334 map<string, size_t> countTrain;
1335 map<string, size_t> countTest;
1336 for (
auto k :
pk) countTrain[k] = 0;
1337 for (
auto k :
pk) countTest[k] = 0;
1339 map<string, ofstream> filesTrain;
1340 map<string, ofstream> filesTest;
1345 for (
auto& m :
p[k].errorTrain) m.second = 0.0;
1346 for (
auto& m :
p[k].errorTest) m.second = 0.0;
1349 for (
auto k : write)
1351 filesTrain[k].open(
strpr(
"%s.%04d",
1352 fileNames.at(k).first.c_str(),
1354 filesTest[k].open(
strpr(
"%s.%04d",
1355 fileNames.at(k).second.c_str(),
1358 vector<string> header;
1361 vector<string> title;
1362 vector<string> colName;
1363 vector<string> colInfo;
1364 vector<size_t> colSize;
1367 title.push_back(
"Energy comparison.");
1368 colSize.push_back(10);
1369 colName.push_back(
"index");
1370 colInfo.push_back(
"Structure index.");
1371 colSize.push_back(16);
1372 colName.push_back(
"Eref");
1373 colInfo.push_back(
"Reference potential energy per atom "
1374 "(training units).");
1375 colSize.push_back(16);
1376 colName.push_back(
"Ennp");
1377 colInfo.push_back(
"NNP potential energy per atom "
1378 "(training units).");
1380 else if (k ==
"force")
1382 title.push_back(
"Force comparison.");
1383 colSize.push_back(10);
1384 colName.push_back(
"index_s");
1385 colInfo.push_back(
"Structure index.");
1386 colSize.push_back(10);
1387 colName.push_back(
"index_a");
1388 colInfo.push_back(
"Atom index (x, y, z components in "
1389 "consecutive lines).");
1390 colSize.push_back(16);
1391 colName.push_back(
"Fref");
1392 colInfo.push_back(
"Reference force (training units).");
1393 colSize.push_back(16);
1394 colName.push_back(
"Fnnp");
1395 colInfo.push_back(
"NNP force (training units).");
1397 else if (k ==
"charge")
1399 title.push_back(
"Charge comparison.");
1400 colSize.push_back(10);
1401 colName.push_back(
"index_s");
1402 colInfo.push_back(
"Structure index.");
1403 colSize.push_back(10);
1404 colName.push_back(
"index_a");
1405 colInfo.push_back(
"Atom index.");
1406 colSize.push_back(16);
1407 colName.push_back(
"Qref");
1408 colInfo.push_back(
"Reference charge.");
1409 colSize.push_back(16);
1410 colName.push_back(
"Qnnp");
1411 colInfo.push_back(
"NNP charge.");
1419 for (vector<Structure>::iterator it =
structures.begin();
1422#ifdef N2P2_NO_SF_GROUPS
1437 if ( !it->hasCharges || (!it->hasAMatrix &&
useForces) )
1457 map<string, double>* error =
nullptr;
1458 size_t* count =
nullptr;
1459 ofstream* file =
nullptr;
1462 error = &(
p[k].errorTrain);
1463 count = &(countTrain.at(k));
1464 if (doWrite(k)) file = &(filesTrain.at(k));
1468 error = &(
p[k].errorTest);
1469 count = &(countTest.at(k));
1470 if (doWrite(k)) file = &(filesTest.at(k));
1473 it->updateError(k, *error, *count);
1476 if (k ==
"energy") (*file) << it->getEnergyLine();
1477 else if (k ==
"force")
1479 for (
auto l : it->getForcesLines()) (*file) << l;
1481 else if (k ==
"charge")
1483 for (
auto l : it->getChargesLines()) (*file) << l;
1488 it->clearElectrostatics();
1497 filesTrain.at(k).close();
1498 filesTest.at(k).close();
1509 omp_set_num_threads(num_threads);
1519 map<string, pair<string, string>> fileNames;
1521 for (
auto const& ip :
p)
1523 string const& k = ip.first;
1525 if (
d.writeCompEvery > 0 &&
1526 (
epoch %
d.writeCompEvery == 0 ||
epoch <=
d.writeCompAlways))
1529 if (k ==
"energy") middle =
"points";
1530 else if (k ==
"force" ) middle =
"forces";
1531 else if (k ==
"charge") middle =
"charges";
1532 fileNames[k] = make_pair(
strpr(
"train%s.%06zu.out",
1533 middle.c_str(),
epoch),
1534 strpr(
"test%s.%06zu.out",
1535 middle.c_str(),
epoch));
1546 Eigen::VectorXd &cVec,
1550 for (
size_t i = 0; i < s.
numAtoms; ++i)
1552 cVec(i) = s.
atoms.at(i).charge - s.
atoms.at(i).chargeRef;
1554 cNorm = cVec.norm();
1560 string metric =
"?";
1561 string peratom =
"";
1563 log <<
"The training loop output covers different errors, update and\n";
1564 log <<
"timing information. The following quantities are organized\n";
1565 log <<
"according to the matrix scheme below:\n";
1566 log <<
"-------------------------------------------------------------------\n";
1567 log <<
"ep ........ Epoch.\n";
1570 string const& pmetric =
p[k].displayMetric;
1571 if (pmetric.find(
"RMSE") != pmetric.npos) metric =
"RMSE";
1572 else if (pmetric.find(
"MAE") != pmetric.npos) metric =
"MAE";
1573 if (pmetric.find(
"pa") != pmetric.npos) peratom =
" per atom";
1575 log <<
p[k].tiny <<
"_count ... Number of " << k <<
" updates.\n";
1576 log <<
p[k].tiny <<
"_train ... " << metric <<
" of training "
1577 <<
p[k].plural << peratom <<
".\n";
1578 log <<
p[k].tiny <<
"_test .... " << metric <<
" of test "
1579 <<
p[k].plural << peratom <<
".\n";
1582 log <<
p[k].tiny <<
"_pt ...... Percentage of time for " << k <<
1583 " updates w.r.t. to t_train.\n";
1585 log <<
"count ..... Total number of updates.\n";
1586 log <<
"train ..... Percentage of time for training.\n";
1587 log <<
"error ..... Percentage of time for error calculation.\n";
1588 log <<
"other ..... Percentage of time for other purposes.\n";
1589 log <<
"epoch ..... Total time for this epoch (seconds).\n";
1590 log <<
"total ..... Total time for all epochs (seconds).\n";
1591 log <<
"-------------------------------------------------------------------\n";
1595 <<
strpr(
" %5s",
"ep")
1596 <<
strpr(
" %7s", (
p[k].tiny +
"_count").c_str())
1597 <<
strpr(
" %11s", (
p[k].tiny +
"_train").c_str())
1598 <<
strpr(
" %11s", (
p[k].tiny +
"_test").c_str())
1599 <<
strpr(
" %5s", (
p[k].tiny +
"_pt").c_str())
1603 <<
strpr(
" %5s",
"ep")
1604 <<
strpr(
" %7s",
"count")
1605 <<
strpr(
" %5s",
"train")
1606 <<
strpr(
" %5s",
"error")
1607 <<
strpr(
" %5s",
"other")
1608 <<
strpr(
" %9s",
"epoch")
1609 <<
strpr(
" %9s",
"total")
1611 log <<
"-------------------------------------------------------------------\n";
1618 double timeLoop =
sw[
"loop"].getLoop();
1619 double timeTrain =
sw[
"train"].getLoop();
1620 size_t totalUpdates = 0;
1623 totalUpdates +=
p[k].countUpdates;
1624 double timeProp =
sw[k].getLoop();
1626 for (
auto& c : caps) c = toupper(c);
1627 log <<
strpr(
"%-6s", caps.c_str());
1633 physical(k,
p[k].errorTrain.at(
p[k].displayMetric)),
1634 physical(k,
p[k].errorTest.at(
p[k].displayMetric)));
1639 p[k].errorTrain.at(
p[k].displayMetric),
1640 p[k].errorTest.at(
p[k].displayMetric));
1643 else log <<
strpr(
" %5.1f", timeProp / timeTrain * 100.0);
1646 double timeOther = timeLoop;
1647 timeOther -=
sw[
"error"].getLoop();
1648 timeOther -=
sw[
"train"].getLoop();
1651 log <<
strpr(
" %7zu", totalUpdates);
1652 log <<
strpr(
" %5.1f",
sw[
"train"].getLoop() / timeLoop * 100.0);
1653 log <<
strpr(
" %5.1f",
sw[
"error"].getLoop() / timeLoop * 100.0);
1654 log <<
strpr(
" %5.1f", timeOther / timeLoop * 100.0);
1655 log <<
strpr(
" %9.2f",
sw[
"loop"].getLoop());
1656 log <<
strpr(
" %9.2f",
sw[
"loop"].getTotal());
1663 string const& fileNameFormat)
const
1669 string fileName =
strpr(fileNameFormat.c_str(),
1671 file.open(fileName.c_str());
1672 elements.at(i).neuralNetworks.at(nnName).writeConnections(file);
1684 string weightFileFormat =
strpr(
".%%03zu.%06d.out",
epoch);
1689 weightFileFormat =
"weights" + weightFileFormat;
1694 weightFileFormat =
"weightse" + weightFileFormat;
1708 string fileName =
strpr(fileNameFormat.c_str(),
1710 file.open(fileName.c_str());
1711 double hardness =
elements.at(i).getHardness();
1713 file << hardness << endl;
1726 string hardnessFileFormat =
strpr(
".%%03zu.%06d.out",
epoch);
1727 hardnessFileFormat =
"hardness" + hardnessFileFormat;
1737 string fileNameActual = fileName;
1741 fileNameActual +=
strpr(
".stage-%zu",
stage);
1744 if (append) file.open(fileNameActual.c_str(), ofstream::app);
1747 file.open(fileNameActual.c_str());
1750 vector<string> title;
1751 vector<string> colName;
1752 vector<string> colInfo;
1753 vector<size_t> colSize;
1754 title.push_back(
"Learning curves for training properties:");
1757 title.push_back(
" * " +
p[k].plural);
1759 colSize.push_back(10);
1760 colName.push_back(
"epoch");
1761 colInfo.push_back(
"Current epoch.");
1763 map<string, string> text;
1764 text[
"RMSEpa"] =
"RMSE of %s %s per atom";
1765 text[
"RMSE"] =
"RMSE of %s %s";
1766 text[
"MAEpa"] =
"MAE of %s %s per atom";
1767 text[
"MAE"] =
"MAE of %s %s";
1771 for (
auto m :
p[k].errorMetrics)
1773 colSize.push_back(16);
1774 colName.push_back(m +
"_" +
p[k].tiny +
"train_pu");
1775 colInfo.push_back(
strpr(
1776 (text[m] +
" (physical units)").c_str(),
1778 p[k].plural.c_str()));
1779 colSize.push_back(16);
1780 colName.push_back(m +
"_" +
p[k].tiny +
"test_pu");
1781 colInfo.push_back(
strpr(
1782 (text[m] +
" (physical units)").c_str(),
1784 p[k].plural.c_str()));
1792 if (!(k ==
"energy" || k ==
"force" || k ==
"charge"))
continue;
1793 for (
auto m :
p[k].errorMetrics)
1795 colSize.push_back(16);
1796 colName.push_back(m +
"_" +
p[k].tiny +
"train_iu");
1797 colInfo.push_back(
strpr(
1798 (text[m] +
" (training units)").c_str(),
1800 p[k].plural.c_str()));
1801 colSize.push_back(16);
1802 colName.push_back(m +
"_" +
p[k].tiny +
"test_iu");
1803 colInfo.push_back(
strpr(
1804 (text[m] +
" (training units)").c_str(),
1806 p[k].plural.c_str()));
1819 if (!(k ==
"energy" || k ==
"force" || k ==
"charge"))
continue;
1820 for (
auto m :
p[k].errorMetrics)
1822 file <<
strpr(
" %16.8E %16.8E",
1830 for (
auto m :
p[k].errorMetrics)
1832 file <<
strpr(
" %16.8E %16.8E",
1833 p[k].errorTrain.at(m),
1834 p[k].errorTest.at(m));
1845 string const& fileName)
const
1850 file.open(fileName.c_str());
1853 vector<string> title;
1854 vector<string> colName;
1855 vector<string> colInfo;
1856 vector<size_t> colSize;
1857 title.push_back(
"Statistics for individual neurons of network \""
1858 +
id +
"\" gathered during RMSE calculation.");
1859 colSize.push_back(10);
1860 colName.push_back(
"element");
1861 colInfo.push_back(
"Element index.");
1862 colSize.push_back(10);
1863 colName.push_back(
"neuron");
1864 colInfo.push_back(
"Neuron number.");
1865 colSize.push_back(10);
1866 colName.push_back(
"count");
1867 colInfo.push_back(
"Number of neuron value computations.");
1868 colSize.push_back(16);
1869 colName.push_back(
"min");
1870 colInfo.push_back(
"Minimum neuron value encounterd.");
1871 colSize.push_back(16);
1872 colName.push_back(
"max");
1873 colInfo.push_back(
"Maximum neuron value encounterd.");
1874 colSize.push_back(16);
1875 colName.push_back(
"mean");
1876 colInfo.push_back(
"Mean neuron value.");
1877 colSize.push_back(16);
1878 colName.push_back(
"sigma");
1879 colInfo.push_back(
"Standard deviation of neuron value.");
1886 size_t n =
elements.at(i).neuralNetworks.at(
id).getNumNeurons();
1887 vector<long> count(n, 0);
1888 vector<double> min(n, 0.0);
1889 vector<double> max(n, 0.0);
1890 vector<double> mean(n, 0.0);
1891 vector<double> sigma(n, 0.0);
1892 elements.at(i).neuralNetworks.at(
id).
1893 getNeuronStatistics(&(count.front()),
1901 MPI_Reduce(MPI_IN_PLACE, &(count.front()), n, MPI_LONG , MPI_SUM, 0,
comm);
1902 MPI_Reduce(MPI_IN_PLACE, &(min.front()) , n, MPI_DOUBLE, MPI_MIN, 0,
comm);
1903 MPI_Reduce(MPI_IN_PLACE, &(max.front()) , n, MPI_DOUBLE, MPI_MAX, 0,
comm);
1904 MPI_Reduce(MPI_IN_PLACE, &(mean.front()) , n, MPI_DOUBLE, MPI_SUM, 0,
comm);
1905 MPI_Reduce(MPI_IN_PLACE, &(sigma.front()), n, MPI_DOUBLE, MPI_SUM, 0,
comm);
1909 MPI_Reduce(&(count.front()), &(count.front()), n, MPI_LONG , MPI_SUM, 0,
comm);
1910 MPI_Reduce(&(min.front()) , &(min.front()) , n, MPI_DOUBLE, MPI_MIN, 0,
comm);
1911 MPI_Reduce(&(max.front()) , &(max.front()) , n, MPI_DOUBLE, MPI_MAX, 0,
comm);
1912 MPI_Reduce(&(mean.front()) , &(mean.front()) , n, MPI_DOUBLE, MPI_SUM, 0,
comm);
1913 MPI_Reduce(&(sigma.front()), &(sigma.front()), n, MPI_DOUBLE, MPI_SUM, 0,
comm);
1917 for (
size_t j = 0; j < n; ++j)
1919 size_t m = count.at(j);
1920 sigma.at(j) = sqrt((m * sigma.at(j) - mean.at(j) * mean.at(j))
1923 file <<
strpr(
"%10d %10d %10d %16.8E %16.8E %16.8E %16.8E\n",
1949 string fileName =
strpr(
"neuron-stats.%s.%06zu.out",
1964 for (vector<Element>::iterator it =
elements.begin();
1967 for (
auto& nn : it->neuralNetworks) nn.second.resetNeuronStatistics();
1973 string const fileNameFormat)
const
1976 string fileNameFormatActual = fileNameFormat;
1980 fileNameFormatActual +=
strpr(
".stage-%zu",
stage);
1988 fileName =
strpr(fileNameFormatActual.c_str(), 0);
1992 fileName =
strpr(fileNameFormatActual.c_str(),
1995 if (append) file.open(fileName.c_str(), ofstream::app);
1998 file.open(fileName.c_str());
2011 for (
auto& uc :
p[property].updateCandidates)
2013 if (property ==
"energy")
2018 else if (property ==
"force")
2022 for (
auto &sc : uc.subCandidates)
2025 sc.error = fabs(ai.
fRef[sc.c] - ai.
f[sc.c]);
2026 uc.error += sc.error;
2028 uc.error /= uc.subCandidates.size();
2030 else if (property ==
"charge")
2034 for (
auto const& a : s.
atoms)
2036 uc.error = fabs(a.chargeRef - a.charge);
2042 sort(
p[property].updateCandidates.begin(),
2043 p[property].updateCandidates.end());
2045 for (
auto &uc :
p[property].updateCandidates)
2047 if (uc.subCandidates.size() > 0)
2048 sort(uc.subCandidates.begin(),
2049 uc.subCandidates.end());
2053 p[property].posUpdateCandidates = 0;
2054 for (
auto &uc :
p[property].updateCandidates)
2056 uc.posSubCandidates = 0;
2064 shuffle(
p[property].updateCandidates.begin(),
2065 p[property].updateCandidates.end(),
2068 for (
auto &uc :
p[property].updateCandidates)
2070 if (uc.subCandidates.size() > 0)
2071 shuffle(uc.subCandidates.begin(),
2072 uc.subCandidates.end(),
2077 p[property].posUpdateCandidates = 0;
2078 for (
auto &uc :
p[property].updateCandidates)
2080 uc.posSubCandidates = 0;
2095 for (
int i =
pk.size() - 1; i >= 0; --i)
2101 if (
pk.size() == 1)
return;
2118 if (
p[k].selectionModeSchedule.find(
epoch)
2119 !=
p[k].selectionModeSchedule.end())
2121 p[k].selectionMode =
p[k].selectionModeSchedule[
epoch];
2124 string message =
"INFO Switching selection mode for "
2125 "property \"" + k +
"\" to ";
2128 message +=
strpr(
"SM_RANDOM (%d).\n",
p[k].selectionMode);
2130 else if (
p[k].selectionMode ==
SM_SORT)
2132 message +=
strpr(
"SM_SORT (%d).\n",
p[k].selectionMode);
2136 message +=
strpr(
"SM_THRESHOLD (%d).\n",
2137 p[k].selectionMode);
2151 log <<
"*** TRAINING LOOP ***********************"
2152 "**************************************\n";
2157 sw[
"error"].start();
2191 for (
auto k :
pk)
p[k].countUpdates = 0;
2207 sw[
"train"].start();
2210 string property =
pk.at(i);
2212 p[property].countUpdates++;
2220 sw[
"error"].start();
2246 log <<
"-----------------------------------------"
2247 "--------------------------------------\n";
2248 log <<
strpr(
"TIMING Training loop finished: %.2f seconds.\n",
2249 sw[
"loop"].getTotal());
2250 log <<
"*****************************************"
2251 "**************************************\n";
2259 string const& k = property;
2264 sw[k].start(newLoop);
2265 sw[k +
"_err"].start(newLoop);
2268 int num_threads = omp_get_max_threads();
2269 omp_set_num_threads(1);
2278 bool derivatives =
false;
2279 if (k ==
"force") derivatives =
true;
2281 vector<size_t> thresholdLoopCount(batchSize, 0);
2282 vector<double> currentRmseFraction(batchSize, 0.0);
2287 vector<size_t> currentUpdateCandidates(batchSize, 0);
2291 fill(pu.
error.at(i).begin(), pu.
error.at(i).end(), 0.0);
2296 for (
size_t b = 0; b < batchSize; ++b)
2300 size_t* posCandidates = NULL;
2301 size_t indexBest = 0;
2302 double rmseFractionBest = 0.0;
2309 for (il = 0; il < trials; ++il)
2323 if (useSubCandidates)
2336 currentUpdateCandidates.at(b) = *posCandidates;
2356 currentRmseFraction.at(b) =
2378 throw runtime_error(
"ERROR: Not implemented.\n");
2381 else if (k ==
"force")
2388 currentRmseFraction.at(b) =
2404 currentRmseFraction.at(b) =
2413 throw runtime_error(
"ERROR: Not implemented.\n");
2416 else if (k ==
"charge")
2423 Eigen::VectorXd QError;
2426 currentRmseFraction.at(b) =
2439 currentRmseFraction.at(b) =
2459 if (currentRmseFraction.at(b) > rmseFractionBest)
2461 rmseFractionBest = currentRmseFraction.at(b);
2462 indexBest = *posCandidates;
2476 thresholdLoopCount.at(b) = il;
2483 currentUpdateCandidates.at(b) = indexBest;
2484 currentRmseFraction.at(b) = rmseFractionBest;
2488#ifdef N2P2_NO_SF_GROUPS
2497 if (useSubCandidates)
2504 MPI_INT, MPI_MIN,
comm);
2518 vector<vector<double>> dXdc;
2521 Eigen::VectorXd QError;
2523 double QErrorNorm = 0;
2527 .getNumConnections();
2528 dXdc.at(i).resize(n, 0.0);
2571 for (vector<Atom>::iterator it = s.
atoms.begin();
2572 it != s.
atoms.end(); ++it)
2574 size_t i = it->element;
2578 for (
size_t j = 0; j < it->G.size(); ++j)
2583 nn.
setInput(it->G.size(), it->charge);
2592 for (
size_t j = 0; j < dXdc.at(i).size(); ++j)
2594 pu.
jacobian.at(iu).at(offset.at(i) + j) +=
2603 throw runtime_error(
"ERROR: Not implemented.\n");
2606 else if (k ==
"force")
2618 vector<size_t>{sC->
c},
2624 for (vector<Atom>::iterator it = s.
atoms.begin();
2625 it != s.
atoms.end(); ++it)
2629#ifndef N2P2_FULL_SFD_MEMORY
2634 double dQdxia = s.
atoms.at(sC->
a).dQdr.at(it->index)[sC->
c];
2641 double dQdxia = s.
atoms.at(sC->
a).dQdr.at(it->index)[sC->
c];
2642 it->dGdxia.resize(it->G.size() + 1);
2643 it->dGdxia.back() = dQdxia;
2646 size_t i = it->element;
2649 for (
size_t j = 0; j < it->G.size(); ++j)
2654 nn.
setInput(it->G.size(), it->charge);
2661#ifndef N2P2_FULL_SFD_MEMORY
2666 &(it->dGdxia.front()));
2671 for (
size_t j = 0; j < dXdc.at(i).size(); ++j)
2673 pu.
jacobian.at(iu).at(offset.at(i) + j) +=
2683 throw runtime_error(
"ERROR: Not implemented.\n");
2686 else if (k ==
"charge")
2693 vector<vector<double>> dChidc;
2695 for (
size_t k = 0; k < s.
numAtoms; ++k)
2699 .getNumConnections();
2700 dChidc.at(k).resize(n, 0.0);
2713 for (
auto& dChidGi : ak.
dChidG)
2715 for (
auto& dChidci : dChidc.at(k))
2722 vector<Eigen::VectorXd> dQdChi;
2724 vector<Eigen::VectorXd> dQdJ;
2729 if (QErrorNorm != 0)
2732 for (
size_t i = 0; i < s.
numAtoms; ++i)
2735 for (
size_t k = 0; k < s.
numAtoms; ++k)
2737 size_t l = s.
atoms.at(k).element;
2738 for (
size_t j = 0; j < dChidc.at(k).size(); ++j)
2741 pu.
jacobian.at(0).at(offset.at(l) + j) +=
2742 1.0 / QErrorNorm * QError(i)
2743 * dQdChi.at(k)(i) * dChidc.at(k).at(j);
2750 .getNumConnections();
2751 pu.
jacobian.at(0).at(offset.at(k) + n) +=
2753 * QError(i) * dQdJ.at(k)(i) * 2
2754 * sqrt(
elements.at(k).getHardness());
2775 for (
size_t j = 0; j < dXdc.at(i).size(); ++j)
2777 pu.
jacobian.at(iu).at(offset.at(i) + j) +=
2791 else if (k ==
"force")
2795 currentRmseFraction.at(b) = fabs(a.
fRef[sC->
c] - a.
f[sC->
c])
2798 else if (k ==
"charge")
2802 currentRmseFraction.at(b) = QErrorNorm / sqrt(s.
numAtoms)
2841 else if (k ==
"force")
2844 pu.
error.at(0).at(offset2) += a.
fRef[sC->
c] - a.
f[sC->
c];
2846 else if (k ==
"charge")
2849 pu.
error.at(0).at(offset2) = -QErrorNorm;
2867 else if (k ==
"force")
2870 pu.
error.at(i).at(offset2) += (a.
fRef[sC->
c] - a.
f[sC->
c])
2874 else if (k ==
"charge")
2878 throw runtime_error(
"ERROR: US_ELEMENT not implemented"
2879 " for HDNNP_4G.\n");
2893 for (
size_t j = 0; j < pu.
error.at(i).size(); ++j)
2897 for (
size_t j = 0; j < pu.
jacobian.at(i).size(); ++j)
2903 sw[k +
"_err"].stop();
2909 sw[k +
"_com"].start(newLoop);
2916 if (
myRank == 0) MPI_Reduce(MPI_IN_PLACE , &(pu.
error.at(i).front()), 1, MPI_DOUBLE, MPI_SUM, 0,
comm);
2917 else MPI_Reduce(&(pu.
error.at(i).front()), &(pu.
error.at(i).front()), 1, MPI_DOUBLE, MPI_SUM, 0,
comm);
2926 MPI_Allreduce(MPI_IN_PLACE, &(pu.
error.at(i).front()), 1, MPI_DOUBLE, MPI_SUM,
comm);
2937 if (
myRank == 0) MPI_Gather(MPI_IN_PLACE , 1, MPI_DOUBLE, &(pu.
error.at(i).front()), 1, MPI_DOUBLE, 0,
comm);
2938 else MPI_Gather(&(pu.
error.at(i).front()), 1, MPI_DOUBLE, NULL , 1, MPI_DOUBLE, 0,
comm);
2947 MPI_Allgather(MPI_IN_PLACE, 1, MPI_DOUBLE, &(pu.
error.at(i).front()), 1, MPI_DOUBLE,
comm);
2973 sw[k +
"_com"].stop();
2979 sw[k +
"_upd"].start(newLoop);
2981 omp_set_num_threads(num_threads);
2984 for (
size_t i = 0; i <
updaters.size(); ++i)
2987 pu.
error.at(i).size());
2989 pu.
error.at(i).size());
3010 sw[k +
"_upd"].stop();
3016 sw[k +
"_log"].start(newLoop);
3019 vector<int> procUpdateCandidate;
3020 vector<size_t> indexStructure;
3021 vector<size_t> indexStructureGlobal;
3022 vector<size_t> indexAtom;
3023 vector<size_t> indexCoordinate;
3025 vector<int> currentUpdateCandidatesPerTask;
3026 vector<int> currentUpdateCandidatesOffset;
3027 int myCurrentUpdateCandidates = currentUpdateCandidates.size();
3031 currentUpdateCandidatesPerTask.resize(
numProcs, 0);
3032 currentUpdateCandidatesPerTask.at(0) = myCurrentUpdateCandidates;
3034 if (
myRank == 0) MPI_Gather(MPI_IN_PLACE , 1, MPI_INT, &(currentUpdateCandidatesPerTask.front()), 1, MPI_INT, 0,
comm);
3035 else MPI_Gather(&(myCurrentUpdateCandidates), 1, MPI_INT, NULL , 1, MPI_INT, 0,
comm);
3039 int totalUpdateCandidates = 0;
3040 for (
size_t i = 0; i < currentUpdateCandidatesPerTask.size(); ++i)
3042 currentUpdateCandidatesOffset.push_back(totalUpdateCandidates);
3043 totalUpdateCandidates += currentUpdateCandidatesPerTask.at(i);
3045 procUpdateCandidate.resize(totalUpdateCandidates, 0);
3046 indexStructure.resize(totalUpdateCandidates, 0);
3047 indexStructureGlobal.resize(totalUpdateCandidates, 0);
3048 indexAtom.resize(totalUpdateCandidates, 0);
3049 indexCoordinate.resize(totalUpdateCandidates, 0);
3051 currentRmseFraction.resize(totalUpdateCandidates, 0.0);
3052 thresholdLoopCount.resize(totalUpdateCandidates, 0.0);
3056 procUpdateCandidate.resize(myCurrentUpdateCandidates, 0);
3057 indexStructure.resize(myCurrentUpdateCandidates, 0);
3058 indexStructureGlobal.resize(myCurrentUpdateCandidates, 0);
3059 indexAtom.resize(myCurrentUpdateCandidates, 0);
3060 indexCoordinate.resize(myCurrentUpdateCandidates, 0);
3062 for (
int i = 0; i < myCurrentUpdateCandidates; ++i)
3064 procUpdateCandidate.at(i) =
myRank;
3067 if (useSubCandidates)
3078 indexStructure.at(i) = c->
s;
3079 indexStructureGlobal.at(i) =
structures.at(c->
s).index;
3080 if (useSubCandidates)
3082 indexAtom.at(i) = sC->
a;
3083 indexCoordinate.at(i) = sC->
c;
3088 MPI_Gatherv(MPI_IN_PLACE, 0, MPI_DOUBLE, &(currentRmseFraction.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()), MPI_DOUBLE, 0,
comm);
3089 MPI_Gatherv(MPI_IN_PLACE, 0,
MPI_SIZE_T, &(thresholdLoopCount.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()),
MPI_SIZE_T, 0,
comm);
3090 MPI_Gatherv(MPI_IN_PLACE, 0, MPI_INT , &(procUpdateCandidate.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()), MPI_INT , 0,
comm);
3091 MPI_Gatherv(MPI_IN_PLACE, 0,
MPI_SIZE_T, &(indexStructure.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()),
MPI_SIZE_T, 0,
comm);
3092 MPI_Gatherv(MPI_IN_PLACE, 0,
MPI_SIZE_T, &(indexStructureGlobal.front()), &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()),
MPI_SIZE_T, 0,
comm);
3093 MPI_Gatherv(MPI_IN_PLACE, 0,
MPI_SIZE_T, &(indexAtom.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()),
MPI_SIZE_T, 0,
comm);
3094 MPI_Gatherv(MPI_IN_PLACE, 0,
MPI_SIZE_T, &(indexCoordinate.front()) , &(currentUpdateCandidatesPerTask.front()), &(currentUpdateCandidatesOffset.front()),
MPI_SIZE_T, 0,
comm);
3098 MPI_Gatherv(&(currentRmseFraction.front()) , myCurrentUpdateCandidates, MPI_DOUBLE, NULL, NULL, NULL, MPI_DOUBLE, 0,
comm);
3099 MPI_Gatherv(&(thresholdLoopCount.front()) , myCurrentUpdateCandidates,
MPI_SIZE_T, NULL, NULL, NULL,
MPI_SIZE_T, 0,
comm);
3100 MPI_Gatherv(&(procUpdateCandidate.front()) , myCurrentUpdateCandidates, MPI_INT , NULL, NULL, NULL, MPI_INT , 0,
comm);
3101 MPI_Gatherv(&(indexStructure.front()) , myCurrentUpdateCandidates,
MPI_SIZE_T, NULL, NULL, NULL,
MPI_SIZE_T, 0,
comm);
3102 MPI_Gatherv(&(indexStructureGlobal.front()), myCurrentUpdateCandidates,
MPI_SIZE_T, NULL, NULL, NULL,
MPI_SIZE_T, 0,
comm);
3104 MPI_Gatherv(&(indexCoordinate.front()) , myCurrentUpdateCandidates,
MPI_SIZE_T, NULL, NULL, NULL,
MPI_SIZE_T, 0,
comm);
3109 for (
size_t i = 0; i < procUpdateCandidate.size(); ++i)
3114 thresholdLoopCount.at(i),
3115 currentRmseFraction.at(i),
3116 indexStructureGlobal.at(i),
3117 indexStructure.at(i));
3119 else if (k ==
"force")
3122 thresholdLoopCount.at(i),
3123 currentRmseFraction.at(i),
3124 indexStructureGlobal.at(i),
3125 indexStructure.at(i),
3127 indexCoordinate.at(i));
3129 else if (k ==
"charge")
3132 thresholdLoopCount.at(i),
3133 currentRmseFraction.at(i),
3134 indexStructureGlobal.at(i),
3135 indexStructure.at(i),
3141 sw[k +
"_log"].stop();
3162 return weights.at(element).at(index);
3167 weights.at(element).at(index) = value;
3177#ifdef N2P2_NO_SF_GROUPS
3183 vector<vector<double> > dEdc;
3184 vector<vector<double> > dedc;
3189 size_t n =
elements.at(i).neuralNetworks.at(
"short")
3190 .getNumConnections();
3191 dEdc.at(i).resize(n, 0.0);
3192 dedc.at(i).resize(n, 0.0);
3194 for (vector<Atom>::iterator it = s.
atoms.begin();
3195 it != s.
atoms.end(); ++it)
3197 size_t i = it->element;
3203 for (
size_t j = 0; j < dedc.at(i).size(); ++j)
3205 dEdc.at(i).at(j) += dedc.at(i).at(j);
3216 std::size_t component)
3219#ifdef N2P2_NO_SF_GROUPS
3226 vector<vector<double> > dFdc;
3227 vector<vector<double> > dfdc;
3232 size_t n =
elements.at(i).neuralNetworks.at(
"short")
3233 .getNumConnections();
3234 dFdc.at(i).resize(n, 0.0);
3235 dfdc.at(i).resize(n, 0.0);
3237 for (vector<Atom>::iterator it = s.
atoms.begin();
3238 it != s.
atoms.end(); ++it)
3240#ifndef N2P2_FULL_SFD_MEMORY
3245 size_t i = it->element;
3251#ifndef N2P2_FULL_SFD_MEMORY
3254 nn.
calculateDFdc(&(dfdc.at(i).front()), &(it->dGdxia.front()));
3256 for (
size_t j = 0; j < dfdc.at(i).size(); ++j)
3258 dFdc.at(i).at(j) += dfdc.at(i).at(j);
3277 n += e.neuralNetworks.at(
id).getNumConnections();
3288 npe.push_back(e.neuralNetworks.at(
id).getNumConnections());
3296 vector<size_t> offset;
3300 offset.push_back(n);
3301 n += e.neuralNetworks.at(
id).getNumConnections();
3309 vector<vector<double>>& dPdc)
3315 if (property ==
"energy")
3319 for (
auto const& a : structure.
atoms)
3321 size_t e = a.element;
3325 vector<double> tmp(npe.at(e), 0.0);
3327 for (
size_t j = 0; j < tmp.size(); ++j)
3329 dPdc.at(0).at(off.at(e) + j) += tmp.at(j);
3333 else if (property ==
"force")
3337 for (
size_t ia = 0; ia < structure.
numAtoms; ++ia)
3339 for (
size_t ixyz = 0; ixyz < 3; ++ixyz)
3342 for (
auto& a : structure.
atoms)
3344#ifndef N2P2_FULL_SFD_MEMORY
3347 a.collectDGdxia(ia, ixyz);
3349 size_t e = a.element;
3355 vector<double> tmp(npe.at(e), 0.0);
3356#ifndef N2P2_FULL_SFD_MEMORY
3361 for (
size_t j = 0; j < tmp.size(); ++j)
3363 dPdc.at(count).at(off.at(e) + j) += tmp.at(j);
3372 throw runtime_error(
"ERROR: Weight derivatives not implemented for "
3373 "property \"" + property +
"\".\n");
3381 vector<vector<double>>& dPdc,
3388 if (property ==
"energy")
3393 for (
size_t ic = 0; ic < npe.at(ie); ++ic)
3395 size_t const o = off.at(ie) + ic;
3396 double const w =
weights.at(0).at(o);
3402 double energyHigh = structure.
energy;
3404 weights.at(0).at(o) -= 2.0 * delta;
3408 double energyLow = structure.
energy;
3410 dPdc.at(0).push_back((energyHigh - energyLow) / (2.0 * delta));
3415 else if (property ==
"force")
3419 for (
size_t ia = 0; ia < structure.
numAtoms; ++ia)
3421 for (
size_t ixyz = 0; ixyz < 3; ++ixyz)
3425 for (
size_t ic = 0; ic < npe.at(ie); ++ic)
3427 size_t const o = off.at(ie) + ic;
3428 double const w =
weights.at(0).at(o);
3434 double forceHigh = structure.
atoms.at(ia).f[ixyz];
3436 weights.at(0).at(o) -= 2.0 * delta;
3440 double forceLow = structure.
atoms.at(ia).f[ixyz];
3442 dPdc.at(count).push_back((forceHigh - forceLow)
3453 throw runtime_error(
"ERROR: Numeric weight derivatives not "
3454 "implemented for property \""
3455 + property +
"\".\n");
3535 string s =
strpr(
" E %5zu %10zu %5d %3zu %10.2E %10zu %5zu\n",
3551 string s =
strpr(
" F %5zu %10zu %5d %3zu %10.2E %10zu %5zu %5zu %2zu\n",
3566 string s =
strpr(
" Q %5zu %10zu %5d %3zu %10.2E %10zu %5zu %5zu\n",
3573#ifndef N2P2_FULL_SFD_MEMORY
3576 size_t indexComponent)
3584 dGdxia.resize(nsf+1, 0.0);
3588 vector<vector<size_t> >
const& tableFull
3594 if (atom.
neighbors[i].index == indexAtom)
3597 vector<size_t>
const& table = tableFull.at(n.
element);
3598 for (
size_t j = 0; j < n.
dGdr.size(); ++j)
3600 dGdxia[table.at(j)] += n.
dGdr[j][indexComponent];
3604 if (atom.
index == indexAtom)
3606 for (
size_t i = 0; i < nsf; ++i)
3618 string keywordNW =
"nguyen_widrow_weights" +
nns.at(
id).keywordSuffix2;
3620 double minWeights = atof(
settings[
"weights_min"].c_str());
3621 double maxWeights = atof(
settings[
"weights_max"].c_str());
3622 log <<
strpr(
"Initial weights selected randomly in interval "
3623 "[%f, %f).\n", minWeights, maxWeights);
3629 for (
size_t j = 0; j < w.size(); ++j)
3631 w.at(j) = minWeights + gsl_rng_uniform(
rngGlobal)
3632 * (maxWeights - minWeights);
3638 log <<
"Weights modified according to Nguyen Widrow scheme.\n";
3639 for (vector<Element>::iterator it =
elements.begin();
3648 throw runtime_error(
"ERROR: Preconditioning of weights not yet"
3653 log <<
"Weights modified accoring to Glorot Bengio scheme.\n";
3655 log <<
"Biases set to zero.\n";
3656 for (vector<Element>::iterator it =
elements.begin();
3671 bool all = (
property ==
"all");
3672 bool isProperty = (find(
pk.begin(),
pk.end(), property) !=
pk.end());
3673 if (!(all || isProperty))
3675 throw runtime_error(
"ERROR: Unknown property for selection mode"
3684 log <<
"Global selection mode settings:\n";
3691 + property)))
return;
3692 log <<
"Selection mode settings specific to property \""
3693 <<
property <<
"\":\n";
3696 if (all) keyword =
"selection_mode";
3697 else keyword =
"selection_mode_" + property;
3701 map<size_t, SelectionMode> schedule;
3703 if (args.size() % 2 != 1)
3705 throw runtime_error(
"ERROR: Incorrect selection mode format.\n");
3708 for (
size_t i = 1; i < args.size(); i = i + 2)
3710 schedule[(size_t)atoi(args.at(i).c_str())] =
3713 for (map<size_t, SelectionMode>::const_iterator it = schedule.begin();
3714 it != schedule.end(); ++it)
3716 log <<
strpr(
"- Selection mode starting with epoch %zu:\n",
3720 log <<
strpr(
" Random selection of update candidates: "
3721 "SelectionMode::SM_RANDOM (%d)\n", it->second);
3723 else if (it->second ==
SM_SORT)
3725 log <<
strpr(
" Update candidates selected according to error: "
3726 "SelectionMode::SM_SORT (%d)\n", it->second);
3730 log <<
strpr(
" Update candidates chosen randomly above RMSE "
3731 "threshold: SelectionMode::SM_THRESHOLD (%d)\n",
3736 throw runtime_error(
"ERROR: Unknown selection mode.\n");
3743 i.second.selectionModeSchedule = schedule;
3744 i.second.selectionMode = schedule[0];
3749 p[property].selectionModeSchedule = schedule;
3750 p[property].selectionMode = schedule[0];
3754 if (all) keyword =
"rmse_threshold";
3755 else keyword =
"rmse_threshold_" + property;
3758 double t = atof(
settings[keyword].c_str());
3759 log <<
strpr(
"- RMSE selection threshold: %.2f * RMSE\n", t);
3760 if (all)
for (
auto& i :
p) i.second.rmseThreshold = t;
3761 else p[property].rmseThreshold = t;
3764 if (all) keyword =
"rmse_threshold_trials";
3765 else keyword =
"rmse_threshold_trials_" + property;
3768 size_t t = atoi(
settings[keyword].c_str());
3769 log <<
strpr(
"- RMSE selection trials : %zu\n", t);
3770 if (all)
for (
auto& i :
p) i.second.rmseThresholdTrials = t;
3771 else p[property].rmseThresholdTrials = t;
3779 string keyword =
"write_";
3780 bool isProperty = (find(
pk.begin(),
pk.end(), type) !=
pk.end());
3781 if (type ==
"energy" ) keyword +=
"trainpoints";
3782 else if (type ==
"force" ) keyword +=
"trainforces";
3783 else if (type ==
"charge" ) keyword +=
"traincharges";
3784 else if (type ==
"weights_epoch") keyword += type;
3785 else if (type ==
"neuronstats" ) keyword += type;
3788 throw runtime_error(
"ERROR: Invalid type for file output setup.\n");
3794 size_t* writeEvery =
nullptr;
3795 size_t* writeAlways =
nullptr;
3799 writeEvery = &(
p[type].writeCompEvery);
3800 writeAlways = &(
p[type].writeCompAlways);
3801 message =
"Property \"" + type +
"\" comparison";
3802 message.at(0) = toupper(message.at(0));
3804 else if (type ==
"weights_epoch")
3810 else if (type ==
"neuronstats")
3814 message =
"Neuron statistics";
3819 if (v.size() == 1) *writeEvery = (size_t)atoi(v.at(0).c_str());
3820 else if (v.size() == 2)
3822 *writeEvery = (size_t)atoi(v.at(0).c_str());
3823 *writeAlways = (size_t)atoi(v.at(1).c_str());
3826 +
" files will be written every %zu epochs.\n").c_str(),
3828 if (*writeAlways > 0)
3831 +
" files will always be written up to epoch "
3832 "%zu.\n").c_str(), *writeAlways);
3841 bool isProperty = (find(
pk.begin(),
pk.end(), property) !=
pk.end());
3844 throw runtime_error(
"ERROR: Unknown property for update plan"
3850 string keyword =
property +
"_fraction";
3853 if (property ==
"force" &&
3857 double const ratio = atof(
settings[
"force_energy_ratio"].c_str());
3860 log <<
"WARNING: Given force fraction is ignored because "
3861 "force/energy ratio is provided.\n";
3863 log <<
strpr(
"Desired force/energy update ratio : %.6f\n",
3865 log <<
"----------------------------------------------\n";
3867 /
p[
"force"].numTrainPatterns;
3875 keyword =
"task_batch_size_" + property;
3901 MPI_Allgather(MPI_IN_PLACE, 1, MPI_INT, &(pa.
errorsPerTask.front()), 1, MPI_INT,
comm);
3930 log <<
"Update plan for property \"" +
property +
"\":\n";
3931 log <<
strpr(
"- Per-task batch size : %zu\n",
3933 log <<
strpr(
"- Fraction of patterns used per epoch : %.6f\n",
3937 log <<
"WARNING: No updates are planned for this property.";
3939 log <<
strpr(
"- Updates per epoch : %zu\n",
3941 log <<
strpr(
"- Patterns used per update (rank %3d / global) : "
3944 log <<
"----------------------------------------------\n";
3951 bool isProperty = (find(
pk.begin(),
pk.end(), property) !=
pk.end());
3954 throw runtime_error(
"ERROR: Unknown property for array allocation.\n");
3957 log <<
"Allocating memory for " +
property +
3958 " error vector and Jacobian.\n";
3976 pa.
error.at(i).resize(size, 0.0);
3978 log <<
strpr(
"Updater %3zu:\n", i);
3979 log <<
strpr(
" - Error size: %zu\n", pa.
error.at(i).size());
3982 log <<
"----------------------------------------------\n";
3990 string fileNameActual = fileName;
3994 fileNameActual +=
strpr(
".stage-%zu",
stage);
3997 vector<string> sub = {
"_err",
"_com",
"_upd",
"_log"};
3998 if (append) file.open(fileNameActual.c_str(), ofstream::app);
4001 file.open(fileNameActual.c_str());
4004 vector<string> title;
4005 vector<string> colName;
4006 vector<string> colInfo;
4007 vector<size_t> colSize;
4008 title.push_back(
"Timing data for training loop.");
4009 colSize.push_back(10);
4010 colName.push_back(
"epoch");
4011 colInfo.push_back(
"Current epoch.");
4012 colSize.push_back(11);
4013 colName.push_back(
"train");
4014 colInfo.push_back(
"Time for training.");
4015 colSize.push_back(7);
4016 colName.push_back(
"ptrain");
4017 colInfo.push_back(
"Time for training (percentage of loop).");
4018 colSize.push_back(11);
4019 colName.push_back(
"error");
4020 colInfo.push_back(
"Time for error calculation.");
4021 colSize.push_back(7);
4022 colName.push_back(
"perror");
4023 colInfo.push_back(
"Time for error calculation (percentage of loop).");
4024 colSize.push_back(11);
4025 colName.push_back(
"epoch");
4026 colInfo.push_back(
"Time for this epoch.");
4027 colSize.push_back(11);
4028 colName.push_back(
"total");
4029 colInfo.push_back(
"Total time for all epochs.");
4032 colSize.push_back(11);
4033 colName.push_back(
p[k].tiny +
"train");
4034 colInfo.push_back(
"");
4035 colSize.push_back(7);
4036 colName.push_back(
p[k].tiny +
"ptrain");
4037 colInfo.push_back(
"");
4043 colSize.push_back(11);
4044 colName.push_back(
p[k].tiny + s);
4045 colInfo.push_back(
"");
4046 colSize.push_back(7);
4047 colName.push_back(
p[k].tiny +
"p" + s);
4048 colInfo.push_back(
"");
4055 double timeLoop =
sw[
"loop"].getLoop();
4057 file <<
strpr(
" %11.3E",
sw[
"train"].getLoop());
4058 file <<
strpr(
" %7.3f",
sw[
"train"].getLoop() / timeLoop);
4059 file <<
strpr(
" %11.3E",
sw[
"error"].getLoop());
4060 file <<
strpr(
" %7.3f",
sw[
"error"].getLoop() / timeLoop);
4061 file <<
strpr(
" %11.3E", timeLoop);
4062 file <<
strpr(
" %11.3E",
sw[
"loop"].getTotal());
4066 file <<
strpr(
" %11.3E",
sw[k].getLoop());
4067 file <<
strpr(
" %7.3f",
sw[k].getLoop() /
sw[
"train"].getLoop());
4073 file <<
strpr(
" %11.3E",
sw[k + s].getLoop());
4074 file <<
strpr(
" %7.3f",
sw[k + s].getLoop() /
sw[k].getLoop());
4086 property (property ),
4087 displayMetric (
"" ),
4090 selectionMode (SM_RANDOM),
4091 numTrainPatterns (0 ),
4092 numTestPatterns (0 ),
4094 writeCompEvery (0 ),
4095 writeCompAlways (0 ),
4096 posUpdateCandidates (0 ),
4097 numGroupedSubCand (0 ),
4098 countGroupedSubCand (0 ),
4099 rmseThresholdTrials (0 ),
4102 patternsPerUpdate (0 ),
4103 patternsPerUpdateGlobal(0 ),
4104 numErrorsGlobal (0 ),
4105 epochFraction (0.0 ),
4106 rmseThreshold (0.0 )
4128 throw runtime_error(
"ERROR: Unknown training property.\n");
Collect and process large data sets.
std::size_t numStructures
Total number of structures in dataset.
gsl_rng * rng
GSL random number generator (different seed for each MPI process).
MPI_Comm comm
Global MPI communicator.
int numProcs
Total number of MPI processors.
std::vector< Structure > structures
All structures in this dataset.
void combineFiles(std::string filePrefix) const
Combine individual MPI proc files to one.
gsl_rng * rngGlobal
Global GSL random number generator (equal seed for each MPI process).
void collectError(std::string const &property, std::map< std::string, double > &error, std::size_t &count) const
Collect error metrics of a property over all MPI procs.
std::size_t atomicNumber(std::size_t index) const
Get atomic number from element index.
Weight updates based on simple gradient descent methods.
void setParametersFixed(double const eta)
Set parameters for fixed step gradient descent algorithm.
void setParametersAdam(double const eta, double const beta1, double const beta2, double const epsilon)
Set parameters for Adam algorithm.
DescentType
Enumerate different gradient descent variants.
@ DT_ADAM
Adaptive moment estimation (Adam).
@ DT_FIXED
Fixed step size.
Implementation of the Kalman filter method.
void setParametersStandard(double const epsilon, double const q0, double const qtau, double const qmin, double const eta0, double const etatau, double const etamax)
Set parameters for standard Kalman filter.
void setParametersFadingMemory(double const epsilon, double const q0, double const qtau, double const qmin, double const lambda, double const nu)
Set parameters for fading memory Kalman filter.
KalmanType
Enumerate different Kalman filter types.
@ KT_STANDARD
Regular Kalman filter.
@ KT_FADINGMEMORY
Kalman filtering with fading memory modification.
void setSizeObservation(std::size_t const sizeObservation)
Set observation vector size.
bool silent
Completely silence output.
std::vector< std::vector< double > > cutoffs
Matrix storing all symmetry function cut-offs for all elements.
void setupNormalization(bool standalone=true)
Set up normalization.
ElementMap elementMap
Global element map, populated by setupElementMap().
void addEnergyOffset(Structure &structure, bool ref=true)
Add atomic energy offsets to reference energy.
void readNeuralNetworkWeights(std::string const &id, std::string const &fileName)
Read in weights for a specific type of neural network.
virtual void setupSymmetryFunctionCache(bool verbose=false)
Set up symmetry function cache.
@ HDNNP_2G
Short range NNP (2G-HDNNP).
@ HDNNP_Q
Short range NNP with charge NN, no electrostatics/Qeq (M.
@ HDNNP_4G
NNP with electrostatics and non-local charge transfer (4G-HDNNP).
ScreeningFunction screeningFunction
void calculateAtomicNeuralNetworks(Structure &structure, bool const derivatives, std::string id="")
Calculate atomic neural networks for all atoms in given structure.
double physical(std::string const &property, double value) const
Undo normalization for a given property.
virtual void setupSymmetryFunctionGroups()
Set up symmetry function groups.
std::vector< Element > elements
double normalized(std::string const &property, double value) const
Apply normalization to given property.
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.
settings::Settings settings
virtual void setupSymmetryFunctionScaling(std::string const &fileName="scaling.data")
Set up symmetry function scaling from file.
void removeEnergyOffset(Structure &structure, bool ref=true)
Remove atomic energy offsets from reference energy.
double getEnergyWithOffset(Structure const &structure, bool ref=true) const
Add atomic energy offsets and return energy.
virtual void setupSymmetryFunctions()
Set up all symmetry functions.
void chargeEquilibration(Structure &structure, bool const derivativesElec)
Perform global charge equilibration method.
void calculateSymmetryFunctions(Structure &structure, bool const derivatives)
Calculate all symmetry functions for all atoms in given structure.
std::map< std::string, NNSetup > nns
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 setupSymmetryFunctionMemory(bool verbose=false)
Extract required memory dimensions for symmetry function derivatives.
void writeSettingsFile(std::ofstream *const &file) const
Write complete settings file.
This class implements a feed-forward neural network.
int getNumConnections() const
Return total number of connections.
void setInput(double const *const &input) const
Set neural network input layer node values.
void modifyConnections(ModificationScheme modificationScheme)
Change connections according to a given modification scheme.
void setConnections(double const *const &connections)
Set neural network weights and biases.
void calculateDFdc(double *dFdc, double const *const &dGdxyz) const
Calculate "second" derivative of output with respect to connections.
@ MS_ZEROBIAS
Set all bias values to zero.
@ MS_GLOROTBENGIO
Normalize connections according to Glorot and Bengio.
@ MS_NGUYENWIDROW
Initialize connections according to Nguyen-Widrow scheme.
void getConnections(double *connections) const
Get neural network weights and biases.
void propagate()
Propagate input information through all layers.
void calculateDEdc(double *dEdc) const
Calculate derivative of output neuron with respect to connections.
void calculateDEdG(double *dEdG) const
Calculate derivative of output neuron with respect to input neurons.
void getOutput(double *output) const
Get neural network output layer node values.
double getOuter() const
Getter for outer.
void setTrainingLogFileName(std::string fileName)
Set training log file name.
UpdaterType updaterType
Updater type used.
std::vector< double > dGdxia
Derivative of symmetry functions with respect to one specific atom coordinate.
std::size_t epoch
Current epoch.
std::vector< std::size_t > numWeightsPerUpdater
Number of weights per updater.
std::size_t writeWeightsAlways
Up to this epoch weights are written every epoch.
std::size_t numWeights
Total number of weights.
bool useForces
Use forces for training.
std::size_t writeWeightsEvery
Write weights every this many epochs.
void setupTraining()
General training settings and setup of weight update routine.
std::mt19937_64 rngGlobalNew
Global random number generator.
std::string nnId
ID of neural network the training is working on.
void writeHardness(std::string const &fileNameFormat) const
Write hardness to files (one file for each element).
void allocateArrays(std::string const &property)
Allocate error and Jacobian arrays for given property.
void writeWeightsEpoch() const
Write weights to files during training loop.
std::ofstream trainingLog
Training log file.
bool hasUpdaters
If this rank performs weight updates.
void checkSelectionMode()
Check if selection mode should be changed.
std::string trainingLogFileName
File name for training log.
JacobianMode jacobianMode
Jacobian mode used.
void selectSets()
Randomly select training and test set structures.
JacobianMode
Jacobian matrix preparation mode.
@ JM_SUM
No Jacobian, sum up contributions from update candidates.
@ JM_TASK
Prepare one Jacobian entry for each task, sum up within tasks.
@ JM_FULL
Prepare full Jacobian matrix.
void randomizeNeuralNetworkWeights(std::string const &id)
Randomly initialize specificy neural network weights.
std::size_t countUpdates
Update counter (for all training quantities together).
std::map< std::string, Stopwatch > sw
Stopwatches for timing overview.
void setupSelectionMode(std::string const &property)
Set selection mode for specific training property.
void dPdc(std::string property, Structure &structure, std::vector< std::vector< double > > &dEdc)
Compute derivatives of property with respect to weights.
void printHeader()
Print training loop header on screen.
std::mt19937_64 rngNew
Per-task random number generator.
PropertyMap p
Actual training properties.
std::vector< Updater * > updaters
Weight updater (combined or for each element).
std::size_t writeNeuronStatisticsAlways
Up to this epoch neuron statistics are written every epoch.
UpdaterType
Type of update routine.
@ UT_GD
Simple gradient descent methods.
@ UT_KF
Kalman filter-based methods.
@ UT_LM
Levenberg-Marquardt algorithm.
double forceWeight
Force update weight.
void update(std::string const &property)
Perform one update.
void dataSetNormalization()
Apply normalization based on initial weights prediction.
ParallelMode
Training parallelization mode.
@ PM_TRAIN_RK0
No training parallelization, only data set distribution.
@ PM_TRAIN_ALL
Parallel gradient computation, update on each task.
bool writeTrainingLog
Whether training log file is written.
std::size_t numUpdaters
Number of updaters (depends on update strategy).
bool advance() const
Check if training loop should be continued.
~Training()
Destructor, updater vector needs to be cleaned.
void sortUpdateCandidates(std::string const &property)
Sort update candidates with descending RMSE.
void writeNeuronStatisticsEpoch() const
Write neuron statistics during training loop.
SelectionMode
How update candidates are selected during Training.
@ SM_RANDOM
Select candidates randomly.
@ SM_THRESHOLD
Select candidates randomly with RMSE above threshold.
@ SM_SORT
Sort candidates according to their RMSE and pick worst first.
std::size_t stage
Training stage.
void collectDGdxia(Atom const &atom, std::size_t indexAtom, std::size_t indexComponent)
Collect derivative of symmetry functions with respect to one atom's coordinate.
void writeHardnessEpoch() const
Write hardness to files during training loop.
void setWeights()
Set weights in neural network class.
void loop()
Execute main training loop.
void printEpoch()
Print preferred error metric and timing information on screen.
std::vector< std::size_t > weightsOffset
Offset of each element's weights in combined array.
std::size_t numEpochs
Number of epochs requested.
void setupFileOutput(std::string const &type)
Set file output intervals for properties and other quantities.
std::size_t getNumConnections(std::string id="short") const
Get total number of NN connections.
void addTrainingLogEntry(int proc, std::size_t il, double f, std::size_t isg, std::size_t is)
Write energy update data to training log file.
void writeTimingData(bool append, std::string const fileName="timing.out")
Write timing data for all clocks.
std::size_t writeNeuronStatisticsEvery
Write neuron statistics every this many epochs.
void getWeights()
Get weights from neural network class.
void dPdcN(std::string property, Structure &structure, std::vector< std::vector< double > > &dEdc, double delta=1.0E-4)
Compute numeric derivatives of property with respect to weights.
std::vector< std::string > setupNumericDerivCheck()
Set up numeric weight derivatives check.
void setSingleWeight(std::size_t element, std::size_t index, double value)
Set a single weight value.
void shuffleUpdateCandidates(std::string const &property)
Shuffle update candidates.
void calculateChargeErrorVec(Structure const &s, Eigen::VectorXd &cVec, double &cNorm)
Calculate vector of charge errors in one structure.
double getSingleWeight(std::size_t element, std::size_t index)
Get a single weight value.
void initializeWeights()
Initialize weights for all elements.
void initializeWeightsMemory(UpdateStrategy updateStrategy=US_COMBINED)
Initialize weights vector according to update strategy.
bool repeatedEnergyUpdates
After force update perform energy update for corresponding structure.
UpdateStrategy updateStrategy
Update strategy used.
void calculateNeighborLists()
Calculate neighbor lists for all structures.
ParallelMode parallelMode
Parallelization mode used.
void writeNeuronStatistics(std::string const &nnName, std::string const &fileName) const
Write neuron statistics collected since last invocation.
std::vector< std::size_t > getConnectionOffsets(std::string id="short") const
Get offsets of NN connections for each element.
std::vector< int > epochSchedule
Update schedule epoch (stage 1: 0 = charge update; stage 2: 0 = energy update, 1 = force update (opti...
void writeUpdaterStatus(bool append, std::string const fileNameFormat="updater.%03zu.out") const
Write updater information to file.
UpdateStrategy
Update strategies available for Training.
@ US_COMBINED
One combined updater for all elements.
@ US_ELEMENT
Separate updaters for individual elements.
std::vector< std::size_t > getNumConnectionsPerElement(std::string id="short") const
Get number of NN connections for each element.
void resetNeuronStatistics()
Reset neuron statistics for all elements.
std::vector< std::vector< double > > calculateWeightDerivatives(Structure *structure)
Calculate derivatives of energy with respect to weights.
void calculateError(std::map< std::string, std::pair< std::string, std::string > > const fileNames)
Calculate error metrics for all structures.
void writeLearningCurve(bool append, std::string const fileName="learning-curve.out") const
Write current RMSEs and epoch information to file.
void calculateErrorEpoch()
Calculate error metrics per epoch for all structures with file names used in training loop.
std::vector< std::string > pk
Vector of actually used training properties.
void setEpochSchedule()
Select energies/forces schedule for one epoch.
bool hasStructures
If this rank holds structure information.
std::vector< std::vector< double > > weights
Neural network weights and biases for each element.
bool freeMemory
Free symmetry function memory after calculation.
void writeSetsToFiles()
Write training and test set to separate files (train.data and test.data, same format as input....
void writeWeights(std::string const &nnName, std::string const &fileNameFormat) const
Write weights to files (one file for each element).
void setStage(std::size_t stage)
Set training stage (if multiple stages are needed for NNP type).
void setupUpdatePlan(std::string const &property)
Set up how often properties are updated.
Base class for different weight update methods.
void writeSettingsFile(std::ofstream *const &file, std::map< std::size_t, std::string > const &replacements={}) const
Write complete settings file.
std::vector< std::string > getSettingsLines() const
Get complete settings file.
bool keywordExists(Key const &key, bool const exact=false) const override
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.
string reduce(string const &line, string const &whitespace, string const &fill)
Replace multiple whitespaces with fill.
void appendLinesToFile(ofstream &file, vector< string > const lines)
Append multiple lines of strings to open file stream.
Struct to store information on neighbor atoms.
std::size_t element
Element index of neighbor atom.
std::vector< Vec3D > dGdr
Derivatives of symmetry functions with respect to neighbor coordinates.
Storage for a single atom.
std::vector< Neighbor > neighbors
Neighbor array (maximum number defined in macros.h.
std::size_t numSymmetryFunctions
Number of symmetry functions used to describe the atom environment.
Vec3D f
Force vector calculated by neural network.
double charge
Atomic charge determined by neural network.
std::size_t index
Index number of this atom.
Vec3D fRef
Reference force vector from data set.
std::vector< Vec3D > dGdr
Derivative of symmetry functions with respect to this atom's coordinates.
double chi
Atomic electronegativity determined by neural network.
std::size_t element
Element index of this atom.
std::vector< double > G
Symmetry function values.
double chargeRef
Atomic reference charge.
std::vector< std::size_t > numNeighborsPerElement
Number of neighbors per element.
std::vector< double > dChidG
Derivative of electronegativity with respect to symmetry functions.
std::size_t numNeighbors
Total number of neighbors.
Setup data for one neural network.
std::string weightFileFormat
Format for weight files.
std::string keywordSuffix2
Suffix for some other keywords (weight file loading related).
std::string name
Description string for log output, e.g. "electronegativity".
Storage for one atomic configuration.
void freeAtoms(bool all, double const maxCutoffRadius=0.0)
Free symmetry function memory for all atoms, see free() in Atom class.
@ ST_TRAINING
Structure is part of the training set.
@ ST_UNKNOWN
Sample type not assigned yet.
@ ST_TEST
Structure is part of the test set.
std::vector< std::size_t > numAtomsPerElement
Number of atoms of each element in this structure.
void calculateDQdChi(std::vector< Eigen::VectorXd > &dQdChi)
Calculates derivative of the charges with respect to electronegativities.
std::size_t index
Index number of this structure.
SampleType sampleType
Sample type (training or test set).
void calculateDQdJ(std::vector< Eigen::VectorXd > &dQdJ)
Calculates derivative of the charges with respect to atomic hardness.
bool hasAMatrix
If A matrix of this structure is currently stored.
void calculateDQdr(std::vector< size_t > const &atomsIndices, std::vector< size_t > const &compIndices, double const maxCutoffRadius, std::vector< Element > const &elements)
Calculates derivative of the charges with respect to the atom's position.
void clearElectrostatics(bool clearDQdr=false)
Clear A-matrix, dAdrQ and optionally dQdr.
double energy
Potential energy determined by neural network.
double energyRef
Reference potential energy.
std::size_t numAtoms
Total number of atoms present in this structure.
std::vector< Atom > atoms
Vector of all atoms in this structure.
bool hasCharges
If all charges of this structure have been calculated (and stay the same, e.g.
bool exists(std::string const &key)
Check if property is present.
Specific training quantity (e.g. energies, forces, charges).
std::size_t countGroupedSubCand
Counter for number of used subcandidates belonging to same update candidate.
std::size_t numErrorsGlobal
Global number of errors per update.
std::size_t numTrainPatterns
Number of training patterns in set.
std::vector< int > errorsPerTask
Errors per task for each update.
std::string property
Copy of identifier within Property map.
SelectionMode selectionMode
Selection mode for update candidates.
std::size_t numUpdates
Number of desired updates per epoch.
double rmseThreshold
RMSE threshold for update candidates.
std::string plural
Plural string of property;.
double epochFraction
Desired update fraction per epoch.
std::size_t patternsPerUpdate
Patterns used per update.
std::size_t taskBatchSize
Batch size for each MPI task.
std::vector< std::vector< double > > error
Global error vector (per updater).
std::string displayMetric
Error metric for display.
std::string tiny
Tiny abbreviation string for property.
std::vector< std::string > errorMetrics
Error metrics available for this property.
Property(std::string const &property)
Constructor.
std::vector< std::vector< int > > offsetJacobian
Stride for Jacobians per task per updater.
std::size_t countUpdates
Counter for updates per epoch.
std::vector< std::vector< int > > weightsPerTask
Weights per task per updater.
std::size_t numGroupedSubCand
Number of subcandidates which are considered before changing the update candidate.
std::vector< int > offsetPerTask
Offset for combined error per task.
std::vector< UpdateCandidate > updateCandidates
Vector with indices of training patterns.
std::map< std::string, double > errorTrain
Current error metrics of training patterns.
std::vector< std::vector< double > > jacobian
Global Jacobian (per updater).
std::size_t patternsPerUpdateGlobal
Patterns used per update (summed over all MPI tasks).
std::map< std::string, double > errorTest
Current error metrics of test patterns.
std::size_t rmseThresholdTrials
Maximum trials for SM_THRESHOLD selection mode.
std::size_t posUpdateCandidates
Current position in update candidate list (SM_SORT/_THRESHOLD).
Contains update candidate which is grouped with others to specific parent update candidate (e....
std::size_t a
Atom index (only used for force candidates).
std::size_t c
Component index (x,y,z -> 0,1,2, only used for force candidates).
Contains location of one update candidate (energy or force).
std::vector< SubCandidate > subCandidates
Vector containing grouped candidates.
std::size_t posSubCandidates
Current position in sub-candidate list (SM_SORT/_THRESHOLD).
std::size_t s
Structure index.