diff --git a/chemtools/__init__.py b/chemtools/__init__.py index fa23ccdc..fa68f85d 100644 --- a/chemtools/__init__.py +++ b/chemtools/__init__.py @@ -23,17 +23,21 @@ # pragma pylint: disable=wildcard-import """The Main ChemTools Package.""" +import sys -from chemtools.wrappers import * from chemtools.toolbox import * from chemtools.conceptual import * from chemtools.denstools import * from chemtools.utils import * from chemtools.outputs import * -import horton +if sys.version_info.major == 2: + from chemtools.wrappers2 import * + import horton + horton.log.head_banner = "" + horton.log.foot_banner = "" +else: + from chemtools.wrappers3 import * -horton.log.head_banner = "" -horton.log.foot_banner = "" __version__ = '0.9.0' diff --git a/chemtools/conceptual/base.py b/chemtools/conceptual/base.py index f2b9b707..a68db84f 100644 --- a/chemtools/conceptual/base.py +++ b/chemtools/conceptual/base.py @@ -402,13 +402,13 @@ def grand_potential_derivative(self, n_elec, order=1): else: # higher-order derivatives are compute with Faa Di Bruno formula # list of hyper-hardneses (derivatives of energy w.r.t. N) - e_deriv = [self.energy_derivative(n_elec, i + 1) for i in xrange(1, order)] - g_deriv = [self.grand_potential_derivative(n_elec, k + 1) for k in xrange(1, order - 1)] + e_deriv = [self.energy_derivative(n_elec, i + 1) for i in range(1, order)] + g_deriv = [self.grand_potential_derivative(n_elec, k + 1) for k in range(1, order - 1)] if any([item is None for item in e_deriv]) or any([item is None for item in g_deriv]): deriv = None else: deriv = 0 - for k in xrange(1, order - 1): + for k in range(1, order - 1): deriv -= g_deriv[k - 1] * sp.bell(order - 1, k, e_deriv[:order - k]) deriv /= sp.bell(order - 1, order - 1, [e_deriv[0]]) return deriv diff --git a/chemtools/conceptual/general.py b/chemtools/conceptual/general.py index 3ec2ce55..8c626b6f 100644 --- a/chemtools/conceptual/general.py +++ b/chemtools/conceptual/general.py @@ -89,7 +89,7 @@ def __init__(self, expr, n0, n_energies, n_symbol=None, n0_symbol=None, guess=No 'the number of electrons.') self._n_symb = n_symbol # store minimum and maximum number of electrons used for interpolation - self._n_min, self._n_max = np.min(n_energies.keys()), np.max(n_energies.keys()) + self._n_min, self._n_max = np.min(list(n_energies.keys())), np.max(list(n_energies.keys())) # substitute N0 in energy expression if n0_symbol: @@ -184,7 +184,7 @@ def _solve_parameters(self, expr, n_energies, guess, opts=None): # construct system of equations to solve system_eqns = [] d_system_eqns = [] - for n, energy in n_energies.iteritems(): + for n, energy in n_energies.items(): eqn = sp.lambdify((params,), expr.subs(self._n_symb, n) - energy, 'numpy') system_eqns.append(eqn) d_eqn_row = [] @@ -234,7 +234,6 @@ def _solve_nmax(self, guess): d_expr = self._expr.diff(self._n_symb) n_max_eqn = sp.lambdify(self._n_symb, d_expr, 'numpy') result = root(n_max_eqn, guess) - print result if result.success: n_max = np.asscalar(result.x) # n_ceil = math.ceil(n_max) diff --git a/chemtools/conceptual/leastnorm.py b/chemtools/conceptual/leastnorm.py index dc1c802a..f5c86d1b 100644 --- a/chemtools/conceptual/leastnorm.py +++ b/chemtools/conceptual/leastnorm.py @@ -26,7 +26,7 @@ """ import numpy as np -from scipy.misc import factorial +from scipy.special import factorial from chemtools.utils.utils import doc_inherit from chemtools.conceptual.base import BaseGlobalTool diff --git a/chemtools/conceptual/utils.py b/chemtools/conceptual/utils.py index c73561cb..8888a71d 100644 --- a/chemtools/conceptual/utils.py +++ b/chemtools/conceptual/utils.py @@ -40,7 +40,7 @@ def check_dict_values(dict_values): # check length of dictionary & its keys if len(dict_values) != 3 or not all([key >= 0 for key in dict_values.keys()]): raise ValueError("The energy model requires 3 keys corresponding to positive " - "number of electrons! Given keys={0}".format(dict_values.keys())) + "number of electrons! Given keys={0}".format(list(dict_values.keys()))) # find reference number of electrons n_ref = sorted(dict_values.keys())[1] if n_ref < 1: @@ -48,7 +48,7 @@ def check_dict_values(dict_values): # check that number of electrons differ by one if sorted(dict_values.keys()) != [n_ref - 1, n_ref, n_ref + 1]: raise ValueError("In current implementation, the number of electrons (keys) should " - "differ by one! Given keys={0}".format(dict_values.keys())) + "differ by one! Given keys={0}".format(list(dict_values.keys()))) # check that all values have the same type if not all([isinstance(value, type(dict_values[n_ref])) for value in dict_values.values()]): raise ValueError("All values in dictionary should be of the same type!") diff --git a/chemtools/outputs/vmd.py b/chemtools/outputs/vmd.py index ca21b130..bd96f3cd 100644 --- a/chemtools/outputs/vmd.py +++ b/chemtools/outputs/vmd.py @@ -592,7 +592,7 @@ def print_vmd_script_multiple_cube(scriptfile, cubes, isosurfs=None, material='O raise TypeError('Each iso-surface value must be a float') if colors is None: - colors = range(len(cubes)) + colors = list(range(len(cubes))) elif not (isinstance(colors, (list, tuple)) and len(colors) == len(cubes)): raise TypeError('The colors must be provided as a list or tuple of the same length as the ' 'number of cube files') diff --git a/chemtools/scripts/chemtools_mot.py b/chemtools/scripts/chemtools_mot.py index 07a1f8df..519d49fe 100644 --- a/chemtools/scripts/chemtools_mot.py +++ b/chemtools/scripts/chemtools_mot.py @@ -25,12 +25,15 @@ """Molecular Orbital Theory (MOT) Script.""" +import sys import numpy as np -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.motbased import OrbPart from chemtools.scripts.common import help_cube, load_molecule_and_grid +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + description_mot = """ Visualize Molecular Orbitals (MO) using VMD package. diff --git a/chemtools/scripts/common.py b/chemtools/scripts/common.py index c4f41fd1..09a7b3eb 100644 --- a/chemtools/scripts/common.py +++ b/chemtools/scripts/common.py @@ -24,11 +24,14 @@ """Common utility for scripts.""" +import sys import numpy as np -from chemtools.wrappers.molecule import Molecule from chemtools.utils.cube import UniformGrid +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + help_cube = """ cubic grid used for evaluation and visualization. diff --git a/chemtools/toolbox/conceptual.py b/chemtools/toolbox/conceptual.py index d8f3ccd5..6e9127d9 100644 --- a/chemtools/toolbox/conceptual.py +++ b/chemtools/toolbox/conceptual.py @@ -27,11 +27,9 @@ compute various conceptual density functional theory (DFT) descriptive tools. """ - +import sys import logging -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.grid import MolecularGrid from chemtools.toolbox.utils import check_arg_molecule, get_matching_attr from chemtools.toolbox.utils import get_dict_energy, get_dict_density, get_dict_population from chemtools.conceptual.linear import LinearGlobalTool, LinearLocalTool, LinearCondensedTool @@ -40,6 +38,14 @@ from chemtools.conceptual.exponential import ExponentialGlobalTool from chemtools.conceptual.rational import RationalGlobalTool from chemtools.conceptual.general import GeneralGlobalTool + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.grid import MolecularGrid +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.grid import MolecularGrid + try: from pathlib2 import Path except ImportError: diff --git a/chemtools/toolbox/densbased.py b/chemtools/toolbox/densbased.py index f6c7ff8f..f846c793 100644 --- a/chemtools/toolbox/densbased.py +++ b/chemtools/toolbox/densbased.py @@ -23,10 +23,15 @@ # pragma pylint: disable=invalid-name """Density-Based Local Tools.""" +import sys -from chemtools.wrappers.molecule import Molecule from chemtools.denstools.densbased import DensGradLapKedTool +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + class DensityLocalTool(DensGradLapKedTool): """Density Local Tool Class.""" diff --git a/chemtools/toolbox/dftbased.py b/chemtools/toolbox/dftbased.py index 391fac8f..cc9b6d69 100644 --- a/chemtools/toolbox/dftbased.py +++ b/chemtools/toolbox/dftbased.py @@ -24,11 +24,15 @@ """Density Functional Theory (DFT) Based Tools.""" +import sys import numpy as np from scipy.optimize import bisect -from chemtools.wrappers.molecule import Molecule +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule class DFTBasedTool(object): diff --git a/chemtools/toolbox/interactions.py b/chemtools/toolbox/interactions.py index bcb97bb5..2df28e92 100644 --- a/chemtools/toolbox/interactions.py +++ b/chemtools/toolbox/interactions.py @@ -24,9 +24,9 @@ """Module for (non)bonding interaction analysis of Quantum Chemistry Output Files.""" +import sys import numpy as np -from chemtools.wrappers.molecule import Molecule from chemtools.denstools.densbased import DensGradTool from chemtools.utils.utils import doc_inherit from chemtools.utils.cube import UniformGrid @@ -35,6 +35,11 @@ from numpy.ma import masked_less +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + class BaseInteraction(object): """Base class for (non)bonding interactions indicators.""" diff --git a/chemtools/toolbox/kinetic.py b/chemtools/toolbox/kinetic.py index 99942e3b..4e28b7f8 100644 --- a/chemtools/toolbox/kinetic.py +++ b/chemtools/toolbox/kinetic.py @@ -24,12 +24,17 @@ """Kinetic Energy Density Module.""" +import sys import numpy as np from chemtools.utils.utils import doc_inherit -from chemtools.wrappers.molecule import Molecule from chemtools.denstools.densbased import DensGradTool, DensGradLapTool, DensGradLapKedTool +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + class KED(object): """Kinetic Energy Density Class.""" diff --git a/chemtools/toolbox/motbased.py b/chemtools/toolbox/motbased.py index e6e53ae9..a2ffb98a 100644 --- a/chemtools/toolbox/motbased.py +++ b/chemtools/toolbox/motbased.py @@ -24,11 +24,16 @@ """Orbital-Based Population Analysis.""" +import sys import numpy as np from chemtools.utils.utils import doc_inherit from chemtools.orbstools.partition import OrbitalPartitionTools -from chemtools.wrappers.molecule import Molecule, MolecularOrbitals + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule, MolecularOrbitals +else: + from chemtools.wrappers3.molecule import Molecule, MolecularOrbitals class OrbPart(object): diff --git a/chemtools/toolbox/oxidation.py b/chemtools/toolbox/oxidation.py index 0e6f8842..b738c2e1 100644 --- a/chemtools/toolbox/oxidation.py +++ b/chemtools/toolbox/oxidation.py @@ -23,11 +23,16 @@ """Module for Oxidation State.""" +import sys import itertools import numpy as np -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.part import DensPart +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.part import DensPart +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.part import DensPart class EOS(object): diff --git a/chemtools/toolbox/test/test_conceptual_condensed.py b/chemtools/toolbox/test/test_conceptual_condensed.py index d7ca619c..eb2fc5f7 100644 --- a/chemtools/toolbox/test/test_conceptual_condensed.py +++ b/chemtools/toolbox/test/test_conceptual_condensed.py @@ -24,13 +24,20 @@ """Test chemtools.analysis.conceptual.CondensedConceptualDFT.""" +import sys import numpy as np -from numpy.testing import assert_raises, assert_equal, assert_almost_equal +from numpy.testing import assert_equal, assert_almost_equal -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.grid import MolecularGrid from chemtools.toolbox.conceptual import CondensedConceptualDFT + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.grid import MolecularGrid +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.grid import MolecularGrid + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_conceptual_global.py b/chemtools/toolbox/test/test_conceptual_global.py index c4f8b823..9bd1a45e 100644 --- a/chemtools/toolbox/test/test_conceptual_global.py +++ b/chemtools/toolbox/test/test_conceptual_global.py @@ -24,12 +24,18 @@ """Test chemtools.analysis.conceptual.GlobalConceptualDFT.""" +import sys import numpy as np from numpy.testing import assert_raises, assert_almost_equal from chemtools.toolbox.conceptual import GlobalConceptualDFT -from chemtools.wrappers.molecule import Molecule + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_conceptual_local.py b/chemtools/toolbox/test/test_conceptual_local.py index 23375400..ef6bcbe0 100644 --- a/chemtools/toolbox/test/test_conceptual_local.py +++ b/chemtools/toolbox/test/test_conceptual_local.py @@ -24,13 +24,20 @@ """Test chemtools.analysis.conceptual.LocalConceptualDFT.""" +import sys import numpy as np from numpy.testing import assert_raises, assert_equal, assert_almost_equal -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.grid import MolecularGrid from chemtools.toolbox.conceptual import LocalConceptualDFT + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.grid import MolecularGrid +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.grid import MolecularGrid + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_densbased.py b/chemtools/toolbox/test/test_densbased.py index d4a37bac..61989bb0 100644 --- a/chemtools/toolbox/test/test_densbased.py +++ b/chemtools/toolbox/test/test_densbased.py @@ -24,12 +24,18 @@ """Test chemtools.toolbox.densbased.""" +import sys import numpy as np from numpy.testing import assert_allclose -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.densbased import DensityLocalTool + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_dftbased.py b/chemtools/toolbox/test/test_dftbased.py index b99579d9..978dd650 100644 --- a/chemtools/toolbox/test/test_dftbased.py +++ b/chemtools/toolbox/test/test_dftbased.py @@ -24,12 +24,18 @@ """Test chemtools.toolbox.dftbased.""" +import sys import numpy as np from numpy.testing import assert_raises, assert_array_almost_equal -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.dftbased import DFTBasedTool + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_elf.py b/chemtools/toolbox/test/test_elf.py index 5493515c..2b687c01 100644 --- a/chemtools/toolbox/test/test_elf.py +++ b/chemtools/toolbox/test/test_elf.py @@ -26,6 +26,7 @@ import numpy as np from numpy.testing import assert_allclose, assert_raises from chemtools.toolbox.interactions import ELF, LOL + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_kinetic.py b/chemtools/toolbox/test/test_kinetic.py index e7ae7ccf..645851d0 100644 --- a/chemtools/toolbox/test/test_kinetic.py +++ b/chemtools/toolbox/test/test_kinetic.py @@ -24,12 +24,18 @@ """Test chemtools.toolbox.kinetic.""" +import sys import numpy as np from numpy.testing import assert_array_almost_equal, assert_allclose, assert_raises -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.kinetic import KED + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_nci.py b/chemtools/toolbox/test/test_nci.py index 9977d8ee..1e6ceb8b 100644 --- a/chemtools/toolbox/test/test_nci.py +++ b/chemtools/toolbox/test/test_nci.py @@ -24,6 +24,7 @@ import os +import sys import shutil import tempfile from contextlib import contextmanager @@ -33,7 +34,12 @@ from chemtools.utils import UniformGrid from chemtools.toolbox.interactions import NCI -from chemtools.wrappers.molecule import Molecule + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_oxidation.py b/chemtools/toolbox/test/test_oxidation.py index b809f6ef..78d8c1e4 100644 --- a/chemtools/toolbox/test/test_oxidation.py +++ b/chemtools/toolbox/test/test_oxidation.py @@ -23,16 +23,22 @@ """Test chemtools.toolbox.oxidation.""" +import sys import glob import numpy as np from numpy.testing import assert_equal, assert_almost_equal, assert_approx_equal -from horton import ProAtomDB -from chemtools.wrappers.grid import MolecularGrid -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.part import DensPart from chemtools.toolbox.oxidation import EOS +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.grid import MolecularGrid + from chemtools.wrappers2.part import DensPart + from horton import ProAtomDB +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.grid import MolecularGrid + from chemtools.wrappers3.part import DensPart try: from importlib_resources import path diff --git a/chemtools/toolbox/test/test_partition.py b/chemtools/toolbox/test/test_partition.py index 6cc0b4e5..41ffeb52 100644 --- a/chemtools/toolbox/test/test_partition.py +++ b/chemtools/toolbox/test/test_partition.py @@ -24,14 +24,20 @@ """Test chemtools.toolbox.motbased.OrbPart""" +import sys import numpy as np from numpy.testing import assert_raises, assert_allclose + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: from importlib.resources import path -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.motbased import OrbPart diff --git a/chemtools/toolbox/test/test_topology.py b/chemtools/toolbox/test/test_topology.py index 4781e2d6..ae726f48 100644 --- a/chemtools/toolbox/test/test_topology.py +++ b/chemtools/toolbox/test/test_topology.py @@ -23,12 +23,17 @@ """Test chemtools.toolbox.topology.""" +import sys import numpy as np from chemtools.toolbox.topology import TopologicalTool -from chemtools.wrappers.molecule import Molecule from chemtools.utils.cube import UniformGrid +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/test/test_utils.py b/chemtools/toolbox/test/test_utils.py index 0440c10a..11de23d4 100644 --- a/chemtools/toolbox/test/test_utils.py +++ b/chemtools/toolbox/test/test_utils.py @@ -22,14 +22,21 @@ # -- """Test chemtools.toolbox.utils.""" + +import sys import numpy as np from numpy.testing import assert_raises from chemtools import UniformGrid -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.utils import get_matching_attr, get_molecular_grid from chemtools.toolbox.utils import get_dict_energy, get_dict_density, get_dict_population + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/toolbox/topology.py b/chemtools/toolbox/topology.py index c4cdfb29..1c7ae6de 100644 --- a/chemtools/toolbox/topology.py +++ b/chemtools/toolbox/topology.py @@ -24,13 +24,18 @@ """Module for topological analysis of scalar fields.""" +import sys import numpy as np -from chemtools.wrappers.molecule import Molecule from chemtools.utils.cube import UniformGrid from chemtools.topology.critical import Topology from chemtools.outputs.vmd import print_vmd_script_topology +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + class TopologicalTool(Topology): """Topological analysis of scalar functions.""" diff --git a/chemtools/toolbox/utils.py b/chemtools/toolbox/utils.py index fa6971c5..c8250828 100644 --- a/chemtools/toolbox/utils.py +++ b/chemtools/toolbox/utils.py @@ -24,13 +24,18 @@ """Utility Functions of Toolbox Module.""" +import sys import numpy as np -from horton import ProAtomDB -from horton.scripts.wpart import wpart_schemes -from chemtools.wrappers.grid import MolecularGrid -from chemtools.wrappers.molecule import Molecule +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.grid import MolecularGrid + from horton import ProAtomDB + from horton.scripts.wpart import wpart_schemes +else: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.grid import MolecularGrid def check_arg_molecule(molecule): diff --git a/chemtools/utils/cube.py b/chemtools/utils/cube.py index e789220f..7d5ada70 100644 --- a/chemtools/utils/cube.py +++ b/chemtools/utils/cube.py @@ -23,10 +23,14 @@ """The Cube Module.""" +import sys import logging import numpy as np -from chemtools.wrappers.molecule import Molecule +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule try: from importlib_resources import path diff --git a/chemtools/utils/test/test_cube.py b/chemtools/utils/test/test_cube.py index 66c12df6..e0472a63 100644 --- a/chemtools/utils/test/test_cube.py +++ b/chemtools/utils/test/test_cube.py @@ -23,15 +23,21 @@ """Test chemtools.utils.cube.""" +import sys import shutil import tempfile from contextlib import contextmanager from numpy.testing import assert_raises, assert_allclose import numpy as np -from chemtools.wrappers.molecule import Molecule from chemtools.toolbox.conceptual import LocalConceptualDFT from chemtools.utils.cube import UniformGrid + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule +else: + from chemtools.wrappers3.molecule import Molecule + try: from importlib_resources import path except ImportError: diff --git a/chemtools/utils/test/test_mesh.py b/chemtools/utils/test/test_mesh.py index 11cae19d..604d5fda 100644 --- a/chemtools/utils/test/test_mesh.py +++ b/chemtools/utils/test/test_mesh.py @@ -1,8 +1,11 @@ """Test chemtools.utils.mesh.""" -from chemtools.utils.mesh import mesh_plane + + import numpy as np from numpy.testing import assert_raises +from chemtools.utils.mesh import mesh_plane + def test_plane_mesh(): """Test chemtools.utils.mesh.plane_mesh.""" diff --git a/chemtools/utils/test/test_utils.py b/chemtools/utils/test/test_utils.py index 110b85e0..fca42bb0 100644 --- a/chemtools/utils/test/test_utils.py +++ b/chemtools/utils/test/test_utils.py @@ -22,7 +22,6 @@ # -- """Test chemtools.utils.utils.""" -import os from numpy.testing import assert_equal, assert_raises diff --git a/chemtools/utils/utils.py b/chemtools/utils/utils.py index 3b259e44..37b85f7a 100644 --- a/chemtools/utils/utils.py +++ b/chemtools/utils/utils.py @@ -22,8 +22,6 @@ # -- """The Utility Module.""" -import os -from glob import glob __all__ = ['doc_inherit'] diff --git a/chemtools/wrappers2/__init__.py b/chemtools/wrappers2/__init__.py new file mode 100644 index 00000000..bb0e102b --- /dev/null +++ b/chemtools/wrappers2/__init__.py @@ -0,0 +1,31 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- +"""The Python 2 Wrappers Module.""" + + +import sys + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import * + from chemtools.wrappers2.grid import * + from chemtools.wrappers2.part import * diff --git a/chemtools/wrappers/grid.py b/chemtools/wrappers2/grid.py similarity index 99% rename from chemtools/wrappers/grid.py rename to chemtools/wrappers2/grid.py index 7ca8f2ff..b5d6e226 100644 --- a/chemtools/wrappers/grid.py +++ b/chemtools/wrappers2/grid.py @@ -26,7 +26,7 @@ import numpy as np from horton import BeckeMolGrid -from chemtools.wrappers.molecule import Molecule +from chemtools.wrappers2.molecule import Molecule __all__ = ['MolecularGrid'] diff --git a/chemtools/wrappers/molecule.py b/chemtools/wrappers2/molecule.py similarity index 100% rename from chemtools/wrappers/molecule.py rename to chemtools/wrappers2/molecule.py diff --git a/chemtools/wrappers/part.py b/chemtools/wrappers2/part.py similarity index 98% rename from chemtools/wrappers/part.py rename to chemtools/wrappers2/part.py index 63c735fd..00329a4c 100644 --- a/chemtools/wrappers/part.py +++ b/chemtools/wrappers2/part.py @@ -28,8 +28,8 @@ from horton import ProAtomDB from horton.scripts.wpart import wpart_schemes -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.grid import MolecularGrid +from chemtools.wrappers2.molecule import Molecule +from chemtools.wrappers2.grid import MolecularGrid __all__ = ['DensPart'] diff --git a/chemtools/wrappers/__init__.py b/chemtools/wrappers2/test/__init__.py similarity index 85% rename from chemtools/wrappers/__init__.py rename to chemtools/wrappers2/test/__init__.py index a44016cf..ec97ea1e 100644 --- a/chemtools/wrappers/__init__.py +++ b/chemtools/wrappers2/test/__init__.py @@ -20,9 +20,3 @@ # along with this program; if not, see # # -- -"""The Wrappers Module.""" - - -from chemtools.wrappers.molecule import * -from chemtools.wrappers.grid import * -from chemtools.wrappers.part import * diff --git a/chemtools/wrappers2/test/test_grid.py b/chemtools/wrappers2/test/test_grid.py new file mode 100644 index 00000000..2531153b --- /dev/null +++ b/chemtools/wrappers2/test/test_grid.py @@ -0,0 +1,81 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 chemtools.wrappers2.grid.""" + + +import sys +import pytest +import numpy as np + +from numpy.testing import assert_raises, assert_allclose + +pytestmark = pytest.mark.skipif(sys.version_info.major != 2, reason="Requires Python 2.") + +if sys.version_info.major == 2: + from chemtools.wrappers2.grid import MolecularGrid + from chemtools.wrappers2.molecule import Molecule + +try: + from importlib_resources import path +except ImportError: + from importlib.resources import path + + +def test_wrapper_grid_raises(): + with path('chemtools.data', 'ch4_uhf_ccpvdz.fchk') as fpath: + assert_raises(TypeError, MolecularGrid.from_molecule, fpath, 'exp:1e-5:20:40:50') + assert_raises(ValueError, MolecularGrid.from_file, fpath, 'ex:1e-5:20:40:50') + assert_raises(ValueError, MolecularGrid.from_file, fpath, 'exp:1e-5:-20:40:50') + assert_raises(ValueError, MolecularGrid.from_file, fpath, 'exp:1e-5:20:40:10') + assert_raises(ValueError, MolecularGrid.from_file, fpath, 'pow:1e-5:20:40:50') + assert_raises(ValueError, MolecularGrid.from_file, fpath, 'veryfin') + with path('chemtools.data', 'ch4_uhf_ccpvdz.fchk') as fpath: + grid = MolecularGrid.from_file(fpath) + assert_raises(ValueError, grid.integrate, np.array([[1., 2., 3.]])) + assert_raises(NotImplementedError, grid.compute_spherical_average, None) + + +def test_wrapper_grid_ch4(): + with path('chemtools.data', 'ch4_uhf_ccpvdz.fchk') as fpath: + mol = Molecule.from_file(fpath) + grid = MolecularGrid(mol.coordinates, mol.numbers, mol.pseudo_numbers, 'exp:1e-5:25:80:230') + assert grid.weights.ndim == 1 + assert grid.points.shape[0] == grid.weights.shape[0] + assert grid.points.shape == (grid.npoints, 3) + # check grid basics + assert_allclose(mol.coordinates, grid.coordinates, rtol=0., atol=1.e-7) + assert_allclose(mol.numbers, grid.numbers, rtol=0., atol=1.e-7) + assert_allclose(mol.pseudo_numbers, grid.pseudo_numbers, rtol=0., atol=1.e-7) + # check integrate + assert_allclose(10., grid.integrate(mol.compute_density(grid.points)), rtol=0., atol=1.e-4) + + +def test_wrapper_grid_from_file_o2(): + with path('chemtools.data', 'o2_uhf.wfn') as fpath: + grid = MolecularGrid.from_file(fpath, 'veryfine') + mol = Molecule.from_file(fpath) + assert grid.weights.ndim == 1 + assert grid.points.shape[0] == grid.weights.shape[0] + assert grid.points.shape == (grid.npoints, 3) + # check integrate + assert_allclose(16., grid.integrate(mol.compute_density(grid.points)), rtol=0., atol=1.e-4) diff --git a/chemtools/wrappers/test/test_molecule.py b/chemtools/wrappers2/test/test_molecule.py similarity index 99% rename from chemtools/wrappers/test/test_molecule.py rename to chemtools/wrappers2/test/test_molecule.py index bef44b29..577fa910 100644 --- a/chemtools/wrappers/test/test_molecule.py +++ b/chemtools/wrappers2/test/test_molecule.py @@ -20,12 +20,19 @@ # along with this program; if not, see # # -- -"""Test chemtools.wrappers.molecule.""" +"""Test chemtools.wrappers2.molecule.""" +import sys +import pytest import numpy as np + from numpy.testing import assert_raises, assert_equal, assert_almost_equal -from chemtools.wrappers import Molecule + +pytestmark = pytest.mark.skipif(sys.version_info.major != 2, reason="Requires Python 2.") + +if sys.version_info.major == 2: + from chemtools.wrappers2 import Molecule try: from importlib_resources import path diff --git a/chemtools/wrappers2/test/test_part.py b/chemtools/wrappers2/test/test_part.py new file mode 100644 index 00000000..019008fc --- /dev/null +++ b/chemtools/wrappers2/test/test_part.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- + + +import sys +import pytest +import numpy as np + +pytestmark = pytest.mark.skipif(sys.version_info.major != 2, reason="Requires Python 2.") + +if sys.version_info.major == 2: + from chemtools.wrappers2.molecule import Molecule + from chemtools.wrappers2.part import DensPart + from chemtools.wrappers2.grid import MolecularGrid + +try: + from importlib_resources import path +except ImportError: + from importlib.resources import path + + +def test_condense_linear_from_file_fmr_h_ch4_fchk(): + # expected populations of CH4 computed with HORTON + with path('chemtools.data', 'ch4_uhf_ccpvdz.fchk') as fname: + part = DensPart.from_file(fname, scheme='h') + expected = np.array([6.11301651, 0.97175462, 0.97175263, 0.9717521, 0.97174353]) + computed = part.numbers - part.charges + assert np.all(abs(expected - computed) < 1.e-3) + assert np.all(abs(part.condense_to_atoms(part.density) - computed) < 1.e-2) + + +def test_condense_quadratic_from_molecule_fmr_mbis_ch4_wfn(): + # expected populations of CH4 computed with HORTON + with path('chemtools.data', 'ch4_uhf_ccpvdz.wfn') as fname: + molecule = Molecule.from_file(fname) + part = DensPart.from_molecule(molecule, scheme='mbis') + expected = np.array([6.46038055, 0.88489494, 0.88492901, 0.88493897, 0.88492396]) + computed = part.numbers - part.charges + assert np.all(abs(expected - computed) < 1.e-2) + assert np.all(abs(part.condense_to_atoms(part.density) - computed) < 1.e-2) + + +def test_condense_quadratic_from_molecule_fmr_mbis_ch4(): + # expected populations of CH4 computed with HORTON + with path('chemtools.data', 'ch4_uhf_ccpvdz.wfn') as fname: + molecule = Molecule.from_file(fname) + grid = MolecularGrid.from_molecule(molecule, 'fine') + part = DensPart.from_molecule(molecule, scheme='mbis', grid=grid) + expected = np.array([6.46038055, 0.88489494, 0.88492901, 0.88493897, 0.88492396]) + computed = part.numbers - part.charges + assert np.all(abs(expected - computed) < 1.e-2) + assert np.all(abs(part.condense_to_atoms(part.density) - computed) < 1.e-2) diff --git a/chemtools/wrappers3/__init__.py b/chemtools/wrappers3/__init__.py new file mode 100644 index 00000000..e1b74ebe --- /dev/null +++ b/chemtools/wrappers3/__init__.py @@ -0,0 +1,31 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- +"""The Python 3 Wrappers Module.""" + + +import sys + +if sys.version_info.major != 2: + from chemtools.wrappers3.molecule import * + from chemtools.wrappers3.grid import * + from chemtools.wrappers3.part import * diff --git a/chemtools/wrappers3/grid.py b/chemtools/wrappers3/grid.py new file mode 100644 index 00000000..03c81eff --- /dev/null +++ b/chemtools/wrappers3/grid.py @@ -0,0 +1,205 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- +"""Grid Wrapper Module.""" + + +import numpy as np + +# from horton import BeckeMolGrid +from chemtools.wrappers3.molecule import Molecule + + +__all__ = ['MolecularGrid'] + + +class MolecularGrid(object): + """Becke-Lebedev molecular grid for numerical integrations.""" + + def __init__(self, coordinates, numbers, pseudo_numbers, specs='medium', k=3, rotate=False): + """Initialize class. + + Parameters + ---------- + coordinates : np.ndarray, shape=(M, 3) + Cartesian coordinates of `M` atoms in the molecule. + numbers : np.ndarray, shape=(M,) + Atomic number of `M` atoms in the molecule. + pseudo_numbers : np.ndarray, shape=(M,) + Pseudo-number of `M` atoms in the molecule. + specs : str, optional + Specification of grid. Either choose from ['coarse', 'medium', 'fine', 'veryfine', + 'ultrafine', 'insane'] or provide a string of 'rname:rmin:rmax:nrad:nang' format. + Here 'rname' denotes the type of radial grid and can be chosen from ['linear', 'exp', + 'power'], 'rmin' and 'rmax' specify the first and last radial grid points in angstrom, + 'nrad' specify the number of radial grid points, and 'nang' specify the number of + angular Lebedev-Laikov grid. The 'nang' can be chosen from (6, 14, 26, 38, 50, 74, 86, + 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, + 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810). + k : int, optional + The order of the switching function in Becke's weighting scheme. + rotate : bool, optional + Whether to randomly rotate spherical grids. + + """ + self._coordinates = coordinates + self._numbers = numbers + self._pseudo_numbers = pseudo_numbers + self._k = k + self._rotate = rotate + self.specs = specs + + self._grid = BeckeMolGrid(self.coordinates, self.numbers, self.pseudo_numbers, + agspec=self.specs, k=k, random_rotate=rotate, mode='keep') + + @classmethod + def from_molecule(cls, molecule, specs='medium', k=3, rotate=False): + """Initialize the class given an instance of Molecule. + + Parameters + ---------- + molecule : instance of Molecule + Instance of Molecule class. + specs : str, optional + Specification of grid. Either choose from ['coarse', 'medium', 'fine', 'veryfine', + 'ultrafine', 'insane'] or provide a string of 'rname:rmin:rmax:nrad:nang' format. + Here 'rname' denotes the type of radial grid and can be chosen from ['linear', 'exp', + 'power'], 'rmin' and 'rmax' specify the first and last radial grid points in angstrom, + 'nrad' specify the number of radial grid points, and 'nang' specify the number of + angular Lebedev-Laikov grid. The 'nang' can be chosen from (6, 14, 26, 38, 50, 74, 86, + 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, + 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810). + k : int, optional + The order of the switching function in Becke's weighting scheme. + rotate : bool, optional + Whether to randomly rotate spherical grids. + + """ + if not isinstance(molecule, Molecule): + raise TypeError('Argument molecule should be an instance of Molecule class.') + coords, nums, pnums = molecule.coordinates, molecule.numbers, molecule.pseudo_numbers + return cls(coords, nums, pnums, specs, k, rotate) + + @classmethod + def from_file(cls, fname, specs='medium', k=3, rotate=False): + """Initialize the class given an instance of Molecule. + + Parameters + ---------- + fname : str + Path to molecule's files. + specs : str, optional + Specification of grid. Either choose from ['coarse', 'medium', 'fine', 'veryfine', + 'ultrafine', 'insane'] or provide a string of 'rname:rmin:rmax:nrad:nang' format. + Here 'rname' denotes the type of radial grid and can be chosen from ['linear', 'exp', + 'power'], 'rmin' and 'rmax' specify the first and last radial grid points in angstrom, + 'nrad' specify the number of radial grid points, and 'nang' specify the number of + angular Lebedev-Laikov grid. The 'nang' can be chosen from (6, 14, 26, 38, 50, 74, 86, + 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, + 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810). + k : int, optional + The order of the switching function in Becke's weighting scheme. + rotate : bool, optional + Whether to randomly rotate spherical grids. + + """ + mol = Molecule.from_file(fname) + return cls.from_molecule(mol, specs, k, rotate) + + def __getattr__(self, item): + return getattr(self._grid, item) + + @property + def center(self): + """Cartesian coordinates of atomic centers.""" + return self._coordinates + + @property + def coordinates(self): + """Cartesian coordinates of atomic centers.""" + return self._coordinates + + @property + def numbers(self): + """Atomic number of atomic centers.""" + return self._numbers + + @property + def pseudo_numbers(self): + """Pseudo atomic number of atomic centers.""" + return self._pseudo_numbers + + @property + def npoints(self): + """Number of grid points.""" + return self._grid.points.shape[0] + + @property + def points(self): + """Cartesian coordinates of grid points.""" + return self._grid.points + + @property + def weights(self): + """Integration weight of grid points.""" + return self._grid.weights + + def integrate(self, *value): + """Integrate the property evaluated on the grid points. + + Parameters + ---------- + value : np.ndarray + Property value evaluated on the grid points. + + """ + # temporary hack because of HORTON + if not isinstance(value, np.ndarray): + temp = value[0] + for item in value[1:]: + if item is not None: + temp = temp * item + value = temp + if value.ndim != 1: + raise ValueError('Argument value should be a 1D array.') + if value.shape != (self.npoints,): + raise ValueError('Argument value should have ({0},) shape!'.format(self.npoints)) + return self._grid.integrate(value) + + def compute_spherical_average(self, value): + """Compute spherical average of given value evaluated on the grid points. + + Note: This method only works for atomic systems with one nuclear center. + + Parameters + ---------- + value : np.ndarray + Property value evaluated on the grid points. + + """ + if len(self.numbers) != 1: + raise NotImplementedError('This method only works for systems with one atom!') + if value.ndim != 1: + raise ValueError('Argument value should be a 1D array.') + if value.shape != (self.npoints,): + raise ValueError('Argument value should have ({0},) shape!'.format(self.npoints)) + return self._grid.get_spherical_average(value) diff --git a/chemtools/wrappers3/molecule.py b/chemtools/wrappers3/molecule.py new file mode 100644 index 00000000..5c1765ee --- /dev/null +++ b/chemtools/wrappers3/molecule.py @@ -0,0 +1,655 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- +"""Wrapper Module.""" + + +import logging +import numpy as np +# from horton import IOData, DenseLinalgFactory + +try: + from importlib_resources import path +except ImportError: + from importlib.resources import path + + +__all__ = ["Molecule"] + + +class Molecule(object): + """Molecule class from HORTON package.""" + + def __init__(self, iodata): + """ + Initialize class. + + Parameters + ---------- + iodata : horton.IOData + An instance of horton.IOData object. + """ + self._iodata = iodata + if hasattr(self._iodata, "obasis"): + self._ao = AtomicOrbitals.from_molecule(self) + else: + self._ao = None + + self._coordinates = self._iodata.coordinates + self._numbers = self._iodata.numbers + + if hasattr(self._iodata, "exp_alpha"): + self._mo = MolecularOrbitals.from_molecule(self) + else: + self._mo = None + + @classmethod + def from_file(cls, fname): + """Initialize class given a file. + + Parameters + ---------- + fname : str + Path to molecule"s files. + + """ + # load molecule + logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s") + try: + iodata = IOData.from_file(str(fname)) + except IOError as _: + try: + with path("chemtools.data.examples", str(fname)) as fname: + logging.info("Loading {0}".format(str(fname))) + iodata = IOData.from_file(str(fname)) + except IOError as error: + logging.info(error) + return cls(iodata) + + def __getattr__(self, attr): + """Return attribute. + + Parameters + ---------- + attr : str + The name of attribute to retrieve. + + """ + value = getattr(self._iodata, attr, None) + return value + + @property + def coordinates(self): + """Cartesian coordinates of atomic centers.""" + return self._coordinates + + @property + def numbers(self): + """Atomic number of atomic centers.""" + return self._numbers + + @property + def ao(self): + """Atomic orbital instance.""" + return self._ao + + @property + def mo(self): + """Molecular orbital instance.""" + return self._mo + + def compute_density_matrix(self, spin="ab", index=None): + """Compute the density matrix array for the specified spin orbitals. + + Parameters + ---------- + spin : str, optional + The type of occupied spin orbitals. Options are "a" (for alpha), "b" (for beta), and + "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + return self.mo.compute_dm(spin, index=index)._array + + def compute_molecular_orbital(self, points, spin="ab", index=None): + """Return molecular orbitals. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + self._check_argument(points) + # assign orbital index (HORTON index the orbitals from 0) + if index is None: + # include all occupied orbitals of specified spin + spin_index = {"a": 0, "b": 1} + index = np.arange(self.mo.homo_index[spin_index[spin]]) + else: + # include specified set of orbitals + index = np.copy(np.asarray(index)) - 1 + if index.ndim == 0: + index = np.array([index]) + if np.any(index < 0): + raise ValueError("Argument index={0} cannot be less than one!".format(index + 1)) + + # get orbital expression of specified spin + if spin == 'b' and not hasattr(self._iodata, "exp_beta"): + exp = self._iodata.exp_alpha + else: + exp = getattr(self._iodata, "exp_" + {'a': 'alpha', 'b': 'beta'}[spin]) + return self._ao.compute_orbitals(exp, points, index) + + def compute_density(self, points, spin="ab", index=None): + r"""Return electron density. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + self._check_argument(points) + + # allocate output array + output = np.zeros((points.shape[0],), float) + + # compute density + if index is None: + # get density matrix corresponding to the specified spin + dm = self.mo.compute_dm(spin) + # include all orbitals + output = self._ao.compute_density(dm, points) + else: + # include subset of molecular orbitals + if spin == "ab": + # compute mo expression of alpha & beta orbitals + mo_a = self.compute_molecular_orbital(points, "a", index) + mo_b = self.compute_molecular_orbital(points, "b", index) + # add density of alpha & beta molecular orbitals + np.sum(mo_a**2, axis=1, out=output) + output += np.sum(mo_b**2, axis=1) + else: + # compute mo expression of specified molecular orbitals + mo = self.compute_molecular_orbital(points, spin, index) + # add density of specified molecular orbitals + np.sum(mo**2, axis=1, out=output) + return output + + def compute_gradient(self, points, spin="ab", index=None): + r"""Return gradient of the electron density. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + self._check_argument(points) + return self._ao.compute_gradient(self.mo.compute_dm(spin, index=index), points) + + def compute_hessian(self, points, spin="ab", index=None): + r"""Return hessian of the electron density. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + self._check_argument(points) + return self._ao.compute_hessian(self.mo.compute_dm(spin, index=index), points) + + def compute_laplacian(self, points, spin="ab", index=None): + r"""Return Laplacian of the electron density. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + hess = self.compute_hessian(points, spin, index) + return np.trace(hess, axis1=1, axis2=2) + + def compute_esp(self, points, spin="ab", index=None, charges=None): + r"""Return molecular electrostatic potential. + + The molecular electrostatic potential at point :math:`\mathbf{r}` is caused by the + electron density :math:`\rho` of the specified spin orbitals and set of point charges + :math:`\{q_A\}_{A=1}^{N_\text{atoms}}` placed at the position of the nuclei. i.e, + + .. math:: + V \left(\mathbf{r}\right) = + \sum_{A=1}^{N_\text{atoms}} \frac{q_A}{\rvert \mathbf{R}_A - \mathbf{r} \lvert} - + \int \frac{\rho \left(\mathbf{r}"\right)}{\rvert \mathbf{r}" - \mathbf{r} \lvert} + d\mathbf{r}" + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + charges : np.ndarray, optional + Array with shape (n,) representing the point charges at the position of the nuclei. + When ``None``, the pseudo numbers are used. + """ + self._check_argument(points) + # assign point charges + if charges is None: + charges = self.pseudo_numbers + elif not isinstance(charges, np.ndarray) or charges.shape != self.numbers.shape: + raise ValueError("Argument charges should be a 1d-array " + "with {0} shape.".format(self.numbers.shape)) + dm = self.mo.compute_dm(spin, index=index) + return self._ao.compute_esp(dm, points, self.coordinates, charges) + + def compute_ked(self, points, spin="ab", index=None): + r"""Return positive definite or Lagrangian kinetic energy density. + + .. math:: + \tau_\text{PD} \left(\mathbf{r}\right) = + \tfrac{1}{2} \sum_i^N n_i \rvert \nabla \phi_i \left(\mathbf{r}\right) \lvert^2 + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + self._check_argument(points) + return self._ao.compute_ked(self.mo.compute_dm(spin, index=index), points) + + def _check_argument(self, points): + """Check given arguments. + + Parameters + ---------- + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + if not isinstance(points, np.ndarray) or points.ndim != 2 or points.shape[1] != 3: + raise ValueError("Argument points should be a 2D-array with 3 columns.") + if not np.issubdtype(points.dtype, np.float64): + raise ValueError("Argument points should be a 2D-array of floats!") + if self._ao is None: + raise AttributeError("Atomic Orbitals information is needed!") + if self._mo is None: + raise AttributeError("Molecular Orbitals information is needed!") + + +class MolecularOrbitals(object): + """Molecular orbital class.""" + + def __init__(self, occs_a, occs_b, energy_a, energy_b, coeffs_a, coeffs_b): + self._occs_a, self._occs_b = occs_a, occs_b + self._energy_a, self._energy_b = energy_a, energy_b + self._coeffs_a, self._coeffs_b = coeffs_a, coeffs_b + + @classmethod + def from_molecule(cls, mol): + """Initialize class given an instance of `Molecule`. + + Parameters + ---------- + mol : Molecule + An instance of `Molecule` class. + + """ + if hasattr(mol, "exp_beta") and mol.exp_beta is not None: + exp_a, exp_b = mol.exp_alpha, mol.exp_beta + else: + exp_a, exp_b = mol.exp_alpha, mol.exp_alpha + occs_a, occs_b = exp_a.occupations, exp_b.occupations + energy_a, energy_b = exp_a.energies, exp_b.energies + coeffs_a, coeffs_b = exp_a.coeffs, exp_b.coeffs + return cls(occs_a, occs_b, energy_a, energy_b, coeffs_a, coeffs_b) + + @classmethod + def from_file(cls, fname): + """Initialize class given a file. + + Parameters + ---------- + fname : str + Path to molecule"s files. + + """ + return cls.from_molecule(Molecule.from_file(fname)) + + @property + def homo_index(self): + """Index of alpha and beta HOMO orbital.""" + # HORTON indexes the orbitals from 0, so 1 is added to get the intuitive index + index_a = np.argwhere(self._occs_a == 0.)[0, 0] + index_b = np.argwhere(self._occs_b == 0.)[0, 0] + return index_a, index_b + + @property + def lumo_index(self): + """Index of alpha and beta LUMO orbital.""" + # HORTON indexes the orbitals from 0, so 1 is added to get the intuitive index + return self.homo_index[0] + 1, self.homo_index[1] + 1 + + @property + def homo_energy(self): + """Energy of alpha and beta HOMO orbital.""" + return self._energy_a[self.homo_index[0] - 1], self._energy_b[self.homo_index[1] - 1] + + @property + def lumo_energy(self): + """Energy of alpha and beta LUMO orbital.""" + return self._energy_a[self.lumo_index[0] - 1], self._energy_b[self.lumo_index[1] - 1] + + @property + def occupation(self): + """Orbital occupation of alpha and beta electrons.""" + return self._occs_a, self._occs_b + + @property + def energy(self): + """Orbital energy of alpha and beta electrons.""" + return self._energy_a, self._energy_b + + @property + def coefficient(self): + """Orbital coefficient of alpha and beta electrons. + + The alpha and beta orbital coefficients are each storied in a 2d-array in which + the columns represent the basis coefficients of each molecular orbital. + """ + return self._coeffs_a, self._coeffs_b + + @property + def nelectrons(self): + """Number of alpha and beta electrons.""" + return np.sum(self._occs_a), np.sum(self._occs_b) + + def compute_overlap(self): + return NotImplementedError + + def compute_dm(self, spin="ab", index=None): + """Return HORTON density matrix object corresponding to the specified spin orbitals. + + Parameters + ---------- + spin : str, optional + Type of occupied spin orbitals which can be either "a" (for alpha), "b" (for + beta), and "ab" (for alpha + beta). + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + # temporary class because of HORTON2 + class DM(object): + def __init__(self, arr): + self._array = arr + + if spin == "ab": + return DM(self.compute_dm("a", index)._array + self.compute_dm("b", index)._array) + elif spin == "a": + arr = np.dot(self._coeffs_a * self._occs_a, self._coeffs_a.T) + elif spin == "b": + arr = np.dot(self._coeffs_b * self._occs_b, self._coeffs_b.T) + else: + raise ValueError("Argument spin={0} is not recognized!".format(spin)) + + if index is not None: + # convert to numpy array + index = np.asarray(index) + # check + if index.ndim == 0: + index = index.reshape(1) + if index.ndim >= 2: + raise ValueError("Indices should be given as a one-dimensional numpy array.") + index -= 1 + if np.any(index < 0): + raise ValueError( + "Indices cannot be less than 1. Note that indices start from 1." + ) + arr = arr[index[:, np.newaxis], index[np.newaxis, :]] + return DM(arr) + + +class AtomicOrbitals(object): + """Gaussian Basis Functions or Atomic Orbitals.""" + + def __init__(self, basis): + self._basis = basis + + @classmethod + def from_molecule(cls, mol): + """Initialize class given an instance of `Molecule`. + + Parameters + ---------- + mol : Molecule + An instance of `Molecule` class. + + """ + basis = mol._iodata.obasis + return cls(basis) + + @classmethod + def from_file(cls, fname): + """Initialize class given a file. + + Parameters + ---------- + fname : str + Path to molecule"s files. + + """ + return cls.from_molecule(Molecule.from_file(fname)) + + @property + def nbasis(self): + """int : number of basis functions.""" + return self._basis.nbasis + + @property + def center_index(self): + # FIXME: following code is pretty hacky. it will be used for the orbital partitioning code + # GOBasis object stores basis set to atom and angular momentum mapping + # by shell and not by contraction. So it needs to be converted + try: + ind_shell_atom = np.array(self._basis.shell_map) + ind_shell_orbtype = np.array(self._basis.shell_types) + + def num_contr(orbtype): + """Return the number of contractions for the given orbital type. + + Parameters + ---------- + orbtype : int + Horton's orbital type scheme. + Positive number corresponds to cartesian type and negative number corresponds to + spherical types. + + Returns + ------- + num_contr : int + Number of contractions for the given orbital type. + + """ + if orbtype < 0: + return 1 + 2 * abs(orbtype) + return 1 + orbtype + sum(i * (i + 1) / 2 for i in range(1, orbtype + 1)) + + numcontr_per_shell = np.array([num_contr(i) for i in ind_shell_orbtype]) + # Duplicate each shell information by the number of contractions in shell + ind_basis_center = np.repeat(ind_shell_atom, numcontr_per_shell) + # self._ind_basis_orbtype = np.repeat(ind_shell_orbtype, numcontr_per_shell) + except AttributeError: + pass + return ind_basis_center + + def compute_overlap(self): + r"""Return overlap matrix :math:`\mathbf{S}` of atomic orbitals. + + .. math:: [\mathbf{S}]_ij = int \phi_i(\mathbf{r}) \phi_j(\mathbf{r}) d\mathbf{r} + + """ + # make linear algebra factory + lf = DenseLinalgFactory(self.nbasis) + # compute overlap matrix + arr = self._basis.compute_overlap(lf)._array + return arr + + def compute_orbitals(self, dm, points, index): + """ + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + index : sequence of int, optional + Sequence of integers representing the occupied spin orbitals which are indexed + from 1 to :attr:`nbasis`. If ``None``, all orbitals of the given spin(s) are included. + + """ + return self._basis.compute_grid_orbitals_exp(dm, points, index) + + def compute_density(self, dm, points): + """Return electron density evaluated on the a set of points. + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + return self._basis.compute_grid_density_dm(dm, points) + + def compute_gradient(self, dm, points): + """Return gradient of the electron density evaluated on the a set of points. + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + return self._basis.compute_grid_gradient_dm(dm, points) + + def compute_hessian(self, dm, points): + """Return hessian of the electron density evaluated on the a set of points. + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + # compute upper triangular elements + output = self._basis.compute_grid_hessian_dm(dm, points) + # convert the (n, 6) shape to (n, 3, 3) + hess = np.zeros((len(points), 9)) + # NOTE: hard coded in the indices of the upper triangular matrix in the flattened form + # in C ordering. Maybe there is a numpy function that does this. This might fail if the + # hessian is not in c-ordering + hess[:, [0, 1, 2, 4, 5, 8]] = output + hess = hess.reshape(len(points), 3, 3) + hess += np.transpose(hess, axes=(0, 2, 1)) + for index in range(len(points)): + hess[index][np.diag_indices(3)] /= 2. + return hess + + def compute_esp(self, dm, points, coordinates, charges): + """Return electrostatic potential evaluated on the a set of points. + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + return self._basis.compute_grid_esp_dm(dm, coordinates, charges, points) + + def compute_ked(self, dm, points): + """Return positive definite kinetic energy density evaluated on the a set of points. + + Parameters + ---------- + dm : ndarray + First order reduced density matrix of B basis sets given as a 2D array of (B, B) shape. + points : ndarray + Cartesian coordinates of N points given as a 2D-array with (N, 3) shape. + + """ + return self._basis.compute_grid_kinetic_dm(dm, points) diff --git a/chemtools/wrappers3/part.py b/chemtools/wrappers3/part.py new file mode 100644 index 00000000..4aa97869 --- /dev/null +++ b/chemtools/wrappers3/part.py @@ -0,0 +1,176 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- +"""Wrapper of Part Module.""" + + +import numpy as np + +# from horton import ProAtomDB +# from horton.scripts.wpart import wpart_schemes + +from chemtools.wrappers3.molecule import Molecule +from chemtools.wrappers3.grid import MolecularGrid + + +__all__ = ['DensPart'] + + +class DensPart(object): + """Density-based Atoms-in-Molecules Partitioning Class.""" + + def __init__(self, coordinates, numbers, pseudo_numbers, density, grid, scheme="h", **kwargs): + """Initialize class. + + Parameters + ---------- + coordinates : np.ndarray, shape=(M, 3) + Cartesian coordinates of `M` atoms in the molecule. + numbers : np.ndarray, shape=(M,) + Atomic number of `M` atoms in the molecule. + pseudo_numbers : np.ndarray, shape=(M,) + Pseudo-number of `M` atoms in the molecule. + density : np.ndarray, shape=(N,) + Total density to be partitioned. + grid : BeckeMolGrid + Instance of BeckeMolGrid numerical integration grid. + scheme : str + Type of atoms-in-molecule partitioning scheme. + + """ + wpart = wpart_schemes[scheme] + # make proatom database + if scheme.lower() not in ["mbis", "b"]: + if "proatomdb" not in kwargs.keys() or kwargs['proatomdb'] is None: + if np.any(numbers > 18): + absent = list(numbers[numbers > 18]) + raise ValueError("Pro-atom for atomic number {} does not exist!".format(absent)) + proatomdb = ProAtomDB.from_refatoms(numbers) + kwargs["proatomdb"] = proatomdb + kwargs["local"] = False + # partition + self.part = wpart(coordinates, numbers, pseudo_numbers, grid, density, **kwargs) + self.part.do_charges() + + self.grid = grid + self.density = density + self.coordines = coordinates + self.numbers = numbers + self.pseudo_numbers = pseudo_numbers + self.charges = self.part['charges'] + + @classmethod + def from_molecule(cls, mol, scheme=None, grid=None, spin="ab", **kwargs): + """Initialize class given a Molecule instance. + + Parameters + ---------- + mol : Molecule + Instance of Molecule class. + scheme : str + Type of atoms-in-molecule partitioning scheme. + grid : MolecularGrid + Instance of MolecularGrid numerical integration grid. + spin : str + Type of occupied spin orbitals; choose either "a" (for alpha), "b" (for beta), + and "ab" (for alpha + beta). + + """ + if grid is None: + grid = MolecularGrid(mol.coordinates, mol.numbers, mol.pseudo_numbers, + specs="fine", rotate=False, k=3) + else: + check_molecule_grid(mol, grid) + # compute molecular electron density + dens = mol.compute_density(grid.points, spin=spin) + if mol.pesudo_numbers is None: + pesudo_numbers = mol.numbers + else: + pesudo_numbers = mol.pesudo_numbers + return cls(mol.coordinates, mol.numbers, pesudo_numbers, dens, grid, scheme, **kwargs) + + @classmethod + def from_file(cls, fname, scheme=None, grid=None, spin="ab", **kwargs): + """Initialize class given a wave-function file. + + Parameters + ---------- + fname : str + Path to molecule's files. + scheme : str + Type of atoms-in-molecule partitioning scheme. + grid : MolecularGrid + Instance of MolecularGrid integration grid. + spin : str + Type of occupied spin orbitals; choose either "a" (for alpha), "b" (for beta), + and "ab" (for alpha + beta). + + """ + mol = Molecule.from_file(fname) + return cls.from_molecule(mol, scheme=scheme, grid=grid, spin=spin, **kwargs) + + def condense_to_atoms(self, value, w_power=1): + condensed = np.zeros(self.part.natom) + for index in range(self.part.natom): + at_grid = self.part.get_grid(index) + at_weight = self.part.cache.load("at_weights", index) + # wcor = self.part.get_wcor(index) + local_prop = self.part.to_atomic_grid(index, value) + condensed[index] = at_grid.integrate(at_weight**w_power * local_prop) + return condensed + + def condense_to_fragments(self, value, fragments=None, w_power=1): + if fragments is None: + fragments = [[index] for index in range(self.part.natom)] + else: + # check fragments + segments = sorted([item for frag in fragments for item in frag]) + segments = np.array(segments) + if segments.size != self.numbers.size: + raise ValueError("Items in Fragments should uniquely represent all atoms.") + condensed = np.zeros(len(fragments)) + for index, frag in enumerate(fragments): + weight = np.zeros(self.grid.points.shape[0]) + for item in frag: + weight += self.part.cache.load("at_weights", item) + share = self.grid.integrate(weight ** w_power, value) + condensed[index] = share + return condensed + + +def check_molecule_grid(mol, grid): + """Check that molecule and grid represent the same system. + + Parameters + ---------- + mol : Molecule + Instance of Molecule class. + grid : MolecularGrid + Instance of MolecularGrid numerical integration grid. + + """ + if not np.max(abs(grid.centers - mol.coordinates)) < 1.e-6: + raise ValueError("Argument molecule & grid should have the same coordinates/centers.") + if not np.max(abs(grid.numbers - mol.numbers)) < 1.e-6: + raise ValueError("Arguments molecule & grid should have the same numbers.") + if not np.max(abs(grid.pseudo_numbers - mol.pseudo_numbers)) < 1.e-6: + raise ValueError("Arguments molecule & grid should have the same pseudo_numbers.") diff --git a/chemtools/wrappers3/test/__init__.py b/chemtools/wrappers3/test/__init__.py new file mode 100644 index 00000000..ec97ea1e --- /dev/null +++ b/chemtools/wrappers3/test/__init__.py @@ -0,0 +1,22 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 +# +# -- diff --git a/chemtools/wrappers/test/test_grid.py b/chemtools/wrappers3/test/test_grid.py similarity index 92% rename from chemtools/wrappers/test/test_grid.py rename to chemtools/wrappers3/test/test_grid.py index ce787c32..e169abcb 100644 --- a/chemtools/wrappers/test/test_grid.py +++ b/chemtools/wrappers3/test/test_grid.py @@ -20,20 +20,25 @@ # along with this program; if not, see # # -- -"""Test chemtools.wrappers.grid.""" +"""Test chemtools.wrappers2.grid.""" +import sys +import pytest import numpy as np + +pytestmark = pytest.mark.skipif(sys.version_info.major == 2, reason="Requires Python 3.") + +from numpy.testing import assert_raises, assert_allclose + +if sys.version_info.major != 2: + from chemtools.wrappers3.grid import MolecularGrid + try: from importlib_resources import path except ImportError: from importlib.resources import path -from numpy.testing import assert_raises, assert_allclose - -from chemtools.wrappers.grid import MolecularGrid -from chemtools.wrappers.molecule import Molecule - def test_wrapper_grid_raises(): with path('chemtools.data', 'ch4_uhf_ccpvdz.fchk') as fpath: diff --git a/chemtools/wrappers3/test/test_molecule.py b/chemtools/wrappers3/test/test_molecule.py new file mode 100644 index 00000000..b1ac6de3 --- /dev/null +++ b/chemtools/wrappers3/test/test_molecule.py @@ -0,0 +1,453 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools 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. +# +# ChemTools 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 chemtools.wrappers2.molecule.""" + + +import sys +import pytest +import numpy as np + +pytestmark = pytest.mark.skipif(sys.version_info.major == 2, reason="Requires Python 3.") + +from numpy.testing import assert_raises, assert_equal, assert_almost_equal + +if sys.version_info.major != 2: + from chemtools.wrappers3.molecule import Molecule + +try: + from importlib_resources import path +except ImportError: + from importlib.resources import path + + +def check_molecule_raises(mol): + """Check expected raised error messages by HortonWaveFunction class.""" + # example point array + points = np.array([[0., 0., 0.], [1., 1., 1.]]) + # check invalid orbital spin argument + assert_raises(ValueError, mol.compute_density_matrix, "alphabeta") + assert_raises(ValueError, mol.compute_density, points, spin="alpha") + assert_raises(ValueError, mol.compute_gradient, points, spin="beta") + assert_raises(ValueError, mol.compute_hessian, points, spin="betaalpha") + assert_raises(ValueError, mol.compute_ked, points, spin="balpha") + # check invalid points argument + assert_raises(ValueError, mol.compute_molecular_orbital, [0.1, 0.5, 0.7], "a") + assert_raises(ValueError, mol.compute_molecular_orbital, np.array([0.1, 0.5, 0.7]), "b") + assert_raises(ValueError, mol.compute_molecular_orbital, np.array([[0.4, 0.2]]), "a") + assert_raises(ValueError, mol.compute_molecular_orbital, np.array([[5, 10, 15]]), "b") + assert_raises(ValueError, mol.compute_density, [0., 0., 0.]) + assert_raises(ValueError, mol.compute_density, np.array([0., 0., 0.])) + assert_raises(ValueError, mol.compute_density, np.array([[0., 0.]])) + assert_raises(ValueError, mol.compute_density, np.array([[0, 0, 0]])) + assert_raises(ValueError, mol.compute_gradient, [.1, .2, .3]) + assert_raises(ValueError, mol.compute_gradient, np.array([.1, .2, .3])) + assert_raises(ValueError, mol.compute_gradient, np.array([[.1, 2., .3, .4]])) + assert_raises(ValueError, mol.compute_gradient, np.array([[1, 2, 3]])) + assert_raises(ValueError, mol.compute_hessian, [.5, .5, .5]) + assert_raises(ValueError, mol.compute_hessian, np.array([.5, .5, .5])) + assert_raises(ValueError, mol.compute_hessian, np.array([[.5, 5.]])) + assert_raises(ValueError, mol.compute_hessian, np.array([[5, 5, 5]])) + assert_raises(ValueError, mol.compute_esp, [1., .5]) + assert_raises(ValueError, mol.compute_esp, np.array([1., .5, .25])) + assert_raises(ValueError, mol.compute_esp, np.array([[1., .25]])) + assert_raises(ValueError, mol.compute_esp, np.array([[1, 25, 10]])) + assert_raises(ValueError, mol.compute_ked, [.5, 0., .2]) + assert_raises(ValueError, mol.compute_ked, np.array([.5, 0., .2])) + assert_raises(ValueError, mol.compute_ked, np.array([[.5, 0., .2, .1, .3]])) + assert_raises(ValueError, mol.compute_ked, np.array([[5, 0, 2]])) + # check invalid charges argument + assert_raises(ValueError, mol.compute_esp, points, charges=np.array([6., 1., 1.])) + assert_raises(ValueError, mol.compute_esp, points, charges=[6., 1., 1., 1., 1.]) + + +def test_molecule_raises_fchk_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_raises(molecule) + + +def test_molecule_raises_fchk_rhf_ch4(): + with path("chemtools.data", "ch4_rhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_raises(molecule) + + +def test_molecule_raises_wfn_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.wfn") as fname: + molecule = Molecule.from_file(fname) + check_molecule_raises(molecule) + + +def check_molecule_basics(mol): + """Check expected basic attributes of HortonWaveFunction class.""" + # check basic numerics + assert_equal(mol.natom, 5) + assert_equal(mol.mo.nelectrons, (5, 5)) + # assert_equal(mol.ao.nbasis, 34) + assert_equal(mol.numbers, [6, 1, 1, 1, 1]) + assert_equal(mol.pseudo_numbers, [6, 1, 1, 1, 1]) + assert_equal(mol.mo.homo_index, (5, 5)) + assert_equal(mol.mo.lumo_index, (6, 6)) + assert_equal(mol.mo.occupation[0], np.array([1] * 5 + [0] * 29)) + assert_equal(mol.mo.occupation[1], np.array([1] * 5 + [0] * 29)) + assert_almost_equal(mol.energy, -4.019868797400735E+01, decimal=8) + # check coordinates + coord = np.array([[-3.77945227E-05, 3.77945227E-05, -1.88972613E-05], + [ 1.04290206E+00, 1.50497789E+00, 9.34507367E-01], + [ 1.28607202E+00, -1.53098052E+00, -4.77307027E-01], + [-1.46467003E+00, -7.02997019E-01, 1.25954026E+00], + [-8.64474117E-01, 7.29131931E-01, -1.71670281E+00]]) + assert_almost_equal(mol.coordinates, coord, decimal=6) + + +def test_molecule_basics_fchk_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + # check basics + check_molecule_basics(molecule) + # check charges + esp = np.array([-0.502277518, 0.125567970, 0.125569655, 0.125566743, 0.125573150]) + assert_almost_equal(molecule.esp_charges, esp, decimal=6) + npa = np.array([-0.791299253, 0.197824989, 0.197825250, 0.197824326, 0.197824689]) + assert_almost_equal(molecule.npa_charges, npa, decimal=6) + mul = np.array([-0.139702704, 0.0349253868, 0.0349266071, 0.0349235395, 0.0349271707]) + assert_almost_equal(molecule.mulliken_charges, mul, decimal=6) + + +def test_molecule_basics_wfn_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.wfn") as fname: + molecule = Molecule.from_file(fname) + check_molecule_basics(molecule) + + +def test_molecule_orbitals_fchk_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + mol = Molecule.from_file(fname) + # check orbital energy + orb_energy = np.array([-1.12152085E+01, -9.42914385E-01, -5.43117091E-01, -5.43114279E-01, + -5.43101269E-01, 1.93295185E-01, 2.74358942E-01, 2.74359310E-01, + 2.74359740E-01, 5.89328697E-01, 5.89342443E-01, 5.89343893E-01, + 8.90214386E-01, 8.90219069E-01, 8.90222480E-01, 9.36275219E-01, + 1.13182813E+00, 1.13184117E+00, 1.25675685E+00, 1.68763897E+00, + 1.68764372E+00, 1.68765502E+00, 1.89570058E+00, 1.89570452E+00, + 1.89571385E+00, 2.21323213E+00, 2.21324619E+00, 2.21328532E+00, + 2.54691042E+00, 2.54694190E+00, 2.75532231E+00, 2.79906776E+00, + 2.79907762E+00, 2.79908651E+00]) + assert_almost_equal(mol.mo.energy[0], orb_energy, decimal=6) + assert_almost_equal(mol.mo.energy[1], orb_energy, decimal=6) + assert_almost_equal(mol.mo.homo_energy[0], orb_energy[4], decimal=6) + assert_almost_equal(mol.mo.homo_energy[1], orb_energy[4], decimal=6) + assert_almost_equal(mol.mo.lumo_energy[0], orb_energy[5], decimal=6) + assert_almost_equal(mol.mo.lumo_energy[1], orb_energy[5], decimal=6) + # check orbital coefficients + orb_coeff = np.array([9.97287609E-01, 1.86004593E-02, -8.24772487E-03]) + assert_almost_equal(mol.mo.coefficient[1][:3, 0], orb_coeff, decimal=6) + assert_almost_equal(mol.mo.coefficient[0][:3, 0], orb_coeff, decimal=6) + assert_almost_equal(mol.mo.coefficient[1][0, 1], -0.188285003, decimal=6) + assert_almost_equal(mol.mo.coefficient[0][0, 1], -0.188285003, decimal=6) + assert_almost_equal(mol.mo.coefficient[1][-1, -1], 1.02960200, decimal=6) + assert_almost_equal(mol.mo.coefficient[0][-1, -1], 1.02960200, decimal=6) + # check overlap matrix + overlap = mol.ao.compute_overlap() + assert_equal(overlap.shape, (34, 34)) + assert_almost_equal(np.diag(overlap), np.ones(34), decimal=6) + assert_almost_equal(overlap, overlap.T, decimal=6) + + +def test_molecule_density_matrix_fchk_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + mol = Molecule.from_file(fname) + # check density matrix against Gaussian (printed in log fname) + expected_diag = np.array([1.03003, 0.13222, 0.05565, 0.17944, 0.17944, 0.17944, 0.03891, + 0.03891, 0.03891, 0.00028, 0.00064, 0.00045, 0.00011, 0.00072, + 0.18575, 0.02710, 0.00044, 0.00072, 0.00038, 0.18575, 0.02710, + 0.00057, 0.00074, 0.00023, 0.18575, 0.02710, 0.00069, 0.00029, + 0.00056, 0.18575, 0.02710, 0.00035, 0.00030, 0.00089]) + # column 5, rows 15-21 + expected_05 = np.array([0.12066, 0.05027, -0.00560, -0.00252, -0.00502, -0.12275, -0.05114]) + # column 15, rows 16-34 + expected_15 = np.array([ 0.06838, -0.00691, -0.00996, -0.00619, -0.01612, -0.01573, 0.00244, + 0.00393, 0.00238, -0.01612, -0.01573, 0.00277, 0.00383, 0.00217, + -0.01612, -0.01573, 0.00270, 0.00366, 0.00253]) + # column 19, rows 20-34 + expected_19 = np.array([-0.00130, 0.00004, 0.00004, -0.00033, 0.00002, 0.00302, 0.00184, + 0.00001, -0.00010, -0.00003, -0.00438, -0.00125, -0.00032, 0.00003, + -0.00035]) + # column 29, rows 30-34 + expected_29 = np.array([-0.00442, -0.00106, -0.00003, 0.00029, -0.00047]) + # check alpha density matrix + dm_array_a = mol.compute_density_matrix(spin="a") + assert_almost_equal(np.diag(dm_array_a), expected_diag, decimal=5) + assert_almost_equal(dm_array_a, dm_array_a.T, decimal=5) + assert_almost_equal(dm_array_a[0, 1:3], np.array([-0.04982, -0.05262]), decimal=5) + assert_almost_equal(dm_array_a[4, 14:21], expected_05, decimal=5) + assert_almost_equal(dm_array_a[14, 15:], expected_15, decimal=5) + assert_almost_equal(dm_array_a[18, 19:], expected_19, decimal=5) + assert_almost_equal(dm_array_a[28, 29:], expected_29, decimal=5) + # check beta density matrix + dm_array_b = mol.compute_density_matrix(spin="b") + assert_almost_equal(np.diag(dm_array_b), expected_diag, decimal=5) + assert_almost_equal(dm_array_b, dm_array_b.T, decimal=5) + assert_almost_equal(dm_array_b[0, 1:3], np.array([-0.04982, -0.05262]), decimal=5) + assert_almost_equal(dm_array_b[4, 14:21], expected_05, decimal=5) + assert_almost_equal(dm_array_b[14, 15:], expected_15, decimal=5) + assert_almost_equal(dm_array_b[18, 19:], expected_19, decimal=5) + assert_almost_equal(dm_array_b[28, 29:], expected_29, decimal=5) + # check total density matrix + dm_array_ab = mol.compute_density_matrix(spin="ab") + assert_almost_equal(np.diag(dm_array_ab), 2 * expected_diag, decimal=5) + assert_almost_equal(dm_array_ab, dm_array_ab.T, decimal=5) + assert_almost_equal(dm_array_ab[0, 1:3], 2*np.array([-0.04982, -0.05262]), decimal=5) + assert_almost_equal(dm_array_ab[4, 14:21], 2*expected_05, decimal=5) + assert_almost_equal(dm_array_ab[14, 15:], 2*expected_15, decimal=5) + assert_almost_equal(dm_array_ab[18, 19:], 2*expected_19, decimal=5) + assert_almost_equal(dm_array_ab[28, 29:], 2*expected_29, decimal=5) + + +def test_molecule_esp_fchk_uhf_ch4(): + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + mol = Molecule.from_file(fname) + # check esp against Gaussian (printed in log file) + # check esp at the position of each nuclei (1.e-14 is added to avoid division by zero) + # excluding the nucleus itself. + point = mol.coordinates[0].reshape(1, 3) + 1.e-14 + charge = np.array([0., 1., 1., 1., 1.]) + assert_almost_equal(mol.compute_esp(point, charges=charge), [-14.745629], decimal=5) + point = mol.coordinates[1].reshape(1, 3) + 1.e-14 + charge = np.array([6., 0., 1., 1., 1.]) + assert_almost_equal(mol.compute_esp(point, charges=charge), [-1.116065], decimal=5) + point = mol.coordinates[2].reshape(1, 3) + 1.e-14 + charge = np.array([6., 1., 0., 1., 1.]) + assert_almost_equal(mol.compute_esp(point, charges=charge), [-1.116065], decimal=5) + point = mol.coordinates[3].reshape(1, 3) + 1.e-14 + charge = np.array([6., 1., 1., 0., 1.]) + assert_almost_equal(mol.compute_esp(point, charges=charge), [-1.116067], decimal=5) + point = mol.coordinates[4].reshape(1, 3) + 1.e-14 + charge = np.array([6., 1., 1., 1., 0.]) + assert_almost_equal(mol.compute_esp(point, charges=charge), [-1.116065], decimal=5) + # check esp at non-nuclei points + points = np.array([[ 0.5, 0.5, 0.5], + [-0.5, -0.5, -0.5], + [-0.5, 0.5, 0.5], + [-0.5, -0.5, 0.5], + [-0.5, 0.5, -0.5], + [ 0.5, -0.5, -0.5], + [ 0.5, -0.5, 0.5], + [ 0.5, 0.5, -0.5]]) / 0.529177 + expected_esp = np.array([0.895650, 0.237257, 0.234243, 0.708301, + 0.499083, 0.479275, 0.241434, 0.235102]) + assert_almost_equal(mol.compute_esp(points), expected_esp, decimal=5) + + +def check_molecule_against_gaussian_ch4(mol): + """Check local properties of HortonWaveFunction class against Gaussian09_C.01"s cubegen.""" + # get expected data computed by Gaussian09_C.01 cubegen + with path("chemtools.data", "data_cubegen_g09_C01_ch4_uhf_ccpvdz.npz") as fname: + data = np.load(str(fname)) + points, dens, grad = data["points"], data["dens"], data["grad"] + lap, hess_xx, esp = data["lap"], data["hess_xx"], data["esp"] + # check density, gradient, esp & hessian + assert_almost_equal(mol.compute_density(points, "ab"), dens, decimal=5) + assert_almost_equal(mol.compute_density(points, "a"), 0.5 * dens, decimal=5) + assert_almost_equal(mol.compute_density(points, "b"), 0.5 * dens, decimal=5) + assert_almost_equal(mol.compute_gradient(points, "ab"), grad, decimal=5) + assert_almost_equal(mol.compute_esp(points, "ab"), esp, decimal=5) + assert_almost_equal(mol.compute_laplacian(points, "ab", None), lap, decimal=5) + assert_almost_equal(mol.compute_hessian(points, "ab", None)[:, 0, 0], hess_xx, decimal=5) + # density computed by summing squared mo expressions + assert_almost_equal(mol.compute_density(points, "ab", range(1, 6)), dens, decimal=5) + assert_almost_equal(mol.compute_density(points, "a", range(1, 6)), 0.5 * dens, decimal=5) + + +def test_molecule_grid_g09_fchk_uhf_ch4(): + # make an instance of molecule from fchk file + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_against_gaussian_ch4(molecule) + + +def test_molecule_grid_g09_fchk_rhf_ch4(): + # make an instance of molecule from fchk file + with path("chemtools.data", "ch4_rhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_against_gaussian_ch4(molecule) + + +def test_molecule_grid_g09_wfn_uhf_ch4(): + # make an instance of molecule from wfn file + with path("chemtools.data", "ch4_uhf_ccpvdz.wfn") as fname: + molecule = Molecule.from_file(fname) + check_molecule_against_gaussian_ch4(molecule) + + +def check_molecule_against_fortran_ch4(mol): + # get expected data computed by Fortran code + with path("chemtools.data", "data_fortran_ch4_uhf_ccpvdz.npz") as fname: + data = np.load(str(fname)) + points, exp8, exp9 = data["points"], data["orb_08"], data["orb_09"] + dens, grad, ke = data["dens"], data["grad"], data["ked_pd"] + # check density & gradient + assert_almost_equal(mol.compute_density(points, "ab", None), dens, decimal=6) + assert_almost_equal(mol.compute_gradient(points, "ab", None), grad, decimal=6) + assert_almost_equal(mol.compute_density(points, "a", None), 0.5 * dens, decimal=6) + assert_almost_equal(mol.compute_density(points, "b", None), 0.5 * dens, decimal=6) + # check density computed by summing squared mo expressions + assert_almost_equal(mol.compute_density(points, "ab", range(1, 6)), dens, decimal=6) + assert_almost_equal(mol.compute_density(points, "a", range(1, 6)), 0.5 * dens, decimal=6) + assert_almost_equal(mol.compute_density(points, "b", range(1, 6)), 0.5 * dens, decimal=6) + # check mo expression + assert_almost_equal(mol.compute_molecular_orbital(points, "a", 8)[:, 0], exp8, decimal=6) + assert_almost_equal(mol.compute_molecular_orbital(points, "b", 8)[:, 0], exp8, decimal=6) + assert_almost_equal(mol.compute_molecular_orbital(points, "a", 9)[:, 0], exp9, decimal=6) + assert_almost_equal(mol.compute_molecular_orbital(points, "b", 9)[:, 0], exp9, decimal=6) + # check positive definite ke + assert_almost_equal(mol.compute_ked(points, "ab"), ke, decimal=6) + + +def test_molecule_grid_fortran_fchk_uhf_ch4(): + # make an instance of molecule + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_against_fortran_ch4(molecule) + + +def test_molecule_grid_fortran_fchk_rhf_ch4(): + # make an instance of molecule + with path("chemtools.data", "ch4_rhf_ccpvdz.fchk") as fname: + molecule = Molecule.from_file(fname) + check_molecule_against_fortran_ch4(molecule) + + +def test_molecule_basic_fchk_uhf_o2(): + with path("chemtools.data", "o2_uhf_virtual.fchk") as fname: + mol = Molecule.from_file(fname) + # print mol.nelectrons + # check basic numerics + assert_equal(mol.natom, 2) + assert_equal(mol.mo.nelectrons, (9, 7)) + assert_equal(mol.ao.nbasis, 44) + assert_equal(mol.numbers, [8, 8]) + assert_equal(mol.pseudo_numbers, [8, 8]) + assert_equal(mol.mo.homo_index, (9, 7)) + assert_equal(mol.mo.lumo_index, (10, 8)) + assert_equal(mol.mo.occupation[0], np.array([1] * 9 + [0] * 35)) + assert_equal(mol.mo.occupation[1], np.array([1] * 7 + [0] * 37)) + assert_almost_equal(mol.energy, -1.496641407696776E+02, decimal=8) + # check coordinates + coord = np.array([[0.0, 0.0, 1.09452498], [0.0, 0.0, -1.09452498]]) + assert_almost_equal(mol.coordinates, coord, decimal=6) + # energy of alpha orbitals + orb_energy_a = np.array([-2.07520003E+01, -2.07511624E+01, -1.77073830E+00, -1.19176563E+00, + -8.67505123E-01, -8.67505123E-01, -7.86590608E-01, -5.41367609E-01, + -5.41367609E-01, 1.75551613E-01, 1.79578059E-01, 2.12136240E-01, + 2.12136240E-01, 2.82746427E-01, 2.82746427E-01, 2.85693824E-01, + 4.25803814E-01, 5.48137441E-01, 1.13330121E+00, 1.13563801E+00, + 1.13563801E+00, 1.21130153E+00, 1.21130153E+00, 1.21869395E+00, + 1.42273629E+00, 1.74301463E+00, 2.54671907E+00, 2.54671907E+00, + 2.83314193E+00, 2.83314193E+00, 3.16996201E+00, 3.16996201E+00, + 3.35513494E+00, 3.91418793E+00, 3.91418793E+00, 4.32509152E+00, + 5.22427895E+00, 5.22427895E+00, 5.43561012E+00, 5.51122695E+00, + 5.51122695E+00, 6.78081131E+00, 5.12932061E+01, 5.15031931E+01]) + # energy of beta orbitals + orb_energy_b = np.array([-2.06970266E+01, -2.06955835E+01, -1.64286825E+00, -9.81871414E-01, + -7.18265821E-01, -6.01903968E-01, -6.01903968E-01, 1.04744787E-01, + 1.04744787E-01, 1.82240025E-01, 1.84775146E-01, 2.25175431E-01, + 2.25175431E-01, 2.81988319E-01, 3.22590360E-01, 3.22590360E-01, + 4.31522630E-01, 6.25012892E-01, 1.14414927E+00, 1.21805822E+00, + 1.21805822E+00, 1.24392742E+00, 1.30587863E+00, 1.30587863E+00, + 1.45650216E+00, 1.79253113E+00, 2.62926234E+00, 2.62926234E+00, + 2.95890359E+00, 2.95890359E+00, 3.32630000E+00, 3.32630000E+00, + 3.42846106E+00, 3.99627997E+00, 3.99627997E+00, 4.36808390E+00, + 5.33007026E+00, 5.33007026E+00, 5.45827702E+00, 5.61601037E+00, + 5.61601037E+00, 6.81582045E+00, 5.13257489E+01, 5.15352582E+01]) + # check orbital energy + assert_almost_equal(mol.mo.energy[0], orb_energy_a, decimal=6) + assert_almost_equal(mol.mo.energy[1], orb_energy_b, decimal=6) + assert_almost_equal(mol.mo.homo_energy[0], orb_energy_a[8], decimal=6) + assert_almost_equal(mol.mo.homo_energy[1], orb_energy_b[6], decimal=6) + assert_almost_equal(mol.mo.lumo_energy[0], orb_energy_a[9], decimal=6) + assert_almost_equal(mol.mo.lumo_energy[1], orb_energy_b[7], decimal=6) + assert_almost_equal(mol.mulliken_charges, 0.0, decimal=6) + # check orbital coefficients + assert_almost_equal(mol.mo.coefficient[0][:3, 0], + np.array([0.389497609, 0.333421243, 0.]), decimal=6) + + +def test_molecule_density_matrix_index_fchk_uhf_ch4(): + # check get_density_matrix for different values of the index/ + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + mol = Molecule.from_file(fname) + dm_full = mol.mo.compute_dm("a")._array + # errors + assert_raises(ValueError, mol.compute_dm, "a", [[1]]) + assert_raises(ValueError, mol.compute_dm, "a", [0]) + # one index + for i in range(1, mol.ao.nbasis + 1): + assert np.allclose(dm_full[i - 1, i - 1], mol.mo.compute_dm("a", i)._array) + # multiple indices + for i in range(1, mol.ao.nbasis + 1): + for j in range(1, mol.ao.nbasis + 1): + # NOTE: indices can be repeated + indices = np.array([i -1, j - 1]) + assert np.allclose(dm_full[indices[:, None], indices[None, :]], + mol.mo.compute_dm("a", [i, j])._array) + + +def test_molecule_horton_h2o(): + with path("chemtools.data", "data_horton_fchk_h2o_ub3lyp_ccpvtz.npz") as fname: + data = np.load(str(fname)) + with path("chemtools.data", "h2o_q+0_ub3lyp_ccpvtz.fchk") as fname: + mol = Molecule.from_file(fname) + # check properties computed at nucleus against HORTON + points = data["coords"] + assert np.allclose(mol.compute_density(points), data["nuc_dens"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_gradient(points), data["nuc_grad"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_hessian(points), data["nuc_hess"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_ked(points), data["nuc_ked_pd"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_esp(points), data["nuc_esp"], rtol=0., atol=1.e-6) + # check properties computed on a grid against HORTON + assert np.allclose(mol.compute_density(data["points"]), data["dens"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_gradient(data["points"]), data["grad"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_hessian(data["points"]), data["hess"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_ked(data["points"]), data["ked_pd"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_esp(data["points"]), data["esp"], rtol=0., atol=1.e-6) + + +def test_molecule_horton_ch4(): + with path("chemtools.data", "data_horton_fchk_ch4_uhf_ccpvdz.npz") as fname: + data = np.load(str(fname)) + with path("chemtools.data", "ch4_uhf_ccpvdz.fchk") as fname: + mol = Molecule.from_file(fname) + # check properties computed at nucleus against HORTON + points = data["coords"] + assert np.allclose(mol.compute_density(points), data["nuc_dens"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_gradient(points), data["nuc_grad"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_hessian(points), data["nuc_hess"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_ked(points), data["nuc_ked_pd"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_esp(points), data["nuc_esp"], rtol=0., atol=1.e-6) + # check properties computed on a grid against HORTON + assert np.allclose(mol.compute_density(data["points"]), data["dens"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_gradient(data["points"]), data["grad"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_hessian(data["points"]), data["hess"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_ked(data["points"]), data["ked_pd"], rtol=0., atol=1.e-6) + assert np.allclose(mol.compute_esp(data["points"]), data["esp"], rtol=0., atol=1.e-6) diff --git a/chemtools/wrappers/test/test_part.py b/chemtools/wrappers3/test/test_part.py similarity index 89% rename from chemtools/wrappers/test/test_part.py rename to chemtools/wrappers3/test/test_part.py index c58b40a4..5ec0d635 100644 --- a/chemtools/wrappers/test/test_part.py +++ b/chemtools/wrappers3/test/test_part.py @@ -21,11 +21,18 @@ # # -- + +import sys +import pytest import numpy as np -from chemtools.wrappers.molecule import Molecule -from chemtools.wrappers.part import DensPart -from chemtools.wrappers.grid import MolecularGrid +pytestmark = pytest.mark.skip(sys.version_info.major == 2, reason="Requires Python 3.") + +if sys.version_info.major != 2: + from chemtools.wrappers3.molecule import Molecule + from chemtools.wrappers3.part import DensPart + from chemtools.wrappers3.grid import MolecularGrid + try: from importlib_resources import path except ImportError: diff --git a/setup.py b/setup.py index 7bf89b2c..c393ba79 100755 --- a/setup.py +++ b/setup.py @@ -46,7 +46,7 @@ install_requires=[ 'numpy>=1.16', 'matplotlib', 'Pillow', 'Image', 'sympy', 'scipy', 'importlib_resources; python_version < "3.7"', - 'nose', + 'nose', 'pytest', ], extras_require={ 'doc': [ diff --git a/tools/conda.recipe/meta.yaml b/tools/conda.recipe/meta.yaml index 0a0e8074..77ee18de 100644 --- a/tools/conda.recipe/meta.yaml +++ b/tools/conda.recipe/meta.yaml @@ -18,6 +18,7 @@ requirements: - horton - matplotlib - nose + - pytest - importlib_resources build: - python @@ -26,6 +27,7 @@ requirements: - sympy - matplotlib - nose + - pytest - horton - git-lfs - sphinxcontrib-bibtex