diff --git a/README.md b/README.md index ce822f43d0..5301a14461 100644 --- a/README.md +++ b/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) diff --git a/pyscf/df/df_jk.py b/pyscf/df/df_jk.py index cf114a5f3e..3ab0f389de 100644 --- a/pyscf/df/df_jk.py +++ b/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()) diff --git a/pyscf/dft/rks.py b/pyscf/dft/rks.py index 522cf89268..62182e5e94 100644 --- a/pyscf/dft/rks.py +++ b/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) diff --git a/pyscf/dft/test/test_h2o.py b/pyscf/dft/test/test_h2o.py index 78959dc4a4..34447b7361 100644 --- a/pyscf/dft/test/test_h2o.py +++ b/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() diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py index 8019a4a01d..78e8f9f4e3 100644 --- a/pyscf/gto/mole.py +++ b/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 diff --git a/pyscf/gto/test/test_mole.py b/pyscf/gto/test/test_mole.py index 1a2f77a4a0..bed6a5ec6d 100644 --- a/pyscf/gto/test/test_mole.py +++ b/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]) diff --git a/pyscf/hessian/rks.py b/pyscf/hessian/rks.py index 1ce83aefad..1f79df7a66 100644 --- a/pyscf/hessian/rks.py +++ b/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 diff --git a/pyscf/hessian/test/test_rks.py b/pyscf/hessian/test/test_rks.py new file mode 100644 index 0000000000..89864e33e7 --- /dev/null +++ b/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() + diff --git a/pyscf/hessian/test/test_uhf.py b/pyscf/hessian/test/test_uhf.py new file mode 100644 index 0000000000..3826515176 --- /dev/null +++ b/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() + diff --git a/pyscf/hessian/test/test_uks.py b/pyscf/hessian/test/test_uks.py new file mode 100644 index 0000000000..b5ecbfc63b --- /dev/null +++ b/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() + diff --git a/pyscf/hessian/uks.py b/pyscf/hessian/uks.py index 1d78cb26ba..4584d6d015 100644 --- a/pyscf/hessian/uks.py +++ b/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 diff --git a/pyscf/scf/_vhf.py b/pyscf/scf/_vhf.py index 37eec9af9a..50568b1da5 100644 --- a/pyscf/scf/_vhf.py +++ b/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) diff --git a/pyscf/scf/hf.py b/pyscf/scf/hf.py index 62197813d4..878fb736c8 100644 --- a/pyscf/scf/hf.py +++ b/pyscf/scf/hf.py @@ -883,6 +883,13 @@ def mulliken_pop(mol, dm, s=None, verbose=logger.DEBUG): .. math:: \delta_i = \sum_j M_{ij} + Returns: + A list : pop, charges + + pop : nparray + Mulliken population on each atomic orbitals + charges : nparray + Mulliken charges ''' if s is None: s = get_ovlp(mol) log = logger.new_logger(mol, verbose) @@ -930,6 +937,13 @@ def mulliken_meta(mol, dm, verbose=logger.DEBUG, | 'minao' : Project GTOs to MINAO basis | 'scf' : Fraction-averaged RHF + Returns: + A list : pop, charges + + pop : nparray + Mulliken population on each atomic orbitals + charges : nparray + Mulliken charges ''' from pyscf.lo import orth if s is None: s = get_ovlp(mol) @@ -1516,11 +1530,7 @@ def _finalize(self): def init_direct_scf(self, mol=None): if mol is None: mol = self.mol - if mol.cart: - intor = 'int2e_cart' - else: - intor = 'int2e_sph' - opt = _vhf.VHFOpt(mol, intor, 'CVHFnrs8_prescreen', + opt = _vhf.VHFOpt(mol, 'int2e', 'CVHFnrs8_prescreen', 'CVHFsetnr_direct_scf', 'CVHFsetnr_direct_scf_dm') opt.direct_scf_tol = self.direct_scf_tol