-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathpolarizability.py
132 lines (107 loc) · 4.5 KB
/
polarizability.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
# static polarizability module
from pyscf import gto, dft, grad, scf
import numpy as np
from pyscf.scf.hf import RHF
from pyscf.dft.rks import RKS
from pyscf.lib import param
# TODO: you can write an iterface code such that you don't
# need to decorate every possible SCF class
# See pyscf/qmmm/itrf.py for an example
# Hartree-Fork SCF Method
def hcore_caused_by_elec_field(mol, e_field=np.array([0,0,0])):
h1 = mol.intor('int1e_r')
return np.einsum('xij,x->ij', h1, e_field)
class ModifiedRHF(RHF):
def __init__(self, mol, e_field=np.zeros(3)):
RHF.__init__(self, mol)
self.e_field = e_field
def get_hcore(self, mol=None):
if mol is None: mol = self.mol
return RHF.get_hcore(self, mol) + hcore_caused_by_elec_field(mol, self.e_field)
def get_dipole_with_field_RHF(mol, e_field=np.zeros(3)):
mf = ModifiedRHF(mol, e_field=e_field)
mf.verbose = 0
mf.stdout = '/dev/null'
mf.conv_tol_grad = 1e-8
mf.max_cycle = 100
mf.kernel()
dm = mf.make_rdm1()
return mf.dip_moment(mol, dm, unit='AU', verbose=0)
def get_dipole_at_grid_RHF(mol, direction, h):
# preparation for central difference
if direction == 'x':
elec_grid = np.array([[i * h, 0, 0] for i in range(-2,3)])
elif direction == 'y':
elec_grid = np.array([[0, i * h, 0] for i in range(-2,3)])
elif direction == 'z':
elec_grid = np.array([[0, 0, i * h] for i in range(-2,3)])
dipole_at_grid = []
for e_field in elec_grid:
dipole_at_grid.append(get_dipole_with_field_RHF(mol, e_field))
return np.array(dipole_at_grid)
def get_polarizability_RHF(mol):
# unit: Angstrom**3
# 4th order central difference
h = 0.001
central_difference_coeff = np.array([1/12, -2/3, 0, 2/3, -1/12])/h
# calculate polarizability
polarizability_caused_by_ex = np.dot(get_dipole_at_grid_RHF(mol, 'x', h).T, central_difference_coeff)
polarizability_caused_by_ey = np.dot(get_dipole_at_grid_RHF(mol, 'y', h).T, central_difference_coeff)
polarizability_caused_by_ez = np.dot(get_dipole_at_grid_RHF(mol, 'z', h).T, central_difference_coeff)
# keep the column vector format
polarizability = np.array([polarizability_caused_by_ex, polarizability_caused_by_ey, polarizability_caused_by_ez]).T
# convert from a.u. to Angstrom**3
polarizability *= param.BOHR**3
return polarizability
def get_beta_RHF(mol):
alpha = get_polarizability_RHF(mol)
a = np.trace(alpha)/3
beta = alpha - a*np.identity(3)
return beta
# DFT Method
class ModifiedRKS(RKS):
def __init__(self, mol, xc='LDA,VWN', e_field=np.zeros(3)):
RKS.__init__(self, mol, xc=xc)
self.e_field = e_field
def get_hcore(self, mol=None):
if mol is None: mol = self.mol
return RKS.get_hcore(self, mol) + hcore_caused_by_elec_field(mol, self.e_field)
def get_dipole_with_field_DFT(mol, e_field=np.zeros(3)):
mf = ModifiedRKS(mol, e_field=e_field, xc='b3lyp')
mf.verbose = 0
mf.stdout = '/dev/null'
mf.conv_tol_grad = 1e-7
mf.max_cycle = 100
mf.kernel()
dm = mf.make_rdm1()
return mf.dip_moment(mol, dm, unit='AU', verbose=0)
def get_dipole_at_grid_DFT(mol, direction, h):
if direction == 'x':
elec_grid = np.array([[i * h, 0, 0] for i in range(-2,3)])
elif direction == 'y':
elec_grid = np.array([[0, i * h, 0] for i in range(-2,3)])
elif direction == 'z':
elec_grid = np.array([[0, 0, i * h] for i in range(-2,3)])
dipole_at_grid = []
for e_field in elec_grid:
dipole_at_grid.append(get_dipole_with_field_DFT(mol, e_field))
return np.array(dipole_at_grid)
def get_polarizability_DFT(mol):
# unit: Angstrom**3
# 4th order central difference
h = 0.001
central_difference_coeff = np.array([1/12, -2/3, 0, 2/3, -1/12])/h
# calculate polarizability
polarizability_caused_by_ex = np.dot(get_dipole_at_grid_DFT(mol, 'x', h).T, central_difference_coeff)
polarizability_caused_by_ey = np.dot(get_dipole_at_grid_DFT(mol, 'y', h).T, central_difference_coeff)
polarizability_caused_by_ez = np.dot(get_dipole_at_grid_DFT(mol, 'z', h).T, central_difference_coeff)
# keep the column vector format
polarizability = np.array([polarizability_caused_by_ex, polarizability_caused_by_ey, polarizability_caused_by_ez]).T
# convert from a.u. to Angstrom**3
polarizability *= param.BOHR**3
return polarizability
def get_beta_DFT(mol):
alpha = get_polarizability_DFT(mol)
a = np.trace(alpha)/3
beta = alpha - a*np.identity(3)
return beta