diff --git a/helpers/Helpers.py b/helpers/Helpers.py index 24298d4..d2ffaef 100644 --- a/helpers/Helpers.py +++ b/helpers/Helpers.py @@ -5,6 +5,7 @@ ''' import os +from copy import deepcopy #=============================================================================== #%% Constants #=============================================================================== @@ -379,6 +380,106 @@ def prettyPrintDictionary(dictionary, level=0, delimiter=' . '): else: print(header+str(k)+':'+str(val)) + +#------------------------------------------------------------------------------ + +def __copyHamiltonian_4keyTo2keys(results_0, reverse=False): + """ Convert hamiltonian from 4-key <*bra, *ket> to 2-key (bra, ket) form""" + results = {} + if not reverse: + for key_, vals in results_0.items(): + bra_, ket_ = key_[:2], key_[2:] + if not bra_ in results: + results[bra_] = {ket_: vals, } + else: + results[bra_][ket_] = vals + else: + for bra_, kvals in results_0.items(): + for ket_, vals in kvals.items(): + results[(*bra_, *ket_)] = vals + return results + +def sortingHamiltonian(results, sorted_2b_comb, is_jt=False, l_ge_10=True): + """ + This function returns a hamiltonian in the order given, applying the phase + changes from the J or JT scheme. + + Note: use it before [recursiveSumOnDictionaries] + Input: + :results {(1, 1): {(101, 1): [{J,T: [values]}] ...}} + :sorted_2b_comb [(1,1), (1,101), (1,10001), ... ] + Returns: + the same result dictionary with the keys in the order of sorted_2b_comb + """ + dict_0 = {} + ## 0. Structure as the imported hamiltonian (keys involve (a,b,c,d):{J}) + if list(results.keys()) == []: return dict_0 + + _transform_to_key4 = False + if list(results.keys())[0].__len__() == 2: + dict_0 = __copyHamiltonian_4keyTo2keys(results, True) + _transform_to_key4 = True + else: + dict_0 = results + + # 1 read all the 4-keys and sort in the order bra-ket (assign permutat and particle order) + dict_1 = {} + for bk, vals in dict_0.items(): + permuts = getAllPermutationsOf2Bstates(bk, l_ge_10, not is_jt) + + for i, item_ in enumerate(permuts): + bk_perm, t_perm, phs = item_ + if phs == 0: continue ## redundant permutation + + not_in_ = not bk_perm in dict_1 + + if is_jt: + dict_1[bk_perm] = {0: {}, 1:{}} + for T in (0,1): + for J in vals: + phs2 = phs * (-1)**(J + T) + if not_in_: + dict_1[bk_perm][T][J] = phs2 * vals[T][J] + else: + # Test if matches + assert abs(dict_1[bk_perm][T][J] + - phs2 * vals[T][J]) < 1.0e-6,\ + "[ERROR]: values of:{} does not match "\ + "\with previous:{}".format(bk_perm, bk) + else: + dict_1[bk_perm] = dict([(J, dict()) for J in vals]) + for J in vals: + phs2 = phs + if i not in (0,4, 3,7): ## double/non exchange has no J dep. + phs2 = phs * (-1)**(J + 1) + for T in range(6): + if not_in_: + dict_1[bk_perm][J][T] = phs2 * vals[J][t_perm[T]] + else: + # Test if matches + # Test if matches + assert abs(dict_1[bk_perm][J][T] + - phs2 * vals[J][t_perm[T]]) < 1.0e-6,\ + "[ERROR]: values of:{} does not match "\ + "\with previous:{}".format(bk_perm, bk) + + # 2 sort in the order of bra (sorting_order) and apply the phs-changes + dict_2 = {} + for i in range(len(sorted_2b_comb)): + bra = sorted_2b_comb[i] + for j in range(i, len(sorted_2b_comb)): + ket = sorted_2b_comb[j] + + srt_key = (*bra, *ket) + if srt_key in dict_1: + dict_2[srt_key] = dict_1[srt_key] + + ## 0.2 In case of modifying to 4-key index: ---------- + if _transform_to_key4: + dict_2 = __copyHamiltonian_4keyTo2keys(dict_2) + return dict_2 + + def recursiveSumOnDictionaries(dict2read, dict2write): """ Given dictionaries with numerical end point values, sum the value to read if @@ -547,6 +648,56 @@ def readAntoine(index, l_ge_10=False): raise Exception("Invalid state index for Antoine Format [{}]".format(index)) +def getAllPermutationsOf2Bstates(tb_states, l_ge_10, is_j_scheme=True): + """ + Given all permutations of the quantum numbers < a,b| c,d> for pnpn states + return in this order: + [(0,1,2,3), (1,0,2,3), (0,1,3,2), (1,0,3,2), :: bra-ket + (2,3,0,1), (3,2,0,1), (2,3,1,0), (3,2,1,0) ] :: ket-bra + with their phases, if the quantum numbers are equal + * If it cannot be permuted, then phs = 0 + + OUTPUT: (all states permutations, particle-label perm associated, phs) + + WARNING: give the tb_states in the correct order (the one for export) + In case of JT-scheme, it will not be considered the t_perm + """ + _return_format2 = False + if len(tb_states) == 2: + _return_format2 = True + tb_states = (tb_states[0][0], tb_states[0][1], + tb_states[1][0], tb_states[1][1]) + + QQNN_PERMUTATIONS = [(0,1,2,3), (1,0,2,3), (0,1,3,2), (1,0,3,2), + (2,3,0,1), (3,2,0,1), (2,3,1,0), (3,2,1,0)] + PART_LABEL_PERMUT = [(1,2,3,4), (3,4,1,2), (2,1,4,3), (4,3,2,1), + (1,3,2,4), (3,1,4,2), (2,4,1,3), (4,2,3,1)] + + permuts = [] + perv_permuts = set() + for i, key_perm in enumerate(QQNN_PERMUTATIONS): + phs = [1, 1] + if i in (1, 3, 6, 7): + phs[0] = (-1)**((readAntoine(tb_states[0], l_ge_10)[2] + + readAntoine(tb_states[1], l_ge_10)[2]) // 2) + if i in (2, 3, 5, 7): + phs[1] = (-1)**((readAntoine(tb_states[2], l_ge_10)[2] + + readAntoine(tb_states[3], l_ge_10)[2]) // 2) + + phs = phs[0] * phs[1] + + st_perm = tuple([tb_states[k] for k in key_perm]) + t_perm = tuple([0, *PART_LABEL_PERMUT[i], 5]) if is_j_scheme else None + if st_perm in perv_permuts: + phs = 0 + perv_permuts.add(st_perm) + permuts.append( (st_perm, t_perm, phs)) + + if _return_format2: + permuts = [] + for bk, t_perm, phs in permuts: + permuts.append( ((bk[0], bk[1]), (bk[2], bk[3]), t_perm, phs) ) + return permuts def shellSHO_Notation(n, l, j=0): diff --git a/helpers/TBME_Runner.py b/helpers/TBME_Runner.py index 3c2f929..2a27d19 100644 --- a/helpers/TBME_Runner.py +++ b/helpers/TBME_Runner.py @@ -20,7 +20,8 @@ QN_2body_jj_J_Coupling from copy import deepcopy from helpers.Helpers import recursiveSumOnDictionaries, getCoreNucleus,\ - Constants, almostEqual, printProgressBar, _LINE_1, _LINE_2 + Constants, almostEqual, printProgressBar, _LINE_1, _LINE_2,\ + sortingHamiltonian import os import numpy as np @@ -375,8 +376,9 @@ def _sortQQNNFromTheValenceSpace(self): q_numbs = map(lambda qn: int(qn), q_numbs) q_numbs = sorted(q_numbs)#, reverse=True) - q_numbs = list(combinations_with_replacement(q_numbs, 2)) + self._valence_space_sorted = q_numbs ## to keep the order of exporting + q_numbs = list(combinations_with_replacement(q_numbs, 2)) self._twoBodyQuantumNumbersSorted = q_numbs all_permut_ = [] @@ -386,7 +388,6 @@ def _sortQQNNFromTheValenceSpace(self): all_permut_.append((b,a)) all_permut_.sort() self._allPermutations_twoBodyQuantumNumbers = all_permut_ - self._allPermutations_twoBodyQuantumNumbers = self._twoBodyQuantumNumbersSorted def _compute1BodyMatrixElements(self, kin=False, force=None): """ @@ -576,7 +577,7 @@ def _readMatrixElementsFromFile(self, force_str, **params): if l_ge_10: l_ge_10 = False if l_ge_10.lower() == 'false' else True - valence_space = list(self.input_obj.Valence_Space.keys()) + valence_space = self._valence_space_sorted data_ = TBME_Reader(filename, ignorelines, constant, valence_space, l_ge_10) @@ -595,12 +596,15 @@ def _readMatrixElementsFromFile(self, force_str, **params): def combineAllResults(self): + ## Get all matrix elements in the correct order. + kin_key = ForceEnum.Kinetic_2Body final_J = {} final_JT = {} for force, results in self.resultsByInteraction.items(): - if self.__class__ == TBME_Runner: - force = self._forceNameFromForceEnum[force] + ## NOTE: Why this was here? + # if self.__class__ == TBME_Runner: + # force = self._forceNameFromForceEnum[force] if self.interactionSchemes[force] == self._Scheme.J: # Kin 2Body is a JT scheme interaction diff --git a/helpers/io_manager.py b/helpers/io_manager.py index 6d670eb..5513a8d 100644 --- a/helpers/io_manager.py +++ b/helpers/io_manager.py @@ -10,7 +10,8 @@ import xml.etree.ElementTree as et from helpers.WaveFunctions import QN_1body_jj from helpers.Helpers import readAntoine, shellSHO_Notation, valenceSpacesDict_l_ge10,\ - getCoreNucleus + getCoreNucleus, getAllPermutationsOf2Bstates, prettyPrintDictionary,\ + sortingHamiltonian from copy import deepcopy @@ -595,6 +596,7 @@ def __init__(self, filename, :constant , Factor to multiply all the matrix elements :val_space , when given, it skips matrix elements that contain other quantum numbers. + :sorting_2b_order , order of two-body wf, [(a,b), ...] """ with open(filename, 'r') as f: data = f.readlines() @@ -606,7 +608,8 @@ def __init__(self, filename, if filename.endswith(OutputFileTypes.sho): self.scheme = self._Scheme.JT - elif filename.endswith(OutputFileTypes.twoBody): + elif (filename.endswith(OutputFileTypes.twoBody) + or filename.endswith(OutputFileTypes.centerOfMass)): self.scheme = self._Scheme.J else: raise ParserException("Unidentified Scheme for importing") @@ -614,6 +617,7 @@ def __init__(self, filename, self._elements = {} self._valence_space = [] self._strict_val_space = False + self._sorting_2b_order = [] if val_space: self._strict_val_space = True self._valence_space = [int(st) for st in val_space] @@ -623,33 +627,17 @@ def __init__(self, filename, def getMatrixElemnts(self, sorting_order=None): """ Get the imported matrix elements, : sorting order: of tuples with the qqnn in """ - + print("\n\n Exporting matrix elements READED:\n\n *****************") if sorting_order: - dict_2 = {} - q_numbs = sorting_order - for i in range(len(q_numbs)): - bra = q_numbs[i] - - if not bra in self._elements: - continue - else: - if not bra in dict_2: - dict_2[q_numbs[i]] = {} + args = ( + self._elements, + sorting_order, + self.scheme == self._Scheme.JT, + self.l_ge_10 + ) + self._elements = sortingHamiltonian(*args) - for j in range(0, len(q_numbs)): # start i = 0 - ket = q_numbs[j] - if not ket in self._elements[bra]: - continue - else: - if not ket in dict_2[bra]: - dict_2[bra][ket] = {} - - block = self._elements[q_numbs[i]][q_numbs[j]] - dict_2[q_numbs[i]][q_numbs[j]] = block - - self._elements = deepcopy(dict_2) - - return deepcopy(self._elements) + return self._elements def _importMatrixElements(self, data): if self.scheme == self._Scheme.JT: @@ -657,53 +645,43 @@ def _importMatrixElements(self, data): elif self.scheme == self._Scheme.J: self._getJSchemeMatrixElements(data) - def _getJTSchemeMatrixElements(self, data): - JT_block = {} bra, ket = None, None - index = 0 - j_min, j_max, T = 0, 0, 0 + index = 0 + j_min, j_max, J = 0, 0, 0 for line in data: line = line.strip() if index == 0: bra, ket, j_min, j_max, skip = self._getBlockQN_HamilType(line) - + j_length = j_max - j_min + 1 if skip: index += 1 continue - # JT_block[bra] = {0: {}, 1: {}} - phs_bra, phs_ket, bra, ket = self._getPermutationPhase(bra, ket) + + j_dict = dict([(j, {}) for j in range(j_min, j_max+1)]) + if bra in JT_block: if ket in JT_block[bra]: print("WARNING, while reading the file found repeated " "matrix element block:", bra, ket) - JT_block[bra][ket] = {0: {}, 1: {}} + JT_block[bra][ket] = deepcopy(j_dict) else: - JT_block[bra] = {ket : {0: {}, 1: {}}} + JT_block[bra] = {ket : deepcopy(j_dict)} elif index > 0 and (not skip): T = index - 1 line = line.split() for i in range(j_max - j_min +1): J = j_min + i - - phs = 1 - if self.bra_exch: - phs *= (-1)**(phs_bra + J + T) - if self.ket_exch: - phs *= (-1)**(phs_ket + J + T) - - JT_block[bra][ket][T][J] = self.const * phs * float(line[i]) - - index = index + 1 if index < 2 else 0 + JT_block[bra][ket][T][J] = self.const * float(line[i]) + + index = index + 1 if index < j_length else 0 self._elements = JT_block - def _getJSchemeMatrixElements(self, data): - J_block = {} bra, ket = None, None index = 0 @@ -720,50 +698,58 @@ def _getJSchemeMatrixElements(self, data): j_dict = dict([(j, {}) for j in range(j_min, j_max+1)]) - phs_bra, phs_ket, bra, ket = self._getPermutationPhase(bra, ket) if bra in J_block: if ket in J_block[bra]: print("WARNING, while reading the file found repeated " "matrix element block:", bra, ket) J_block[bra][ket] = deepcopy(j_dict) else: - J_block[bra] = {ket : deepcopy(j_dict)} - - t_indexing = self._getParticleLabelIndexingAfterPermutation() - + J_block[bra] = {ket : deepcopy(j_dict)} elif index > 0 and (not skip): J = j_min + index - 1 - line = line.split() - phs = 1 - if self.bra_exch: - phs *= (-1)**(phs_bra + J + 1) - if self.ket_exch: - phs *= (-1)**(phs_ket + J + 1) - for T in t_indexing: - J_block[bra][ket][J][T] = self.const * phs * float(line[T]) - + for i in range(6): + J_block[bra][ket][J][i] = self.const * float(line[i]) + index = index + 1 if index < j_length else 0 self._elements = J_block def _getParticleLabelIndexingAfterPermutation(self): """ When bra or ket are exchanged, the index of the particle - lables have to be reorganized """ + lables have to be reorganized + NOTE: the permutation only considers the bra or ket exchange + but cannot tell the switch ket-bra for the v.s. order. + this permutation does not affect the value. + """ + t_indexing = [0, 1, 2, 3, 4, 5] - # pppp pnpn pnnp nppn npnp nnnn - if self.bra_exch: + # pppp pnpn pnnp nppn npnp nnnn + if self.bra_ket_switch: + t_indexing = [0, 1, 3, 2, 4, 5] + # pppp pnpn nppn pnnp npnp nnnn + if self.bra_exch: + if self.ket_exch: + # pppp npnp pnnp nppn pnpn nnnn + t_indexing = [0, 4, 2, 3, 1, 5] + else: + # pppp pnnp npnp pnpn nppn nnnn + t_indexing = [0, 2, 4, 1, 3, 5] if self.ket_exch: - # pppp npnp nppn pnnp pnpn nnnn - t_indexing = [0, 4, 3, 2, 1, 5] - else: - # pppp nppn npnp pnpn pnnp nnnn - t_indexing = [0, 3, 4, 1, 2, 5] - if self.ket_exch: - # pppp pnnp pnpn npnp nppn nnnn - t_indexing = [0, 2, 1, 4, 3, 5] - + # pppp nppn pnpn npnp pnnp nnnn + t_indexing = [0, 3, 1, 4, 2, 5] + else: + if self.bra_exch: + if self.ket_exch: + # pppp npnp nppn pnnp pnpn nnnn + t_indexing = [0, 4, 3, 2, 1, 5] + else: + # pppp nppn npnp pnpn pnnp nnnn + t_indexing = [0, 3, 4, 1, 2, 5] + if self.ket_exch: + # pppp pnnp pnpn npnp nppn nnnn + t_indexing = [0, 2, 1, 4, 3, 5] return t_indexing def _getBlockQN_HamilType(self, header): @@ -785,39 +771,12 @@ def _getBlockQN_HamilType(self, header): skip = True else: self._valence_space.append(st) - - bra = min(tuple(spss[0:2]), tuple(spss[2:4])) - ket = max(tuple(spss[0:2]), tuple(spss[2:4])) - + + bra, ket = tuple(spss[0:2]), tuple(spss[2:4]) + ## NOTE: do not sort here bra-ket, since it has a particle-label in the + ## J-scheme related to the possition in the imported line. return bra, ket, int(header[-2]), int(header[-1]), skip - def _getPermutationPhase(self, bra, ket): - """ - return the bra and the ket in increasing order (avoids repetition). - - phases do not have the (J + 1) or (J + T) part (append afterwards), - just j1 + j2 - """ - - self.bra_exch = False - self.ket_exch = False - phs_bra, phs_ket = 0, 0 - - if bra[0] > bra[1]: - self.bra_exch = True - bra = (bra[1], bra[0]) - phs_bra = (readAntoine(bra[0], self.l_ge_10)[2] + - readAntoine(bra[1], self.l_ge_10)[2]) // 2 - if ket[0] > ket[1]: - self.ket_exch = True - ket = (ket[1], ket[0]) - phs_ket = (readAntoine(ket[0], self.l_ge_10)[2] + - readAntoine(ket[1], self.l_ge_10)[2]) // 2 - - return phs_bra, phs_ket, bra, ket - - - #=============================================================================== # AUXIALIRY METHODS FOR EXPORTATION, #===============================================================================