diff --git a/NOTICES b/NOTICES index 94ff94a..4c9ee9f 100644 --- a/NOTICES +++ b/NOTICES @@ -1,37 +1,57 @@ -The Molecular Design Toolkit incorporates code from the following sources: +This project incorporates code from the following sources: - Transformations.py ------------------------------- -Copyright (c) 2006-2015, Christoph Gohlke -Copyright (c) 2006-2015, The Regents of the University of California -SOURCE CODE: moldesign/external/transformations.py -DESCRIPTION: Included in its entirety, with minor modifications to imports -LICENSE: moldesign/external/transformations.py -WEBSITE: http://www.lfd.uci.edu/~gohlke/ + Python +-------------------------------- +Copyright (c) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, +2011, 2012, 2013, 2014, 2015, 2016 Python Software Foundation; All Rights +Reserved +SOURCE CODE: moldesign/utils/utils.py (`methodcaller` class) +DESCRIPTION: unmodified copy of the pure-python `methodcaller` class +WEBSITE: https://www.python.org +LICENSE: moldesign/external/licenses/python + + + PyQuante2 +-------------------------------- +By Rick Muller +Copyright (c) 2013 PyQuante Development Team +SOURCE CODE: moldesign/external/pyquante +WEBSITE: http://pyquante.sourceforge.net/ +DESCRIPTION: Code for evaulating one-electron integrals of spherical gaussian basis +functions. Except for the imports, the files are unmodified. +LICENSE: moldesign/external/pyquante/LICENSE - Sphinx + Sphinx ------------------------- Copyright (c) 2007-2016 by the Sphinx team SOURCE CODE: moldesign/utils/docparsers/google.py DESCRIPTION: Includes code derived from the Napoleon Google-style docstring parsers. Note that these routines have been heavily modified as described in the file. -LICENSE: moldesign/utils/docparsers/sphinxlicense WEBSITE: http://sphinx-doc.org +LICENSE: moldesign/utils/docparsers/sphinxlicense - Python --------------------------------- -Copyright (c) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, -2011, 2012, 2013, 2014, 2015, 2016 Python Software Foundation; All Rights -Reserved -SOURCE CODE: moldesign/utils/utils.py (`methodcaller` class) -DESCRIPTION: unmodified copy of the pure-python `methodcaller` class -LICENSE: moldesign/external/licenses/python -WEBSITE: https://www.python.org + Symmol.f +------------------------------------ +By Tullio Pilati and Alessandra Forni +Released by the CCP14 Project +SOURCE CODE: DockerMakefiles/buildfiles/symmol/symmol.f +WEBSITE: http://www.ccp14.ac.uk +LICENSE: None available + + + Transformations.py +------------------------------ +Copyright (c) 2006-2015, Christoph Gohlke +Copyright (c) 2006-2015, The Regents of the University of California +SOURCE CODE: moldesign/external/transformations.py +DESCRIPTION: Included in its entirety, with minor modifications to imports +LICENSE: moldesign/external/transformations.py +WEBSITE: http://www.lfd.uci.edu/~gohlke/ - Versioneer + Versioneer ------------------------- By Brian Warner SOURCE CODE: versioneer.py, moldesign/_version.py @@ -40,11 +60,3 @@ LICENSE: To make Versioneer easier to embed, all its code is dedicated to the pu The _version.py that it creates is also in the public domain. Specifically, both are released under the Creative Commons "Public Domain Dedication" license (CC0-1.0), as described in https://creativecommons.org/publicdomain/zero/1.0/ - - Symmol.f ------------------------------------- -By Tullio Pilati and Alessandra Forni -Released by the CCP14 Project -SOURCE CODE: docker_images/symmol/symmol.f -WEBSITE: http://www.ccp14.ac.uk -LICENSE: None available diff --git a/moldesign/_static_data/pint_atomic_units.txt b/moldesign/_static_data/pint_atomic_units.txt index b08f425..52c83a1 100644 --- a/moldesign/_static_data/pint_atomic_units.txt +++ b/moldesign/_static_data/pint_atomic_units.txt @@ -3,11 +3,13 @@ atomic_time = hbar / hartree = t0 kcalpermol = kcal / 6.0221412927e23 kjpermol = kJ / 6.0221412927e23 gpermol = g / 6.0221412927e23 -fs = femtosecond = 1.0 * femtosecond -ps = picosecond = 1.0 * piscosecond -amu = 1.0 * atomic_mass_units -ang = 1.0 * angstrom -eV = electron_volt debye = 0.393430307 * elementary_charge * bohr electron_charge = elementary_charge -q_e = elementary_charge \ No newline at end of file +q_e = elementary_charge + +@system nano using international + nanometers + nanoseconds + atomic_mass_units +@end + diff --git a/moldesign/_tests/conftest.py b/moldesign/_tests/conftest.py index 108f711..70cf2d6 100644 --- a/moldesign/_tests/conftest.py +++ b/moldesign/_tests/conftest.py @@ -1,8 +1,13 @@ def pytest_itemcollected(item): - if hasattr(item.module, '__PYTEST_MARK__'): - item.add_marker(item.module.__PYTEST_MARK__) + marks = getattr(item.module, '__PYTEST_MARK__', None) + if marks is None: + return + if isinstance(marks, str): + marks = [marks] + for mark in marks: + item.add_marker(mark) # TODO: nicer output strings for git commit status diff --git a/moldesign/_tests/helpers.py b/moldesign/_tests/helpers.py index e89adc5..07fc108 100644 --- a/moldesign/_tests/helpers.py +++ b/moldesign/_tests/helpers.py @@ -229,3 +229,20 @@ def assert_almost_equal(actual, desired, **kwargs): reason='Network connection timed out') """ Decorator to disable tests that need internet if we can't connect to the network """ + +def generate_grid(g, g2=None): + from moldesign import mathutils + + if g2 is None: g2 = g + if not hasattr(g, 'alpha'): + alpha = min(min(p.alpha for p in g), min(p.alpha for p in g2)) + else: + alpha = min(g.alpha, g2.alpha) + + width = np.sqrt(1.0/alpha) + ranges = np.ones((3, 2))*5.0*width + ranges[:, 0] *= -1 + ranges += ((g.center + g2.center)/2.0)[:, None] + grid = mathutils.VolumetricGrid(*ranges, npoints=64) + allpoints = grid.allpoints() + return allpoints, grid diff --git a/moldesign/_tests/molecule_fixtures.py b/moldesign/_tests/molecule_fixtures.py index c94bdbd..d2ff3b5 100644 --- a/moldesign/_tests/molecule_fixtures.py +++ b/moldesign/_tests/molecule_fixtures.py @@ -4,10 +4,11 @@ import moldesign as mdt import moldesign.units as u +from moldesign.utils import exports from .helpers import get_data_path -from moldesign.interfaces.openbabel import force_remote as openbabel_missing +__all__ = ['molecule_standards'] molecule_standards = {} @@ -16,6 +17,7 @@ def typedfixture(*types, **kwargs): We'll later use this type to determine what tests to run on the result""" def fixture_wrapper(func): + __all__.append(func.__name__) for t in types: molecule_standards.setdefault(t, []).append(func.__name__) return pytest.fixture(**kwargs)(func) @@ -56,12 +58,14 @@ def ethylene_waterbox_2na_2cl(): return solvated +@exports @pytest.fixture def random_atoms_from_3aid(pdb3aid): atoms = mdt.molecules.atomcollections.AtomList(random.sample(pdb3aid.atoms, 10)) return atoms +@exports @pytest.fixture(scope='session') def cached_small_molecule(): mol = mdt.from_smiles('CNCOS(=O)C') @@ -69,16 +73,19 @@ def cached_small_molecule(): return mol +@exports @pytest.fixture def small_molecule(cached_small_molecule): return cached_small_molecule.copy() +@exports @pytest.fixture(scope='session') def cached_benzene(): return mdt.from_smiles('c1ccccc1') +@exports @pytest.fixture def benzene(cached_benzene): return cached_benzene.copy() @@ -104,25 +111,30 @@ def heh_plus(): return mol +@exports @pytest.fixture(scope='session') def cached_ethylene(): return mdt.from_smiles('C=C') +@exports @pytest.fixture def ethylene(cached_ethylene): return cached_ethylene.copy() +@exports @pytest.fixture(scope='session') def cached_pdb1yu8(): return mdt.read(get_data_path('1yu8.pdb')) +@exports @pytest.fixture def pdb1yu8(): return mdt.read(get_data_path('1yu8.pdb')) +@exports @pytest.fixture(scope='session') def cached_mol_from_xyz(): return mdt.read("""43 @@ -173,11 +185,13 @@ def cached_mol_from_xyz(): """, format='xyz') +@exports @pytest.fixture def mol_from_xyz(cached_mol_from_xyz): return cached_mol_from_xyz.copy() +@exports @pytest.fixture(scope='session') def cached_mol_from_sdf(): return mdt.read(""" @@ -209,6 +223,7 @@ def cached_mol_from_sdf(): $$$$ """, format='sdf') +@exports @pytest.fixture def mol_from_sdf(cached_mol_from_sdf): return cached_mol_from_sdf.copy() @@ -224,6 +239,7 @@ def nucleic(): ######################################################################################## # Molecules with forcefields assigned - these use a session-scoped constructor w/ a copy factory +@exports @pytest.fixture(scope='session') def cached_mol_parameterized_with_zeros(cached_small_molecule): return _param_small_mol(cached_small_molecule.copy(), 'zero') @@ -234,6 +250,7 @@ def mol_with_zerocharge_params(cached_mol_parameterized_with_zeros): return cached_mol_parameterized_with_zeros.copy() +@exports @pytest.fixture(scope='session') def cached_mol_parameterized_with_am1bcc(cached_small_molecule): """ We don't use this fixture directly, rather use another fixture that copies these results @@ -247,6 +264,7 @@ def mol_with_am1bcc_params(cached_mol_parameterized_with_am1bcc): return cached_mol_parameterized_with_am1bcc.copy() +@exports @pytest.fixture(scope='session') def cached_mol_parameterized_with_gasteiger(cached_small_molecule): """ We don't use this fixture directly, rather use another fixture that copies these results @@ -271,6 +289,7 @@ def _param_small_mol(cached_small_molecule, chargemodel): return mol +@exports @pytest.fixture(scope='session') def cached_protein_with_default_amber_ff(cached_pdb1yu8): """ We don't use this fixture directly, rather use another fixture that copies these results @@ -288,14 +307,46 @@ def protein_default_amber_forcefield(cached_protein_with_default_amber_ff): return cached_protein_with_default_amber_ff.copy() +@exports @pytest.fixture(scope='session') -def cached_h2_rhfwfn(): +def cached_h2_rhf_sto3g(): mol = h2() # fixture is not cached, so just call it directly mol.set_energy_model(mdt.models.PySCFPotential, basis='sto-3g', theory='rhf') mol.calculate(requests=['forces']) return mol +@exports @pytest.fixture -def h2_rhfwfn(cached_h2_rhfwfn): - return cached_h2_rhfwfn.copy() +def h2_rhf_sto3g(cached_h2_rhf_sto3g): + return cached_h2_rhf_sto3g.copy() + + +@exports +@pytest.fixture(scope='session') +def cached_h2_rhf_augccpvdz(): + mol = h2() + mol.set_energy_model(mdt.models.RHF, basis='aug-cc-pvdz') + mol.calculate() + return mol + + +@exports +@pytest.fixture +def h2_rhf_augccpvdz(cached_h2_rhf_augccpvdz): + return cached_h2_rhf_augccpvdz.copy() + + +@exports +@pytest.fixture(scope='session') +def cached_acetylene_dft_631g(): + mol = mdt.from_smiles('C#C') + mol.set_energy_model(mdt.models.B3LYP, basis='6-31g') + mol.calculate() + return mol + + +@exports +@pytest.fixture +def acetylene_dft_631g(cached_acetylene_dft_631g): + return cached_acetylene_dft_631g.copy() diff --git a/moldesign/_tests/object_fixtures.py b/moldesign/_tests/object_fixtures.py index 54f991f..aa4098b 100644 --- a/moldesign/_tests/object_fixtures.py +++ b/moldesign/_tests/object_fixtures.py @@ -18,11 +18,15 @@ registered_types = {key:val[:] for key,val in molecule_standards.items()} +__all__ = ['registered_types'] + + def typedfixture(*types, **kwargs): """This is a decorator that lets us associate fixtures with one or more arbitrary types. We'll later use this type to determine what tests to run on the result""" def fixture_wrapper(func): + __all__.append(func.__name__) for t in types: registered_types.setdefault(t, []).append(func.__name__) return pytest.fixture(**kwargs)(func) @@ -58,7 +62,7 @@ def simple_unit_array(): @typedfixture('pickleable', 'equality') def unit_number(): - return 391.23948 * u.ureg.kg * u.ang / u.alpha + return 391.23948 * u.ureg.kg * u.angstrom / u.alpha ###################################### @@ -104,13 +108,13 @@ def mol_bond_graph(h2): @typedfixture('pickleable') -def mol_wfn(h2_rhfwfn): - return h2_rhfwfn.copy().wfn +def mol_wfn(h2_rhf_sto3g): + return h2_rhf_sto3g.copy().wfn @typedfixture('pickleable') -def mol_properties(h2_rhfwfn): - return h2_rhfwfn.copy().properties +def mol_properties(h2_rhf_sto3g): + return h2_rhf_sto3g.copy().properties @typedfixture('trajectory') @@ -153,3 +157,5 @@ def h2_harmonic_atoms(h2_harmonic): registered_types['trajectory'] + registered_types['atom']) pickleable = registered_types['pickleable'] + moldesign_objects + +__all__.extend(['moldesign_objects', 'pickleable']) diff --git a/moldesign/_tests/test_atom_bond_computed_properties.py b/moldesign/_tests/test_atom_bond_computed_properties.py index 0bde8b3..d261cdb 100644 --- a/moldesign/_tests/test_atom_bond_computed_properties.py +++ b/moldesign/_tests/test_atom_bond_computed_properties.py @@ -46,22 +46,22 @@ def test_bond_ffterms_returns_none_if_no_ff(h2): assert bond.ff is None -def test_basis_function_atom_access(h2_rhfwfn): - mol = h2_rhfwfn +def test_basis_function_atom_access(h2_rhf_sto3g): + mol = h2_rhf_sto3g for atom in mol.atoms: assert len(atom.basis_functions) == 1 # good ol' sto-3g assert len(atom.basis_functions[0].primitives) == 3 @pytest.mark.screening -def test_atom_property_access_to_mulliken_charges(h2_rhfwfn): - mol = h2_rhfwfn +def test_atom_property_access_to_mulliken_charges(h2_rhf_sto3g): + mol = h2_rhf_sto3g for atom in mol.atoms: assert abs(atom.properties.mulliken) <= 1e-5 * u.q_e -def test_atomic_forces(h2_rhfwfn): - mol = h2_rhfwfn +def test_atomic_forces(h2_rhf_sto3g): + mol = h2_rhf_sto3g helpers.assert_almost_equal(mol.atoms[0].force, -mol.atoms[1].force) diff --git a/moldesign/_tests/test_atoms.py b/moldesign/_tests/test_atoms.py index 8efc4d5..4e9c617 100644 --- a/moldesign/_tests/test_atoms.py +++ b/moldesign/_tests/test_atoms.py @@ -1,6 +1,9 @@ import pytest +import numpy as np import moldesign as mdt +from moldesign import units as u + from .molecule_fixtures import * from . import helpers diff --git a/moldesign/_tests/test_copies.py b/moldesign/_tests/test_copies.py index 2ef0efd..45bb5bf 100644 --- a/moldesign/_tests/test_copies.py +++ b/moldesign/_tests/test_copies.py @@ -233,11 +233,11 @@ def test_constraints_copied_with_molecule(mol_with_zerocharge_params): assert atom.molecule is mcpy -def test_properties_copied_with_molecule(cached_h2_rhfwfn): - original = cached_h2_rhfwfn +def test_properties_copied_with_molecule(cached_h2_rhf_sto3g): + original = cached_h2_rhf_sto3g assert original.potential_energy is not None # sanity check - mol = cached_h2_rhfwfn.copy() + mol = cached_h2_rhf_sto3g.copy() assert mol is not original assert mol.properties is not original.properties @@ -255,8 +255,8 @@ def test_properties_copied_with_molecule(cached_h2_rhfwfn): assert mol.properties[prop] is not val -def test_wfn_copied_with_molecule(cached_h2_rhfwfn): - original = cached_h2_rhfwfn +def test_wfn_copied_with_molecule(cached_h2_rhf_sto3g): + original = cached_h2_rhf_sto3g assert original.wfn is not None # sanity check mol = original.copy() @@ -269,8 +269,8 @@ def test_wfn_copied_with_molecule(cached_h2_rhfwfn): assert mol.wfn.aobasis.fock is not original.wfn.aobasis.fock -def test_wfn_copy(cached_h2_rhfwfn): - original = cached_h2_rhfwfn +def test_wfn_copy(cached_h2_rhf_sto3g): + original = cached_h2_rhf_sto3g wfn = original.wfn.copy() assert wfn.mol is original diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 77bb9e5..db80898 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -3,22 +3,25 @@ import random +import itertools import numpy as np -import numpy.testing as npt import pytest import moldesign from moldesign import units as u +from moldesign.mathutils import spherical_harmonics as harmonics + +from . import helpers +from .molecule_fixtures import * registered_types = {} +__PYTEST_MARK__ = ['math', 'gaussians'] -# TODO: implement and test spherical gaussians (not yet implemented) def typedfixture(*types, **kwargs): """This is a decorator that lets us associate fixtures with one or more arbitrary types. We'll later use this type to determine what tests to run on the result""" - def fixture_wrapper(func): for t in types: registered_types.setdefault(t, []).append(func.__name__) @@ -27,73 +30,277 @@ def fixture_wrapper(func): return fixture_wrapper -@typedfixture('gaussian') +@pytest.fixture def std_1d_gaussian(): g = moldesign.orbitals.gaussians.Gaussian([0.0]*u.angstrom, 1.0/u.angstrom ** 2) return g -@typedfixture('gaussian') +@typedfixture('basis_fn') def std_3d_gaussian(): g = moldesign.orbitals.gaussians.Gaussian([0.0, 0.0, 0.0]*u.angstrom, 1.0/u.angstrom ** 2) return g -@typedfixture('gaussian') -def std_10d_gaussian(): - g = moldesign.orbitals.gaussians.Gaussian([0.0 for i in range(10)]*u.angstrom, - 1.0/u.angstrom ** 2) +@typedfixture('basis_fn') +def cartesian_3d_gaussian(): + g = moldesign.orbitals.CartesianGaussian( + center=[random.random() for i in range(3)]*u.angstrom, + powers=[1, 3, 0], + alpha=1.1/u.angstrom ** 2, + coeff=1.0) return g -@typedfixture('cartesiangaussian') -def cart_10d_gaussian(): - g = moldesign.orbitals.gaussians.CartesianGaussian( - center=[random.random() for i in range(10)]*u.angstrom, - powers=[random.choice([0, 1, 2, 3]) for i in range(10)], - exp=1.1/u.angstrom ** 2) - return g +@typedfixture('basis_fn') +def spherical_3d_gaussian(): + g = moldesign.orbitals.SphericalGaussian( + center=[random.random() for i in range(3)]*u.angstrom, + l=3, m=-2, + alpha=1.1/u.angstrom ** 2, + coeff=1.0) + return g -@pytest.mark.parametrize('objkey', registered_types['gaussian']) +@pytest.mark.parametrize('objkey', ['std_1d_gaussian','std_3d_gaussian']) @pytest.mark.screening def test_gaussian_integral_and_dimensionality(objkey, request): - g = request.getfuncargvalue(objkey) - assert g.ndims == len(g.center) + g = request.getfixturevalue(objkey) + assert g.ndim == len(g.center) intval = g.integral - expectval = g.coeff*(np.pi/g.exp) ** (g.ndims/2.0) + expectval = g.coeff*(np.pi/g.alpha) ** (g.ndim/2.0) _assert_almost_equal(intval, expectval, decimal=10) -@pytest.mark.parametrize('objkey', - registered_types['gaussian'] + registered_types['cartesiangaussian']) +@pytest.fixture +def linear_combination(): + return _make_rando_linear_combination(True) + + +def _make_rando_gaussian(withunits=True): + if withunits: + length = u.angstrom + else: + length = 1.0 + return moldesign.orbitals.Gaussian((np.random.rand(3)-0.5)*1.0 * length, + (random.random()*5)/(length ** 2), + coeff=random.random()) + + +def _make_rando_cartesian_gaussian(powers, withunits=True): + if withunits: + length = u.angstrom + else: + length = 1.0 + return moldesign.orbitals.CartesianGaussian((np.random.rand(3)-0.5)*1.0 * length, + (random.random()*5)/(length ** 2), + powers=powers, + coeff=random.random()) + + +def _make_rando_spherical_gaussian(l,m, withunits=True): + if withunits: + length = u.angstrom + else: + length = 1.0 + return moldesign.orbitals.SphericalGaussian((np.random.rand(3)-0.5)*1.0 * length, + (random.random()*5)/(length ** 2), + l,m, + coeff=random.random()) + + +def _make_rando_linear_combination(withunits=True): + gaussians = [] + if withunits: + length = u.angstrom + else: + length = 1.0 + + center = (np.random.rand(3)-0.5)*1.0 * length + + for pwr in [(0,0,0), (1,1,1), (3,2,1)]: + gaussians.append( + moldesign.orbitals.CartesianGaussian( + center=center, + powers=pwr, + alpha=(10.0 * (random.random()+3))/(length**2), + coeff=1/(np.sqrt(3.0)))) + lc = moldesign.orbitals.PrimitiveSum(gaussians) + lc.ndims = 3 # so it works with the test suite + return lc + + +@pytest.mark.parametrize('withunits', [True, False], ids=['quantity','number']) +def test_numerical_vs_analytical_overlap_gauss(withunits): + p1 = _make_rando_gaussian(withunits) + p2 = _make_rando_gaussian(withunits) + _assert_numerical_analytical_overlaps_match(p1, p2) + + +@pytest.mark.parametrize('withunits', [True, False], ids=['quantity','number']) +def test_numerical_vs_analytical_overlap_cartesian(withunits): + p1 = _make_rando_cartesian_gaussian((1,2,3), withunits) + p2 = _make_rando_cartesian_gaussian((1,0,1), withunits) + _assert_numerical_analytical_overlaps_match(p1, p2) + + +@pytest.mark.parametrize('withunits', [True, False], ids=['quantity','number']) +def test_numerical_vs_analytical_overlap_spherical(withunits): + p1 = _make_rando_spherical_gaussian(1,-1, withunits) + p2 = _make_rando_spherical_gaussian(2,0, withunits) + _assert_numerical_analytical_overlaps_match(p1, p2) + + +@pytest.mark.parametrize('withunits', [True, False], ids=['quantity','number']) +def test_numerical_vs_analytical_overlap_linear_combination(withunits): + p1 = _make_rando_linear_combination(withunits) + p2 = _make_rando_linear_combination(withunits) + _assert_numerical_analytical_overlaps_match(p1, p2) + + +def _assert_numerical_analytical_overlaps_match(g1, g2): + olap = g1.overlap(g2) + try: + prod = g1*g2 + except NotImplementedError: + assert isinstance(g1, moldesign.orbitals.SphericalGaussian) + assert isinstance(g2, moldesign.orbitals.SphericalGaussian) + else: + helpers.assert_almost_equal(prod.integral, olap) + + allpoints, grid = helpers.generate_grid(g1, g2) + with np.errstate(under='ignore'): + prodvals = g1(allpoints) * g2(allpoints) + numsum = prodvals.sum() * grid.dx * grid.dy * grid.dz + helpers.assert_almost_equal(numsum, olap, decimal=5) + + +@pytest.mark.parametrize('withunits', [False, True]) +def test_gaussian_multiplication_amplitudes(withunits): + g1 = _make_rando_gaussian(withunits) + g2 = _make_rando_gaussian(withunits) + _assert_same_function_values(g1, g2, withunits) + + +# parameterizations across a sample of cartesian gaussians +test_powers = ((0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0), (1,1,1), (2,0,2), (4,1,1)) +cartesian_test_suite = list(itertools.product(test_powers, test_powers, [True, False])) +cartesian_test_ids = ['[%d%d%d]*[%d%d%d]/%s' % (p[0] + p[1] + ('units' if p[2] else 'c-num',)) + for p in cartesian_test_suite] + + +@pytest.mark.parametrize('p1,p2,withunits', + cartesian_test_suite, + ids=cartesian_test_ids) +def test_cartesian_gaussian_multiplication_amplitudes(p1, p2, withunits): + """ Tests that ``g1(x) * g2(x) == (g1 * g2)(x)`` + """ + g1 = _make_rando_cartesian_gaussian(p1, withunits) + g2 = _make_rando_cartesian_gaussian(p2, withunits) + _assert_same_function_values(g1, g2, withunits) + + +def _assert_same_function_values(g1, g2, withunits): + testcoords = 6.0*(np.random.rand(50, 3)-0.5) + if withunits: + testcoords = testcoords*u.angstrom + g1g2 = g1*g2 + gvals = g1g2(testcoords) + g1vals = g1(testcoords) + g2vals = g2(testcoords) + prodvals = g1vals*g2vals + helpers.assert_almost_equal(prodvals, gvals) + + +def test_initial_gaussian_normalization_gaussian(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + g2 = moldesign.orbitals.Gaussian(center, exp, normalized=True) + helpers.assert_almost_equal(1.0, _numerical_norm(g2), decimal=3) + helpers.assert_almost_equal(1.0, g2.norm) + + +def test_initial_gaussian_normalization_with_prefactor(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + g1 = moldesign.orbitals.Gaussian(center, exp, coeff=3.0*u.angstrom, normalized=True) + helpers.assert_almost_equal(3.0*u.angstrom, _numerical_norm(g1), decimal=3) + helpers.assert_almost_equal(3.0*u.angstrom, g1.norm) + + +def test_initial_normalization_cartesian(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + for powers in itertools.product(range(4), range(4), range(4)): + g2 = moldesign.orbitals.CartesianGaussian(center, exp, powers, normalized=True) + helpers.assert_almost_equal(1.0, _numerical_norm(g2), decimal=3) + helpers.assert_almost_equal(1.0, g2.norm) + + +def test_initial_normalization_cartesian_with_prefactor(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + for powers in itertools.product(range(4), range(4), range(4)): + g1 = moldesign.orbitals.CartesianGaussian(center, exp, powers, coeff=3.0, normalized=True) + helpers.assert_almost_equal(3.0, _numerical_norm(g1), decimal=3) + helpers.assert_almost_equal(3.0, g1.norm) + + +def test_initial_normalization_spherical(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + for l in range(5): + for m in range(-l, l+1): + g2 = moldesign.orbitals.SphericalGaussian(center, exp, l, m, normalized=True) + helpers.assert_almost_equal(1.0, _numerical_norm(g2), decimal=3) + helpers.assert_almost_equal(1.0, g2.norm) + + +def test_initial_normalization_spherical_with_prefactor(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + for l in range(5): + for m in range(-l, l+1): + g1 = moldesign.orbitals.SphericalGaussian(center, exp, l, m, + coeff=3.0 * u.angstrom, normalized=True) + helpers.assert_almost_equal(3.0 * u.angstrom, _numerical_norm(g1), decimal=3) + helpers.assert_almost_equal(3.0 * u.angstrom, g1.norm) + + +def _numerical_norm(g): + allpoints, grid = helpers.generate_grid(g) + with np.errstate(under='ignore'): + vals = g(allpoints) + numnorm = np.sqrt(grid.dx * grid.dy * grid.dz * (vals**2).sum()) + return numnorm + + +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_gaussian_function_values(objkey, request): - g = request.getfuncargvalue(objkey) + g = request.getfixturevalue(objkey) for idim in range(g.ndims): coord = g.center.copy() - randoffset = 4.0 * (random.random() - 0.5) * g.exp**-0.5 + randoffset = 4.0 * (random.random() - 0.5) * g.alpha**-0.5 coord[idim] += randoffset funcval = _gfuncval(g, coord) retval = g(coord) _assert_almost_equal(funcval, retval) -@pytest.mark.parametrize('objkey', - registered_types['gaussian'] + registered_types['cartesiangaussian']) -@pytest.mark.screening +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_vectorized_gaussian_function_evaluations(objkey, request): - g = request.getfuncargvalue(objkey) + g = request.getfixturevalue(objkey) coords = np.zeros((5, g.ndims)) * g.center.units for i in range(5): coords[i] = g.center - randoffset = 4.0 * (random.random() - 0.5) * g.exp**-0.5 + randoffset = 4.0 * (random.random() - 0.5) * g.alpha**-0.5 idim = random.randrange(g.ndims) coords[i, idim] += randoffset @@ -105,24 +312,75 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): _assert_almost_equal(vector_results, expected, decimal=8) -def test_gaussian_self_overlap_is_unity(): - g1 = moldesign.orbitals.gaussians.Gaussian([0.0, 0.0, 0.0]*u.angstrom, - 1.0/u.angstrom ** 2, - coeff=-1.0) - g2 = moldesign.orbitals.gaussians.Gaussian([0.0, 0.0, 0.0]*u.angstrom, - 1.0/u.angstrom ** 2, - coeff=123.456) - npt.assert_almost_equal(-1.0, - g1.overlap(g2, normalized=True)) +@pytest.mark.parametrize('objkey', registered_types['basis_fn'] + ['linear_combination']) +def test_gaussian_str_and_repr_works(objkey, request): + g1 = request.getfixturevalue(objkey) + str(g1) + repr(g1) + + +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) +def test_normalized_gaussian_self_overlap_is_unity(objkey, request): + g1 = request.getfixturevalue(objkey) + g2 = g1.copy() + g1.coeff = -10.0 + g2.coeff = 12341.1832 + olap = g1.overlap(g2, normalized=True) + assert abs(-1.0 - olap) < 1e-12 + + g1.coeff = 10.0 + olap = g1.overlap(g2, normalized=True) + assert abs(1.0 - olap) < 1e-12 + + +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) +def test_normalization(objkey, request): + g1 = request.getfixturevalue(objkey) + oldnorm = g1.norm + + g1.coeff = (random.random() - 0.5) * 428.23 + try: + assert g1.norm != oldnorm + except u.DimensionalityError: + pass # this is a reasonable thing to happen too + + g1.normalize() + assert abs(g1.norm - 1.0) < 1e-12 + + +def test_linear_combination_normalization(linear_combination): + g1 = linear_combination + oldnorm = g1.norm + + prefactor = (random.random() - 0.5) * 428.23 + for prim in g1: + prim.coeff *= prefactor + + try: + assert g1.norm != oldnorm + except u.DimensionalityError: + pass # this is a reasonable thing to happen too + + g1.normalize() + assert abs(g1.norm - 1.0) < 1e-12 def _gfuncval(g, coord): - fv = g.coeff * np.exp(-g.exp * np.sum((g.center - coord)**2)) - for r, r0, pow in zip(coord, g.center, g.powers): - if pow != 0: - fv *= (r-r0)**pow + r = g.center - coord + if len(coord.shape) > 1: + r2 = np.sum(r**2, axis=1) + else: + r2 = np.sum(r**2) + fv = g.coeff * np.exp(-g.alpha * r2) + if isinstance(g, moldesign.orbitals.SphericalGaussian): + fv *= r2**(g.l/2.0) * harmonics.Y(g.l, g.m)(coord - g.center) + elif isinstance(g, moldesign.orbitals.CartesianGaussian): # assume cartesian + for r, r0, pow in zip(coord, g.center, g.powers): + if pow != 0: + fv *= (r-r0)**pow return fv + def _assert_almost_equal(a, b, *args, **kwargs): a_is_quantity = hasattr(a,'units') b_is_quantity = hasattr(b,'units') @@ -138,8 +396,102 @@ def _assert_almost_equal(a, b, *args, **kwargs): *args, **kwargs) -def test_cart_to_powers(): +def test_convert_cartesian_label_to_array_of_integer_powers(): from moldesign.orbitals.gaussians import cart_to_powers assert cart_to_powers('y') == [0, 1, 0] assert cart_to_powers('xxyz') == [2, 1, 1] - assert cart_to_powers('zx^3') == [3,0,1] \ No newline at end of file + assert cart_to_powers('zx^3') == [3,0,1] + + +@pytest.mark.parametrize('key', registered_types['basis_fn'] + ['linear_combination']) +def test_numerical_vs_analytical_norm(key, request): + g = request.getfixturevalue(key) + numnorm = _numerical_norm(g) + helpers.assert_almost_equal(g.norm, numnorm) + + +@pytest.mark.screening +def test_s_orbitals_equivalent_among_gaussian_types(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + + g_bare = moldesign.orbitals.Gaussian(center, exp) + g_cart = moldesign.orbitals.CartesianGaussian(center, exp, [0,0,0]) + g_sphr = moldesign.orbitals.SphericalGaussian(center, exp, 0, 0) + + for gauss in (g_bare, g_cart, g_sphr): + # normalize to amplitude of 1.0 at center + gauss.coeff = gauss.coeff / gauss(center) + assert gauss(center) == 1.0 + + _assert_orbitals_equivalent(g_bare, g_cart) + _assert_orbitals_equivalent(g_bare, g_sphr) + + +# real spherical harmonics that can be represented as a single cartesian term: +LM_TO_CART = {(1,-1): (0,1,0), + (1,0): (0,0,1), + (1,1): (1,0,0), + (2,-2): (1,1,0), + (2,-1): (0,1,1), + (2,1): (1,0,1), + (3,-2): (1,1,1)} + + +@pytest.mark.parametrize('lm,powers', + LM_TO_CART.items(), + ids=['lm:%d,%d, xyz:%d%d%d' % (args[0] + args[1]) + for args in LM_TO_CART.items()]) +def test_orbitals_same_in_cartesian_and_spherical(lm, powers): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + + g_cart = moldesign.orbitals.CartesianGaussian(center, exp, powers) + g_sphr = moldesign.orbitals.SphericalGaussian(center, exp, *lm) + + _assert_orbitals_equivalent(g_cart, g_sphr) + + +@pytest.mark.parametrize('l', range(4), ids=lambda x:'l=%s' % x) +def test_spherical_to_cartesian(l): + for m in range(-l,l+1): + center = np.random.rand(3)*u.angstrom + exp = random.random()*2.0/u.angstrom ** 2 + + bf = moldesign.orbitals.SphericalGaussian(center, exp, l, m, normalized=True) + _assert_orbitals_equivalent(bf, bf.to_cart()) + + +def _assert_orbitals_equivalent(g1, g2): + helpers.assert_almost_equal(g1.norm, + g2.norm) + + testcoords = 6.0*(np.random.rand(50, 3)-0.5)*u.angstrom + helpers.assert_almost_equal(g1(testcoords), + g2(testcoords)) + + +def test_pyscf_and_mdt_norms_are_the_same(h2_rhf_augccpvdz): + mol = h2_rhf_augccpvdz + basis = mol.wfn.aobasis + + for bf in basis: + assert abs(bf.norm - 1.0) < 1e-12 + + +def test_pyscf_and_mdt_overlaps_are_the_same(h2_rhf_augccpvdz): + mol = h2_rhf_augccpvdz + basis = mol.wfn.aobasis + + calc_overlap_mat = [] + for i in range(len(basis)): + calc_overlap_mat.append( + [basis[i].overlap(basis[j]) for j in range(len(basis))] + ) + + overlaps = u.array(calc_overlap_mat) + assert isinstance(overlaps, np.ndarray) or overlaps.units == u.dimensionless + + np.testing.assert_allclose(mol.wfn.aobasis.overlaps, + overlaps, + atol=5.0e-7) diff --git a/moldesign/_tests/test_mathutils.py b/moldesign/_tests/test_mathutils.py index de63195..6d13dd2 100644 --- a/moldesign/_tests/test_mathutils.py +++ b/moldesign/_tests/test_mathutils.py @@ -1,13 +1,16 @@ +import itertools import pytest import numpy as np from moldesign import mathutils from moldesign import units as u +from moldesign.mathutils import spherical_harmonics as sh + from . import helpers -__PYTEST_MARK__ = 'internal' +__PYTEST_MARK__ = ['internal', 'math', 'mathutils'] CONSTRUCTORS = [] IDS = [] @@ -97,6 +100,145 @@ def test_perpendicular(setupfn): assert np.abs(1.0 - mathutils.norm(perpvec)) < 1e-12 +def test_spherical_angles(): + assert sh.cart_to_polar_angles(np.zeros(3)) == (0.0, 0.0) + + oneangles = sh.cart_to_polar_angles(np.array([1.0, 0.0, 1.0])) + np.testing.assert_allclose(np.array(oneangles), + np.array([np.pi/4, 0.0])) + + oneangles = sh.cart_to_polar_angles(np.array([0.0, 1.0, 1.0])) + np.testing.assert_allclose(np.array(oneangles), + np.array([np.pi/4, np.pi/2.0])) + + np.testing.assert_allclose(np.array(sh.cart_to_polar_angles(RANDARRAY)), + np.array(RAND_ARRAY_ANGLES).T) + + +@pytest.mark.parametrize('shell', range(4)) +def test_cart_and_spherical_equivalence(shell): + l = shell + + testpoints = (np.random.rand(10, 3)-0.5)*10.0 + + for m in range(-l,l+1): + real_y = sh.Y(l, m) + cart_y = sh.SPHERE_TO_CART[l,m] + + # single point + np.testing.assert_allclose(real_y(testpoints[0]), + cart_y(testpoints[0]), + err_msg="Failed at %s" % ((l,m),)) + + # vectorized + np.testing.assert_allclose(real_y(testpoints), + cart_y(testpoints), + err_msg="Failed at %s" % ((l,m),)) + + +@pytest.mark.parametrize('l1,m1,l2,m2', + [(0,0,0,0), + (1,0,1,0), + (3,-2,2,-1), + (1,0,1,1), + (1,-1,1,1), + (10,5,10,5)]) +def test_spherical_harmonics_orthonormal(l1, m1, l2, m2): + y1 = sh.Y(l1, m1) + y2 = sh.Y(l2, m2) + + NPHI = 1000 + NTHETA = 700 + + thetas = np.linspace(0, np.pi, num=NTHETA) + dtheta = np.pi/NTHETA + phis = np.linspace(0, 2.0 * np.pi, num=NPHI) + cosphi = np.cos(phis) + sinphi = np.sin(phis) + dphi = 2.0 * np.pi/NPHI + + phibuffer = np.zeros((NPHI, 3)) + integ = 0.0 + for theta in thetas: + sintheta = np.sin(theta) + phibuffer[:,0] = cosphi * sintheta + phibuffer[:,1] = sinphi * sintheta + phibuffer[:,2] = np.cos(theta) + values = y1(phibuffer) * y2(phibuffer) + integ += dphi * dtheta * sintheta * values.sum() + + if l1 == l2 and m1 == m2: + assert abs(integ - 1.0) < 1e-2 + else: + assert abs(integ) < 1e-2 + + +@pytest.fixture +def ndarray_ranges(): + ranges = np.array([[1,4], + [10, 40], + [100, 400]]) + return ranges + + +@pytest.fixture +def ranges_with_units(ndarray_ranges): + return ndarray_ranges * u.angstrom + + +@pytest.mark.parametrize('key', ['ndarray_ranges', 'ranges_with_units']) +def test_volumetric_grid_point_list(key, request): + ranges = request.getfixturevalue(key) + grid = mathutils.VolumetricGrid(*ranges, xpoints=3, ypoints=4, zpoints=5) + assert (grid.xpoints, grid.ypoints, grid.zpoints) == (3,4,5) + pl = list(grid.iter_points()) + pa = grid.allpoints() + assert (u.array(pl) == pa).all() + + for i in range(3): + assert pa[:,i].min() == ranges[i][0] + assert pa[:,i].max() == ranges[i][1] + + assert grid.npoints == 3*4*5 + assert len(pl) == grid.npoints + + for idim in range(3): + for ipoint in range(1,grid.points[idim]): + helpers.assert_almost_equal(grid.spaces[idim][ipoint] - grid.spaces[idim][ipoint-1], + grid.deltas[idim]) + + +@pytest.mark.parametrize('key', ['ndarray_ranges', 'ranges_with_units']) +def test_volumetric_iteration(key, request): + ranges = request.getfixturevalue(key) + grid = mathutils.VolumetricGrid(*ranges, npoints=4) + + grid_iterator = grid.iter_points() + assert (grid.xpoints, grid.ypoints, grid.zpoints) == (4,4,4) + assert grid.npoints == 4**3 + + for ix,x in enumerate(grid.xspace): + for iy, y in enumerate(grid.yspace): + for iz, z in enumerate(grid.zspace): + gridpoint = next(grid_iterator) + if ix == iy == iz == 0: + assert (gridpoint == grid.origin).all() + assert x == gridpoint[0] + assert y == gridpoint[1] + assert z == gridpoint[2] + + +RANDARRAY = np.array([[ 0.03047527, -4.49249374, 5.83154878], + [-7.38885486, 8.30076569, -7.97671261], + [ 6.79409582, -4.07664079, 8.45983794], + [-4.25045695, 0.03411114, 9.02986275]]) + +RAND_ARRAY_ANGLES = [(0.656426762, -1.564012833), + (2.1933588598416, 2.2981378887583), + (0.75266068344496, -0.54043933754409), + (0.43995560467481, 3.1335675381452)] + + TESTARRAY = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], diff --git a/moldesign/_tests/test_minimizers.py b/moldesign/_tests/test_minimizers.py index 2977c44..03cb913 100644 --- a/moldesign/_tests/test_minimizers.py +++ b/moldesign/_tests/test_minimizers.py @@ -8,6 +8,8 @@ from moldesign import units as u +__PYTEST_MARK__ = 'minimization' + @pytest.fixture(scope='function') def harmonic_atom(): mol = mdt.Molecule([mdt.Atom(1)]) diff --git a/moldesign/_tests/test_molecules.py b/moldesign/_tests/test_molecules.py index b50355a..133cf8c 100644 --- a/moldesign/_tests/test_molecules.py +++ b/moldesign/_tests/test_molecules.py @@ -3,9 +3,11 @@ import pickle from past.builtins import basestring -import moldesign.exceptions -import moldesign.molecules.atomcollections -import moldesign.utils.classes +import pytest +import numpy as np + +from moldesign import units as u +import moldesign as mdt from .object_fixtures import * from .molecule_fixtures import * @@ -15,7 +17,6 @@ __PYTEST_MARK__ = 'internal' # mark all tests in this module with this label (see ./conftest.py) - def test_h2_protected_atom_arrays(h2): atom1, atom2 = h2.atoms with pytest.raises(TypeError): @@ -86,7 +87,7 @@ def test_h2_cache_flush(h2_harmonic): h2 = h2_harmonic pe = h2.calc_potential_energy() f = h2.forces - h2.atoms[0].x += 0.1285*u.ang + h2.atoms[0].x += 0.1285*u.angstrom pe2 = h2.calc_potential_energy() f2 = h2.forces assert pe != pe2 @@ -95,19 +96,19 @@ def test_h2_cache_flush(h2_harmonic): def test_h2_not_calculated_yet(h2_harmonic): h2_harmonic.calculate() - h2_harmonic.atoms[1].x += 0.3*u.ang - with pytest.raises(moldesign.exceptions.NotCalculatedError): + h2_harmonic.atoms[1].x += 0.3*u.angstrom + with pytest.raises(mdt.NotCalculatedError): h2_harmonic.forces - with pytest.raises(moldesign.exceptions.NotCalculatedError): + with pytest.raises(mdt.NotCalculatedError): h2_harmonic.potential_energy def h2_properties_raises_not_calculated_yet(h2_harmonic): h2_harmonic.calculate() - h2_harmonic.atoms[1].x += 0.3*u.ang - with pytest.raises(moldesign.exceptions.NotCalculatedError): + h2_harmonic.atoms[1].x += 0.3*u.angstrom + with pytest.raises(mdt.NotCalculatedError): h2_harmonic.properties.forces - with pytest.raises(moldesign.exceptions.NotCalculatedError): + with pytest.raises(mdt.NotCalculatedError): h2_harmonic.properties.potential_energy @@ -219,7 +220,7 @@ def test_molecule_bonds(molkey, request): @pytest.mark.parametrize('molkey', registered_types['molecule']) def test_molecule_types(molkey, request): mol = request.getfixturevalue(molkey) - assert issubclass(type(mol.atoms), moldesign.molecules.atomcollections.AtomList) + assert issubclass(type(mol.atoms), mdt.AtomList) for atom in mol.atoms: assert issubclass(type(atom), mdt.Atom) for residue in mol.residues: @@ -408,7 +409,7 @@ def test_initialization_charges(): assert mol.charge == -1 * u.q_e with pytest.raises(TypeError): - mdt.Atom('H', charge=3) # it needs to be "formal_charge", to distinguish from partial charge + mdt.Atom('H', charge=3) # it needs to be "formal_charge" to distinguish from partial charge m2 = mdt.Molecule([a1], charge=-1) assert m2.charge == -1 * u.q_e diff --git a/moldesign/_tests/test_objects.py b/moldesign/_tests/test_objects.py index 81ce74f..945e5ef 100644 --- a/moldesign/_tests/test_objects.py +++ b/moldesign/_tests/test_objects.py @@ -1,8 +1,12 @@ import pickle +import pytest +import numpy as np + from moldesign.utils import Alias from .object_fixtures import * from .molecule_fixtures import * +from .object_fixtures import TESTDICT __PYTEST_MARK__ = 'internal' # mark all tests in this module with this label (see ./conftest.py) diff --git a/moldesign/_tests/test_qm_xfaces.py b/moldesign/_tests/test_qm_xfaces.py index decdcb0..f01a53d 100644 --- a/moldesign/_tests/test_qm_xfaces.py +++ b/moldesign/_tests/test_qm_xfaces.py @@ -49,26 +49,6 @@ def test_minimization_trajectory(h2_with_model): helpers.assert_something_resembling_minimization_happened(p1, e1, traj, mol) -@typedfixture('model') -def pyscf_rhf(): - return mdt.models.PySCFPotential(basis='sto-3g', theory='rhf') - - -@typedfixture('model') -def pyscf_dft(): - return mdt.models.PySCFPotential(basis='sto-3g', theory='rks', functional='b3lyp') - - -@typedfixture('model') -def nwchem_rhf(): - return mdt.models.NWChemQM(basis='sto-3g', theory='rhf', functional='b3lyp') - - -@typedfixture('model') -def nwchem_dft(): - return mdt.models.NWChemQM(basis='sto-3g', theory='rks', functional='b3lyp') - - @pytest.mark.parametrize('objkey', TESTSET) def test_pyscf_rhf_sto3g_properties(objkey, request): mol = request.getfixturevalue(objkey) @@ -139,9 +119,9 @@ def test_calc_eri_tensor(h2): @pytest.mark.screening -def test_aobasis(h2_rhfwfn): +def test_aobasis(h2_rhf_sto3g): # it's sto-3g, so structure is simple - aobasis = h2_rhfwfn.wfn.aobasis + aobasis = h2_rhf_sto3g.wfn.aobasis assert aobasis.basisname == 'sto-3g' assert (aobasis.coeffs == np.identity(2)).all() np.testing.assert_allclose(aobasis.fock, aobasis.fock.T) @@ -149,7 +129,7 @@ def test_aobasis(h2_rhfwfn): assert aobasis.fock.dimensionality == u.eV.dimensionality - for orb in h2_rhfwfn.wfn.aobasis: + for orb in h2_rhf_sto3g.wfn.aobasis: assert orb.aotype == '1s' assert orb.orbtype == 's' assert len(orb.primitives) == 3 diff --git a/moldesign/_tests/test_trajectory.py b/moldesign/_tests/test_trajectory.py index c1e34b2..b8d7e9e 100644 --- a/moldesign/_tests/test_trajectory.py +++ b/moldesign/_tests/test_trajectory.py @@ -141,7 +141,7 @@ def h2_wfn_traj(h2): def test_align_phases(h2_wfn_traj): - # TODO: actually check results + # TODO: actually check that the orbitals were aligned h2_wfn_traj.align_orbital_phases() h2_wfn_traj.align_orbital_phases(reference_frame=1) diff --git a/moldesign/_tests/test_units.py b/moldesign/_tests/test_units.py index 1b69521..172aa7d 100644 --- a/moldesign/_tests/test_units.py +++ b/moldesign/_tests/test_units.py @@ -6,7 +6,8 @@ from moldesign import units from moldesign import units as u -__PYTEST_MARK__ = 'internal' # mark all tests in this module with this label (see ./conftest.py) +# mark all tests in this module with these label (see ./conftest.py) +__PYTEST_MARK__ = ['internal', 'units'] @pytest.mark.screening @@ -57,6 +58,41 @@ def test_default_units(): assert units.default.energy == units.eV +@pytest.fixture +def si_units(): + system = units.UnitSystem(length=u.ureg.meters, mass=u.ureg.kg, + time=u.ureg.seconds, energy=u.ureg.joules) + return system + + +@pytest.mark.parametrize(['dimension', 'expected', 'weird_units'], + [['length', u.ureg.meters, u.ureg.lightyear], + ['mass', u.ureg.kg, u.ureg.stone], + ['time', u.ureg.seconds, u.ureg.year], + ['momentum', u.ureg.kg * u.ureg.meter / u.ureg.second, + u.ureg.lightyear * u.ureg.stone / u.ureg.year], + ['force', u.ureg.joule / u.ureg.meter, u.ureg.force_ton], + ['energy', u.ureg.joule, u.ureg.horsepower * u.ureg.seconds]], + ids='length mass time momentum force energy'.split()) +def test_si_unit_system(si_units, dimension, expected, weird_units): + assert getattr(si_units, dimension) == expected + + default_unit = getattr(u.default, dimension) + my_quantity = 3.141 * default_unit + si_quantity = si_units.convert(my_quantity) + + assert my_quantity.units != expected + assert si_quantity.units == expected + assert u.default.convert(abs(my_quantity - si_quantity)) < 1e-15 * default_unit + + setattr(si_units, dimension, weird_units) + assert getattr(si_units, dimension) == weird_units + + weird_quantity = si_units.convert(my_quantity) + assert abs(my_quantity-weird_quantity) < 1e-15 * default_unit + assert weird_quantity.units == weird_units + + @pytest.mark.screening def test_default_unit_conversions(): my_list = [1.0*units.angstrom, 1.0*units.nm, 1.0*units.a0] @@ -68,6 +104,47 @@ def test_default_unit_conversions(): np.testing.assert_almost_equal(result[2], 0.05291772, 7) +def test_dimensional_array_from_mixed_data(): + data = [np.arange(3) * u.angstrom, + [1.0 * u.nm, 1.0 * u.nm**2/u.angstrom, 3.0*u.ureg.km], + [1,1,1] * u.bohr] + arr = u.array(data) + + bohr_in_ang = (1.0 * u.bohr).value_in(u.angstrom) + expected = u.angstrom * [[0,1,2], + [10.0, 100.0, 3.0e13], + [bohr_in_ang, bohr_in_ang, bohr_in_ang]] + np.testing.assert_allclose(arr, expected) + assert isinstance(arr, units.MdtQuantity) + + +def test_dimensionless_array_from_mixed_data(): + data = [np.arange(3), + [1.0 * u.dimensionless, 1.0 * u.nm/u.angstrom, 3.0], + [1,1,1] * u.ureg.kg / u.ureg.g] + arr = u.array(data) + np.testing.assert_allclose(arr, + np.array([[0,1,2], + [1,10,3], + [1000,1000,1000]],dtype='float') ) + assert isinstance(arr, np.ndarray) + + +def test_inconsistent_array_units_raises_dimensionality_error(): + with pytest.raises(units.DimensionalityError): + u.array([[1,2,3] * u.angstrom, [3*u.bohr, 4*u.eV, 5*u.bohr]]) + + +def test_mixed_units_and_nonunits_raises_dimensionality_error(): + with pytest.raises(units.DimensionalityError): + u.array([[1,2,3] * u.angstrom, [3*u.bohr, 4, 5*u.bohr]]) + + +def test_no_nonsquare_arrays(): + with pytest.raises(ValueError): + u.array([[1,2],[3]]) + + def test_scalar_comparisons_to_zero_ignore_units(): num = 1.0*units.angstrom assert num > 0.0 @@ -169,7 +246,7 @@ def test_default_unit_conversions(): def test_getunits_doctests(): - assert units.get_units(1.0*units.angstrom) == units.MdtUnit('ang') + assert units.get_units(1.0*units.angstrom) == units.MdtUnit('angstrom') assert units.get_units(np.array([1.0, 2, 3.0])) == units.MdtUnit('dimensionless') assert units.get_units([[1.0*units.dalton, 3.0*units.eV], ['a'], 'gorilla']) == units.MdtUnit('amu') diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py new file mode 100644 index 0000000..7239a0c --- /dev/null +++ b/moldesign/_tests/test_wfn.py @@ -0,0 +1,78 @@ +import pytest +import numpy as np + +import moldesign as mdt +from moldesign import units as u +from moldesign.interfaces.pyscf_interface import basis_values + +from .molecule_fixtures import * +from . import helpers + + +__PYTEST_MARK__ = 'wfn' + +TESTSYSTEMS = ['h2_rhf_augccpvdz', 'h2_rhf_sto3g', 'acetylene_dft_631g'] + + +@pytest.mark.parametrize('molkey', TESTSYSTEMS) +def test_pyscf_orbital_grid_works(molkey, request): + """ Tests the basic input/output of the pyscf basis_values function + + Doesn't actually test the values directly - just that the answers are mathematically consistent + """ + mol = request.getfixturevalue(molkey) + wfn = mol.wfn + nbasis = len(wfn.aobasis) + + coords = u.array([mol.com, + np.zeros(3)*u.angstrom, + 10.0 * np.ones(3) * u.angstrom, + np.ones(3)*u.nm]) + + # First - check that the shape is appropriate if called without orbital coefficients + values_nocoeffs = basis_values(mol, wfn.aobasis, coords) + assert values_nocoeffs.shape == (len(coords), nbasis) + assert (values_nocoeffs[-1] == values_nocoeffs[-2]).all() # these 2 coordinates are the same + + # Second - explicitly send orbital coefficients for first 2 basis functions + coeffs = np.zeros((2, nbasis)) + coeffs[:2, :2] = np.identity(2) + vals_b0 = basis_values(mol, wfn.aobasis, coords, coeffs=coeffs) + assert vals_b0.shape == (len(coords), len(coeffs)) + np.testing.assert_allclose(values_nocoeffs[:,:2], vals_b0) + + # Third - send symmetric and anti-symmetric combinations of basis functions and check answers + plusminus = np.zeros((2, nbasis)) + plusminus[:2, :2] = 1.0 / np.sqrt(2) + plusminus[1,1] = -1.0 / np.sqrt(2) + vals_plusminus = basis_values(mol, wfn.aobasis, coords, coeffs=plusminus) + assert vals_plusminus.shape == (len(coords), len(coeffs)) + helpers.assert_almost_equal(vals_plusminus[:,0], + (vals_b0[:,0] + vals_b0[:,1])/np.sqrt(2)) + helpers.assert_almost_equal(vals_plusminus[:,1], + (vals_b0[:,0] - vals_b0[:,1])/np.sqrt(2)) + + +@pytest.mark.parametrize('molkey', TESTSYSTEMS) +def test_basis_function_3d_grids_same_in_pyscf_and_mdt(molkey, request): + mol = request.getfixturevalue(molkey) + randocoords = 6.0 * u.angstrom * (np.random.rand(200, 3) - 0.5) + pyscf_vals = basis_values(mol, mol.wfn.aobasis, randocoords) + with np.errstate(under='ignore'): + mdt_vals = mol.wfn.aobasis(randocoords) + helpers.assert_almost_equal(mdt_vals, pyscf_vals, decimal=6) + + +@pytest.mark.parametrize('molkey', ['h2_rhf_augccpvdz', 'h2_rhf_sto3g']) +@pytest.mark.screening +def test_pyscf_basis_function_space_integral_normalized(molkey, request): + mol = request.getfixturevalue(molkey) + grid = mdt.mathutils.padded_grid(mol.positions, 8.0 * u.angstrom, npoints=150) + points = grid.allpoints() + + pyscf_vals = basis_values(mol, mol.wfn.aobasis, points) + assert pyscf_vals.shape == (len(points), len(mol.wfn.aobasis)) + + pyscf_normsqr = (pyscf_vals ** 2).sum(axis=0) * grid.dx * grid.dy * grid.dz + helpers.assert_almost_equal(pyscf_normsqr, 1.0, + decimal=4) diff --git a/moldesign/data/data.py b/moldesign/data/data.py index d6a37e4..f6853a0 100644 --- a/moldesign/data/data.py +++ b/moldesign/data/data.py @@ -38,6 +38,14 @@ def __getitem__(self, item): DEFAULT_FORCE_TOLERANCE = (0.0001 * u.hartree / u.bohr).defunits() # taken from GAMESS OPTTOL keyword +# Precompute odd factorial values (N!!) +ODD_FACTORIAL = {0: 1} # by convention +_ofact = 1 +for _i in range(1, 20, 2): + _ofact *= _i + ODD_FACTORIAL[_i] = float(_ofact) + + def print_environment(): """For reporting bugs - spits out the user's environment""" import sys diff --git a/moldesign/external/pyquante2/LICENSE b/moldesign/external/pyquante2/LICENSE new file mode 100644 index 0000000..719a0a0 --- /dev/null +++ b/moldesign/external/pyquante2/LICENSE @@ -0,0 +1,28 @@ +Copyright (c) 2013 PyQuante Development Team + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + a. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + b. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + c. Neither the name of SymPy nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. \ No newline at end of file diff --git a/moldesign/external/pyquante2/__init__.py b/moldesign/external/pyquante2/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/moldesign/external/pyquante2/one.py b/moldesign/external/pyquante2/one.py new file mode 100644 index 0000000..ce4904d --- /dev/null +++ b/moldesign/external/pyquante2/one.py @@ -0,0 +1,253 @@ +""" +One electron integrals. +""" + +from numpy import pi,exp,floor,array,isclose +from math import factorial + +## BEGIN MODIFIED CODE: +#from pyquante2.utils import binomial, fact2, Fgamma, norm2 +from .utils import binomial, fact2, Fgamma, norm2 +#END MODIFIED CODE + +# Notes: +# The versions S,T,V include the normalization constants +# The version overlap,kinetic,nuclear_attraction do not. +# This is so, for example, the kinetic routines can call the potential routines +# without the normalization constants getting in the way. + +def S(a,b): + """ + Simple interface to the overlap function. + >>> from pyquante2 import pgbf,cgbf + >>> s = pgbf(1) + >>> isclose(S(s,s),1.0) + True + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(S(sc,sc),1.0) + True + + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(S(sc,s),1.0) + True + >>> isclose(S(s,sc),1.0) + True + + """ + if b.contracted: + return sum(cb*S(pb,a) for (cb,pb) in b) + elif a.contracted: + return sum(ca*S(b,pa) for (ca,pa) in a) + return a.norm*b.norm*overlap(a.exponent,a.powers, + a.origin,b.exponent,b.powers,b.origin) + +def T(a,b): + """ + Simple interface to the kinetic function. + >>> from pyquante2 import pgbf,cgbf + >>> from pyquante2.basis.pgbf import pgbf + >>> s = pgbf(1) + >>> isclose(T(s,s),1.5) + True + + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(T(sc,sc),1.5) + True + + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(T(sc,s),1.5) + True + >>> isclose(T(s,sc),1.5) + True + + """ + if b.contracted: + return sum(cb*T(pb,a) for (cb,pb) in b) + elif a.contracted: + return sum(ca*T(b,pa) for (ca,pa) in a) + return a.norm*b.norm*kinetic(a.exponent,a.powers,a.origin, + b.exponent,b.powers,b.origin) + +def V(a,b,C): + """ + Simple interface to the nuclear attraction function. + >>> from pyquante2 import pgbf,cgbf + >>> s = pgbf(1) + >>> isclose(V(s,s,(0,0,0)),-1.595769) + True + + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(V(sc,sc,(0,0,0)),-1.595769) + True + + >>> sc = cgbf(exps=[1],coefs=[1]) + >>> isclose(V(sc,s,(0,0,0)),-1.595769) + True + + >>> isclose(V(s,sc,(0,0,0)),-1.595769) + True + + """ + if b.contracted: + return sum(cb*V(pb,a,C) for (cb,pb) in b) + elif a.contracted: + return sum(ca*V(b,pa,C) for (ca,pa) in a) + return a.norm*b.norm*nuclear_attraction(a.exponent,a.powers,a.origin, + b.exponent,b.powers,b.origin,C) + +def overlap(alpha1,lmn1,A,alpha2,lmn2,B): + """ + Full form of the overlap integral. Taken from THO eq. 2.12 + >>> isclose(overlap(1,(0,0,0),array((0,0,0),'d'),1,(0,0,0),array((0,0,0),'d')),1.968701) + True + """ + l1,m1,n1 = lmn1 + l2,m2,n2 = lmn2 + rab2 = norm2(A-B) + gamma = alpha1+alpha2 + P = gaussian_product_center(alpha1,A,alpha2,B) + + pre = pow(pi/gamma,1.5)*exp(-alpha1*alpha2*rab2/gamma) + + wx = overlap1d(l1,l2,P[0]-A[0],P[0]-B[0],gamma) + wy = overlap1d(m1,m2,P[1]-A[1],P[1]-B[1],gamma) + wz = overlap1d(n1,n2,P[2]-A[2],P[2]-B[2],gamma) + return pre*wx*wy*wz + +def overlap1d(l1,l2,PAx,PBx,gamma): + """ + The one-dimensional component of the overlap integral. Taken from THO eq. 2.12 + >>> isclose(overlap1d(0,0,0,0,1),1.0) + True + """ + total = 0 + for i in range(1+int(floor(0.5*(l1+l2)))): + total += binomial_prefactor(2*i,l1,l2,PAx,PBx)* \ + fact2(2*i-1)/pow(2*gamma,i) + return total + +def gaussian_product_center(alpha1,A,alpha2,B): + """ + The center of the Gaussian resulting from the product of two Gaussians: + >>> gaussian_product_center(1,array((0,0,0),'d'),1,array((0,0,0),'d')) + array([ 0., 0., 0.]) + """ + return (alpha1*A+alpha2*B)/(alpha1+alpha2) + +def binomial_prefactor(s,ia,ib,xpa,xpb): + """ + The integral prefactor containing the binomial coefficients from Augspurger and Dykstra. + >>> binomial_prefactor(0,0,0,0,0) + 1 + """ + total= 0 + for t in range(s+1): + if s-ia <= t <= ib: + total += binomial(ia,s-t)*binomial(ib,t)* \ + pow(xpa,ia-s+t)*pow(xpb,ib-t) + return total + +def kinetic(alpha1,lmn1,A,alpha2,lmn2,B): + """ + The full form of the kinetic energy integral + >>> isclose(kinetic(1,(0,0,0),array((0,0,0),'d'),1,(0,0,0),array((0,0,0),'d')),2.953052) + True + """ + l1,m1,n1 = lmn1 + l2,m2,n2 = lmn2 + term0 = alpha2*(2*(l2+m2+n2)+3)*\ + overlap(alpha1,(l1,m1,n1),A,\ + alpha2,(l2,m2,n2),B) + term1 = -2*pow(alpha2,2)*\ + (overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2+2,m2,n2),B) + + overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2,m2+2,n2),B) + + overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2,m2,n2+2),B)) + term2 = -0.5*(l2*(l2-1)*overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2-2,m2,n2),B) + + m2*(m2-1)*overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2,m2-2,n2),B) + + n2*(n2-1)*overlap(alpha1,(l1,m1,n1),A, + alpha2,(l2,m2,n2-2),B)) + return term0+term1+term2 + +def nuclear_attraction(alpha1,lmn1,A,alpha2,lmn2,B,C): + """ + Full form of the nuclear attraction integral + >>> isclose(nuclear_attraction(1,(0,0,0),array((0,0,0),'d'),1,(0,0,0),array((0,0,0),'d'),array((0,0,0),'d')),-3.141593) + True + """ + l1,m1,n1 = lmn1 + l2,m2,n2 = lmn2 + gamma = alpha1+alpha2 + + P = gaussian_product_center(alpha1,A,alpha2,B) + rab2 = norm2(A-B) + rcp2 = norm2(C-P) + + dPA = P-A + dPB = P-B + dPC = P-C + + Ax = A_array(l1,l2,dPA[0],dPB[0],dPC[0],gamma) + Ay = A_array(m1,m2,dPA[1],dPB[1],dPC[1],gamma) + Az = A_array(n1,n2,dPA[2],dPB[2],dPC[2],gamma) + + total = 0. + for I in range(l1+l2+1): + for J in range(m1+m2+1): + for K in range(n1+n2+1): + total += Ax[I]*Ay[J]*Az[K]*Fgamma(I+J+K,rcp2*gamma) + + val= -2*pi/gamma*exp(-alpha1*alpha2*rab2/gamma)*total + return val + +def A_term(i,r,u,l1,l2,PAx,PBx,CPx,gamma): + """ + THO eq. 2.18 + + >>> A_term(0,0,0,0,0,0,0,0,1) + 1.0 + >>> A_term(0,0,0,0,1,1,1,1,1) + 1.0 + >>> A_term(1,0,0,0,1,1,1,1,1) + -1.0 + >>> A_term(0,0,0,1,1,1,1,1,1) + 1.0 + >>> A_term(1,0,0,1,1,1,1,1,1) + -2.0 + >>> A_term(2,0,0,1,1,1,1,1,1) + 1.0 + >>> A_term(2,0,1,1,1,1,1,1,1) + -0.5 + >>> A_term(2,1,0,1,1,1,1,1,1) + 0.5 + """ + return pow(-1,i)*binomial_prefactor(i,l1,l2,PAx,PBx)*\ + pow(-1,u)*factorial(i)*pow(CPx,i-2*r-2*u)*\ + pow(0.25/gamma,r+u)/factorial(r)/factorial(u)/factorial(i-2*r-2*u) + +def A_array(l1,l2,PA,PB,CP,g): + """ + THO eq. 2.18 and 3.1 + + >>> A_array(0,0,0,0,0,1) + [1.0] + >>> A_array(0,1,1,1,1,1) + [1.0, -1.0] + >>> A_array(1,1,1,1,1,1) + [1.5, -2.5, 1.0] + """ + Imax = l1+l2+1 + A = [0]*Imax + for i in range(Imax): + for r in range(int(floor(i/2)+1)): + for u in range(int(floor((i-2*r)/2)+1)): + I = i-2*r-u + A[I] = A[I] + A_term(i,r,u,l1,l2,PA,PB,CP,g) + return A + +if __name__ == '__main__': + import doctest; doctest.testmod() diff --git a/moldesign/external/pyquante2/utils.py b/moldesign/external/pyquante2/utils.py new file mode 100644 index 0000000..bb9ad86 --- /dev/null +++ b/moldesign/external/pyquante2/utils.py @@ -0,0 +1,242 @@ +""" +utils.py - Simple utilility funtions used in pyquante2. +""" +import numpy as np +from math import factorial,lgamma +from itertools import combinations_with_replacement,combinations +from functools import reduce + +def pairs(it): return combinations_with_replacement(it,2) +def upairs(it): return combinations(it,2) + +def fact2(n): + """ + fact2(n) - n!!, double factorial of n + >>> fact2(0) + 1 + >>> fact2(1) + 1 + >>> fact2(3) + 3 + >>> fact2(8) + 384 + >>> fact2(-1) + 1 + """ + return reduce(int.__mul__,range(n,0,-2),1) + +def norm2(a): return np.dot(a,a) + +def binomial(n,k): + """ + Binomial coefficient + >>> binomial(5,2) + 10 + >>> binomial(10,5) + 252 + """ + if n == k: return 1 + assert n>k, "Attempting to call binomial(%d,%d)" % (n,k) + return factorial(n)//(factorial(k)*factorial(n-k)) + +def Fgamma(m,x): + """ + Incomplete gamma function + >>> np.isclose(Fgamma(0,0),1.0) + True + """ + SMALL=1e-12 + x = max(x,SMALL) + return 0.5*pow(x,-m-0.5)*gamm_inc(m+0.5,x) + +# def gamm_inc_scipy(a,x): +# """ +# Demonstration on how to replace the gamma calls with scipy.special functions. +# By default, pyquante only requires numpy, but this may change as scipy +# builds become more stable. +# >>> np.isclose(gamm_inc_scipy(0.5,1),1.49365) +# True +# >>> np.isclose(gamm_inc_scipy(1.5,2),0.6545103) +# True +# >>> np.isclose(gamm_inc_scipy(2.5,1e-12),0) +# True +# """ +# from scipy.special import gamma,gammainc +# return gamma(a)*gammainc(a,x) + +def gamm_inc(a,x): + """ + Incomple gamma function \gamma; computed from NumRec routine gammp. + >>> np.isclose(gamm_inc(0.5,1),1.49365) + True + >>> np.isclose(gamm_inc(1.5,2),0.6545103) + True + >>> np.isclose(gamm_inc(2.5,1e-12),0) + True + """ + assert (x > 0 and a >= 0), "Invalid arguments in routine gamm_inc: %s,%s" % (x,a) + + if x < (a+1.0): #Use the series representation + gam,gln = _gser(a,x) + else: #Use continued fractions + gamc,gln = _gcf(a,x) + gam = 1-gamc + return np.exp(gln)*gam + +def _gser(a,x): + "Series representation of Gamma. NumRec sect 6.1." + ITMAX=100 + EPS=3.e-7 + + gln=lgamma(a) + assert(x>=0),'x < 0 in gser' + if x == 0 : return 0,gln + + ap = a + delt = sum = 1./a + for i in range(ITMAX): + ap=ap+1. + delt=delt*x/ap + sum=sum+delt + if abs(delt) < abs(sum)*EPS: break + else: + print('a too large, ITMAX too small in gser') + gamser=sum*np.exp(-x+a*np.log(x)-gln) + return gamser,gln + +def _gcf(a,x): + "Continued fraction representation of Gamma. NumRec sect 6.1" + ITMAX=100 + EPS=3.e-7 + FPMIN=1.e-30 + + gln=lgamma(a) + b=x+1.-a + c=1./FPMIN + d=1./b + h=d + for i in range(1,ITMAX+1): + an=-i*(i-a) + b=b+2. + d=an*d+b + if abs(d) < FPMIN: d=FPMIN + c=b+an/c + if abs(c) < FPMIN: c=FPMIN + d=1./d + delt=d*c + h=h*delt + if abs(delt-1.) < EPS: break + else: + print('a too large, ITMAX too small in gcf') + gammcf=np.exp(-x+a*np.log(x)-gln)*h + return gammcf,gln + +def trace2(A,B): + "Return trace(AB) of matrices A and B" + return np.sum(A*B) + +def dmat(c,nclosed,nopen=0): + """Form the density matrix from the first nclosed orbitals of c. If nopen != 0, + add in half the density matrix from the next nopen orbitals. + """ + d = np.dot(c[:,:nclosed],c[:,:nclosed].T) + if nopen > 0: + d += 0.5*np.dot(c[:,nclosed:(nclosed+nopen)],c[:,nclosed:(nclosed+nopen)].T) + return d + +def symorth(S): + "Symmetric orthogonalization" + E,U = np.linalg.eigh(S) + n = len(E) + Shalf = np.identity(n,'d') + for i in range(n): + Shalf[i,i] /= np.sqrt(E[i]) + return simx(Shalf,U,True) + +def canorth(S): + "Canonical orthogonalization U/sqrt(lambda)" + E,U = np.linalg.eigh(S) + for i in range(len(E)): + U[:,i] = U[:,i] / np.sqrt(E[i]) + return U + +def cholorth(S): + "Cholesky orthogonalization" + return np.linalg.inv(np.linalg.cholesky(S)).T + +def simx(A,B,transpose=False): + "Similarity transform B^T(AB) or B(AB^T) (if transpose)" + if transpose: + return np.dot(B,np.dot(A,B.T)) + return np.dot(B.T,np.dot(A,B)) + +def ao2mo(H,C): return simx(H,C) +def mo2ao(H,C,S): return simx(H,np.dot(S,C),transpose=True) + +def geigh(H,S): + "Solve the generalized eigensystem Hc = ESc" + A = cholorth(S) + E,U = np.linalg.eigh(simx(H,A)) + return E,np.dot(A,U) + +def parseline(line,format): + """\ + Given a line (a string actually) and a short string telling + how to format it, return a list of python objects that result. + + The format string maps words (as split by line.split()) into + python code: + x -> Nothing; skip this word + s -> Return this word as a string + i -> Return this word as an int + d -> Return this word as an int + f -> Return this word as a float + + Basic parsing of strings: + >>> parseline('Hello, World','ss') + ['Hello,', 'World'] + + You can use 'x' to skip a record; you also don't have to parse + every record: + >>> parseline('1 2 3 4','xdd') + [2, 3] + + >>> parseline('C1 0.0 0.0 0.0','sfff') + ['C1', 0.0, 0.0, 0.0] + + Should this return an empty list? + >>> parseline('This line wont be parsed','xx') + """ + xlat = {'x':None,'s':str,'f':float,'d':int,'i':int} + result = [] + words = line.split() + for i in range(len(format)): + f = format[i] + trans = xlat.get(f,None) + if trans: result.append(trans(words[i])) + if len(result) == 0: return None + if len(result) == 1: return result[0] + return result + +def colorscale(mag, cmin, cmax): + """ + Return a tuple of floats between 0 and 1 for R, G, and B. + From Python Cookbook (9.11?) + """ + # Normalize to 0-1 + try: + x = float(mag-cmin)/(cmax-cmin) + except ZeroDivisionError: + x = 0.5 # cmax == cmin + blue = min((max((4*(0.75-x), 0.)), 1.)) + red = min((max((4*(x-0.25), 0.)), 1.)) + green = min((max((4*abs(x-0.5)-1., 0.)), 1.)) + return red, green, blue + +#Todo: replace with np.isclose +#def isnear(a,b,tol=1e-6): return abs(a-b) < tol + +if __name__ == '__main__': + import doctest + doctest.testmod() + diff --git a/moldesign/helpers/helpers.py b/moldesign/helpers/helpers.py index 645d427..f01e9c4 100644 --- a/moldesign/helpers/helpers.py +++ b/moldesign/helpers/helpers.py @@ -22,40 +22,6 @@ import collections -import numpy as np - -from moldesign import units as u - - -class VolumetricGrid(object): - """ - Helper object for preparing gaussian CUBE files - """ - UNITS = u.angstrom - def __init__(self, positions, padding=2.5*u.angstrom, npoints=25): - mins = positions.min(axis=0) - padding - maxes = positions.max(axis=0) + padding - self.npoints = npoints - self.xr = (mins[0], maxes[0]) - self.yr = (mins[1], maxes[1]) - self.zr = (mins[2], maxes[2]) - self._origin = mins.value_in(self.UNITS) - self.dx = (self.xr[1] - self.xr[0]).value_in(self.UNITS) / (float(npoints) - 1) - self.dy = (self.yr[1] - self.yr[0]).value_in(self.UNITS) / (float(npoints) - 1) - self.dz = (self.zr[1] - self.zr[0]).value_in(self.UNITS) / (float(npoints) - 1) - self.fxyz = None - - def xyzlist(self): - stride = self.npoints * 1j - grids = np.mgrid[self.xr[0]:self.xr[1]:stride, - self.yr[0]:self.yr[1]:stride, - self.zr[0]:self.zr[1]:stride] - return grids * self.UNITS - - def origin(self): - return tuple(self._origin) - - def get_all_atoms(*objects): """ Given Atoms, AtomContainers, lists of Atoms, and lists of AtomContainers, return a flat list of all atoms contained therein. diff --git a/moldesign/interfaces/pyscf_interface.py b/moldesign/interfaces/pyscf_interface.py index efc0258..7250da8 100644 --- a/moldesign/interfaces/pyscf_interface.py +++ b/moldesign/interfaces/pyscf_interface.py @@ -99,7 +99,7 @@ def __call__(self, info): @compute.runsremotely(enable=force_remote) def get_eris_in_basis(basis, orbs): - """ Get electron repulsion integrals transformed in (in form eri[i,j,k,l] = (ij|kl)) + """ Get electron repulsion integrals transformed into this basis (in form eri[i,j,k,l] = (ij|kl)) """ from pyscf import ao2mo @@ -121,17 +121,17 @@ def basis_values(mol, basis, coords, coeffs=None, positions=None): values are returned) Returns: - Array[length]: if ``coeffs`` is not passed, an array of basis fn values at each + Array[length**(-1.5)]: if ``coeffs`` is not passed, an array of basis fn values at each coordinate. Otherwise, a list of orbital values at each coordinate """ from pyscf.dft import numint # TODO: more than just create the basis by name ... pmol = mol_to_pyscf(mol, basis=basis.basisname, positions=positions) - aovals = numint.eval_ao(pmol, np.ascontiguousarray(coords.value_in(u.bohr))) + aovals = numint.eval_ao(pmol, np.ascontiguousarray(coords.value_in(u.bohr))) * (u.a0**-1.5) if coeffs is None: return aovals else: - return aovals.dot(coeffs) + return aovals.dot(coeffs.T) diff --git a/moldesign/mathutils/__init__.py b/moldesign/mathutils/__init__.py index 2736341..fb97fc3 100644 --- a/moldesign/mathutils/__init__.py +++ b/moldesign/mathutils/__init__.py @@ -1,2 +1,3 @@ from .vectormath import * from .eigen import * +from .grids import * \ No newline at end of file diff --git a/moldesign/mathutils/grids.py b/moldesign/mathutils/grids.py new file mode 100644 index 0000000..e49c1a8 --- /dev/null +++ b/moldesign/mathutils/grids.py @@ -0,0 +1,163 @@ +from __future__ import print_function, absolute_import, division +from future.builtins import * +from future import standard_library +standard_library.install_aliases() + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import itertools + +import numpy as np + +from .. import utils +from .. import units as u + + +@utils.exports +class VolumetricGrid(object): + """ Creates a 3D, rectangular grid of points + + Args: + xrange (Tuple[len=2]): (min,max) in x direction + yrange (Tuple[len=2]): (min,max) in y direction + zrange (Tuple[len=2]): (min,max) in z direction + xpoints (int): number of grid lines in x direction (default: npoints) + ypoints (int): number of grid lines in y direction (default: npoints) + zpoints (int): number of grid lines in z direction (default: npoints) + npoints (int): synonym for "xpoints" + """ + dx, dy, dz = (utils.IndexView('deltas', i) for i in range(3)) + xr, yr, zr = (utils.IndexView('ranges', i) for i in range(3)) + xpoints, ypoints, zpoints = (utils.IndexView('points', i) for i in range(3)) + xspace, yspace, zspace = (utils.IndexView('spaces', i) for i in range(3)) + + def __init__(self, xrange, yrange, zrange, + xpoints=None, ypoints=None, zpoints=None, + npoints=32): + + if xpoints is not None: + npoints = xpoints + + self.points = np.array([xpoints if xpoints is not None else npoints, + ypoints if ypoints is not None else npoints, + zpoints if zpoints is not None else npoints], + dtype='int') + + self.ranges = u.array([xrange, yrange, zrange]) + self.deltas = (self.ranges[:,1] - self.ranges[:,0]) / (self.points - 1.0) + self.spaces = [u.linspace(*r, num=p) for r,p in zip(self.ranges, self.points)] + + @property + def ndims(self): + return len(self.ranges) + ndim = num_dims = ndims + + @property + def origin(self): + """ Vector[len=3]: the origin of the grid (the lowermost corner in each dimension) + """ + origin = [r[0] for r in (self.ranges)] + try: + return u.array(origin) + except u.DimensionalityError: + return origin + + @property + def npoints(self): + """ int: total number of grid points in this grid + """ + return np.product(self.points) + + def iter_points(self): + """ Iterate through every point on the grid. + + Always iterates in the same order. The x-index is the most slowly varying, + the z-index is the fastest. + + Yields: + Vector[len=3]: x,y,z coordinate of each point on the grid + """ + for i,j,k in itertools.product(*map(range, self.points)): + yield self.origin + self.deltas * [i,j,k] + + def allpoints(self): + """ Return an array of all coordinates on the grid. + + This obviously takes a lot of memory, but is useful for evaluating vectorized functions + on this grid. Points are returned in the same order as ``iter_points``. + + Yields: + Matrix[shape=(self.npoints**3,3)]: x,y,z coordinate of each point on the grid + """ + return _cartesian_product(self.spaces) + + +def _cartesian_product(arrays): + """ Fast grid creation routine from @senderle on Stack Overflow. This is awesome. + + Modifications: + 1. Module import names + 2. Unit handling + + References: + https://stackoverflow.com/a/11146645/1958900 + """ + import functools + + # AMV mod for units handling + units = getattr(arrays[0], 'units', u.ureg.dimensionless) + + # original code + broadcastable = np.ix_(*arrays) + broadcasted = np.broadcast_arrays(*broadcastable) + rows, cols = functools.reduce(np.multiply, broadcasted[0].shape), len(broadcasted) + out = np.empty(rows * cols, dtype=broadcasted[0].dtype) + start, end = 0, rows + for a in broadcasted: + out[start:end] = a.reshape(-1) + start, end = end, end + rows + + result = out.reshape(cols, rows).T + + # AMV mod for units handling + if units == u.ureg.dimensionless: + return result + else: + # do it this way to avoid copying a huge array + quantity = u.MdtQuantity(result, units=units) + return quantity + + +@utils.exports +def padded_grid(positions, padding, npoints=25): + """ Creates a 3D, rectangular grid of points surrounding a set of positions + The points in the grid are evenly spaced in each dimension, but spacing may differ between + in different dimensions. + + Args: + positions (Matrix[shape=(*,3)]): positions to create the grid around + padding (Scalar): how far to extend the grid past the positions in each dimension + npoints (int): number of points in each direction (total number of points is npoints**3) + + Returns: + VolumetricGrid: grid object + """ + mins = positions.min(axis=0)-padding + maxes = positions.max(axis=0)+padding + + xr = (mins[0], maxes[0]) + yr = (mins[1], maxes[1]) + zr = (mins[2], maxes[2]) + + return VolumetricGrid(xr, yr, zr, npoints) diff --git a/moldesign/mathutils/spherical_harmonics.py b/moldesign/mathutils/spherical_harmonics.py new file mode 100644 index 0000000..96cb60e --- /dev/null +++ b/moldesign/mathutils/spherical_harmonics.py @@ -0,0 +1,199 @@ +from __future__ import print_function, absolute_import, division +from future import standard_library +standard_library.install_aliases() +from future.builtins import * + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import numpy as np + +__all__ = ['SPHERE_TO_CART'] + + +def cart_to_polar_angles(coords): + if len(coords.shape) == 2: + r_xy_2 = np.sum(coords[:,:2]**2, axis=1) + theta = np.arctan2(np.sqrt(r_xy_2), coords[:,2]) + phi = np.arctan2(coords[:,1], coords[:,0]) + return theta, phi + else: + assert len(coords) == 3 and len(coords.shape) == 1 + r_xy_2 = np.sum(coords[:2]**2) + theta = np.arctan2(np.sqrt(r_xy_2), coords[2]) + phi = np.arctan2(coords[1], coords[0]) + return theta, phi + + +class Y(object): + r""" A real-valued spherical harmonic function + + These functions are orthonormalized over the sphere such that + + .. math:: + \int d^2\Omega \: Y_{lm}(\theta, \phi) Y_{l'm'}(\theta, \phi) = \delta_{l,l'} \delta_{m,m'} + + References: + https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics + """ + def __init__(self, l, m): + self.l = l + self.m = m + + self._posneg = -1 + if self.m < 0: + self._posneg *= -1 + if self.m % 2 == 0: + self._posneg *= -1 + + def __call__(self, coords): + from scipy.special import sph_harm + + theta, phi = cart_to_polar_angles(coords) + if self.m == 0: + return (sph_harm(self.m, self.l, phi, theta)).real + + vplus = sph_harm(abs(self.m), self.l, phi, theta) + vminus = sph_harm(-abs(self.m), self.l, phi, theta) + value = np.sqrt(1/2.0) * (self._posneg*vplus + vminus) + + if self.m < 0: + return -value.imag + else: + return value.real + + +class Cart(object): + def __init__(self, px, py, pz, coeff): + self.powers = np.array([px, py, pz], dtype='int') + self.coeff = coeff + self._l = self.powers.sum() + + def __call__(self, coords): + """ Evaluate this function at the given list of coordinates + + Args: + coords (Vector[len=3] or Matrix[matrix=shape(*,3)): coordinate(s) at which to evaluate + + Returns: + Scalar or Vector: value of the function at these coordinates + """ + + c = coords**self.powers + + if len(coords.shape) == 2: + r_l = (coords**2).sum(axis=1) ** (-self._l / 2.0) + return self.coeff * np.product(c, axis=1) * r_l + else: + r_l = (coords**2).sum() ** (-self._l / 2.0) + return self.coeff * np.product(c) * r_l + + def __iter__(self): + yield self # for interface compatibility with the CartSum class + + +class CartSum(object): + def __init__(self, coeff, carts): + self.coeff = coeff + + prefacs = [] + powers = [] + self.l = None + for cart in carts: + powers.append(np.array(cart[:3], dtype='int')) + prefacs.append(np.array(cart[3])) + if self.l is None: + self.l = sum(cart[:3]) + else: + assert self.l == sum(cart[:3]) + + self.prefactors = np.array(prefacs) + self.powers = np.array(powers) + + def __call__(self, coords): + # likely can speed this up a lot using clever numpy broadcasting + + if len(coords.shape) == 2: + c = np.zeros((len(coords), )) + axis = 1 + r_l = (coords**2).sum(axis=1) ** (-self.l / 2.0) + else: + c = 0.0 + axis = None + r_l = (coords**2).sum() ** (-self.l / 2.0) + + for factor, power in zip(self.prefactors, self.powers): + c += factor * np.product(coords**power, axis=axis) + + return c * self.coeff * r_l + + def __iter__(self): + for pf, (px, py, pz) in zip(self.prefactors, self.powers): + yield Cart(px, py, pz, coeff=self.coeff*pf) + + + +def sqrt_x_over_pi(num, denom): + return np.sqrt(num / (denom*np.pi)) + + + ############ s ############ +SPHERE_TO_CART = {(0, 0): Cart(0, 0, 0, sqrt_x_over_pi(1, 4)), + + ############ p ############ + (1, -1): Cart(0, 1, 0, sqrt_x_over_pi(3, 4)), + (1, 0): Cart(0, 0, 1, sqrt_x_over_pi(3, 4)), + (1, 1): Cart(1, 0, 0, sqrt_x_over_pi(3, 4)), + + ############ d ############ + (2, -2): Cart(1, 1, 0, sqrt_x_over_pi(15, 4)), + (2, -1): Cart(0, 1, 1, sqrt_x_over_pi(15, 4)), + (2, 0): CartSum(sqrt_x_over_pi(5, 16), + [(2, 0, 0, -1.0), + (0, 2, 0, -1.0), + (0, 0, 2, 2.0)]), + (2, 1): Cart(1, 0, 1, sqrt_x_over_pi(15, 4)), + (2, 2): CartSum(sqrt_x_over_pi(15, 16), + [(2, 0, 0, 1.0), + (0, 2, 0, -1.0)]), + + ############ f ############ + (3, -3): CartSum(sqrt_x_over_pi(35, 32), + [(2, 1, 0, 3.0), + (0, 3, 0, -1.0)]), + + (3, -2): Cart(1, 1, 1, sqrt_x_over_pi(105, 4)), + + (3, -1): CartSum(sqrt_x_over_pi(21, 32), + [(0, 1, 2, 4.0), + (2, 1, 0, -1.0), + (0, 3, 0, -1.0)]), + + (3, 0): CartSum(sqrt_x_over_pi(7, 16), + [(0, 0, 3, 2.0), + (2, 0, 1, -3.0), + (0, 2, 1, -3.0)]), + + (3, 1): CartSum(sqrt_x_over_pi(21, 32), + [(1, 0, 2, 4.0), + (3, 0, 0, -1.0), + (1, 2, 0, -1.0)]), + + (3, 2): CartSum(sqrt_x_over_pi(105, 16), + [(2, 0, 1, 1.0), + (0, 2, 1, -1.0)]), + + (3, 3): CartSum(sqrt_x_over_pi(35, 32), + [(3, 0, 0, 1.0), + (1, 2, 0, -3.0)]), + } \ No newline at end of file diff --git a/moldesign/min/smart.py b/moldesign/min/smart.py index ba88cd6..44b0414 100644 --- a/moldesign/min/smart.py +++ b/moldesign/min/smart.py @@ -60,8 +60,8 @@ def _run(self): forces = self.mol.calculate_forces() if abs(forces).max() <= self.gd_threshold: self._spmin = self._make_quadratic_method() + self._spmin.traj = self.traj self._spmin._run() - self.traj = self._spmin.traj return # Otherwise, remove large forces with gradient descent; exit if we pass the cycle limit diff --git a/moldesign/models/pyscf.py b/moldesign/models/pyscf.py index 60f3cda..08c1113 100644 --- a/moldesign/models/pyscf.py +++ b/moldesign/models/pyscf.py @@ -311,7 +311,7 @@ def _converge(self, method, dm0=None): # fallback 2: level shift, slower convergence # this probably won't converge, but is intended to get us in the right basin for the next step - # NEWFEATURE: should dynamically adjust level shift instead of hardcoded cycles + # TODO: should dynamically adjust level shift instead of hardcoded cycles self.logger.handled('SCF failed to converge. Performing %d iterations with level shift of -0.5 hartree' % (old_div(method.max_cycle, 2))) failed.append(method) @@ -440,7 +440,7 @@ def _get_ao_basis_functions(self): atom = self.mol.atoms[pmol.bas_atom(ishell)] angular = pmol.bas_angular(ishell) num_momentum_states = angular*2 + 1 - exps = pmol.bas_exp(ishell) + exps = pmol.bas_exp(ishell) / (u.bohr**2) num_contractions = pmol.bas_nctr(ishell) coeffs = pmol.bas_ctr_coeff(ishell) @@ -450,12 +450,12 @@ def _get_ao_basis_functions(self): sphere_label = label[3] l, m = SPHERICAL_NAMES[sphere_label] assert l == angular - # TODO: This is not really the principal quantum number + # Note: this is for metadata only, should not be used in any calculations n = int(''.join(x for x in label[2] if x.isdigit())) - primitives = [orbitals.SphericalGaussian(atom.position.copy(), - exp, n, l, m, - coeff=coeff[ictr]) + exp, l, m, + coeff=coeff[ictr], + normalized=True) for exp, coeff in zip(exps, coeffs)] bfs.append(orbitals.AtomicBasisFunction(atom, n=n, l=angular, m=m, diff --git a/moldesign/molecules/trajectory.py b/moldesign/molecules/trajectory.py index 1a1fbdc..f853354 100644 --- a/moldesign/molecules/trajectory.py +++ b/moldesign/molecules/trajectory.py @@ -370,7 +370,7 @@ def _new_property(self, key, value): else: try: proplist = self.unit_system.convert(u.array([value])) - except (TypeError, u.UndefinedUnitError): + except (ValueError, TypeError, u.UndefinedUnitError): proplist = [value] else: proplist.make_resizable() diff --git a/moldesign/orbitals/__init__.py b/moldesign/orbitals/__init__.py index 43f0b03..4d3ae95 100644 --- a/moldesign/orbitals/__init__.py +++ b/moldesign/orbitals/__init__.py @@ -3,7 +3,11 @@ def toplevel(o): return o __all__ = [] +from .primitives import * from .gaussians import * from .orbitals import * +from .cartesian import * +from .spherical import * +from .atomic_basis_fn import * from .basis import * from .wfn import * diff --git a/moldesign/orbitals/atomic_basis_fn.py b/moldesign/orbitals/atomic_basis_fn.py new file mode 100644 index 0000000..a182d16 --- /dev/null +++ b/moldesign/orbitals/atomic_basis_fn.py @@ -0,0 +1,120 @@ +from __future__ import print_function, absolute_import, division +from future.builtins import * +from future import standard_library +standard_library.install_aliases() + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from . import PrimitiveSum, SHELLS, SPHERICALNAMES + + +class AtomicBasisFunction(PrimitiveSum): + """ Stores an atomic basis function. + + Note: + Either l and m should be passed, or cart, but not both. + + Args: + atom (moldesign.Atom): The atom this basis function belongs to + index (int): the index of this basis function (it is stored as + ``wfn.basis[self.index]``) + n (int): principal quantum number (``n>=1``) - shell metadata (used for labeling only) + l (int): total angular momentum quantum number (``l<=n-1``) + m (int): z-angular momentum quantum number (optional - + for spherical sets only; ``|m|<=l``) + cart (str): cartesian component (optional; for cartesian sets only) + primitives (List[PrimitiveBase]): List of primitives, if available + """ + def __init__(self, atom, n=None, l=None, m=None, cart=None, primitives=None): + super().__init__(primitives) + self._atom = atom + self.atom_index = atom.index + self.n = n + self.l = l + self.m = m + if cart is not None: + assert self.m is None, 'Both cartesian and spherical components passed!' + assert len(cart) == self.l, \ + 'Angular momentum does not match specified component %s' % cart + for e in cart: + assert e in 'xyz' + self.cart = ''.join(sorted(cart)) + else: + self.cart = None + + # These quantities can't be defined until we assemble the entire basis + self.coeffs = None + self.molecule = atom.molecule + self.basis = None + self.occupation = None + self.wfn = None + + @property + def atom(self): + """ moldesign.Atom: the atom this basis function belongs to + + We get the atom via an indirect reference, making it easier to copy the wavefunction + """ + if self.wfn is not None: + return self.wfn.mol.atoms[self.atom_index] + else: + return self._atom + + @property + def orbtype(self): + """ A string describing the orbital's angular momentum state. + + Examples: + >>> AtomicBasisFunction(n=1, l=0).orbtype + 's' + >>> AtomicBasisFunction(n=2, l=1, cart='y').orbtype + 'py' + >>> AtomicBasisFunction(n=3, l=2, m=0).orbtype + 'd(z^2)' + """ + if self.l == 0: t = 's' + elif self.cart is not None: t = SHELLS[self.l] + self.cart + else: t = SPHERICALNAMES[self.l, self.m] + return t + + @property + def aotype(self): + """ A string describing the orbital's state. + + Examples: + >>> AtomicBasisFunction(n=1, l=0).aotype + '1s' + >>> AtomicBasisFunction(n=2, l=1, cart='y').aotype + '2py' + >>> AtomicBasisFunction(n=3, l=2, m=0).aotype + '3d(z^2)' + """ + t = self.orbtype + if self.n: + return '%s%s' % (self.n, t) + else: + return t + + def __str__(self): + return 'AO ' + self.name + + @property + def name(self): + try: + return '%s on atom %s' % (self.aotype, self.atom.name) + except: + return 'Basis Fn' + + def __repr__(self): + return '<%s %s>' % (self.__class__.__name__, self.name) \ No newline at end of file diff --git a/moldesign/orbitals/basis.py b/moldesign/orbitals/basis.py index a65913a..68d4951 100644 --- a/moldesign/orbitals/basis.py +++ b/moldesign/orbitals/basis.py @@ -23,21 +23,6 @@ from . import toplevel, MolecularOrbitals -class PrimitiveBase(object): - def __init__(self, *args, **kwargs): - raise NotImplementedError() - - def __call__(self, position): - raise NotImplementedError() - - def overlap(self, other): - raise NotImplementedError() - - @property - def norm(self, other): - raise NotImplementedError() - - @toplevel class BasisSet(MolecularOrbitals): """ @@ -85,26 +70,41 @@ def __init__(self, mol, orbitals, name=None, self._basis_fn_by_atom.setdefault(fn.atom.index, []).append(fn) def __call__(self, coords, coeffs=None): - """Calculate the value of each basis function at the specified coordinates. + """Calculate the amplitude of a list of orbitals at the specified coordinates. - Note: - This is just a pass-through to a specific implementation - PYSCF's eval_ao routine - for now. + Returns an array of orbital amplitudes. The amplitude of the _n_th orbital at the + _j_th coordinate is stored at position ``(j, n)`` in the returned array. Args: - coords (Array[length]): List of coordinates. - coeffs (Vector): List of ao coefficients (optional; if not passed, all basis fn - values are returned) + coords (Matrix[shape=(*,3)]): List of coordinates. + coeffs (Matrix[shape=(*, nbasis)]): List of ao coefficients (optional; + if not passed, the amplitudes of all basis functions will be returned Returns: - Array[length]: if ``coeffs`` is not passed, an array of basis fn values at each - coordinate. Otherwise, a list of orbital values at each coordinate + Matrix: + - if ``coeffs`` is passed, an array of orbital amplitudes at the given positions + of size ``(len(coords), len(coeffs))``. + - if ``coeffs`` is NOT passed, an array of basis function amplitudes + of size ``(len(coords), len(aobasis))``. """ - from moldesign.interfaces.pyscf_interface import basis_values - return basis_values(self.wfn.mol, self, coords, coeffs=coeffs, - positions=self.wfn.positions) + basis_vals = np.zeros((len(coords), len(self))) * self.orbitals[0]._get_wfn_units() + for ibf, bf in enumerate(self.orbitals): + basis_vals[:, ibf] = bf(coords) + + if coeffs is None: + return basis_vals + else: + return np.dot(basis_vals, coeffs.T) def get_basis_functions_on_atom(self, atom): + """ Return a list of basis functions on this atom + + Args: + atom (moldesign.Atom): query atom + + Returns: + List[AtomicBasisFunction]: basis functions centered on this atom + """ return self._basis_fn_by_atom[atom.index] def __repr__(self): diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py new file mode 100644 index 0000000..d3bc709 --- /dev/null +++ b/moldesign/orbitals/cartesian.py @@ -0,0 +1,281 @@ +from __future__ import print_function, absolute_import, division +from future.builtins import * +from future import standard_library +standard_library.install_aliases() + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import itertools + +import numpy as np +from scipy.special import binom + +from .. import units as u +from .. import utils +from . import Primitive, PrimitiveSum, Gaussian + + +class CartesianGaussian(Primitive): + r""" Stores an N-dimensional gaussian function of the form: + + .. math:: + G(\mathbf r) = C \times \left( \prod_{i=1}^N{{r_i}^{p_i} } \right) + e^{-\alpha |\mathbf r - \mathbf{r}_0|^2} + + For a three-dimensional gaussian, this is + + ..math:: + G(x,y,z) = C \times x^{p_1} y^{p_2} z^{p_3} e^{-\alpha |\mathbf r - \mathbf{r}_0|^2} + + where *C* is ``self.coeff``, :math:`\alpha` is ``self.alpha``, *r0* is ``self.center``, and + :math:`p_1, p_2, ...` are given in the array ``self.powers`` + + This function is evaluated as the the product of three terms that can be accessed separately: + - ``self.coeff`` - the constant prefactor _C_ + - ``self.radial_part(r)`` - :math:`e^{-\alpha \left| \mathbf r - \mathbf r_0 \right|^2}` + - ``self.angular_part(r)`` - :math:`\prod_{i=1}^N{{r_i}^{p_i} }` + + Args: + center (Vector[length]): location of the gaussian's centroid + powers (List[int]): cartesian powers in each dimension (see + equations in :class:`CartesianGaussian` docs) + alpha (Scalar[1/length**2]): gaussian width parameter + coeff (Scalar): multiplicative coefficient (if None, gaussian will be automatically + normalized) + + Note: + The dimensionality of the gaussian is determined by the dimensionality + of the centroid location vector and the power vector. So, if scalars are passed for the + ``center`` and ``powers``, it's 1-D. If length-3 vectors are passed for ``center`` + and ``powers``, it's 3D. + + References: + Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. + """ + center = utils.Alias('radial_part.center') + alpha = utils.Alias('radial_part.alpha') + ndims = ndim = num_dims = utils.Alias('radial_part.ndim') + + def __init__(self, center, alpha, powers, coeff=None, normalized=True): + assert len(powers) == len(center), "Inconsistent dimensionality - number of cartesian " \ + "powers must match dimensionality of centroid vector" + + self.powers = np.array(powers) + self.shell = sum(self.powers) + self.radial_part = Gaussian(center, alpha, coeff=1.0, normalized=False) + super().__init__(coeff=coeff, normalized=normalized) + + def __repr__(self): + return ("<{ndim}-D cartesian gaussian (norm: {norm:4.2f}, " + "cartesian powers: {powers}, " + "alpha: {exp:4.2f}, " + "center: {center}>").format( + ndim=self.ndim, + center=self.center, exp=self.alpha, + powers=tuple(self.powers), norm=self.norm) + + def _get_wfn_units(self): + return u.MdtQuantity(self.coeff * u.default.length**self.angular).units + + @property + def angular(self): + """ Angular momentum of this function (sum of cartesian powers) + """ + return self.powers.sum() + + def overlap(self, other, normalized=False): + r""" Overlap of this function with another basis function + + .. math:: + \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r + + Args: + other (Primitive): + normalized (bool): If True, return the overlap of the two NORMALIZED functions. + + Returns: + Scalar: value of the overlap + """ + olap = 0.0 + for other_prim in other.iterprimitives(): + if hasattr(other_prim, 'to_cart'): + other_prim = other_prim.to_cart() + + for p in other_prim.iterprimitives(): + olap += self._overlap_cart_cart(p) + + if normalized: + olap = olap / (self.norm * other.norm) + + return olap + + def _overlap_cart_cart(self, other): + """ Overlap of two cartesian functions. Uses the PyQuante2 python implementation + + Users shouldn't need to call directly, this should be automatically called when necessary + """ + from ..external.pyquante2.one import overlap + lengthunits = u.default.get_default(self.center) + alphaunits = (1/lengthunits**2).units + + alphas = [alphaunits.value_of(g.alpha) for g in (self, other)] + centers = [lengthunits.value_of(g.center) for g in (self, other)] + + olap_base = overlap(alphas[0], self.powers, centers[0], + alphas[1], other.powers, centers[1]) + + olap = self.coeff * other.coeff * olap_base + + if lengthunits != u.dimensionless: + olap *= lengthunits ** (3 + sum(self.powers) + sum(other.powers)) + + return olap + + def __call__(self, coords): + """ Evaluate this function at the given coordinates. + + Can be called either with a 1D column (e.g., ``[1,2,3]*u.angstrom ``) or + an ARRAY of coordinates (``[[0,0,0],[1,1,1]] * u.angstrom``) + + Args: + coords (Vector[length, len=3] or Matrix[length, shape=(*,3)]): Coordinates or + list of 3D coordinates + + Examples: + >>> g = CartesianGaussian([0,0,0]*u.angstrom, alpha=1.0/u.angstrom**2, powers=(0,0,0)) + >>> # Value at a single coordinate: + >>> g([0,0,0] * u.angstrom) + 1.0 + >>> # Values at a list of coordinates + >>> g[[0,0,0], [0,0,1], [0.5,0.5,0.5] * u.angstrom] + array([ 1.0, 0.36787944, 0.47236655]) + + Returns: + Scalar or Vector: function value(s) at the passed coordinates + """ + result = self.radial_part(coords) + + if self.shell > 0: + result *= self.angular_part(coords) + return result * self.coeff + + def angular_part(self, coords): + if self.shell == 0: + return 1.0 + + if len(coords.shape) > 1: + axis = 1 + else: + axis = None + + disp = coords - self.center + if hasattr(disp, 'units'): + mag = disp.magnitude + units = disp.units ** self.powers.sum() + else: + mag = disp + units = 1.0 + return np.product(mag**self.powers, axis=axis) * units + + def __mul__(self, other): + """ Returns product of two cartiesian gaussian functions as a list of product gaussians + + The returned gaussians all have the same center and width, and differ only in their + angular parts. + + Args: + other (CartesianGaussian): other gaussian basis function + + Returns: + CartesianGaussian or PrimitiveSum[CartesianGaussian]: product functions + """ + + if (self.center == other.center).all(): # this case is much easier than the general one + cpy = self.copy() + cpy.radial_part.alpha = self.alpha + other.alpha + cpy.powers = self.powers + other.powers + cpy.coeff = self.coeff * other.coeff + return cpy + + newcenter = self.radial_part * other.radial_part + + partial_coeffs = [{} for idim in range(self.ndim)] + for idim in range(self.ndim): + r_self = self.center[idim] + r_other = other.center[idim] + r_new = newcenter.center[idim] + + for m in range(self.powers[idim]+1): + for k in range(other.powers[idim]+1): + powercoeff = (binom(self.powers[idim], m) * binom(other.powers[idim], k) * + ((r_new - r_self) ** (self.powers[idim]-m)) * + ((r_new - r_other) ** (other.powers[idim]-k))) + if powercoeff == 0.0: + continue + newpower = m+k + if newpower not in partial_coeffs[idim]: + partial_coeffs[idim][newpower] = powercoeff + else: + partial_coeffs[idim][newpower] += powercoeff + + new_gaussians = [] + for powers in itertools.product(*[x.keys() for x in partial_coeffs]): + final_coeff = self.coeff * other.coeff * newcenter.coeff + for idim, p in enumerate(powers): + final_coeff *= partial_coeffs[idim][p] + new_gaussians.append(CartesianGaussian(newcenter.center, newcenter.alpha, + powers=powers, coeff=final_coeff, + normalized=False)) + + if len(new_gaussians) == 0: # no non-zero components, just return a zeroed gaussian + return CartesianGaussian(newcenter, newcenter.alpha, + self.powers + other.powers, coeff=0.0, + normalized=False) + elif len(new_gaussians) == 1: + return new_gaussians[0] + else: + return PrimitiveSum(new_gaussians) + + @property + def integral(self): + r"""Integral of this this gaussian over all N-dimensional space. + + This is implemented only for 0 and positive integer cartesian powers. + The integral is 0 if any of the powers are odd. Otherwise, the integral + is given by: + .. math:: + \int G(\mathbf r) d^N \mathbf r & = c \int d^N e^{-a x^2} \mathbf r + \prod_{i=1}^N{{r_i}^{p_i} } \\ + &= (2a)^{-\sum_i p_i} \left( \frac{\pi}{2 a} \right) ^ {N/2} \prod_{i=1}^N{(p_i-1)!!} + + where *N* is the dimensionality of the gaussian, :math:`p_i` are the cartesian powers, + and _!!_ is the "odd factorial" (:math:`n!!=1\times 3\times 5 \times ... \times n`) + + References: + Dwight, Herbert B. Tables of Integrals and other Mathematical Data, 3rd ed. + Macmillan 1957. 201. + """ + from ..data import ODD_FACTORIAL + + integ = self.coeff * self.radial_part.integral + for p in self.powers: + if p == 0: # no contribution + continue + elif p % 2 == 1: # integral of odd function is exactly 0 + return 0.0 + elif p < 0: + raise ValueError('Powers must be positive or 0') + else: + integ *= ODD_FACTORIAL[p-1]/((2.0 * self.alpha) ** (p//2)) + return integ diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 4971b3f..fe38be0 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -16,142 +16,52 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -from past.utils import old_div import numpy as np -from .orbitals import Orbital, SHELLS, SPHERICALNAMES +from .. import units as u +from . import Primitive -class AbstractFunction(object): - """ Abstract base class for basis functions - """ - def normalize(self): - """ Give this function unit norm by adjusting its coefficient - """ - self.coeff /= np.sqrt(self.norm) - - def overlap(self, other, normalized=False): - r""" Overlap of this function with another: - - .. math:: - \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r - - Args: - other (AbstractFunction): - normalized (bool): If True, return the overlap of the two NORMALIZED functions. - - Returns: - Scalar: value of the overlap - """ - newfn = self * other - integral = newfn.integral - if normalized: - integral /= np.sqrt(self.norm*other.norm) - return integral - - @property - def norm(self): - r""" The L2-Norm of this gaussian: - - .. math:: - \int \left| G(\mathbf r) \right|^2 d^N \mathbf r - """ - return self.overlap(self) - - def __str__(self): - return "%d-D %s with norm %s"%(self.ndims, self.__class__, self.norm) - - -class CartesianGaussian(AbstractFunction): - r""" Stores an N-dimensional gaussian function of the form: +class Gaussian(Primitive): + r""" N-dimensional gaussian function. + The function is given by: .. math:: - G(\mathbf r) = C \times \left( \prod_{i=1}^N{{r_i}^{p_i} } \right) - e^{-a |\mathbf r - \mathbf{r}_0|^2} - - For a three-dimensional gaussian, this is - - ..math:: - G(x,y,z) = C \times x^{p_1} y^{p_2} z^{p_3} e^{-a |\mathbf r - \mathbf{r}_0|^2} + G(\mathbf r) = C e^{-\alpha\left| \mathbf r - \mathbf r_0 \right|^2} - where *C* is ``self.coeff``, *a* is ``self.exp``, *r0* is ``self.center``, and - :math:`p_1, p_2, ...` are given in the array ``self.powers`` + where *C* is ``self.coeff``, *a* is ``self.alpha``, and :math:`\mathbf r_0` is ``self.center``. - References: - Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. + Args: + center (Vector): centroid of this function (r_0 in the equation above) + alpha (Scalar): constant in the exponential (alpha in the equation above) + coeff (Scalar): coefficient for this function. If normalized=True, the coefficient + multiplies the _normalized_ gaussian. If not passed, the gaussian will automatically + be normalized + normalized (bool): if True, normalize the function before applying the coefficient """ - def __init__(self, center, exp, powers, coeff=None): - """ Initialization: - - Args: - center (Vector[length]): location of the gaussian's centroid - powers (List[int]): cartesian powers in each dimension (see - equations in :class:`CartesianGaussian` docs) - exp (Scalar[1/length**2]): gaussian width parameter - coeff (Scalar): multiplicative coefficient (if None, gaussian will be automatically - normalized) - - Note: - The dimensionality of the gaussian is determined by the dimensionality - of the centroid location vector and the power vector. So, if scalars are passed for the - ``center`` and ``powers``, it's 1-D. If length-3 vectors are passed for ``center`` - and ``powers``, it's 3D. - """ - assert len(powers) == len(center), "Inconsistent dimensionality - number of cartesian " \ - "powers must match dimensionality of centroid vector" - self.center = center - self.exp = exp - self.powers = np.array(powers) - if coeff is None: - self.coeff = 1.0 # dummy value overwritten by self.normalize() - self.normalize() - else: - self.coeff = coeff - - self._cartesian = (self.powers != 0).any() - - def __repr__(self): - return ("<{ndim}-D Gaussian (cart) (coeff: {coeff:4.2f}, " - "cartesian powers: {powers}, " - "exponent: {exp:4.2f}, " - "center: {center}>").format( - ndim=self.ndim, - center=self.center, exp=self.exp, - powers=tuple(self.powers), coeff=self.coeff) + def __init__(self, center, alpha, coeff=None, normalized=True): + self.center = u.default.convert(center) + self.alpha = u.default.convert(alpha) + super().__init__(coeff=coeff, normalized=normalized) - @property - def ndim(self): - return len(self.powers) - num_dimensions = ndims = ndim - - @property - def angular(self): - """ Angular momentum of this function (sum of cartesian powers) - """ - return self.powers.sum() + def _get_wfn_units(self): + return u.MdtQuantity(self.coeff).units - def __call__(self, coords, _include_cartesian=True): + def __call__(self, coords, _getvals=False): """ Evaluate this function at the given coordinates. Can be called either with a 1D column (e.g., ``[1,2,3]*u.angstrom ``) or an ARRAY of coordinates (``[[0,0,0],[1,1,1]] * u.angstrom``) Args: - _include_cartesian (bool): include the contribution from the cartesian components - (for computational efficiency, this can sometimes omited now and included later) - - Example: - >>> g = CartesianGaussian([0,0,0]*u.angstrom, exp=1.0/u.angstrom**2, powers=(0,0,0)) - >>> g([0,0,0] * u.angstrom) - 1.0 - >>> g[[0,0,0], [0,0,1], [0.5,0.5,0.5] * u.angstrom] - array([ 1.0, 0.36787944, 0.47236655]) + _include_angular (bool): include the contribution from the non-exponential parts + (for computational efficiency, this can be omitted now and included later) Args: - coords (Vector[length]): Coordinates or list of coordinates + coords (Vector[length]): 3D Coordinates or list of 3D coordinates Returns: - Scalar: function value(s) at the passed coordinates + Scalar or Vector: function value(s) at the passed coordinates """ if len(coords.shape) > 1: axis = 1 @@ -162,266 +72,51 @@ def __call__(self, coords, _include_cartesian=True): prd = disp*disp # don't use np.dot - allow multiple coords at once r2 = prd.sum(axis=axis) - result = self.coeff * np.exp(-self.exp * r2) - if self._cartesian and _include_cartesian: - result *= (np.product(disp.magnitude**self.powers, axis=axis) - * disp.units**self.powers.sum() ) - return result - - def __mul__(self, other): - """ Returns product of two gaussian functions, which is also a gaussian - - Args: - other (CartesianGaussian): other gaussian wavepacket - - Returns: - CartesianGaussian: product gaussian - """ - - # convert widths to prefactor form - a = self.exp - b = other.exp - exp = a + b - center = old_div((a*self.center + b*other.center),(a+b)) - powers = self.powers + other.powers - return CartesianGaussian(center=center, exp=exp, - powers=powers, coeff=self.coeff*other.coeff) - - @property - def integral(self): - r"""Integral of this this gaussian over all N-dimensional space. - - This is implemented only for 0 and positive integer cartesian powers. - The integral is 0 if any of the powers are odd. Otherwise, the integral - is given by: - .. math:: - \int G(\mathbf r) d^N \mathbf r & = c \int d^N e^{-a x^2} \mathbf r - \prod_{i=1}^N{{r_i}^{p_i} } \\ - &= (2a)^{-\sum_i p_i} \left( \frac{\pi}{2 a} \right) ^ {N/2} \prod_{i=1}^N{(p_i-1)!!} - - where *N* is the dimensionality of the gaussian, :math:`p_i` are the cartesian powers, - and _!!_ is the "odd factorial" (:math:`n!!=1\times 3\times 5 \times ... \times n`) - - References: - Dwight, Herbert B. Tables of Integrals and other Mathematical Data, 3rd ed. - Macmillan 1957. 201. - """ - integ = (old_div(np.pi,self.exp))**(old_div(self.ndim,2.0)) - for p in self.powers: - if p == 0: # no contribution - continue - elif p % 2 == 1: # integral of odd function is exactly 0 - return 0.0 - elif p < 0: - raise ValueError('Powers must be positive or 0') - else: - integ *= old_div(_ODD_FACTORIAL[p-1],(2 ** (p+1))) - return self.coeff * integ - - -def Gaussian(center, exp, coeff=1.0): - r""" Constructor for an N-dimensional gaussian function. - - The function is given by: - .. math:: - G(\mathbf r) = C e^{-a\left| \mathbf r - \mathbf r_0 \right|^2} - - where *C* is ``self.coeff``, *a* is ``self.exp``, and :math:`\mathbf r_0` is ``self.center``. - - Note: - This is just a special case of a cartesian gaussian where all the powers are 0. - """ - return CartesianGaussian(center, exp, - powers=[0 for x in center], - coeff=coeff) - - -class SphericalGaussian(AbstractFunction): - r""" Stores a 3-dimensional spherical gaussian function: - - .. math:: - G_{nlm}(\mathbf r) = C Y^l_m(\mathbf r - \mathbf r_0) r^n e^{-a\left| \mathbf r - \mathbf r_0 \right|^2} - - where *C* is ``self.coeff``, *a* is ``self.exp``, and :math:`\mathbf r_0` is ``self.center``, - *(n,l,m)* are given by ``(self.n, self.l, self.m)``, and :math:`Y^l_m(\mathbf{r})` is a - spherical harmonic. - - Note: - self.SPHERE_TO_CART stores expansion coefficients for spherical gaussians in terms of - cartesian gaussians. They are taken from the reference below. - - References: - Schlegel and Frisch. Transformation between cartesian and pure spherical harmonic gaussians. - Int J Quantum Chem 54, 83-87 (1995). doi:10.1002/qua.560540202 - """ - - ndims = 3 - - def __init__(self, center, exp, n, l, m, coeff=None): - self.center = center - assert len(self.center) == 3 - self.exp = exp - self.n = n - self.l = l - self.m = m - - if coeff is None: - self.coeff = 1.0 - self.normalize() + result = self.coeff * np.exp(-self.alpha * r2) + if _getvals: + return result, disp, r2 # extra quantities so that they don't need to be recomputed else: - self.coeff = coeff + return result def __repr__(self): - return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " + return ("<{ndim}-D gaussian (coeff: {coeff:4.2f}, " "exponent: {exp:4.2f}, " - "(n,l,m) = {qnums}").format( - center=self.center, exp=self.exp, coeff=self.coeff, - qnums=(self.n, self.l, self.m)) - - def __call__(self, coords, _include_spherical=True): - """ Evaluate this function at the given coordinates. - - Can be called either with a 1D column (e.g., ``[1,2,3]*u.angstrom ``) or - an ARRAY of coordinates (``[[0,0,0],[1,1,1]] * u.angstrom``) - - Args: - _include_spherical (bool): include the contribution from the non-exponential parts - (for computational efficiency, this can sometimes omitted now and included later) - - Args: - coords (Vector[length]): 3D Coordinates or list of 3D coordinates - - Returns: - Scalar: function value(s) at the passed coordinates - """ - raise NotImplementedError - - -class AtomicBasisFunction(Orbital): - def __init__(self, atom, n=None, l=None, m=None, cart=None, primitives=None): - """ Initialization: - - Args: - atom (moldesign.Atom): The atom this basis function belongs to - index (int): the index of this basis function (it is stored as - ``wfn.basis[self.index]``) - n (int): principal quantum number (``n>=1``) - l (int): total angular momentum quantum number (``l<=n-1``) - m (int): z-angular momentum quantum number (optional - - for spherical sets only; ``|m|<=l``) - cart (str): cartesian component (optional; for cartesian sets only) - primitives (List[PrimitiveBase]): List of primitives, if available - """ - self._atom = atom - self.atom_index = atom.index - self.n = n - self.l = l - self.m = m - self.primitives = primitives - if cart is not None: - assert self.m is None, 'Both cartesian and spherical components passed!' - assert len(cart) == self.l, \ - 'Angular momentum does not match specified component %s' % cart - for e in cart: - assert e in 'xyz' - self.cart = ''.join(sorted(cart)) - else: - self.cart = None - - # These quantities can't be defined until we assemble the entire basis - self.coeffs = None - self.molecule = atom.molecule - self.basis = None - self.occupation = None - self.wfn = None - - @property - def atom(self): - """ moldesign.Atom: the atom this basis function belongs to - - We get the atom via an indirect reference, making it easier to copy the wavefunction - """ - if self.wfn is not None: - return self.wfn.mol.atoms[self.atom_index] - else: - return self._atom - - @property - def num_primitives(self): - return len(self.primitives) - - @property - def norm(self): - """ Calculate this orbital's norm - - Returns: - float: norm :math:`` - """ - norm = 0.0 - for p1 in self.primitives: - for p2 in self.primitives: - norm += p1.overlap(p2) - return norm - - def normalize(self): - """ Scale primitive coefficients to normalize this basis function - """ - prefactor = old_div(1.0,np.sqrt(self.norm)) - for primitive in self.primitives: - primitive *= prefactor + "center: {center}>").format( + ndim=self.ndim, + center=self.center, exp=self.alpha, + coeff=self.coeff) @property - def orbtype(self): - """ A string describing the orbital's angular momentum state. + def ndim(self): + return len(self.center) + num_dimensions = ndims = ndim - Examples: - >>> AtomicBasisFunction(n=1, l=0).orbtype - 's' - >>> AtomicBasisFunction(n=2, l=1, cart='y').orbtype - 'py' - >>> AtomicBasisFunction(n=3, l=2, m=0).orbtype - 'd(z^2)' + def __mul__(self, other): + """ Returns a new gaussian centroid. """ - if self.l == 0: t = 's' - elif self.cart is not None: t = SHELLS[self.l] + self.cart - else: t = SPHERICALNAMES[self.l, self.m] - return t + a1 = self.alpha + a2 = other.alpha + x1 = self.center + x2 = other.center + alpha = a1 + a2 + center = (a1*x1 + a2*x2)/(a1+a2) + prefactor = self.coeff * other.coeff + for i in range(self.ndim): + prefactor *= np.exp(-(a1*x1[i]**2 + a2*x2[i]**2) + + alpha * center[i]**2) + + return Gaussian(center, alpha, coeff=prefactor, normalized=False) @property - def aotype(self): - """ A string describing the orbital's state. - - Examples: - >>> AtomicBasisFunction(n=1, l=0).aotype - '1s' - >>> AtomicBasisFunction(n=2, l=1, cart='y').aotype - '2py' - >>> AtomicBasisFunction(n=3, l=2, m=0).aotype - '3d(z^2)' + def integral(self): + r"""Integral of this this gaussian over all N-dimensional space """ - t = self.orbtype - if self.n: return '%s%s' % (self.n, t) - else: return t - - def __str__(self): - return 'AO ' + self.name - - @property - def name(self): - try: - return '%s on atom %s' % (self.aotype, self.atom.name) - except: - return 'Basis Fn' - - def __repr__(self): - return '<%s %s>' % (self.__class__.__name__, self.name) + return self.coeff * (np.pi/self.alpha)**(self.ndim/2.0) -# Precompute odd factorial values (N!!) -_ODD_FACTORIAL = {0: 1} # by convention -_ofact = 1 -for _i in range(1, 20, 2): - _ofact *= _i - _ODD_FACTORIAL[_i] = float(_ofact) + def to_cart(self): + from . import CartesianGaussian + return CartesianGaussian(self.center, self.alpha, (0,0,0), + self.coeff, normalized=False) def cart_to_powers(s): diff --git a/moldesign/orbitals/orbitals.py b/moldesign/orbitals/orbitals.py index 99f71be..d703805 100644 --- a/moldesign/orbitals/orbitals.py +++ b/moldesign/orbitals/orbitals.py @@ -44,7 +44,7 @@ (3, -3): 'f(3yx^2-y^3)', (3, -2): 'f(xyz)', (3, -1): 'f(yz^2)', (3, 0): 'f(z^3)', (3, 1): 'f(xz^2)', (3, 2): 'f(zx^2-zy^2)', (3, 3): 'f(x^3-3xy^2)'} -ANGULAR_NAME_TO_COMPONENT = {'': (0, 0), 'x': (1, -1), 'y': (1, 1), 'z': (1, 0), +ANGULAR_NAME_TO_COMPONENT = {'': (0, 0), 'x': (1, 1), 'y': (1, -1), 'z': (1, 0), 'z^2': (2, 0), 'xy': (2, -2), 'yz': (2, -1), 'xz': (2, 1), 'x^2-y^2': (2, 2), 'zx^2-zy^2': (3, 2), 'xyz': (3, -2), 'z^3': (3, 0), diff --git a/moldesign/orbitals/primitives.py b/moldesign/orbitals/primitives.py new file mode 100644 index 0000000..34eceb4 --- /dev/null +++ b/moldesign/orbitals/primitives.py @@ -0,0 +1,164 @@ +from __future__ import print_function, absolute_import, division +from future.builtins import * +from future import standard_library +standard_library.install_aliases() + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import copy +import itertools + +import numpy as np +from .. import utils + + +class Primitive(object): + """ Abstract base class for basis function primitives. All functions are assumed real. + """ + def __init__(self, coeff=None, normalized=False): + self.coeff = 1.0 # initialize with a dummy value + + if normalized or (coeff is None): + self.normalize() + + if coeff is not None: + self.coeff *= coeff + + def _get_wfn_units(self): + return NotImplementedError() + + def copy(self): + return copy.deepcopy(self) + + def normalize(self): + """ Give this function unit norm by adjusting its coefficient + """ + self.coeff /= self.norm + + def overlap(self, other, normalized=False): + r""" Overlap of this function with another: + + .. math:: + \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r + + Args: + other (Primitive): + normalized (bool): If True, return the overlap of the two NORMALIZED functions. + + Note: + This implementation computes overlap by computing the integral of the product of the + two functions. Subclasses may override this method if more efficient techniques + are applicable + + Returns: + Scalar: value of the overlap + """ + newfn = self * other + integral = newfn.integral + if normalized: + integral /= (self.norm*other.norm) + return integral + + @property + def norm(self): + r""" The L2-Norm of this object, calculated as the square root of its self overlap. + + .. math:: + \sqrt{\int \left| G(\mathbf r) \right|^2 d^N \mathbf r} + """ + return np.sqrt(self.overlap(self)) + + def iterprimitives(self): + yield self + + +class PrimitiveSum(object): + """ Stores a linear combination of primitive functions. + + Args: + primitives (List[Primitives]): List of primitives, if available + """ + __getitem__ = utils.Alias('primitives.__getitem__') + __len__ = utils.Alias('primitives.__len__') + __iter__ = utils.Alias('primitives.__iter__') + iterprimitives = __iter__ + + def __init__(self, primitives): + self.primitives = primitives + self._wfn_units = self.primitives[0]._get_wfn_units() + + def _get_wfn_units(self): + return self._wfn_units + + def __call__(self, coords): + outvals = np.zeros(len(coords)) * self._get_wfn_units() + for primitive in self.primitives: + outvals += primitive(coords) + return outvals + + @property + def num_primitives(self): + return len(self.primitives) + + @property + def center(self): + c = self.primitives[0].center.copy() + for p in self: + if (p.center != c).any(): + raise ValueError("This is a sum over multiple centers") + else: + return c + + @property + def norm(self): + """ Scalar: :math:`\sqrt{}` + """ + norm = 0.0 + for p1 in self.primitives: + for p2 in self.primitives: + norm += p1.overlap(p2) + return np.sqrt(norm) + + def normalize(self): + """ Scale primitive coefficients to normalize this basis function + """ + prefactor = 1.0 / self.norm + for primitive in self.primitives: + primitive.coeff *= prefactor + + def overlap(self, other): + """ Calculate orbital overlap with another object + + Args: + other (AbstractFunction or Orbital): object to calculate overlaps with + """ + olap = sum(p1.overlap(p2) + for p1, p2 in itertools.product(self.iterprimitives(), + other.iterprimitives())) + return olap + + @property + def integral(self): + return sum(x.integral for x in self) + + def __str__(self): + return '%s with %d primitives' % (self.__class__.__name__, self.num_primitives) + + def __repr__(self): + return '<%s>' % self + + def __mul__(self, other): + newprims = [p1*p2 + for p1, p2 in itertools.product(self.iterprimitives(), other.iterprimitives())] + return PrimitiveSum(newprims) diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py new file mode 100644 index 0000000..f53a827 --- /dev/null +++ b/moldesign/orbitals/spherical.py @@ -0,0 +1,160 @@ +from __future__ import print_function, absolute_import, division +from future.builtins import * +from future import standard_library +standard_library.install_aliases() + +# Copyright 2017 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import numpy as np +from scipy.special import factorial + +from ..mathutils import spherical_harmonics +from .. import units as u +from .. import utils +from . import Primitive, Gaussian, CartesianGaussian, PrimitiveSum + + +class SphericalGaussian(Primitive): + r""" Stores a 3-dimensional real spherical gaussian function: + + .. math:: + G_{nlm}(\mathbf r) = C Y^l_m(\mathbf r - \mathbf r_0) r^l e^{-\alpha \left| \mathbf r - \mathbf r_0 \right|^2} + + where: + - *C* is ``self.coeff``, + - :math:`\alpha` is ``self.alpha``, + - :math:`\mathbf r_0` is ``self.center``, + - _(l,m)_ are given by ``(self.l, self.m)``, and + - :math:`Y^l_m(\mathbf{r})` are the orthonormal real-valued spherical harmonics. + + This function is evaluated as the the product of three terms that can be accessed separately: + - ``self.coeff`` - the constant prefactor _C_ + - ``self.radial_part(r)`` - :math:`e^{-a\left| \mathbf r - \mathbf r_0 \right|^2}` + - ``self.angular_part(r)`` - :math:`Y^l_m(\mathbf r - \mathbf r_0) r^l` + + References: + Schlegel and Frisch. Transformation between cartesian and pure spherical harmonic gaussians. + Int J Quantum Chem 54, 83-87 (1995). doi:10.1002/qua.560540202 + + https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics + """ + center = utils.Alias('radial_part.center') + alpha = utils.Alias('radial_part.alpha') + ndims = ndim = num_dims = 3 + + def __init__(self, center, alpha, l, m, coeff=None, normalized=True): + assert len(center) == 3, "Spherical gaussians must be 3-dimensional" + self.l = l + self.m = m + self.radial_part = Gaussian(center, alpha, coeff=1.0, normalized=False) + self._spherical_harmonic = spherical_harmonics.Y(l, m) + super().__init__(coeff=coeff, normalized=normalized) + + def __repr__(self): + return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " + "alpha: {alpha:4.2f}, " + "(l,m) = {qnums}>").format( + center=self.center, alpha=self.alpha, coeff=self.coeff, + qnums=(self.l, self.m)) + + def __mul__(self, other): + raise NotImplementedError("Cannot yet multiply spherical gaussian functions") + + def _get_wfn_units(self): + return u.MdtQuantity(self.coeff * u.default.length**self.l).units + + def __call__(self, coords): + """ Evaluate this function at the given coordinates. + + Can be called either with a 1D column (e.g., ``[1,2,3]*u.angstrom ``) or + an ARRAY of coordinates (``[[0,0,0],[1,1,1]] * u.angstrom``) + + Args: + coords (Vector[length]): 3D Coordinates or list of 3D coordinates + + Returns: + Scalar or Vector: function value(s) at the passed coordinates + """ + result, disp, r2 = self.radial_part(coords, _getvals=True) + result *= r2**(self.l/2.0) * self._spherical_harmonic(disp) + return result * self.coeff + + def overlap(self, other, normalized=False): + r""" Overlap of this function with another basis function + + .. math:: + \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r + + Args: + other (Primitive): + normalized (bool): If True, return the overlap of the two NORMALIZED functions. + + Returns: + Scalar: value of the overlap + """ + from ..data import ODD_FACTORIAL + + if (self.center != other.center).any(): + return self.to_cart().overlap(other) + + elif self.l != other.l or self.m != other.m: + return 0.0 + + else: + newexp = self.alpha + other.alpha + power = self.l + other.l + 2 + + # In this case, radial part integrates to 1, so we just + # integrate r^(2+l) exp(-(a1+a2) r^2) from 0 to infinity + if power % 2 == 0: + n = power // 2 + val = np.sqrt(np.pi/newexp) * ODD_FACTORIAL[2*n-1] / (newexp**n * 2**(n+1)) + else: + n = (power - 1) // 2 + val = factorial(n, True) / (2 * newexp**(n+1)) + + val *= self.coeff * other.coeff + + if normalized: + return val / (self.norm * other.norm) + else: + return val + + def angular_part(self, coords): + if len(coords.shape) > 1: + axis = 1 + else: + axis = None + + disp = coords - self.center + prd = disp*disp + r2 = prd.sum(axis=axis) + return r2**(self.l/2.0) * self._spherical_harmonic(disp) + + def to_cart(self): + """ Convert this function to a linear combination of cartesian functions + + Returns: + List[CartesianGaussian]: cartesian components of this function + """ + carts = [CartesianGaussian(self.center, self.alpha, + cart_fn.powers, self.coeff * cart_fn.coeff, + normalized=False) + for cart_fn in spherical_harmonics.SPHERE_TO_CART[self.l, self.m]] + return PrimitiveSum(carts) + + @property + def integral(self): + return self.to_cart().integral diff --git a/moldesign/orbitals/wfn.py b/moldesign/orbitals/wfn.py index f7ae967..c7daa14 100644 --- a/moldesign/orbitals/wfn.py +++ b/moldesign/orbitals/wfn.py @@ -109,7 +109,9 @@ def align_orbital_phases(self, other, assert_same=True): """ for orbtype in self.orbitals: if orbtype not in other.orbitals: - if assert_same: assert False, '%s has orbital type %s, but %s does not.' % (self, orbtype, other) + if assert_same: + assert False, '%s has orbital type %s, but %s does not.' % ( + self, orbtype, other) else: continue self.orbitals[orbtype].align_phases(other.orbitals[orbtype]) diff --git a/moldesign/units/constants.py b/moldesign/units/constants.py index 31cd5ce..c8f2ab3 100644 --- a/moldesign/units/constants.py +++ b/moldesign/units/constants.py @@ -42,8 +42,9 @@ electron_charge = q_e = ureg.elementary_charge # useful units -fs = femtoseconds = ureg.fs -ps = picoseconds = ureg.ps +fs = femtoseconds = ureg.femtosecond +ps = picoseconds = ureg.picosecond +ns = nanoseconds = ureg.nanosecond eV = electronvolts = ureg.eV kcalpermol = ureg.kcalpermol gpermol = ureg.gpermol @@ -53,7 +54,7 @@ amu = da = dalton = ureg.amu kelvin = ureg.kelvin nm = ureg.nanometers -ang = angstrom = ureg.ang +angstrom = ureg.angstrom molar = ureg.mole / ureg.liter debye = ureg.debye @@ -64,4 +65,4 @@ def_mass = amu def_momentum = def_mass * def_vel def_force = def_momentum / def_time -def_energy = eV \ No newline at end of file +def_energy = eV diff --git a/moldesign/units/quantity.py b/moldesign/units/quantity.py index 2e5f98f..1c10e2f 100644 --- a/moldesign/units/quantity.py +++ b/moldesign/units/quantity.py @@ -30,8 +30,9 @@ # Set up pint's unit definitions ureg = UnitRegistry() -unit_def_file = join(abspath(dirname(__file__)), '../_static_data/pint_atomic_units.txt') +unit_def_file = join(abspath(dirname(__file__)), '..', '_static_data','pint_atomic_units.txt') ureg.load_definitions(unit_def_file) +ureg.default_system = 'nano' set_application_registry(ureg) @@ -42,6 +43,10 @@ class MdtUnit(ureg.Unit): def __reduce__(self): return _get_unit, (str(self),) + @property + def units(self): + return self + def convert(self, value): """ Returns quantity converted to these units diff --git a/moldesign/units/tools.py b/moldesign/units/tools.py index 9367839..c0d82cb 100644 --- a/moldesign/units/tools.py +++ b/moldesign/units/tools.py @@ -17,9 +17,12 @@ # See the License for the specific language governing permissions and # limitations under the License. +from .. import utils + from .quantity import * + def unitsum(iterable): """ Faster method to compute sums of iterables if they're all in the right units @@ -55,6 +58,18 @@ def dot(a1, a2): return a1.dot(a2) +@utils.args_from(np.linspace) +def linspace(start, stop, **kwargs): + u1 = getattr(start, 'units', ureg.dimensionless) + u2 = getattr(stop, 'units', ureg.dimensionless) + if u1 == u2 == ureg.dimensionless: + return np.linspace(start, stop, **kwargs) + else: + q1mag = start.magnitude + q2mag = stop.value_in(start.units) + return np.linspace(q1mag, q2mag, **kwargs) * start.units + + def arrays_almost_equal(a1, a2): """ Return true if arrays are almost equal up to numerical noise @@ -131,7 +146,7 @@ def get_units(q): Examples: >>> from moldesign import units >>> units.get_units(1.0 * units.angstrom) - + >>> units.get_units(np.array([1.0, 2, 3.0])) >>> # We dive on the first element of each iterable until we can determine a unit system: @@ -153,36 +168,68 @@ def get_units(q): else: if isinstance(x, str): raise TypeError('Found string data while trying to determine units') - return MdtQuantity(x).units + + q = MdtQuantity(x) + if q.dimensionless: + return ureg.dimensionless + else: + return q.units -def array(qlist, baseunit=None): +def array(qlist, defunits=False, _baseunit=None): """ Facilitates creating an array with units - like numpy.array, but it also checks - units for all components of the array + units for all components of the array. + + Note: + Unlike numpy.array, these arrays must have numeric type - this routine will + raise a ValueError if a non-square array is passed. Args: qlist (List[MdtQuantity]): List-like object of quantity objects - baseunit (MdtUnit) unit to standardize with + defunits (bool): if True, convert the array to the default units Returns: - MdtQuantity: array with standardized units + MdtQuantity: ndarray-like object with standardized units + + Raises: + DimensionalityError: if the array has inconsistent units + ValueError: if the array could not be converted to a square numpy array """ + from . import default + if hasattr(qlist, 'units') and hasattr(qlist, 'magnitude'): return MdtQuantity(qlist) - if baseunit is None: - baseunit = get_units(qlist) - try: - if baseunit == 1.0: - return np.array(qlist) - except DimensionalityError: - pass + if _baseunit is None: + _baseunit = get_units(qlist) + if _baseunit.dimensionless: + return _make_nparray(qlist) + if defunits: + _baseunit = default.get_default(_baseunit) + + if hasattr(qlist, 'to'): # if already a quantity, just convert and return + return qlist.to(_baseunit) + try: # try to create a quantity + return _baseunit * [array(item, _baseunit=_baseunit).value_in(_baseunit) for item in qlist] + except TypeError: # if here, one or more objects cannot be converted to quantities + raise DimensionalityError(_baseunit, ureg.dimensionless, + extra_msg='Object "%s" does not have units' % qlist) + + +def _make_nparray(q): + """ Turns a list of dimensionless numbers into a numpy array. Does not permit object arrays + """ + if hasattr(q, 'units'): + return q.value_in(ureg.dimensionless) try: - newlist = [array(item, baseunit=baseunit).value_in(baseunit) for item in qlist] - return baseunit * newlist - except TypeError as exc: - return qlist.to(baseunit) + arr = np.array([_make_nparray(x) for x in q]) + if arr.dtype == 'O': + raise ValueError("Could not create numpy array of numeric data - is your input square?") + else: + return arr + except TypeError: + return q #@utils.args_from(np.broadcast_to) diff --git a/moldesign/units/unitsystem.py b/moldesign/units/unitsystem.py index d57c2d7..46413e0 100644 --- a/moldesign/units/unitsystem.py +++ b/moldesign/units/unitsystem.py @@ -25,6 +25,17 @@ class UnitSystem(object): In MDT, many methods will automatically convert output using the UnitSystem at ``moldesign.units.default`` + + Args: + length (MdtUnit): length units + mass (MdtUnit): mass units + time (MdtUnit): time units + energy (MdtUnit): energy units + temperature (MdtUnit): temperature units (default: kelvin) + force (MdtUnit): force units (default: energy/length) + momentum (MdtUnit): momentum units (default: mass * length / time) + angle (MdtUnit): angle units (default: radians) + charge (MdtUnit): charge units (default: fundamental charge) """ def __init__(self, length, mass, time, energy, temperature=kelvin, @@ -75,16 +86,26 @@ def convert(self, quantity): """ Convert a quantity into this unit system. Args: - quantity (MdtQuantity): quantity to convert + quantity (MdtQuantity or MdtUnit): quantity to convert """ baseunit = self.get_baseunit(quantity) - if isinstance(baseunit, int): - assert baseunit == 1 + if baseunit == ureg.dimensionless: return quantity * ureg.dimensionless else: result = quantity.to(baseunit) return result + def get_default(self, q): + """ Return the default unit system for objects with these dimensions + + Args: + q (MdtQuantity or MdtUnit): quantity to get default units for + + Returns: + MdtUnit: Proper units for this quantity + """ + return self.get_baseunit(1.0 * q).units + def convert_if_possible(self, quantity): if isinstance(quantity, MdtQuantity): return self.convert(quantity) @@ -109,16 +130,27 @@ def get_baseunit(self, quantity): try: q = quantity[0] except (TypeError, StopIteration): + if isinstance(quantity, (int, float, complex)): + return ureg.dimensionless raise TypeError('This type of object cannot have physical units') if isinstance(q, str): raise TypeError('This type of object cannot have physical units') try: return self.get_baseunit(q) except (IndexError, TypeError): # Assume dimensionless - return 1 - baseunit = 1 - - # Factor out energy units (this doesn't work for things like energy/length) + return ureg.dimensionless + baseunit = ureg.dimensionless + + # Factor out force units + if self._force: + if '[length]' in dims and '[mass]' in dims and '[time]' in dims: + while dims['[length]'] >= 1 and dims['[mass]'] >= 1 and dims['[time]'] <= -2: + baseunit *= self['force'] + dims['[length]'] -= 1 + dims['[mass]'] -= 1 + dims['[time]'] += 2 + + # Factor out energy units if '[length]' in dims and '[mass]' in dims and '[time]' in dims: while dims['[length]'] >= 1 and dims['[mass]'] >= 1 and dims['[time]'] <= -2: baseunit *= self['energy'] @@ -126,6 +158,15 @@ def get_baseunit(self, quantity): dims['[mass]'] -= 1 dims['[time]'] += 2 + # Factor out momentum units + if self._momentum: + if '[length]' in dims and '[mass]' in dims and '[time]' in dims: + while dims['[length]'] >= 1 and dims['[mass]'] >= 1 and dims['[time]'] <= -1: + baseunit *= self['momentum'] + dims['[length]'] -= 1 + dims['[mass]'] -= 1 + dims['[time]'] += 1 + if '[current]' in dims: dims.setdefault('[charge]', 0) dims.setdefault('[time]', 0) @@ -141,7 +182,7 @@ def get_baseunit(self, quantity): baseunit *= self[unit]**dims[unit] except AttributeError: baseunit *= ureg[unit]**dims[unit] - return baseunit + return baseunit.units default = UnitSystem(length=angstrom, mass=amu, time=fs, energy=eV) diff --git a/moldesign/utils/descriptors.py b/moldesign/utils/descriptors.py index 16d3bc1..206ef6a 100644 --- a/moldesign/utils/descriptors.py +++ b/moldesign/utils/descriptors.py @@ -52,6 +52,14 @@ def _method_getter(s, *args, **kwargs): return _method_getter +class IndexView(object): + def __init__(self, attr, index): + self.attr = attr + self.index = index + + def __get__(self, instance, owner): + return getattr(instance, self.attr)[self.index] + class Synonym(object): """ An attribute (class or intance) that is just a synonym for another. diff --git a/requirements.txt b/requirements.txt index 29b2964..fb48154 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,7 +6,7 @@ funcsigs numpy >= 1.12 pathlib2 ; python_version < '3.3' parmed >= 2.7.3 -pint >= 0.8 +pint >= 0.8.1 pyccc >= 0.7.9 pyyaml requests