Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix frequency-dependent iterative CPHF #91

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Open

Fix frequency-dependent iterative CPHF #91

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
f2e1391
Tighten ground-state SCF convergence thresholds
berquist Jan 15, 2020
9ee0ead
Make sure the RHS and response vectors from a previous run are cleare…
berquist Jan 15, 2020
bcd7f0d
Initial rework of dynamic iterative CPHF
berquist Jan 15, 2020
7e4ed94
Fix storage of RHS and response vectors for dynamic iterative solver
berquist Jan 18, 2020
7a092fc
Replace loop update of AU with block assignment
berquist Jan 19, 2020
b04fdc5
Fully vectorize dynamic iterative solver
berquist Jan 19, 2020
1d4f23b
Remove explicit counts for number of perturbation matrices
berquist Jan 19, 2020
a98461a
Update code documentation
berquist Jan 19, 2020
9a4309e
Move JK object initialization up
berquist Jan 19, 2020
18bcb33
Remove duplicate resets of x abd rhsvecs
berquist Jan 19, 2020
bb825c1
Switch from matmul operator to np.dot
berquist Jan 22, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Response-Theory/Self-Consistent-Field/CPHF.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,8 +42,8 @@
psi4.set_options({"basis": "aug-cc-pVDZ",
"scf_type": "direct",
"df_scf_guess": False,
"e_convergence": 1e-9,
"d_convergence": 1e-9,
"e_convergence": 1e-11,
"d_convergence": 1e-11,
"cphf_tasks": ['polarizability']})

# Set defaults
Expand Down
138 changes: 74 additions & 64 deletions Response-Theory/Self-Consistent-Field/helper_CPHF.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,9 +64,6 @@ def __init__(self, mol, numpy_memory=2):
Fia *= -2
self.dipoles_xyz.append(Fia)

self.x = None
self.rhsvecs = None

def run(self, method='direct', omega=None):
self.method = method
if self.method == 'direct':
Expand Down Expand Up @@ -273,97 +270,110 @@ def solve_static_iterative(self, maxiter=20, conv=1.e-9, use_diis=True):
print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' %
(CPHF_ITER, avg_RMS, max_RMS))

def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-9, use_diis=True):
def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=True):

# Init JK object
jk = psi4.core.JK.build(self.scf_wfn.basisset())
jk.initialize()

# Add blank matrices to the JK object and NumPy hooks to
# C_right; there are 6 sets of matrices to account for X and Y
# vectors separately.
npC_right = []
for xyz in range(6):
jk.C_left_add(self.Co)
mC = psi4.core.Matrix(self.nbf, self.nocc)
npC_right.append(np.asarray(mC))
jk.C_right_add(mC)

# Build initial guess, previous vectors, diis object, and C_left updates
x_l, x_r = [], []
x_l_old, x_r_old = [], []
diis_l, diis_r = [], []
ia_denom_l = self.epsilon[self.nocc:] - self.epsilon[:self.nocc].reshape(-1, 1) - omega
ia_denom_r = self.epsilon[self.nocc:] - self.epsilon[:self.nocc].reshape(-1, 1) + omega
for xyz in range(3):
x_l.append(self.dipoles_xyz[xyz] / ia_denom_l)
x_r.append(self.dipoles_xyz[xyz] / ia_denom_r)
x_l_old.append(np.zeros(ia_denom_l.shape))
x_r_old.append(np.zeros(ia_denom_r.shape))
diis_l.append(DIIS_helper())
diis_r.append(DIIS_helper())
# Build initial guess, previous vectors, and DIIS objects
norb = self.scf_wfn.nmo()
C = np.asarray(self.C)
rhsmats = [C.T @ np.asarray(dipmat) @ C for dipmat in self.tmp_dipoles]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you use np.dot here. I do like the @ syntax better, just not sure if people are ready for it yet.

x = []
x_old = []
diis = []
ia_denom = self.epsilon[self.nocc:] - self.epsilon[:self.nocc].reshape(-1, 1) - omega
all_denom = self.epsilon.reshape(-1, 1) - self.epsilon - omega
for rhsmat in rhsmats:
U = np.zeros((norb, norb))
U[:self.nocc, self.nocc:] = rhsmat[:self.nocc, self.nocc:] / ia_denom
U[self.nocc:, :self.nocc] = rhsmat[self.nocc:, :self.nocc] / -ia_denom.T
x.append(U)
x_old.append(np.zeros_like(U))
diis.append(DIIS_helper())

