Merge pull request #5 from deepmodeling/devel

devel, up to date

Author: denghuilu

README.md

@@ -32,6 +32,7 @@
 	    - [Include deepmd in the pair style](#include-deepmd-in-the-pair-style)
 	    - [Long-range interaction](#long-range-interaction)
 	- [Run path-integral MD with i-PI](#run-path-integral-md-with-i-pi)
+	- [Use deep potential with ASE](#use-deep-potential-with-ase)
 - [Troubleshooting](#troubleshooting)
 
 # About DeePMD-kit
@@ -553,6 +554,30 @@ The option **`graph_file`** provides the file name of the frozen model.
 
 The `dp_ipi` gets the atom names from an [XYZ file](https://en.wikipedia.org/wiki/XYZ_file_format) provided by **`coord_file`** (meanwhile ignores all coordinates in it), and translates the names to atom types by rules provided by **`atom_type`**.
 
+## Use deep potential with ASE
+
+Deep potential can be set up as a calculator with ASE to obtain potential energies and forces.
+```python
+from ase import Atoms
+from deepmd.calculator import DP
+
+water = Atoms('H2O',
+              positions=[(0.7601, 1.9270, 1),
+                         (1.9575, 1, 1),
+                         (1., 1., 1.)],
+              cell=[100, 100, 100],
+              calculator=DP(model="frozen_model.pb"))
+print(water.get_potential_energy())
+print(water.get_forces())
+```
+
+Optimization is also available:
+```python
+from ase.optimize import BFGS
+dyn = BFGS(water)
+dyn.run(fmax=1e-6)
+print(water.get_positions())
+```
 
 # Troubleshooting
 In consequence of various differences of computers or systems, problems may occur. Some common circumstances are listed as follows. 

source/lib/include/NNPInter.h

@@ -200,6 +200,9 @@ class NNPInterModelDevi
   void compute_std_f (vector<VALUETYPE> &		std,
 		      const vector<VALUETYPE> &		avg,
 		      const vector<vector<VALUETYPE> >& xx);
+  void compute_relative_std_f (vector<VALUETYPE> &		std,
+		      const vector<VALUETYPE> &		avg,
+		      const VALUETYPE eps);
 private:
   unsigned numb_models;
   vector<Session*> sessions;

source/lib/src/NNPInter.cc

@@ -1690,3 +1690,22 @@ compute_std_f (vector<VALUETYPE> &		std,
   }
 }
 
+void
+NNPInterModelDevi::
+compute_relative_std_f (vector<VALUETYPE> &std,
+						const vector<VALUETYPE> &avg,
+						const VALUETYPE eps)
+{
+  unsigned nloc = std.size();
+  for (unsigned ii = 0; ii < nloc; ++ii){
+      const VALUETYPE * tmp_avg = &(avg[ii*3]);
+      VALUETYPE vdiff[3];
+      vdiff[0] = tmp_avg[0];
+      vdiff[1] = tmp_avg[1];
+      vdiff[2] = tmp_avg[2];
+      VALUETYPE f_norm = sqrt(MathUtilities::dot(vdiff, vdiff));
+      // relative std = std/(abs(f)+eps)
+      std[ii] /= f_norm + eps;
+  }
+}
+

source/lmp/pair_nnp.cpp

@@ -450,6 +450,9 @@ void PairNNP::compute(int eflag, int vflag)
 	vector<double> tmp_avg_f;
 	nnp_inter_model_devi.compute_avg (tmp_avg_f, all_force);  
 	nnp_inter_model_devi.compute_std_f (std_f, tmp_avg_f, all_force);
+	if (out_rel == 1){
+	    nnp_inter_model_devi.compute_relative_std_f (std_f, tmp_avg_f, eps);
+	}
 #else 
 	vector<float> tmp_avg_f_, std_f_;
 	for (unsigned ii = 0; ii < all_force_.size(); ++ii){
@@ -461,6 +464,9 @@ void PairNNP::compute(int eflag, int vflag)
 	nnp_inter_model_devi.compute_std_f (std_f_, tmp_avg_f_, all_force_);
 	std_f.resize(std_f_.size());
 	for (int dd = 0; dd < std_f_.size(); ++dd) std_f[dd] = std_f_[dd];
+	if (out_rel == 1){
+	    nnp_inter_model_devi.compute_relative_std_f (std_f, tmp_avg_f_, eps);
+	}
 #endif
 	double min = numeric_limits<double>::max(), max = 0, avg = 0;
 	ana_st(max, min, avg, std_f, nlocal);
@@ -514,16 +520,7 @@ void PairNNP::compute(int eflag, int vflag)
 	      // 1. If the atom_style is not atomic (e.g. charge), the order of std_f is different from that of atom ids.
               // 2. std_f is not gathered by MPI.
 	      for (int dd = 0; dd < all_nlocal; ++dd) {
-          if (out_rel == 1){
-            // relative std = std/(abs(f)+1)
-#ifdef HIGH_PREC
-            fp << " " << setw(18) << std_f[dd] / (fabs(tmp_avg_f[dd]) + eps);
-#else
-            fp << " " << setw(18) << std_f[dd] / (fabsf(tmp_avg_f_[dd]) + eps);
-#endif
-          } else {
             fp << " " << setw(18) << std_f[dd];	
-          }
         }
 	  }
 	  fp << endl;

source/train/CMakeLists.txt

@@ -2,7 +2,7 @@
 
 configure_file("RunOptions.py.in" "${CMAKE_CURRENT_BINARY_DIR}/RunOptions.py" @ONLY)
 
-file(GLOB LIB_PY main.py common.py env.py compat.py Network.py Deep*.py Data.py DataSystem.py Model*.py Descrpt*.py Fitting.py Loss.py LearningRate.py Trainer.py TabInter.py ${CMAKE_CURRENT_BINARY_DIR}/RunOptions.py)
+file(GLOB LIB_PY main.py common.py env.py compat.py calculator.py Network.py Deep*.py Data.py DataSystem.py Model*.py Descrpt*.py Fitting.py Loss.py LearningRate.py Trainer.py TabInter.py ${CMAKE_CURRENT_BINARY_DIR}/RunOptions.py)
 
 file(GLOB CLS_PY  Local.py Slurm.py)
 

source/train/DescrptSeAR.py

@@ -72,7 +72,7 @@ def build (self,
                reuse = None):
         # dout
         self.dout_a = self.descrpt_a.build(coord_, atype_, natoms, box, mesh, davg[0], dstd[0], suffix=suffix+'_a', reuse=reuse)
-        self.dout_r = self.descrpt_r.build(coord_, atype_, natoms, box, mesh, davg[1], dstd[1], suffix=suffix+'_r', reuse=reuse)
+        self.dout_r = self.descrpt_r.build(coord_, atype_, natoms, box, mesh, davg[1], dstd[1], suffix=suffix,      reuse=reuse)
         self.dout_a = tf.reshape(self.dout_a, [-1, self.descrpt_a.get_dim_out()])
         self.dout_r = tf.reshape(self.dout_r, [-1, self.descrpt_r.get_dim_out()])
         self.dout = tf.concat([self.dout_a, self.dout_r], axis = 1)

source/train/calculator.py

@@ -0,0 +1,49 @@
+"""An ASE calculator interface.
+
+Example:
+```python
+from ase import Atoms
+from deepmd.calculator import DP
+
+water = Atoms('H2O',
+              positions=[(0.7601, 1.9270, 1),
+                         (1.9575, 1, 1),
+                         (1., 1., 1.)],
+              cell=[100, 100, 100],
+              calculator=DP(model="frozen_model.pb"))
+print(water.get_potential_energy())
+print(water.get_forces())
+```
+
+Optimization is also available:
+```python
+from ase.optimize import BFGS
+dyn = BFGS(water)
+dyn.run(fmax=1e-6)
+print(water.get_positions())
+```
+"""
+
+from ase.calculators.calculator import Calculator, all_changes
+import deepmd.DeepPot as DeepPot
+
+
+class DP(Calculator):
+    name = "DP"
+    implemented_properties = ["energy", "forces", "stress"]
+
+    def __init__(self, model, label="DP", **kwargs):
+        Calculator.__init__(self, label=label, **kwargs)
+        self.dp = DeepPot(model)
+        self.type_dict = dict(zip(self.dp.get_type_map(), range(self.dp.get_ntypes())))
+
+    def calculate(self, atoms=None, properties=["energy", "forces", "stress"], system_changes=all_changes):
+        coord = atoms.get_positions().reshape([1, -1])
+        cell = atoms.get_cell().reshape([1, -1])
+        symbols = atoms.get_chemical_symbols()
+        atype = [self.type_dict[k] for k in symbols]
+        e, f, v = self.dp.eval(coord, cell, atype)
+        self.results['energy'] = e[0]
+        self.results['forces'] = f[0]
+        self.results['stress'] = v[0]
+

source/train/common.py

@@ -1,3 +1,5 @@
+import warnings
+
 data_requirement = {}
 
 def add_data_requirement(key,