Source code for pyscf.solvent.ddpcm

#!/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.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

'''
domain decomposition PCM (In testing)

See also
JCP, 144, 054101
JCP, 144, 160901
'''

import warnings
import numpy
from pyscf import lib
from pyscf.lib import logger
from pyscf import gto
from pyscf.dft import gen_grid
from pyscf.data import radii
from pyscf.solvent import ddcosmo
from pyscf.symm import sph
from pyscf.solvent import _attach_solvent

warnings.warn('Module ddPCM is under testing')


[docs] @lib.with_doc(_attach_solvent._for_scf.__doc__) def ddpcm_for_scf(mf, solvent_obj=None, dm=None): if solvent_obj is None: solvent_obj = DDPCM(mf.mol) return _attach_solvent._for_scf(mf, solvent_obj, dm)
[docs] @lib.with_doc(_attach_solvent._for_casscf.__doc__) def ddpcm_for_casscf(mc, solvent_obj=None, dm=None): if solvent_obj is None: if isinstance(getattr(mc._scf, 'with_solvent', None), DDPCM): solvent_obj = mc._scf.with_solvent else: solvent_obj = DDPCM(mc.mol) return _attach_solvent._for_casscf(mc, solvent_obj, dm)
[docs] @lib.with_doc(_attach_solvent._for_casci.__doc__) def ddpcm_for_casci(mc, solvent_obj=None, dm=None): if solvent_obj is None: if isinstance(getattr(mc._scf, 'with_solvent', None), DDPCM): solvent_obj = mc._scf.with_solvent else: solvent_obj = DDPCM(mc.mol) return _attach_solvent._for_casci(mc, solvent_obj, dm)
[docs] @lib.with_doc(_attach_solvent._for_post_scf.__doc__) def ddpcm_for_post_scf(method, solvent_obj=None, dm=None): if solvent_obj is None: if isinstance(getattr(method._scf, 'with_solvent', None), DDPCM): solvent_obj = method._scf.with_solvent else: solvent_obj = DDPCM(method.mol) return _attach_solvent._for_post_scf(method, solvent_obj, dm)
[docs] @lib.with_doc(_attach_solvent._for_tdscf.__doc__) def ddpcm_for_tdscf(method, solvent_obj=None, dm=None): scf_solvent = getattr(method._scf, 'with_solvent', None) assert scf_solvent is None or isinstance(scf_solvent, DDPCM) if solvent_obj is None: solvent_obj = DDPCM(method.mol) return _attach_solvent._for_tdscf(method, solvent_obj, dm)
# Inject ddPCM to other methods from pyscf import scf from pyscf import mcscf from pyscf import mp, ci, cc from pyscf import tdscf scf.hf.SCF.ddPCM = scf.hf.SCF.DDPCM = ddpcm_for_scf mp.mp2.MP2.ddPCM = mp.mp2.MP2.DDPCM = ddpcm_for_post_scf ci.cisd.CISD.ddPCM = ci.cisd.CISD.DDPCM = ddpcm_for_post_scf cc.ccsd.CCSDBase.ddPCM = cc.ccsd.CCSDBase.DDPCM = ddpcm_for_post_scf tdscf.rhf.TDBase.ddPCM = tdscf.rhf.TDBase.DDPCM = ddpcm_for_tdscf mcscf.casci.CASCI.ddPCM = mcscf.casci.CASCI.DDPCM = ddpcm_for_casci mcscf.mc1step.CASSCF.ddPCM = mcscf.mc1step.CASSCF.DDPCM = ddpcm_for_casscf
[docs] def gen_ddpcm_solver(pcmobj, verbose=None): mol = pcmobj.mol if pcmobj.grids.coords is None: pcmobj.grids.build(with_non0tab=True) 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 logger.debug(pcmobj, 'Num points exposed %d', nexposed) logger.debug(pcmobj, 'Num points buried %d', nbury) logger.debug(pcmobj, 'Num points on shell %d', on_shell) nlm = (lmax+1)**2 Lmat = ddcosmo.make_L(pcmobj, r_vdw, ylm_1sph, fi) Lmat = Lmat.reshape(natm*nlm,-1) Amat = make_A(pcmobj, r_vdw, ylm_1sph, ui).reshape(natm*nlm,-1) fac = 2*numpy.pi * (pcmobj.eps+1) / (pcmobj.eps-1) A_diele = Amat + fac * numpy.eye(natm*nlm) A_inf = Amat + 2*numpy.pi * numpy.eye(natm*nlm) cached_pol = ddcosmo.cache_fake_multipoles(pcmobj.grids, r_vdw, lmax) def gen_vind(dm): phi = ddcosmo.make_phi(pcmobj, dm, r_vdw, ui) phi = numpy.linalg.solve(A_diele, A_inf.dot(phi.ravel())) Xvec = numpy.linalg.solve(Lmat, phi.ravel()).reshape(natm,-1) psi, vmat = ddcosmo.make_psi_vmat(pcmobj, dm, r_vdw, ui, pcmobj.grids, ylm_1sph, cached_pol, Xvec, Lmat)[:2] dielectric = pcmobj.eps f_epsilon = (dielectric-1.)/dielectric epcm = .5 * f_epsilon * numpy.einsum('jx,jx', psi, Xvec) vpcm = .5 * f_epsilon * vmat return epcm, vpcm return gen_vind
[docs] def energy(pcmobj, dm): r''' ddPCM energy Es = 1/2 f(eps) \int rho(r) W(r) dr ''' epcm = gen_ddpcm_solver(pcmobj, pcmobj.verbose)(dm)[0] return epcm
[docs] def regularize_xt(t, eta): xt = numpy.zeros_like(t) inner = t <= 1-eta on_shell = (1-eta < t) & (t < 1) xt[inner] = 1 ti = t[on_shell] - eta*.5 # JCP, 144, 054101 xt[on_shell] = 1./eta**4 * (1-ti)**2 * (ti-1+2*eta)**2 return xt
[docs] def make_A(pcmobj, r_vdw, ylm_1sph, ui): # Part of A matrix defined in JCP, 144, 054101, Eq (43), (44) mol = pcmobj.mol natm = mol.natm lmax = pcmobj.lmax # eta = pcmobj.eta nlm = (lmax+1)**2 coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(pcmobj.lebedev_order) ngrid_1sph = weights_1sph.size atom_coords = mol.atom_coords() ylm_1sph = ylm_1sph.reshape(nlm,ngrid_1sph) Amat = numpy.zeros((natm,nlm,natm,nlm)) for ja in range(natm): # w_u = precontract w_n U_j w_u = weights_1sph * ui[ja] p1 = 0 for l in range(lmax+1): fac = 2*numpy.pi/(l*2+1) p0, p1 = p1, p1 + (l*2+1) a = numpy.einsum('xn,n,mn->xm', ylm_1sph, w_u, ylm_1sph[p0:p1]) Amat[ja,:,ja,p0:p1] += -fac * a for ka in ddcosmo.atoms_with_vdw_overlap(ja, atom_coords, r_vdw): vjk = r_vdw[ja] * coords_1sph + atom_coords[ja] - atom_coords[ka] rjk = lib.norm(vjk, axis=1) pol = sph.multipoles(vjk, lmax) p1 = 0 weights = w_u / rjk**(l*2+1) for l in range(lmax+1): fac = 4*numpy.pi*l/(l*2+1) * r_vdw[ka]**(l+1) p0, p1 = p1, p1 + (l*2+1) a = numpy.einsum('xn,n,mn->xm', ylm_1sph, weights, pol[l]) Amat[ja,:,ka,p0:p1] += -fac * a return Amat
[docs] class ddPCM(ddcosmo.DDCOSMO): def __init__(self, mol): ddcosmo.DDCOSMO.__init__(self, mol) self.method = 'ddPCM'
[docs] def dump_flags(self, verbose=None): logger.info(self, '******** %s (In testing) ********', self.__class__) logger.warn(self, 'ddPCM is an experimental feature. It is ' 'still in testing.\nFeatures and APIs may be changed ' 'in the future.') logger.info(self, 'lebedev_order = %s (%d grids per sphere)', self.lebedev_order, gen_grid.LEBEDEV_ORDER[self.lebedev_order]) logger.info(self, 'lmax = %s' , self.lmax) logger.info(self, 'eta = %s' , self.eta) logger.info(self, 'eps = %s' , self.eps) logger.info(self, 'frozen = %s' , self.frozen) logger.info(self, 'equilibrium_solvation = %s', self.equilibrium_solvation) logger.debug2(self, 'radii_table %s', self.radii_table) if self.atom_radii: logger.info(self, 'User specified atomic radii %s', str(self.atom_radii)) self.grids.dump_flags(verbose) return self
[docs] def build(self): if self.grids.coords is None: self.grids.build(with_non0tab=True) mol = self.mol natm = mol.natm lmax = self.lmax r_vdw = ddcosmo.get_atomic_radii(self) coords_1sph, weights_1sph = ddcosmo.make_grids_one_sphere(self.lebedev_order) ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, lmax, True)) fi = ddcosmo.make_fi(self, 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 logger.debug(self, 'Num points exposed %d', nexposed) logger.debug(self, 'Num points buried %d', nbury) logger.debug(self, 'Num points on shell %d', on_shell) nlm = (lmax+1)**2 Lmat = ddcosmo.make_L(self, r_vdw, ylm_1sph, fi) Lmat = Lmat.reshape(natm*nlm,-1) Amat = make_A(self, r_vdw, ylm_1sph, ui).reshape(natm*nlm,-1) fac = 2*numpy.pi * (self.eps+1) / (self.eps-1) A_diele = Amat + fac * numpy.eye(natm*nlm) A_inf = Amat + 2*numpy.pi * numpy.eye(natm*nlm) cached_pol = ddcosmo.cache_fake_multipoles(self.grids, r_vdw, lmax) self._intermediates = { 'r_vdw': r_vdw, 'ylm_1sph': ylm_1sph, 'ui': ui, 'Lmat': Lmat, 'A_diele': A_diele, 'A_inf': A_inf, 'cached_pol': cached_pol, }
def _get_vind(self, dm): '''A single shot solvent effects for given density matrix. ''' if not self._intermediates or self.grids.coords is None: self.build() mol = self.mol r_vdw = self._intermediates['r_vdw' ] ylm_1sph = self._intermediates['ylm_1sph' ] ui = self._intermediates['ui' ] Lmat = self._intermediates['Lmat' ] A_diele = self._intermediates['A_diele' ] A_inf = self._intermediates['A_inf' ] cached_pol = self._intermediates['cached_pol'] if not (isinstance(dm, numpy.ndarray) and dm.ndim == 2): # spin-traced DM for UHF or ROHF dm = dm[0] + dm[1] phi = ddcosmo.make_phi(self, dm, r_vdw, ui, ylm_1sph) phi = numpy.linalg.solve(A_diele, A_inf.dot(phi.ravel())) Xvec = numpy.linalg.solve(Lmat, phi.ravel()).reshape(mol.natm,-1) psi, vmat = ddcosmo.make_psi_vmat(self, dm, r_vdw, ui, ylm_1sph, cached_pol, Xvec, Lmat)[:2] dielectric = self.eps f_epsilon = (dielectric-1.)/dielectric epcm = .5 * f_epsilon * numpy.einsum('jx,jx', psi, Xvec) vpcm = .5 * f_epsilon * vmat return epcm, vpcm def _B_dot_x(self, dm): ''' Compute the matrix-vector product B * x. The B matrix, as defined in the paper R. Cammi, JPCA, 104, 5631 (2000), is the second order derivatives of E_solvation wrt density matrices. Note: In ddCOSMO, strictly, B is not symmetric. To make it compatible with the CIS framework, it is symmetrized in current implementation. ''' if not self._intermediates or self.grids.coords is None: self.build() mol = self.mol r_vdw = self._intermediates['r_vdw' ] ylm_1sph = self._intermediates['ylm_1sph' ] ui = self._intermediates['ui' ] Lmat = self._intermediates['Lmat' ] A_diele = self._intermediates['A_diele' ] A_inf = self._intermediates['A_inf' ] cached_pol = self._intermediates['cached_pol'] natm = mol.natm nlm = (self.lmax+1)**2 phi = ddcosmo.make_phi(self, dm, r_vdw, ui, ylm_1sph, with_nuc=False) phi = numpy.linalg.solve(A_diele, A_inf.dot(phi.reshape(-1,natm*nlm).T)) Xvec = numpy.linalg.solve(Lmat, phi) Xvec = Xvec.reshape(natm,nlm,-1).transpose(2,0,1) vmat = ddcosmo.make_psi_vmat(self, dm, r_vdw, ui, ylm_1sph, cached_pol, Xvec, Lmat, with_nuc=False)[1] dielectric = self.eps f_epsilon = (dielectric-1.)/dielectric return .5 * f_epsilon * vmat gen_solver = as_solver = gen_ddpcm_solver
[docs] def regularize_xt(self, t, eta, scale=1): return regularize_xt(t, eta)
[docs] def nuc_grad_method(self, grad_method): raise NotImplementedError
DDPCM = ddPCM if __name__ == '__main__': from pyscf import scf mol = gto.M(atom='H 0 0 0; H 0 1 1.2; H 1. .1 0; H .5 .5 1') numpy.random.seed(1) nao = mol.nao_nr() dm = numpy.random.random((nao,nao)) dm = dm + dm.T #dm = scf.RHF(mol).run().make_rdm1() e, vmat = DDPCM(mol).kernel(dm) print(e + 1.2446306643473923) print(lib.fp(vmat) - 0.77873361914445294) 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' #cc-pvdz' mol.build() cm = DDPCM(mol) cm.verbose = 4 mf = ddpcm_for_scf(scf.RHF(mol), cm)#.newton() mf.verbose = 4 mf.kernel() # -75.5697645601958