diff --git a/gpu4pyscf/properties/__init__.py b/gpu4pyscf/properties/__init__.py
new file mode 100644
index 00000000..cf49865c
--- /dev/null
+++ b/gpu4pyscf/properties/__init__.py
@@ -0,0 +1,16 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+from gpu4pyscf.properties import polarizability, ir, shielding
\ No newline at end of file
diff --git a/gpu4pyscf/properties/ir.py b/gpu4pyscf/properties/ir.py
new file mode 100644
index 00000000..9f813df9
--- /dev/null
+++ b/gpu4pyscf/properties/ir.py
@@ -0,0 +1,129 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+from functools import reduce
+from pyscf.hessian import thermo
+import numpy as np
+import cupy
+from pyscf.data import elements, nist
+from scipy.constants import physical_constants
+from gpu4pyscf.lib import logger
+from gpu4pyscf.lib.cupy_helper import contract
+
+LINDEP_THRESHOLD = 1e-7
+
+
+def eval_ir_freq_intensity(mf, hessian_obj):
+ """main function to calculate the polarizability
+
+ Args:
+ mf: mean field object
+ unit (str, optional): the unit of the polarizability. Defaults to 'au'.
+
+ Returns:
+ polarizability (numpy.array): polarizability
+ """
+ log = logger.new_logger(hessian_obj, mf.mol.verbose)
+ hessian = hessian_obj.kernel()
+ hartree_kj = nist.HARTREE2J*1e3
+ unit2cm = ((hartree_kj * nist.AVOGADRO)**.5 / (nist.BOHR*1e-10)
+ / (2*np.pi*nist.LIGHT_SPEED_SI) * 1e-2)
+ natm = mf.mol.natm
+ nao = mf.mol.nao
+ dm0 = mf.make_rdm1()
+
+ atom_charges = mf.mol.atom_charges()
+ mass = cupy.array([elements.MASSES[atom_charges[i]] for i in range(natm)])
+ # hessian_mass = contract('ijkl,i,j->ijkl', hessian,
+ # 1/cupy.sqrt(mass), 1/cupy.sqrt(mass))
+ hessian_mass = contract('ijkl,i->ijkl', cupy.array(hessian), 1/cupy.sqrt(mass))
+ hessian_mass = contract('ijkl,j->ijkl', hessian_mass, 1/cupy.sqrt(mass))
+
+ TR = thermo._get_TR(mass.get(), mf.mol.atom_coords())
+ TRspace = []
+ TRspace.append(TR[:3])
+
+ rot_const = thermo.rotation_const(mass.get(), mf.mol.atom_coords())
+ rotor_type = thermo._get_rotor_type(rot_const)
+ if rotor_type == 'ATOM':
+ pass
+ elif rotor_type == 'LINEAR': # linear molecule
+ TRspace.append(TR[3:5])
+ else:
+ TRspace.append(TR[3:])
+
+ if TRspace:
+ TRspace = cupy.vstack(TRspace)
+ q, r = cupy.linalg.qr(TRspace.T)
+ P = cupy.eye(natm * 3) - q.dot(q.T)
+ w, v = cupy.linalg.eigh(P)
+ bvec = v[:,w > LINDEP_THRESHOLD]
+ h = reduce(cupy.dot, (bvec.T, hessian_mass.transpose(0,2,1,3).reshape(3*natm,3*natm), bvec))
+ e, mode = cupy.linalg.eigh(h)
+ mode = bvec.dot(mode)
+
+ c = contract('ixn,i->ixn', mode.reshape(natm, 3, -1),
+ 1/np.sqrt(mass)).reshape(3*natm, -1)
+ freq = cupy.sign(e)*cupy.sqrt(cupy.abs(e))*unit2cm
+
+ mo_coeff = mf.mo_coeff
+ mo_occ = mf.mo_occ
+ mo_energy = mf.mo_energy
+ mo_coeff = cupy.array(mo_coeff)
+ mo_occ = cupy.array(mo_occ)
+ mo_energy = cupy.array(mo_energy)
+ mocc = mo_coeff[:, mo_occ > 0]
+ mocc = cupy.array(mocc)
+
+ atmlst = range(natm)
+ h1ao = hessian_obj.make_h1(mo_coeff, mo_occ, None, atmlst)
+ # TODO: compact with hessian method, which can save one time cphf solve.
+ # ! Different from PySCF, mo1 is all in mo!
+ mo1, mo_e1 = hessian_obj.solve_mo1(mo_energy, mo_coeff, mo_occ, h1ao,
+ None, atmlst, hessian_obj.max_memory, log)
+
+
+ tmp = cupy.empty((3, 3, natm)) # dipole moment, x,y,z
+ aoslices = mf.mol.aoslice_by_atom()
+ with mf.mol.with_common_orig((0, 0, 0)):
+ hmuao = mf.mol.intor('int1e_r') # mu
+ hmuao11 = -mf.mol.intor('int1e_irp').reshape(3, 3, nao, nao)
+ hmuao = cupy.array(hmuao)
+ hmuao11 = cupy.array(hmuao11)
+ for i0, ia in enumerate(atmlst):
+ shl0, shl1, p0, p1 = aoslices[ia]
+ h11ao = cupy.zeros((3, 3, nao, nao))
+
+ h11ao[:, :, :, p0:p1] += hmuao11[:, :, :, p0:p1]
+ h11ao[:, :, p0:p1] += hmuao11[:, :, :, p0:p1].transpose(0, 1, 3, 2)
+
+ tmp0 = contract('ypi,vi->ypv', mo1[ia], mocc) # nabla
+ dm1 = contract('ypv,up->yuv', tmp0, mo_coeff)
+ tmp[:, :, ia] = -contract('xuv,yuv->xy', hmuao, dm1) * 4 #the minus means the density should be negative, but calculate it is positive.
+ tmp[:, :, ia] -= contract('xyuv,vu->xy', h11ao, dm0)
+ tmp[:, :, ia] += mf.mol.atom_charge(ia)*cupy.eye(3)
+
+ alpha = physical_constants["fine-structure constant"][0]
+ amu = physical_constants["atomic mass constant"][0]
+ m_e = physical_constants["electron mass"][0]
+ N_A = physical_constants["Avogadro constant"][0]
+ a_0 = physical_constants["Bohr radius"][0]
+ unit_kmmol = alpha**2 * (1e-3 / amu) * m_e * N_A * np.pi * a_0 / 3
+
+ intensity = contract('xym,myn->xn', tmp, c.reshape(natm, 3, -1))
+ intensity = contract('xn,xn->n', intensity, intensity)
+ intensity *= unit_kmmol
+
+ return freq, intensity
diff --git a/gpu4pyscf/properties/polarizability.py b/gpu4pyscf/properties/polarizability.py
new file mode 100644
index 00000000..567f71c1
--- /dev/null
+++ b/gpu4pyscf/properties/polarizability.py
@@ -0,0 +1,92 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import numpy as np
+from gpu4pyscf.scf import cphf
+import cupy
+from gpu4pyscf.lib.cupy_helper import contract
+
+
+def gen_vind(mf, mo_coeff, mo_occ):
+ """get the induced potential. This is the same as contract the mo1 with the kernel.
+
+ Args:
+ mf: mean field object
+ mo_coeff (numpy.array): mo coefficients
+ mo_occ (numpy.array): mo_coefficients
+
+ Returns:
+ fx (function): a function to calculate the induced potential with the input as the mo1.
+ """
+ nao, nmo = mo_coeff.shape
+ mocc = mo_coeff[:, mo_occ > 0]
+ mvir = mo_coeff[:, mo_occ == 0]
+ nocc = mocc.shape[1]
+ nvir = nmo - nocc
+ vresp = mf.gen_response(mo_coeff, mo_occ, hermi=1)
+
+ def fx(mo1):
+ mo1 = mo1.reshape(-1, nvir, nocc) # * the saving pattern
+ mo1_mo_real = contract('nai,ua->nui', mo1, mvir)
+ dm1 = 2*contract('nui,vi->nuv', mo1_mo_real, mocc.conj())
+ dm1+= dm1.transpose(0,2,1)
+
+ v1 = vresp(dm1) # (nset, nao, nao)
+ tmp = contract('nuv,vi->nui', v1, mocc)
+ v1vo = contract('nui,ua->nai', tmp, mvir.conj())
+
+ return v1vo
+ return fx
+
+
+def eval_polarizability(mf, unit='au'):
+ """main function to calculate the polarizability
+
+ Args:
+ mf: mean field object
+ unit (str, optional): the unit of the polarizability. Defaults to 'au'.
+
+ Returns:
+ polarizability (numpy.array): polarizability
+ """
+
+ polarizability = np.empty((3, 3))
+
+ mo_coeff = mf.mo_coeff
+ mo_occ = mf.mo_occ
+ mo_energy = mf.mo_energy
+ mo_coeff = cupy.array(mo_coeff)
+ mo_occ = cupy.array(mo_occ)
+ mo_energy = cupy.array(mo_energy)
+ fx = gen_vind(mf, mo_coeff, mo_occ)
+ mocc = mo_coeff[:, mo_occ > 0]
+ mvir = mo_coeff[:, mo_occ == 0]
+
+ with mf.mol.with_common_orig((0, 0, 0)):
+ h1 = mf.mol.intor('int1e_r')
+ h1 = cupy.array(h1)
+ for idirect in range(3):
+ h1ai = -mvir.T.conj()@h1[idirect]@mocc
+ mo1 = cphf.solve(fx, mo_energy, mo_occ, h1ai, max_cycle=20, tol=1e-10)[0]
+ for jdirect in range(idirect, 3):
+ p10 = np.trace(mo1.conj().T@mvir.conj().T@h1[jdirect]@mocc)*2
+ p01 = np.trace(mocc.conj().T@h1[jdirect]@mvir@mo1)*2
+ polarizability[idirect, jdirect] = p10+p01
+ polarizability[1, 0] = polarizability[0, 1]
+ polarizability[2, 0] = polarizability[0, 2]
+ polarizability[2, 1] = polarizability[1, 2]
+
+ return polarizability
+
diff --git a/gpu4pyscf/properties/shielding.py b/gpu4pyscf/properties/shielding.py
new file mode 100644
index 00000000..dbbeccc6
--- /dev/null
+++ b/gpu4pyscf/properties/shielding.py
@@ -0,0 +1,263 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+from gpu4pyscf.scf import cphf
+import numpy as np
+from pyscf.data import nist
+import cupy
+from pyscf.scf import _vhf, jk
+from gpu4pyscf.dft import numint
+from gpu4pyscf.lib.cupy_helper import contract, take_last2d, add_sparse
+
+
+def gen_vind(mf, mo_coeff, mo_occ):
+ """get the induced potential. This is the same as contract the mo1 with the kernel.
+
+ Args:
+ mf: mean field object
+ mo_coeff (numpy.array): mo coefficients
+ mo_occ (numpy.array): mo_coefficients
+
+ Returns:
+ fx (function): a function to calculate the induced potential with the input as the mo1.
+ """
+ nao, nmo = mo_coeff.shape
+ mocc = mo_coeff[:, mo_occ > 0]
+ mvir = mo_coeff[:, mo_occ == 0]
+ nocc = mocc.shape[1]
+ nvir = nmo - nocc
+ omega, alpha, hyb = mf._numint.rsh_and_hybrid_coeff(
+ mf.xc, spin=mf.mol.spin)
+
+ def fx(mo1):
+ mo1 = mo1.reshape(-1, nvir, nocc) # * the saving pattern
+ mo1_mo_real = contract('nai,ua->nui', mo1, mvir)
+ dm1 = 2*contract('nui,vi->nuv', mo1_mo_real, mocc.conj())
+ dm1 -= dm1.transpose(0, 2, 1)
+ if hasattr(mf,'with_df'):
+ v1 = cupy.empty((3, nao, nao))
+ for i in range(3):
+ v1[i] =+mf.get_jk(mf.mol, dm1[i], hermi=2, with_j=False)[1]*0.5*hyb
+ else:
+ v1 = np.empty((3, nao, nao))
+ for i in range(3):
+ v1[i] = -jk.get_jk(mf.mol, dm1[i].get(), 'ijkl,jk->il')*0.5*hyb
+ v1 = cupy.array(v1)
+ tmp = contract('nuv,vi->nui', v1, mocc)
+ v1vo = contract('nui,ua->nai', tmp, mvir.conj())
+
+ return v1vo
+
+ return fx
+
+
+def nr_rks(ni, mol, grids, xc_code, dms):
+
+ xctype = ni._xc_type(xc_code)
+ mo_coeff = getattr(dms, 'mo_coeff', None)
+ mo_occ = getattr(dms, 'mo_occ', None)
+ nao = mo_coeff.shape[1]
+
+ opt = getattr(ni, 'gdftopt', None)
+ if opt is None:
+ ni.build(mol, grids.coords)
+ opt = ni.gdftopt
+ coeff = cupy.asarray(opt.coeff)
+ nao, nao0 = coeff.shape
+ dms = cupy.asarray(dms).reshape(-1,nao0,nao0)
+ dms = take_last2d(dms, opt.ao_idx)
+ mo_coeff = mo_coeff[opt.ao_idx]
+
+ vmat = cupy.zeros((3, nao, nao))
+ if xctype == 'LDA':
+ ao_deriv = 0
+ else:
+ ao_deriv = 1
+
+ for ao, index, weight, coords in ni.block_loop(opt.mol, grids, nao, ao_deriv):
+ mo_coeff_mask = mo_coeff[index,:]
+ rho = numint.eval_rho2(opt.mol, ao, mo_coeff_mask, mo_occ, None, xctype)
+ vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype)[1]
+ if xctype == 'LDA':
+ wv = weight * vxc[0]
+ giao = opt.mol.eval_gto('GTOval_ig', coords.get(), comp=3) # (#C(0 1) g) |AO>
+ giao = cupy.array(giao)
+ giao_aux = giao[:,:,index]
+ for idirect in range(3):
+ vtmp = contract('pu,p,vp->uv', giao_aux[idirect], wv, ao)
+ vtmp = cupy.ascontiguousarray(vtmp)
+ add_sparse(vmat[idirect], vtmp, index)
+
+ elif xctype == 'GGA':
+ wv = vxc * weight
+ giao = opt.mol.eval_gto('GTOval_ig', coords.get(), comp=3)
+ giao_nabla = opt.mol.eval_gto('GTOval_ipig', coords.get()).reshape(3, 3, -1, nao)
+ giao = cupy.array(giao)
+ giao_nabla = cupy.array(giao_nabla)
+ giao_aux = giao[:,:,index]
+ giao_nabla_aux = giao_nabla[:,:,:,index]
+ for idirect in range(3):
+ # * write like the gpu4pyscf numint part, but explicitly use the einsum
+ aow = contract('pn,p->pn', giao_aux[idirect], wv[0])
+ aow += contract('xpn,xp->pn', giao_nabla_aux[:, idirect, :, :], wv[1:4])
+ vtmp = contract('pn,mp->nm', aow, ao[0])
+ vtmp = cupy.ascontiguousarray(vtmp)
+ add_sparse(vmat[idirect], vtmp, index)
+ aow = contract('pn,xp->xpn', giao_aux[idirect], wv[1:4])
+ vtmp = contract('xpn,xmp->nm', aow, ao[1:4])
+ vtmp = cupy.ascontiguousarray(vtmp)
+ add_sparse(vmat[idirect], vtmp, index)
+
+ elif xctype == 'NLC':
+ raise NotImplementedError('NLC')
+ elif xctype == 'MGGA':
+ raise NotImplementedError('Meta-GGA')
+ elif xctype == 'HF':
+ pass
+ else:
+ raise NotImplementedError(
+ f'numint.nr_rks for functional {xc_code}')
+
+ ao = None
+
+ # vmat = contract('pi,npq->niq', coeff, vmat)
+ # vmat = contract('qj,niq->nij', coeff, vmat)
+ vmat = take_last2d(vmat, opt.rev_ao_idx)
+
+ if numint.FREE_CUPY_CACHE:
+ dms = None
+ cupy.get_default_memory_pool().free_all_blocks()
+
+ return (vmat - vmat.transpose(0, 2, 1))
+
+
+# TODO: this can be modified into jk, not _vhf
+def get_jk(mol, dm0):
+ # J = Im[(i i|\mu g\nu) + (i gi|\mu \nu)] = -i (i i|\mu g\nu)
+ # K = Im[(\mu gi|i \nu) + (\mu i|i g\nu)]
+ # = [-i (\mu g i|i \nu)] - h.c. (-h.c. for anti-symm because of the factor -i)
+ intor = mol._add_suffix('int2e_ig1')
+ vj, vk = _vhf.direct_mapdm(intor, # (g i,j|k,l)
+ 'a4ij', ('lk->s1ij', 'jk->s1il'),
+ dm0.get(), 3, # xyz, 3 components
+ mol._atm, mol._bas, mol._env)
+ vk = vk - np.swapaxes(vk, -1, -2)
+ return -cupy.array(vj), -cupy.array(vk)
+
+
+def get_vxc(mf, dm0):
+ vxc = nr_rks(mf._numint, mf.mol, mf.grids, mf.xc, mf.make_rdm1())
+ # ! imaginary part
+ vj, vk = get_jk(mf.mol, dm0)
+ if not mf._numint.libxc.is_hybrid_xc(mf.xc):
+ vk = None
+ vxc += vj
+ else:
+ omega, alpha, hyb = mf._numint.rsh_and_hybrid_coeff(
+ mf.xc, spin=mf.mol.spin)
+ vxc += vj - vk*hyb*0.5
+ return vxc
+
+
+def get_h1ao(mf):
+ dm0 = mf.make_rdm1()
+ # ! imaginary part
+ h1ao = -0.5*mf.mol.intor('int1e_giao_irjxp')
+ h1ao += -mf.mol.intor('int1e_igkin')
+ h1ao += -mf.mol.intor('int1e_ignuc')
+ h1ao = cupy.array(h1ao)
+ h1ao += get_vxc(mf, dm0)
+
+ return h1ao
+
+
+def eval_shielding(mf):
+
+ mo_coeff = mf.mo_coeff
+ mo_occ = mf.mo_occ
+ mo_energy = mf.mo_energy
+ mo_coeff = cupy.array(mo_coeff)
+ mo_occ = cupy.array(mo_occ)
+ mo_energy = cupy.array(mo_energy)
+
+ nao = mo_energy.shape[0]
+ idx_occ = mo_occ > 0
+ idx_vir = mo_occ == 0
+ fx = gen_vind(mf, mo_coeff, mo_occ)
+ mocc = mo_coeff[:, idx_occ]
+ mvir = mo_coeff[:, idx_vir]
+ s1ao = -mf.mol.intor('int1e_igovlp')
+ s1ao = cupy.array(s1ao)
+ dm0 = mf.make_rdm1()
+ natom = mf.mol.natm
+
+ h1ao = get_h1ao(mf)
+ tmp = contract('xuv,ua->xav', s1ao, mvir)
+ s1ai = contract('xav,vi->xai', tmp, mocc)
+ tmp = contract('ned,ue->nud', s1ao, dm0)
+ s1dm = contract('nud,dv->nuv', tmp, dm0)*0.5
+ tmp = contract('xpq,pi->xiq', s1ao, mocc)
+ s1jk = -contract('xiq,qj->xij', tmp, mocc)*0.5
+ tmp = contract('nai,ua->nui', s1jk, mocc)
+ s1jkdm1 = contract('nui,vi->nuv', tmp, mocc.conj())*2
+ s1jkdm1 -= s1jkdm1.transpose(0, 2, 1)
+ omega, alpha, hyb = mf._numint.rsh_and_hybrid_coeff(
+ mf.xc, spin=mf.mol.spin)
+ if hasattr(mf,'with_df'):
+ vk2 = cupy.empty((3, nao, nao))
+ for i in range(3):
+ vk2[i] = +mf.get_jk(mf.mol, s1jkdm1[i], hermi=2, with_j=False)[1]*0.5*hyb
+
+ else:
+ vk2 = np.empty((3, nao, nao))
+ for i in range(3):
+ vk2[i] = -jk.get_jk(mf.mol, s1jkdm1[i].get(), 'ijkl,jk->il')*0.5*hyb
+ vk2 = cupy.array(vk2)
+ h1ao += vk2
+ tmp = contract('xuv,ua->xav', h1ao, mvir)
+ Veff_ai = contract('xav,vi->xai', tmp, mocc)
+ Veff_ai -= contract('xai,i->xai', s1ai, mo_energy[idx_occ])
+ Veff_ai = cupy.array(Veff_ai)
+ mo1 = cphf.solve(fx, mo_energy, mo_occ, Veff_ai, max_cycle=20, tol=1e-10)[0]
+
+ shielding_d = cupy.empty((natom, 3, 3))
+ shielding_p = cupy.empty((natom, 3, 3))
+
+ for atm_id in range(natom):
+ mf.mol.set_rinv_origin(mf.mol.atom_coord(atm_id))
+ # ! imaginary part (p* | rinv cross p | )
+ int_h01 = -mf.mol.intor('int1e_prinvxp')
+ # ! imaginary part (g | nabla-rinv cross p |)
+ int_h01_giao = mf.mol.intor('int1e_a01gp').reshape(
+ 3, 3, nao, nao)
+ int_h11 = mf.mol.intor('int1e_giao_a11part').reshape(
+ 3, 3, nao, nao) # ! (-.5 | nabla-rinv | r)
+ int_h11 = cupy.array(int_h11)
+ int_h01 = cupy.array(int_h01)
+ int_h01_giao = cupy.array(int_h01_giao)
+ int_h11_diag = int_h11[0, 0] + int_h11[1, 1] + int_h11[2, 2]
+ int_h11[0, 0] -= int_h11_diag
+ int_h11[1, 1] -= int_h11_diag
+ int_h11[2, 2] -= int_h11_diag
+
+ tmp = contract('ua,xai->xui', mvir, mo1)
+ dm1 = contract('xui,vi->xvu', tmp, mocc)
+ dm1 *= 2
+ shielding_p[atm_id] = contract('yuv,xvu->xy', int_h01, dm1)*2.0
+ shielding_p[atm_id] += contract('yuv,xvu->xy', int_h01, s1dm)
+ shielding_d[atm_id] = contract('xyuv,vu->xy', int_h11, dm0)
+ shielding_d[atm_id] += contract('xyuv,vu->xy', int_h01_giao, dm0)
+ ppm = nist.ALPHA**2 * 1e6
+ return shielding_d*ppm, shielding_p*ppm
diff --git a/gpu4pyscf/properties/tests/test_ir_intensity.py b/gpu4pyscf/properties/tests/test_ir_intensity.py
new file mode 100644
index 00000000..8614b588
--- /dev/null
+++ b/gpu4pyscf/properties/tests/test_ir_intensity.py
@@ -0,0 +1,101 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import unittest
+import numpy as np
+import pyscf
+from pyscf import lib
+from gpu4pyscf.dft import rks, uks
+from gpu4pyscf.properties import ir
+
+lib.num_threads(8)
+
+atom = '''
+O 0.0000000000 -0.0000000000 0.1174000000
+H -0.7570000000 -0.0000000000 -0.4696000000
+H 0.7570000000 0.0000000000 -0.4696000000
+'''
+
+bas='def2tzvpp'
+grids_level = 6
+
+def setUpModule():
+ global mol
+ mol = pyscf.M(atom=atom, basis=bas, max_memory=32000)
+ mol.output = '/dev/null'
+ mol.build()
+ mol.verbose = 1
+
+def tearDownModule():
+ global mol
+ mol.stdout.close()
+ del mol
+
+def run_dft_df_if(xc):
+ mf = rks.RKS(mol, xc=xc).density_fit(auxbasis='def2-universal-jkfit')
+ mf.grids.level = grids_level
+ mf.conv_tol_cpscf = 1e-3
+ e_dft = mf.kernel()
+ h = mf.Hessian()
+ h.auxbasis_response = 2
+ freq, intensity = ir.eval_ir_freq_intensity(mf, h)
+ return e_dft, freq, intensity
+
+
+class KnownValues(unittest.TestCase):
+ '''
+ known values are obtained by Q-Chem
+
+ $rem
+ MEM_STATIC 2000
+ JOBTYPE freq
+ METHOD b3lyp
+ BASIS def2-tzvpp
+ SCF_CONVERGENCE 12
+ XC_GRID 000099000590
+ SYMMETRY FALSE
+ SYM_IGNORE = TRUE
+ ri_j True
+ ri_k True
+ aux_basis RIJK-def2-tzvpp
+ THRESH 14
+ $end
+
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ Mode: 1 2 3
+ Frequency: 1630.86 3850.08 3949.35
+ Force Cnst: 1.6970 9.1261 9.9412
+ Red. Mass: 1.0829 1.0450 1.0818
+ IR Active: YES YES YES
+ IR Intens: 69.334 4.675 47.943
+ Raman Active: YES YES YES
+
+ '''
+
+ def test_rks_b3lyp_df(self):
+ print('-------- RKS density fitting B3LYP -------------')
+ e_tot, freq, intensity = run_dft_df_if('B3LYP')
+ assert np.allclose(e_tot, -76.4666819537)
+ qchem_freq = np.array([1630.86, 3850.08, 3949.35])
+ qchem_intensity = np.array([69.334, 4.675, 47.943])
+ print(freq, intensity)
+ assert np.allclose(freq, qchem_freq, rtol=1e-03)
+ assert np.allclose(intensity, qchem_intensity, rtol=1e-02)
+
+
+
+if __name__ == "__main__":
+ print("Full Tests for ir intensity")
+ unittest.main()
\ No newline at end of file
diff --git a/gpu4pyscf/properties/tests/test_polarizability.py b/gpu4pyscf/properties/tests/test_polarizability.py
new file mode 100644
index 00000000..06c6ed15
--- /dev/null
+++ b/gpu4pyscf/properties/tests/test_polarizability.py
@@ -0,0 +1,128 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import unittest
+import numpy as np
+import pyscf
+from pyscf import lib
+from gpu4pyscf.dft import rks, uks
+from gpu4pyscf.properties import polarizability
+
+lib.num_threads(8)
+
+atom = '''
+O 0.0000000000 -0.0000000000 0.1174000000
+H -0.7570000000 -0.0000000000 -0.4696000000
+H 0.7570000000 0.0000000000 -0.4696000000
+'''
+
+bas='def2tzvpp'
+grids_level = 6
+
+def setUpModule():
+ global mol
+ mol = pyscf.M(atom=atom, basis=bas, max_memory=32000)
+ mol.output = '/dev/null'
+ mol.build()
+ mol.verbose = 1
+
+def tearDownModule():
+ global mol
+ mol.stdout.close()
+ del mol
+
+def run_dft_polarizability(xc):
+ mf = rks.RKS(mol, xc=xc)
+ mf.grids.level = grids_level
+ e_dft = mf.kernel()
+ polar = polarizability.eval_polarizability(mf)
+ return e_dft, polar
+
+def run_dft_df_polarizability(xc):
+ mf = rks.RKS(mol, xc=xc).density_fit(auxbasis='def2-universal-jkfit')
+ mf.grids.level = grids_level
+ e_dft = mf.kernel()
+ polar = polarizability.eval_polarizability(mf)
+ return e_dft, polar
+
+
+class KnownValues(unittest.TestCase):
+ '''
+ known values are obtained by Q-Chem
+
+ $rem
+ MEM_STATIC 2000
+ JOBTYPE POLARIZABILITY
+ METHOD b3lyp
+ RESPONSE_POLAR 0
+ BASIS def2-tzvpp
+ SCF_CONVERGENCE 12
+ XC_GRID 000099000590
+ SYMMETRY FALSE
+ SYM_IGNORE = TRUE
+ $end
+
+ -----------------------------------------------------------------------------
+ Polarizability tensor [a.u.]
+ 8.5899463 -0.0000000 -0.0000000
+ -0.0000000 6.0162267 -0.0000000
+ -0.0000000 -0.0000000 7.5683123
+
+ $rem
+ MEM_STATIC 2000
+ JOBTYPE POLARIZABILITY
+ METHOD b3lyp
+ RESPONSE_POLAR -1
+ BASIS def2-tzvpp
+ SCF_CONVERGENCE 12
+ XC_GRID 000099000590
+ SYMMETRY FALSE
+ SYM_IGNORE = TRUE
+ ri_j True
+ ri_k True
+ aux_basis RIJK-def2-tzvpp
+ THRESH 14
+ $end
+
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ Polarizability tensor [a.u.]
+ 8.5902540 -0.0000000 -0.0000000
+ 0.0000000 6.0167648 0.0000000
+ -0.0000000 -0.0000000 7.5688173
+
+ '''
+ def test_rks_b3lyp(self):
+ print('-------- RKS B3LYP -------------')
+ e_tot, polar = run_dft_polarizability('B3LYP')
+ assert np.allclose(e_tot, -76.4666494276)
+ qchem_polar = np.array([ [ 8.5899463, -0.0000000, -0.0000000],
+ [-0.0000000, 6.0162267, -0.0000000],
+ [-0.0000000, -0.0000000, 7.5683123]])
+ assert np.allclose(polar, qchem_polar)
+
+ def test_rks_b3lyp_df(self):
+ print('-------- RKS density fitting B3LYP -------------')
+ e_tot, polar = run_dft_df_polarizability('B3LYP')
+ assert np.allclose(e_tot, -76.4666819553)
+ qchem_polar = np.array([ [ 8.5902540, -0.0000000, -0.0000000],
+ [ 0.0000000, 6.0167648, -0.0000000],
+ [ -0.0000000, -0.0000000, 7.5688173]])
+ assert np.allclose(polar, qchem_polar)
+
+
+
+if __name__ == "__main__":
+ print("Full Tests for polarizabillity")
+ unittest.main()
\ No newline at end of file
diff --git a/gpu4pyscf/properties/tests/test_shielding.py b/gpu4pyscf/properties/tests/test_shielding.py
new file mode 100644
index 00000000..04352ad1
--- /dev/null
+++ b/gpu4pyscf/properties/tests/test_shielding.py
@@ -0,0 +1,140 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+import unittest
+import numpy as np
+import pyscf
+from pyscf import lib
+from gpu4pyscf.dft import rks, uks
+from gpu4pyscf.properties import shielding
+
+lib.num_threads(8)
+
+atom = '''
+O 0.0000000000 -0.0000000000 0.1174000000
+H -0.7570000000 -0.0000000000 -0.4696000000
+H 0.7570000000 0.0000000000 -0.4696000000
+'''
+
+bas='def2tzvpp'
+grids_level = 6
+
+def setUpModule():
+ global mol
+ mol = pyscf.M(atom=atom, basis=bas, max_memory=32000)
+ mol.output = '/dev/null'
+ mol.build()
+ mol.verbose = 1
+
+def tearDownModule():
+ global mol
+ mol.stdout.close()
+ del mol
+
+def run_dft_nmr_shielding(xc):
+ mf = rks.RKS(mol, xc=xc)
+ mf.grids.level = grids_level
+ e_dft = mf.kernel()
+ tensor = shielding.eval_shielding(mf)
+ return e_dft, tensor
+
+def run_dft_df_nmr_shielding(xc):
+ mf = rks.RKS(mol, xc=xc).density_fit(auxbasis='def2-universal-jkfit')
+ mf.grids.level = grids_level
+ e_dft = mf.kernel()
+ tensor = shielding.eval_shielding(mf)
+ return e_dft, tensor
+
+
+class KnownValues(unittest.TestCase):
+ '''
+ known values are obtained by Q-Chem
+
+ $rem
+ METHOD b3lyp
+ BASIS def2-tzvpp
+ SCF_CONVERGENCE 12
+ XC_GRID 000099000590
+ SYMMETRY FALSE
+ SYM_IGNORE TRUE
+ MOPROP 1
+ MOPROP_PERTNUM 0 ! do all perturbations at once
+ MOPROP_CONV_1ST 10 ! sets the CPSCF convergence threshold
+ MOPROP_DIIS_DIM_SS 4 ! no. of DIIS subspace vectors
+ MOPROP_MAXITER_1ST 100 ! max iterations
+ MOPROP_DIIS 5 ! turns on DIIS (=0 to turn off)
+ MOPROP_DIIS_THRESH 1
+ MOPROP_DIIS_SAVE 0
+ $end
+
+ ATOM ISOTROPIC ANISOTROPIC REL. SHIFTS
+------------------------------------------------------------------------
+ Atom O 1 332.07586444 45.66898820
+ Atom H 2 31.39150070 19.36432602
+ Atom H 3 31.39060707 19.36602790
+
+ $rem
+ METHOD b3lyp
+ BASIS def2-tzvpp
+ SCF_CONVERGENCE 10
+ XC_GRID 000099000590
+ ri_j True
+ ri_k True
+ aux_basis RIJK-def2-tzvpp
+ THRESH 14
+ SYMMETRY FALSE
+ SYM_IGNORE TRUE
+ MOPROP 1
+ MOPROP_PERTNUM 0 ! do all perturbations at once
+ MOPROP_CONV_1ST 10 ! sets the CPSCF convergence threshold
+ MOPROP_DIIS_DIM_SS 4 ! no. of DIIS subspace vectors
+ MOPROP_MAXITER_1ST 100 ! max iterations
+ MOPROP_DIIS 5 ! turns on DIIS (=0 to turn off)
+ MOPROP_DIIS_THRESH 1
+ MOPROP_DIIS_SAVE 0
+ $end
+
+
+ ATOM ISOTROPIC ANISOTROPIC REL. SHIFTS
+------------------------------------------------------------------------
+ Atom O 1 332.07961083 45.66777197
+ Atom H 2 31.39250594 19.36477246
+ Atom H 3 31.39160966 19.36647972
+
+ '''
+ def test_rks_b3lyp(self):
+ print('-------- RKS B3LYP -------------')
+ e_tot, tensor = run_dft_nmr_shielding('B3LYP')
+ nmr_total = tensor[0].get()+tensor[1].get()
+ isotropic_pyscf = np.array([nmr_total[0].trace()/3, nmr_total[1].trace()/3, nmr_total[2].trace()/3])
+ isotropic_qchem = np.array([332.07586444, 31.39150070, 31.39060707])
+
+ assert np.allclose(e_tot, -76.4666494276)
+ assert np.allclose(isotropic_pyscf, isotropic_qchem, rtol=1.0E-4)
+
+ def test_rks_b3lyp_df(self):
+ print('-------- RKS density fitting B3LYP -------------')
+ e_tot, tensor = run_dft_df_nmr_shielding('B3LYP')
+ nmr_total = tensor[0].get()+tensor[1].get()
+ isotropic_pyscf = np.array([nmr_total[0].trace()/3, nmr_total[1].trace()/3, nmr_total[2].trace()/3])
+ isotropic_qchem = np.array([332.07961083, 31.39250594, 31.39160966])
+ assert np.allclose(e_tot, -76.4666819553)
+ assert np.allclose(isotropic_pyscf, isotropic_qchem, rtol=1.0E-4)
+
+
+
+if __name__ == "__main__":
+ print("Full Tests for nmr shielding constants")
+ unittest.main()
\ No newline at end of file