diff --git a/iodata/formats/mol.py b/iodata/formats/mol.py new file mode 100644 index 000000000..76a56e5a3 --- /dev/null +++ b/iodata/formats/mol.py @@ -0,0 +1,120 @@ +# IODATA is an input and output module for quantum chemistry. +# Copyright (C) 2011-2019 The IODATA Development Team +# +# This file is part of IODATA. +# +# IODATA is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# IODATA is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# -- +"""MOL File Format. + +This script handles either individual molecules or a molecular database +originating from V2000 compliant *.mol molecules. + +More information on *.mol files can be found here +https://en.wikipedia.org/wiki/Chemical_table_file +""" + +from typing import TextIO, Iterator + +import numpy as np + +from ..docstrings import (document_load_one, document_load_many) +from ..docstrings import (document_dump_one, document_dump_many) +from ..iodata import IOData +from ..periodic import sym2num, num2sym +from ..utils import angstrom, LineIterator + + +__all__ = [] + +PATTERNS = ['*.mol'] + + +@document_load_one("MOL", ['atcoords', 'atnums', 'bonds', 'title']) +def load_one(lit: LineIterator) -> dict: + """Do not edit this docstring. It will be overwritten.""" + title = next(lit).strip() # mol title row 1 + next(lit) # mol timestamp row 2 + next(lit) # mol blank line row 3 + words = next(lit).split() + if words[-1].upper() != "V2000": + lit.error("Only V2000 SDF files are supported.") + natom = int(words[0]) # number atoms col0 row 4 + nbond = int(words[1]) # number bonds from col1 row4 + atcoords = np.empty((natom, 3), float) # xyz coords with shape 3 + atnums = np.empty(natom, int) # each atom as iterated + bonds = np.empty((nbond, 3), int) # each bond as iterated + for iatom in range(natom): # iterate entire atom block + words = next(lit).split() # split each atom line + atcoords[iatom, 0] = float(words[0]) * angstrom # coord angstrom + atcoords[iatom, 1] = float(words[1]) * angstrom # coord angstrom + atcoords[iatom, 2] = float(words[2]) * angstrom # coord angstrom + atnums[iatom] = sym2num.get(words[3].title()) # atomic symbol + for ibond in range(nbond): # iterate entire bond block + words = next(lit).split() # each element in bond line + bonds[ibond, 0] = int(words[0]) - 1 # first atom for bond + bonds[ibond, 1] = int(words[1]) - 1 # second atom for bond + bonds[ibond, 2] = int(words[2]) # bond type (single,double,triple) + while True: # iterate each commentary until EOF + try: + words = next(lit) # skip each commentary line + except StopIteration: # what if file has wrong termination? + lit.error("WARNING: Molecule specification did not end properly with M END !") + if words == "M END\n": # what if file has right termination? + break # terminate while loop + return{ + 'title': title, + 'atcoords': atcoords, + 'atnums': atnums, + 'bonds': bonds} + + +@document_load_many("MOL", ['atcoords', 'atnums', 'bonds', 'title']) +def load_many(lit: LineIterator) -> Iterator[dict]: + """Do not edit this docstring. It will be overwritten.""" + while True: + try: + yield load_one(lit) + except StopIteration: + return + + +@document_dump_one("MOL", ['atcoords', 'atnums', 'bonds', 'title']) +def dump_one(f: TextIO, data: IOData): + """Do not edit this docstring. It will be overwritten.""" + print(data.title or 'Created with IOData', file=f) # row 1, title + print('', file=f) # row 2, empty + print('', file=f) # row 3, empty + if data.bonds is None: # What if no bond block data? + nbond = 0 + else: # What if bond block data? + nbond = len(data.bonds) # num bonds is length of bonds array + # populate line 4 with atom number and bond number from file + print("{:3d}{:3d} 0 0 0 0 0 0 0999 V2000".format(data.natom, nbond), file=f) + for beginatom in range(data.natom): # populate entire atom block + n = num2sym[data.atnums[beginatom]] + x, y, z = data.atcoords[beginatom] / angstrom # convert to arb units + print(f'{x:10.4f}{y:10.4f}{z:10.4f} {n:<3s} 0 0 0 0 0 0 0 0 0 0 0 0', file=f) + if data.bonds is not None: # populate bond block + for beginatom, endatom, bondtype in data.bonds: + print('{:3d}{:3d}{:3d} 0 0 0 0'.format(beginatom + 1, endatom + 1, + bondtype), file=f) + print("M END", file=f) + + +@document_dump_many("MOL", ['atcoords', 'atnums', 'bonds', 'title']) +def dump_many(f: TextIO, datas: IOData): + """Do not edit this docstring. It will be overwritten.""" + for data in datas: + dump_one(f, data) diff --git a/iodata/test/data/benzoicacid.mol b/iodata/test/data/benzoicacid.mol new file mode 100644 index 000000000..35689506f --- /dev/null +++ b/iodata/test/data/benzoicacid.mol @@ -0,0 +1,35 @@ +C7H6O2 +APtclcactv06032219373D 0 0.00000 0.00000 + + 15 15 0 0 0 0 0 0 0 0999 V2000 + 1.6471 0.0768 -0.0012 C 0 0 0 0 0 0 0 0 0 0 0 0 + 0.1706 0.0297 -0.0002 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.5689 1.2141 0.0004 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.9473 1.1628 0.0009 C 0 0 0 0 0 0 0 0 0 0 0 0 + -2.5974 -0.0586 0.0013 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.8708 -1.2362 0.0011 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.4918 -1.1994 -0.0057 C 0 0 0 0 0 0 0 0 0 0 0 0 + 2.2206 1.1476 -0.0015 O 0 0 0 0 0 0 0 0 0 0 0 0 + 2.3575 -1.0677 0.0038 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.0627 2.1682 0.0009 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.5205 2.0782 0.0014 H 0 0 0 0 0 0 0 0 0 0 0 0 + -3.6768 -0.0931 0.0020 H 0 0 0 0 0 0 0 0 0 0 0 0 + -2.3844 -2.1862 0.0021 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.0741 -2.1193 -0.0061 H 0 0 0 0 0 0 0 0 0 0 0 0 + 3.3211 -0.9864 0.0029 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 2 3 2 0 0 0 0 + 3 4 1 0 0 0 0 + 4 5 2 0 0 0 0 + 5 6 1 0 0 0 0 + 6 7 2 0 0 0 0 + 2 7 1 0 0 0 0 + 1 8 2 0 0 0 0 + 1 9 1 0 0 0 0 + 3 10 1 0 0 0 0 + 4 11 1 0 0 0 0 + 5 12 1 0 0 0 0 + 6 13 1 0 0 0 0 + 7 14 1 0 0 0 0 + 9 15 1 0 0 0 0 +M END diff --git a/iodata/test/data/database.mol b/iodata/test/data/database.mol new file mode 100644 index 000000000..4fed70487 --- /dev/null +++ b/iodata/test/data/database.mol @@ -0,0 +1,20 @@ +H2O +APtclcactv06032219383D 0 0.00000 0.00000 + + 3 2 0 0 0 0 0 0 0 0999 V2000 + -0.0000 -0.0589 -0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 1 3 1 0 0 0 0 +M END +H2O +APtclcactv06032219383D 0 0.00000 0.00000 + + 3 2 0 0 0 0 0 0 0 0999 V2000 + -0.0000 -0.0589 -0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 1 3 1 0 0 0 0 +M END diff --git a/iodata/test/data/water.mol b/iodata/test/data/water.mol new file mode 100644 index 000000000..2cb441342 --- /dev/null +++ b/iodata/test/data/water.mol @@ -0,0 +1,10 @@ +H2O +APtclcactv06032219383D 0 0.00000 0.00000 + + 3 2 0 0 0 0 0 0 0 0999 V2000 + -0.0000 -0.0589 -0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 1 3 1 0 0 0 0 +M END diff --git a/iodata/test/data/waterv3000.mol b/iodata/test/data/waterv3000.mol new file mode 100644 index 000000000..e043c5aa1 --- /dev/null +++ b/iodata/test/data/waterv3000.mol @@ -0,0 +1,10 @@ +H2O +APtclcactv06032219383D 0 0.00000 0.00000 + + 3 2 0 0 0 0 0 0 0 0999 V3000 + -0.0000 -0.0589 -0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 0.8110 0.4677 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 0 0 0 + 1 3 1 0 0 0 0 +M END diff --git a/iodata/test/test_mol.py b/iodata/test/test_mol.py new file mode 100644 index 000000000..4d7976eee --- /dev/null +++ b/iodata/test/test_mol.py @@ -0,0 +1,160 @@ +# IODATA is an input and output module for quantum chemistry. +# Copyright (C) 2011-2019 The IODATA Development Team +# +# This file is part of IODATA. +# +# IODATA is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# IODATA is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# -- +"""Test iodata.formats.mol module.""" + +import os + +import pytest +from numpy.testing import assert_equal, assert_allclose + +from .common import truncated_file +from ..api import load_one, load_many, dump_one, dump_many +from ..utils import angstrom, FileFormatError + +try: + from importlib_resources import path +except ImportError: + from importlib.resources import path + + +def test_mol_load_one_water(): + # test mol one structure + with path('iodata.test.data', 'water.mol') as fn_mol: + mol = load_one(str(fn_mol)) + assert mol.title == 'H2O' + assert mol.natom == 3 + assert len(mol.bonds) == 2 + assert_equal(mol.atnums, [8, 1, 1]) + # check coordinates + atcoords_ang = mol.atcoords / angstrom + assert_allclose(atcoords_ang[0], [-0.0000, -0.0589, -0.0000]) + assert_allclose(atcoords_ang[1], [-0.8110, 0.4677, 0.0000]) + assert_allclose(atcoords_ang[2], [0.8110, 0.4677, 0.0000]) + assert_equal(mol.bonds[0], [0, 1, 1]) + assert_equal(mol.bonds[1], [0, 2, 1]) + + +def test_mol_load_one_benzoicacid(): + # test mol one structure + with path('iodata.test.data', 'benzoicacid.mol') as fn_mol: + mol = load_one(str(fn_mol)) + assert mol.title == 'C7H6O2' + assert mol.natom == 15 + assert len(mol.bonds) == 15 + assert_equal(mol.atnums, [6, 6, 6, 6, 6, 6, 6, 8, 8, 1, 1, 1, 1, 1, 1]) + atcoords_ang = mol.atcoords / angstrom + # test atoms + assert_allclose(atcoords_ang[0], [1.6471, 0.0768, -0.0012]) + assert_allclose(atcoords_ang[1], [0.1706, 0.0297, -0.0002]) + assert_allclose(atcoords_ang[2], [-0.5689, 1.2141, 0.0004]) + assert_allclose(atcoords_ang[3], [-1.9473, 1.1628, 0.0009]) + assert_allclose(atcoords_ang[4], [-2.5974, -0.0586, 0.0013]) + assert_allclose(atcoords_ang[5], [-1.8708, -1.2362, 0.0011]) + assert_allclose(atcoords_ang[6], [-0.4918, -1.1994, -0.0057]) + assert_allclose(atcoords_ang[7], [2.2206, 1.1476, -0.0015]) + assert_allclose(atcoords_ang[8], [2.3575, -1.0677, 0.0038]) + assert_allclose(atcoords_ang[9], [-0.0627, 2.1682, 0.0009]) + assert_allclose(atcoords_ang[10], [-2.5205, 2.0782, 0.0014]) + assert_allclose(atcoords_ang[11], [-3.6768, -0.0931, 0.0020]) + assert_allclose(atcoords_ang[12], [-2.3844, -2.1862, 0.0021]) + assert_allclose(atcoords_ang[13], [0.0741, -2.1193, -0.0061]) + assert_allclose(atcoords_ang[14], [3.3211, -0.9864, 0.0029]) + # test bonds + assert_equal(mol.bonds[0], [0, 1, 1]) + assert_equal(mol.bonds[1], [1, 2, 2]) + assert_equal(mol.bonds[2], [2, 3, 1]) + assert_equal(mol.bonds[3], [3, 4, 2]) + assert_equal(mol.bonds[4], [4, 5, 1]) + assert_equal(mol.bonds[5], [5, 6, 2]) + assert_equal(mol.bonds[6], [1, 6, 1]) + assert_equal(mol.bonds[7], [0, 7, 2]) + assert_equal(mol.bonds[8], [0, 8, 1]) + assert_equal(mol.bonds[9], [2, 9, 1]) + assert_equal(mol.bonds[10], [3, 10, 1]) + assert_equal(mol.bonds[11], [4, 11, 1]) + assert_equal(mol.bonds[12], [5, 12, 1]) + assert_equal(mol.bonds[13], [6, 13, 1]) + assert_equal(mol.bonds[14], [8, 14, 1]) + + +def test_mol_formaterror(tmpdir): + # test mol file if it has an ending other than 'M END' + with path('iodata.test.data', 'benzoicacid.mol') as fn_test: + with truncated_file(fn_test, 15, 0, tmpdir) as fn: + with pytest.raises(IOError): + load_one(str(fn)) + + +def check_load_dump_consistency(tmpdir, fn): + """Check if dumping and loading an MOL file results in the same data.""" + mol0 = load_one(str(fn)) + fn_tmp = os.path.join(tmpdir, 'test.mol') + dump_one(mol0, fn_tmp) + mol1 = load_one(fn_tmp) + assert mol0.title == mol1.title + assert_equal(mol0.atnums, mol1.atnums) + assert_allclose(mol0.atcoords, mol1.atcoords, atol=1.e-5) + assert_equal(mol0.bonds, mol1.bonds) + + +def test_dump_consistency(tmpdir): + with path('iodata.test.data', 'water.mol') as fn_mol: + check_load_dump_consistency(tmpdir, fn_mol) + with path('iodata.test.data', 'benzoicacid.mol') as fn_mol: + check_load_dump_consistency(tmpdir, fn_mol) + with path('iodata.test.data', 'benzene.mol2') as fn_mol: + check_load_dump_consistency(tmpdir, fn_mol) + + +def test_load_many(): + with path('iodata.test.data', 'database.mol') as fn_mol: + mols = list(load_many(str(fn_mol))) + assert len(mols) == 2 + assert mols[1].title == 'H2O' + assert mols[0].title == 'H2O' + assert_equal(mols[0].bonds[0], [0, 1, 1]) + assert_equal(mols[1].bonds[0], [0, 1, 1]) + assert_equal(mols[0].bonds[1], [0, 2, 1]) + assert_equal(mols[1].bonds[1], [0, 2, 1]) + assert_equal(mols[0].natom, mols[1].natom) + assert_equal(mols[1].natom, 3) + assert_allclose(mols[0].atcoords[0] / angstrom, [-0.0000, -0.0589, -0.0000]) + assert_allclose(mols[1].atcoords[0] / angstrom, [-0.0000, -0.0589, -0.0000]) + assert_allclose(mols[0].atcoords[1] / angstrom, [-0.8110, 0.4677, 0.0000]) + assert_allclose(mols[1].atcoords[1] / angstrom, [-0.8110, 0.4677, 0.0000]) + + +def test_load_dump_many_consistency(tmpdir): + with path('iodata.test.data', 'water.mol') as fn_mol: + mols0 = list(load_many(str(fn_mol))) + fn_tmp = os.path.join(tmpdir, 'test') + dump_many(mols0, fn_tmp, fmt='mol') + mols1 = list(load_many(fn_tmp, fmt='mol')) + assert len(mols0) == len(mols1) + for mol0, mol1 in zip(mols0, mols1): + assert mol0.title == mol1.title + assert_equal(mol0.atnums, mol1.atnums) + assert_allclose(mol0.atcoords, mol1.atcoords, atol=1.e-5) + assert_equal(mol0.bonds, mol1.bonds) + + +def test_v2000_check(): + with path('iodata.test.data', 'waterv3000.mol') as fn_mol: + with pytest.raises(FileFormatError): + load_one(fn_mol)