From f2e139185eea0003a0b851e44fc8c4fc3d09a83f Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 14 Jan 2020 22:49:22 -0500 Subject: [PATCH 01/11] Tighten ground-state SCF convergence thresholds --- Response-Theory/Self-Consistent-Field/CPHF.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/CPHF.py b/Response-Theory/Self-Consistent-Field/CPHF.py index 2769633d..f23f10b0 100644 --- a/Response-Theory/Self-Consistent-Field/CPHF.py +++ b/Response-Theory/Self-Consistent-Field/CPHF.py @@ -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 From 9ee0eadf69ef42a7c4c9e735a78683e42c4b8146 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 14 Jan 2020 22:57:16 -0500 Subject: [PATCH 02/11] Make sure the RHS and response vectors from a previous run are cleared out --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 15e5ca3a..a059386e 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -84,6 +84,8 @@ def run(self, method='direct', omega=None): self.form_polarizability() def solve_static_direct(self): + self.x = None + self.rhsvecs = None # Run a quick check to make sure everything will fit into memory I_Size = (self.nbf ** 4) * 8.e-9 oNNN_Size = (self.nocc * self.nbf ** 3) * 8.e-9 @@ -137,6 +139,8 @@ def solve_static_direct(self): self.rhsvecs.append(rhsvec) def solve_dynamic_direct(self, omega=0.0): + self.x = None + self.rhsvecs = None # Adapted completely from TDHF.py eps_v = self.epsilon[self.nocc:] @@ -197,6 +201,8 @@ def solve_dynamic_direct(self, omega=0.0): self.x.append(xcomp) def solve_static_iterative(self, maxiter=20, conv=1.e-9, use_diis=True): + self.x = None + self.rhsvecs = None # Init JK object jk = psi4.core.JK.build(self.scf_wfn.basisset()) @@ -274,6 +280,8 @@ def solve_static_iterative(self, maxiter=20, conv=1.e-9, use_diis=True): (CPHF_ITER, avg_RMS, max_RMS)) def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-9, use_diis=True): + self.x = None + self.rhsvecs = None # Init JK object jk = psi4.core.JK.build(self.scf_wfn.basisset()) From bcd7f0dc18784ba98480651794c8bda03354926e Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 14 Jan 2020 23:09:48 -0500 Subject: [PATCH 03/11] Initial rework of dynamic iterative CPHF --- .../Self-Consistent-Field/helper_CPHF.py | 118 ++++++++++-------- 1 file changed, 65 insertions(+), 53 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index a059386e..2a7d5e9e 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -279,50 +279,55 @@ 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): self.x = None self.rhsvecs = None + # Convert Co and Cv to numpy arrays + Co = np.asarray(self.Co) + Cv = np.asarray(self.Cv) + + # 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] + x = [] + x_old = [] + diis = [] + for rhsmat in rhsmats: + U = np.zeros((norb, norb)) + for i in range(self.nocc): + for a in range(self.nocc, norb): + U[i, a] = rhsmat[i, a] / (self.epsilon[i] - self.epsilon[a] - omega) + U[a, i] = rhsmat[a, i] / (self.epsilon[a] - self.epsilon[i] - omega) + x.append(U) + x_old.append(np.zeros_like(U)) + diis.append(DIIS_helper()) + # 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): + npCx_l = [] + npCx_r = [] + for _ in range(len(rhsmats)): + mCx_l = psi4.core.Matrix(self.nbf, self.nocc) + npCx_l.append(np.asarray(mCx_l)) 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()) - - # Convert Co and Cv to numpy arrays - Co = np.asarray(self.Co) - Cv = np.asarray(self.Cv) + jk.C_right_add(mCx_l) + for _ in range(len(rhsmats)): + mCx_r = psi4.core.Matrix(self.nbf, self.nocc) + npCx_r.append(np.asarray(mCx_r)) + jk.C_left_add(mCx_r) + 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) + npCx_l[xyz][:] = C.dot(x[xyz].T)[:, :self.nocc] + npCx_r[xyz][:] = (-C).dot(x[xyz])[:, :self.nocc] # Perform generalized J/K build jk.compute() @@ -336,30 +341,27 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-9, use_diis=Tr 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 - + U = x[xyz].copy() + GPxMO = (C.T).dot(2 * J_l - K_l + 2 * J_r - K_r).dot(C) + AU = np.zeros_like(GPxMO) + for p in range(norb): + for q in range(norb): + AU[p, q] = x[xyz][p, q] * (self.epsilon[p] - self.epsilon[q] - omega) - GPxMO[p, q] + for i in range(self.nocc): + for a in range(self.nocc, norb): + U[i, a] += (rhsmats[xyz][i, a] - AU[i, a]) / (self.epsilon[i] - self.epsilon[a] - omega) + U[a, i] += (rhsmats[xyz][a, i] - AU[a, i]) / (self.epsilon[a] - self.epsilon[i] - omega) # 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) @@ -367,11 +369,21 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-9, use_diis=Tr if max_RMS < conv: print('CPHF converged in %d iterations and %.2f seconds.' % (CPHF_ITER, time.time() - t)) self.rhsvecs = [] - 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.x = [] + # do one last perturbed density build in the AO basis + for i, rhsmat in enumerate(rhsmats): + l = C[:, :self.nocc] @ C.dot(x[i].T)[:, :self.nocc].T - C.dot(x[i])[:, :self.nocc] @ C[:, :self.nocc].T + r = self.tmp_dipoles[i].np + self.rhsvecs.append(-2 * l.reshape(-1)) + self.x.append(r.reshape(-1)) + # res = l.reshape(-1).dot(r.reshape(-1)) * -2 + # 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)))) break print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' % From 7e4ed94cf6d79e0b31d0172b27feb40d50709546 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sat, 18 Jan 2020 18:40:00 -0500 Subject: [PATCH 04/11] Fix storage of RHS and response vectors for dynamic iterative solver --- .../Self-Consistent-Field/helper_CPHF.py | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 2a7d5e9e..6118935b 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -370,24 +370,26 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T print('CPHF converged in %d iterations and %.2f seconds.' % (CPHF_ITER, time.time() - t)) self.rhsvecs = [] self.x = [] - # do one last perturbed density build in the AO basis - for i, rhsmat in enumerate(rhsmats): - l = C[:, :self.nocc] @ C.dot(x[i].T)[:, :self.nocc].T - C.dot(x[i])[:, :self.nocc] @ C[:, :self.nocc].T - r = self.tmp_dipoles[i].np - self.rhsvecs.append(-2 * l.reshape(-1)) - self.x.append(r.reshape(-1)) - # res = l.reshape(-1).dot(r.reshape(-1)) * -2 - # 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)))) + for numx in range(3): + 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' % - (CPHF_ITER, avg_RMS, max_RMS)) + print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' % (CPHF_ITER, avg_RMS, max_RMS)) def form_polarizability(self): self.polar = np.empty((3, 3)) From 7a092fc8bf9dbf69fbb12da084cfae92bf126985 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sat, 18 Jan 2020 19:21:25 -0500 Subject: [PATCH 05/11] Replace loop update of AU with block assignment --- .../Self-Consistent-Field/helper_CPHF.py | 21 +++++++------------ 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 6118935b..6c67277d 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -283,10 +283,6 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T self.x = None self.rhsvecs = None - # Convert Co and Cv to numpy arrays - Co = np.asarray(self.Co) - Cv = np.asarray(self.Cv) - # Build initial guess, previous vectors, and DIIS objects norb = self.scf_wfn.nmo() C = np.asarray(self.C) @@ -294,12 +290,12 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T 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)) - for i in range(self.nocc): - for a in range(self.nocc, norb): - U[i, a] = rhsmat[i, a] / (self.epsilon[i] - self.epsilon[a] - omega) - U[a, i] = rhsmat[a, i] / (self.epsilon[a] - self.epsilon[i] - omega) + 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()) @@ -342,11 +338,7 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T # Bulid new guess U = x[xyz].copy() - GPxMO = (C.T).dot(2 * J_l - K_l + 2 * J_r - K_r).dot(C) - AU = np.zeros_like(GPxMO) - for p in range(norb): - for q in range(norb): - AU[p, q] = x[xyz][p, q] * (self.epsilon[p] - self.epsilon[q] - omega) - GPxMO[p, q] + AU = x[xyz] * all_denom - (C.T).dot(2 * J_l - K_l + 2 * J_r - K_r).dot(C) for i in range(self.nocc): for a in range(self.nocc, norb): U[i, a] += (rhsmats[xyz][i, a] - AU[i, a]) / (self.epsilon[i] - self.epsilon[a] - omega) @@ -389,7 +381,8 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T ) break - print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' % (CPHF_ITER, avg_RMS, max_RMS)) + print('CPHF Iteration %3d: Average RMS = %3.8f Maximum RMS = %3.8f' % + (CPHF_ITER, avg_RMS, max_RMS)) def form_polarizability(self): self.polar = np.empty((3, 3)) From b04fdc5a4dd84a55e4ebc97c861101779b0fa779 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sat, 18 Jan 2020 22:36:49 -0500 Subject: [PATCH 06/11] Fully vectorize dynamic iterative solver --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 6c67277d..7eed6674 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -338,11 +338,9 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T # Bulid new guess U = x[xyz].copy() - AU = x[xyz] * all_denom - (C.T).dot(2 * J_l - K_l + 2 * J_r - K_r).dot(C) - for i in range(self.nocc): - for a in range(self.nocc, norb): - U[i, a] += (rhsmats[xyz][i, a] - AU[i, a]) / (self.epsilon[i] - self.epsilon[a] - omega) - U[a, i] += (rhsmats[xyz][a, i] - AU[a, i]) / (self.epsilon[a] - self.epsilon[i] - omega) + upd = (rhsmats[xyz] - (x[xyz] * all_denom - (C.T).dot(2 * J_l - K_l + 2 * J_r - K_r).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[xyz].add(U, U - x_old[xyz]) From 1d4f23bfc997844ea75ac684e5c965fc9fb4516c Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 19 Jan 2020 13:39:44 -0500 Subject: [PATCH 07/11] Remove explicit counts for number of perturbation matrices --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 7eed6674..780e89e0 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -321,7 +321,7 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T t = time.time() for CPHF_ITER in range(1, maxiter + 1): - for xyz in range(3): + for xyz in range(len(rhsmats)): npCx_l[xyz][:] = C.dot(x[xyz].T)[:, :self.nocc] npCx_r[xyz][:] = (-C).dot(x[xyz])[:, :self.nocc] @@ -329,12 +329,12 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T 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]) + J_r = np.asarray(jk.J()[xyz + len(rhsmats)]) + K_r = np.asarray(jk.K()[xyz + len(rhsmats)]) # Bulid new guess U = x[xyz].copy() From a98461a04e253143114cd0531c5274984979e78e Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 19 Jan 2020 16:43:54 -0500 Subject: [PATCH 08/11] Update code documentation --- .../Self-Consistent-Field/helper_CPHF.py | 44 +++++++++++-------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 780e89e0..a92dcdff 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -304,17 +304,23 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T jk = psi4.core.JK.build(self.scf_wfn.basisset()) jk.initialize() - npCx_l = [] - npCx_r = [] + # 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_l = psi4.core.Matrix(self.nbf, self.nocc) - npCx_l.append(np.asarray(mCx_l)) + 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_l) + jk.C_right_add(mCx_1) for _ in range(len(rhsmats)): - mCx_r = psi4.core.Matrix(self.nbf, self.nocc) - npCx_r.append(np.asarray(mCx_r)) - jk.C_left_add(mCx_r) + 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:') @@ -322,23 +328,25 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T for CPHF_ITER in range(1, maxiter + 1): for xyz in range(len(rhsmats)): - npCx_l[xyz][:] = C.dot(x[xyz].T)[:, :self.nocc] - npCx_r[xyz][:] = (-C).dot(x[xyz])[:, :self.nocc] + 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(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 + len(rhsmats)]) - K_r = np.asarray(jk.K()[xyz + len(rhsmats)]) - - # Bulid new guess + 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_l - K_l + 2 * J_r - K_r).dot(C))) / all_denom + 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 From 9a4309e9a75d706839ce7bb9945a2e6d78845ec0 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 19 Jan 2020 16:46:06 -0500 Subject: [PATCH 09/11] Move JK object initialization up --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index a92dcdff..a1d7c25f 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -283,6 +283,10 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T self.x = None self.rhsvecs = None + # Init JK object + jk = psi4.core.JK.build(self.scf_wfn.basisset()) + jk.initialize() + # Build initial guess, previous vectors, and DIIS objects norb = self.scf_wfn.nmo() C = np.asarray(self.C) @@ -300,10 +304,6 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T x_old.append(np.zeros_like(U)) diis.append(DIIS_helper()) - # 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, # which will contain MO coefficients contracted with the response # matrix. 1 and 2 refer to the first and second terms in the perturbed From 18bcb332bb582981163cfca133bb4762b7cb4ef6 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 19 Jan 2020 17:37:59 -0500 Subject: [PATCH 10/11] Remove duplicate resets of x abd rhsvecs --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index a1d7c25f..10c851e5 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -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': @@ -84,8 +81,6 @@ def run(self, method='direct', omega=None): self.form_polarizability() def solve_static_direct(self): - self.x = None - self.rhsvecs = None # Run a quick check to make sure everything will fit into memory I_Size = (self.nbf ** 4) * 8.e-9 oNNN_Size = (self.nocc * self.nbf ** 3) * 8.e-9 @@ -139,8 +134,6 @@ def solve_static_direct(self): self.rhsvecs.append(rhsvec) def solve_dynamic_direct(self, omega=0.0): - self.x = None - self.rhsvecs = None # Adapted completely from TDHF.py eps_v = self.epsilon[self.nocc:] @@ -201,8 +194,6 @@ def solve_dynamic_direct(self, omega=0.0): self.x.append(xcomp) def solve_static_iterative(self, maxiter=20, conv=1.e-9, use_diis=True): - self.x = None - self.rhsvecs = None # Init JK object jk = psi4.core.JK.build(self.scf_wfn.basisset()) @@ -280,8 +271,6 @@ def solve_static_iterative(self, maxiter=20, conv=1.e-9, use_diis=True): (CPHF_ITER, avg_RMS, max_RMS)) def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=True): - self.x = None - self.rhsvecs = None # Init JK object jk = psi4.core.JK.build(self.scf_wfn.basisset()) From bb825c13dfa201cf4add691dda759fc83ca7e068 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 21 Jan 2020 22:24:56 -0500 Subject: [PATCH 11/11] Switch from matmul operator to np.dot --- Response-Theory/Self-Consistent-Field/helper_CPHF.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Response-Theory/Self-Consistent-Field/helper_CPHF.py b/Response-Theory/Self-Consistent-Field/helper_CPHF.py index 10c851e5..8886ebaf 100644 --- a/Response-Theory/Self-Consistent-Field/helper_CPHF.py +++ b/Response-Theory/Self-Consistent-Field/helper_CPHF.py @@ -279,7 +279,7 @@ def solve_dynamic_iterative(self, omega=0.0, maxiter=20, conv=1.e-10, use_diis=T # 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] + rhsmats = [(C.T).dot(np.asarray(dipmat)).dot(C) for dipmat in self.tmp_dipoles] x = [] x_old = [] diis = []