diff --git a/src/ML-PACE/compute_pace.cpp b/src/ML-PACE/compute_pace.cpp index 070be9ad50..4689106452 100644 --- a/src/ML-PACE/compute_pace.cpp +++ b/src/ML-PACE/compute_pace.cpp @@ -36,6 +36,8 @@ ComputePaceAtom::ComputePaceAtom(class LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 3) error->all(FLERR, "Illegal compute pace/atom command"); peratom_flag = 1; + size_peratom_cols = 0; + scalar_flag = 1; // compute max of gamma } ComputePaceAtom::~ComputePaceAtom() { @@ -43,7 +45,7 @@ ComputePaceAtom::~ComputePaceAtom() { void ComputePaceAtom::init() { if (force->pair == nullptr) - error->all(FLERR, "Compute pace/atom requires a pair style be defined"); + error->all(FLERR, "Compute pace/atom requires a pair style pace/al be defined"); int count = 0; for (int i = 0; i < modify->ncompute; i++) @@ -53,14 +55,23 @@ void ComputePaceAtom::init() { pair_pace_al = force->pair_match("pace/al", 1); if (!pair_pace_al) error->all(FLERR, "Compute pace/atom requires a `pace/al` pair style"); +} +double ComputePaceAtom::compute_scalar() { + invoked_peratom = update->ntimestep; + auto pair = (PairPACEActiveLearning *) pair_pace_al; + if (invoked_peratom != pair->bevaluator_timestep) + error->all(FLERR, "PACE/gamma was not computed on needed timestep.\nIncrease `freq` in pair_style pace/al [gamma_lower_bound] [gamma_upper_bound] [freq] or reset timestep to 0"); + + scalar = pair->max_gamma_grade_per_structure; + return scalar; } void ComputePaceAtom::compute_peratom() { invoked_peratom = update->ntimestep; auto pair = (PairPACEActiveLearning *) pair_pace_al; if (invoked_peratom != pair->bevaluator_timestep) - error->all(FLERR, "PACE/gamma was not computed on needed timestep"); + error->all(FLERR, "PACE/gamma was not computed on needed timestep.\nIncrease `freq` in pair_style pace/al [gamma_lower_bound] [gamma_upper_bound] [freq] or reset timestep to 0"); vector_atom = pair->extrapolation_grade_gamma; - scalar = pair->max_gamma_grade_per_structure; + } \ No newline at end of file diff --git a/src/ML-PACE/compute_pace.h b/src/ML-PACE/compute_pace.h index 1d2f130470..0830ffadcd 100644 --- a/src/ML-PACE/compute_pace.h +++ b/src/ML-PACE/compute_pace.h @@ -35,6 +35,7 @@ namespace LAMMPS_NS { ComputePaceAtom(class LAMMPS *, int, char **); ~ComputePaceAtom() override; void init() override; + double compute_scalar() override; void compute_peratom() override; private: int nmax; diff --git a/src/ML-PACE/pair_pace_al.cpp b/src/ML-PACE/pair_pace_al.cpp index 3718586648..baa14c8c08 100644 --- a/src/ML-PACE/pair_pace_al.cpp +++ b/src/ML-PACE/pair_pace_al.cpp @@ -71,6 +71,8 @@ using namespace MathConst; #define MAXLINE 1024 #define DELTA 4 #define PACE_AL_EXTRAPOLATION_GRADE_FNAME "grade.dat" +#define EXTRAPOLATIVE_STRUCTURES_FNAME "extrapolative_structures.dat" +#define SPECIES_TYPE_FNAME "species_types.dat" //added YL @@ -144,8 +146,12 @@ PairPACEActiveLearning::~PairPACEActiveLearning() { memory->destroy(extrapolation_grade_gamma); } - if (dump) + if (dump) { delete dump; + dump = nullptr; + } + + //computePaceAtom will be deleted by lmp->modify, as it was registered there } /* ---------------------------------------------------------------------- */ @@ -315,8 +321,9 @@ void PairPACEActiveLearning::compute(int eflag, int vflag) { MPI_Comm_rank(world, &mpi_rank); if (max_gamma_grade_per_structure > gamma_upper_bound) { - if (mpi_rank == 0) dump_extrapolation_grade(update->ntimestep, max_gamma_grade_per_structure); - dump->write(); + if (mpi_rank == 0) + dump_extrapolation_grade(update->ntimestep, max_gamma_grade_per_structure); + if (is_dump_extrapolative_structures) dump->write(); MPI_Barrier(world); if (mpi_rank == 0) { @@ -326,8 +333,10 @@ void PairPACEActiveLearning::compute(int eflag, int vflag) { MPI_Abort(world, 1); //abort properly with error code '1' if not using 4 processes exit(EXIT_FAILURE); } else if (max_gamma_grade_per_structure > gamma_lower_bound) { - if (mpi_rank == 0) dump_extrapolation_grade(update->ntimestep, max_gamma_grade_per_structure); - dump->write(); + if (mpi_rank == 0) + dump_extrapolation_grade(update->ntimestep, max_gamma_grade_per_structure); + if (is_dump_extrapolative_structures) + dump->write(); } } @@ -351,9 +360,10 @@ void PairPACEActiveLearning::allocate() { ------------------------------------------------------------------------- */ void PairPACEActiveLearning::settings(int narg, char **arg) { - if (narg > 3) { + //TODO: make keyword base interface: pair_style pace/al freq 100 gamma_lo 1.5 gamma_hi 10 dump yes + if (narg > 4) { error->all(FLERR, - "Illegal pair_style command. Correct form:\n\tpair_style pace/al [gamma_lower_bound] [gamma_upper_bound]"); + "Illegal pair_style command. Correct form:\n\tpair_style pace/al [gamma_lower_bound] [gamma_upper_bound] [freq] [nodump]"); } if (narg > 0) { @@ -381,12 +391,23 @@ void PairPACEActiveLearning::settings(int narg, char **arg) { "Illegal gamma_grade_eval_freq value: it should be integer number >= 1"); } + if (narg > 3) { + if (strcmp(arg[3], "nodump") == 0) is_dump_extrapolative_structures = false; + if (strcmp(arg[3], "dump") == 0) is_dump_extrapolative_structures = true; + } + if (comm->me == 0) { if (screen) { utils::logmesg(lmp, "ACE/AL version: {}.{}.{}\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY); utils::logmesg(lmp, "Extrapolation grade thresholds (lower/upper): {}/{}\n", gamma_lower_bound, gamma_upper_bound); utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq); + if (is_dump_extrapolative_structures) + utils::logmesg(lmp, "Extrapolative structures will be dumped into {}\n", + EXTRAPOLATIVE_STRUCTURES_FNAME); + else + utils::logmesg(lmp, "No extrapolative structures will be dumped\n"); + } } } @@ -448,7 +469,10 @@ void PairPACEActiveLearning::coeff(int narg, char **arg) { aceimpl->rec_ace->element_type_mapping.init(atom->ntypes + 1); aceimpl->rec_ace->element_type_mapping.fill(-1); //-1 means atom not included into potential - FILE *species_type_file = fopen("species_types.dat", "w"); + FILE *species_type_file = nullptr; + if (is_dump_extrapolative_structures) + species_type_file = fopen(SPECIES_TYPE_FNAME, "w"); + const int n = atom->ntypes; for (int i = 1; i <= n; i++) { char *elemname = elemtypes[i - 1]; @@ -459,9 +483,8 @@ void PairPACEActiveLearning::coeff(int narg, char **arg) { map[i] = -1; if (comm->me == 0) utils::logmesg(lmp, "Skipping LAMMPS atom type #{}(NULL)\n", i); } else { - // dump species types for reconstruction of atomic configurations - fprintf(species_type_file, "%s ", elemname); + if (is_dump_extrapolative_structures) fprintf(species_type_file, "%s ", elemname); int atomic_number = AtomicNumberByName_pace_al(elemname); if (atomic_number == -1) error->all(FLERR, "'{}' is not a valid element\n", elemname); SPECIES_TYPE mu = aceimpl->basis_set->get_species_index_by_name(elemname); @@ -479,7 +502,7 @@ void PairPACEActiveLearning::coeff(int narg, char **arg) { } } } - fclose(species_type_file); + if (is_dump_extrapolative_structures) fclose(species_type_file); aceimpl->ace->set_basis(*aceimpl->basis_set); aceimpl->rec_ace->set_basis(*aceimpl->ctilde_basis_set); @@ -508,14 +531,14 @@ void PairPACEActiveLearning::coeff(int narg, char **arg) { } // prepare dump class - if (!dump) { + if (is_dump_extrapolative_structures && !dump) { // dump pace all custom freq extrapolation.dat id type mass x y z char **dumpargs = new char *[12]; dumpargs[0] = (char *) "pace"; // dump id dumpargs[1] = (char *) "all"; // group dumpargs[2] = (char *) "custom"; // dump style dumpargs[3] = (char *) "1"; // dump frequency - dumpargs[4] = (char *) "extrapolative_structures.dat"; // fname + dumpargs[4] = (char *) EXTRAPOLATIVE_STRUCTURES_FNAME; // fname dumpargs[5] = (char *) "id"; dumpargs[6] = (char *) "type"; dumpargs[7] = (char *) "mass"; diff --git a/src/ML-PACE/pair_pace_al.h b/src/ML-PACE/pair_pace_al.h index 47616cd0b1..7e66603092 100644 --- a/src/ML-PACE/pair_pace_al.h +++ b/src/ML-PACE/pair_pace_al.h @@ -59,6 +59,7 @@ namespace LAMMPS_NS { struct ACEALImpl *aceimpl; int gamma_grade_eval_freq = 1; + bool is_dump_extrapolative_structures = true; DumpCustom *dump = nullptr; Compute *computePaceAtom = nullptr; int natoms; //total number of atoms