Merge pull request #304 from sunqm/master

Handle Cartesian GTOs in RSH integrals

README.md

@@ -6,7 +6,7 @@ Python-based Simulations of Chemistry Framework
 ===============================================
 [![Build Status](https://travis-ci.org/pyscf/pyscf.svg?branch=master)](https://travis-ci.org/pyscf/pyscf)
 
-2018-12-31
+2019-03-15
 
 * [Stable release 1.6.1](https://github.com/pyscf/pyscf/releases/tag/v1.6.1)
 * [1.7 alpha](https://github.com/pyscf/pyscf/tree/dev)

pyscf/df/df_jk.py

@@ -148,18 +148,13 @@ def _cderi(self, x):
         def nuc_grad_method(self):
             raise NotImplementedError
 
-    for k, v in _METHODS_TO_ASSIGN.items():
-        setattr(DFHF, k, v)
     return DFHF(mf)
 
-# A tag to label the derived SCF class
+# 1. A tag to label the derived SCF class
+# 2. A hook to register DF specific methods, such as nuc_grad_method.
 class _DFHF:
     pass
 
-# Some methods should be attached to the DF-SCF class in density_fit function
-# above. They can be registered in _METHODS_TO_ASSIGN.
-_METHODS_TO_ASSIGN = {}
-
 
 def get_jk(dfobj, dm, hermi=1, vhfopt=None, with_j=True, with_k=True):
     t0 = t1 = (time.clock(), time.time())

pyscf/dft/rks.py

@@ -154,7 +154,8 @@ def _get_k_lr(mol, dm, omega=0, hermi=0):
     with mol.with_range_coulomb(omega):
         # Compute the long range part of ERIs temporarily with omega. Restore
         # the original omega when the block ends
-        vklr = jk.get_jk(mol, dms, ['ijkl,jk->il']*len(dms))
+        intor = mol._add_suffix('int2e')
+        vklr = jk.get_jk(mol, dms, ['ijkl,jk->il']*len(dms), intor=intor)
     return numpy.asarray(vklr).reshape(dm.shape)
 
 

pyscf/dft/test/test_h2o.py

@@ -389,6 +389,16 @@ def test_nr_rks_rsh(self):
         vxc = method.get_veff(h2o, dm, dm, vxc)
         self.assertAlmostEqual(lib.finger(vxc), 23.058813088809824, 8)
 
+    def test_nr_rks_rsh_cart(self):
+        mol1 = h2o.copy()
+        mol1.basis = 'ccpvdz'
+        mol1.cart = True
+        mol1.build(0, 0)
+        method = dft.RKS(mol1)
+        method.xc = 'B97M_V'
+        method.grids.atom_grid = {"H": (50, 194), "O": (50, 194),}
+        self.assertAlmostEqual(method.kernel(), -76.44022393692919, 9)
+
     def test_nr_uks_rsh(self):
         method = dft.UKS(h2o)
         dm = method.get_init_guess()

pyscf/gto/mole.py

@@ -2463,8 +2463,8 @@ def set_rinv_zeta(self, zeta):
     set_rinv_zeta_ = set_rinv_zeta  # for backward compatibility
 
     def with_rinv_zeta(self, zeta):
-        '''Retuen a temporary mol context which has the rquired charge
-        distribution on the "rinv_origin": rho(r) = Norm * exp(-zeta * r^2).
+        '''Retuen a temporary mol context which has the rquired Gaussian charge
+        distribution placed at "rinv_origin": rho(r) = Norm * exp(-zeta * r^2).
         See also :func:`mol.set_rinv_zeta`
 
         Examples:
@@ -2476,8 +2476,9 @@ def with_rinv_zeta(self, zeta):
         return _TemporaryMoleContext(self.set_rinv_zeta, (zeta,), (zeta0,))
 
     def with_rinv_as_nucleus(self, atm_id):
-        '''Retuen a temporary mol context which has the rquired origin of 1/r
-        operator and the required nuclear charge distribution on 1/r.
+        '''Retuen a temporary mol context in which the rinv operator (1/r) is
+        treated like the Coulomb potential of a Gaussian charge distribution
+        rho(r) = Norm * exp(-zeta * r^2) at the place of the input atm_id.
 
         Examples:
 
@@ -2592,6 +2593,11 @@ def atom_pure_symbol(self, atm_id):
         '''
         return _std_symbol(self._atom[atm_id][0])
 
+    @property
+    def elements(self):
+        '''A list of elements in the molecule'''
+        return [self.atom_pure_symbol(i) for i in range(self.natm)]
+
     def atom_charge(self, atm_id):
         r'''Nuclear effective charge of the given atom id
         Note "atom_charge /= charge(atom_symbol)" when ECP is enabled.
@@ -2618,7 +2624,7 @@ def atom_nelec_core(self, atm_id):
         '''
         return charge(self.atom_symbol(atm_id)) - self.atom_charge(atm_id)
 
-    def atom_coord(self, atm_id):
+    def atom_coord(self, atm_id, unit='Bohr'):
         r'''Coordinates (ndarray) of the given atom id
 
         Args:
@@ -2632,12 +2638,18 @@ def atom_coord(self, atm_id):
         [ 0.          0.          2.07869874]
         '''
         ptr = self._atm[atm_id,PTR_COORD]
-        return self._env[ptr:ptr+3]
+        if unit[:3].upper() == 'ANG':
+            return self._env[ptr:ptr+3] * param.BOHR
+        else:
+            return self._env[ptr:ptr+3]
 
-    def atom_coords(self):
+    def atom_coords(self, unit='Bohr'):
         '''np.asarray([mol.atom_coords(i) for i in range(mol.natm)])'''
         ptr = self._atm[:,PTR_COORD]
-        return self._env[numpy.vstack((ptr,ptr+1,ptr+2)).T]
+        c = self._env[numpy.vstack((ptr,ptr+1,ptr+2)).T]
+        if unit[:3].upper() == 'ANG':
+            c *= param.BOHR
+        return c
 
     atom_mass_list = atom_mass_list
 

pyscf/gto/test/test_mole.py

@@ -558,6 +558,7 @@ def test_atom_method(self):
             self.assertEqual(len(shls), nshls)
             self.assertEqual(mol0.atom_nshells(i), nshls)
         aoslice = mol0.aoslice_2c_by_atom()
+        mol0.elements
         self.assertEqual([x[2] for x in aoslice], [0, 8, 56])
         self.assertEqual([x[3] for x in aoslice], [8, 56, 64])
 

pyscf/hessian/rks.py

@@ -489,6 +489,7 @@ class Hessian(rhf_hess.Hessian):
     def __init__(self, mf):
         rhf_hess.Hessian.__init__(self, mf)
         self.grids = None
+        self.grid_response = False
         self._keys = self._keys.union(['grids'])
 
     partial_hess_elec = partial_hess_elec

pyscf/hessian/test/test_rks.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
+#
+# 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 unittest
+import numpy
+from pyscf import gto, dft, lib
+from pyscf import grad, hessian
+
+mol = gto.Mole()
+mol.verbose = 5
+mol.output = '/dev/null'
+mol.atom.extend([
+    ["O" , (0. , 0.     , 0.)],
+    [1   , (0. , -0.757 , 0.587)],
+    [1   , (0. , 0.757  , 0.587)] ])
+mol.basis = '6-31g'
+mol.build()
+
+def tearDownModule():
+    global mol
+    del mol
+
+class KnownValues(unittest.TestCase):
+    def test_finite_diff_lda_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.7828771346902333, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    def test_finite_diff_b3lyp_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.7590878171493624, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    def test_finite_diff_wb97x_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'wb97x'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.7637876979690904, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 2)
+
+
+if __name__ == "__main__":
+    print("Full Tests for RKS Hessian")
+    unittest.main()
+

pyscf/hessian/test/test_uhf.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
+#
+# 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 unittest
+import numpy
+from pyscf import gto, scf, lib
+from pyscf import grad, hessian
+
+mol = gto.Mole()
+mol.verbose = 5
+mol.output = '/dev/null'
+mol.atom.extend([
+    ["O" , (0. , 0.     , 0.)],
+    [1   , (0. , -0.757 , 0.587)],
+    [1   , (0. , 0.757  , 0.587)] ])
+mol.basis = '6-31g'
+mol.spin = 2
+mol.build()
+
+def tearDownModule():
+    global mol
+    del mol
+
+class KnownValues(unittest.TestCase):
+    def test_finite_diff_rhf_hess(self):
+        mf = scf.UHF(mol)
+        mf.conv_tol = 1e-14
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.20243405976628576, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 5)
+
+
+if __name__ == "__main__":
+    print("Full Tests for UHF Hessian")
+    unittest.main()
+

pyscf/hessian/test/test_uks.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+# Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
+#
+# 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 unittest
+import numpy
+from pyscf import gto, dft, lib
+from pyscf import grad, hessian
+
+mol = gto.Mole()
+mol.verbose = 5
+mol.output = '/dev/null'
+mol.atom.extend([
+    ["O" , (0. , 0.     , 0.)],
+    [1   , (0. , -0.757 , 0.587)],
+    [1   , (0. , 0.757  , 0.587)] ])
+mol.basis = '6-31g'
+mol.charge = 1
+mol.spin = 1
+mol.build()
+
+def tearDownModule():
+    global mol
+    del mol
+
+class KnownValues(unittest.TestCase):
+    def test_finite_diff_lda_hess(self):
+        mf = dft.UKS(mol)
+        mf.conv_tol = 1e-14
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.8503072107510495, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    def test_finite_diff_b3lyp_hess(self):
+        mf = dft.UKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.8208641727673912, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    def test_finite_diff_wb97x_hess(self):
+        mf = dft.UKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'wb97x'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.finger(hess), -0.8207572641132195, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 2)
+
+
+if __name__ == "__main__":
+    print("Full Tests for UKS Hessian")
+    unittest.main()
+

pyscf/hessian/uks.py

@@ -575,6 +575,7 @@ class Hessian(uhf_hess.Hessian):
     def __init__(self, mf):
         uhf_hess.Hessian.__init__(self, mf)
         self.grids = None
+        self.grid_response = False
         self._keys = self._keys.union(['grids'])
 
     partial_hess_elec = partial_hess_elec

pyscf/scf/_vhf.py

@@ -280,7 +280,10 @@ def direct_mapdm(intor, aosym, jkdescript,
         vshape = (n_dm,ncomp) + get_dims(vsym[-2:], shls_slice, ao_loc)
         vjk.append(numpy.empty(vshape))
         for j in range(n_dm):
-            assert(dms[j].shape == get_dims(dmsym, shls_slice, ao_loc))
+            if dms[j].shape != get_dims(dmsym, shls_slice, ao_loc):
+                raise RuntimeError('dm[%d] shape %s is inconsistent with the '
+                                   'shls_slice shape %s' %
+                                   (j, dms[j].shape, get_dims(dmsym, shls_slice, ao_loc)))
             dmsptr[i*n_dm+j] = dms[j].ctypes.data_as(ctypes.c_void_p)
             vjkptr[i*n_dm+j] = vjk[i][j].ctypes.data_as(ctypes.c_void_p)
             fjk[i*n_dm+j] = f1
@@ -356,7 +359,10 @@ def direct_bindm(intor, aosym, jkdescript,
     for i, (dmsym, vsym) in enumerate(descr_sym):
         f1 = _fpointer('CVHFnr%s_%s_%s'%(aosym, dmsym, vsym))
 
-        assert(dms[i].shape == get_dims(dmsym, shls_slice, ao_loc))
+        if dms[i].shape != get_dims(dmsym, shls_slice, ao_loc):
+            raise RuntimeError('dm[%d] shape %s is inconsistent with the '
+                               'shls_slice shape %s' %
+                               (i, dms[i].shape, get_dims(dmsym, shls_slice, ao_loc)))
         vshape = (ncomp,) + get_dims(vsym[-2:], shls_slice, ao_loc)
         vjk.append(numpy.empty(vshape))
         dmsptr[i] = dms[i].ctypes.data_as(ctypes.c_void_p)