From 5eadf397a0d8a8ece79a0872886bf145dd77713b Mon Sep 17 00:00:00 2001 From: Daniel Claudino Date: Fri, 9 Aug 2024 13:21:26 +0000 Subject: [PATCH] Added alternate spins ordering to chemistry Observable plugins Signed-off-by: Daniel Claudino --- python/plugins/observables/psi4_observable.py | 34 +++++++- .../plugins/observables/pyscf_observable.py | 82 +++++++++++++------ 2 files changed, 88 insertions(+), 28 deletions(-) diff --git a/python/plugins/observables/psi4_observable.py b/python/plugins/observables/psi4_observable.py index 669ca97b2..aa43c7cb2 100644 --- a/python/plugins/observables/psi4_observable.py +++ b/python/plugins/observables/psi4_observable.py @@ -13,6 +13,15 @@ def __init__(self): xacc.Observable.__init__(self) self.observable = None self.asPauli = None + self.index_map = None + + def _build_index_map(self, n_bits): + if self.index_map is None: + self.index_map = {} + for i in range(n_bits // 2): + self.index_map[i] = 2 * i + self.index_map[i + n_bits // 2] = 2 * i + 1 + return def toString(self): return self.observable.toString() @@ -137,6 +146,12 @@ def spin_block_tei(I): pos_or_neg = lambda x : ' + ' if x > 0. else ' - ' + alternate_spins = False + if 'alternate-spins' in inputParams: + alternate_spins = inputParams['alternate-spins'] + if alternate_spins: + self._build_index_map(n_active) + # --- 1-body frozen-core: Hamiltonian_fc_1body = np.zeros((n_active, n_active)) Hamiltonian_fc_1body_tmp = np.zeros((n_active, n_active)) @@ -155,7 +170,12 @@ def spin_block_tei(I): ia = MSO_frozen_list[a] Hamiltonian_fc_1body[p, q] += gmo[ia, ip, ia, iq] if abs(Hamiltonian_fc_1body[p,q]) > 1e-12: - f_str += pos_or_neg(Hamiltonian_fc_1body[p,q]) + str(abs(Hamiltonian_fc_1body[p,q])) + ' ' + str(p) + '^ ' + str(q) + f_str += pos_or_neg(Hamiltonian_fc_1body[p, q]) + str( + abs(Hamiltonian_fc_1body[p, q])) + ' ' + if alternate_spins: + f_str += str(self.index_map[p]) + '^ ' + str(self.index_map[q]) + else: + f_str += str(p) + '^ ' + str(q) # ------- 2-body frozen-core: @@ -242,7 +262,15 @@ def spin_block_tei(I): for r in range(n_active): ir = MSO_active_list[r] for ss in range(n_active): - if abs(Hamiltonian_fc_2body_tmp[p,q,r,ss]) > 1e-12: - f_str += pos_or_neg(Hamiltonian_fc_2body_tmp[p,q,r,ss]) + str(abs(Hamiltonian_fc_2body_tmp[p,q,r,ss])) + ' ' + str(p) + '^ ' + str(q) + '^ ' + str(r) + ' ' + str(ss) + if abs(Hamiltonian_fc_2body_tmp[p, q, r, ss]) > 1e-12: + f_str += pos_or_neg(Hamiltonian_fc_2body_tmp[p, q, r, ss]) + str( + abs(Hamiltonian_fc_2body_tmp[p, q, r, ss])) + ' ' + if alternate_spins: + f_str += str(self.index_map[p]) + '^ ' + f_str += str(self.index_map[q]) + '^ ' + f_str += str(self.index_map[r]) + ' ' + f_str += str(self.index_map[ss]) + else: + f_str += str(p) + '^ ' + str(q) + '^ ' + str(r) + ' ' + str(ss) self.observable = xacc.getObservable('fermion', f_str) self.asPauli = xacc.transformToPauli('jw', self.observable) \ No newline at end of file diff --git a/python/plugins/observables/pyscf_observable.py b/python/plugins/observables/pyscf_observable.py index 556bd2d7d..1aed12936 100644 --- a/python/plugins/observables/pyscf_observable.py +++ b/python/plugins/observables/pyscf_observable.py @@ -13,6 +13,15 @@ def __init__(self): xacc.Observable.__init__(self) self.observable = None self.asPauli = None + self.index_map = None + + def _build_index_map(self, n_bits): + if self.index_map is None: + self.index_map = {} + for i in range(n_bits // 2): + self.index_map[i] = 2 * i + self.index_map[i + n_bits // 2] = 2 * i + 1 + return def toString(self): return self.observable.toString() @@ -33,7 +42,7 @@ def __iter__(self): return self.asPauli.__iter__() def fromOptions(self, inputParams): - import numpy as np, sys, io, os + import numpy as np, sys from pyscf import gto, scf, dft, tddft from pyscf.lib import logger @@ -46,20 +55,20 @@ def fromOptions(self, inputParams): else: mol.build(verbose=logger.QUIET) if 'spin' in inputParams and inputParams['spin'] != 0: - mol.spin = inputParams['spin'] - scf_wfn = scf.ROHF(mol) + mol.spin = inputParams['spin'] + scf_wfn = scf.ROHF(mol) else: - scf_wfn = scf.RHF(mol) # needs to be changed for open-shells + scf_wfn = scf.RHF(mol) # needs to be changed for open-shells scf_wfn.conv_tol = 1e-8 - scf_wfn.kernel() # runs RHF calculations + scf_wfn.kernel() # runs RHF calculations E_nucl = mol.energy_nuc() # Get orbital coefficients: Ca = scf_wfn.mo_coeff Cb = scf_wfn.mo_coeff C = np.block([ - [Ca,np.zeros_like(Cb)], - [np.zeros_like(Ca),Cb] + [Ca, np.zeros_like(Cb)], + [np.zeros_like(Ca), Cb] ]) # Get the two electron integrals using MintsHelper @@ -109,7 +118,8 @@ def spin_block_tei(I): # --- th is version is safer than above (H_1b is permutted correctly if eigvecs are permutted) Hamiltonian_1body_ao = np.block([[H_core_ao, np.zeros_like(H_core_ao)], [np.zeros_like(H_core_ao), H_core_ao]]) - Hamiltonian_1body = np.einsum('ij, jk, kl -> il', C.T, Hamiltonian_1body_ao, C) + Hamiltonian_1body = np.einsum( + 'ij, jk, kl -> il', C.T, Hamiltonian_1body_ao, C) Hamiltonian_2body = gmo if 'frozen-spin-orbitals' in inputParams and 'active-spin-orbitals' in inputParams: @@ -137,7 +147,13 @@ def spin_block_tei(I): f_str = str(Hamiltonian_fc_0body) - pos_or_neg = lambda x : ' + ' if x > 0. else ' - ' + def pos_or_neg(x): return ' + ' if x > 0. else ' - ' + + alternate_spins = False + if 'alternate-spins' in inputParams: + alternate_spins = inputParams['alternate-spins'] + if alternate_spins: + self._build_index_map(n_active) # --- 1-body frozen-core: Hamiltonian_fc_1body = np.zeros((n_active, n_active)) @@ -156,13 +172,18 @@ def spin_block_tei(I): ia = MSO_frozen_list[a] Hamiltonian_fc_1body[p, q] += gmo[ia, ip, ia, iq] - if abs(Hamiltonian_fc_1body[p,q]) > 1e-12: - f_str += pos_or_neg(Hamiltonian_fc_1body[p,q]) + str(abs(Hamiltonian_fc_1body[p,q])) + ' ' + str(p) + '^ ' + str(q) - + if abs(Hamiltonian_fc_1body[p, q]) > 1e-12: + f_str += pos_or_neg(Hamiltonian_fc_1body[p, q]) + str( + abs(Hamiltonian_fc_1body[p, q])) + ' ' + if alternate_spins: + f_str += str(self.index_map[p]) + '^ ' + str(self.index_map[q]) + else: + f_str += str(p) + '^ ' + str(q) # ------- 2-body frozen-core: - Hamiltonian_fc_2body = np.zeros((n_active, n_active, n_active, n_active)) + Hamiltonian_fc_2body = np.zeros( + (n_active, n_active, n_active, n_active)) for p in range(n_active): ip = MSO_active_list[p] @@ -179,7 +200,7 @@ def spin_block_tei(I): iss = MSO_active_list[ss] #Hamiltonian_fc_2body[p,q,r,ss]= 0.25* gmo[ip,iq,ir,iss] - Hamiltonian_fc_2body[p, q, r,ss] = gmo[ip, iq, ir, iss] + Hamiltonian_fc_2body[p, q, r, ss] = gmo[ip, iq, ir, iss] #Hamiltonian_fc_2body[p,q,r,ss]= 0.25* gmo[ip,iq,iss,ir] # checking whether to reduce Hamiltonian @@ -190,15 +211,16 @@ def spin_block_tei(I): if reduce_hamiltonian: if Hamiltonian_fc_1body.shape[0] != 4: - raise NotImplementedError("Reduction only implemented for Hamiltonians with 4 spin orbitals / 2 spatial orbitals.") + raise NotImplementedError( + "Reduction only implemented for Hamiltonians with 4 spin orbitals / 2 spatial orbitals.") # Gather pieces of the Hamiltonian - h11 = Hamiltonian_fc_1body[0,0] - h22 = Hamiltonian_fc_1body[1,1] - J11 = Hamiltonian_fc_2body[0,2,0,2] - J12 = Hamiltonian_fc_2body[1,2,1,2] - J22 = Hamiltonian_fc_2body[1,3,1,3] - K12 = Hamiltonian_fc_2body[0,2,1,3] + h11 = Hamiltonian_fc_1body[0, 0] + h22 = Hamiltonian_fc_1body[1, 1] + J11 = Hamiltonian_fc_2body[0, 2, 0, 2] + J12 = Hamiltonian_fc_2body[1, 2, 1, 2] + J22 = Hamiltonian_fc_2body[1, 3, 1, 3] + K12 = Hamiltonian_fc_2body[0, 2, 1, 3] e1 = h11 + J11 e2 = h22 + 2 * J12 - K12 @@ -213,7 +235,8 @@ def spin_block_tei(I): # <10|H|01> = K12 = g4 from numpy import linalg - M = np.array([[1, 1, 1, 1], [1, -1, 1, -1], [1, 1, -1, -1], [1, -1, -1, 1]]) + M = np.array([[1, 1, 1, 1], [1, -1, 1, -1], + [1, 1, -1, -1], [1, -1, -1, 1]]) Minv = linalg.inv(M) b = np.zeros(4) @@ -235,7 +258,8 @@ def spin_block_tei(I): else: - Hamiltonian_fc_2body_tmp = 0.25 * Hamiltonian_fc_2body.transpose(0, 1, 3, 2) + Hamiltonian_fc_2body_tmp = 0.25 * \ + Hamiltonian_fc_2body.transpose(0, 1, 3, 2) for p in range(n_active): ip = MSO_active_list[p] for q in range(n_active): @@ -243,7 +267,15 @@ def spin_block_tei(I): for r in range(n_active): ir = MSO_active_list[r] for ss in range(n_active): - if abs(Hamiltonian_fc_2body_tmp[p,q,r,ss]) > 1e-12: - f_str += pos_or_neg(Hamiltonian_fc_2body_tmp[p,q,r,ss]) + str(abs(Hamiltonian_fc_2body_tmp[p,q,r,ss])) + ' ' + str(p) + '^ ' + str(q) + '^ ' + str(r) + ' ' + str(ss) + if abs(Hamiltonian_fc_2body_tmp[p, q, r, ss]) > 1e-12: + f_str += pos_or_neg(Hamiltonian_fc_2body_tmp[p, q, r, ss]) + str( + abs(Hamiltonian_fc_2body_tmp[p, q, r, ss])) + ' ' + if alternate_spins: + f_str += str(self.index_map[p]) + '^ ' + f_str += str(self.index_map[q]) + '^ ' + f_str += str(self.index_map[r]) + ' ' + f_str += str(self.index_map[ss]) + else: + f_str += str(p) + '^ ' + str(q) + '^ ' + str(r) + ' ' + str(ss) self.observable = xacc.getObservable('fermion', f_str) self.asPauli = xacc.transformToPauli('jw', self.observable) \ No newline at end of file