support for reading in xyz files (#151)

* support for reading in xyz files

tangelo/toolboxes/molecular_computation/molecule.py

@@ -18,8 +18,9 @@
 
 import copy
 from dataclasses import dataclass, field
+
 import numpy as np
-from pyscf import gto, scf, ao2mo, symm
+from pyscf import gto, scf, ao2mo, symm, lib
 import openfermion
 import openfermion.ops.representations as reps
 from openfermion.chem.molecular_data import spinorb_from_spatial
@@ -93,7 +94,6 @@ class Molecule:
     n_min_orbitals: int = field(init=False)
 
     def __post_init__(self):
-        self.xyz = atom_string_to_list(self.xyz) if isinstance(self.xyz, str) else self.xyz
         mol = self.to_pyscf(basis="sto-3g", symmetry=False)
         self.n_atoms = mol.natm
         self.n_electrons = mol.nelectron
@@ -117,16 +117,29 @@ def to_pyscf(self, basis="sto-3g", symmetry=False):
             pyscf.gto.Mole: PySCF compatible object.
         """
 
-        mol = gto.Mole()
-        mol.atom = self.xyz
+        mol = gto.Mole(atom=self.xyz)
         mol.basis = basis
         mol.charge = self.q
         mol.spin = self.spin
         mol.symmetry = symmetry
         mol.build()
 
+        self.xyz = list()
+        for sym, xyz in mol._atom:
+            self.xyz += [tuple([sym, tuple([x*lib.parameters.BOHR for x in xyz])])]
+
         return mol
 
+    def to_file(self, filename, format=None):
+        """Write molecule geometry to filename in specified format
+
+        Args:
+            filename (str): The name of the file to output the geometry.
+            format (str): The output type of "raw", "xyz", or "zmat". If None, will be inferred by the filename
+        """
+        mol = self.to_pyscf(basis="sto-3g")
+        mol.tofile(filename, format)
+
     def to_openfermion(self, basis="sto-3g"):
         """Method to return a openfermion.MolecularData object.
 

tangelo/toolboxes/molecular_computation/tests/data/h2.xyz

@@ -0,0 +1,4 @@
+2
+test h2
+H 0.0 0.0 0.0
+H 0.0 0.0 0.7414

tangelo/toolboxes/molecular_computation/tests/test_molecule.py

@@ -13,6 +13,7 @@
 # limitations under the License.
 
 import unittest
+import os
 
 from tangelo import SecondQuantizedMolecule
 from tangelo.molecule_library import mol_H2_sto3g, xyz_H2O
@@ -24,12 +25,23 @@
 H       0.0        0.0        0.0
 H       0.0        0.0        0.7414
 """
+H2_file_xyz = os.path.dirname(os.path.abspath(__file__))+"/data/h2.xyz"
 
 H2O_list = [("O", (0., 0., 0.11779)),
             ("H", (0., 0.75545, -0.47116)),
             ("H", (0., -0.75545, -0.47116))]
 
 
+def atom_list_close(atom1, atom2, atol):
+
+    for (a0, xyz0), (a1, xyz1) in zip(atom1, atom2):
+        if a0 != a1:
+            raise AssertionError(f"Atoms are not the same {a0} != {a1}")
+        for x0, x1 in zip(xyz0, xyz1):
+            if abs(x0 - x1) > atol:
+                raise AssertionError(f"geometries for atom {a0} are different. {x0} != {x1}")
+
+
 class CoordsTest(unittest.TestCase):
 
     def test_atoms_string_to_list(self):
@@ -52,6 +64,14 @@ def test_instantiate_H2(self):
         assert(molecule.q == 0)
         assert(molecule.n_electrons == 2)
 
+        molecule = SecondQuantizedMolecule(H2_file_xyz, 0, 0, "sto-3g")
+        assert(molecule.elements == ["H"]*2)
+        assert(molecule.basis == "sto-3g")
+        assert(molecule.spin == 0)
+        assert(molecule.q == 0)
+        assert(molecule.n_electrons == 2)
+        atom_list_close(molecule.xyz, H2_list, 1.e-14)
+
     def test_freezing_orbitals(self):
         """Verify freezing orbitals functionalities."""