Major update

examples/interface-LAMMPS/H2O_RPBE-D3_Committee/md.lmp

@@ -9,7 +9,7 @@ clear
 # Configuration files
 variable cfgFile         string "h2o_8640_liquid_NpT_RPBE-D3.data"
 # Timesteps
-variable numSteps        equal  5
+variable numSteps        equal  100
 variable dt              equal  0.0005
 # NN
 variable nnpCutoff       equal  6.36
@@ -30,12 +30,12 @@ mass 2 ${mass_O}
 timestep ${dt}
 thermo 1
 thermo_modify format float %14.6f
-dump 4a all custom 5 dump.myforce.* id type fx fy fz
+dump 4a all custom 1 dump.myforce.* id type fx fy fz
 
 ###############################################################################
 # NN
 ###############################################################################
-pair_style nnp dir ${nnpDir} showew no showewsum 10 resetew no maxew 100 maxcd 0.001 cflength 1.8897261328 cfenergy 0.0367493254 emap "1:H,2:O"
+pair_style nnp dir ${nnpDir} showew no showewsum 10 resetew no maxew 100 calcd 1 wricde 1.6 maxcde 3.0 wricdf 0.04 maxcdf 0.10 cflength 1.8897261328 cfenergy 0.0367493254 emap "1:H,2:O"
 pair_coeff * * ${nnpCutoff}
 
 ###############################################################################

src/.vscode/settings.json

@@ -0,0 +1,69 @@
+{
+    "files.associations": {
+        "iosfwd": "cpp",
+        "array": "cpp",
+        "atomic": "cpp",
+        "bit": "cpp",
+        "*.tcc": "cpp",
+        "cctype": "cpp",
+        "cinttypes": "cpp",
+        "clocale": "cpp",
+        "cmath": "cpp",
+        "complex": "cpp",
+        "cstdarg": "cpp",
+        "cstddef": "cpp",
+        "cstdint": "cpp",
+        "cstdio": "cpp",
+        "cstdlib": "cpp",
+        "cstring": "cpp",
+        "ctime": "cpp",
+        "cwchar": "cpp",
+        "cwctype": "cpp",
+        "deque": "cpp",
+        "map": "cpp",
+        "set": "cpp",
+        "unordered_map": "cpp",
+        "vector": "cpp",
+        "exception": "cpp",
+        "algorithm": "cpp",
+        "functional": "cpp",
+        "iterator": "cpp",
+        "memory": "cpp",
+        "memory_resource": "cpp",
+        "numeric": "cpp",
+        "optional": "cpp",
+        "random": "cpp",
+        "string": "cpp",
+        "string_view": "cpp",
+        "system_error": "cpp",
+        "tuple": "cpp",
+        "type_traits": "cpp",
+        "utility": "cpp",
+        "fstream": "cpp",
+        "initializer_list": "cpp",
+        "iomanip": "cpp",
+        "iostream": "cpp",
+        "istream": "cpp",
+        "limits": "cpp",
+        "new": "cpp",
+        "ostream": "cpp",
+        "sstream": "cpp",
+        "stdexcept": "cpp",
+        "streambuf": "cpp",
+        "typeinfo": "cpp",
+        "regex": "cpp",
+        "orderingmethods": "cpp",
+        "bitset": "cpp",
+        "chrono": "cpp",
+        "condition_variable": "cpp",
+        "csetjmp": "cpp",
+        "csignal": "cpp",
+        "list": "cpp",
+        "unordered_set": "cpp",
+        "ratio": "cpp",
+        "mutex": "cpp",
+        "shared_mutex": "cpp",
+        "thread": "cpp",
+        "cfenv": "cpp"
+    }
+}

src/interface/LAMMPS/src/USER-NNP/pair_nnp.cpp

@@ -58,23 +58,30 @@ void PairNNP::compute(int eflag, int vflag)
     handleExtrapolationWarnings();
   }
 
-  // Calculate forces of local and ghost atoms.
-  interface.getForces(atom->f);
-
+  // Get committee size from setting file
   comSize = interface.getComSize();
+
+  // Calculate forces of local and ghost atoms and committee disagreement.
+  if (comSize > 1 && calcd > 0 && update->ntimestep % calcd == 0 && maxcdf > 0.000001)
+  {
+    globalNumAtoms = interface.reduceNumAtoms();
+    interface.getForcesCom(atom->f, globalNumAtoms);
+    std::vector<double> globalForceCom = interface.reduceForceCom(globalNumAtoms);
+    maxForceCom = interface.getComDisForce(globalForceCom, comSize, wricdf, update->ntimestep);
+  }
+  else {
+    interface.getForces(atom->f);
+  }
   // Add energy contribution to total energy.
   if (eflag_global)
   {
     ev_tally(0,0,atom->nlocal,1,interface.getEnergy(),0.0,0.0,0.0,0.0,0.0);
-    if (comSize > 1)
+    if (maxcde > 0.000001 && comSize > 1 && calcd > 0 && update->ntimestep % calcd == 0)
     { 
       std::vector<double> globalEnergyCom = interface.reduceEnergyCom();
       globalNumAtoms = interface.reduceNumAtoms();
-      if(comm->me == 0)
-      {
-        committeeDisagreement = interface.getCommitteeDisagreement(globalEnergyCom, globalNumAtoms);
-        interface.writeCommittee();
-      }
+      if(comm->me == 0){
+        maxEnergyCom = interface.getComDisEnergy(globalEnergyCom, globalNumAtoms, wricde, update->ntimestep);}
     }
   }
 
