1#!/usr/bin/env python 2# Copyright 2014-2020 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: Yang Gao <younggao1994@gmail.com> 17 18# 19''' 20Analytical electron-phonon matrix for unrestricted kohn sham 21''' 22import time 23import numpy as np 24from pyscf import lib 25from pyscf.hessian import uks as uks_hess 26from pyscf.hessian import rks as rks_hess 27from pyscf.hessian import rhf as rhf_hess 28from pyscf.grad import rks as rks_grad 29from pyscf.dft import numint 30from pyscf.eph import rhf as rhf_eph 31from pyscf.eph.uhf import uhf_deriv_generator 32from pyscf.data.nist import MP_ME 33 34CUTOFF_FREQUENCY = rhf_eph.CUTOFF_FREQUENCY 35KEEP_IMAG_FREQUENCY = rhf_eph.KEEP_IMAG_FREQUENCY 36 37def _get_vxc_deriv1(hessobj, mo_coeff, mo_occ, max_memory): 38 mol = hessobj.mol 39 mf = hessobj.base 40 if hessobj.grids is not None: 41 grids = hessobj.grids 42 else: 43 grids = mf.grids 44 if grids.coords is None: 45 grids.build(with_non0tab=True) 46 47 nao, nmo = mo_coeff[0].shape 48 ni = mf._numint 49 xctype = ni._xc_type(mf.xc) 50 aoslices = mol.aoslice_by_atom() 51 shls_slice = (0, mol.nbas) 52 ao_loc = mol.ao_loc_nr() 53 dm0a, dm0b = mf.make_rdm1(mo_coeff, mo_occ) 54 55 vmata = np.zeros((mol.natm,3,nao,nao)) 56 vmatb = np.zeros((mol.natm,3,nao,nao)) 57 max_memory = max(2000, max_memory-(vmata.size+vmatb.size)*8/1e6) 58 if xctype == 'LDA': 59 ao_deriv = 1 60 for ao, mask, weight, coords \ 61 in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): 62 rhoa = ni.eval_rho2(mol, ao[0], mo_coeff[0], mo_occ[0], mask, xctype) 63 rhob = ni.eval_rho2(mol, ao[0], mo_coeff[1], mo_occ[1], mask, xctype) 64 vxc, fxc = ni.eval_xc(mf.xc, (rhoa,rhob), 1, deriv=2)[1:3] 65 u_u, u_d, d_d = fxc[0].T 66 67 ao_dm0a = numint._dot_ao_dm(mol, ao[0], dm0a, mask, shls_slice, ao_loc) 68 ao_dm0b = numint._dot_ao_dm(mol, ao[0], dm0b, mask, shls_slice, ao_loc) 69 for ia in range(mol.natm): 70 p0, p1 = aoslices[ia][2:] 71 # First order density = rho1 * 2. *2 is not applied because + c.c. in the end 72 rho1a = np.einsum('xpi,pi->xp', ao[1:,:,p0:p1], ao_dm0a[:,p0:p1]) 73 rho1b = np.einsum('xpi,pi->xp', ao[1:,:,p0:p1], ao_dm0b[:,p0:p1]) 74 75 wv = u_u * rho1a + u_d * rho1b 76 wv *= weight 77 aow = np.einsum('pi,xp->xpi', ao[0], wv) 78 rks_grad._d1_dot_(vmata[ia], mol, aow, ao[0], mask, ao_loc, True) 79 80 wv = u_d * rho1a + d_d * rho1b 81 wv *= weight 82 aow = np.einsum('pi,xp->xpi', ao[0], wv) 83 rks_grad._d1_dot_(vmatb[ia], mol, aow, ao[0], mask, ao_loc, True) 84 ao_dm0a = ao_dm0b = aow = None 85 86 for ia in range(mol.natm): 87 vmata[ia] = -vmata[ia] - vmata[ia].transpose(0,2,1) 88 vmatb[ia] = -vmatb[ia] - vmatb[ia].transpose(0,2,1) 89 90 elif xctype == 'GGA': 91 ao_deriv = 2 92 for ao, mask, weight, coords \ 93 in ni.block_loop(mol, grids, nao, ao_deriv, max_memory): 94 rhoa = ni.eval_rho2(mol, ao[:4], mo_coeff[0], mo_occ[0], mask, xctype) 95 rhob = ni.eval_rho2(mol, ao[:4], mo_coeff[1], mo_occ[1], mask, xctype) 96 vxc, fxc = ni.eval_xc(mf.xc, (rhoa,rhob), 1, deriv=2)[1:3] 97 98 wva, wvb = numint._uks_gga_wv0((rhoa,rhob), vxc, weight) 99 100 ao_dm0a = [numint._dot_ao_dm(mol, ao[i], dm0a, mask, shls_slice, ao_loc) 101 for i in range(4)] 102 ao_dm0b = [numint._dot_ao_dm(mol, ao[i], dm0b, mask, shls_slice, ao_loc) 103 for i in range(4)] 104 for ia in range(mol.natm): 105 wva = dR_rho1a = rks_hess._make_dR_rho1(ao, ao_dm0a, ia, aoslices) 106 wvb = dR_rho1b = rks_hess._make_dR_rho1(ao, ao_dm0b, ia, aoslices) 107 wva[0], wvb[0] = numint._uks_gga_wv1((rhoa,rhob), (dR_rho1a[0],dR_rho1b[0]), vxc, fxc, weight) 108 wva[1], wvb[1] = numint._uks_gga_wv1((rhoa,rhob), (dR_rho1a[1],dR_rho1b[1]), vxc, fxc, weight) 109 wva[2], wvb[2] = numint._uks_gga_wv1((rhoa,rhob), (dR_rho1a[2],dR_rho1b[2]), vxc, fxc, weight) 110 111 aow = np.einsum('npi,Xnp->Xpi', ao[:4], wva) 112 rks_grad._d1_dot_(vmata[ia], mol, aow, ao[0], mask, ao_loc, True) 113 aow = np.einsum('npi,Xnp->Xpi', ao[:4], wvb) 114 rks_grad._d1_dot_(vmatb[ia], mol, aow, ao[0], mask, ao_loc, True) 115 ao_dm0a = ao_dm0b = aow = None 116 117 for ia in range(mol.natm): 118 vmata[ia] = -vmata[ia] - vmata[ia].transpose(0,2,1) 119 vmatb[ia] = -vmatb[ia] - vmatb[ia].transpose(0,2,1) 120 121 elif xctype == 'MGGA': 122 raise NotImplementedError('meta-GGA') 123 124 return vmata, vmatb 125 126def get_eph(ephobj, mo1, omega, vec, mo_rep): 127 if isinstance(mo1, str): 128 mo1 = lib.chkfile.load(mo1, 'scf_mo1') 129 mo1a = mo1['0'] 130 mo1b = mo1['1'] 131 mo1a = dict([(int(k), mo1a[k]) for k in mo1a]) 132 mo1b = dict([(int(k), mo1b[k]) for k in mo1b]) 133 134 mol = ephobj.mol 135 mf = ephobj.base 136 ni = mf._numint 137 ni.libxc.test_deriv_order(mf.xc, 2, raise_error=True) 138 omg, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin) 139 140 vnuc_deriv = ephobj.vnuc_generator(mol) 141 aoslices = mol.aoslice_by_atom() 142 143 mo_coeff, mo_occ = mf.mo_coeff, mf.mo_occ 144 vind = uhf_deriv_generator(mf, mf.mo_coeff, mf.mo_occ) 145 mem_now = lib.current_memory()[0] 146 max_memory = max(2000, mf.max_memory*.9-mem_now) 147 vxc1aoa, vxc1aob = _get_vxc_deriv1(ephobj, mf.mo_coeff, mf.mo_occ, max_memory) 148 nao, nmo = mo_coeff[0].shape 149 mocca = mo_coeff[0][:,mo_occ[0]>0] 150 moccb = mo_coeff[1][:,mo_occ[1]>0] 151 dm0a = np.dot(mocca, mocca.T) 152 dm0b = np.dot(moccb, moccb.T) 153 154 natoms = mol.natm 155 vcorea = [] 156 vcoreb = [] 157 158 for ia in range(natoms): 159 h1 = vnuc_deriv(ia) 160 moia = np.hstack((mo1a[ia], mo1b[ia])) 161 v1 = vind(moia) 162 shl0, shl1, p0, p1 = aoslices[ia] 163 shls_slice = (shl0, shl1) + (0, mol.nbas)*3 164 if abs(hyb)>1e-10: 165 vja, vjb, vka, vkb = \ 166 rhf_hess._get_jk(mol, 'int2e_ip1', 3, 's2kl', 167 ['ji->s2kl', -dm0a[:,p0:p1], #vja 168 'ji->s2kl', -dm0b[:,p0:p1], #vjb 169 'li->s1kj', -dm0a[:,p0:p1], 170 'li->s1kj', -dm0b[:,p0:p1]], #vka 171 shls_slice=shls_slice) 172 vhfa = vja + vjb - hyb * vka 173 vhfb = vjb + vja - hyb * vkb 174 if abs(omg) > 1e-10: 175 with mol.with_range_coulomb(omg): 176 vka, vkb = \ 177 rhf_hess._get_jk(mol, 'int2e_ip1', 3, 's2kl', 178 ['li->s1kj', -dm0a[:,p0:p1], 179 'li->s1kj', -dm0b[:,p0:p1]], # vk1 180 shls_slice=shls_slice) 181 vhfa -= (alpha-hyb) * vka 182 vhfb -= (alpha-hyb) * vkb 183 else: 184 vja, vjb = rhf_hess._get_jk(mol, 'int2e_ip1', 3, 's2kl', 185 ['ji->s2kl', -dm0a[:,p0:p1], 186 'ji->s2kl', -dm0b[:,p0:p1]], # vj1 187 shls_slice=shls_slice) 188 vhfa = vhfb = vja + vjb 189 vtota = h1 + v1[0] + vxc1aoa[ia] + vhfa + vhfa.transpose(0,2,1) 190 vtotb = h1 + v1[1] + vxc1aob[ia] + vhfb + vhfb.transpose(0,2,1) 191 vcorea.append(vtota) 192 vcoreb.append(vtotb) 193 194 vcorea = np.asarray(vcorea).reshape(-1,nao,nao) 195 vcoreb = np.asarray(vcoreb).reshape(-1,nao,nao) 196 197 mass = mol.atom_mass_list() * MP_ME 198 vec = rhf_eph._freq_mass_weighted_vec(vec, omega, mass) 199 mata = np.einsum('xJ,xuv->Juv', vec, vcorea) 200 matb = np.einsum('xJ,xuv->Juv', vec, vcoreb) 201 if mo_rep: 202 mata = np.einsum('Juv,up,vq->Jpq', mata, mf.mo_coeff[0].conj(), mf.mo_coeff[0], optimize=True) 203 matb = np.einsum('Juv,up,vq->Jpq', matb, mf.mo_coeff[1].conj(), mf.mo_coeff[1], optimize=True) 204 return np.asarray([mata,matb]) 205 206 207class EPH(uks_hess.Hessian): 208 '''EPH for unrestricted DFT 209 210 Attributes: 211 cutoff_frequency : float or int 212 cutoff frequency in cm-1. Default is 80 213 keep_imag_frequency : bool 214 Whether to keep imaginary frequencies in the output. Default is False 215 216 Saved results 217 218 omega : numpy.ndarray 219 Vibrational frequencies in au. 220 vec : numpy.ndarray 221 Polarization vectors of the vibration modes 222 eph : numpy.ndarray 223 Electron phonon matrix eph[spin,j,a,b] (j in nmodes, a,b in norbs) 224 ''' 225 226 def __init__(self, scf_method, cutoff_frequency=CUTOFF_FREQUENCY, 227 keep_imag_frequency=KEEP_IMAG_FREQUENCY): 228 uks_hess.Hessian.__init__(self, scf_method) 229 self.cutoff_frequency = cutoff_frequency 230 self.keep_imag_frequency = keep_imag_frequency 231 232 get_mode = rhf_eph.get_mode 233 get_eph = get_eph 234 vnuc_generator = rhf_eph.vnuc_generator 235 kernel = rhf_eph.kernel 236 237if __name__ == '__main__': 238 from pyscf import gto, dft 239 mol = gto.M() 240 mol.atom = [['O', [0.000000000000, 0.000000002577,0.868557119905]], 241 ['H', [0.000000000000,-1.456050381698,2.152719488376]], 242 ['H', [0.000000000000, 1.456050379121,2.152719486067]]] 243 244 mol.unit = 'Bohr' 245 mol.basis = 'sto3g' 246 mol.verbose=4 247 mol.build() # this is a pre-computed relaxed geometry 248 249 mf = dft.UKS(mol) 250 mf.grids.level=6 251 mf.xc = 'b3lyp' 252 mf.conv_tol = 1e-16 253 mf.conv_tol_grad = 1e-10 254 mf.kernel() 255 256 grad = mf.nuc_grad_method().kernel() 257 print("Force on the atoms/au:") 258 print(grad) 259 260 261 myeph = EPH(mf) 262 (epha, ephb), omega = myeph.kernel() 263 264 from pyscf.eph.rks import EPH as REPH 265 mf = dft.RKS(mol) 266 mf.grids.level=6 267 mf.xc = 'b3lyp' 268 mf.conv_tol = 1e-16 269 mf.conv_tol_grad = 1e-10 270 mf.kernel() 271 272 myeph = REPH(mf) 273 rmat, omega = myeph.kernel() 274 print(np.linalg.norm(epha-ephb)) 275 for i in range(len(rmat)): 276 print(min(np.linalg.norm(epha[i]-rmat[i]), np.linalg.norm(epha[i]+rmat[i]))) 277