1#!/usr/bin/env python 2# Copyright 2014-2019 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18 19'''Non-relativistic RKS analytical nuclear gradients''' 20 21 22import numpy 23from pyscf import gto 24from pyscf import lib 25from pyscf.lib import logger 26from pyscf.grad import rhf as rhf_grad 27from pyscf.dft import numint, radi, gen_grid 28from pyscf import __config__ 29 30 31def get_veff(ks_grad, mol=None, dm=None): 32 ''' 33 First order derivative of DFT effective potential matrix (wrt electron coordinates) 34 35 Args: 36 ks_grad : grad.uhf.Gradients or grad.uks.Gradients object 37 ''' 38 if mol is None: mol = ks_grad.mol 39 if dm is None: dm = ks_grad.base.make_rdm1() 40 t0 = (logger.process_clock(), logger.perf_counter()) 41 42 mf = ks_grad.base 43 ni = mf._numint 44 if ks_grad.grids is not None: 45 grids = ks_grad.grids 46 else: 47 grids = mf.grids 48 if grids.coords is None: 49 grids.build(with_non0tab=True) 50 51 if mf.nlc != '': 52 raise NotImplementedError 53 #enabling range-separated hybrids 54 omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin) 55 56 mem_now = lib.current_memory()[0] 57 max_memory = max(2000, ks_grad.max_memory*.9-mem_now) 58 if ks_grad.grid_response: 59 exc, vxc = get_vxc_full_response(ni, mol, grids, mf.xc, dm, 60 max_memory=max_memory, 61 verbose=ks_grad.verbose) 62 logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0)) 63 else: 64 exc, vxc = get_vxc(ni, mol, grids, mf.xc, dm, 65 max_memory=max_memory, verbose=ks_grad.verbose) 66 t0 = logger.timer(ks_grad, 'vxc', *t0) 67 68 if abs(hyb) < 1e-10 and abs(alpha) < 1e-10: 69 vj = ks_grad.get_j(mol, dm) 70 vxc += vj 71 else: 72 vj, vk = ks_grad.get_jk(mol, dm) 73 vk *= hyb 74 if abs(omega) > 1e-10: # For range separated Coulomb operator 75 with mol.with_range_coulomb(omega): 76 vk += ks_grad.get_k(mol, dm) * (alpha - hyb) 77 vxc += vj - vk * .5 78 79 return lib.tag_array(vxc, exc1_grid=exc) 80 81 82def get_vxc(ni, mol, grids, xc_code, dms, relativity=0, hermi=1, 83 max_memory=2000, verbose=None): 84 xctype = ni._xc_type(xc_code) 85 make_rho, nset, nao = ni._gen_rho_evaluator(mol, dms, hermi) 86 ao_loc = mol.ao_loc_nr() 87 88 vmat = numpy.zeros((nset,3,nao,nao)) 89 if xctype == 'LDA': 90 ao_deriv = 1 91 for ao, mask, weight, coords \ 92 in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): 93 for idm in range(nset): 94 rho = make_rho(idm, ao[0], mask, 'LDA') 95 vxc = ni.eval_xc(xc_code, rho, 0, relativity, 1, 96 verbose=verbose)[1] 97 vrho = vxc[0] 98 aow = numpy.einsum('pi,p->pi', ao[0], weight*vrho) 99 _d1_dot_(vmat[idm], mol, ao[1:4], aow, mask, ao_loc, True) 100 rho = vxc = vrho = aow = None 101 elif xctype == 'GGA': 102 ao_deriv = 2 103 for ao, mask, weight, coords \ 104 in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): 105 for idm in range(nset): 106 rho = make_rho(idm, ao[:4], mask, 'GGA') 107 vxc = ni.eval_xc(xc_code, rho, 0, relativity, 1, 108 verbose=verbose)[1] 109 wv = numint._rks_gga_wv0(rho, vxc, weight) 110 _gga_grad_sum_(vmat[idm], mol, ao, wv, mask, ao_loc) 111 rho = vxc = vrho = wv = None 112 elif xctype == 'NLC': 113 raise NotImplementedError('NLC') 114 elif xctype == 'MGGA': 115 raise NotImplementedError('meta-GGA') 116 117 exc = None 118 if nset == 1: 119 vmat = vmat.reshape(3,nao,nao) 120 # - sign because nabla_X = -nabla_x 121 return exc, -vmat 122 123def _make_dR_dao_w(ao, wv): 124 aow = numpy.einsum('npi,p->npi', ao[1:4], wv[0]) 125 # XX, XY, XZ = 4, 5, 6 126 # YX, YY, YZ = 5, 7, 8 127 # ZX, ZY, ZZ = 6, 8, 9 128 aow[0] += numpy.einsum('pi,p->pi', ao[4], wv[1]) # dX nabla_x 129 aow[0] += numpy.einsum('pi,p->pi', ao[5], wv[2]) # dX nabla_y 130 aow[0] += numpy.einsum('pi,p->pi', ao[6], wv[3]) # dX nabla_z 131 aow[1] += numpy.einsum('pi,p->pi', ao[5], wv[1]) # dY nabla_x 132 aow[1] += numpy.einsum('pi,p->pi', ao[7], wv[2]) # dY nabla_y 133 aow[1] += numpy.einsum('pi,p->pi', ao[8], wv[3]) # dY nabla_z 134 aow[2] += numpy.einsum('pi,p->pi', ao[6], wv[1]) # dZ nabla_x 135 aow[2] += numpy.einsum('pi,p->pi', ao[8], wv[2]) # dZ nabla_y 136 aow[2] += numpy.einsum('pi,p->pi', ao[9], wv[3]) # dZ nabla_z 137 return aow 138 139def _d1_dot_(vmat, mol, ao1, ao2, mask, ao_loc, dR1_on_bra=True): 140 shls_slice = (0, mol.nbas) 141 if dR1_on_bra: 142 vmat[0] += numint._dot_ao_ao(mol, ao1[0], ao2, mask, shls_slice, ao_loc) 143 vmat[1] += numint._dot_ao_ao(mol, ao1[1], ao2, mask, shls_slice, ao_loc) 144 vmat[2] += numint._dot_ao_ao(mol, ao1[2], ao2, mask, shls_slice, ao_loc) 145 else: 146 vmat[0] += numint._dot_ao_ao(mol, ao1, ao2[0], mask, shls_slice, ao_loc) 147 vmat[1] += numint._dot_ao_ao(mol, ao1, ao2[1], mask, shls_slice, ao_loc) 148 vmat[2] += numint._dot_ao_ao(mol, ao1, ao2[2], mask, shls_slice, ao_loc) 149 150def _gga_grad_sum_(vmat, mol, ao, wv, mask, ao_loc): 151 aow = numpy.einsum('npi,np->pi', ao[:4], wv) 152 _d1_dot_(vmat, mol, ao[1:4], aow, mask, ao_loc, True) 153 aow = _make_dR_dao_w(ao, wv) 154 _d1_dot_(vmat, mol, aow, ao[0], mask, ao_loc, True) 155 return vmat 156 157 158def get_vxc_full_response(ni, mol, grids, xc_code, dms, relativity=0, hermi=1, 159 max_memory=2000, verbose=None): 160 '''Full response including the response of the grids''' 161 xctype = ni._xc_type(xc_code) 162 make_rho, nset, nao = ni._gen_rho_evaluator(mol, dms, hermi) 163 ao_loc = mol.ao_loc_nr() 164 165 excsum = 0 166 vmat = numpy.zeros((3,nao,nao)) 167 if xctype == 'LDA': 168 ao_deriv = 1 169 vtmp = numpy.empty((3,nao,nao)) 170 for atm_id, (coords, weight, weight1) in enumerate(grids_response_cc(grids)): 171 mask = gen_grid.make_mask(mol, coords) 172 ao = ni.eval_ao(mol, coords, deriv=ao_deriv, non0tab=mask) 173 rho = make_rho(0, ao[0], mask, 'LDA') 174 exc, vxc = ni.eval_xc(xc_code, rho, 0, relativity, 1, 175 verbose=verbose)[:2] 176 vrho = vxc[0] 177 178 vtmp = numpy.zeros((3,nao,nao)) 179 aow = numpy.einsum('pi,p->pi', ao[0], weight*vrho) 180 _d1_dot_(vtmp, mol, ao[1:4], aow, mask, ao_loc, True) 181 vmat += vtmp 182 183 # response of weights 184 excsum += numpy.einsum('r,r,nxr->nx', exc, rho, weight1) 185 # response of grids coordinates 186 excsum[atm_id] += numpy.einsum('xij,ji->x', vtmp, dms) * 2 187 rho = vxc = vrho = aow = None 188 189 elif xctype == 'GGA': 190 ao_deriv = 2 191 for atm_id, (coords, weight, weight1) in enumerate(grids_response_cc(grids)): 192 mask = gen_grid.make_mask(mol, coords) 193 ao = ni.eval_ao(mol, coords, deriv=ao_deriv, non0tab=mask) 194 rho = make_rho(0, ao[:4], mask, 'GGA') 195 exc, vxc = ni.eval_xc(xc_code, rho, 0, relativity, 1, 196 verbose=verbose)[:2] 197 198 vtmp = numpy.zeros((3,nao,nao)) 199 wv = numint._rks_gga_wv0(rho, vxc, weight) 200 _gga_grad_sum_(vtmp, mol, ao, wv, mask, ao_loc) 201 vmat += vtmp 202 203 # response of weights 204 excsum += numpy.einsum('r,r,nxr->nx', exc, rho[0], weight1) 205 # response of grids coordinates 206 excsum[atm_id] += numpy.einsum('xij,ji->x', vtmp, dms) * 2 207 rho = vxc = vrho = wv = None 208 209 elif xctype == 'NLC': 210 raise NotImplementedError('NLC') 211 elif xctype == 'MGGA': 212 raise NotImplementedError('meta-GGA') 213 214 # - sign because nabla_X = -nabla_x 215 return excsum, -vmat 216 217 218# JCP 98, 5612 (1993); DOI:10.1063/1.464906 219def grids_response_cc(grids): 220 mol = grids.mol 221 atom_grids_tab = grids.gen_atomic_grids(mol, grids.atom_grid, 222 grids.radi_method, 223 grids.level, grids.prune) 224 atm_coords = numpy.asarray(mol.atom_coords() , order='C') 225 atm_dist = gto.inter_distance(mol, atm_coords) 226 227 def _radii_adjust(mol, atomic_radii): 228 charges = mol.atom_charges() 229 if grids.radii_adjust == radi.treutler_atomic_radii_adjust: 230 rad = numpy.sqrt(atomic_radii[charges]) + 1e-200 231 elif grids.radii_adjust == radi.becke_atomic_radii_adjust: 232 rad = atomic_radii[charges] + 1e-200 233 else: 234 fadjust = lambda i, j, g: g 235 gadjust = lambda *args: 1 236 return fadjust, gadjust 237 238 rr = rad.reshape(-1,1) * (1./rad) 239 a = .25 * (rr.T - rr) 240 a[a<-.5] = -.5 241 a[a>0.5] = 0.5 242 243 def fadjust(i, j, g): 244 return g + a[i,j]*(1-g**2) 245 246 #: d[g + a[i,j]*(1-g**2)] /dg = 1 - 2*a[i,j]*g 247 def gadjust(i, j, g): 248 return 1 - 2*a[i,j]*g 249 return fadjust, gadjust 250 251 fadjust, gadjust = _radii_adjust(mol, grids.atomic_radii) 252 253 def gen_grid_partition(coords, atom_id): 254 ngrids = coords.shape[0] 255 grid_dist = [] 256 grid_norm_vec = [] 257 for ia in range(mol.natm): 258 v = (atm_coords[ia] - coords).T 259 normv = numpy.linalg.norm(v,axis=0) + 1e-200 260 v /= normv 261 grid_dist.append(normv) 262 grid_norm_vec.append(v) 263 264 def get_du(ia, ib): # JCP 98, 5612 (1993); (B10) 265 uab = atm_coords[ia] - atm_coords[ib] 266 duab = 1./atm_dist[ia,ib] * grid_norm_vec[ia] 267 duab-= uab[:,None]/atm_dist[ia,ib]**3 * (grid_dist[ia]-grid_dist[ib]) 268 return duab 269 270 pbecke = numpy.ones((mol.natm,ngrids)) 271 dpbecke = numpy.zeros((mol.natm,mol.natm,3,ngrids)) 272 for ia in range(mol.natm): 273 for ib in range(ia): 274 g = 1/atm_dist[ia,ib] * (grid_dist[ia]-grid_dist[ib]) 275 p0 = fadjust(ia, ib, g) 276 p1 = (3 - p0**2) * p0 * .5 277 p2 = (3 - p1**2) * p1 * .5 278 p3 = (3 - p2**2) * p2 * .5 279 t_uab = 27./16 * (1-p2**2) * (1-p1**2) * (1-p0**2) * gadjust(ia, ib, g) 280 281 s_uab = .5 * (1 - p3 + 1e-200) 282 s_uba = .5 * (1 + p3 + 1e-200) 283 pbecke[ia] *= s_uab 284 pbecke[ib] *= s_uba 285 pt_uab =-t_uab / s_uab 286 pt_uba = t_uab / s_uba 287 288# * When grid is on atom ia/ib, ua/ub == 0, d_uba/d_uab may have huge error 289# How to remove this error? 290 duab = get_du(ia, ib) 291 duba = get_du(ib, ia) 292 if ia == atom_id: 293 dpbecke[ia,ia] += pt_uab * duba 294 dpbecke[ia,ib] += pt_uba * duba 295 else: 296 dpbecke[ia,ia] += pt_uab * duab 297 dpbecke[ia,ib] += pt_uba * duab 298 299 if ib == atom_id: 300 dpbecke[ib,ib] -= pt_uba * duab 301 dpbecke[ib,ia] -= pt_uab * duab 302 else: 303 dpbecke[ib,ib] -= pt_uba * duba 304 dpbecke[ib,ia] -= pt_uab * duba 305 306# * JCP 98, 5612 (1993); (B8) (B10) miss many terms 307 if ia != atom_id and ib != atom_id: 308 ua_ub = grid_norm_vec[ia] - grid_norm_vec[ib] 309 ua_ub /= atm_dist[ia,ib] 310 dpbecke[atom_id,ia] -= pt_uab * ua_ub 311 dpbecke[atom_id,ib] -= pt_uba * ua_ub 312 313 for ia in range(mol.natm): 314 dpbecke[:,ia] *= pbecke[ia] 315 316 return pbecke, dpbecke 317 318 natm = mol.natm 319 for ia in range(natm): 320 coords, vol = atom_grids_tab[mol.atom_symbol(ia)] 321 coords = coords + atm_coords[ia] 322 pbecke, dpbecke = gen_grid_partition(coords, ia) 323 z = 1./pbecke.sum(axis=0) 324 w1 = dpbecke[:,ia] * z 325 w1 -= pbecke[ia] * z**2 * dpbecke.sum(axis=1) 326 w1 *= vol 327 w0 = vol * pbecke[ia] * z 328 yield coords, w0, w1 329 330 331class Gradients(rhf_grad.Gradients): 332 333 # This parameter has no effects for HF gradients. Add this attribute so that 334 # the kernel function can be reused in the DFT gradients code. 335 grid_response = getattr(__config__, 'grad_rks_Gradients_grid_response', False) 336 337 def __init__(self, mf): 338 rhf_grad.Gradients.__init__(self, mf) 339 self.grids = None 340 # This parameter has no effects for HF gradients. Add this attribute so that 341 # the kernel function can be reused in the DFT gradients code. 342 self.grid_response = False 343 self._keys = self._keys.union(['grid_response', 'grids']) 344 345 def dump_flags(self, verbose=None): 346 rhf_grad.Gradients.dump_flags(self, verbose) 347 logger.info(self, 'grid_response = %s', self.grid_response) 348 #if callable(self.base.grids.prune): 349 # logger.info(self, 'Grid pruning %s may affect DFT gradients accuracy.' 350 # 'Call mf.grids.run(prune=False) to mute grid pruning', 351 # self.base.grids.prune) 352 return self 353 354 get_veff = get_veff 355 356 def extra_force(self, atom_id, envs): 357 '''Hook for extra contributions in analytical gradients. 358 359 Contributions like the response of auxiliary basis in density fitting 360 method, the grid response in DFT numerical integration can be put in 361 this function. 362 ''' 363 if self.grid_response: 364 vhf = envs['vhf'] 365 log = envs['log'] 366 log.debug('grids response for atom %d %s', 367 atom_id, vhf.exc1_grid[atom_id]) 368 return vhf.exc1_grid[atom_id] 369 else: 370 return 0 371 372Grad = Gradients 373 374from pyscf import dft 375dft.rks.RKS.Gradients = dft.rks_symm.RKS.Gradients = lib.class_as_method(Gradients) 376 377 378if __name__ == '__main__': 379 from pyscf import dft 380 381 mol = gto.Mole() 382 mol.atom = [ 383 ['O' , (0. , 0. , 0.)], 384 [1 , (0. , -0.757 , 0.587)], 385 [1 , (0. , 0.757 , 0.587)] ] 386 mol.basis = '631g' 387 mol.build() 388 mf = dft.RKS(mol) 389 mf.conv_tol = 1e-14 390 #mf.grids.atom_grid = (20,86) 391 e0 = mf.scf() 392 g = mf.Gradients() 393 print(lib.finger(g.kernel()) - -0.049887865971659243) 394#[[ -4.20040265e-16 -6.59462771e-16 2.10150467e-02] 395# [ 1.42178271e-16 2.81979579e-02 -1.05137653e-02] 396# [ 6.34069238e-17 -2.81979579e-02 -1.05137653e-02]] 397 g.grid_response = True 398 print(lib.finger(g.kernel()) - -0.049891265876709084) 399# O 0.0000000000 -0.0000000000 0.0210225191 400# H 0.0000000000 0.0281984036 -0.0105112595 401# H -0.0000000000 -0.0281984036 -0.0105112595 402 403 mf.xc = 'b88,p86' 404 e0 = mf.scf() 405 g = Gradients(mf) 406 print(lib.finger(g.kernel()) - -0.050382923259300716) 407#[[ -8.20194970e-16 -2.04319288e-15 2.44405835e-02] 408# [ 4.36709255e-18 2.73690416e-02 -1.22232039e-02] 409# [ 3.44483899e-17 -2.73690416e-02 -1.22232039e-02]] 410 g.grid_response = True 411 print(lib.finger(g.kernel()) - -0.05036316927480719) 412 413 mf.xc = 'b3lypg' 414 e0 = mf.scf() 415 g = Gradients(mf) 416 print(lib.finger(g.kernel()) - -0.035613964330885352) 417#[[ -3.59411142e-16 -2.68753987e-16 1.21557501e-02] 418# [ 4.04977877e-17 2.11112794e-02 -6.08181640e-03] 419# [ 1.52600378e-16 -2.11112794e-02 -6.08181640e-03]] 420 421 422 mol = gto.Mole() 423 mol.atom = [ 424 ['H' , (0. , 0. , 1.804)], 425 ['F' , (0. , 0. , 0. )], ] 426 mol.unit = 'B' 427 mol.basis = '631g' 428 mol.build() 429 430 mf = dft.RKS(mol) 431 mf.conv_tol = 1e-14 432 mf.kernel() 433 print(lib.finger(Gradients(mf).kernel()) - 0.0018831588319051444) 434# sum over z direction non-zero, due to meshgrid response 435#[[ 0 0 -2.68934738e-03] 436# [ 0 0 2.69333577e-03]] 437 mf = dft.RKS(mol) 438 mf.grids.prune = None 439 mf.grids.level = 6 440 mf.conv_tol = 1e-14 441 mf.kernel() 442 print(lib.finger(Gradients(mf).kernel()) - 0.0018819497229394144) 443#[[ 0 0 -2.68931547e-03] 444# [ 0 0 2.68911282e-03]] 445