@@ -106,11 +113,16 @@ void PairNNP::settings(int narg, char **arg)
   maxew = 0;
   globalNumAtoms = 1;
   comSize = 1;
-  maxcd = 0.0;
+  calcd = 1;
+  wricdf = 1.0;
+  wricde = 1.0;
+  maxcdf = 0.0;
+  maxcde = 0.0;
+  maxForceCom = 0.0;
+  maxEnergyCom = 0.0;
   resetew = false;
   cflength = 1.0;
   cfenergy = 1.0;
-  committeeDisagreement = 0.0;
   len = strlen("") + 1;
   emap = new char[len];
   strcpy(emap,"");
@@ -182,11 +194,35 @@ void PairNNP::settings(int narg, char **arg)
         error->all(FLERR,"Illegal pair_style command");
       cfenergy = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
-    // committee disagreement threshold
-    } else if (strcmp(arg[iarg],"maxcd") == 0) {
+    // calculate committee disagreement for each calcd step
+    } else if (strcmp(arg[iarg],"calcd") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal pair_style command");
+      calcd = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    // write to file force committee disagreement if wricd is exceeded
+    } else if (strcmp(arg[iarg],"wricdf") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal pair_style command");
+      wricdf = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    // write to file energy committee disagreement if wricd is exceeded
+      } else if (strcmp(arg[iarg],"wricde") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal pair_style command");
+      wricde = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    // force committee disagreement threshold for stopping calculation
+    } else if (strcmp(arg[iarg],"maxcdf") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal pair_style command");
-      maxcd = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      maxcdf = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    // energy committee disagreement threshold for stopping calculation
+    } else if (strcmp(arg[iarg],"maxcde") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal pair_style command");
+      maxcde = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
     } else error->all(FLERR,"Illegal pair_style command");
   }
@@ -432,10 +468,19 @@ void PairNNP::handleExtrapolationWarnings()
     error->one(FLERR,"Too many extrapolation warnings");
   }
 
+ // Stop if force committee disagreement is exceeded.
+  if (comSize > 1 && calcd > 0 && update->ntimestep % calcd == 0 && maxcdf > 0.000001) {
+    if (maxForceCom > maxcdf) {
+      error->one(FLERR,"Force committee disagreement exceeded set threshold (check maxcdf)");
+    }
+  }
+
  // Stop if energy committee disagreement is exceeded.
-  if (comSize > 1) {
-    if (committeeDisagreement/globalNumAtoms > maxcd) {
-      error->one(FLERR,"Committee energy per atom exceeded set threshold");
+  if(comm->me == 0){
+    if (maxcde > 0.000001 && comSize > 1 && calcd > 0 && update->ntimestep % calcd == 0) {
+      if (maxEnergyCom > maxcde) {
+        error->one(FLERR,"Energy committee disagreement exceeded set threshold (check maxcde)");
+      }
     }
   }
 

src/interface/LAMMPS/src/USER-NNP/pair_nnp.h

@@ -46,15 +46,21 @@ class PairNNP : public Pair {
   int maxew;
   int globalNumAtoms;
   int comSize;
+  int calcd;
   long numExtrapolationWarningsTotal;
   long numExtrapolationWarningsSummary;
   double cflength;
   double cfenergy;
   double maxCutoffRadius;
-  double committeeDisagreement;
-  double maxcd;
+  double wricdf;
+  double wricde;
+  double maxcdf;
+  double maxcde;
+  double maxForceCom;
+  double maxEnergyCom;
   char* directory;
   char* emap;
+
   nnp::InterfaceLammps interface;
 };
 

src/libnnp/Structure.cpp

@@ -908,3 +908,15 @@ double Structure::calcDisagreement() const
 
     return sqrt(energySum / energyCom.size());
 } 
+
+void Structure::addForceLocal(  std::size_t const& alpha,
+                                std::size_t const& c,
+                                std::size_t const& atom,
+                                std::size_t const& comSize, 
+                                double const& fLoc ) const
+{
+    std::size_t index = 3*(atom-1)*comSize + 3*c + alpha;
+    forceLoc[index] -= fLoc;
+
+    return;
+}

src/libnnp/Structure.h

@@ -94,6 +94,8 @@ struct Structure
     Vec3D                    invbox[3];
     /// Same as #energy but for all committee members
     std::vector<double>      energyCom;
+    /// Force vector for all committee members
+    mutable std::vector<double>      forceLoc;
     /// Number of atoms of each element in this structure.
     std::vector<std::size_t> numAtomsPerElement;
     /// Vector of all atoms in this structure.
@@ -341,7 +343,19 @@ struct Structure
      * @return Return committee disagreement for energy.
      */
     double                  calcDisagreement() const;
-
+    /** Add local force contribution to forceLoc vector.
+     * 
+     * @param[in] alpha Coordinate index.
+     * @param[in] c Committee index.
+     * @param[in] atom Atom index.
+     * @param[in] comSize Committee size.
+     * @param[in] fLoc Force contributation added to forceLoc.
+     */
+    void                    addForceLocal(  std::size_t const& alpha,
+                                            std::size_t const& c,
+                                            std::size_t const& atom,
+                                            std::size_t const& comSize,
+                                            double const& fLoc ) const;
 };
 
 }