Read xyz from orca input

Author: briling

qstack/compound.py

@@ -18,15 +18,17 @@ def xyz_to_mol(fin, basis="def2-svp", charge=0, spin=0, ignore=False, unit='ANG'
         basis (str or dict): Basis set.
         charge (int): Charge of the molecule.
         spin (int): Spin of the molecule (alpha electrons - beta electrons).
+        ignore (bool): If assume molecule closed-shell an assign charge either 0 or -1
+        unit (str): units (Ang or Bohr)
+        ecp (str) : ECP to use
 
     Returns:
         A pyscf Mole object containing the molecule information.
     """
 
     # Open and read the file
-    f = open(fin, "r")
-    molxyz = "\n".join(f.read().split("\n")[2:])
-    f.close()
+    with open(fin, "r") as f:
+        molxyz = "\n".join(f.read().split("\n")[2:])
 
     # Define attributes to the Mole object and build it
     mol = gto.Mole()

qstack/orcaio.py

@@ -6,6 +6,46 @@
 from qstack.tools import reorder_ao
 
 
+def read_input(fname, basis, ecp=None):
+    """Read the structure from an Orca input (XYZ coordinates in simple format only)
+
+    Note: we do not read basis set info from the file.
+    TODO: read also %coords block?
+
+    Args:
+        fname (str) : path to file
+        basis (str/dict) : basis name, path to file, or dict in the pyscf format
+    Kwargs:
+        ecp (str) : ECP to use
+
+    Returns:
+        pyscf Mole object.
+    """
+
+    with open(fname, "r") as f:
+        lines = [x.strip() for x in f.readlines()]
+
+    command_line = '\n'.join(y[1:] for y in filter(lambda x: x.startswith('!'), lines)).lower().split()
+    if 'bohrs' in command_line:
+        unit = 'Bohr'
+    else:
+        unit = 'Angstrom'
+
+    idx_xyz_0, idx_xyz_1 = [y[0] for y in filter(lambda x: x[1].startswith('*'), enumerate(lines))][:2]
+    charge, mult = map(int, lines[idx_xyz_0][1:].split()[1:])
+    molxyz = '\n'.join(lines[idx_xyz_0+1:idx_xyz_1])
+
+    mol = pyscf.gto.Mole()
+    mol.atom = molxyz
+    mol.charge = charge
+    mol.spin = mult-1
+    mol.unit = unit
+    mol.basis = basis
+    mol.ecp = ecp
+    mol.build()
+    return mol
+
+
 def read_density(mol, basename, directory='./', version=500, openshell=False, reorder_dest='pyscf'):
     """Read densities from an ORCA output.
 

tests/test_orca.py

@@ -75,7 +75,18 @@ def test_orca_gbw_reader_def2tzvp():
     assert np.linalg.norm(mol_dip-mol_dip_true) < 1e-5
 
 
+def test_orca_input_reader():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol0 = qstack.compound.xyz_to_mol(path+'/data/orca/H2O.xyz', 'sto3g', charge=1, spin=1)
+    mol = qstack.orcaio.read_input(path+'/data/orca/H2O.orca504.inp', 'sto3g')
+    assert mol.natm == mol0.natm
+    assert mol.nelectron == mol0.nelectron
+    assert np.all(mol.elements==mol0.elements)
+    assert np.allclose(mol.atom_coords(), mol0.atom_coords())
+
+
 if __name__ == '__main__':
+    test_orca_input_reader()
     test_orca_density_reader()
     test_orca_gbw_reader()
     test_orca_gbw_reader_def2tzvp()