update frame parameter according to LAMMPS compute command

doc/third-party/lammps-command.md

@@ -29,14 +29,16 @@ pair_style deepmd models ... keyword value ...
 - models = frozen model(s) to compute the interaction. 
 If multiple models are provided, then only the first model serves to provide energy and force prediction for each timestep of molecular dynamics, 
 and the model deviation will be computed among all models every `out_freq` timesteps.
-- keyword = *out_file* or *out_freq* or *fparam* or *atomic* or *relative* or *relative_v* or *aparam* or *ttm*
+- keyword = *out_file* or *out_freq* or *fparam* or *fparam_from_compute* or *atomic* or *relative* or *relative_v* or *aparam* or *ttm*
 <pre>
     <i>out_file</i> value = filename
         filename = The file name for the model deviation output. Default is model_devi.out
     <i>out_freq</i> value = freq
         freq = Frequency for the model deviation output. Default is 100.
     <i>fparam</i> value = parameters
         parameters = one or more frame parameters required for model evaluation.
+    <i>fparam_from_compute</i> value = id
+        id = compute id used to update the frame parameter.
     <i>atomic</i> = no value is required. 
         If this keyword is set, the model deviation of each atom will be output.
     <i>relative</i> value = level
@@ -55,6 +57,9 @@ pair_style deepmd graph.pb
 pair_style deepmd graph.pb fparam 1.2
 pair_style deepmd graph_0.pb graph_1.pb graph_2.pb out_file md.out out_freq 10 atomic relative 1.0
 pair_coeff * * O H
