#!/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: Yang Gao <younggao1994@gmail.com>
#
'''
Analytical electron-phonon matrix for restricted hartree fock
'''
import numpy as np
import scipy.linalg
from pyscf.hessian import rhf
from pyscf.lib import logger, chkfile
from pyscf.scf._response_functions import _gen_rhf_response
from pyscf import __config__
from pyscf.data.nist import HARTREE2WAVENUMBER, MP_ME
CUTOFF_FREQUENCY = getattr(__config__, 'eph_cutoff_frequency', 80) # 80 cm-1
KEEP_IMAG_FREQUENCY = getattr(__config__, 'eph_keep_imaginary_frequency', False)
IMAG_CUTOFF_FREQUENCY = getattr(__config__, 'eph_imag_cutoff_frequency', 1e-4)
[docs]
def kernel(ephobj, mo_energy=None, mo_coeff=None, mo_occ=None, mo_rep=False):
if mo_energy is None: mo_energy = ephobj.base.mo_energy
if mo_coeff is None: mo_coeff = ephobj.base.mo_coeff
if mo_occ is None: mo_occ = ephobj.base.mo_occ
# chkfile is used to pass first orbitals from hessian methods to eph methods
# TODO: Remove the dependence to chfile and return first orbitals from a function
assert ephobj.chkfile is not None, 'chkfile is required to save first order orbitals'
de = ephobj.hess_elec(mo_energy, mo_coeff, mo_occ)
ephobj.de = de + ephobj.hess_nuc(ephobj.mol)
omega, vec = ephobj.get_mode(ephobj.mol, ephobj.de)
ephobj.omega, ephobj.vec = omega, vec
ephobj.eph = ephobj.get_eph(ephobj.chkfile, omega, vec, mo_rep)
return ephobj.eph, ephobj.omega
[docs]
def solve_hmat(mol, hmat, cutoff_frequency=CUTOFF_FREQUENCY,
keep_imag_frequency=KEEP_IMAG_FREQUENCY):
log = logger.new_logger(mol, mol.verbose)
mass = mol.atom_mass_list() * MP_ME
natom = len(mass)
h = np.empty_like(hmat) #[atom, axis, atom, axis]
for i in range(natom):
for j in range(natom):
h[i,j] = hmat[i,j] / np.sqrt(mass[i]*mass[j])
forcemat = h.transpose(0,2,1,3).reshape(natom*3, natom*3)
forcemat[abs(forcemat)<1e-12]=0 #improve stability
w, c = scipy.linalg.eig(forcemat)
idx = np.argsort(w.real)[::-1] # sort the mode of the frequency
w = w[idx]
c = c[:,idx]
w_au = w**0.5
w_cm = w_au * HARTREE2WAVENUMBER
log.info('****Eigenmodes(cm-1)****')
for i, omega in enumerate(w_cm):
if abs(omega.imag) < IMAG_CUTOFF_FREQUENCY:
w_au[i] = w_au[i].real
w_cm[i] = w_cm[i].real
if omega.real > cutoff_frequency:
log.info("Mode %i Omega=%.4f", i, omega.real)
else:
log.info("Mode %i Omega=%.4f, mode filtered", i, omega.real)
else:
log.info("Mode %i Omega=%.4fj, imaginary mode", i, omega.imag)
if KEEP_IMAG_FREQUENCY:
idx_real = np.where(w_cm.real>cutoff_frequency)[0]
idx_imag = np.where(abs(w_cm.imag)>IMAG_CUTOFF_FREQUENCY)[0]
idx = np.concatenate([idx_real, idx_imag])
else:
w_au = w_au.real
idx = np.where(w_cm.real>cutoff_frequency)[0]
w_new = w_au[idx]
c_new = c[:,idx]
log.info('****Remaining Eigenmodes(cm-1)****')
for i, omega in enumerate(w_cm[idx]):
if omega.imag == 0:
log.info("Mode %i Omega=%.4f", i, omega.real)
else:
log.info("Mode %i Omega=%.4fj", i, omega.imag)
return w_new, c_new
[docs]
def get_mode(ephobj, mol=None, de=None):
if mol is None: mol = ephobj.mol
if de is None:
if ephobj.de is None:
de = ephobj.hess_elec() + ephobj.hess_nuc()
else:
de = ephobj.de
return solve_hmat(mol, de, ephobj.cutoff_frequency, ephobj.keep_imag_frequency)
[docs]
def rhf_deriv_generator(mf, mo_coeff, mo_occ):
nao, nmo = mo_coeff.shape
mocc = mo_coeff[:,mo_occ>0]
nocc = mocc.shape[1]
vresp = _gen_rhf_response(mf, mo_coeff, mo_occ, hermi=1)
def fx(mo1):
mo1 = mo1.reshape(-1,nmo,nocc)
nset = len(mo1)
dm1 = np.empty((nset,nao,nao))
for i, x in enumerate(mo1):
dm = np.dot(x*2, mocc.T) # *2 for double occupancy
dm1[i] = dm + dm.T
v1 = vresp(dm1)
return v1
return fx
[docs]
def vnuc_generator(ephobj, mol):
if mol is None: mol = ephobj.mol
aoslices = mol.aoslice_by_atom()
def vnuc_deriv(atm_id):
shl0, shl1, p0, p1 = aoslices[atm_id]
with mol.with_rinv_at_nucleus(atm_id):
vrinv = mol.intor('int1e_iprinv', comp=3) # <\nabla|1/r|>
vrinv *= -mol.atom_charge(atm_id)
return vrinv + vrinv.transpose(0,2,1)
return vnuc_deriv
def _freq_mass_weighted_vec(vec, omega, mass):
nmodes, natoms = len(omega), len(mass)
dtype = np.result_type(omega, vec)
vec = vec.reshape(natoms,3,nmodes).astype(dtype)
for i in range(natoms):
for j in range(nmodes):
vec[i,:,j] /= np.sqrt(2*mass[i]*omega[j])
vec = vec.reshape(3*natoms,nmodes)
return vec
[docs]
def get_eph(ephobj, mo1, omega, vec, mo_rep):
if isinstance(mo1, str):
mo1 = chkfile.load(mo1, 'scf_mo1')
mo1 = dict([(int(k), mo1[k]) for k in mo1])
mol = ephobj.mol
mf = ephobj.base
vnuc_deriv = ephobj.vnuc_generator(mol)
aoslices = mol.aoslice_by_atom()
vind = rhf_deriv_generator(mf, mf.mo_coeff, mf.mo_occ)
mocc = mf.mo_coeff[:,mf.mo_occ>0]
dm0 = np.dot(mocc, mocc.T) * 2
natoms = mol.natm
nao = mol.nao_nr()
vcore = []
for ia in range(natoms):
h1 = vnuc_deriv(ia)
v1 = vind(mo1[ia])
shl0, shl1, p0, p1 = aoslices[ia]
shls_slice = (shl0, shl1) + (0, mol.nbas)*3
vj1, vk1= rhf._get_jk(mol, 'int2e_ip1', 3, 's2kl',
['ji->s2kl', -dm0[:,p0:p1], # vj1
'li->s1kj', -dm0[:,p0:p1]], # vk1
shls_slice=shls_slice)
vhf = vj1 - vk1*.5
vtot = h1 + v1 + vhf + vhf.transpose(0,2,1)
vcore.append(vtot)
vcore = np.asarray(vcore).reshape(-1,nao,nao)
mass = mol.atom_mass_list() * MP_ME
vec = _freq_mass_weighted_vec(vec, omega, mass)
mat = np.einsum('xJ,xuv->Juv', vec, vcore)
if mo_rep:
mat = np.einsum('Juv,up,vq->Jpq', mat, mf.mo_coeff.conj(), mf.mo_coeff, optimize=True)
return mat
[docs]
class EPH(rhf.Hessian):
'''EPH for restricted Hartree Fock
Attributes:
cutoff_frequency : float or int
cutoff frequency in cm-1. Default is 80
keep_imag_frequency : bool
Whether to keep imaginary frequencies in the output. Default is False
Saved results
omega : numpy.ndarray
Vibrational frequencies in au.
vec : numpy.ndarray
Polarization vectors of the vibration modes
eph : numpy.ndarray
Electron phonon matrix eph[j,a,b] (j in nmodes, a,b in norbs)
'''
def __init__(self, scf_method, cutoff_frequency=CUTOFF_FREQUENCY,
keep_imag_frequency=KEEP_IMAG_FREQUENCY):
rhf.Hessian.__init__(self, scf_method)
self.cutoff_frequency = cutoff_frequency
self.keep_imag_frequency = keep_imag_frequency
get_mode = get_mode
get_eph = get_eph
vnuc_generator = vnuc_generator
kernel = kernel
if __name__ == '__main__':
from pyscf import gto, scf
mol = gto.M()
mol.atom = [['O', [0.000000000000, -0.000000000775, 0.923671924285]],
['H', [-0.000000000000, -1.432564848017, 2.125164039823]],
['H', [0.000000000000, 1.432564848792, 2.125164035930]]]
mol.unit = 'Bohr'
mol.basis = 'sto3g'
mol.verbose=4
mol.build() # this is a pre-computed relaxed geometry
mf = scf.RHF(mol)
mf.conv_tol = 1e-16
mf.conv_tol_grad = 1e-10
mf.kernel()
myeph = EPH(mf)
grad = mf.nuc_grad_method().kernel()
print("Force on the atoms/au:")
print(grad)
eph, omega = myeph.kernel(mo_rep=True)
print(np.amax(eph))