From 8c63e57d621221332b1d0998b29ada2e73a10b13 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 27 Jun 2017 19:58:43 -0700 Subject: [PATCH 01/33] Stash changes --- moldesign/interfaces/pyscf_interface.py | 2 +- moldesign/orbitals/gaussians.py | 46 +++++++++++++------------ 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/moldesign/interfaces/pyscf_interface.py b/moldesign/interfaces/pyscf_interface.py index efc0258..5655380 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 diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index d66cd2d..7e3fc16 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -242,10 +242,10 @@ 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} + G_{nlm}(\mathbf r) = C Y^l_m(\mathbf r - \mathbf r_0) r^l 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 + *(l,m)* are given by ``(self.l, self.m)``, and :math:`Y^l_m(\mathbf{r})` is a spherical harmonic. Note: @@ -259,11 +259,10 @@ class SphericalGaussian(AbstractFunction): ndims = 3 - def __init__(self, center, exp, n, l, m, coeff=None): + def __init__(self, center, exp, l, m, coeff=None): self.center = center assert len(self.center) == 3 self.exp = exp - self.n = n self.l = l self.m = m @@ -276,9 +275,9 @@ def __init__(self, center, exp, n, l, m, coeff=None): def __repr__(self): return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " "exponent: {exp:4.2f}, " - "(n,l,m) = {qnums}").format( + "(l,m) = {qnums}").format( center=self.center, exp=self.exp, coeff=self.coeff, - qnums=(self.n, self.l, self.m)) + qnums=(self.l, self.m)) def __call__(self, coords, _include_spherical=True): """ Evaluate this function at the given coordinates. @@ -300,20 +299,23 @@ def __call__(self, coords, _include_spherical=True): class AtomicBasisFunction(Orbital): - def __init__(self, atom, n=None, l=None, m=None, cart=None, primitives=None): - """ Initialization: + """ Stores an atomic basis function. - 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 - """ + 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``) - this is useful only as a metadata object + 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): self.atom = atom self.n = n self.l = l @@ -343,18 +345,18 @@ def norm(self): """ Calculate this orbital's norm Returns: - float: norm :math:`` + float: norm :math:`\sqrt{}` """ norm = 0.0 for p1 in self.primitives: for p2 in self.primitives: norm += p1.overlap(p2) - return norm + return np.sqrt(norm) def normalize(self): """ Scale primitive coefficients to normalize this basis function """ - prefactor = old_div(1.0,np.sqrt(self.norm)) + prefactor = 1.0 / self.norm for primitive in self.primitives: primitive *= prefactor From efaa71b4a1fc01e7edf944b04ebce5c1ec2a92b9 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 27 Jun 2017 23:30:03 -0700 Subject: [PATCH 02/33] Add spherical harmonics up to f functions --- moldesign/_tests/test_mathutils.py | 49 ++++++ moldesign/mathutils/spherical_harmonics.py | 187 +++++++++++++++++++++ 2 files changed, 236 insertions(+) create mode 100644 moldesign/mathutils/spherical_harmonics.py diff --git a/moldesign/_tests/test_mathutils.py b/moldesign/_tests/test_mathutils.py index de63195..230a2df 100644 --- a/moldesign/_tests/test_mathutils.py +++ b/moldesign/_tests/test_mathutils.py @@ -4,6 +4,8 @@ from moldesign import mathutils from moldesign import units as u +from moldesign.mathutils import spherical_harmonics as sh + from . import helpers @@ -97,6 +99,53 @@ 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),)) + + +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/mathutils/spherical_harmonics.py b/moldesign/mathutils/spherical_harmonics.py new file mode 100644 index 0000000..d208e7f --- /dev/null +++ b/moldesign/mathutils/spherical_harmonics.py @@ -0,0 +1,187 @@ +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'] + +SQ2INV = np.sqrt(1/2.0) + + +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): + """ A *real* spherical harmonic + + 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 + else: + vplus = sph_harm(self.m, self.l, phi, theta) + vminus = sph_harm(-self.m, self.l, phi, theta) + + if self.m < 0: + return -(SQ2INV * (self._posneg * vplus + vminus)).imag + else: + return (SQ2INV * (self._posneg * vplus + vminus)).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 + + +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 sqrt_x_over_pi(num, denom): + return np.sqrt(num / (denom*np.pi)) + + +SPHERE_TO_CART = {(0, 0): Cart(0, 0, 0, sqrt_x_over_pi(1, 4)), + + ############ l=1 ############ + (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)), + + ############ l=2 ############ + (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)]), + + ############ l=3 ############ + (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 From 54f29c81898d9fc417f60959ebc3163f3f68b05c Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Wed, 28 Jun 2017 00:16:55 -0700 Subject: [PATCH 03/33] Hook up SphericalGaussians, fix real harmonics signs --- moldesign/mathutils/spherical_harmonics.py | 15 ++++---- moldesign/orbitals/gaussians.py | 44 +++++++++++++++------- 2 files changed, 38 insertions(+), 21 deletions(-) diff --git a/moldesign/mathutils/spherical_harmonics.py b/moldesign/mathutils/spherical_harmonics.py index d208e7f..cf92b09 100644 --- a/moldesign/mathutils/spherical_harmonics.py +++ b/moldesign/mathutils/spherical_harmonics.py @@ -61,8 +61,8 @@ def __call__(self, coords): if self.m == 0: return (sph_harm(self.m, self.l, phi, theta)).real else: - vplus = sph_harm(self.m, self.l, phi, theta) - vminus = sph_harm(-self.m, self.l, phi, theta) + vplus = sph_harm(abs(self.m), self.l, phi, theta) + vminus = sph_harm(-abs(self.m), self.l, phi, theta) if self.m < 0: return -(SQ2INV * (self._posneg * vplus + vminus)).imag @@ -136,15 +136,16 @@ 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)), - ############ l=1 ############ + ############ 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)), - ############ l=2 ############ - (2, -2): Cart(1, 1, 0, -sqrt_x_over_pi(15, 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), @@ -155,12 +156,12 @@ def sqrt_x_over_pi(num, denom): [(2, 0, 0, 1.0), (0, 2, 0, -1.0)]), - ############ l=3 ############ + ############ 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, -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), diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 7e3fc16..e0cbee7 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -18,7 +18,7 @@ # 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 collections import numpy as np from .orbitals import Orbital, SHELLS, SPHERICALNAMES @@ -184,7 +184,7 @@ def __mul__(self, other): a = self.exp b = other.exp exp = a + b - center = old_div((a*self.center + b*other.center),(a+b)) + center = (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) @@ -208,7 +208,7 @@ def integral(self): 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)) + integ = (np.pi/self.exp)**(self.ndim/2.0) for p in self.powers: if p == 0: # no contribution continue @@ -217,7 +217,7 @@ def integral(self): elif p < 0: raise ValueError('Powers must be positive or 0') else: - integ *= old_div(_ODD_FACTORIAL[p-1],(2 ** (p+1))) + integ *= _ODD_FACTORIAL[p-1]/(2 ** (p+1)) return self.coeff * integ @@ -239,27 +239,27 @@ def Gaussian(center, exp, coeff=1.0): class SphericalGaussian(AbstractFunction): - r""" Stores a 3-dimensional spherical gaussian function: + 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^{-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``, - *(l,m)* are given by ``(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. + *(l,m)* are given by ``(self.l, self.m)``, and :math:`Y^l_m(\mathbf{r})` are the _real_ + spherical harmonics. 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 """ ndims = 3 def __init__(self, center, exp, l, m, coeff=None): + from moldesign.mathutils import spherical_harmonics + self.center = center assert len(self.center) == 3 self.exp = exp @@ -272,6 +272,8 @@ def __init__(self, center, exp, l, m, coeff=None): else: self.coeff = coeff + self._spherical_harmonic = spherical_harmonics.Y(l, m) + def __repr__(self): return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " "exponent: {exp:4.2f}, " @@ -295,7 +297,19 @@ def __call__(self, coords, _include_spherical=True): Returns: Scalar: function value(s) at the passed coordinates """ - raise NotImplementedError + if len(coords.shape) > 1: + axis = 1 + else: + axis = None + + disp = coords - self.center + prd = disp*disp + r2 = prd.sum(axis=axis) + + result = self.coeff * np.exp(-self.exp * r2) + if _include_spherical: + result *= r2**(self.l/2.0) * self._spherical_harmonic(coords) + return result class AtomicBasisFunction(Orbital): @@ -390,8 +404,10 @@ def aotype(self): '3d(z^2)' """ t = self.orbtype - if self.n: return '%s%s' % (self.n, t) - else: return t + if self.n: + return '%s%s' % (self.n, t) + else: + return t def __str__(self): return 'AO ' + self.name From cc2f5d1fd3be53122588c68a7779d9bcc49d5e14 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 30 Jun 2017 23:20:29 -0700 Subject: [PATCH 04/33] Add more electronic wavefunction test fixtures --- moldesign/_tests/molecule_fixtures.py | 32 +++++++++++++++++-- moldesign/_tests/object_fixtures.py | 8 ++--- .../test_atom_bond_computed_properties.py | 12 +++---- moldesign/_tests/test_copies.py | 14 ++++---- moldesign/_tests/test_qm_xfaces.py | 6 ++-- moldesign/_tests/test_trajectory.py | 2 +- 6 files changed, 50 insertions(+), 24 deletions(-) diff --git a/moldesign/_tests/molecule_fixtures.py b/moldesign/_tests/molecule_fixtures.py index c94bdbd..1f7baf4 100644 --- a/moldesign/_tests/molecule_fixtures.py +++ b/moldesign/_tests/molecule_fixtures.py @@ -289,7 +289,7 @@ def protein_default_amber_forcefield(cached_protein_with_default_amber_ff): @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']) @@ -297,5 +297,31 @@ def cached_h2_rhfwfn(): @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() + + +@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 + + +@pytest.fixture +def h2_rhf_augccpvdz(cached_h2_rhf_augccpvdz): + return cached_h2_rhf_augccpvdz.copy() + + +@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 + + +@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..eef3e7d 100644 --- a/moldesign/_tests/object_fixtures.py +++ b/moldesign/_tests/object_fixtures.py @@ -104,13 +104,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') 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_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_qm_xfaces.py b/moldesign/_tests/test_qm_xfaces.py index decdcb0..3037bde 100644 --- a/moldesign/_tests/test_qm_xfaces.py +++ b/moldesign/_tests/test_qm_xfaces.py @@ -139,9 +139,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 +149,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) From b9249bb1cd0e812b5fe6f00d38d11c454d0c13a5 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 30 Jun 2017 23:21:13 -0700 Subject: [PATCH 05/33] MDT and PySCF amplitude grids The two programs more or less agree up to a multiplicative constant, currently. Big disagreements around small values, which is probably fine --- moldesign/interfaces/pyscf_interface.py | 2 +- moldesign/models/pyscf.py | 6 +- moldesign/orbitals/__init__.py | 1 + moldesign/orbitals/atomic_basis_fn.py | 152 +++++++++++++++++++ moldesign/orbitals/basis.py | 54 +++---- moldesign/orbitals/gaussians.py | 194 ++++++------------------ 6 files changed, 229 insertions(+), 180 deletions(-) create mode 100644 moldesign/orbitals/atomic_basis_fn.py diff --git a/moldesign/interfaces/pyscf_interface.py b/moldesign/interfaces/pyscf_interface.py index 5655380..cd6843c 100644 --- a/moldesign/interfaces/pyscf_interface.py +++ b/moldesign/interfaces/pyscf_interface.py @@ -132,6 +132,6 @@ def basis_values(mol, basis, coords, coeffs=None, positions=None): if coeffs is None: return aovals else: - return aovals.dot(coeffs) + return aovals.dot(coeffs.T) diff --git a/moldesign/models/pyscf.py b/moldesign/models/pyscf.py index 60f3cda..728c97b 100644 --- a/moldesign/models/pyscf.py +++ b/moldesign/models/pyscf.py @@ -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) # very unsure if these units are right! num_contractions = pmol.bas_nctr(ishell) coeffs = pmol.bas_ctr_coeff(ishell) @@ -450,11 +450,11 @@ 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, + exp, l, m, coeff=coeff[ictr]) for exp, coeff in zip(exps, coeffs)] bfs.append(orbitals.AtomicBasisFunction(atom, diff --git a/moldesign/orbitals/__init__.py b/moldesign/orbitals/__init__.py index 43f0b03..f589c6e 100644 --- a/moldesign/orbitals/__init__.py +++ b/moldesign/orbitals/__init__.py @@ -5,5 +5,6 @@ def toplevel(o): from .gaussians import * from .orbitals 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..d388ac6 --- /dev/null +++ b/moldesign/orbitals/atomic_basis_fn.py @@ -0,0 +1,152 @@ +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 .orbitals import Orbital, SHELLS, SPHERICALNAMES + + +class AtomicBasisFunction(Orbital): + """ 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): + 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 + + def __call__(self, coords): + outvals = np.zeros(len(coords)) + for primitive in self.primitives: + outvals += primitive(coords) + return outvals + + @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:`\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 *= prefactor + + @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..f437348 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))) + 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): + + 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/gaussians.py b/moldesign/orbitals/gaussians.py index 5e4daa3..2f795b6 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -16,11 +16,8 @@ # 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 collections import numpy as np -from .orbitals import Orbital, SHELLS, SPHERICALNAMES - class AbstractFunction(object): """ Abstract base class for basis functions @@ -79,24 +76,22 @@ class CartesianGaussian(AbstractFunction): References: Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. + + 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. """ 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 @@ -108,7 +103,7 @@ def __init__(self, center, exp, powers, coeff=None): else: self.coeff = coeff - self._cartesian = (self.powers != 0).any() + self.shell = sum(self.powers) def __repr__(self): return ("<{ndim}-D Gaussian (cart) (coeff: {coeff:4.2f}, " @@ -130,7 +125,7 @@ def angular(self): """ return self.powers.sum() - def __call__(self, coords, _include_cartesian=True): + def __call__(self, coords, _include_angular=True): """ Evaluate this function at the given coordinates. Can be called either with a 1D column (e.g., ``[1,2,3]*u.angstrom ``) or @@ -163,11 +158,22 @@ def __call__(self, coords, _include_cartesian=True): 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() ) + if self.shell > 0 and _include_angular: + result *= self.angular_part(coords) return result + 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 + return np.product(disp.magnitude**self.powers, axis=axis)*disp.units ** self.powers.sum() + def __mul__(self, other): """ Returns product of two gaussian functions, which is also a gaussian @@ -177,6 +183,8 @@ def __mul__(self, other): Returns: CartesianGaussian: product gaussian """ + if self.center != other.center: + raise NotImplementedError() # convert widths to prefactor form a = self.exp @@ -279,15 +287,15 @@ def __repr__(self): center=self.center, exp=self.exp, coeff=self.coeff, qnums=(self.l, self.m)) - def __call__(self, coords, _include_spherical=True): + def __call__(self, coords, _include_angular=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) + _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]): 3D Coordinates or list of 3D coordinates @@ -305,134 +313,22 @@ def __call__(self, coords, _include_spherical=True): r2 = prd.sum(axis=axis) result = self.coeff * np.exp(-self.exp * r2) - if _include_spherical: - result *= r2**(self.l/2.0) * self._spherical_harmonic(coords) + if _include_angular: + result *= r2**(self.l/2.0) * self._spherical_harmonic(disp) return result - -class AtomicBasisFunction(Orbital): - """ 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``) - this is useful only as a metadata object - 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): - 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:`\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 *= prefactor - - @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) + def angular_part(self, coords): + if len(coords.shape) > 1: + axis = 1 else: - return t + axis = None - def __str__(self): - return 'AO ' + self.name + disp = coords - self.center + prd = disp*disp + r2 = prd.sum(axis=axis) - @property - def name(self): - try: - return '%s on atom %s' % (self.aotype, self.atom.name) - except: - return 'Basis Fn' + return r2**(self.l/2.0) * self._spherical_harmonic(disp) - def __repr__(self): - return '<%s %s>' % (self.__class__.__name__, self.name) # Precompute odd factorial values (N!!) _ODD_FACTORIAL = {0: 1} # by convention From a5b4e7d41841ffa265578828de67a71a243c7d61 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Sat, 1 Jul 2017 00:45:22 -0700 Subject: [PATCH 06/33] Wfn grid tests --- moldesign/_tests/test_wfn.py | 61 ++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 moldesign/_tests/test_wfn.py diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py new file mode 100644 index 0000000..b321a3c --- /dev/null +++ b/moldesign/_tests/test_wfn.py @@ -0,0 +1,61 @@ +import pytest +import numpy as np + +from .molecule_fixtures import * + +import moldesign as mdt + +from moldesign.interfaces.pyscf_interface import basis_values + + +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)) + np.testing.assert_allclose(vals_plusminus[:,0], + (vals_b0[:,0] + vals_b0[:,1])/np.sqrt(2)) + np.testing.assert_allclose(vals_plusminus[:,1], + (vals_b0[:,0] - vals_b0[:,1])/np.sqrt(2)) + + +@pytest.mark.parametrize('molkey', TESTSYSTEMS) +def test_pyscf_and_mdt_grids_are_the_same(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) + mdt_vals = mol.wfn.aobasis(randocoords) + np.testing.assert_allclose(mdt_vals, pyscf_vals) + + + From dd164f797b6dab81bc369c5ad8f1dfc33e2b5b07 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 7 Jul 2017 12:39:56 -0700 Subject: [PATCH 07/33] Limit the test fixture `import *` side effects --- moldesign/_tests/molecule_fixtures.py | 27 ++++++++++++++++++++++++++- moldesign/_tests/object_fixtures.py | 6 ++++++ moldesign/_tests/test_molecules.py | 21 +++++++++++---------- moldesign/_tests/test_objects.py | 4 ++++ moldesign/_tests/test_qm_xfaces.py | 20 -------------------- 5 files changed, 47 insertions(+), 31 deletions(-) diff --git a/moldesign/_tests/molecule_fixtures.py b/moldesign/_tests/molecule_fixtures.py index 1f7baf4..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,6 +307,7 @@ 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_rhf_sto3g(): mol = h2() # fixture is not cached, so just call it directly @@ -296,11 +316,13 @@ def cached_h2_rhf_sto3g(): return mol +@exports @pytest.fixture 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() @@ -309,11 +331,13 @@ def cached_h2_rhf_augccpvdz(): 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') @@ -322,6 +346,7 @@ def cached_acetylene_dft_631g(): 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 eef3e7d..08725d5 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) @@ -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_molecules.py b/moldesign/_tests/test_molecules.py index b50355a..7902dee 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): @@ -96,18 +97,18 @@ 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): + 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): + 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 3037bde..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) From 82e4f143ddb5ba0670166f2310a5eac0832cd4b3 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 7 Jul 2017 14:23:53 -0700 Subject: [PATCH 08/33] Test and stabilize the UnitSystem object --- moldesign/_tests/test_atoms.py | 3 +++ moldesign/_tests/test_units.py | 35 +++++++++++++++++++++++++++ moldesign/units/constants.py | 1 + moldesign/units/tools.py | 24 ++++++++++++------- moldesign/units/unitsystem.py | 44 ++++++++++++++++++++++++++++++++-- 5 files changed, 97 insertions(+), 10 deletions(-) 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_units.py b/moldesign/_tests/test_units.py index 1b69521..28051a2 100644 --- a/moldesign/_tests/test_units.py +++ b/moldesign/_tests/test_units.py @@ -57,6 +57,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] diff --git a/moldesign/units/constants.py b/moldesign/units/constants.py index 31cd5ce..13fe61c 100644 --- a/moldesign/units/constants.py +++ b/moldesign/units/constants.py @@ -44,6 +44,7 @@ # useful units fs = femtoseconds = ureg.fs ps = picoseconds = ureg.ps +ns = nanoseconds = ureg.ns eV = electronvolts = ureg.eV kcalpermol = ureg.kcalpermol gpermol = ureg.gpermol diff --git a/moldesign/units/tools.py b/moldesign/units/tools.py index 9367839..1fd5142 100644 --- a/moldesign/units/tools.py +++ b/moldesign/units/tools.py @@ -156,33 +156,41 @@ def get_units(q): return MdtQuantity(x).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 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 """ + from . import default + if hasattr(qlist, 'units') and hasattr(qlist, 'magnitude'): return MdtQuantity(qlist) - if baseunit is None: - baseunit = get_units(qlist) + if _baseunit is None: + _baseunit = get_units(qlist) try: - if baseunit == 1.0: + if _baseunit == 1.0: + # It's dimensionless, return an ndarray return np.array(qlist) except DimensionalityError: pass + if defunits and isinstance(_baseunit, MdtUnit): + _baseunit = default.get_default(_baseunit) + try: - newlist = [array(item, baseunit=baseunit).value_in(baseunit) for item in qlist] - return baseunit * newlist + # Convert all sublists to the common "baseunit" + newlist = [array(item, _baseunit=_baseunit).value_in(_baseunit) for item in qlist] + return _baseunit*newlist except TypeError as exc: - return qlist.to(baseunit) + # Otherwise, + return qlist.to(_baseunit) #@utils.args_from(np.broadcast_to) diff --git a/moldesign/units/unitsystem.py b/moldesign/units/unitsystem.py index d57c2d7..f86b1e0 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,7 +86,7 @@ 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): @@ -85,6 +96,17 @@ def convert(self, quantity): 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) @@ -118,7 +140,16 @@ def get_baseunit(self, quantity): return 1 baseunit = 1 - # Factor out energy units (this doesn't work for things like energy/length) + # 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 +157,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) From a0c9ca712e590f17d3c7cb71a082281122898895 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 7 Jul 2017 14:57:08 -0700 Subject: [PATCH 09/33] Rework and expose the VolumetricGrid class --- moldesign/_tests/test_wfn.py | 63 +++++++++++++++- moldesign/helpers/helpers.py | 34 --------- moldesign/mathutils/__init__.py | 1 + moldesign/mathutils/grids.py | 124 ++++++++++++++++++++++++++++++++ moldesign/units/tools.py | 15 ++++ moldesign/utils/descriptors.py | 8 +++ 6 files changed, 208 insertions(+), 37 deletions(-) create mode 100644 moldesign/mathutils/grids.py diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py index b321a3c..4753a92 100644 --- a/moldesign/_tests/test_wfn.py +++ b/moldesign/_tests/test_wfn.py @@ -1,15 +1,17 @@ import pytest import numpy as np -from .molecule_fixtures import * - 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 + 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 @@ -58,4 +60,59 @@ def test_pyscf_and_mdt_grids_are_the_same(molkey, request): np.testing.assert_allclose(mdt_vals, pyscf_vals) +@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 = mdt.mathutils.VolumetricGrid(*ranges, 3, 4, 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 = mdt.mathutils.VolumetricGrid(*ranges, 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] + + 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/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..c4200c9 --- /dev/null +++ b/moldesign/mathutils/grids.py @@ -0,0 +1,124 @@ +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: 25) + ypoints (int): number of grid lines in y direction (default: xpoints) + zpoints (int): number of grid lines in z direction (default: 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=25, ypoints=None, zpoints=None): + + self.points = np.array([xpoints, + ypoints if ypoints is not None else xpoints, + zpoints if zpoints is not None else xpoints], + 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 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, dtype='float'): + """ 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 + """ + ap = np.empty((self.npoints, 3), dtype=dtype) + if hasattr(self.xr, 'units'): + ap = ap * self.xr.units + + for ip, point in enumerate(self.iter_points()): + ap[ip] = point + + return ap + + +@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/units/tools.py b/moldesign/units/tools.py index 1fd5142..62c5061 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 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. From 4a72172af4ca076c617e4c5640cc112fc38aec7e Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 7 Jul 2017 16:00:10 -0700 Subject: [PATCH 10/33] Analytical and numerical norms agree for basic gaussians --- moldesign/_tests/test_gaussian_math.py | 33 +++++++++++++++++++++----- moldesign/mathutils/grids.py | 5 ++-- moldesign/orbitals/gaussians.py | 17 ++++++------- 3 files changed, 39 insertions(+), 16 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 77bb9e5..b3268b5 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -10,10 +10,10 @@ import moldesign from moldesign import units as u -registered_types = {} +from . import helpers +registered_types = {} -# 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. @@ -60,7 +60,7 @@ def cart_10d_gaussian(): @pytest.mark.parametrize('objkey', registered_types['gaussian']) @pytest.mark.screening def test_gaussian_integral_and_dimensionality(objkey, request): - g = request.getfuncargvalue(objkey) + g = request.getfixturevalue(objkey) assert g.ndims == len(g.center) intval = g.integral @@ -73,7 +73,7 @@ def test_gaussian_integral_and_dimensionality(objkey, request): @pytest.mark.parametrize('objkey', registered_types['gaussian'] + registered_types['cartesiangaussian']) 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() @@ -88,7 +88,7 @@ def test_gaussian_function_values(objkey, request): registered_types['gaussian'] + registered_types['cartesiangaussian']) @pytest.mark.screening 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): @@ -123,6 +123,7 @@ def _gfuncval(g, coord): 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') @@ -142,4 +143,24 @@ def test_cart_to_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] + + +def test_numerical_vs_analytical_norm(std_3d_gaussian): + g = std_3d_gaussian + ananorm = g.norm + + width = np.sqrt(1/g.exp) + ranges = np.ones((3,2)) * 5.0 * width + ranges[:,0] *= -1 + ranges += g.center[:, None] + grid = moldesign.mathutils.VolumetricGrid(*ranges, 25) + allpoints = grid.allpoints() + + with np.errstate(under='ignore'): + vals = g(allpoints) + + numnorm = np.sqrt(grid.dx**3 * (vals**2).sum()) + helpers.assert_almost_equal(ananorm, numnorm) + + diff --git a/moldesign/mathutils/grids.py b/moldesign/mathutils/grids.py index c4200c9..57b9d6d 100644 --- a/moldesign/mathutils/grids.py +++ b/moldesign/mathutils/grids.py @@ -91,12 +91,13 @@ def allpoints(self, dtype='float'): Matrix[shape=(self.npoints**3,3)]: x,y,z coordinate of each point on the grid """ ap = np.empty((self.npoints, 3), dtype=dtype) - if hasattr(self.xr, 'units'): - ap = ap * self.xr.units for ip, point in enumerate(self.iter_points()): ap[ip] = point + if hasattr(self.xr, 'units'): + ap = ap * self.xr.units + return ap diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 2f795b6..3c2f4a4 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -51,12 +51,12 @@ def norm(self): r""" The L2-Norm of this gaussian: .. math:: - \int \left| G(\mathbf r) \right|^2 d^N \mathbf r + \sqrt{\int \left| G(\mathbf r) \right|^2 d^N \mathbf r} """ - return self.overlap(self) + return np.sqrt(self.overlap(self)) def __str__(self): - return "%d-D %s with norm %s"%(self.ndims, self.__class__, self.norm) + return "%d-D %s with norm %s" % (self.ndims, self.__class__, self.norm) class CartesianGaussian(AbstractFunction): @@ -132,19 +132,20 @@ def __call__(self, coords, _include_angular=True): 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 _include_cartesian (bool): include the contribution from the cartesian components (for computational efficiency, this can sometimes omited now and included later) - Example: + Examples: >>> g = CartesianGaussian([0,0,0]*u.angstrom, exp=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]) - Args: - coords (Vector[length]): Coordinates or list of coordinates - Returns: Scalar: function value(s) at the passed coordinates """ @@ -183,7 +184,7 @@ def __mul__(self, other): Returns: CartesianGaussian: product gaussian """ - if self.center != other.center: + if (self.center != other.center).any(): raise NotImplementedError() # convert widths to prefactor form From ef5d5a67deb4dfc8826247c6a06ec68ac3ea3927 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 7 Jul 2017 16:58:43 -0700 Subject: [PATCH 11/33] Analytical cartesian gaussian norms agree with numerical ones --- moldesign/_tests/test_gaussian_math.py | 31 ++++++++++++++++++++++---- moldesign/orbitals/gaussians.py | 5 +++-- 2 files changed, 30 insertions(+), 6 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index b3268b5..cce164f 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -41,12 +41,23 @@ def std_3d_gaussian(): return g +@typedfixture('cartesiangaussian') +def cart_3d_gaussian(): + g = moldesign.orbitals.CartesianGaussian( + center=[random.random() for i in range(3)]*u.angstrom, + powers=[1, 3, 0], + exp=1.1/u.angstrom ** 2, + coeff=1.0) + 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) return g + @typedfixture('cartesiangaussian') def cart_10d_gaussian(): g = moldesign.orbitals.gaussians.CartesianGaussian( @@ -56,7 +67,6 @@ def cart_10d_gaussian(): return g - @pytest.mark.parametrize('objkey', registered_types['gaussian']) @pytest.mark.screening def test_gaussian_integral_and_dimensionality(objkey, request): @@ -70,6 +80,18 @@ def test_gaussian_integral_and_dimensionality(objkey, request): decimal=10) +def test_squared_gaussian_values(cart_3d_gaussian): + g = cart_3d_gaussian + randocoords = 6.0 * u.angstrom * (np.random.rand(200, 3) - 0.5) + + g2 = g * g + gvals = g(randocoords) + g2vals = g2(randocoords) + + helpers.assert_almost_equal(gvals**2, g2vals) + + + @pytest.mark.parametrize('objkey', registered_types['gaussian'] + registered_types['cartesiangaussian']) def test_gaussian_function_values(objkey, request): @@ -105,7 +127,7 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): _assert_almost_equal(vector_results, expected, decimal=8) -def test_gaussian_self_overlap_is_unity(): +def test_normalized_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) @@ -146,8 +168,9 @@ def test_cart_to_powers(): assert cart_to_powers('zx^3') == [3,0,1] -def test_numerical_vs_analytical_norm(std_3d_gaussian): - g = std_3d_gaussian +@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian']) +def test_numerical_vs_analytical_norm(key, request): + g = request.getfixturevalue(key) ananorm = g.norm width = np.sqrt(1/g.exp) diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 3c2f4a4..5b73b91 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -55,6 +55,7 @@ def norm(self): """ return np.sqrt(self.overlap(self)) + def __str__(self): return "%d-D %s with norm %s" % (self.ndims, self.__class__, self.norm) @@ -173,7 +174,7 @@ def angular_part(self, coords): axis = None disp = coords - self.center - return np.product(disp.magnitude**self.powers, axis=axis)*disp.units ** self.powers.sum() + return np.product(disp.magnitude**self.powers, axis=axis)*(disp.units ** self.powers.sum()) def __mul__(self, other): """ Returns product of two gaussian functions, which is also a gaussian @@ -224,7 +225,7 @@ def integral(self): elif p < 0: raise ValueError('Powers must be positive or 0') else: - integ *= _ODD_FACTORIAL[p-1]/(2 ** (p+1)) + integ *= _ODD_FACTORIAL[p-1]/((2.0 * self.exp) ** (p//2)) return self.coeff * integ From 72a8ebff9f69b6221d52d5e2cef24b94d15f18e4 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 14:19:06 -0700 Subject: [PATCH 12/33] Implement multiplication for arbitrary cartesian GBFs --- moldesign/_tests/test_gaussian_math.py | 76 ++++++++++++++++++-------- moldesign/orbitals/gaussians.py | 74 +++++++++++++++++++++---- 2 files changed, 116 insertions(+), 34 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index cce164f..0a42f50 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -3,6 +3,7 @@ import random +import itertools import numpy as np import numpy.testing as npt import pytest @@ -41,7 +42,7 @@ def std_3d_gaussian(): return g -@typedfixture('cartesiangaussian') +@typedfixture('gaussian') def cart_3d_gaussian(): g = moldesign.orbitals.CartesianGaussian( center=[random.random() for i in range(3)]*u.angstrom, @@ -52,18 +53,12 @@ def cart_3d_gaussian(): @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) - 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) +def spherical_3d_gaussian(): + g = moldesign.orbitals.SphericalGaussian( + center=[random.random() for i in range(3)]*u.angstrom, + l=3, m=-2, + exp=1.1/u.angstrom ** 2, + coeff=1.0) return g @@ -80,20 +75,53 @@ def test_gaussian_integral_and_dimensionality(objkey, request): decimal=10) -def test_squared_gaussian_values(cart_3d_gaussian): - g = cart_3d_gaussian - randocoords = 6.0 * u.angstrom * (np.random.rand(200, 3) - 0.5) - - g2 = g * g - gvals = g(randocoords) - g2vals = g2(randocoords) +def _make_rando_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()) + +test_powers = ((0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0), (0,2,0), (0,0,2), + (1,1,1), (2,2,2), (1,2,3), (0,1,4)) +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_gaussian_multiplication_amplitudes(p1, p2, withunits): + """ Tests that ``g1(x) * g2(x) == (g1 * g2)(x)`` + """ + gauss1 = _make_rando_gaussian(p1, withunits) + gauss2 = _make_rando_gaussian(p2, withunits) + + testcoords = 6.0*(np.random.rand(50, 3)-0.5) + if withunits: testcoords = testcoords * u.angstrom + + # test the product of the two AND the squares of each + #for g1, g2 in itertools.product((gauss1, gauss2), (gauss1, gauss2)): + g1, g2 = gauss1, gauss2 # temp + g1g2 = g1 * g2 + if isinstance(g1g2, moldesign.orbitals.AbstractFunction): + gvals = g1g2(testcoords) + assert g1g2.coeff != 0.0 + else: + gvals = sum(gg(testcoords) for gg in g1g2) + g1vals = g1(testcoords) + g2vals = g2(testcoords) - helpers.assert_almost_equal(gvals**2, g2vals) + prodvals = g1vals * g2vals + helpers.assert_almost_equal(prodvals, gvals) @pytest.mark.parametrize('objkey', - registered_types['gaussian'] + registered_types['cartesiangaussian']) + registered_types['gaussian']) def test_gaussian_function_values(objkey, request): g = request.getfixturevalue(objkey) @@ -107,7 +135,7 @@ def test_gaussian_function_values(objkey, request): @pytest.mark.parametrize('objkey', - registered_types['gaussian'] + registered_types['cartesiangaussian']) + registered_types['gaussian']) @pytest.mark.screening def test_vectorized_gaussian_function_evaluations(objkey, request): g = request.getfixturevalue(objkey) @@ -168,7 +196,7 @@ def test_cart_to_powers(): assert cart_to_powers('zx^3') == [3,0,1] -@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian']) +@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian', 'cart_3d_gaussian']) def test_numerical_vs_analytical_norm(key, request): g = request.getfixturevalue(key) ananorm = g.norm diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 5b73b91..ead962e 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -16,7 +16,11 @@ # 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 moldesign import units as u class AbstractFunction(object): @@ -55,7 +59,6 @@ def norm(self): """ return np.sqrt(self.overlap(self)) - def __str__(self): return "%d-D %s with norm %s" % (self.ndims, self.__class__, self.norm) @@ -95,7 +98,7 @@ class CartesianGaussian(AbstractFunction): def __init__(self, center, exp, powers, coeff=None): assert len(powers) == len(center), "Inconsistent dimensionality - number of cartesian " \ "powers must match dimensionality of centroid vector" - self.center = center + self.center = u.array(center) self.exp = exp self.powers = np.array(powers) if coeff is None: @@ -174,28 +177,79 @@ def angular_part(self, coords): axis = None disp = coords - self.center - return np.product(disp.magnitude**self.powers, axis=axis)*(disp.units ** self.powers.sum()) + 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 gaussian functions, which is also a gaussian + """ 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 wavepacket Returns: - CartesianGaussian: product gaussian + CartesianGaussian or List[CartesianGaussian]: product gaussian(s) + + TODO: + - Should return a composite gaussian object with the same interface as regular gaussians, + rather than either an object OR a list """ - if (self.center != other.center).any(): - raise NotImplementedError() + prefactor, newcenter, newexp = self._mul_gaussians(other) + new_prefactor = prefactor * self.coeff * other.coeff + + coeffs = [{}, {}, {}] + for idim in range(3): + r_self = self.center[idim] + r_other = other.center[idim] + r_new = newcenter[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 coeffs[idim]: + coeffs[idim][newpower] = powercoeff + else: + coeffs[idim][newpower] += powercoeff + + new_gaussians = [] + for (l, c0), (m, c1), (n,c2) in itertools.product(*[x.items() for x in coeffs]): + final_coeff = c0 * c1 * c2 * new_prefactor + new_gaussians.append(CartesianGaussian(newcenter, newexp, powers=[l,m,n], + coeff=final_coeff)) + + if len(new_gaussians) == 0: # no non-zero components, just return a zeroed gaussian + return CartesianGaussian(newcenter, newexp, [0,0,0], coeff=0.0) + elif len(new_gaussians) == 1: + return new_gaussians[0] + else: + return new_gaussians + def _mul_gaussians(self, other): + """ Returns the new centroid and width of THE EXPONENTIAL PART ONLY of two gaussians. + """ # convert widths to prefactor form a = self.exp b = other.exp exp = a + b center = (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) + prefactor = 1.0 + for i in range(3): + prefactor *= np.exp(-(self.exp*self.center[i]**2 + other.exp * other.center[i]**2) + + exp * center[i]**2) + + return prefactor, center, exp @property def integral(self): From f48a8e9585848d028cc8407898eea8df179b54f8 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 15:51:53 -0700 Subject: [PATCH 13/33] Clean up gaussian math and tests --- moldesign/_tests/test_gaussian_math.py | 87 ++++++++++++++++++-------- moldesign/orbitals/gaussians.py | 39 ++++++++---- 2 files changed, 89 insertions(+), 37 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 0a42f50..20fae85 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -10,6 +10,7 @@ import moldesign from moldesign import units as u +from moldesign.mathutils import spherical_harmonics as harmonics from . import helpers @@ -42,7 +43,7 @@ def std_3d_gaussian(): return g -@typedfixture('gaussian') +@typedfixture('cart-gaussian') def cart_3d_gaussian(): g = moldesign.orbitals.CartesianGaussian( center=[random.random() for i in range(3)]*u.angstrom, @@ -52,7 +53,7 @@ def cart_3d_gaussian(): return g -@typedfixture('gaussian') +@typedfixture('spherical-gaussian') def spherical_3d_gaussian(): g = moldesign.orbitals.SphericalGaussian( center=[random.random() for i in range(3)]*u.angstrom, @@ -85,27 +86,26 @@ def _make_rando_gaussian(powers, withunits=True): powers=powers, coeff=random.random()) -test_powers = ((0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0), (0,2,0), (0,0,2), - (1,1,1), (2,2,2), (1,2,3), (0,1,4)) +# 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_gaussian_multiplication_amplitudes(p1, p2, withunits): +def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): """ Tests that ``g1(x) * g2(x) == (g1 * g2)(x)`` """ - gauss1 = _make_rando_gaussian(p1, withunits) - gauss2 = _make_rando_gaussian(p2, withunits) + g1 = _make_rando_gaussian(p1, withunits) + g2 = _make_rando_gaussian(p2, withunits) testcoords = 6.0*(np.random.rand(50, 3)-0.5) - if withunits: testcoords = testcoords * u.angstrom + if withunits: + testcoords = testcoords * u.angstrom - # test the product of the two AND the squares of each - #for g1, g2 in itertools.product((gauss1, gauss2), (gauss1, gauss2)): - g1, g2 = gauss1, gauss2 # temp g1g2 = g1 * g2 if isinstance(g1g2, moldesign.orbitals.AbstractFunction): gvals = g1g2(testcoords) @@ -121,7 +121,9 @@ def test_gaussian_multiplication_amplitudes(p1, p2, withunits): @pytest.mark.parametrize('objkey', - registered_types['gaussian']) + registered_types['gaussian'] + + registered_types['cart-gaussian'] + + registered_types['spherical-gaussian']) def test_gaussian_function_values(objkey, request): g = request.getfixturevalue(objkey) @@ -155,22 +157,57 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): _assert_almost_equal(vector_results, expected, decimal=8) -def test_normalized_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) +@pytest.mark.parametrize('objkey', + registered_types['gaussian'] + + registered_types['cart-gaussian']) +def test_gaussian_str_and_repr_works(objkey, request): + g1 = request.getfixturevalue(objkey) + str(g1) + repr(g1) + + +@pytest.mark.parametrize('objkey', + registered_types['gaussian'] + + registered_types['cart-gaussian'] + + registered_types['spherical-gaussian']) +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 npt.assert_almost_equal(-1.0, g1.overlap(g2, normalized=True)) + g1.coeff = 10.0 + npt.assert_almost_equal(1.0, + g1.overlap(g2, normalized=True)) + + +@pytest.mark.parametrize('objkey', + registered_types['gaussian'] + + registered_types['cart-gaussian'] + + registered_types['spherical-gaussian']) +def test_normalization(objkey, request): + g1 = request.getfixturevalue(objkey) + oldnorm = g1.norm + + g1.coeff = (random.random() - 0.5) * 428.23 + assert g1.norm != oldnorm + + 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 + r2 = np.sum(r**2) + fv = g.coeff * np.exp(-g.exp * r2) + if isinstance(g, moldesign.orbitals.SphericalGaussian): + fv *= r2**(g.l/2.0) * harmonics.Y(g.l, g.m)(coord - g.center) + else: # assume cartesian + for r, r0, pow in zip(coord, g.center, g.powers): + if pow != 0: + fv *= (r-r0)**pow return fv @@ -189,14 +226,14 @@ 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] -@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian', 'cart_3d_gaussian']) +@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian', 'spherical_3d_gaussian']) def test_numerical_vs_analytical_norm(key, request): g = request.getfixturevalue(key) ananorm = g.norm diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index ead962e..9a0ca4d 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -17,6 +17,7 @@ # See the License for the specific language governing permissions and # limitations under the License. import itertools +import copy import numpy as np from scipy.special import binom @@ -26,10 +27,14 @@ class AbstractFunction(object): """ Abstract base class for basis functions """ + + def copy(self): + return copy.deepcopy(self) + def normalize(self): """ Give this function unit norm by adjusting its coefficient """ - self.coeff /= np.sqrt(self.norm) + self.coeff /= self.norm def overlap(self, other, normalized=False): r""" Overlap of this function with another: @@ -47,7 +52,7 @@ def overlap(self, other, normalized=False): newfn = self * other integral = newfn.integral if normalized: - integral /= np.sqrt(self.norm*other.norm) + integral /= (self.norm*other.norm) return integral @property @@ -201,11 +206,19 @@ def __mul__(self, other): - Should return a composite gaussian object with the same interface as regular gaussians, rather than either an object OR a list """ + + if (self.center == other.center).all(): # this case is much easier than the general one + cpy = self.copy() + cpy.exp = self.exp + other.exp + cpy.powers = self.powers + other.powers + cpy.coeff = self.coeff * other.coeff + return cpy + prefactor, newcenter, newexp = self._mul_gaussians(other) new_prefactor = prefactor * self.coeff * other.coeff - coeffs = [{}, {}, {}] - for idim in range(3): + 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[idim] @@ -218,19 +231,21 @@ def __mul__(self, other): if powercoeff == 0.0: continue newpower = m+k - if newpower not in coeffs[idim]: - coeffs[idim][newpower] = powercoeff + if newpower not in partial_coeffs[idim]: + partial_coeffs[idim][newpower] = powercoeff else: - coeffs[idim][newpower] += powercoeff + partial_coeffs[idim][newpower] += powercoeff new_gaussians = [] - for (l, c0), (m, c1), (n,c2) in itertools.product(*[x.items() for x in coeffs]): - final_coeff = c0 * c1 * c2 * new_prefactor - new_gaussians.append(CartesianGaussian(newcenter, newexp, powers=[l,m,n], + for powers in itertools.product(*[x.keys() for x in partial_coeffs]): + final_coeff = 1.0 * new_prefactor + for idim, p in enumerate(powers): + final_coeff *= partial_coeffs[idim][p] + new_gaussians.append(CartesianGaussian(newcenter, newexp, powers=powers, coeff=final_coeff)) if len(new_gaussians) == 0: # no non-zero components, just return a zeroed gaussian - return CartesianGaussian(newcenter, newexp, [0,0,0], coeff=0.0) + return CartesianGaussian(newcenter, newexp, self.powers + other.powers, coeff=0.0) elif len(new_gaussians) == 1: return new_gaussians[0] else: @@ -245,7 +260,7 @@ def _mul_gaussians(self, other): exp = a + b center = (a*self.center + b*other.center)/(a+b) prefactor = 1.0 - for i in range(3): + for i in range(self.ndim): prefactor *= np.exp(-(self.exp*self.center[i]**2 + other.exp * other.center[i]**2) + exp * center[i]**2) From d0a9d8e5918d9211f6f6db4864f49856823e5f48 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 19:17:27 -0700 Subject: [PATCH 14/33] Refactor gaussian subpackage and get spherical overlaps working --- moldesign/orbitals/__init__.py | 2 + moldesign/orbitals/cartesian.py | 227 ++++++++++++++++++++++ moldesign/orbitals/gaussians.py | 335 ++++---------------------------- moldesign/orbitals/spherical.py | 146 ++++++++++++++ 4 files changed, 413 insertions(+), 297 deletions(-) create mode 100644 moldesign/orbitals/cartesian.py create mode 100644 moldesign/orbitals/spherical.py diff --git a/moldesign/orbitals/__init__.py b/moldesign/orbitals/__init__.py index f589c6e..65e8d7f 100644 --- a/moldesign/orbitals/__init__.py +++ b/moldesign/orbitals/__init__.py @@ -4,6 +4,8 @@ def toplevel(o): __all__ = [] from .gaussians import * +from .cartesian import * +from .spherical import * from .orbitals import * from .atomic_basis_fn import * from .basis import * diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py new file mode 100644 index 0000000..fc7363d --- /dev/null +++ b/moldesign/orbitals/cartesian.py @@ -0,0 +1,227 @@ +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 Gaussian + + +class CartesianGaussian(Gaussian): + 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^{-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} + + 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`` + + References: + Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. + + 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. + """ + def __init__(self, center, exp, powers, coeff=None): + assert len(powers) == len(center), "Inconsistent dimensionality - number of cartesian " \ + "powers must match dimensionality of centroid vector" + + super().__init__(center, exp, coeff=1.0) # temporarily set coefficient + self.powers = np.array(powers) + self.shell = sum(self.powers) + if coeff is None: + self.normalize() + else: + self.coeff = coeff + + def __repr__(self): + return ("<{ndim}-D cartesian gaussian (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) + + @property + def angular(self): + """ Angular momentum of this function (sum of cartesian powers) + """ + return self.powers.sum() + + def __call__(self, coords, _include_angular=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: + coords (Vector[length, len=3] or Matrix[length, shape=(*,3)]): Coordinates or + list of 3D coordinates + _include_cartesian (bool): include the contribution from the cartesian components + (for computational efficiency, this can sometimes omited now and included later) + + Examples: + >>> g = CartesianGaussian([0,0,0]*u.angstrom, exp=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 = super().__call__(coords) + + if self.shell > 0 and _include_angular: + result *= self.angular_part(coords) + return result + + 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 wavepacket + + Returns: + CartesianGaussian or List[CartesianGaussian]: product gaussian(s) + + TODO: + - Should return a composite gaussian object with the same interface as regular gaussians, + rather than either an object OR a list + """ + + if (self.center == other.center).all(): # this case is much easier than the general one + cpy = self.copy() + cpy.exp = self.exp + other.exp + cpy.powers = self.powers + other.powers + cpy.coeff = self.coeff * other.coeff + return cpy + + newcenter = super().__mul__(other) + + 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 = 1.0 * newcenter.coeff + for idim, p in enumerate(powers): + final_coeff *= partial_coeffs[idim][p] + new_gaussians.append(CartesianGaussian(newcenter.center, newcenter.exp, + powers=powers, coeff=final_coeff)) + + if len(new_gaussians) == 0: # no non-zero components, just return a zeroed gaussian + return CartesianGaussian(newcenter, newcenter.exp, + self.powers + other.powers, coeff=0.0) + elif len(new_gaussians) == 1: + return new_gaussians[0] + else: + return 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 = super().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.exp) ** (p//2)) + return integ diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 9a0ca4d..66a760b 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -16,12 +16,10 @@ # 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 copy - import numpy as np -from scipy.special import binom -from moldesign import units as u + +from .. import units as u class AbstractFunction(object): @@ -57,7 +55,7 @@ def overlap(self, other, normalized=False): @property def norm(self): - r""" The L2-Norm of this gaussian: + 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} @@ -68,95 +66,40 @@ 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(AbstractFunction): + 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} - - 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`` - - References: - Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. - - 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) + G(\mathbf r) = C e^{-a\left| \mathbf r - \mathbf r_0 \right|^2} - 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. + where *C* is ``self.coeff``, *a* is ``self.exp``, and :math:`\mathbf r_0` is ``self.center``. """ - def __init__(self, center, exp, powers, coeff=None): - assert len(powers) == len(center), "Inconsistent dimensionality - number of cartesian " \ - "powers must match dimensionality of centroid vector" + def __init__(self, center, exp, coeff=None): self.center = u.array(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.shell = sum(self.powers) - - 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) - - @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 __call__(self, coords, _include_angular=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: - coords (Vector[length, len=3] or Matrix[length, shape=(*,3)]): Coordinates or - list of 3D coordinates - _include_cartesian (bool): include the contribution from the cartesian components - (for computational efficiency, this can sometimes omited now and included later) - - Examples: - >>> g = CartesianGaussian([0,0,0]*u.angstrom, exp=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]) + _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]): 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 @@ -168,245 +111,43 @@ def __call__(self, coords, _include_angular=True): r2 = prd.sum(axis=axis) result = self.coeff * np.exp(-self.exp * r2) - if self.shell > 0 and _include_angular: - result *= self.angular_part(coords) - return result - - 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() + if _getvals: + return result, disp, r2 # extra quantities so that they don't need to be recomputed 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 wavepacket + return result - Returns: - CartesianGaussian or List[CartesianGaussian]: product gaussian(s) - - TODO: - - Should return a composite gaussian object with the same interface as regular gaussians, - rather than either an object OR a list - """ + def __repr__(self): + return ("<{ndim}-D gaussian (coeff: {coeff:4.2f}, " + "exponent: {exp:4.2f}, " + "center: {center}>").format( + ndim=self.ndim, + center=self.center, exp=self.exp, + coeff=self.coeff) - if (self.center == other.center).all(): # this case is much easier than the general one - cpy = self.copy() - cpy.exp = self.exp + other.exp - cpy.powers = self.powers + other.powers - cpy.coeff = self.coeff * other.coeff - return cpy - - prefactor, newcenter, newexp = self._mul_gaussians(other) - new_prefactor = prefactor * self.coeff * other.coeff - - 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[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 = 1.0 * new_prefactor - for idim, p in enumerate(powers): - final_coeff *= partial_coeffs[idim][p] - new_gaussians.append(CartesianGaussian(newcenter, newexp, powers=powers, - coeff=final_coeff)) - - if len(new_gaussians) == 0: # no non-zero components, just return a zeroed gaussian - return CartesianGaussian(newcenter, newexp, self.powers + other.powers, coeff=0.0) - elif len(new_gaussians) == 1: - return new_gaussians[0] - else: - return new_gaussians + @property + def ndim(self): + return len(self.center) + num_dimensions = ndims = ndim - def _mul_gaussians(self, other): - """ Returns the new centroid and width of THE EXPONENTIAL PART ONLY of two gaussians. + def __mul__(self, other): + """ Returns a new gaussian centroid. """ - # convert widths to prefactor form a = self.exp b = other.exp exp = a + b center = (a*self.center + b*other.center)/(a+b) - prefactor = 1.0 + prefactor = self.coeff * other.coeff for i in range(self.ndim): - prefactor *= np.exp(-(self.exp*self.center[i]**2 + other.exp * other.center[i]**2) + + prefactor *= np.exp(-(a*self.center[i]**2 + b*other.center[i]**2) + exp * center[i]**2) - return prefactor, center, exp + return Gaussian(center, exp, coeff=prefactor) @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 = (np.pi/self.exp)**(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 *= _ODD_FACTORIAL[p-1]/((2.0 * self.exp) ** (p//2)) - 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 real spherical gaussian function: - - .. math:: - G_{nlm}(\mathbf r) = C Y^l_m(\mathbf r - \mathbf r_0) r^l 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``, - *(l,m)* are given by ``(self.l, self.m)``, and :math:`Y^l_m(\mathbf{r})` are the _real_ - spherical harmonics. - - 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 - """ - - ndims = 3 - - def __init__(self, center, exp, l, m, coeff=None): - from moldesign.mathutils import spherical_harmonics - - self.center = center - assert len(self.center) == 3 - self.exp = exp - self.l = l - self.m = m - - if coeff is None: - self.coeff = 1.0 - self.normalize() - else: - self.coeff = coeff - - self._spherical_harmonic = spherical_harmonics.Y(l, m) - - def __repr__(self): - return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " - "exponent: {exp:4.2f}, " - "(l,m) = {qnums}").format( - center=self.center, exp=self.exp, coeff=self.coeff, - qnums=(self.l, self.m)) - - def __call__(self, coords, _include_angular=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_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]): 3D Coordinates or list of 3D coordinates - - Returns: - Scalar: function value(s) at the passed coordinates + r"""Integral of this this gaussian over all N-dimensional space """ - if len(coords.shape) > 1: - axis = 1 - else: - axis = None - - disp = coords - self.center - prd = disp*disp - r2 = prd.sum(axis=axis) - - result = self.coeff * np.exp(-self.exp * r2) - if _include_angular: - result *= r2**(self.l/2.0) * self._spherical_harmonic(disp) - return result - - 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) - - -# 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) + return self.coeff * (np.pi/self.exp)**(self.ndim/2.0) def cart_to_powers(s): diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py new file mode 100644 index 0000000..96ddb39 --- /dev/null +++ b/moldesign/orbitals/spherical.py @@ -0,0 +1,146 @@ +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 . import Gaussian + + +class SphericalGaussian(Gaussian): + 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^{-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``, + *(l,m)* are given by ``(self.l, self.m)``, and :math:`Y^l_m(\mathbf{r})` are the _real_ + spherical harmonics. + + 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 + """ + + ndims = 3 + + def __init__(self, center, exp, l, m, coeff=None): + from moldesign.mathutils import spherical_harmonics + + super().__init__(center, exp, coeff=1.0) # temporarily set coefficient + + self.center = center + assert len(self.center) == 3 + self.exp = exp + self.l = l + self.m = m + + if coeff is None: + self.coeff = 1.0 + self.normalize() + else: + self.coeff = coeff + + self._spherical_harmonic = spherical_harmonics.Y(l, m) + + def __repr__(self): + return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " + "exponent: {exp:4.2f}, " + "(l,m) = {qnums}").format( + center=self.center, exp=self.exp, coeff=self.coeff, + qnums=(self.l, self.m)) + + def __mul__(self, other): + raise NotImplementedError("Cannot multiply spherical gaussian functions") + + def __call__(self, coords, _include_angular=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_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]): 3D Coordinates or list of 3D coordinates + + Returns: + Scalar or Vector: function value(s) at the passed coordinates + """ + result, disp, r2 = super().__call__(coords, _getvals=True) + + if _include_angular: + result *= r2**(self.l/2.0) * self._spherical_harmonic(disp) + return result + + def overlap(self, other, normalized=False): + r""" Overlap of this function with another spherical basis function: + + .. 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 + """ + from ..data import ODD_FACTORIAL + + if (self.center != other.center).any(): + raise NotImplementedError("Two-center spherical harmonic overlaps not yet supported") + + elif self.l != other.l or self.m != other.m: + return 0.0 + + else: + newexp = self.exp + other.exp + 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) \ No newline at end of file From 15797d3b49191f6179a3554f7e6cfba0a047037a Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 19:34:51 -0700 Subject: [PATCH 15/33] Update tests for spherical basis --- moldesign/_tests/test_gaussian_math.py | 101 +++++++++++++++++---- moldesign/_tests/test_mathutils.py | 38 ++++++++ moldesign/data/data.py | 8 ++ moldesign/mathutils/spherical_harmonics.py | 25 ++--- 4 files changed, 143 insertions(+), 29 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 20fae85..696cd5d 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -76,7 +76,7 @@ def test_gaussian_integral_and_dimensionality(objkey, request): decimal=10) -def _make_rando_gaussian(powers, withunits=True): +def _make_rando_cart_gaussian(powers, withunits=True): if withunits: length = u.angstrom else: @@ -86,6 +86,17 @@ def _make_rando_gaussian(powers, withunits=True): 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()) + # 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])) @@ -99,14 +110,13 @@ def _make_rando_gaussian(powers, withunits=True): def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): """ Tests that ``g1(x) * g2(x) == (g1 * g2)(x)`` """ - g1 = _make_rando_gaussian(p1, withunits) - g2 = _make_rando_gaussian(p2, withunits) + g1 = _make_rando_cart_gaussian(p1, withunits) + g2 = _make_rando_cart_gaussian(p2, withunits) testcoords = 6.0*(np.random.rand(50, 3)-0.5) if withunits: - testcoords = testcoords * u.angstrom - - g1g2 = g1 * g2 + testcoords = testcoords*u.angstrom + g1g2 = g1*g2 if isinstance(g1g2, moldesign.orbitals.AbstractFunction): gvals = g1g2(testcoords) assert g1g2.coeff != 0.0 @@ -114,9 +124,7 @@ def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): gvals = sum(gg(testcoords) for gg in g1g2) g1vals = g1(testcoords) g2vals = g2(testcoords) - - prodvals = g1vals * g2vals - + prodvals = g1vals*g2vals helpers.assert_almost_equal(prodvals, gvals) @@ -138,7 +146,6 @@ def test_gaussian_function_values(objkey, request): @pytest.mark.parametrize('objkey', registered_types['gaussian']) -@pytest.mark.screening def test_vectorized_gaussian_function_evaluations(objkey, request): g = request.getfixturevalue(objkey) @@ -159,7 +166,8 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): @pytest.mark.parametrize('objkey', registered_types['gaussian'] + - registered_types['cart-gaussian']) + registered_types['cart-gaussian'] + + registered_types['spherical-gaussian']) def test_gaussian_str_and_repr_works(objkey, request): g1 = request.getfixturevalue(objkey) str(g1) @@ -175,12 +183,12 @@ def test_normalized_gaussian_self_overlap_is_unity(objkey, request): g2 = g1.copy() g1.coeff = -10.0 g2.coeff = 12341.1832 - npt.assert_almost_equal(-1.0, - g1.overlap(g2, normalized=True)) + olap = g1.overlap(g2, normalized=True) + assert abs(-1.0 - olap) < 1e-12 g1.coeff = 10.0 - npt.assert_almost_equal(1.0, - g1.overlap(g2, normalized=True)) + olap = g1.overlap(g2, normalized=True) + assert abs(1.0 - olap) < 1e-12 @pytest.mark.parametrize('objkey', @@ -192,7 +200,10 @@ def test_normalization(objkey, request): oldnorm = g1.norm g1.coeff = (random.random() - 0.5) * 428.23 - assert g1.norm != oldnorm + 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 @@ -200,11 +211,14 @@ def test_normalization(objkey, request): def _gfuncval(g, coord): r = g.center - coord - r2 = np.sum(r**2) + if len(coord.shape) > 1: + r2 = np.sum(r**2, axis=1) + else: + r2 = np.sum(r**2) fv = g.coeff * np.exp(-g.exp * r2) if isinstance(g, moldesign.orbitals.SphericalGaussian): fv *= r2**(g.l/2.0) * harmonics.Y(g.l, g.m)(coord - g.center) - else: # assume cartesian + 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 @@ -252,3 +266,54 @@ def test_numerical_vs_analytical_norm(key, request): helpers.assert_almost_equal(ananorm, 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 + testcoords = 6.0*(np.random.rand(50, 3)-0.5) * u.angstrom + + 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 + + barevals = g_bare(testcoords) + cartvals = g_cart(testcoords) + spherevals = g_sphr(testcoords) + + helpers.assert_almost_equal(barevals, cartvals) + helpers.assert_almost_equal(barevals, spherevals) + + helpers.assert_almost_equal(g_bare.norm, g_cart.norm) + helpers.assert_almost_equal(g_sphr.norm, g_cart.norm) + + +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 + testcoords = 6.0*(np.random.rand(50, 3)-0.5) * u.angstrom + + g_cart = moldesign.orbitals.CartesianGaussian(center, exp, powers) + g_sphr = moldesign.orbitals.SphericalGaussian(center, exp, *lm) + + helpers.assert_almost_equal(g_cart(testcoords), + g_sphr(testcoords)) + + helpers.assert_almost_equal(g_sphr.norm, g_cart.norm) diff --git a/moldesign/_tests/test_mathutils.py b/moldesign/_tests/test_mathutils.py index 230a2df..84b5ddd 100644 --- a/moldesign/_tests/test_mathutils.py +++ b/moldesign/_tests/test_mathutils.py @@ -1,3 +1,4 @@ +import itertools import pytest import numpy as np @@ -135,6 +136,43 @@ def test_cart_and_spherical_equivalence(shell): 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 + + RANDARRAY = np.array([[ 0.03047527, -4.49249374, 5.83154878], [-7.38885486, 8.30076569, -7.97671261], [ 6.79409582, -4.07664079, 8.45983794], 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/mathutils/spherical_harmonics.py b/moldesign/mathutils/spherical_harmonics.py index cf92b09..047e43d 100644 --- a/moldesign/mathutils/spherical_harmonics.py +++ b/moldesign/mathutils/spherical_harmonics.py @@ -20,8 +20,6 @@ __all__ = ['SPHERE_TO_CART'] -SQ2INV = np.sqrt(1/2.0) - def cart_to_polar_angles(coords): if len(coords.shape) == 2: @@ -38,7 +36,12 @@ def cart_to_polar_angles(coords): class Y(object): - """ A *real* spherical harmonic + r""" A *real*-valued spherical harmonic function + + These functions are orthonormalized over the sphere such that + + .. math:: + \int d\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 @@ -57,17 +60,17 @@ 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 - else: - vplus = sph_harm(abs(self.m), self.l, phi, theta) - vminus = sph_harm(-abs(self.m), self.l, phi, theta) - if self.m < 0: - return -(SQ2INV * (self._posneg * vplus + vminus)).imag - else: - return (SQ2INV * (self._posneg * vplus + vminus)).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): From 44a7a785a2a7412587f8873c68f72f358de09ae4 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 21:24:49 -0700 Subject: [PATCH 16/33] Factor generalized contraction class out of AtomicBasisFunction --- moldesign/_tests/test_gaussian_math.py | 7 +---- moldesign/orbitals/__init__.py | 3 +- moldesign/orbitals/atomic_basis_fn.py | 38 ++------------------------ moldesign/orbitals/cartesian.py | 10 ++----- 4 files changed, 9 insertions(+), 49 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 696cd5d..7294146 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -5,7 +5,6 @@ import itertools import numpy as np -import numpy.testing as npt import pytest import moldesign @@ -117,11 +116,7 @@ def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): if withunits: testcoords = testcoords*u.angstrom g1g2 = g1*g2 - if isinstance(g1g2, moldesign.orbitals.AbstractFunction): - gvals = g1g2(testcoords) - assert g1g2.coeff != 0.0 - else: - gvals = sum(gg(testcoords) for gg in g1g2) + gvals = g1g2(testcoords) g1vals = g1(testcoords) g2vals = g2(testcoords) prodvals = g1vals*g2vals diff --git a/moldesign/orbitals/__init__.py b/moldesign/orbitals/__init__.py index 65e8d7f..dec4e13 100644 --- a/moldesign/orbitals/__init__.py +++ b/moldesign/orbitals/__init__.py @@ -4,9 +4,10 @@ def toplevel(o): __all__ = [] from .gaussians import * +from .orbitals import * +from .contraction import * from .cartesian import * from .spherical import * -from .orbitals 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 index d388ac6..ae31b66 100644 --- a/moldesign/orbitals/atomic_basis_fn.py +++ b/moldesign/orbitals/atomic_basis_fn.py @@ -16,12 +16,10 @@ # 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 BasisContraction, SHELLS, SPHERICALNAMES -import numpy as np -from .orbitals import Orbital, SHELLS, SPHERICALNAMES - -class AtomicBasisFunction(Orbital): +class AtomicBasisFunction(BasisContraction): """ Stores an atomic basis function. Note: @@ -39,12 +37,12 @@ class AtomicBasisFunction(Orbital): 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 - self.primitives = primitives if cart is not None: assert self.m is None, 'Both cartesian and spherical components passed!' assert len(cart) == self.l, \ @@ -62,12 +60,6 @@ def __init__(self, atom, n=None, l=None, m=None, cart=None, primitives=None): self.occupation = None self.wfn = None - def __call__(self, coords): - outvals = np.zeros(len(coords)) - for primitive in self.primitives: - outvals += primitive(coords) - return outvals - @property def atom(self): """ moldesign.Atom: the atom this basis function belongs to @@ -79,30 +71,6 @@ def atom(self): 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:`\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 *= prefactor - @property def orbtype(self): """ A string describing the orbital's angular momentum state. diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index fc7363d..a9d264e 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -22,7 +22,7 @@ from scipy.special import binom from .. import units as u -from . import Gaussian +from . import Gaussian, BasisContraction class CartesianGaussian(Gaussian): @@ -142,11 +142,7 @@ def __mul__(self, other): other (CartesianGaussian): other gaussian wavepacket Returns: - CartesianGaussian or List[CartesianGaussian]: product gaussian(s) - - TODO: - - Should return a composite gaussian object with the same interface as regular gaussians, - rather than either an object OR a list + CartesianGaussian or BasisContraction: product functions """ if (self.center == other.center).all(): # this case is much easier than the general one @@ -191,7 +187,7 @@ def __mul__(self, other): elif len(new_gaussians) == 1: return new_gaussians[0] else: - return new_gaussians + return BasisContraction(new_gaussians) @property def integral(self): From a17886e6953abf01d7ccc9cd2a18cd4d910e179f Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 21:25:56 -0700 Subject: [PATCH 17/33] Add the actual contractoin source code ... --- moldesign/orbitals/contraction.py | 73 +++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 moldesign/orbitals/contraction.py diff --git a/moldesign/orbitals/contraction.py b/moldesign/orbitals/contraction.py new file mode 100644 index 0000000..bba75d9 --- /dev/null +++ b/moldesign/orbitals/contraction.py @@ -0,0 +1,73 @@ +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 .orbitals import Orbital + + +class BasisContraction(Orbital): + """ Stores a linear combination of primitive functions. + + Args: + primitives (List[PrimitiveBase]): List of primitives, if available + """ + def __init__(self, primitives): + self.primitives = primitives + + def __call__(self, coords): + outvals = np.zeros(len(coords)) + for primitive in self.primitives: + outvals += primitive(coords) + return outvals + + @property + def num_primitives(self): + return len(self.primitives) + + @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 *= 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(other) for p1 in self.primitives) + return olap + + def __str__(self): + return '%s with %d primitives' % (self.__class__.__name__, self.num_primitives) + + def __repr__(self): + return '<%s>' % self \ No newline at end of file From 02aa93e5e51484600361a2a80dfb2f7dd240ce86 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 10 Jul 2017 23:00:17 -0700 Subject: [PATCH 18/33] Explicitly factor radial and angular parts --- moldesign/orbitals/cartesian.py | 64 ++++++++++++++-------------- moldesign/orbitals/gaussians.py | 55 ++++++++++++++---------- moldesign/orbitals/spherical.py | 74 ++++++++++++++++----------------- 3 files changed, 102 insertions(+), 91 deletions(-) diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index a9d264e..14dff6a 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -22,10 +22,11 @@ from scipy.special import binom from .. import units as u -from . import Gaussian, BasisContraction +from .. import utils +from . import AbstractFunction, BasisContraction, Gaussian -class CartesianGaussian(Gaussian): +class CartesianGaussian(AbstractFunction): r""" Stores an N-dimensional gaussian function of the form: .. math:: @@ -35,9 +36,9 @@ class CartesianGaussian(Gaussian): 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(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``, *a* is ``self.exp``, *r0* is ``self.center``, and + where *C* is ``self.coeff``, *alpha* is ``self.alpha``, *r0* is ``self.center``, and :math:`p_1, p_2, ...` are given in the array ``self.powers`` References: @@ -47,7 +48,7 @@ class CartesianGaussian(Gaussian): 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 + alpha (Scalar[1/length**2]): gaussian width parameter coeff (Scalar): multiplicative coefficient (if None, gaussian will be automatically normalized) @@ -57,26 +58,27 @@ class CartesianGaussian(Gaussian): ``center`` and ``powers``, it's 1-D. If length-3 vectors are passed for ``center`` and ``powers``, it's 3D. """ - def __init__(self, center, exp, powers, coeff=None): + 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" - super().__init__(center, exp, coeff=1.0) # temporarily set coefficient self.powers = np.array(powers) self.shell = sum(self.powers) - if coeff is None: - self.normalize() - else: - self.coeff = coeff + 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 (coeff: {coeff:4.2f}, " + return ("<{ndim}-D cartesian gaussian (norm: {norm:4.2f}, " "cartesian powers: {powers}, " - "exponent: {exp:4.2f}, " + "width: {exp:4.2f}, " "center: {center}>").format( ndim=self.ndim, - center=self.center, exp=self.exp, - powers=tuple(self.powers), coeff=self.coeff) + center=self.center, exp=self.alpha, + powers=tuple(self.powers), norm=self.norm) @property def angular(self): @@ -84,7 +86,7 @@ def angular(self): """ return self.powers.sum() - def __call__(self, coords, _include_angular=True): + 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 @@ -93,11 +95,9 @@ def __call__(self, coords, _include_angular=True): Args: coords (Vector[length, len=3] or Matrix[length, shape=(*,3)]): Coordinates or list of 3D coordinates - _include_cartesian (bool): include the contribution from the cartesian components - (for computational efficiency, this can sometimes omited now and included later) Examples: - >>> g = CartesianGaussian([0,0,0]*u.angstrom, exp=1.0/u.angstrom**2, powers=(0,0,0)) + >>> 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 @@ -108,11 +108,11 @@ def __call__(self, coords, _include_angular=True): Returns: Scalar or Vector: function value(s) at the passed coordinates """ - result = super().__call__(coords) + result = self.radial_part(coords) - if self.shell > 0 and _include_angular: + if self.shell > 0: result *= self.angular_part(coords) - return result + return result * self.coeff def angular_part(self, coords): if self.shell == 0: @@ -147,12 +147,12 @@ def __mul__(self, other): if (self.center == other.center).all(): # this case is much easier than the general one cpy = self.copy() - cpy.exp = self.exp + other.exp + cpy.radial_part.alpha = self.alpha + other.alpha cpy.powers = self.powers + other.powers cpy.coeff = self.coeff * other.coeff return cpy - newcenter = super().__mul__(other) + newcenter = self.radial_part * other.radial_part partial_coeffs = [{} for idim in range(self.ndim)] for idim in range(self.ndim): @@ -175,15 +175,17 @@ def __mul__(self, other): new_gaussians = [] for powers in itertools.product(*[x.keys() for x in partial_coeffs]): - final_coeff = 1.0 * newcenter.coeff + 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.exp, - powers=powers, coeff=final_coeff)) + 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.exp, - self.powers + other.powers, coeff=0.0) + return CartesianGaussian(newcenter, newcenter.alpha, + self.powers + other.powers, coeff=0.0, + normalized=False) elif len(new_gaussians) == 1: return new_gaussians[0] else: @@ -210,7 +212,7 @@ def integral(self): """ from ..data import ODD_FACTORIAL - integ = super().integral + integ = self.coeff * self.radial_part.integral for p in self.powers: if p == 0: # no contribution continue @@ -219,5 +221,5 @@ def integral(self): elif p < 0: raise ValueError('Powers must be positive or 0') else: - integ *= ODD_FACTORIAL[p-1]/((2.0 * self.exp) ** (p//2)) + 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 66a760b..ab0425d 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -25,6 +25,13 @@ class AbstractFunction(object): """ Abstract base class for basis functions """ + def __init__(self, coeff=None, normalized=False): + if normalized or coeff is None: + self.coeff = 1.0 # dummy value overwritten by self.normalize() + self.normalize() + + if coeff is not None: + self.coeff = coeff def copy(self): return copy.deepcopy(self) @@ -62,28 +69,28 @@ def norm(self): """ return np.sqrt(self.overlap(self)) - def __str__(self): - return "%d-D %s with norm %s" % (self.ndims, self.__class__, self.norm) - class Gaussian(AbstractFunction): r""" N-dimensional gaussian function. The function is given by: .. math:: - G(\mathbf r) = C e^{-a\left| \mathbf r - \mathbf r_0 \right|^2} + G(\mathbf r) = C e^{-\alpha\left| \mathbf r - \mathbf r_0 \right|^2} + + where *C* is ``self.coeff``, *a* is ``self.alpha``, and :math:`\mathbf r_0` is ``self.center``. - where *C* is ``self.coeff``, *a* is ``self.exp``, and :math:`\mathbf r_0` is ``self.center``. + 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, coeff=None): + def __init__(self, center, alpha, coeff=None, normalized=True): self.center = u.array(center) - self.exp = exp - - if coeff is None: - self.coeff = 1.0 # dummy value overwritten by self.normalize() - self.normalize() - else: - self.coeff = coeff + self.alpha = alpha + super().__init__(coeff=coeff, normalized=normalized) def __call__(self, coords, _getvals=False): """ Evaluate this function at the given coordinates. @@ -110,7 +117,7 @@ def __call__(self, coords, _getvals=False): 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) + 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: @@ -121,7 +128,7 @@ def __repr__(self): "exponent: {exp:4.2f}, " "center: {center}>").format( ndim=self.ndim, - center=self.center, exp=self.exp, + center=self.center, exp=self.alpha, coeff=self.coeff) @property @@ -132,22 +139,24 @@ def ndim(self): def __mul__(self, other): """ Returns a new gaussian centroid. """ - a = self.exp - b = other.exp - exp = a + b - center = (a*self.center + b*other.center)/(a+b) + 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(-(a*self.center[i]**2 + b*other.center[i]**2) + - exp * center[i]**2) + prefactor *= np.exp(-(a1*x1[i]**2 + a2*x2[i]**2) + + alpha * center[i]**2) - return Gaussian(center, exp, coeff=prefactor) + return Gaussian(center, alpha, coeff=prefactor, normalized=False) @property def integral(self): r"""Integral of this this gaussian over all N-dimensional space """ - return self.coeff * (np.pi/self.exp)**(self.ndim/2.0) + return self.coeff * (np.pi/self.alpha)**(self.ndim/2.0) def cart_to_powers(s): diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py index 96ddb39..0bcf4ee 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -19,16 +19,19 @@ import numpy as np from scipy.special import factorial -from . import Gaussian +from ..mathutils import spherical_harmonics +from .. import utils +from . import AbstractFunction, Gaussian, CartesianGaussian -class SphericalGaussian(Gaussian): + +class SphericalGaussian(AbstractFunction): 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^{-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``, + where *C* is ``self.coeff``, *a* is ``self.alpha``, and :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 _real_ spherical harmonics. @@ -38,62 +41,46 @@ class SphericalGaussian(Gaussian): 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 - ndims = 3 - - def __init__(self, center, exp, l, m, coeff=None): - from moldesign.mathutils import spherical_harmonics - - super().__init__(center, exp, coeff=1.0) # temporarily set coefficient - - self.center = center - assert len(self.center) == 3 - self.exp = exp + 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 - - if coeff is None: - self.coeff = 1.0 - self.normalize() - else: - self.coeff = coeff - + 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}, " - "exponent: {exp:4.2f}, " + "width: {alpha:4.2f}, " "(l,m) = {qnums}").format( - center=self.center, exp=self.exp, coeff=self.coeff, + center=self.center, alpha=self.alpha, coeff=self.coeff, qnums=(self.l, self.m)) def __mul__(self, other): - raise NotImplementedError("Cannot multiply spherical gaussian functions") + raise NotImplementedError("Cannot yet multiply spherical gaussian functions") - def __call__(self, coords, _include_angular=True): + 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: - _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]): 3D Coordinates or list of 3D coordinates Returns: Scalar or Vector: function value(s) at the passed coordinates """ - result, disp, r2 = super().__call__(coords, _getvals=True) - - if _include_angular: - result *= r2**(self.l/2.0) * self._spherical_harmonic(disp) - return result + 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 spherical basis function: + r""" Overlap of this function with another gaussian basis function .. math:: \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r @@ -108,13 +95,15 @@ def overlap(self, other, normalized=False): from ..data import ODD_FACTORIAL if (self.center != other.center).any(): - raise NotImplementedError("Two-center spherical harmonic overlaps not yet supported") + # work in the cartesian basis + carts = self.to_cart() + raise NotImplementedError() elif self.l != other.l or self.m != other.m: return 0.0 else: - newexp = self.exp + other.exp + newexp = self.alpha + other.alpha power = self.l + other.l + 2 # In this case, radial part integrates to 1, so we just @@ -142,5 +131,16 @@ def angular_part(self, coords): 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) + for cart_fn in spherical_harmonics.SPHERE_TO_CART[self.l, self.m]] + return carts - return r2**(self.l/2.0) * self._spherical_harmonic(disp) \ No newline at end of file From fcb2b51b65344a696126ddab6eeff1d2a0eb2033 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 09:55:42 -0700 Subject: [PATCH 19/33] Rename AbstractFn -> Primitve, BasisContraction -> PrimitiveSum --- moldesign/_tests/test_gaussian_math.py | 72 +++++++++++++------ moldesign/orbitals/__init__.py | 2 +- moldesign/orbitals/atomic_basis_fn.py | 4 +- moldesign/orbitals/cartesian.py | 19 +++-- moldesign/orbitals/gaussians.py | 52 +------------- .../{contraction.py => primitives.py} | 54 +++++++++++++- moldesign/orbitals/spherical.py | 24 ++++--- 7 files changed, 134 insertions(+), 93 deletions(-) rename moldesign/orbitals/{contraction.py => primitives.py} (59%) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 7294146..3d77713 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -19,7 +19,6 @@ 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__) @@ -47,7 +46,7 @@ def cart_3d_gaussian(): g = moldesign.orbitals.CartesianGaussian( center=[random.random() for i in range(3)]*u.angstrom, powers=[1, 3, 0], - exp=1.1/u.angstrom ** 2, + alpha=1.1/u.angstrom ** 2, coeff=1.0) return g @@ -57,7 +56,7 @@ def spherical_3d_gaussian(): g = moldesign.orbitals.SphericalGaussian( center=[random.random() for i in range(3)]*u.angstrom, l=3, m=-2, - exp=1.1/u.angstrom ** 2, + alpha=1.1/u.angstrom ** 2, coeff=1.0) return g @@ -66,15 +65,25 @@ def spherical_3d_gaussian(): @pytest.mark.screening def test_gaussian_integral_and_dimensionality(objkey, request): g = request.getfixturevalue(objkey) - assert g.ndims == len(g.center) + 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) +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_cart_gaussian(powers, withunits=True): if withunits: length = u.angstrom @@ -96,6 +105,14 @@ def _make_rando_spherical_gaussian(l,m, withunits=True): l,m, coeff=random.random()) + +@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])) @@ -111,7 +128,10 @@ def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): """ g1 = _make_rando_cart_gaussian(p1, withunits) g2 = _make_rando_cart_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 @@ -132,7 +152,7 @@ def test_gaussian_function_values(objkey, request): 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) @@ -147,7 +167,7 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): 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 @@ -210,7 +230,7 @@ def _gfuncval(g, coord): r2 = np.sum(r**2, axis=1) else: r2 = np.sum(r**2) - fv = g.coeff * np.exp(-g.exp * r2) + 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 @@ -247,7 +267,7 @@ def test_numerical_vs_analytical_norm(key, request): g = request.getfixturevalue(key) ananorm = g.norm - width = np.sqrt(1/g.exp) + width = np.sqrt(1/g.alpha) ranges = np.ones((3,2)) * 5.0 * width ranges[:,0] *= -1 ranges += g.center[:, None] @@ -265,7 +285,6 @@ def test_numerical_vs_analytical_norm(key, request): def test_s_orbitals_equivalent_among_gaussian_types(): center = np.random.rand(3) * u.angstrom exp = 5.12 / u.angstrom**2 - testcoords = 6.0*(np.random.rand(50, 3)-0.5) * u.angstrom g_bare = moldesign.orbitals.Gaussian(center, exp) g_cart = moldesign.orbitals.CartesianGaussian(center, exp, [0,0,0]) @@ -276,15 +295,9 @@ def test_s_orbitals_equivalent_among_gaussian_types(): gauss.coeff = gauss.coeff / gauss(center) assert gauss(center) == 1.0 - barevals = g_bare(testcoords) - cartvals = g_cart(testcoords) - spherevals = g_sphr(testcoords) + _assert_orbitals_equivalent(g_bare, g_cart) + _assert_orbitals_equivalent(g_bare, g_sphr) - helpers.assert_almost_equal(barevals, cartvals) - helpers.assert_almost_equal(barevals, spherevals) - - helpers.assert_almost_equal(g_bare.norm, g_cart.norm) - helpers.assert_almost_equal(g_sphr.norm, g_cart.norm) LM_TO_CART = {(1,-1): (0,1,0), @@ -303,12 +316,27 @@ def test_s_orbitals_equivalent_among_gaussian_types(): def test_orbitals_same_in_cartesian_and_spherical(lm, powers): center = np.random.rand(3) * u.angstrom exp = 5.12 / u.angstrom**2 - testcoords = 6.0*(np.random.rand(50, 3)-0.5) * u.angstrom g_cart = moldesign.orbitals.CartesianGaussian(center, exp, powers) g_sphr = moldesign.orbitals.SphericalGaussian(center, exp, *lm) - helpers.assert_almost_equal(g_cart(testcoords), - g_sphr(testcoords)) + _assert_orbitals_equivalent(g_cart, g_sphr) + + +@pytest.mark.parametrize('l', range(4)) +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) - helpers.assert_almost_equal(g_sphr.norm, g_cart.norm) + testcoords = 6.0*(np.random.rand(50, 3)-0.5)*u.angstrom + helpers.assert_almost_equal(g1(testcoords), + g2(testcoords)) \ No newline at end of file diff --git a/moldesign/orbitals/__init__.py b/moldesign/orbitals/__init__.py index dec4e13..4d3ae95 100644 --- a/moldesign/orbitals/__init__.py +++ b/moldesign/orbitals/__init__.py @@ -3,9 +3,9 @@ def toplevel(o): return o __all__ = [] +from .primitives import * from .gaussians import * from .orbitals import * -from .contraction import * from .cartesian import * from .spherical import * from .atomic_basis_fn import * diff --git a/moldesign/orbitals/atomic_basis_fn.py b/moldesign/orbitals/atomic_basis_fn.py index ae31b66..a182d16 100644 --- a/moldesign/orbitals/atomic_basis_fn.py +++ b/moldesign/orbitals/atomic_basis_fn.py @@ -16,10 +16,10 @@ # 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 BasisContraction, SHELLS, SPHERICALNAMES +from . import PrimitiveSum, SHELLS, SPHERICALNAMES -class AtomicBasisFunction(BasisContraction): +class AtomicBasisFunction(PrimitiveSum): """ Stores an atomic basis function. Note: diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index 14dff6a..ab9a25d 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -23,26 +23,28 @@ from .. import units as u from .. import utils -from . import AbstractFunction, BasisContraction, Gaussian +from . import Primitive, PrimitiveSum, Gaussian -class CartesianGaussian(AbstractFunction): +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^{-a |\mathbf r - \mathbf{r}_0|^2} + 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``, *alpha* is ``self.alpha``, *r0* is ``self.center``, and + 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`` - References: - Levine, Ira N. Quantum Chemistry, 5th ed. Prentice Hall, 2000. 486-94. + 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 @@ -57,6 +59,9 @@ class CartesianGaussian(AbstractFunction): 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') @@ -189,7 +194,7 @@ def __mul__(self, other): elif len(new_gaussians) == 1: return new_gaussians[0] else: - return BasisContraction(new_gaussians) + return PrimitiveSum(new_gaussians) @property def integral(self): diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index ab0425d..7878c34 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -16,61 +16,13 @@ # 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 numpy as np from .. import units as u +from . import Primitive -class AbstractFunction(object): - """ Abstract base class for basis functions - """ - def __init__(self, coeff=None, normalized=False): - if normalized or coeff is None: - self.coeff = 1.0 # dummy value overwritten by self.normalize() - self.normalize() - - if coeff is not None: - self.coeff = coeff - - 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 (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 /= (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)) - - -class Gaussian(AbstractFunction): +class Gaussian(Primitive): r""" N-dimensional gaussian function. The function is given by: diff --git a/moldesign/orbitals/contraction.py b/moldesign/orbitals/primitives.py similarity index 59% rename from moldesign/orbitals/contraction.py rename to moldesign/orbitals/primitives.py index bba75d9..6f34bc6 100644 --- a/moldesign/orbitals/contraction.py +++ b/moldesign/orbitals/primitives.py @@ -16,12 +16,60 @@ # 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 numpy as np -from .orbitals import Orbital -class BasisContraction(Orbital): +class Primitive(object): + """ Abstract base class for basis function primitives + """ + def __init__(self, coeff=None, normalized=False): + if normalized or coeff is None: + self.coeff = 1.0 # dummy value overwritten by self.normalize() + self.normalize() + + if coeff is not None: + self.coeff = coeff + + 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. + + 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)) + + +class PrimitiveSum(object): """ Stores a linear combination of primitive functions. Args: @@ -70,4 +118,4 @@ def __str__(self): return '%s with %d primitives' % (self.__class__.__name__, self.num_primitives) def __repr__(self): - return '<%s>' % self \ No newline at end of file + return '<%s>' % self diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py index 0bcf4ee..9207338 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -22,18 +22,26 @@ from ..mathutils import spherical_harmonics from .. import utils -from . import AbstractFunction, Gaussian, CartesianGaussian +from . import Primitive, Gaussian, CartesianGaussian, PrimitiveSum -class SphericalGaussian(AbstractFunction): +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^{-a\left| \mathbf r - \mathbf r_0 \right|^2} + 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``, *a* is ``self.alpha``, and :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 _real_ - spherical harmonics. + 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. @@ -86,7 +94,7 @@ def overlap(self, other, normalized=False): \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r Args: - other (AbstractFunction): + other (Primitive): normalized (bool): If True, return the overlap of the two NORMALIZED functions. Returns: @@ -142,5 +150,5 @@ def to_cart(self): carts = [CartesianGaussian(self.center, self.alpha, cart_fn.powers, self.coeff * cart_fn.coeff) for cart_fn in spherical_harmonics.SPHERE_TO_CART[self.l, self.m]] - return carts + return PrimitiveSum(carts) From 79509f948a3360580721374577151c6f8664b178 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 13:48:01 -0700 Subject: [PATCH 20/33] Vendor PyQuante 1-el integral libs --- NOTICES | 70 ++++--- moldesign/external/pyquante2/LICENSE | 28 +++ moldesign/external/pyquante2/__init__.py | 0 moldesign/external/pyquante2/one.py | 253 +++++++++++++++++++++++ moldesign/external/pyquante2/utils.py | 242 ++++++++++++++++++++++ 5 files changed, 564 insertions(+), 29 deletions(-) create mode 100644 moldesign/external/pyquante2/LICENSE create mode 100644 moldesign/external/pyquante2/__init__.py create mode 100644 moldesign/external/pyquante2/one.py create mode 100644 moldesign/external/pyquante2/utils.py 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/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() + From fd4c7eb951270add178d9e73c33d28fcd79641f7 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 13:49:36 -0700 Subject: [PATCH 21/33] Add multiple __PYTEST_MARK__ marks --- moldesign/_tests/conftest.py | 9 +++- moldesign/_tests/test_gaussian_math.py | 67 ++++++++++++++++++++++---- moldesign/_tests/test_mathutils.py | 59 ++++++++++++++++++++++- moldesign/_tests/test_units.py | 3 +- moldesign/_tests/test_wfn.py | 58 ---------------------- 5 files changed, 124 insertions(+), 72 deletions(-) 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/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 3d77713..8cdc2b7 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -15,6 +15,8 @@ registered_types = {} +__PYTEST_MARK__ = ['math', 'gaussians'] + def typedfixture(*types, **kwargs): """This is a decorator that lets us associate fixtures with one or more arbitrary types. @@ -106,6 +108,44 @@ def _make_rando_spherical_gaussian(l,m, withunits=True): coeff=random.random()) +@pytest.mark.parametrize('withunits', [True, False], ids=['units','c-num']) +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=['units','c-num']) +def test_numerical_vs_analytical_overlap_cart(withunits): + p1 = _make_rando_cart_gaussian((1,2,3), withunits) + p2 = _make_rando_cart_gaussian((1,0,1), withunits) + _assert_numerical_analytical_overlaps_match(p1, p2) + + +@pytest.mark.parametrize('withunits', [True, False], ids=['units','c-num']) +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) + + +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 = _generate_grid(g1, g2) + with np.errstate(under='ignore'): + prodvals = g1(allpoints) * g2(allpoints) + numsum = prodvals.sum() * grid.dx ** 3 + 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) @@ -160,7 +200,9 @@ def test_gaussian_function_values(objkey, request): @pytest.mark.parametrize('objkey', - registered_types['gaussian']) + registered_types['gaussian']+ + registered_types['cart-gaussian'] + + registered_types['spherical-gaussian']) def test_vectorized_gaussian_function_evaluations(objkey, request): g = request.getfixturevalue(objkey) @@ -267,12 +309,7 @@ def test_numerical_vs_analytical_norm(key, request): g = request.getfixturevalue(key) ananorm = g.norm - width = np.sqrt(1/g.alpha) - ranges = np.ones((3,2)) * 5.0 * width - ranges[:,0] *= -1 - ranges += g.center[:, None] - grid = moldesign.mathutils.VolumetricGrid(*ranges, 25) - allpoints = grid.allpoints() + allpoints, grid = _generate_grid(g) with np.errstate(under='ignore'): vals = g(allpoints) @@ -281,6 +318,17 @@ def test_numerical_vs_analytical_norm(key, request): helpers.assert_almost_equal(ananorm, numnorm) +def _generate_grid(g, g2=None): + if g2 is None: g2 = g + width = max(np.sqrt(1/g.alpha), np.sqrt(1/g2.alpha)) + ranges = np.ones((3, 2))*5.0*width + ranges[:, 0] *= -1 + ranges += ((g.center + g2.center)/2.0)[:, None] + grid = moldesign.mathutils.VolumetricGrid(*ranges, 64) + allpoints = grid.allpoints() + return allpoints, grid + + @pytest.mark.screening def test_s_orbitals_equivalent_among_gaussian_types(): center = np.random.rand(3) * u.angstrom @@ -299,7 +347,6 @@ def test_s_orbitals_equivalent_among_gaussian_types(): _assert_orbitals_equivalent(g_bare, g_sphr) - LM_TO_CART = {(1,-1): (0,1,0), (1,0): (0,0,1), (1,1): (1,0,0), @@ -323,7 +370,7 @@ def test_orbitals_same_in_cartesian_and_spherical(lm, powers): _assert_orbitals_equivalent(g_cart, g_sphr) -@pytest.mark.parametrize('l', range(4)) +@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 @@ -339,4 +386,4 @@ def _assert_orbitals_equivalent(g1, g2): testcoords = 6.0*(np.random.rand(50, 3)-0.5)*u.angstrom helpers.assert_almost_equal(g1(testcoords), - g2(testcoords)) \ No newline at end of file + g2(testcoords)) diff --git a/moldesign/_tests/test_mathutils.py b/moldesign/_tests/test_mathutils.py index 84b5ddd..de63414 100644 --- a/moldesign/_tests/test_mathutils.py +++ b/moldesign/_tests/test_mathutils.py @@ -10,7 +10,7 @@ from . import helpers -__PYTEST_MARK__ = 'internal' +__PYTEST_MARK__ = ['internal', 'math', 'mathutils'] CONSTRUCTORS = [] IDS = [] @@ -173,6 +173,63 @@ def test_spherical_harmonics_orthonormal(l1, m1, l2, m2): 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, 3, 4, 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, 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], diff --git a/moldesign/_tests/test_units.py b/moldesign/_tests/test_units.py index 28051a2..2dc1a1a 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 labels (see ./conftest.py) +__PYTEST_MARK__ = ['internal', 'units'] @pytest.mark.screening diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py index 4753a92..cd4fdb8 100644 --- a/moldesign/_tests/test_wfn.py +++ b/moldesign/_tests/test_wfn.py @@ -58,61 +58,3 @@ def test_pyscf_and_mdt_grids_are_the_same(molkey, request): pyscf_vals = basis_values(mol, mol.wfn.aobasis, randocoords) mdt_vals = mol.wfn.aobasis(randocoords) np.testing.assert_allclose(mdt_vals, pyscf_vals) - - -@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 = mdt.mathutils.VolumetricGrid(*ranges, 3, 4, 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 = mdt.mathutils.VolumetricGrid(*ranges, 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] - - - From 0449dcc0a863cbacfc0a58a562cdda293a452838 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 13:56:38 -0700 Subject: [PATCH 22/33] Much faster grid generation --- moldesign/mathutils/grids.py | 49 ++++++++++++++++++++++++++++++------ 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/moldesign/mathutils/grids.py b/moldesign/mathutils/grids.py index 57b9d6d..d1ee73c 100644 --- a/moldesign/mathutils/grids.py +++ b/moldesign/mathutils/grids.py @@ -42,7 +42,7 @@ class VolumetricGrid(object): xspace, yspace, zspace = (utils.IndexView('spaces', i) for i in range(3)) def __init__(self, xrange, yrange, zrange, - xpoints=25, ypoints=None, zpoints=None): + xpoints=32, ypoints=None, zpoints=None): self.points = np.array([xpoints, ypoints if ypoints is not None else xpoints, @@ -53,6 +53,11 @@ def __init__(self, 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) @@ -81,7 +86,7 @@ def iter_points(self): for i,j,k in itertools.product(*map(range, self.points)): yield self.origin + self.deltas * [i,j,k] - def allpoints(self, dtype='float'): + 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 @@ -90,15 +95,43 @@ def allpoints(self, dtype='float'): Yields: Matrix[shape=(self.npoints**3,3)]: x,y,z coordinate of each point on the grid """ - ap = np.empty((self.npoints, 3), dtype=dtype) + return _cartesian_product(self.spaces) - for ip, point in enumerate(self.iter_points()): - ap[ip] = point - if hasattr(self.xr, 'units'): - ap = ap * self.xr.units +def _cartesian_product(arrays): + """ Fast grid creation routine from @senderle on Stack Overflow. This is awesome. - return ap + 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 From 66fc58f30f9c7ace38427c5ab2fbd23ed9b91298 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 13:57:41 -0700 Subject: [PATCH 23/33] Delegate to PyQuante overlap routines in most cases --- moldesign/mathutils/spherical_harmonics.py | 12 ++++- moldesign/orbitals/cartesian.py | 52 +++++++++++++++++++++- moldesign/orbitals/gaussians.py | 5 +++ moldesign/orbitals/primitives.py | 31 +++++++++++-- moldesign/orbitals/spherical.py | 11 ++--- 5 files changed, 99 insertions(+), 12 deletions(-) diff --git a/moldesign/mathutils/spherical_harmonics.py b/moldesign/mathutils/spherical_harmonics.py index 047e43d..96cb60e 100644 --- a/moldesign/mathutils/spherical_harmonics.py +++ b/moldesign/mathutils/spherical_harmonics.py @@ -36,12 +36,12 @@ def cart_to_polar_angles(coords): class Y(object): - r""" A *real*-valued spherical harmonic function + r""" A real-valued spherical harmonic function These functions are orthonormalized over the sphere such that .. math:: - \int d\Omega \: Y_{lm}(\theta, \phi) Y_{l'm'}(\theta, \phi) = \delta_{l,l'} \delta_{m,m'} + \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 @@ -98,6 +98,9 @@ def __call__(self, coords): 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): @@ -134,6 +137,11 @@ def __call__(self, coords): 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)) diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index ab9a25d..0d0d84c 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -91,6 +91,54 @@ def angular(self): """ 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. @@ -144,10 +192,10 @@ def __mul__(self, other): angular parts. Args: - other (CartesianGaussian): other gaussian wavepacket + other (CartesianGaussian): other gaussian basis function Returns: - CartesianGaussian or BasisContraction: product functions + CartesianGaussian or PrimitiveSum[CartesianGaussian]: product functions """ if (self.center == other.center).all(): # this case is much easier than the general one diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index 7878c34..e9ec74e 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -110,6 +110,11 @@ def integral(self): """ return self.coeff * (np.pi/self.alpha)**(self.ndim/2.0) + 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): """ Convert a string to a list of cartesian powers diff --git a/moldesign/orbitals/primitives.py b/moldesign/orbitals/primitives.py index 6f34bc6..fd51ab9 100644 --- a/moldesign/orbitals/primitives.py +++ b/moldesign/orbitals/primitives.py @@ -17,12 +17,14 @@ # 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 + """ Abstract base class for basis function primitives. All functions are assumed real. """ def __init__(self, coeff=None, normalized=False): if normalized or coeff is None: @@ -50,6 +52,11 @@ def overlap(self, other, normalized=False): 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 """ @@ -68,13 +75,21 @@ def norm(self): """ return np.sqrt(self.overlap(self)) + def iterprimitives(self): + yield self + class PrimitiveSum(object): """ Stores a linear combination of primitive functions. Args: - primitives (List[PrimitiveBase]): List of primitives, if available + 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 @@ -111,11 +126,21 @@ def overlap(self, other): Args: other (AbstractFunction or Orbital): object to calculate overlaps with """ - olap = sum(p1.overlap(other) for p1 in self.primitives) + 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 index 9207338..b555a6a 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -88,7 +88,7 @@ def __call__(self, coords): return result * self.coeff def overlap(self, other, normalized=False): - r""" Overlap of this function with another gaussian basis function + r""" Overlap of this function with another basis function .. math:: \int f_1(\mathbf r) f_2(\mathbf r) d^N \mathbf r @@ -102,10 +102,8 @@ def overlap(self, other, normalized=False): """ from ..data import ODD_FACTORIAL - if (self.center != other.center).any(): - # work in the cartesian basis - carts = self.to_cart() - raise NotImplementedError() + 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 @@ -152,3 +150,6 @@ def to_cart(self): 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 From 635290360d50bba4dc23f7f97d29d4946b70da84 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 16:23:04 -0700 Subject: [PATCH 24/33] Handle more units corner cases --- moldesign/_tests/test_units.py | 43 +++++++++++++++++++++++++++- moldesign/units/quantity.py | 4 +++ moldesign/units/tools.py | 51 +++++++++++++++++++++++----------- moldesign/units/unitsystem.py | 11 ++++---- 4 files changed, 87 insertions(+), 22 deletions(-) diff --git a/moldesign/_tests/test_units.py b/moldesign/_tests/test_units.py index 2dc1a1a..1786fc3 100644 --- a/moldesign/_tests/test_units.py +++ b/moldesign/_tests/test_units.py @@ -6,7 +6,7 @@ from moldesign import units from moldesign import units as u -# mark all tests in this module with these labels (see ./conftest.py) +# mark all tests in this module with these label (see ./conftest.py) __PYTEST_MARK__ = ['internal', 'units'] @@ -104,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 diff --git a/moldesign/units/quantity.py b/moldesign/units/quantity.py index 2e5f98f..ca6f243 100644 --- a/moldesign/units/quantity.py +++ b/moldesign/units/quantity.py @@ -42,6 +42,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 62c5061..46273f3 100644 --- a/moldesign/units/tools.py +++ b/moldesign/units/tools.py @@ -173,14 +173,22 @@ def get_units(q): 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 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 @@ -189,24 +197,35 @@ def array(qlist, defunits=False, _baseunit=None): if _baseunit is None: _baseunit = get_units(qlist) - try: - if _baseunit == 1.0: - # It's dimensionless, return an ndarray - return np.array(qlist) - except DimensionalityError: - pass - - if defunits and isinstance(_baseunit, MdtUnit): + if _baseunit.dimensionless: + return _make_nparray(qlist) + if defunits: _baseunit = default.get_default(_baseunit) - try: - # Convert all sublists to the common "baseunit" - newlist = [array(item, _baseunit=_baseunit).value_in(_baseunit) for item in qlist] - return _baseunit*newlist - except TypeError as exc: - # Otherwise, + 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: + 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) def broadcast_to(arr, *args, **kwargs): diff --git a/moldesign/units/unitsystem.py b/moldesign/units/unitsystem.py index f86b1e0..46413e0 100644 --- a/moldesign/units/unitsystem.py +++ b/moldesign/units/unitsystem.py @@ -89,8 +89,7 @@ def convert(self, quantity): 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) @@ -131,14 +130,16 @@ 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 + return ureg.dimensionless + baseunit = ureg.dimensionless # Factor out force units if self._force: @@ -181,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) From 24c3dd5778363a7d5c25eb13e5e0604cc51e2c87 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 22:52:32 -0700 Subject: [PATCH 25/33] Create a "nano" base unit system in pint Unclear if this actually does any good, but it doesn't appear to _hurt_ anything --- moldesign/_static_data/pint_atomic_units.txt | 14 ++++++++------ moldesign/_tests/object_fixtures.py | 2 +- moldesign/_tests/test_units.py | 2 +- moldesign/units/constants.py | 10 +++++----- moldesign/units/quantity.py | 3 ++- moldesign/units/tools.py | 2 +- requirements.txt | 2 +- 7 files changed, 19 insertions(+), 16 deletions(-) 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/object_fixtures.py b/moldesign/_tests/object_fixtures.py index 08725d5..aa4098b 100644 --- a/moldesign/_tests/object_fixtures.py +++ b/moldesign/_tests/object_fixtures.py @@ -62,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 ###################################### diff --git a/moldesign/_tests/test_units.py b/moldesign/_tests/test_units.py index 1786fc3..172aa7d 100644 --- a/moldesign/_tests/test_units.py +++ b/moldesign/_tests/test_units.py @@ -246,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/units/constants.py b/moldesign/units/constants.py index 13fe61c..c8f2ab3 100644 --- a/moldesign/units/constants.py +++ b/moldesign/units/constants.py @@ -42,9 +42,9 @@ electron_charge = q_e = ureg.elementary_charge # useful units -fs = femtoseconds = ureg.fs -ps = picoseconds = ureg.ps -ns = nanoseconds = ureg.ns +fs = femtoseconds = ureg.femtosecond +ps = picoseconds = ureg.picosecond +ns = nanoseconds = ureg.nanosecond eV = electronvolts = ureg.eV kcalpermol = ureg.kcalpermol gpermol = ureg.gpermol @@ -54,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 @@ -65,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 ca6f243..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) diff --git a/moldesign/units/tools.py b/moldesign/units/tools.py index 46273f3..933b43f 100644 --- a/moldesign/units/tools.py +++ b/moldesign/units/tools.py @@ -146,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: 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 From e38cb8a7e5aac42d1a12de777f103749800a7d04 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 22:53:14 -0700 Subject: [PATCH 26/33] Correct normalization of basis functions from PySCF --- moldesign/models/pyscf.py | 7 ++++--- moldesign/orbitals/primitives.py | 10 ++++++---- moldesign/orbitals/spherical.py | 5 +++-- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/moldesign/models/pyscf.py b/moldesign/models/pyscf.py index 728c97b..7d6668a 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) / (u.bohr**2) # very unsure if these units are right! + exps = pmol.bas_exp(ishell) / (u.bohr**2) num_contractions = pmol.bas_nctr(ishell) coeffs = pmol.bas_ctr_coeff(ishell) @@ -455,7 +455,8 @@ def _get_ao_basis_functions(self): primitives = [orbitals.SphericalGaussian(atom.position.copy(), exp, l, m, - coeff=coeff[ictr]) + 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/orbitals/primitives.py b/moldesign/orbitals/primitives.py index fd51ab9..d909337 100644 --- a/moldesign/orbitals/primitives.py +++ b/moldesign/orbitals/primitives.py @@ -27,12 +27,13 @@ class Primitive(object): """ Abstract base class for basis function primitives. All functions are assumed real. """ def __init__(self, coeff=None, normalized=False): - if normalized or coeff is None: - self.coeff = 1.0 # dummy value overwritten by self.normalize() + 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 + self.coeff *= coeff def copy(self): return copy.deepcopy(self) @@ -127,7 +128,8 @@ def overlap(self, other): other (AbstractFunction or Orbital): object to calculate overlaps with """ olap = sum(p1.overlap(p2) - for p1, p2 in itertools.product(self.iterprimitives(), other.iterprimitives())) + for p1, p2 in itertools.product(self.iterprimitives(), + other.iterprimitives())) return olap @property diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py index b555a6a..9d5e414 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -102,7 +102,7 @@ def overlap(self, other, normalized=False): """ from ..data import ODD_FACTORIAL - if (self.center - other.center).any(): + if (self.center != other.center).any(): return self.to_cart().overlap(other) elif self.l != other.l or self.m != other.m: @@ -146,7 +146,8 @@ def to_cart(self): List[CartesianGaussian]: cartesian components of this function """ carts = [CartesianGaussian(self.center, self.alpha, - cart_fn.powers, self.coeff * cart_fn.coeff) + 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) From 52fc707acff85cecc2ea845cb0178f6fd679de5e Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 23:04:22 -0700 Subject: [PATCH 27/33] Immediately convert gaussians to default units for stability --- moldesign/_tests/test_gaussian_math.py | 61 ++++++++++++++++++++++++++ moldesign/_tests/test_molecules.py | 6 +-- moldesign/orbitals/gaussians.py | 4 +- 3 files changed, 66 insertions(+), 5 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 8cdc2b7..6d1a538 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -12,6 +12,7 @@ from moldesign.mathutils import spherical_harmonics as harmonics from . import helpers +from .molecule_fixtures import * registered_types = {} @@ -183,6 +184,40 @@ def _assert_same_function_values(g1, g2, withunits): helpers.assert_almost_equal(prodvals, gvals) +def test_initial_normalization_gaussian(): + center = np.random.rand(3) * u.angstrom + exp = 5.12 / u.angstrom**2 + g1 = moldesign.orbitals.Gaussian(center, exp, coeff=3.0, normalized=True) + assert abs(3.0 - g1.norm) < 1e-12 + + g1 = moldesign.orbitals.Gaussian(center, exp, normalized=True) + assert abs(1.0 - g1.norm) < 1e-12 + + +def test_initial_normalization_cart(): + 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) + assert abs(3.0 - g1.norm) < 1e-12 + + g1 = moldesign.orbitals.CartesianGaussian(center, exp, powers, normalized=True) + assert abs(1.0 - g1.norm) < 1e-12 + + +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): + g1 = moldesign.orbitals.SphericalGaussian(center, exp, l, m, + coeff=3.0, normalized=True) + assert abs(3.0 - g1.norm) < 1e-12 + + g1 = moldesign.orbitals.SphericalGaussian(center, exp, l, m, normalized=True) + assert abs(1.0 - g1.norm) < 1e-12 + + @pytest.mark.parametrize('objkey', registered_types['gaussian'] + registered_types['cart-gaussian'] + @@ -347,6 +382,7 @@ def test_s_orbitals_equivalent_among_gaussian_types(): _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), @@ -387,3 +423,28 @@ def _assert_orbitals_equivalent(g1, g2): 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) diff --git a/moldesign/_tests/test_molecules.py b/moldesign/_tests/test_molecules.py index 7902dee..133cf8c 100644 --- a/moldesign/_tests/test_molecules.py +++ b/moldesign/_tests/test_molecules.py @@ -87,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 @@ -96,7 +96,7 @@ 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 + h2_harmonic.atoms[1].x += 0.3*u.angstrom with pytest.raises(mdt.NotCalculatedError): h2_harmonic.forces with pytest.raises(mdt.NotCalculatedError): @@ -105,7 +105,7 @@ def test_h2_not_calculated_yet(h2_harmonic): def h2_properties_raises_not_calculated_yet(h2_harmonic): h2_harmonic.calculate() - h2_harmonic.atoms[1].x += 0.3*u.ang + h2_harmonic.atoms[1].x += 0.3*u.angstrom with pytest.raises(mdt.NotCalculatedError): h2_harmonic.properties.forces with pytest.raises(mdt.NotCalculatedError): diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index e9ec74e..b8aecb9 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -40,8 +40,8 @@ class Gaussian(Primitive): normalized (bool): if True, normalize the function before applying the coefficient """ def __init__(self, center, alpha, coeff=None, normalized=True): - self.center = u.array(center) - self.alpha = alpha + self.center = u.default.convert(center) + self.alpha = u.default.convert(alpha) super().__init__(coeff=coeff, normalized=normalized) def __call__(self, coords, _getvals=False): From ed785cabd816cb1d0a3de325edf30034ce2ec86b Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Tue, 11 Jul 2017 23:16:27 -0700 Subject: [PATCH 28/33] Fix PySCF label parsing --- moldesign/_tests/test_gaussian_math.py | 3 ++- moldesign/orbitals/orbitals.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 6d1a538..6d719ce 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -447,4 +447,5 @@ def test_pyscf_and_mdt_overlaps_are_the_same(h2_rhf_augccpvdz): assert isinstance(overlaps, np.ndarray) or overlaps.units == u.dimensionless np.testing.assert_allclose(mol.wfn.aobasis.overlaps, - overlaps) + overlaps, + atol=5.0e-7) 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), From 0df1a41362bdd771b7fbfc076685aa357d206508 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Fri, 14 Jul 2017 22:39:14 -0700 Subject: [PATCH 29/33] Fix a couple minimization trajectory bugs --- moldesign/_tests/test_minimizers.py | 2 ++ moldesign/min/smart.py | 2 +- moldesign/molecules/trajectory.py | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) 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/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/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() From 1dd9619ebccbd7febe73cb084292b4ff38945080 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 17 Jul 2017 13:04:45 -0700 Subject: [PATCH 30/33] Fix units on wavefunction values --- moldesign/_tests/helpers.py | 17 +++ moldesign/_tests/test_gaussian_math.py | 175 ++++++++++++++---------- moldesign/_tests/test_wfn.py | 31 ++++- moldesign/interfaces/pyscf_interface.py | 4 +- moldesign/orbitals/basis.py | 4 +- moldesign/orbitals/cartesian.py | 3 + moldesign/orbitals/gaussians.py | 3 + moldesign/orbitals/primitives.py | 18 ++- moldesign/orbitals/spherical.py | 4 + moldesign/units/tools.py | 7 +- 10 files changed, 180 insertions(+), 86 deletions(-) diff --git a/moldesign/_tests/helpers.py b/moldesign/_tests/helpers.py index e89adc5..86788b9 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, 64) + allpoints = grid.allpoints() + return allpoints, grid diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 6d719ce..02ec734 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -30,22 +30,22 @@ 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('cart-gaussian') -def cart_3d_gaussian(): +@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], @@ -54,7 +54,7 @@ def cart_3d_gaussian(): return g -@typedfixture('spherical-gaussian') +@typedfixture('basis_fn') def spherical_3d_gaussian(): g = moldesign.orbitals.SphericalGaussian( center=[random.random() for i in range(3)]*u.angstrom, @@ -64,7 +64,7 @@ def spherical_3d_gaussian(): 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.getfixturevalue(objkey) @@ -77,6 +77,11 @@ def test_gaussian_integral_and_dimensionality(objkey, request): decimal=10) +@pytest.fixture +def linear_combination(): + return _make_rando_linear_combination(True) + + def _make_rando_gaussian(withunits=True): if withunits: length = u.angstrom @@ -87,7 +92,7 @@ def _make_rando_gaussian(withunits=True): coeff=random.random()) -def _make_rando_cart_gaussian(powers, withunits=True): +def _make_rando_cartesian_gaussian(powers, withunits=True): if withunits: length = u.angstrom else: @@ -109,27 +114,55 @@ def _make_rando_spherical_gaussian(l,m, withunits=True): coeff=random.random()) -@pytest.mark.parametrize('withunits', [True, False], ids=['units','c-num']) +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()/(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=['units','c-num']) -def test_numerical_vs_analytical_overlap_cart(withunits): - p1 = _make_rando_cart_gaussian((1,2,3), withunits) - p2 = _make_rando_cart_gaussian((1,0,1), withunits) +@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=['units','c-num']) +@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: @@ -140,10 +173,10 @@ def _assert_numerical_analytical_overlaps_match(g1, g2): else: helpers.assert_almost_equal(prod.integral, olap) - allpoints, grid = _generate_grid(g1, g2) + allpoints, grid = helpers.generate_grid(g1, g2) with np.errstate(under='ignore'): prodvals = g1(allpoints) * g2(allpoints) - numsum = prodvals.sum() * grid.dx ** 3 + numsum = prodvals.sum() * grid.dx * grid.dy * grid.dz helpers.assert_almost_equal(numsum, olap, decimal=5) @@ -164,11 +197,11 @@ def test_gaussian_multiplication_amplitudes(withunits): @pytest.mark.parametrize('p1,p2,withunits', cartesian_test_suite, ids=cartesian_test_ids) -def test_cart_gaussian_multiplication_amplitudes(p1, p2, withunits): +def test_cartesian_gaussian_multiplication_amplitudes(p1, p2, withunits): """ Tests that ``g1(x) * g2(x) == (g1 * g2)(x)`` """ - g1 = _make_rando_cart_gaussian(p1, withunits) - g2 = _make_rando_cart_gaussian(p2, withunits) + g1 = _make_rando_cartesian_gaussian(p1, withunits) + g2 = _make_rando_cartesian_gaussian(p2, withunits) _assert_same_function_values(g1, g2, withunits) @@ -184,44 +217,70 @@ def _assert_same_function_values(g1, g2, withunits): helpers.assert_almost_equal(prodvals, gvals) -def test_initial_normalization_gaussian(): +def test_initial_gaussian_normalization_gaussian(): center = np.random.rand(3) * u.angstrom exp = 5.12 / u.angstrom**2 - g1 = moldesign.orbitals.Gaussian(center, exp, coeff=3.0, normalized=True) - assert abs(3.0 - g1.norm) < 1e-12 + 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) - g1 = moldesign.orbitals.Gaussian(center, exp, normalized=True) - assert abs(1.0 - g1.norm) < 1e-12 +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_cart(): + +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)): - g1 = moldesign.orbitals.CartesianGaussian(center, exp, powers, coeff=3.0, normalized=True) - assert abs(3.0 - g1.norm) < 1e-12 + 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) - g1 = moldesign.orbitals.CartesianGaussian(center, exp, powers, normalized=True) - assert abs(1.0 - g1.norm) < 1e-12 + +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, normalized=True) - assert abs(3.0 - g1.norm) < 1e-12 + 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) - g1 = moldesign.orbitals.SphericalGaussian(center, exp, l, m, normalized=True) - assert abs(1.0 - g1.norm) < 1e-12 + +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['gaussian'] + - registered_types['cart-gaussian'] + - registered_types['spherical-gaussian']) +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_gaussian_function_values(objkey, request): g = request.getfixturevalue(objkey) @@ -234,10 +293,7 @@ def test_gaussian_function_values(objkey, request): _assert_almost_equal(funcval, retval) -@pytest.mark.parametrize('objkey', - registered_types['gaussian']+ - registered_types['cart-gaussian'] + - registered_types['spherical-gaussian']) +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_vectorized_gaussian_function_evaluations(objkey, request): g = request.getfixturevalue(objkey) @@ -256,20 +312,14 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): _assert_almost_equal(vector_results, expected, decimal=8) -@pytest.mark.parametrize('objkey', - registered_types['gaussian'] + - registered_types['cart-gaussian'] + - registered_types['spherical-gaussian']) +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_gaussian_str_and_repr_works(objkey, request): g1 = request.getfixturevalue(objkey) str(g1) repr(g1) -@pytest.mark.parametrize('objkey', - registered_types['gaussian'] + - registered_types['cart-gaussian'] + - registered_types['spherical-gaussian']) +@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() @@ -283,10 +333,7 @@ def test_normalized_gaussian_self_overlap_is_unity(objkey, request): assert abs(1.0 - olap) < 1e-12 -@pytest.mark.parametrize('objkey', - registered_types['gaussian'] + - registered_types['cart-gaussian'] + - registered_types['spherical-gaussian']) +@pytest.mark.parametrize('objkey', registered_types['basis_fn']) def test_normalization(objkey, request): g1 = request.getfixturevalue(objkey) oldnorm = g1.norm @@ -339,29 +386,11 @@ def test_convert_cartesian_label_to_array_of_integer_powers(): assert cart_to_powers('zx^3') == [3,0,1] -@pytest.mark.parametrize('key', ['std_3d_gaussian', 'cart_3d_gaussian', 'spherical_3d_gaussian']) +@pytest.mark.parametrize('key', registered_types['basis_fn'] + ['linear_combination']) def test_numerical_vs_analytical_norm(key, request): g = request.getfixturevalue(key) - ananorm = g.norm - - allpoints, grid = _generate_grid(g) - - with np.errstate(under='ignore'): - vals = g(allpoints) - - numnorm = np.sqrt(grid.dx**3 * (vals**2).sum()) - helpers.assert_almost_equal(ananorm, numnorm) - - -def _generate_grid(g, g2=None): - if g2 is None: g2 = g - width = max(np.sqrt(1/g.alpha), np.sqrt(1/g2.alpha)) - ranges = np.ones((3, 2))*5.0*width - ranges[:, 0] *= -1 - ranges += ((g.center + g2.center)/2.0)[:, None] - grid = moldesign.mathutils.VolumetricGrid(*ranges, 64) - allpoints = grid.allpoints() - return allpoints, grid + numnorm = _numerical_norm(g) + helpers.assert_almost_equal(g.norm, numnorm) @pytest.mark.screening diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py index cd4fdb8..f94d1e3 100644 --- a/moldesign/_tests/test_wfn.py +++ b/moldesign/_tests/test_wfn.py @@ -9,6 +9,8 @@ from . import helpers +__PYTEST_MARK__ = 'wfn' + TESTSYSTEMS = ['h2_rhf_augccpvdz', 'h2_rhf_sto3g', 'acetylene_dft_631g'] @@ -45,16 +47,31 @@ def test_pyscf_orbital_grid_works(molkey, request): 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)) - np.testing.assert_allclose(vals_plusminus[:,0], - (vals_b0[:,0] + vals_b0[:,1])/np.sqrt(2)) - np.testing.assert_allclose(vals_plusminus[:,1], - (vals_b0[:,0] - vals_b0[:,1])/np.sqrt(2)) + 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_pyscf_and_mdt_grids_are_the_same(molkey, request): +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) - mdt_vals = mol.wfn.aobasis(randocoords) - np.testing.assert_allclose(mdt_vals, pyscf_vals) + 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']) +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/interfaces/pyscf_interface.py b/moldesign/interfaces/pyscf_interface.py index cd6843c..7250da8 100644 --- a/moldesign/interfaces/pyscf_interface.py +++ b/moldesign/interfaces/pyscf_interface.py @@ -121,14 +121,14 @@ 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: diff --git a/moldesign/orbitals/basis.py b/moldesign/orbitals/basis.py index f437348..68d4951 100644 --- a/moldesign/orbitals/basis.py +++ b/moldesign/orbitals/basis.py @@ -87,7 +87,7 @@ def __call__(self, coords, coeffs=None): - if ``coeffs`` is NOT passed, an array of basis function amplitudes of size ``(len(coords), len(aobasis))``. """ - basis_vals = np.zeros((len(coords), len(self))) + 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) @@ -100,7 +100,7 @@ def get_basis_functions_on_atom(self, atom): """ Return a list of basis functions on this atom Args: - atom (moldesign.Atom): + atom (moldesign.Atom): query atom Returns: List[AtomicBasisFunction]: basis functions centered on this atom diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index 0d0d84c..641bed3 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -85,6 +85,9 @@ def __repr__(self): 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) diff --git a/moldesign/orbitals/gaussians.py b/moldesign/orbitals/gaussians.py index b8aecb9..fe38be0 100644 --- a/moldesign/orbitals/gaussians.py +++ b/moldesign/orbitals/gaussians.py @@ -44,6 +44,9 @@ def __init__(self, center, alpha, coeff=None, normalized=True): self.alpha = u.default.convert(alpha) super().__init__(coeff=coeff, normalized=normalized) + def _get_wfn_units(self): + return u.MdtQuantity(self.coeff).units + def __call__(self, coords, _getvals=False): """ Evaluate this function at the given coordinates. diff --git a/moldesign/orbitals/primitives.py b/moldesign/orbitals/primitives.py index d909337..72eff7c 100644 --- a/moldesign/orbitals/primitives.py +++ b/moldesign/orbitals/primitives.py @@ -35,6 +35,9 @@ def __init__(self, coeff=None, normalized=False): if coeff is not None: self.coeff *= coeff + def _get_wfn_units(self): + return NotImplementedError() + def copy(self): return copy.deepcopy(self) @@ -93,9 +96,13 @@ class PrimitiveSum(object): 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)) + outvals = np.zeros(len(coords)) * self._get_wfn_units() for primitive in self.primitives: outvals += primitive(coords) return outvals @@ -104,6 +111,15 @@ def __call__(self, coords): 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{}` diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py index 9d5e414..b402eaa 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -21,6 +21,7 @@ from scipy.special import factorial from ..mathutils import spherical_harmonics +from .. import units as u from .. import utils from . import Primitive, Gaussian, CartesianGaussian, PrimitiveSum @@ -71,6 +72,9 @@ def __repr__(self): 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. diff --git a/moldesign/units/tools.py b/moldesign/units/tools.py index 933b43f..c0d82cb 100644 --- a/moldesign/units/tools.py +++ b/moldesign/units/tools.py @@ -168,7 +168,12 @@ 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, defunits=False, _baseunit=None): From 1a03b447af51d57fa95e3e54fa4f30f7b06729a1 Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 17 Jul 2017 14:11:50 -0700 Subject: [PATCH 31/33] Code cleanup --- moldesign/_tests/helpers.py | 2 +- moldesign/_tests/test_gaussian_math.py | 2 +- moldesign/mathutils/grids.py | 15 ++++++++------- moldesign/models/pyscf.py | 1 - moldesign/orbitals/cartesian.py | 2 +- moldesign/orbitals/spherical.py | 4 ++-- moldesign/orbitals/wfn.py | 4 +++- 7 files changed, 16 insertions(+), 14 deletions(-) diff --git a/moldesign/_tests/helpers.py b/moldesign/_tests/helpers.py index 86788b9..07fc108 100644 --- a/moldesign/_tests/helpers.py +++ b/moldesign/_tests/helpers.py @@ -243,6 +243,6 @@ def generate_grid(g, g2=None): ranges = np.ones((3, 2))*5.0*width ranges[:, 0] *= -1 ranges += ((g.center + g2.center)/2.0)[:, None] - grid = mathutils.VolumetricGrid(*ranges, 64) + grid = mathutils.VolumetricGrid(*ranges, npoints=64) allpoints = grid.allpoints() return allpoints, grid diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index 02ec734..d568f22 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -128,7 +128,7 @@ def _make_rando_linear_combination(withunits=True): moldesign.orbitals.CartesianGaussian( center=center, powers=pwr, - alpha=10.0 * random.random()/(length**2), + 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 diff --git a/moldesign/mathutils/grids.py b/moldesign/mathutils/grids.py index d1ee73c..6c260e2 100644 --- a/moldesign/mathutils/grids.py +++ b/moldesign/mathutils/grids.py @@ -32,9 +32,10 @@ class VolumetricGrid(object): 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: 25) - ypoints (int): number of grid lines in y direction (default: xpoints) - zpoints (int): number of grid lines in z direction (default: xpoints) + npoints (int): default number of grid lines in each direction (default: 32) + 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) """ dx, dy, dz = (utils.IndexView('deltas', i) for i in range(3)) xr, yr, zr = (utils.IndexView('ranges', i) for i in range(3)) @@ -42,11 +43,11 @@ class VolumetricGrid(object): xspace, yspace, zspace = (utils.IndexView('spaces', i) for i in range(3)) def __init__(self, xrange, yrange, zrange, - xpoints=32, ypoints=None, zpoints=None): + npoints=32, xpoints=None, ypoints=None, zpoints=None): - self.points = np.array([xpoints, - ypoints if ypoints is not None else xpoints, - zpoints if zpoints is not None else 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]) diff --git a/moldesign/models/pyscf.py b/moldesign/models/pyscf.py index 7d6668a..08c1113 100644 --- a/moldesign/models/pyscf.py +++ b/moldesign/models/pyscf.py @@ -452,7 +452,6 @@ def _get_ao_basis_functions(self): assert l == angular # 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, l, m, coeff=coeff[ictr], diff --git a/moldesign/orbitals/cartesian.py b/moldesign/orbitals/cartesian.py index 641bed3..d3bc709 100644 --- a/moldesign/orbitals/cartesian.py +++ b/moldesign/orbitals/cartesian.py @@ -79,7 +79,7 @@ def __init__(self, center, alpha, powers, coeff=None, normalized=True): def __repr__(self): return ("<{ndim}-D cartesian gaussian (norm: {norm:4.2f}, " "cartesian powers: {powers}, " - "width: {exp:4.2f}, " + "alpha: {exp:4.2f}, " "center: {center}>").format( ndim=self.ndim, center=self.center, exp=self.alpha, diff --git a/moldesign/orbitals/spherical.py b/moldesign/orbitals/spherical.py index b402eaa..f53a827 100644 --- a/moldesign/orbitals/spherical.py +++ b/moldesign/orbitals/spherical.py @@ -64,8 +64,8 @@ def __init__(self, center, alpha, l, m, coeff=None, normalized=True): def __repr__(self): return ("<3D Gaussian (Spherical) (coeff: {coeff:4.2f}, " - "width: {alpha:4.2f}, " - "(l,m) = {qnums}").format( + "alpha: {alpha:4.2f}, " + "(l,m) = {qnums}>").format( center=self.center, alpha=self.alpha, coeff=self.coeff, qnums=(self.l, self.m)) 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]) From a16d8a4feaf5263d3a21e81bc4217659f776f74b Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 17 Jul 2017 14:49:45 -0700 Subject: [PATCH 32/33] Fix grid initialization --- moldesign/_tests/test_mathutils.py | 6 ++---- moldesign/mathutils/grids.py | 8 ++++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/moldesign/_tests/test_mathutils.py b/moldesign/_tests/test_mathutils.py index de63414..6d13dd2 100644 --- a/moldesign/_tests/test_mathutils.py +++ b/moldesign/_tests/test_mathutils.py @@ -173,8 +173,6 @@ def test_spherical_harmonics_orthonormal(l1, m1, l2, m2): assert abs(integ) < 1e-2 - - @pytest.fixture def ndarray_ranges(): ranges = np.array([[1,4], @@ -191,7 +189,7 @@ def ranges_with_units(ndarray_ranges): @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, 3, 4, 5) + 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() @@ -213,7 +211,7 @@ def test_volumetric_grid_point_list(key, request): @pytest.mark.parametrize('key', ['ndarray_ranges', 'ranges_with_units']) def test_volumetric_iteration(key, request): ranges = request.getfixturevalue(key) - grid = mathutils.VolumetricGrid(*ranges, 4) + grid = mathutils.VolumetricGrid(*ranges, npoints=4) grid_iterator = grid.iter_points() assert (grid.xpoints, grid.ypoints, grid.zpoints) == (4,4,4) diff --git a/moldesign/mathutils/grids.py b/moldesign/mathutils/grids.py index 6c260e2..e49c1a8 100644 --- a/moldesign/mathutils/grids.py +++ b/moldesign/mathutils/grids.py @@ -32,10 +32,10 @@ class VolumetricGrid(object): 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 - npoints (int): default number of grid lines in each direction (default: 32) 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)) @@ -43,7 +43,11 @@ class VolumetricGrid(object): xspace, yspace, zspace = (utils.IndexView('spaces', i) for i in range(3)) def __init__(self, xrange, yrange, zrange, - npoints=32, xpoints=None, ypoints=None, zpoints=None): + 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, From 9c7ef1176921e2258bb6b564358610170d66206e Mon Sep 17 00:00:00 2001 From: Aaron Virshup Date: Mon, 17 Jul 2017 15:37:58 -0700 Subject: [PATCH 33/33] Add test for linear combination normalization --testall --- moldesign/_tests/test_gaussian_math.py | 19 ++++++++++++++++++- moldesign/_tests/test_wfn.py | 1 + moldesign/orbitals/primitives.py | 2 +- 3 files changed, 20 insertions(+), 2 deletions(-) diff --git a/moldesign/_tests/test_gaussian_math.py b/moldesign/_tests/test_gaussian_math.py index d568f22..db80898 100644 --- a/moldesign/_tests/test_gaussian_math.py +++ b/moldesign/_tests/test_gaussian_math.py @@ -312,7 +312,7 @@ def test_vectorized_gaussian_function_evaluations(objkey, request): _assert_almost_equal(vector_results, expected, decimal=8) -@pytest.mark.parametrize('objkey', registered_types['basis_fn']) +@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) @@ -348,6 +348,23 @@ def test_normalization(objkey, request): 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): r = g.center - coord if len(coord.shape) > 1: diff --git a/moldesign/_tests/test_wfn.py b/moldesign/_tests/test_wfn.py index f94d1e3..7239a0c 100644 --- a/moldesign/_tests/test_wfn.py +++ b/moldesign/_tests/test_wfn.py @@ -64,6 +64,7 @@ def test_basis_function_3d_grids_same_in_pyscf_and_mdt(molkey, request): @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) diff --git a/moldesign/orbitals/primitives.py b/moldesign/orbitals/primitives.py index 72eff7c..34eceb4 100644 --- a/moldesign/orbitals/primitives.py +++ b/moldesign/orbitals/primitives.py @@ -135,7 +135,7 @@ def normalize(self): """ prefactor = 1.0 / self.norm for primitive in self.primitives: - primitive *= prefactor + primitive.coeff *= prefactor def overlap(self, other): """ Calculate orbital overlap with another object