+
+pair_style deepmd cp.pb fparam_from_compute TEMP
+compute    TEMP all temp
 ```
 
 ### Description
@@ -74,6 +79,7 @@ If the keyword `relative_v` is set, then the relative model deviation of the vir
 $$E_{v_i}=\frac{\left|D_{v_i}\right|}{\left|v_i\right|+l}$$
 
 If the keyword `fparam` is set, the given frame parameter(s) will be fed to the model.
+If the keyword `fparam_from_compute` is set, the global parameter(s) from compute command (e.g., temperature from [compute temp command](https://docs.lammps.org/compute_temp.html)) will be fed to the model as the frame parameter(s).
 If the keyword `aparam` is set, the given atomic parameter(s) will be fed to the model, where each atom is assumed to have the same atomic parameter(s). 
 If the keyword `ttm` is set, electronic temperatures from [fix ttm command](https://docs.lammps.org/fix_ttm.html) will be fed to the model as the atomic parameters.
 

source/lmp/pair_deepmd.cpp

@@ -7,6 +7,7 @@
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
+#include "compute.h"
 #include "memory.h"
 #include "update.h"
 #include "output.h"
@@ -179,6 +180,81 @@ make_uniform_aparam(
   }
 }
 
+void PairDeepMD::make_fparam_from_compute(
+#ifdef HIGH_PREC
+  vector<double > & fparam
+#else
+  vector<float > & fparam
+#endif
+)
+{
+  assert(do_compute);
+
+  int icompute = modify->find_compute(compute_id);
+  Compute *compute = modify->compute[icompute];
+
+  assert(compute);
+  fparam.resize(dim_fparam);
+
+  if (dim_fparam == 1){
+    compute->compute_scalar();
+    fparam[0] = compute->scalar;
+  }
+  else if (dim_fparam >1){
+    compute->compute_vector();
+    double *cvector = compute->vector;
+    for (int jj = 0; jj < dim_aparam; ++jj){
+    fparam[jj] =cvector[jj];
+    }
+  }
+}
+
+#ifdef USE_TTM
+void PairDeepMD::make_ttm_fparam(
+#ifdef HIGH_PREC
+    vector<double > & fparam
+#else
+    vector<float > & fparam
+#endif
+    )
+{
+  assert(do_ttm);
+  // get ttm_fix
+  const FixTTMDP * ttm_fix = NULL;
+  for (int ii = 0; ii < modify->nfix; ii++) {
+    if (string(modify->fix[ii]->id) == ttm_fix_id){
+      ttm_fix = dynamic_cast<FixTTMDP*>(modify->fix[ii]);
+    }
+  }
+  assert(ttm_fix);
+
+  fparam.resize(dim_fparam);
+
+  vector<int> nnodes = ttm_fix->get_nodes();
+  int nxnodes = nnodes[0];
+  int nynodes = nnodes[1];
+  int nznodes = nnodes[2];
+  double *** const T_electron = ttm_fix->get_T_electron();
+
+  int numb_effective_nodes = 0;
+  double total_Te = 0;
+
+  // loop over grids to get average electron temperature 
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++) 
+        {
+          if (T_electron[ixnode][iynode][iznode] != 0)
+          {
+            numb_effective_nodes += 1;
+            total_Te += T_electron[ixnode][iynode][iznode];
+          }
+        }
+
+  fparam[0] = total_Te/numb_effective_nodes;
+}
+#endif
+
 #ifdef USE_TTM
 void PairDeepMD::make_ttm_aparam(
 #ifdef HIGH_PREC
@@ -363,9 +439,18 @@ void PairDeepMD::compute(int eflag, int vflag)
   }
   else if (do_ttm) {
 #ifdef USE_TTM
+    if (dim_aparam > 0){
     make_ttm_aparam(daparam);
+    }
+    else if (dim_fparam > 0){
+    make_ttm_fparam(fparam);
+    }
 #endif
   }
+
+  if (do_compute){
+    make_fparam_from_compute(fparam);
+  }
 
   // int ago = numb_models > 1 ? 0 : neighbor->ago;
   int ago = neighbor->ago;
@@ -806,6 +891,7 @@ is_key (const string& input)
   keys.push_back("out_file");
   keys.push_back("fparam");
   keys.push_back("aparam");
+  keys.push_back("fparam_from_compute");
   keys.push_back("ttm");
   keys.push_back("atomic");
   keys.push_back("relative");
@@ -921,6 +1007,22 @@ void PairDeepMD::settings(int narg, char **arg)
       error->all(FLERR, "The deepmd-kit was compiled without support for TTM, please rebuild it with LAMMPS version >=20210831");
 #endif      
     }
+
+    ///////////////////////////////////////////////
+    // pair_style     deepmd cp.pb fparam_from_compute TEMP
+    // compute        TEMP all temp
+    //////////////////////////////////////////////
+    else if (string(arg[iarg]) == string("fparam_from_compute")) {
+      for (int ii = 0; ii < 1; ++ii){
+        if (iarg+1+ii >= narg || is_key(arg[iarg+1+ii])) {
+          error->all(FLERR, "invalid fparam_from_compute key: should be fparam_from_compute compute_id(str)");
+        }
+      }	
+      do_compute = true;
+      compute_id = arg[iarg+1];
+      iarg += 1 + 1;
+    }
+
     else if (string(arg[iarg]) == string("atomic")) {
       out_each = 1;
       iarg += 1;
@@ -948,6 +1050,9 @@ void PairDeepMD::settings(int narg, char **arg)
   if (do_ttm && aparam.size() > 0) {
     error->all(FLERR,"aparam and ttm should NOT be set simultaneously");
   }
+  if (do_compute && fparam.size() > 0) {
+    error->all(FLERR,"fparam and fparam_from_compute should NOT be set simultaneously");
+  }
 
   if (comm->me == 0){
     if (numb_models > 1 && out_freq > 0){
@@ -982,20 +1087,34 @@ void PairDeepMD::settings(int narg, char **arg)
     cout << endl
 	 << pre << "rcut in model:      " << cutoff << endl
 	 << pre << "ntypes in model:    " << numb_types << endl;
-    if (dim_fparam > 0) {
+    if (fparam.size() > 0) {
       cout << pre << "using fparam(s):    " ;
       for (int ii = 0; ii < dim_fparam; ++ii){
 	cout << fparam[ii] << "  " ;
       }
       cout << endl;
     }
+    if (do_compute){
+      cout << pre << "using compute id:      " ;
+      cout << compute_id << "  " << endl;
+    }
     if (aparam.size() > 0) {
       cout << pre << "using aparam(s):    " ;
       for (int ii = 0; ii < aparam.size(); ++ii){
 	cout << aparam[ii] << "  " ;
       }
       cout << endl;
     }
+    if (do_ttm){
+      cout << pre << "using ttm fix:      " ;
+      cout << ttm_fix_id << "  " ;
+      if (dim_fparam > 0){
+        cout << "(fparam)" << endl;
+      }
+      else if (dim_aparam > 0){
+        cout << "(aparam)" << endl;
+      }
+    }  
   }
 
   comm_reverse = numb_models * 3;

source/lmp/pair_deepmd.h.in

@@ -97,6 +97,25 @@ private:
   float eps;
   float eps_v;
 #endif
+
+  void make_fparam_from_compute(
+#ifdef HIGH_PREC
+      std::vector<double > & fparam
+#else
+      std::vector<float > & fparam
+#endif
+      );
+  bool do_compute;
+  std::string compute_id;
+
+    void make_ttm_fparam(
+#ifdef HIGH_PREC
+      std::vector<double > & fparam
+#else
+      std::vector<float > & fparam
+#endif
+      );
+
   void make_ttm_aparam(
 #ifdef HIGH_PREC
       std::vector<double > & dparam