fix an error in stress by ase interface

njzjz · Jan 29, 2021 · a24971f · a24971f
1 parent 1a6cd1c
commit a24971f
Showing 1 changed file with 17 additions and 5 deletions.
diff --git a/source/train/calculator.py b/source/train/calculator.py
@@ -24,13 +24,15 @@
 ```
 """
 
-from ase.calculators.calculator import Calculator, all_changes
+from ase.calculators.calculator import (
+    Calculator, all_changes, PropertyNotImplementedError
+)
 import deepmd.DeepPot as DeepPot
 
 
 class DP(Calculator):
     name = "DP"
-    implemented_properties = ["energy", "forces", "stress"]
+    implemented_properties = ["energy", "forces", "virial", "stress"]
 
     def __init__(self, model, label="DP", type_dict=None, **kwargs):
         Calculator.__init__(self, label=label, **kwargs)
@@ -40,7 +42,7 @@ def __init__(self, model, label="DP", type_dict=None, **kwargs):
         else:
             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):
+    def calculate(self, atoms=None, properties=["energy", "forces", "virial"], system_changes=all_changes):
         coord = atoms.get_positions().reshape([1, -1])
         if sum(atoms.get_pbc())>0:
            cell = atoms.get_cell().reshape([1, -1])
@@ -49,7 +51,17 @@ def calculate(self, atoms=None, properties=["energy", "forces", "stress"], syste
         symbols = atoms.get_chemical_symbols()
         atype = [self.type_dict[k] for k in symbols]
         e, f, v = self.dp.eval(coords=coord, cells=cell, atom_types=atype)
-        self.results['energy'] = e[0]
+        self.results['energy'] = e[0][0]
         self.results['forces'] = f[0]
-        self.results['stress'] = v[0]
+        self.results['virial'] = v[0].reshape(3,3)
 
+        # convert virial into stress for lattice relaxation
+        if "stress" in properties:
+            if sum(atoms.get_pbc()) > 0:
+                # the usual convention (tensile stress is positive)
+                # stress = -virial / volume
+                stress = -0.5*(v[0].copy()+v[0].copy().T) / atoms.get_volume()
+                # Voigt notation 
+                self.results['stress'] = stress.flat[[0,4,8,5,2,1]] 
+            else:
+                raise PropertyNotImplementedError