From bb89e51cb8d208ad223e69c9f589f4166f9dcf38 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Tue, 19 Sep 2023 22:09:23 -0700 Subject: [PATCH] Debugging --- .github/workflows/ci_linux/python_deps.sh | 16 +- pyscf/lib/linalg_helper.py | 1 + pyscf/scf/_response_functions.py | 1 + pyscf/solvent/test/test_ddcosmo.py | 507 -------- pyscf/solvent/test/test_ddcosmo_grad.py | 1033 ----------------- pyscf/solvent/test/test_ddcosmo_tdscf_grad.py | 14 +- pyscf/solvent/test/test_ddpcm.py | 40 - pyscf/solvent/test/test_pol_embed.py | 322 ----- pyscf/tdscf/rhf.py | 2 + 9 files changed, 19 insertions(+), 1917 deletions(-) delete mode 100644 pyscf/solvent/test/test_ddcosmo.py delete mode 100644 pyscf/solvent/test/test_ddcosmo_grad.py delete mode 100644 pyscf/solvent/test/test_ddpcm.py delete mode 100644 pyscf/solvent/test/test_pol_embed.py diff --git a/.github/workflows/ci_linux/python_deps.sh b/.github/workflows/ci_linux/python_deps.sh index 1511ded008..90938a1015 100755 --- a/.github/workflows/ci_linux/python_deps.sh +++ b/.github/workflows/ci_linux/python_deps.sh @@ -1,11 +1,11 @@ #!/usr/bin/env bash python -m pip install --upgrade pip pip install "numpy!=1.16,!=1.17" "scipy!=1.5" h5py pytest pytest-cov pytest-timer -pip install pyberny geometric -pip install spglib - -#cppe -version=$(python -c 'import sys; version=sys.version_info[:2]; print("{0}.{1}".format(*version))') -if [ $version != '2.7' ] && [ $version != '3.5' ]; then - pip install cppe -fi +#pip install pyberny geometric +#pip install spglib +# +##cppe +#version=$(python -c 'import sys; version=sys.version_info[:2]; print("{0}.{1}".format(*version))') +#if [ $version != '2.7' ] && [ $version != '3.5' ]; then +# pip install cppe +#fi diff --git a/pyscf/lib/linalg_helper.py b/pyscf/lib/linalg_helper.py index fe2892364c..8df836b95e 100644 --- a/pyscf/lib/linalg_helper.py +++ b/pyscf/lib/linalg_helper.py @@ -465,6 +465,7 @@ def davidson1(aop, x0, precond, tol=1e-12, max_cycle=50, max_space=12, fill_heff(heff, xs, ax, xt, axt, dot) xt = axt = None w, v = scipy.linalg.eigh(heff[:space,:space]) + print('linalg.v.dtype=', v.dtype) if callable(pick): w, v, idx = pick(w, v, nroots, locals()) if len(w) == 0: diff --git a/pyscf/scf/_response_functions.py b/pyscf/scf/_response_functions.py index 213770e2d6..6b2f5e7e19 100644 --- a/pyscf/scf/_response_functions.py +++ b/pyscf/scf/_response_functions.py @@ -79,6 +79,7 @@ def vind(dm1): if hermi == 2: v1 = numpy.zeros_like(dm1) else: + print(hermi, 'hermi') v1 = ni.nr_rks_fxc(mol, mf.grids, mf.xc, dm0, dm1, 0, hermi, rho0, vxc, fxc, max_memory=max_memory) if hybrid: diff --git a/pyscf/solvent/test/test_ddcosmo.py b/pyscf/solvent/test/test_ddcosmo.py deleted file mode 100644 index be13e22f20..0000000000 --- a/pyscf/solvent/test/test_ddcosmo.py +++ /dev/null @@ -1,507 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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 -from functools import reduce -import numpy -import scipy.special -from pyscf import lib, gto, scf, dft, ao2mo, df -from pyscf.solvent import ddcosmo -from pyscf.solvent import _attach_solvent -from pyscf.symm import sph - - -def make_v_phi(mol, dm, r_vdw, lebedev_order): - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - natm = mol.natm - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(lebedev_order) - - pmol = mol.copy() - v_phi = [] - for ia in range(natm): - for i,c in enumerate(coords_1sph): - r = atom_coords[ia] + r_vdw[ia] * c - dr = atom_coords - r - v_nuc = (atom_charges / numpy.linalg.norm(dr, axis=1)).sum() - pmol.set_rinv_orig(r) - v_e = numpy.einsum('ij,ji', pmol.intor('int1e_rinv'), dm) - v_phi.append(v_nuc - v_e) - v_phi = numpy.array(v_phi).reshape(natm,-1) - return v_phi - -def make_L(pcmobj, r_vdw, lebedev_order, lmax, eta=0.1): - mol = pcmobj.mol - natm = mol.natm - nlm = (lmax+1)**2 - - leb_coords, leb_weights = ddcosmo.make_grids_one_sphere(lebedev_order) - nleb_grid = leb_weights.size - atom_coords = mol.atom_coords() - Ylm_sphere = numpy.vstack(sph.real_sph_vec(leb_coords, lmax, True)) - fi = ddcosmo.make_fi(pcmobj, r_vdw) - - L_diag = numpy.zeros((natm,nlm)) - p1 = 0 - for l in range(lmax+1): - p0, p1 = p1, p1 + (l*2+1) - L_diag[:,p0:p1] = 4*numpy.pi/(l*2+1) - L_diag /= r_vdw.reshape(-1,1) - L = numpy.diag(L_diag.ravel()).reshape(natm,nlm,natm,nlm) - for ja in range(natm): - for ka in range(natm): - if ja == ka: - continue - vjk = r_vdw[ja] * leb_coords + atom_coords[ja] - atom_coords[ka] - v = lib.norm(vjk, axis=1) - tjk = v / r_vdw[ka] - sjk = vjk / v.reshape(-1,1) - Ys = sph.real_sph_vec(sjk, lmax, True) - # scale the weight, see JCTC 9, 3637, Eq (16) - wjk = pcmobj.regularize_xt(tjk, eta) - wjk[fi[ja]>1] /= fi[ja,fi[ja]>1] - tt = numpy.ones_like(wjk) - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) / r_vdw[ka] - p0, p1 = p1, p1 + (l*2+1) - val = numpy.einsum('n,xn,n,mn->xm', leb_weights, Ylm_sphere, wjk*tt, Ys[l]) - L[ja,:,ka,p0:p1] += -fac * val - tt *= tjk - return L.reshape(natm*nlm,natm*nlm) - -def make_psi(mol, dm, r_vdw, lmax): - grids = ddcosmo.Grids(mol) - atom_grids_tab = grids.gen_atomic_grids(mol) - grids.build() - - ao = dft.numint.eval_ao(mol, grids.coords) - den = dft.numint.eval_rho(mol, ao, dm) - den *= grids.weights - natm = mol.natm - nlm = (lmax+1)**2 - psi = numpy.empty((natm,nlm)) - i1 = 0 - for ia in range(natm): - xnj, w = atom_grids_tab[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + w.size - r = lib.norm(xnj, axis=1) - snj = xnj/r.reshape(-1,1) - Ys = sph.real_sph_vec(snj, lmax, True) - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - rr = numpy.zeros_like(r) - rr[r<=r_vdw[ia]] = r[r<=r_vdw[ia]]**l / r_vdw[ia]**(l+1) - rr[r> r_vdw[ia]] = r_vdw[ia]**l / r[r>r_vdw[ia]]**(l+1) - psi[ia,p0:p1] = -fac * numpy.einsum('n,n,mn->m', den[i0:i1], rr, Ys[l]) - psi[ia,0] += numpy.sqrt(4*numpy.pi)/r_vdw[ia] * mol.atom_charge(ia) - return psi - -def make_vmat(pcm, r_vdw, lebedev_order, lmax, LX, LS): - mol = pcm.mol - grids = ddcosmo.Grids(mol) - atom_grids_tab = grids.gen_atomic_grids(mol) - grids.build() - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(lebedev_order) - ao = dft.numint.eval_ao(mol, grids.coords) - nao = ao.shape[1] - vmat = numpy.zeros((nao,nao)) - i1 = 0 - for ia in range(mol.natm): - xnj, w = atom_grids_tab[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + w.size - r = lib.norm(xnj, axis=1) - Ys = sph.real_sph_vec(xnj/r.reshape(-1,1), lmax, True) - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - rr = numpy.zeros_like(r) - rr[r<=r_vdw[ia]] = r[r<=r_vdw[ia]]**l / r_vdw[ia]**(l+1) - rr[r> r_vdw[ia]] = r_vdw[ia]**l / r[r>r_vdw[ia]]**(l+1) - eta_nj = fac * numpy.einsum('n,mn,m->n', rr, Ys[l], LX[ia,p0:p1]) - vmat -= numpy.einsum('n,np,nq->pq', grids.weights[i0:i1] * eta_nj, - ao[i0:i1], ao[i0:i1]) - - atom_coords = mol.atom_coords() - Ylm_sphere = numpy.vstack(sph.real_sph_vec(coords_1sph, lmax, True)) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - xi_nj = numpy.einsum('n,jn,xn,jx->jn', weights_1sph, ui, Ylm_sphere, LS) - - pmol = mol.copy() - for ia in range(mol.natm): - for i,c in enumerate(coords_1sph): - r = atom_coords[ia] + r_vdw[ia] * c - pmol.set_rinv_orig(r) - vmat += pmol.intor('int1e_rinv') * xi_nj[ia,i] - return vmat - - -def make_B(pcmobj, r_vdw, ui, ylm_1sph, cached_pol, L): - mol = pcmobj.mol - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - mol = pcmobj.mol - natm = mol.natm - nao = mol.nao - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - grids = pcmobj.grids - - extern_point_idx = ui > 0 - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - - max_memory = pcmobj.max_memory - lib.current_memory()[0] - blksize = int(max(max_memory*.9e6/8/nao**2, 400)) - - cav_coords = cav_coords[extern_point_idx] - int3c2e = mol._add_suffix('int3c2e') - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, - mol._env, int3c2e) - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s2ij', cintopt=cintopt) - nao_pair = v_nj.shape[0] - v_phi = numpy.zeros((nao_pair, natm, ngrid_1sph)) - v_phi[:,extern_point_idx] += v_nj - - phi = numpy.einsum('n,xn,jn,ijn->ijx', weights_1sph, ylm_1sph, ui, v_phi) - - Xvec = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi.reshape(-1,natm*nlm).T) - Xvec = Xvec.reshape(natm,nlm,nao_pair) - - ao = mol.eval_gto('GTOval', grids.coords) - aow = numpy.einsum('gi,g->gi', ao, grids.weights) - aopair = lib.pack_tril(numpy.einsum('gi,gj->gij', ao, aow)) - - psi = numpy.zeros((nao_pair, natm, nlm)) - i1 = 0 - for ia in range(natm): - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + fak_pol[0].shape[1] - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - psi[:,ia,p0:p1] = -fac * numpy.einsum('mn,ni->im', fak_pol[l], aopair[i0:i1]) - - B = lib.einsum('pnl,nlq->pq', psi, Xvec) - B = B + B.T - B = ao2mo.restore(1, B, nao) - return B - - -def setUpModule(): - global mol - mol = gto.Mole() - mol.atom = ''' O 0.00000000 0.00000000 -0.11081188 - H -0.00000000 -0.84695236 0.59109389 - H -0.00000000 0.89830571 0.52404783 ''' - mol.basis = '3-21g' - mol.verbose = 5 - mol.output = '/dev/null' - mol.build() - -def tearDownModule(): - global mol - mol.stdout.close() - del mol - -class KnownValues(unittest.TestCase): - def test_ddcosmo_scf(self): - mol = gto.M(atom=''' H 0 0 0 ''', charge=1, basis='sto3g', verbose=7, - output='/dev/null') - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 10 - pcm.lebedev_order = 29 - mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol), pcm) - mf.init_guess = '1e' - mf.run() - self.assertAlmostEqual(mf.e_tot, -0.1645636146393864, 9) - - mol = gto.M(atom=''' - 6 0.000000 0.000000 -0.542500 - 8 0.000000 0.000000 0.677500 - 1 0.000000 0.935307 -1.082500 - 1 0.000000 -0.935307 -1.082500 - ''', basis='sto3g', verbose=7, - output='/dev/null') - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 6 - pcm.lebedev_order = 17 - mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol), pcm).run() - self.assertAlmostEqual(mf.e_tot, -112.35463433688, 9) - - def test_ddcosmo_scf_with_overwritten_attributes(self): - mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol)) - mf.kernel() - self.assertAlmostEqual(mf.e_tot, -75.57006258287, 9) - - mf.with_solvent.lebedev_order = 15 - mf.with_solvent.lmax = 5 - mf.with_solvent.eps = .5 - mf.with_solvent.conv_tol = 1e-8 - mf.kernel() - self.assertAlmostEqual(mf.e_tot, -75.55351392557, 9) - - mf.with_solvent.grids.radi_method = dft.mura_knowles - mf.with_solvent.grids.atom_grid = {"H": (8, 50), "O": (8, 50),} - mf.kernel() - self.assertAlmostEqual(mf.e_tot, -75.55237426980, 9) - - def test_make_ylm(self): - numpy.random.seed(1) - lmax = 6 - r = numpy.random.random((100,3)) - numpy.ones(3)*.5 - r = r / lib.norm(r,axis=1).reshape(-1,1) - - ngrid = r.shape[0] - cosphi = r[:,2] - sinphi = (1-cosphi**2)**.5 - costheta = numpy.ones(ngrid) - sintheta = numpy.zeros(ngrid) - costheta[sinphi!=0] = r[sinphi!=0,0] / sinphi[sinphi!=0] - sintheta[sinphi!=0] = r[sinphi!=0,1] / sinphi[sinphi!=0] - costheta[costheta> 1] = 1 - costheta[costheta<-1] =-1 - sintheta[sintheta> 1] = 1 - sintheta[sintheta<-1] =-1 - varphi = numpy.arccos(cosphi) - theta = numpy.arccos(costheta) - theta[sintheta<0] = 2*numpy.pi - theta[sintheta<0] - ylmref = [] - for l in range(lmax+1): - ylm = numpy.empty((l*2+1,ngrid)) - ylm[l] = scipy.special.sph_harm(0, l, theta, varphi).real - for m in range(1, l+1): - f1 = scipy.special.sph_harm(-m, l, theta, varphi) - f2 = scipy.special.sph_harm( m, l, theta, varphi) - # complex to real spherical functions - if m % 2 == 1: - ylm[l-m] = (-f1.imag - f2.imag) / numpy.sqrt(2) - ylm[l+m] = ( f1.real - f2.real) / numpy.sqrt(2) - else: - ylm[l-m] = (-f1.imag + f2.imag) / numpy.sqrt(2) - ylm[l+m] = ( f1.real + f2.real) / numpy.sqrt(2) - if l == 1: - ylm = ylm[[2,0,1]] - ylmref.append(ylm) - ylmref = numpy.vstack(ylmref) - ylm = numpy.vstack(sph.real_sph_vec(r, lmax, True)) - self.assertTrue(abs(ylmref - ylm).max() < 1e-14) - - def test_L_x(self): - pcm = ddcosmo.DDCOSMO(mol) - r_vdw = ddcosmo.get_atomic_radii(pcm) - n = mol.natm * (pcm.lmax+1)**2 - Lref = make_L(pcm, r_vdw, pcm.lebedev_order, pcm.lmax, pcm.eta).reshape(n,n) - - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - fi = ddcosmo.make_fi(pcm, r_vdw) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi).reshape(n,n) - - numpy.random.seed(1) - x = numpy.random.random(n) - self.assertTrue(abs(Lref.dot(n)-L.dot(n)).max() < 1e-12) - - def test_phi(self): - pcm = ddcosmo.DDCOSMO(mol) - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - - numpy.random.seed(1) - nao = mol.nao_nr() - dm = numpy.random.random((nao,nao)) - dm = dm + dm.T - - v_phi = make_v_phi(mol, dm, r_vdw, pcm.lebedev_order) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - phi = -numpy.einsum('n,xn,jn,jn->jx', weights_1sph, ylm_1sph, ui, v_phi) - phi1 = ddcosmo.make_phi(pcm, dm, r_vdw, ui, ylm_1sph) - self.assertTrue(abs(phi - phi1).max() < 1e-12) - - def test_psi_vmat(self): - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 2 - pcm.eps = 0 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - grids = ddcosmo.Grids(mol).build() - pcm.grids = grids - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - - numpy.random.seed(1) - nao = mol.nao_nr() - dm = numpy.random.random((nao,nao)) - dm = dm + dm.T - natm = mol.natm - nlm = (pcm.lmax+1)**2 - LX = numpy.random.random((natm,nlm)) - - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - psi, vmat = ddcosmo.make_psi_vmat(pcm, dm, r_vdw, ui, - ylm_1sph, cached_pol, LX, L)[:2] - psi_ref = make_psi(pcm.mol, dm, r_vdw, pcm.lmax) - self.assertAlmostEqual(abs(psi_ref - psi).max(), 0, 12) - - LS = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, - psi_ref.ravel()).reshape(natm,nlm) - vmat_ref = make_vmat(pcm, r_vdw, pcm.lebedev_order, pcm.lmax, LX, LS) - self.assertAlmostEqual(abs(vmat_ref - vmat).max(), 0, 12) - - def test_B_dot_x(self): - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 2 - pcm.eps = 0 - natm = mol.natm - nao = mol.nao - nlm = (pcm.lmax+1)**2 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - grids = ddcosmo.Grids(mol).run(level=0) - pcm.grids = grids - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - B = make_B(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - - numpy.random.seed(19) - dm = numpy.random.random((2,nao,nao)) - Bx = numpy.einsum('ijkl,xkl->xij', B, dm) - - phi = ddcosmo.make_phi(pcm, dm, r_vdw, ui, ylm_1sph, with_nuc=False) - Xvec = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi.reshape(-1,natm*nlm).T) - Xvec = Xvec.reshape(natm,nlm,-1).transpose(2,0,1) - psi, vref, LS = ddcosmo.make_psi_vmat(pcm, dm, r_vdw, ui, ylm_1sph, - cached_pol, Xvec, L, with_nuc=False) - self.assertAlmostEqual(abs(Bx - vref).max(), 0, 12) - e1 = numpy.einsum('nij,nij->n', psi, Xvec) - e2 = numpy.einsum('nij,nij->n', phi, LS) - e3 = numpy.einsum('nij,nij->n', dm, vref) * .5 - self.assertAlmostEqual(abs(e1-e2).max(), 0, 12) - self.assertAlmostEqual(abs(e1-e3).max(), 0, 12) - - vmat = pcm._B_dot_x(dm) - self.assertEqual(vmat.shape, (2,nao,nao)) - self.assertAlmostEqual(abs(vmat-vref*.5).max(), 0, 12) - - def test_vmat(self): - mol = gto.M(atom='H 0 0 0; H 0 1 1.2; H 1. .1 0; H .5 .5 1', verbose=0) - pcmobj = ddcosmo.DDCOSMO(mol) - f = pcmobj.as_solver() - nao = mol.nao_nr() - numpy.random.seed(1) - dm1 = numpy.random.random((nao,nao)) - dm1 = dm1 + dm1.T - e0, vmat0 = f(dm1) - dx = 0.0001 - vmat1 = numpy.zeros_like(dm1) - for i in range(nao): - for j in range(i): - dm1[i,j] += dx - dm1[j,i] += dx - e1 = f(dm1)[0] - vmat1[i,j] = vmat1[j,i] = (e1 - e0) / (dx*2) - dm1[i,j] -= dx - dm1[j,i] -= dx - dm1[i,i] += dx - e1 = f(dm1)[0] - vmat1[i,i] = (e1 - e0) / dx - dm1[i,i] -= dx - self.assertAlmostEqual(abs(vmat0-vmat1).max(), 0, 4) - - def test_as_scanner(self): - mol = gto.M(atom=''' - 6 0.000000 0.000000 -0.542500 - 8 0.000000 0.000000 0.677500 - 1 0.000000 0.935307 -1.082500 - 1 0.000000 -0.935307 -1.082500 - ''', basis='sto3g', verbose=7, - output='/dev/null') - mf_scanner = ddcosmo.ddcosmo_for_scf(scf.RHF(mol)).as_scanner() - mf_scanner(mol) - self.assertEqual(mf_scanner.with_solvent.grids.coords.shape, (48212, 3)) - mf_scanner('H 0. 0. 0.; H 0. 0. .9') - self.assertEqual(mf_scanner.with_solvent.grids.coords.shape, (20048, 3)) - - h2 = gto.M(atom='H 0. 0. 0.; H 0. 0. .9', basis='sto3g', verbose=7, - output='/dev/null') - mf_h2 = ddcosmo.ddcosmo_for_scf(scf.RHF(h2)).run() - self.assertAlmostEqual(mf_h2.e_tot, mf_scanner.e_tot, 9) - - def test_newton_rohf(self): - mf = mol.ROHF(max_memory=0).ddCOSMO() - mf = mf.newton() - e = mf.kernel() - self.assertAlmostEqual(e, -75.57006258287, 9) - - mf = mol.RHF().ddCOSMO() - e = mf.kernel() - self.assertAlmostEqual(e, -75.57006258287, 9) - - def test_convert_scf(self): - mf = mol.RHF().ddCOSMO() - mf = mf.to_uhf() - self.assertTrue(isinstance(mf, scf.uhf.UHF)) - self.assertTrue(isinstance(mf, _attach_solvent._Solvation)) - - def test_reset(self): - mol1 = gto.M(atom='H 0 0 0; H 0 0 .9', basis='cc-pvdz') - mf = scf.RHF(mol).density_fit().ddCOSMO().newton() - mf.reset(mol1) - self.assertTrue(mf.mol is mol1) - self.assertTrue(mf.with_df.mol is mol1) - self.assertTrue(mf.with_solvent.mol is mol1) - self.assertTrue(mf._scf.with_df.mol is mol1) - self.assertTrue(mf._scf.with_solvent.mol is mol1) - - def test_rhf_tda(self): - # TDA with equilibrium_solvation - mf = mol.RHF().ddCOSMO().run(conv_tol=1e-10) - td = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) - ref = numpy.array([0.30124900879, 0.358722766464, 0.3950184783571]) - self.assertAlmostEqual(abs(ref - td.e).max(), 0, 7) - - - # TDA without equilibrium_solvation - mf = mol.RHF().ddCOSMO().run(conv_tol=1e-10) - td = mf.TDA().ddCOSMO().run() - ref = numpy.array([0.301421953639, 0.358782851661, 0.400409174628]) - self.assertAlmostEqual(abs(ref - td.e).max(), 0, 7) - -# TODO: add tests for direct-scf, ROHF, ROKS, .newton(), and their mixes - - -if __name__ == "__main__": - print("Full Tests for ddcosmo") - unittest.main() diff --git a/pyscf/solvent/test/test_ddcosmo_grad.py b/pyscf/solvent/test/test_ddcosmo_grad.py deleted file mode 100644 index 08ccc9fd8d..0000000000 --- a/pyscf/solvent/test/test_ddcosmo_grad.py +++ /dev/null @@ -1,1033 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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 -from functools import reduce -import numpy -from pyscf import lib -from pyscf import gto -from pyscf import scf -from pyscf import ao2mo -from pyscf import mcscf -from pyscf import cc -from pyscf import tdscf -from pyscf import dft -from pyscf import df -from pyscf import solvent -from pyscf.scf import cphf -from pyscf.grad import rhf as rhf_grad -from pyscf.grad import rks as rks_grad -from pyscf.solvent import ddcosmo -from pyscf.solvent import ddcosmo_grad -from pyscf.solvent import _ddcosmo_tdscf_grad -from pyscf.symm import sph - -def tda_grad(td, z): - '''ddcosmo TDA gradients''' - mol = td.mol - mf = td._scf - mo_coeff = mf.mo_coeff - mo_energy = mf.mo_energy - mo_occ = mf.mo_occ - nao, nmo = mo_coeff.shape - nocc = (mo_occ>0).sum() - nvir = nmo - nocc - z = z[0].reshape(nocc,nvir).T * numpy.sqrt(2) - orbv = mo_coeff[:,nocc:] - orbo = mo_coeff[:,:nocc] - - r_vdw = ddcosmo.get_atomic_radii(td.with_solvent) - fi = ddcosmo.make_fi(td.with_solvent, r_vdw) - ui = 1 - fi - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(td.with_solvent.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, td.with_solvent.lmax, True)) - grids = td.with_solvent.grids - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, td.with_solvent.lmax) - L = ddcosmo.make_L(td.with_solvent, r_vdw, ylm_1sph, fi) - - def fvind(x): - v_mo = numpy.einsum('iabj,xai->xbj', g[:nocc,nocc:,nocc:,:nocc], x) - v_mo += numpy.einsum('aibj,xai->xbj', g[nocc:,:nocc,nocc:,:nocc], x) - return v_mo - - h1 = rhf_grad.get_hcore(mol) - s1 = rhf_grad.get_ovlp(mol) - - eri1 = -mol.intor('int2e_ip1', aosym='s1', comp=3) - eri1 = eri1.reshape(3,nao,nao,nao,nao) - eri0 = ao2mo.kernel(mol, mo_coeff) - eri0 = ao2mo.restore(1, eri0, nmo).reshape(nmo,nmo,nmo,nmo) - g = eri0 * 2 - eri0.transpose(0,3,2,1) - zeta = lib.direct_sum('i+j->ij', mo_energy, mo_energy) * .5 - zeta[nocc:,:nocc] = mo_energy[:nocc] - zeta[:nocc,nocc:] = mo_energy[nocc:] - - dielectric = td.with_solvent.eps - if dielectric > 0: - f_epsilon = (dielectric-1.)/dielectric - else: - f_epsilon = 1 - pcm_nuc = .5 * f_epsilon * nuc_part1(td.with_solvent, r_vdw, ui, ylm_1sph, cached_pol, L) - B0 = .5 * f_epsilon * make_B(td.with_solvent, r_vdw, ui, ylm_1sph, cached_pol, L) - B0 = lib.einsum('pqrs,pi,qj,rk,sl->ijkl', B0, mo_coeff, mo_coeff, mo_coeff, mo_coeff) - g += B0 * 2 - B1 = .5 * f_epsilon * make_B1(td.with_solvent, r_vdw, ui, ylm_1sph, cached_pol, L) - - offsetdic = mol.offset_nr_by_atom() - de = numpy.zeros((mol.natm,3)) - for ia in range(mol.natm): - shl0, shl1, p0, p1 = offsetdic[ia] - - mol.set_rinv_origin(mol.atom_coord(ia)) - h1ao = -mol.atom_charge(ia) * mol.intor('int1e_iprinv', comp=3) - h1ao[:,p0:p1] += h1[:,p0:p1] - h1ao = h1ao + h1ao.transpose(0,2,1) - h1ao += pcm_nuc[ia] - h1mo = numpy.einsum('pi,xpq,qj->xij', mo_coeff, h1ao, mo_coeff) - s1mo = numpy.einsum('pi,xpq,qj->xij', mo_coeff[p0:p1], s1[:,p0:p1], mo_coeff) - s1mo = s1mo + s1mo.transpose(0,2,1) - - f1 = h1mo - numpy.einsum('xpq,pq->xpq', s1mo, zeta) - f1-= numpy.einsum('klpq,xlk->xpq', g[:nocc,:nocc], s1mo[:,:nocc,:nocc]) - - eri1a = eri1.copy() - eri1a[:,:p0] = 0 - eri1a[:,p1:] = 0 - eri1a = eri1a + eri1a.transpose(0,2,1,3,4) - eri1a = eri1a + eri1a.transpose(0,3,4,1,2) - g1 = lib.einsum('xpqrs,pi,qj,rk,sl->xijkl', eri1a, mo_coeff, mo_coeff, mo_coeff, mo_coeff) - tmp1 = lib.einsum('xpqrs,pi,qj,rk,sl->xijkl', B1[ia], mo_coeff, mo_coeff, mo_coeff, mo_coeff) - g1 = g1 * 2 - g1.transpose(0,1,4,3,2) - g1 += tmp1 * 2 - f1 += numpy.einsum('xkkpq->xpq', g1[:,:nocc,:nocc]) - f1ai = f1[:,nocc:,:nocc].copy() - - c1 = s1mo * -.5 - c1vo = cphf.solve(fvind, mo_energy, mo_occ, f1ai, max_cycle=50)[0] - c1[:,nocc:,:nocc] = c1vo - c1[:,:nocc,nocc:] = -(s1mo[:,nocc:,:nocc]+c1vo).transpose(0,2,1) - f1 += numpy.einsum('kapq,xak->xpq', g[:nocc,nocc:], c1vo) - f1 += numpy.einsum('akpq,xak->xpq', g[nocc:,:nocc], c1vo) - - e1 = numpy.einsum('xaijb,ai,bj->x', g1[:,nocc:,:nocc,:nocc,nocc:], z, z) - e1 += numpy.einsum('xab,ai,bi->x', f1[:,nocc:,nocc:], z, z) - e1 -= numpy.einsum('xij,ai,aj->x', f1[:,:nocc,:nocc], z, z) - - g1 = numpy.einsum('pjkl,xpi->xijkl', g, c1) - g1 += numpy.einsum('ipkl,xpj->xijkl', g, c1) - g1 += numpy.einsum('ijpl,xpk->xijkl', g, c1) - g1 += numpy.einsum('ijkp,xpl->xijkl', g, c1) - e1 += numpy.einsum('xaijb,ai,bj->x', g1[:,nocc:,:nocc,:nocc,nocc:], z, z) - - de[ia] = e1 - - return de - -def nuc_part(pcmobj, r_vdw, ui, ylm_1sph, cached_pol, L): - '''0th order''' - mol = pcmobj.mol - natm = mol.natm - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - nao = mol.nao - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - grids = pcmobj.grids - - extern_point_idx = ui > 0 - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - - v_phi = numpy.zeros((natm, ngrid_1sph)) - for ia in range(natm): -# Note (-) sign is not applied to atom_charges, because (-) is explicitly -# included in rhs and L matrix - d_rs = atom_coords.reshape(-1,1,3) - cav_coords[ia] - v_phi[ia] = numpy.einsum('z,zp->p', atom_charges, 1./lib.norm(d_rs,axis=2)) - phi = -numpy.einsum('n,xn,jn,jn->jx', weights_1sph, ylm_1sph, ui, v_phi) - - Xvec = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi.ravel()) - Xvec = Xvec.reshape(natm,nlm) - - i1 = 0 - scaled_weights = numpy.empty((grids.weights.size)) - for ia in range(natm): - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + fak_pol[0].shape[1] - eta_nj = 0 - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - eta_nj += fac * numpy.einsum('mn,m->n', fak_pol[l], Xvec[ia,p0:p1]) - scaled_weights[i0:i1] = eta_nj - scaled_weights *= grids.weights - - ao = mol.eval_gto('GTOval', grids.coords) - vmat = -lib.einsum('g,gi,gj->ij', scaled_weights, ao, ao) - -# Contribution of nuclear charges to the total density -# The factor numpy.sqrt(4*numpy.pi) is due to the product of 4*pi * Y_0^0 - psi = numpy.zeros((natm, nlm)) - for ia in range(natm): - psi[ia,0] += numpy.sqrt(4*numpy.pi)/r_vdw[ia] * mol.atom_charge(ia) - - # -> Psi = SL the adjoint equation to LX = g - L_S = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi.ravel()) - L_S = L_S.reshape(natm,nlm) - xi_jn = numpy.einsum('n,jn,xn,jx->jn', weights_1sph, ui, ylm_1sph, L_S) - cav_coords = cav_coords[extern_point_idx] - xi_jn = xi_jn[extern_point_idx] - - max_memory = pcmobj.max_memory - lib.current_memory()[0] - blksize = int(max(max_memory*.9e6/8/nao**2, 400)) - - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, 'int3c2e') - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor='int3c2e', aosym='s1', cintopt=cintopt) - vmat += numpy.einsum('ijn,n->ij', v_nj, xi_jn) - return vmat - -def nuc_part1(pcmobj, r_vdw, ui, ylm_1sph, cached_pol, L): - '''1st order''' - mol = pcmobj.mol - natm = mol.natm - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - nao = mol.nao - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - grids = pcmobj.grids - aoslices = mol.aoslice_by_atom() - - extern_point_idx = ui > 0 - fi0 = ddcosmo.make_fi(pcmobj, r_vdw) - fi1 = ddcosmo_grad.make_fi1(pcmobj, pcmobj.get_atomic_radii()) - fi1[:,:,ui==0] = 0 - ui1 = -fi1 - - vmat1 = numpy.zeros((natm,3,nao,nao)) - - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - - v_phi = numpy.zeros((natm, ngrid_1sph)) - for ia in range(natm): -# Note (-) sign is not applied to atom_charges, because (-) is explicitly -# included in rhs and L matrix - d_rs = atom_coords.reshape(-1,1,3) - cav_coords[ia] - v_phi[ia] = numpy.einsum('z,zp->p', atom_charges, 1./lib.norm(d_rs,axis=2)) - phi0 = -numpy.einsum('n,xn,jn,jn->jx', weights_1sph, ylm_1sph, ui, v_phi) - - Xvec0 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi0.ravel()) - Xvec0 = Xvec0.reshape(natm,nlm) - - ngrid_1sph = weights_1sph.size - v_phi0 = numpy.empty((natm,ngrid_1sph)) - for ia in range(natm): - cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph - d_rs = atom_coords.reshape(-1,1,3) - cav_coords - v_phi0[ia] = numpy.einsum('z,zp->p', atom_charges, 1./lib.norm(d_rs,axis=2)) - phi1 = -numpy.einsum('n,ln,azjn,jn->azjl', weights_1sph, ylm_1sph, ui1, v_phi0) - - for ia in range(natm): - cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph - for ja in range(natm): - rs = atom_coords[ja] - cav_coords - d_rs = lib.norm(rs, axis=1) - v_phi = atom_charges[ja] * numpy.einsum('px,p->px', rs, 1./d_rs**3) - tmp = numpy.einsum('n,ln,n,nx->xl', weights_1sph, ylm_1sph, ui[ia], v_phi) - phi1[ja,:,ia] += tmp # response of the other atoms - phi1[ia,:,ia] -= tmp # response of cavity grids - - L1 = ddcosmo_grad.make_L1(pcmobj, r_vdw, ylm_1sph, fi0) - - phi1 -= lib.einsum('aziljm,jm->azil', L1, Xvec0) - Xvec1 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi1.reshape(-1,natm*nlm).T) - Xvec1 = Xvec1.T.reshape(natm,3,natm,nlm) - - i1 = 0 - for ia, (coords, weight, weight1) in enumerate(rks_grad.grids_response_cc(grids)): - i0, i1 = i1, i1 + weight.size - ao = mol.eval_gto('GTOval_sph_deriv1', coords) - aow = numpy.einsum('gi,g->gi', ao[0], weight) - aopair1 = lib.einsum('xgi,gj->xgij', ao[1:], aow) - aow = numpy.einsum('gi,zxg->zxgi', ao[0], weight1) - aopair0 = lib.einsum('zxgi,gj->zxgij', aow, ao[0]) - - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - fac_pol = ddcosmo._vstack_factor_fak_pol(fak_pol, lmax) - vmat1 -= numpy.einsum('m,mn,zxnij->zxij', Xvec0[ia], fac_pol, aopair0) - vtmp = numpy.einsum('m,mn,xnij->xij', Xvec0[ia],fac_pol, aopair1) - vmat1[ia,:] -= vtmp - vmat1[ia,:] -= vtmp.transpose(0,2,1) - - for ja in range(natm): - shl0, shl1, p0, p1 = aoslices[ja] - vmat1[ja,:,p0:p1,:] += vtmp[:,p0:p1] - vmat1[ja,:,:,p0:p1] += vtmp[:,p0:p1].transpose(0,2,1) - - scaled_weights = lib.einsum('azm,mn->azn', Xvec1[:,:,ia], fac_pol) - scaled_weights *= weight - aow = numpy.einsum('gi,azg->azgi', ao[0], scaled_weights) - vmat1 -= numpy.einsum('gi,azgj->azij', ao[0], aow) - - psi0 = numpy.zeros((natm, nlm)) - for ia in range(natm): - psi0[ia,0] += numpy.sqrt(4*numpy.pi)/r_vdw[ia] * mol.atom_charge(ia) - - LS0 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi0.ravel()) - LS0 = LS0.reshape(natm,nlm) - - LS1 = numpy.einsum('il,aziljm->azjm', LS0, L1) - LS1 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, LS1.reshape(-1,natm*nlm).T) - LS1 = LS1.T.reshape(natm,3,natm,nlm) - - int3c2e = mol._add_suffix('int3c2e') - int3c2e_ip1 = mol._add_suffix('int3c2e_ip1') - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e) - for ia in range(natm): - cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph - #fakemol = gto.fakemol_for_charges(cav_coords[ui[ia]>0]) - fakemol = gto.fakemol_for_charges(cav_coords) - wtmp = lib.einsum('l,n,ln->ln', LS0[ia], weights_1sph, ylm_1sph) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt) - vmat1 -= numpy.einsum('azl,n,ln,n,pqn->azpq', LS1[:,:,ia], weights_1sph, ylm_1sph, ui[ia], v_nj) - vmat1 += lib.einsum('ln,azn,ijn->azij', wtmp, ui1[:,:,ia], v_nj) - - v_e1_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1, comp=3, aosym='s1') - vtmp = lib.einsum('ln,n,xijn->xij', wtmp, ui[ia], v_e1_nj) - vmat1[ia] += vtmp - vmat1[ia] += vtmp.transpose(0,2,1) - - for ja in range(natm): - shl0, shl1, p0, p1 = aoslices[ja] - vmat1[ja,:,p0:p1,:] -= vtmp[:,p0:p1] - vmat1[ja,:,:,p0:p1] -= vtmp[:,p0:p1].transpose(0,2,1) - - return vmat1 - -def make_B(pcmobj, r_vdw, ui, ylm_1sph, cached_pol, L): - '''0th order''' - mol = pcmobj.mol - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - mol = pcmobj.mol - natm = mol.natm - nao = mol.nao - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - grids = pcmobj.grids - - extern_point_idx = ui > 0 - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - - cav_coords = cav_coords[extern_point_idx] - int3c2e = mol._add_suffix('int3c2e') - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e) - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt) - nao_pair = v_nj.shape[0] - v_phi = numpy.zeros((natm, ngrid_1sph, nao, nao)) - v_phi[extern_point_idx] += v_nj.transpose(2,0,1) - - phi = numpy.einsum('n,xn,jn,jnpq->jxpq', weights_1sph, ylm_1sph, ui, v_phi) - - Xvec = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi.reshape(natm*nlm,-1)) - Xvec = Xvec.reshape(natm,nlm,nao,nao) - - ao = mol.eval_gto('GTOval', grids.coords) - aow = numpy.einsum('gi,g->gi', ao, grids.weights) - aopair = numpy.einsum('gi,gj->gij', ao, aow) - - psi = numpy.zeros((natm, nlm, nao, nao)) - i1 = 0 - for ia in range(natm): - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + fak_pol[0].shape[1] - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - psi[ia,p0:p1] = -fac * numpy.einsum('mn,nij->mij', fak_pol[l], aopair[i0:i1]) - - B = lib.einsum('nlpq,nlrs->pqrs', psi, Xvec) - B = B + B.transpose(2,3,0,1) - return B - -def make_B1(pcmobj, r_vdw, ui, ylm_1sph, cached_pol, L): - '''1st order''' - mol = pcmobj.mol - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - mol = pcmobj.mol - natm = mol.natm - nao = mol.nao - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - grids = pcmobj.grids - - extern_point_idx = ui > 0 - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - - cav_coords = cav_coords[extern_point_idx] - int3c2e = mol._add_suffix('int3c2e') - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, - mol._env, int3c2e) - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt) - nao_pair = v_nj.shape[0] - v_phi = numpy.zeros((natm, ngrid_1sph, nao, nao)) - v_phi[extern_point_idx] += v_nj.transpose(2,0,1) - phi0 = numpy.einsum('n,xn,jn,jnpq->jxpq', weights_1sph, ylm_1sph, ui, v_phi) - - Xvec0 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi0.reshape(natm*nlm,-1)) - Xvec0 = Xvec0.reshape(natm,nlm,nao,nao) - - ao = mol.eval_gto('GTOval', grids.coords) - aow = numpy.einsum('gi,g->gi', ao, grids.weights) - aopair = numpy.einsum('gi,gj->gij', ao, aow) - - psi0 = numpy.zeros((natm, nlm, nao, nao)) - i1 = 0 - for ia in range(natm): - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - i0, i1 = i1, i1 + fak_pol[0].shape[1] - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - psi0[ia,p0:p1] = -fac * numpy.einsum('mn,nij->mij', fak_pol[l], aopair[i0:i1]) - - fi0 = ddcosmo.make_fi(pcmobj, r_vdw) - fi1 = ddcosmo_grad.make_fi1(pcmobj, pcmobj.get_atomic_radii()) - fi1[:,:,ui==0] = 0 - ui1 = -fi1 - - phi1 = numpy.zeros(ui1.shape[:3] + (nlm,nao,nao)) - int3c2e = mol._add_suffix('int3c2e') - int3c2e_ip1 = mol._add_suffix('int3c2e_ip1') - aoslices = mol.aoslice_by_atom() - for ia in range(natm): - cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph - #fakemol = gto.fakemol_for_charges(cav_coords[ui[ia]>0]) - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1') - phi1[:,:,ia] += lib.einsum('n,ln,azn,ijn->azlij', weights_1sph, ylm_1sph, ui1[:,:,ia], v_nj) - - v_e1_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1, comp=3, aosym='s1') - v_e2_nj = v_e1_nj + v_e1_nj.transpose(0,2,1,3) - phi1[ia,:,ia] += lib.einsum('n,ln,n,xijn->xlij', weights_1sph, ylm_1sph, ui[ia], v_e2_nj) - - for ja in range(natm): - shl0, shl1, p0, p1 = aoslices[ja] - v = numpy.einsum('n,ln,n,xijn->xlij', weights_1sph, ylm_1sph, ui[ia], v_e1_nj[:,p0:p1]) - phi1[ja,:,ia,:,p0:p1,:] -= v - phi1[ja,:,ia,:,:,p0:p1] -= v.transpose(0,1,3,2) - - psi1 = numpy.zeros((natm,3,natm,nlm,nao,nao)) - i1 = 0 - for ia, (coords, weight, weight1) in enumerate(rks_grad.grids_response_cc(grids)): - i0, i1 = i1, i1 + weight.size - ao = mol.eval_gto('GTOval_sph_deriv1', coords) - aow = numpy.einsum('gi,g->gi', ao[0], weight) - aopair1 = lib.einsum('xgi,gj->xgij', ao[1:], aow) - aow = numpy.einsum('gi,zxg->zxgi', ao[0], weight1) - aopair0 = lib.einsum('zxgi,gj->zxgij', aow, ao[0]) - - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - p1 = 0 - for l in range(lmax+1): - fac = 4*numpy.pi/(l*2+1) - p0, p1 = p1, p1 + (l*2+1) - psi1[: ,:,ia,p0:p1] -= fac*numpy.einsum('mn,zxnij->zxmij', fak_pol[l], aopair0) - vtmp = fac*numpy.einsum('mn,xnij->xmij', fak_pol[l], aopair1) - psi1[ia,:,ia,p0:p1] -= vtmp - psi1[ia,:,ia,p0:p1] -= vtmp.transpose(0,1,3,2) - - for ja in range(natm): - shl0, shl1, q0, q1 = aoslices[ja] - psi1[ja,:,ia,p0:p1,q0:q1,:] += vtmp[:,:,q0:q1] - psi1[ja,:,ia,p0:p1,:,q0:q1] += vtmp[:,:,q0:q1].transpose(0,1,3,2) - - L1 = ddcosmo_grad.make_L1(pcmobj, r_vdw, ylm_1sph, fi0) - - Xvec1 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi1.transpose(2,3,0,1,4,5).reshape(natm*nlm,-1)) - Xvec1 = Xvec1.reshape(natm,nlm,natm,3,nao,nao).transpose(2,3,0,1,4,5) - LS0 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi0.reshape(natm*nlm,-1)) - LS0 = LS0.reshape(natm,nlm,nao,nao) - - B = lib.einsum('ixnlpq,nlrs->ixpqrs', psi1, Xvec0) - B+= lib.einsum('nlpq,ixnlrs->ixpqrs', psi0, Xvec1) - B-= lib.einsum('ilpq,aziljm,jmrs->azpqrs', LS0, L1, Xvec0) - B = B + B.transpose(0,1,4,5,2,3) - return B - -def B1_dot_x(pcmobj, dm, r_vdw, ui, ylm_1sph, cached_pol, L): - mol = pcmobj.mol - mol = pcmobj.mol - natm = mol.natm - nao = mol.nao - lmax = pcmobj.lmax - nlm = (lmax+1)**2 - - dms = numpy.asarray(dm) - is_single_dm = dms.ndim == 2 - dms = dms.reshape(-1,nao,nao) - n_dm = dms.shape[0] - - atom_coords = mol.atom_coords() - atom_charges = mol.atom_charges() - aoslices = mol.aoslice_by_atom() - grids = pcmobj.grids - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ngrid_1sph = coords_1sph.shape[0] - - extern_point_idx = ui > 0 - fi0 = ddcosmo.make_fi(pcmobj, r_vdw) - fi1 = ddcosmo_grad.make_fi1(pcmobj, pcmobj.get_atomic_radii()) - fi1[:,:,ui==0] = 0 - ui1 = -fi1 - - Bx = numpy.zeros((natm,3,nao,nao)) - - ao = mol.eval_gto('GTOval', grids.coords) - aow = numpy.einsum('gi,g->gi', ao, grids.weights) - aopair = numpy.einsum('gi,gj->gij', ao, aow) - den = numpy.einsum('gij,ij->g', aopair, dm) - psi0 = numpy.zeros((natm, nlm)) - i1 = 0 - for ia in range(natm): - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - fac_pol = ddcosmo._vstack_factor_fak_pol(fak_pol, lmax) - i0, i1 = i1, i1 + fac_pol.shape[1] - psi0[ia] = -numpy.einsum('mn,n->m', fac_pol, den[i0:i1]) - LS0 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi0.ravel()) - LS0 = LS0.reshape(natm,nlm) - - phi0 = numpy.zeros((natm,nlm)) - phi1 = numpy.zeros((natm,3,natm,nlm)) - int3c2e = mol._add_suffix('int3c2e') - int3c2e_ip1 = mol._add_suffix('int3c2e_ip1') - cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e) - cintopt_ip1 = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1) - for ia in range(natm): - cav_coords = atom_coords[ia] + r_vdw[ia] * coords_1sph - #fakemol = gto.fakemol_for_charges(cav_coords[ui[ia]>0]) - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt) - - v_phi = numpy.einsum('pqg,pq->g', v_nj, dm) - phi0[ia] = numpy.einsum('n,ln,n,n->l', weights_1sph, ylm_1sph, ui[ia], v_phi) - phi1[:,:,ia] += lib.einsum('n,ln,azn,n->azl', weights_1sph, ylm_1sph, ui1[:,:,ia], v_phi) - Bx += lib.einsum('l,n,ln,azn,ijn->azij', LS0[ia], weights_1sph, ylm_1sph, ui1[:,:,ia], v_nj) - - wtmp = lib.einsum('n,ln,n->ln', weights_1sph, ylm_1sph, ui[ia]) - - v_e1_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1, comp=3, - aosym='s1', cintopt=cintopt_ip1) - vtmp = lib.einsum('l,ln,xijn->xij', LS0[ia], wtmp, v_e1_nj) - Bx[ia] += vtmp - Bx[ia] += vtmp.transpose(0,2,1) - - for ja in range(natm): - shl0, shl1, p0, p1 = aoslices[ja] - Bx[ja,:,p0:p1,:] -= vtmp[:,p0:p1] - Bx[ja,:,:,p0:p1] -= vtmp[:,p0:p1].transpose(0,2,1) - tmp = numpy.einsum('xijn,ij->xn', v_e1_nj[:,p0:p1], dm[p0:p1]) - tmp += numpy.einsum('xijn,ji->xn', v_e1_nj[:,p0:p1], dm[:,p0:p1]) - phitmp = numpy.einsum('ln,xn->xl', wtmp, tmp) - phi1[ja,:,ia] -= phitmp - phi1[ia,:,ia] += phitmp - - Xvec0 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi0.ravel()) - Xvec0 = Xvec0.reshape(natm,nlm) - - L1 = ddcosmo_grad.make_L1(pcmobj, r_vdw, ylm_1sph, fi0) - - phi1 -= lib.einsum('aziljm,jm->azil', L1, Xvec0) - Xvec1 = numpy.linalg.solve(L.reshape(natm*nlm,-1), phi1.reshape(-1,natm*nlm).T) - Xvec1 = Xvec1.T.reshape(natm,3,natm,nlm) - - psi1 = numpy.zeros((natm,3,natm,nlm)) - i1 = 0 - for ia, (coords, weight, weight1) in enumerate(rks_grad.grids_response_cc(grids)): - i0, i1 = i1, i1 + weight.size - ao = mol.eval_gto('GTOval_sph_deriv1', coords) - aow = numpy.einsum('gi,g->gi', ao[0], weight) - aopair1 = lib.einsum('xgi,gj->xgij', ao[1:], aow) - aow = numpy.einsum('gi,zxg->zxgi', ao[0], weight1) - aopair0 = lib.einsum('zxgi,gj->zxgij', aow, ao[0]) - den0 = numpy.einsum('zxgij,ij->zxg', aopair0, dm) - den1 = numpy.empty((natm,3,weight.size)) - for ja in range(natm): - shl0, shl1, p0, p1 = aoslices[ja] - den1[ja] = numpy.einsum('xgij,ij->xg', aopair1[:,:,p0:p1], dm[p0:p1,:]) - den1[ja]+= numpy.einsum('xgij,ji->xg', aopair1[:,:,p0:p1], dm[:,p0:p1]) - - fak_pol, leak_idx = cached_pol[mol.atom_symbol(ia)] - fac_pol = ddcosmo._vstack_factor_fak_pol(fak_pol, lmax) - scaled_weights = lib.einsum('azm,mn->azn', Xvec1[:,:,ia], fac_pol) - scaled_weights *= weight - aow = numpy.einsum('gi,azg->azgi', ao[0], scaled_weights) - Bx -= numpy.einsum('gi,azgj->azij', ao[0], aow) - - tmp = numpy.einsum('mn,zxn->zxm', fac_pol, den1) - psi1[: ,:,ia] -= numpy.einsum('mn,zxn->zxm', fac_pol, den0) - psi1[ia,:,ia] -= tmp.sum(axis=0) - for ja in range(natm): - psi1[ja,:,ia] += tmp[ja] - - eta_nj = lib.einsum('mn,m->n', fac_pol, Xvec0[ia]) - Bx -= lib.einsum('n,zxnpq->zxpq', eta_nj, aopair0) - vtmp = lib.einsum('n,xnpq->xpq', eta_nj, aopair1) - Bx[ia] -= vtmp - Bx[ia] -= vtmp.transpose(0,2,1) - for ja in range(natm): - shl0, shl1, q0, q1 = aoslices[ja] - Bx[ja,:,q0:q1,:] += vtmp[:,q0:q1] - Bx[ja,:,:,q0:q1] += vtmp[:,q0:q1].transpose(0,2,1) - - psi1 -= numpy.einsum('il,aziljm->azjm', LS0, L1) - LS1 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi1.reshape(-1,natm*nlm).T) - LS1 = LS1.T.reshape(natm,3,natm,nlm) - - cav_coords = (atom_coords.reshape(natm,1,3) - + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph)) - cav_coords = cav_coords[extern_point_idx] - fakemol = gto.fakemol_for_charges(cav_coords) - v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt) - v_phi = numpy.zeros((natm, ngrid_1sph, nao, nao)) - v_phi[extern_point_idx] += v_nj.transpose(2,0,1) - Bx += lib.einsum('azjx,n,xn,jn,jnpq->azpq', LS1, weights_1sph, ylm_1sph, ui, v_phi) - - return Bx - - -def setUpModule(): - global dx, mol0, mol1, mol2, nao, dm - dx = 0.0001 - mol0 = gto.M(atom='H 0 0 0; H 0 1 1.2; H 1. .1 0; H .5 .5 1', unit='B') - mol1 = gto.M(atom='H 0 0 %g; H 0 1 1.2; H 1. .1 0; H .5 .5 1'%(-dx), unit='B') - mol2 = gto.M(atom='H 0 0 %g; H 0 1 1.2; H 1. .1 0; H .5 .5 1'%dx, unit='B') - dx = dx * 2 - nao = mol0.nao_nr() - numpy.random.seed(1) - dm = numpy.random.random((nao,nao)) - dm = dm + dm.T - -def tearDownModule(): - global dx, mol0, mol1, mol2, nao, dm - del dx, mol0, mol1, mol2, nao, dm - -class KnownValues(unittest.TestCase): - - def test_e_psi1(self): - def get_e_psi1(pcmobj): - pcmobj.grids.build() - mol = pcmobj.mol - natm = mol.natm - lmax = pcmobj.lmax - - r_vdw = ddcosmo.get_atomic_radii(pcmobj) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, lmax, True)) - - fi = ddcosmo.make_fi(pcmobj, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - nexposed = numpy.count_nonzero(ui==1) - nbury = numpy.count_nonzero(ui==0) - on_shell = numpy.count_nonzero(ui>0) - nexposed - - nlm = (lmax+1)**2 - Lmat = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi) - Lmat = Lmat.reshape(natm*nlm,-1) - - cached_pol = ddcosmo.cache_fake_multipoles(pcmobj.grids, r_vdw, lmax) - - phi = ddcosmo.make_phi(pcmobj, dm, r_vdw, ui, ylm_1sph) - L_X = numpy.linalg.solve(Lmat, phi.ravel()).reshape(natm,-1) - psi, vmat, L_S = \ - ddcosmo.make_psi_vmat(pcmobj, dm, r_vdw, ui, ylm_1sph, - cached_pol, L_X, Lmat) - psi1 = ddcosmo_grad.make_e_psi1(pcmobj, dm, r_vdw, ui, ylm_1sph, - cached_pol, L_X, Lmat) - return L_X, psi, psi1 - - pcmobj = ddcosmo.DDCOSMO(mol0) - L_X, psi0, psi1 = get_e_psi1(pcmobj) - - pcmobj = ddcosmo.DDCOSMO(mol1) - L_X1, psi = get_e_psi1(pcmobj)[:2] - e1 = numpy.einsum('jx,jx', psi, L_X) - - pcmobj = ddcosmo.DDCOSMO(mol2) - L_X2, psi = get_e_psi1(pcmobj)[:2] - e2 = numpy.einsum('jx,jx', psi, L_X) - self.assertAlmostEqual(abs((e2-e1)/dx - psi1[0,2]).max(), 0, 7) - - def test_phi(self): - def get_phi1(pcmojb): - pcmobj.grids.build() - mol = pcmobj.mol - natm = mol.natm - lmax = pcmobj.lmax - - r_vdw = pcmobj.get_atomic_radii() - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, lmax, True)) - - fi = ddcosmo.make_fi(pcmobj, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - nexposed = numpy.count_nonzero(ui==1) - nbury = numpy.count_nonzero(ui==0) - on_shell = numpy.count_nonzero(ui>0) - nexposed - - nlm = (lmax+1)**2 - Lmat = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi) - Lmat = Lmat.reshape(natm*nlm,-1) - - cached_pol = ddcosmo.cache_fake_multipoles(pcmobj.grids, r_vdw, lmax) - - phi = ddcosmo.make_phi(pcmobj, dm, r_vdw, ui, ylm_1sph) - L_X = numpy.linalg.solve(Lmat, phi.ravel()).reshape(natm,-1) - psi, vmat, L_S = \ - ddcosmo.make_psi_vmat(pcmobj, dm, r_vdw, ui, ylm_1sph, - cached_pol, L_X, Lmat) - phi1 = ddcosmo_grad.make_phi1(pcmobj, dm, r_vdw, ui, ylm_1sph) - phi1 = numpy.einsum('izjx,jx->iz', phi1, L_S) - return L_S, phi, phi1 - - pcmobj = ddcosmo.DDCOSMO(mol0) - L_S, phi0, phi1 = get_phi1(pcmobj) - - pcmobj = ddcosmo.DDCOSMO(mol1) - L_S1, phi = get_phi1(pcmobj)[:2] - e1 = numpy.einsum('jx,jx', phi, L_S) - - pcmobj = ddcosmo.DDCOSMO(mol2) - L_S2, phi = get_phi1(pcmobj)[:2] - e2 = numpy.einsum('jx,jx', phi, L_S) - self.assertAlmostEqual(abs((e2-e1)/dx - phi1[0,2]).max(), 0, 6) - - def test_fi(self): - pcmobj = ddcosmo.DDCOSMO(mol0) - fi1 = ddcosmo_grad.make_fi1(pcmobj, pcmobj.get_atomic_radii()) - ui1 = -fi1 - fi = ddcosmo.make_fi(pcmobj, pcmobj.get_atomic_radii()) - ui = 1 - fi - ui1[:,:,ui<0] = 0 - - pcmobj = ddcosmo.DDCOSMO(mol1) - fi_1 = ddcosmo.make_fi(pcmobj, pcmobj.get_atomic_radii()) - ui_1 = 1 - fi_1 - ui_1[ui_1<0] = 0 - - pcmobj = ddcosmo.DDCOSMO(mol2) - fi_2 = ddcosmo.make_fi(pcmobj, pcmobj.get_atomic_radii()) - ui_2 = 1 - fi_2 - ui_2[ui_2<0] = 0 - self.assertAlmostEqual(abs((fi_2-fi_1)/dx - fi1[0,2]).max(), 0, 5) - self.assertAlmostEqual(abs((ui_2-ui_1)/dx - ui1[0,2]).max(), 0, 5) - - def test_L1(self): - pcmobj = ddcosmo.DDCOSMO(mol0) - r_vdw = pcmobj.get_atomic_radii() - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcmobj.lmax, True)) - - fi = ddcosmo.make_fi(pcmobj, r_vdw) - L1 = ddcosmo_grad.make_L1(pcmobj, r_vdw, ylm_1sph, fi) - - pcmobj = ddcosmo.DDCOSMO(mol1) - fi = ddcosmo.make_fi(pcmobj, r_vdw) - L_1 = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi) - - pcmobj = ddcosmo.DDCOSMO(mol2) - fi = ddcosmo.make_fi(pcmobj, r_vdw) - L_2 = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi) - self.assertAlmostEqual(abs((L_2-L_1)/dx - L1[0,2]).max(), 0, 6) - - def test_e_cosmo_grad(self): - pcmobj = ddcosmo.DDCOSMO(mol0) - de = ddcosmo_grad.kernel(pcmobj, dm) - pcmobj = ddcosmo.DDCOSMO(mol1) - e1 = pcmobj.energy(dm) - pcmobj = ddcosmo.DDCOSMO(mol2) - e2 = pcmobj.energy(dm) - self.assertAlmostEqual(abs((e2-e1)/dx - de[0,2]).max(), 0, 6) - - def test_scf_grad(self): - mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol0)).run() - # solvent only - de_cosmo = ddcosmo_grad.kernel(mf.with_solvent, mf.make_rdm1()) - self.assertAlmostEqual(lib.fp(de_cosmo), 0.000902640319, 6) - # solvent + solute - de = mf.nuc_grad_method().kernel() - self.assertAlmostEqual(lib.fp(de), -0.191856565, 6) - - dm1 = mf.make_rdm1() - - mf1 = ddcosmo.ddcosmo_for_scf(scf.RHF(mol1)).run() - e1 = mf1.e_tot - e1_cosmo = mf1.with_solvent.energy(dm1) - - mf2 = ddcosmo.ddcosmo_for_scf(scf.RHF(mol2)).run() - e2 = mf2.e_tot - e2_cosmo = mf2.with_solvent.energy(dm1) - self.assertAlmostEqual(abs((e2-e1)/dx - de[0,2]).max(), 0, 7) - self.assertAlmostEqual(abs((e2_cosmo-e1_cosmo)/dx - de_cosmo[0,2]).max(), 0, 7) - - sc = mf.nuc_grad_method().as_scanner() - e, g = sc('H 0 1 0; H 0 1 1.2; H 1. 0 0; H .5 .5 0') - self.assertAlmostEqual(e, -0.83152362, 8) - self.assertAlmostEqual(lib.fp(g), 0.068317954, 6) - - mol3 = gto.M(atom='H 0 1 0; H 0 1 1.2; H 1. 0 0; H .5 .5 0', unit='B') - mf = ddcosmo.ddcosmo_for_scf(scf.RHF(mol3)).run() - de = mf.nuc_grad_method().kernel() - self.assertAlmostEqual(lib.fp(de), 0.0683179013, 6) - - def test_casci_grad(self): - mf = scf.RHF(mol0).ddCOSMO().run() - mc = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2)) - e, de = mc.nuc_grad_method().as_scanner()(mol0) - self.assertAlmostEqual(e, -1.18433554, 7) - self.assertAlmostEqual(lib.fp(de), -0.18543118, 5) - - mf = scf.RHF(mol1).run() - mc1 = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2)).run() - e1 = mc1.e_tot - mf = scf.RHF(mol2).run() - mc2 = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2)).run() - e2 = mc2.e_tot - self.assertAlmostEqual((e2-e1)/dx, de[0,2], 3) - - ## FIXME: seems working? - ## frozen dm in CASCI - #mf = scf.RHF(mol0).ddCOSMO().run() - #mc = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2), dm=mf.make_rdm1()) - #e, de = mc.nuc_grad_method().as_scanner()(mol0) - #self.assertAlmostEqual(e, -1.1845042661517311, 7) - #self.assertAlmostEqual(lib.fp(de), -0.18563349186388467, 5) - - #mf = scf.RHF(mol1).run() - #mc1 = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2), dm=mf.make_rdm1()).run() - #e1 = mc1.e_tot - #mf = scf.RHF(mol2).run() - #mc2 = solvent.ddCOSMO(mcscf.CASCI(mf, 2, 2), dm=mf.make_rdm1()).run() - #e2 = mc2.e_tot - #self.assertAlmostEqual((e2-e1)/dx, de[0,2], 4) - - def test_casscf_grad(self): - mf = scf.RHF(mol0).ddCOSMO().run() - mc = solvent.ddCOSMO(mcscf.CASSCF(mf, 2, 2)).set(conv_tol=1e-9) - mc_g = mc.nuc_grad_method().as_scanner() - e, de = mc_g(mol0) - self.assertAlmostEqual(e, -1.19627418, 5) - self.assertAlmostEqual(lib.fp(de), -0.1831184, 4) - - mf = scf.RHF(mol1).run() - mc1 = solvent.ddCOSMO(mcscf.CASSCF(mf, 2, 2)).run(conv_tol=1e-9) - e1 = mc1.e_tot - mf = scf.RHF(mol2).run() - mc2 = solvent.ddCOSMO(mcscf.CASSCF(mf, 2, 2)).run(conv_tol=1e-9) - e2 = mc2.e_tot - # ddcosmo-CASSCF is not fully variational. Errors will be found large - # in this test. - self.assertAlmostEqual((e2-e1)/dx, de[0,2], 2) - - def test_ccsd_grad(self): - mf = scf.RHF(mol0).ddCOSMO().run() - mycc = cc.CCSD(mf).ddCOSMO() - e, de = mycc.nuc_grad_method().as_scanner()(mol0) - self.assertAlmostEqual(e, -1.2060391657, 7) - self.assertAlmostEqual(lib.fp(de), -0.1794318433, 5) - - mf = scf.RHF(mol1).run() - mycc1 = solvent.ddCOSMO(cc.CCSD(mf)).run() - e1 = mycc1.e_tot - mf = scf.RHF(mol2).run() - mycc2 = solvent.ddCOSMO(cc.CCSD(mf)).run() - e2 = mycc2.e_tot - self.assertAlmostEqual((e2-e1)/dx, de[0,2], 4) - - def test_tda_grad(self): - mol0 = gto.M(atom='H 0 0 0 ; H .5 .5 .1', unit='B', basis='321g') - mol1 = gto.M(atom='H 0 0 -.001; H .5 .5 .1', unit='B', basis='321g') - mol2 = gto.M(atom='H 0 0 0.001; H .5 .5 .1', unit='B', basis='321g') - mf = scf.RHF(mol0).ddCOSMO().run() - td = solvent.ddCOSMO(tdscf.TDA(mf)).run(equilibrium_solvation=True) - ref = tda_grad(td, td.xy[0]) + mf.nuc_grad_method().kernel() - - e, de = td.nuc_grad_method().as_scanner(state=1)(mol0) - de = td.nuc_grad_method().kernel() - self.assertAlmostEqual(abs(ref - de).max(), 0, 12) - - td1 = mol1.RHF().ddCOSMO().run().TDA().ddCOSMO().run(equilibrium_solvation=True) - td2 = mol2.RHF().ddCOSMO().run().TDA().ddCOSMO().run(equilibrium_solvation=True) - e1 = td1.e_tot[0] - e2 = td2.e_tot[0] - self.assertAlmostEqual((e2-e1)/0.002, de[0,2], 5) - - def test_solvent_nuc(self): - def get_nuc(mol): - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 2 - pcm.eps = 0 - natm = mol.natm - nao = mol.nao - nlm = (pcm.lmax+1)**2 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - pcm.grids = grids = ddcosmo.Grids(mol).run(level=0) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - return nuc_part(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - - pcm = ddcosmo.DDCOSMO(mol0) - pcm.lmax = 2 - pcm.eps = 0 - natm = mol0.natm - nao = mol0.nao - nlm = (pcm.lmax+1)**2 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - pcm.grids = grids = ddcosmo.Grids(mol0).run(level=0) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - dvmat = nuc_part1(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - - vmat1 = get_nuc(mol1) - vmat2 = get_nuc(mol2) - self.assertAlmostEqual(abs((vmat2-vmat1)/dx - dvmat[0,2]).max(), 0, 8) - - nao = mol0.nao - numpy.random.seed(19) - dm = numpy.random.random((nao,nao)) - vref = pcm._get_vind(dm)[1] - vmat = 0.5 * get_nuc(mol0) - vmat += pcm._B_dot_x(dm) - self.assertAlmostEqual(abs(vmat-vref).max(), 0, 14) - - dm1 = numpy.random.random((2,nao,nao)) - de = _ddcosmo_tdscf_grad._grad_ne(pcm, dm1, r_vdw, ui, ylm_1sph, cached_pol, L) - ref = numpy.einsum('azij,nij->naz', dvmat, dm1) - self.assertAlmostEqual(abs(de - ref).max(), 0, 12) - - def test_B1(self): - def getB(mol): - pcm = ddcosmo.DDCOSMO(mol) - pcm.lmax = 2 - pcm.eps = 0 - natm = mol.natm - nao = mol.nao - nlm = (pcm.lmax+1)**2 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - pcm.grids = grids = ddcosmo.Grids(mol).run(level=0) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - return make_B(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - - pcm = ddcosmo.DDCOSMO(mol0) - pcm.lmax = 2 - pcm.eps = 0 - natm = mol0.natm - nao = mol0.nao - nlm = (pcm.lmax+1)**2 - r_vdw = ddcosmo.get_atomic_radii(pcm) - fi = ddcosmo.make_fi(pcm, r_vdw) - ui = 1 - fi - ui[ui<0] = 0 - pcm.grids = grids = ddcosmo.Grids(mol0).run(level=0) - coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcm.lebedev_order) - ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcm.lmax, True)) - cached_pol = ddcosmo.cache_fake_multipoles(grids, r_vdw, pcm.lmax) - L = ddcosmo.make_L(pcm, r_vdw, ylm_1sph, fi) - dB = make_B1(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - - B1 = getB(mol1) - B2 = getB(mol2) - self.assertAlmostEqual(abs((B2-B1)/dx - dB[0,2]).max(), 0, 8) - - nao = mol0.nao - numpy.random.seed(1) - dm1 = numpy.random.random((2,nao,nao)) - dm2 = numpy.random.random((2,nao,nao)) - dm = dm1[0] - ref = numpy.einsum('azpqrs,npq->nazrs', dB, dm1) - v = B1_dot_x(pcm, dm, r_vdw, ui, ylm_1sph, cached_pol, L) - self.assertAlmostEqual(abs(v-ref[0]).max(), 0, 12) - - de = _ddcosmo_tdscf_grad._grad_ee(pcm, dm1, dm2, r_vdw, ui, ylm_1sph, cached_pol, L) - ref = numpy.einsum('nazij,nij->naz', ref, dm2) - self.assertAlmostEqual(abs(de - ref).max(), 0, 12) - - numpy.random.seed(1) - dm = numpy.random.random((nao,nao)) - dm = dm + dm.T - ref = ddcosmo_grad.kernel(pcm, dm) - dielectric = pcm.eps - if dielectric > 0: - f_epsilon = (dielectric-1.)/dielectric - else: - f_epsilon = 1 - de = _ddcosmo_tdscf_grad._grad_nn(pcm, r_vdw, ui, ylm_1sph, cached_pol, L) - de+= _ddcosmo_tdscf_grad._grad_ne(pcm, dm, r_vdw, ui, ylm_1sph, cached_pol, L) - de+= .5*_ddcosmo_tdscf_grad._grad_ee(pcm, dm, dm, r_vdw, ui, ylm_1sph, cached_pol, L) - de *= .5 * f_epsilon - self.assertAlmostEqual(abs(de-ref).max(), 0, 12) - - def test_regularize_xt(self): - pcmobj = ddcosmo.DDCOSMO(mol0) - numpy.random.seed(2) - t = numpy.random.rand(4) - eta = 0.8 - L1 = ddcosmo_grad.regularize_xt1(t, eta) - L_1 = ddcosmo.regularize_xt(t-1e-4, eta) - L_2 = ddcosmo.regularize_xt(t+1e-4, eta) - self.assertAlmostEqual(abs((L_2-L_1)/2e-4 - L1).max(), 0, 6) - - -if __name__ == "__main__": - print("Full Tests for ddcosmo gradients") - unittest.main() diff --git a/pyscf/solvent/test/test_ddcosmo_tdscf_grad.py b/pyscf/solvent/test/test_ddcosmo_tdscf_grad.py index b54fdf7ee5..a63ea02e31 100644 --- a/pyscf/solvent/test/test_ddcosmo_tdscf_grad.py +++ b/pyscf/solvent/test/test_ddcosmo_tdscf_grad.py @@ -171,15 +171,15 @@ class KnownValues(unittest.TestCase): def test_b3lyp_tda(self): # TDA gga with equilibrium_solvation - mf = mol0.RKS().ddCOSMO().run(xc='b3lyp') - td = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) - g1 = td.nuc_grad_method().kernel() - - mf = mol1.RKS().ddCOSMO().run(xc='b3lyp') - td1 = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) +# mf = mol0.RKS().ddCOSMO().run(xc='b3lyp') +# td = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) +# g1 = td.nuc_grad_method().kernel() +# +# mf = mol1.RKS().ddCOSMO().run(xc='b3lyp') +# td1 = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) mf = mol2.RKS().ddCOSMO().run(xc='b3lyp') td2 = mf.TDA().ddCOSMO().run(equilibrium_solvation=True) - self.assertAlmostEqual((td2.e_tot[0]-td1.e_tot[0])/0.002, g1[0,2], 4) +# self.assertAlmostEqual((td2.e_tot[0]-td1.e_tot[0])/0.002, g1[0,2], 4) # # TDA gga without equilibrium_solvation # mf = mol0.RKS().ddCOSMO().run(xc='b3lyp') diff --git a/pyscf/solvent/test/test_ddpcm.py b/pyscf/solvent/test/test_ddpcm.py deleted file mode 100644 index 4496bc2bff..0000000000 --- a/pyscf/solvent/test/test_ddpcm.py +++ /dev/null @@ -1,40 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2020 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 -from pyscf import gto, scf -from pyscf.solvent import ddpcm - - -class KnownValues(unittest.TestCase): - def test_ddpcm_scf(self): - mol = gto.M(atom=''' - 6 0.000000 0.000000 -0.542500 - 8 0.000000 0.000000 0.677500 - 1 0.000000 0.935307 -1.082500 - 1 0.000000 -0.935307 -1.082500 - ''', basis='sto3g', verbose=7, - output='/dev/null') - pcm = ddpcm.DDPCM(mol) - pcm.lmax = 6 - pcm.lebedev_order = 17 - mf = scf.RHF(mol).ddPCM(pcm).run() - self.assertAlmostEqual(mf.e_tot, -112.3544929827, 8) - - -if __name__ == "__main__": - print("Full Tests for ddpcm") - unittest.main() - diff --git a/pyscf/solvent/test/test_pol_embed.py b/pyscf/solvent/test/test_pol_embed.py deleted file mode 100644 index 24d6f0151c..0000000000 --- a/pyscf/solvent/test/test_pol_embed.py +++ /dev/null @@ -1,322 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2019 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 os -import tempfile -import numpy -from numpy.testing import assert_allclose -from pyscf import lib, gto, scf - -have_pe = False -try: - import cppe - import pyscf.solvent as solvent - from pyscf.solvent import pol_embed - have_pe = True -except ImportError: - pass - - -dname = os.path.dirname(__file__) - -def setUpModule(): - global potf, potf2, mol, mol2, potfile, potfile2 - potf = tempfile.NamedTemporaryFile() - potf.write(b'''! -@COORDINATES -3 -AA -O 3.53300000 2.99600000 0.88700000 1 -H 4.11100000 3.13200000 1.63800000 2 -H 4.10500000 2.64200000 0.20600000 3 -@MULTIPOLES -ORDER 0 -3 -1 -0.67444000 -2 0.33722000 -3 0.33722000 -@POLARIZABILITIES -ORDER 1 1 -3 -1 5.73935000 0.00000000 0.00000000 5.73935000 0.00000000 5.73935000 -2 2.30839000 0.00000000 0.00000000 2.30839000 0.00000000 2.30839000 -3 2.30839000 0.00000000 0.00000000 2.30839000 0.00000000 2.30839000 -EXCLISTS -3 3 -1 2 3 -2 1 3 -3 1 2''') - potf.flush() - potfile = potf.name - - mol = gto.M(atom=''' - 6 0.000000 0.000000 -0.542500 - 8 0.000000 0.000000 0.677500 - 1 0.000000 0.935307 -1.082500 - 1 0.000000 -0.935307 -1.082500 - ''', basis='sto3g', verbose=7, - output='/dev/null') - - potf2 = tempfile.NamedTemporaryFile() - potf2.write(b'''! water molecule + a large, positive charge to force electron spill-out -@COORDINATES -4 -AA -O 21.41500000 20.20300000 28.46300000 1 -H 20.84000000 20.66500000 29.07300000 2 -H 22.26700000 20.19000000 28.89900000 3 -X 21.41500000 20.20300000 28.46300000 4 -@MULTIPOLES -ORDER 0 -4 -1 -0.68195912 -2 0.34097956 -3 0.34097956 -4 8.00000000 -ORDER 1 -4 -1 0.07364864 0.11937993 0.27811003 -2 0.13201293 -0.10341957 -0.13476149 -3 -0.19289612 0.00473166 -0.09514399 -4 0 0 0 -ORDER 2 -4 -1 -3.33216702 -0.28299356 -0.01420064 -4.16411766 0.21213644 -3.93180619 -2 -0.26788067 -0.15634676 -0.21817674 -0.30867972 0.17884149 -0.20228345 -3 -0.01617821 0.00575875 0.24171375 -0.44448719 -0.00422271 -0.31817844 -4 0 0 0 0 0 0 -@POLARIZABILITIES -ORDER 1 1 -4 -1 2.98942355 -0.37779481 0.04141173 1.82815813 0.40319255 2.35026343 -2 1.31859676 -0.60024842 -0.89358696 1.20607500 0.56785158 1.40652272 -3 2.28090986 0.01949959 0.86466278 0.68674835 -0.13216306 0.96339257 -4 0 0 0 0 0 0 -EXCLISTS -4 4 -1 2 3 4 -2 1 3 4 -3 1 2 4 -4 1 2 3''') - potf2.flush() - potfile2 = potf2.name - - mol2 = gto.M( - atom=""" - O 22.931000 21.390000 23.466000 - C 22.287000 21.712000 22.485000 - N 22.832000 22.453000 21.486000 - H 21.242000 21.408000 22.312000 - H 23.729000 22.867000 21.735000 - H 22.234000 23.026000 20.883000 - """, - basis="3-21++G", - verbose=7, - output='/dev/null' - ) - - -def tearDownModule(): - global potf, potf2, mol, mol2 - mol.stdout.close() - mol2.stdout.close() - del potf, potf2, mol, mol2 - - -def _compute_multipole_potential_integrals(pe, site, order, moments): - if order > 2: - raise NotImplementedError("""Multipole potential integrals not - implemented for order > 2.""") - pe.mol.set_rinv_orig(site) - # TODO: only calculate up to requested order! - integral0 = pe.mol.intor("int1e_rinv") - integral1 = pe.mol.intor("int1e_iprinv") + pe.mol.intor("int1e_iprinv").transpose(0, 2, 1) - integral2 = pe.mol.intor("int1e_ipiprinv") + pe.mol.intor("int1e_ipiprinv").transpose(0, 2, 1) + 2.0 * pe.mol.intor("int1e_iprinvip") - - # k = 2: 0,1,2,4,5,8 = XX, XY, XZ, YY, YZ, ZZ - # add the lower triangle to the upper triangle, i.e., - # XY += YX : 1 + 3 - # XZ += ZX : 2 + 6 - # YZ += ZY : 5 + 7 - # and divide by 2 - integral2[1] += integral2[3] - integral2[2] += integral2[6] - integral2[5] += integral2[7] - integral2[1] *= 0.5 - integral2[2] *= 0.5 - integral2[5] *= 0.5 - - op = integral0 * moments[0] * cppe.prefactors(0) - if order > 0: - op += numpy.einsum('aij,a->ij', integral1, - moments[1] * cppe.prefactors(1)) - if order > 1: - op += numpy.einsum('aij,a->ij', - integral2[[0, 1, 2, 4, 5, 8], :, :], - moments[2] * cppe.prefactors(2)) - - return op - -def _compute_field_integrals(pe, site, moment): - pe.mol.set_rinv_orig(site) - integral = pe.mol.intor("int1e_iprinv") + pe.mol.intor("int1e_iprinv").transpose(0, 2, 1) - op = numpy.einsum('aij,a->ij', integral, -1.0*moment) - return op - -def _compute_field(pe, site, D): - pe.mol.set_rinv_orig(site) - integral = pe.mol.intor("int1e_iprinv") + pe.mol.intor("int1e_iprinv").transpose(0, 2, 1) - return numpy.einsum('ij,aij->a', D, integral) - -def _exec_cppe(pe, dm, elec_only=False): - V_es = numpy.zeros((pe.mol.nao, pe.mol.nao), dtype=numpy.float64) - for p in pe.potentials: - moments = [] - for m in p.multipoles: - m.remove_trace() - moments.append(m.values) - V_es += _compute_multipole_potential_integrals(pe, p.position, m.k, moments) - - pe.cppe_state.energies["Electrostatic"]["Electronic"] = ( - numpy.einsum('ij,ij->', V_es, dm) - ) - - n_sitecoords = 3 * pe.cppe_state.get_polarizable_site_number() - V_ind = numpy.zeros((pe.mol.nao, pe.mol.nao), dtype=numpy.float64) - if n_sitecoords: - # TODO: use list comprehensions - current_polsite = 0 - elec_fields = numpy.zeros(n_sitecoords, dtype=numpy.float64) - for p in pe.potentials: - if not p.is_polarizable: - continue - elec_fields_s = _compute_field(pe, p.position, dm) - elec_fields[3*current_polsite:3*current_polsite + 3] = elec_fields_s - current_polsite += 1 - pe.cppe_state.update_induced_moments(elec_fields, elec_only) - induced_moments = numpy.array(pe.cppe_state.get_induced_moments()) - current_polsite = 0 - for p in pe.potentials: - if not p.is_polarizable: - continue - site = p.position - V_ind += _compute_field_integrals(pe, site=site, moment=induced_moments[3*current_polsite:3*current_polsite + 3]) - current_polsite += 1 - e = pe.cppe_state.total_energy - if not elec_only: - vmat = V_es + V_ind - else: - vmat = V_ind - e = pe.cppe_state.energies["Polarization"]["Electronic"] - return e, vmat - -@unittest.skipIf(not have_pe, "CPPE library not found.") -class TestPolEmbed(unittest.TestCase): - def test_exec_cppe(self): - pe = solvent.PE(mol, os.path.join(dname, "pna_6w.potential")) - numpy.random.seed(2) - nao = mol.nao - dm = numpy.random.random((2,nao,nao)) - dm = dm + dm.transpose(0,2,1) - - eref, vref = _exec_cppe(pe, dm[1], elec_only=False) - e, v = pe._exec_cppe(dm, elec_only=False) - self.assertAlmostEqual(abs(vref - v[1]).max(), 0, 10) - - eref, vref = _exec_cppe(pe, dm[0], elec_only=True) - e, v = pe._exec_cppe(dm, elec_only=True) - self.assertAlmostEqual(eref, e[0], 9) - self.assertAlmostEqual(abs(vref - v[0]).max(), 0, 9) - - def test_pol_embed_scf(self): - mol = gto.Mole() - mol.atom = ''' - C 8.64800 1.07500 -1.71100 - C 9.48200 0.43000 -0.80800 - C 9.39600 0.75000 0.53800 - C 8.48200 1.71200 0.99500 - C 7.65300 2.34500 0.05500 - C 7.73200 2.03100 -1.29200 - H 10.18300 -0.30900 -1.16400 - H 10.04400 0.25200 1.24700 - H 6.94200 3.08900 0.38900 - H 7.09700 2.51500 -2.01800 - N 8.40100 2.02500 2.32500 - N 8.73400 0.74100 -3.12900 - O 7.98000 1.33100 -3.90100 - O 9.55600 -0.11000 -3.46600 - H 7.74900 2.71100 2.65200 - H 8.99100 1.57500 2.99500 - ''' - mol.basis = "STO-3G" - mol.build() - pe_options = {"potfile": os.path.join(dname, "pna_6w.potential")} - pe = pol_embed.PolEmbed(mol, pe_options) - mf = solvent.PE(scf.RHF(mol), pe) - mf.conv_tol = 1e-10 - mf.kernel() - ref_pe_energy = -0.03424830892844 - ref_scf_energy = -482.9411084900 - assert_allclose(ref_pe_energy, mf.with_solvent.e, atol=1e-6) - assert_allclose(ref_scf_energy, mf.e_tot, atol=1e-6) - - def test_pe_scf(self): - pe = solvent.PE(mol, potfile) - mf = solvent.PE(mol.RHF(), pe).run(conv_tol=1e-10) - self.assertAlmostEqual(mf.e_tot, -112.35232445743728, 9) - self.assertAlmostEqual(mf.with_solvent.e, 0.00020182314249546455, 9) - - def test_pe_scf_ecp(self): - pe = solvent.PE(mol2, {"potfile": potfile2, "ecp": True}) - mf = solvent.PE(mol2.RHF(), pe).run(conv_tol=1e-10) - self.assertAlmostEqual(mf.e_tot, -168.147494986446, 8) - - def test_as_scanner(self): - mf = mol.RHF(chkfile=tempfile.NamedTemporaryFile().name) - mf_scanner = solvent.PE(mf, potfile).as_scanner() - mf_scanner(mol) - self.assertAlmostEqual(mf_scanner.with_solvent.e, 0.00020182314249546455, 9) - # Change solute. cppe may not support this - mf_scanner('H 0. 0. 0.; H 0. 0. .9') - self.assertAlmostEqual(mf_scanner.with_solvent.e, 5.2407234004672825e-05, 9) - - def test_newton_rohf(self): - mf = solvent.PE(mol.ROHF(max_memory=0), potfile) - mf = mf.newton() - e = mf.kernel() - self.assertAlmostEqual(e, -112.35232445745123, 9) - - mf = solvent.PE(mol.ROHF(max_memory=0), potfile) - e = mf.kernel() - self.assertAlmostEqual(e, -112.35232445745123, 9) - - def test_rhf_tda(self): - # TDA with equilibrium_solvation - mf = solvent.PE(mol.RHF(), potfile).run(conv_tol=1e-10) - td = solvent.PE(mf.TDA(), potfile).run(equilibrium_solvation=True) - ref = numpy.array([0.1506426609354755, 0.338251407831332, 0.4471267328974609]) - self.assertAlmostEqual(abs(ref - td.e).max(), 0, 7) - - # TDA without equilibrium_solvation - mf = solvent.PE(mol.RHF(), potfile).run(conv_tol=1e-10) - td = solvent.PE(mf.TDA(), potfile).run() - ref = numpy.array([0.1506431269137912, 0.338254809044639, 0.4471487090255076]) - self.assertAlmostEqual(abs(ref - td.e).max(), 0, 7) - - -if __name__ == "__main__": - print("Full Tests for pol_embed") - unittest.main() diff --git a/pyscf/tdscf/rhf.py b/pyscf/tdscf/rhf.py index dadfc73b7f..e7cd0f5de0 100644 --- a/pyscf/tdscf/rhf.py +++ b/pyscf/tdscf/rhf.py @@ -864,10 +864,12 @@ def pickeig(w, v, nroots, envs): idx = numpy.where(w > self.positive_eig_threshold)[0] return w[idx], v[:,idx], idx + print('tdscf x0', x0 is None) if x0 is None: x0 = self.init_guess(self._scf, self.nstates) x0 = self.trunc_workspace(vind, x0, nstates=self.nstates, pick=pickeig)[1] + print('tdscf.initial guess', x0.dtype) self.converged, self.e, x1 = \ lib.davidson1(vind, x0, precond, tol=self.conv_tol,