Skip to content

Commit

Permalink
Merge pull request #34 from danclaudino/alternate_spins
Browse files Browse the repository at this point in the history 
Added alternate spins ordering to chemistry Observable plugins
  • Loading branch information
@danclaudino
danclaudino authored Aug 9, 2024
2 parents ef919eb + 5eadf39 commit d6b9249
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 28 deletions.
34 changes: 31 additions & 3 deletions python/plugins/observables/psi4_observable.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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))
Expand All @@ -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:
Expand Down Expand Up @@ -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)
82 changes: 57 additions & 25 deletions python/plugins/observables/pyscf_observable.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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))
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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)
Expand All @@ -235,15 +258,24 @@ 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):
iq = MSO_active_list[q]
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)

0 comments on commit d6b9249

Please sign in to comment.