fix psi4 length conversion (#571)

The length conversion in Psi4 is not fixed to Bohr. Instead, the output
unit follows the input unit and is written in the output file.

---------

Signed-off-by: Jinzhe Zeng <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: njzjz

dpdata/plugins/psi4.py

@@ -3,9 +3,8 @@
 from dpdata.format import Format
 from dpdata.psi4.input import write_psi4_input
 from dpdata.psi4.output import read_psi4_output
-from dpdata.unit import EnergyConversion, ForceConversion, LengthConversion
+from dpdata.unit import EnergyConversion, ForceConversion
 
-length_convert = LengthConversion("bohr", "angstrom").value()
 energy_convert = EnergyConversion("hartree", "eV").value()
 force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
 
@@ -44,7 +43,7 @@ def from_labeled_system(self, file_name: str, **kwargs) -> dict:
             "atom_types": atom_types,
             "atom_names": list(atom_names),
             "atom_numbs": list(atom_numbs),
-            "coords": (coord * length_convert).reshape((1, natoms, 3)),
+            "coords": (coord).reshape((1, natoms, 3)),
             "energies": np.array([energy * energy_convert]),
             "forces": (forces * force_convert).reshape((1, natoms, 3)),
             "cells": np.zeros((1, 3, 3)),

dpdata/psi4/output.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 
+from dpdata.unit import LengthConversion
+
 
 def read_psi4_output(fn: str) -> Tuple[str, np.ndarray, float, np.ndarray]:
     """Read from Psi4 output.
@@ -28,6 +30,7 @@ def read_psi4_output(fn: str) -> Tuple[str, np.ndarray, float, np.ndarray]:
     symbols = None
     forces = None
     energy = None
+    length_unit = None
     with open(fn) as f:
         flag = 0
         for line in f:
@@ -53,14 +56,19 @@ def read_psi4_output(fn: str) -> Tuple[str, np.ndarray, float, np.ndarray]:
                 flag = 1
                 coord = []
                 symbols = []
+            elif line.startswith("    Geometry (in "):
+                # remove ),
+                length_unit = line.split()[2][:-2].lower()
             elif line.startswith("  ## Total Gradient"):
                 flag = 3
                 forces = []
             elif line.startswith("    Total Energy ="):
                 energy = float(line.split()[-1])
+    assert length_unit is not None
+    length_convert = LengthConversion(length_unit, "angstrom").value()
     symbols = np.array(symbols)
     forces = -np.array(forces)
-    coord = np.array(coord)
+    coord = np.array(coord) * length_convert
     assert coord.shape == forces.shape
 
     return symbols, coord, energy, forces