# Convert Co and Cv to numpy arrays
Co = np.asarray(self.Co)
Cv = np.asarray(self.Cv)
# Add blank matrices to the jk object and NumPy hooks to C_right,
# which will contain MO coefficients contracted with the response
# matrix. 1 and 2 refer to the first and second terms in the perturbed
# density build:
# P^{x} = C @ U^{x}.T @ C.T - C @ U^{x} @ C.T
# = C @ C^{x}_1.T + C^{x}_2 @ C.T
npCx_1 = []
npCx_2 = []
for _ in range(len(rhsmats)):
mCx_1 = psi4.core.Matrix(self.nbf, self.nocc)
npCx_1.append(np.asarray(mCx_1))
jk.C_left_add(self.Co)
jk.C_right_add(mCx_1)
for _ in range(len(rhsmats)):
mCx_2 = psi4.core.Matrix(self.nbf, self.nocc)
npCx_2.append(np.asarray(mCx_2))
jk.C_left_add(mCx_2)
jk.C_right_add(self.Co)

print('\nStarting CPHF iterations:')
t = time.time()
for CPHF_ITER in range(1, maxiter + 1):

# Update jk's C_right; ordering is Xx, Xy, Xz, Yx, Yy, Yz
for xyz in range(3):
npC_right[xyz][:] = Cv.dot(x_l[xyz].T)
npC_right[xyz + 3][:] = Cv.dot(x_r[xyz].T)
for xyz in range(len(rhsmats)):
npCx_1[xyz][:] = C.dot(x[xyz].T)[:, :self.nocc]
npCx_2[xyz][:] = (-C).dot(x[xyz])[:, :self.nocc]

# Perform generalized J/K build
jk.compute()

# Update amplitudes
for xyz in range(3):
for xyz in range(len(rhsmats)):
# Build J and K objects
J_l = np.asarray(jk.J()[xyz])
K_l = np.asarray(jk.K()[xyz])
J_r = np.asarray(jk.J()[xyz + 3])
K_r = np.asarray(jk.K()[xyz + 3])

# Bulid new guess
X_l = self.dipoles_xyz[xyz].copy()
X_r = self.dipoles_xyz[xyz].copy()
X_l -= (Co.T).dot(2 * J_l - K_l).dot(Cv)
X_r -= (Co.T).dot(2 * J_r - K_r).dot(Cv)
X_l /= ia_denom_l
X_r /= ia_denom_r

J_1 = np.asarray(jk.J()[xyz])
K_1 = np.asarray(jk.K()[xyz])
J_2 = np.asarray(jk.J()[xyz + len(rhsmats)])
K_2 = np.asarray(jk.K()[xyz + len(rhsmats)])

# Build new guess: work in the full [norb, norb] space, then
# select only those parameters corresponding to
# occ->virt/virt->occ rotation, leaving the occ-occ and
# virt-virt parameters zero.
U = x[xyz].copy()
upd = (rhsmats[xyz] - (x[xyz] * all_denom - (C.T).dot(2 * J_1 - K_1 + 2 * J_2 - K_2).dot(C))) / all_denom
U[:self.nocc, self.nocc:] += upd[:self.nocc, self.nocc:]
U[self.nocc:, :self.nocc] += upd[self.nocc:, :self.nocc]
# DIIS for good measure
if use_diis:
diis_l[xyz].add(X_l, X_l - x_l_old[xyz])
X_l = diis_l[xyz].extrapolate()
diis_r[xyz].add(X_r, X_r - x_r_old[xyz])
X_r = diis_r[xyz].extrapolate()
x_l[xyz] = X_l.copy()
x_r[xyz] = X_r.copy()
diis[xyz].add(U, U - x_old[xyz])
U = diis[xyz].extrapolate()
x[xyz] = U.copy()

# Check for convergence
rms = []
for xyz in range(3):
rms_l = np.max((x_l[xyz] - x_l_old[xyz]) ** 2)
rms_r = np.max((x_r[xyz] - x_r_old[xyz]) ** 2)
rms.append(max(rms_l, rms_r))
x_l_old[xyz] = x_l[xyz]
x_r_old[xyz] = x_r[xyz]
rms.append(np.max((x[xyz] - x_old[xyz]) ** 2))
x_old[xyz] = x[xyz]

avg_RMS = sum(rms) / 3
max_RMS = max(rms)

if max_RMS < conv:
print('CPHF converged in %d iterations and %.2f seconds.' % (CPHF_ITER, time.time() - t))
self.rhsvecs = []
self.x = []
for numx in range(3):
rhsvec = self.dipoles_xyz[numx].reshape(-1)
self.rhsvecs.append(np.concatenate((rhsvec, -rhsvec)))
self.x.append(np.concatenate((x_l[numx].reshape(-1),
x_r[numx].reshape(-1))))
self.rhsvecs.append(
np.concatenate(
(
self.dipoles_xyz[numx].reshape(-1),
-self.dipoles_xyz[numx].T.reshape(-1),
)
)
)
self.x.append(
np.concatenate(
(
x[numx][: self.nocc, self.nocc :].reshape(-1),
x[numx][self.nocc :, : self.nocc].reshape(-1),
)
)
)
break

print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' %
Expand Down