Skip to content

Commit c5b8222

Browse files
committed
add in k vec solving routine
1 parent c80528d commit c5b8222

File tree

3 files changed

+14444
-218
lines changed

3 files changed

+14444
-218
lines changed

GSASII/GSASIIpwdGUI.py

Lines changed: 7 additions & 218 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@
3939
from . import G2shapes
4040
from . import SUBGROUPS as kSUB
4141
from . import k_vector_search as kvs
42+
from . import k_vec_solve as kvsolve
4243
from GSASII.imports.G2phase_CIF import CIFPhaseReader as CIFpr
4344
from . import GSASIIscriptable as G2sc
4445
from . import GSASIImiscGUI as G2IO
@@ -5215,38 +5216,9 @@ def grab_all_kvecs(out_html):
52155216

52165217
all_kvecs = {}
52175218
for key in final_kvec_keys:
5218-
pattern_str = key.split("(")[1].split(")")[0]
5219-
a_val = pattern_str.split(",")[0]
5220-
b_val = pattern_str.split(",")[1]
5221-
c_val = pattern_str.split(",")[2]
5222-
if a_val == "1/2":
5223-
a_val = Fraction(1, 2)
5224-
elif a_val == "1/3":
5225-
a_val = Fraction(1, 3)
5226-
elif a_val == "1":
5227-
a_val = 1
5228-
elif a_val == "0":
5229-
a_val = 0
5230-
5231-
if b_val == "1/2":
5232-
b_val = Fraction(1, 2)
5233-
elif b_val == "1/3":
5234-
b_val = Fraction(1, 3)
5235-
elif b_val == "1":
5236-
b_val = 1
5237-
elif b_val == "0":
5238-
b_val = 0
5239-
5240-
if c_val == "1/2":
5241-
c_val = Fraction(1, 2)
5242-
elif c_val == "1/3":
5243-
c_val = Fraction(1, 3)
5244-
elif c_val == "1":
5245-
c_val = 1
5246-
elif c_val == "0":
5247-
c_val = 0
5248-
5249-
all_kvecs[key] = (a_val, b_val, c_val)
5219+
pattern_str = f"({key.split("(")[1]}"
5220+
k_form = kvsolve.parse_expression_string(pattern_str)
5221+
all_kvecs[key] = k_form
52505222

52515223
return all_kvecs
52525224
else:
@@ -5265,138 +5237,7 @@ def setup_kvec_input(k_vec, k_vec_dict, symmetry=None):
52655237
dict: New entries and those need to be corrected in the data
52665238
to be used in the post request.
52675239
"""
5268-
from itertools import product
5269-
5270-
def extract_coeff_and_var(s):
5271-
"""Extract coefficient and var from strings like '2/3a', '1/3b', '1/2', 'a', etc.
5272-
"""
5273-
s = s.strip()
5274-
5275-
if not s:
5276-
return Fraction(1), ''
5277-
5278-
# Split into coefficient and character parts
5279-
# Find the last alphabetic character
5280-
char_match = re.search(r'[a-zA-Z]', s)
5281-
5282-
if char_match:
5283-
# There's a character
5284-
char_pos = char_match.start()
5285-
coeff_part = s[:char_pos]
5286-
char_part = s[char_pos:]
5287-
5288-
# Validate that character part is just one letter
5289-
if not re.match(r'^[a-zA-Z]$', char_part):
5290-
print(f"Error: Invalid character part: {char_part}")
5291-
else:
5292-
# No character
5293-
coeff_part = s
5294-
char_part = ''
5295-
5296-
# Parse coefficient
5297-
if not coeff_part or coeff_part in ['+', '']:
5298-
coefficient = Fraction(1)
5299-
elif coeff_part == '-':
5300-
coefficient = Fraction(-1)
5301-
else:
5302-
try:
5303-
coefficient = Fraction(coeff_part)
5304-
except ValueError:
5305-
try:
5306-
coefficient = Fraction(float(coeff_part)).limit_denominator()
5307-
except ValueError as e:
5308-
print(f"Error converting coefficient: {e}")
5309-
5310-
return coefficient, char_part
5311-
5312-
def generate_all_combinations(expressions, max_var_value=100):
5313-
"""Generate all possible combinations of variable assignments and calculate sums.
5314-
5315-
Args:
5316-
expressions (list): List of strings like ['2/3a', '1/3b', '1/2', 'a']
5317-
max_var_value (int): Maximum value to assign to variables (1 to max_var_value)
5318-
5319-
Returns:
5320-
list: List of tuples (sum_value, variable_assignments)
5321-
"""
5322-
parsed_terms = []
5323-
variables = set()
5324-
5325-
for expr in expressions:
5326-
coeff, var = extract_coeff_and_var(expr)
5327-
parsed_terms.append((coeff, var))
5328-
if var:
5329-
variables.add(var)
5330-
5331-
variables = sorted(list(variables))
5332-
5333-
if not variables:
5334-
total_sum = sum(coeff for coeff, var in parsed_terms if not var)
5335-
yield (total_sum, {})
5336-
return
5337-
5338-
# Generate combinations one at a time
5339-
var_ranges = [range(1, max_var_value + 1) for _ in variables]
5340-
5341-
for var_values in product(*var_ranges):
5342-
var_assignment = dict(zip(variables, var_values))
5343-
5344-
total_sum = Fraction(0)
5345-
for coeff, var in parsed_terms:
5346-
if var:
5347-
total_sum += coeff * var_assignment[var]
5348-
else:
5349-
total_sum += coeff
5350-
5351-
yield (total_sum, var_assignment)
5352-
5353-
def get_unique_sums(expressions, max_var_value=100):
5354-
"""
5355-
Get unique sum values and count their occurrences
5356-
5357-
Returns:
5358-
dict: {sum_value: count}
5359-
"""
5360-
sum_counts = {}
5361-
5362-
for total_sum, var_assignment in generate_all_combinations(expressions, max_var_value):
5363-
if total_sum in sum_counts:
5364-
sum_counts[total_sum] += 1
5365-
else:
5366-
sum_counts[total_sum] = 1
5367-
5368-
return sum_counts
5369-
5370-
def match_vector_pattern(k_vec, k_vec_dict, symmetry=None):
5371-
"""Check the k-vector against the standard form in isodistort.
5372-
5373-
Args:
5374-
k_vec (str): The k-vector to use for the isodistort request.
5375-
k_vec_dict (dict): The standard k-vector form in isodistort.
5376-
symmetry (str, optional): The crystal system.
5377-
5378-
Returns:
5379-
str: The standard k-vector form in isodistort.
5380-
"""
5381-
from itertools import permutations
5382-
5383-
all_matches = {}
5384-
5385-
for desc, pattern in k_vec_dict.items():
5386-
if len(pattern) != len(k_vec):
5387-
continue
5388-
5389-
if symmetry in ['cubic', 'rhombohedral']:
5390-
k_vec_sequences = list(permutations(k_vec))
5391-
elif symmetry in ['hexagonal', 'trigonal', 'tetragonal']:
5392-
k_vec_sequences = [k_vec, (k_vec[1], k_vec[0], k_vec[2])]
5393-
else:
5394-
k_vec_sequences = [k_vec]
5395-
5396-
for k_vec_seq in k_vec_sequences:
5397-
pass
5398-
5399-
return
5240+
k_forms = k_vec_dict.items()
54005241

54015242
k_vec_form = match_vector_pattern(k_vec, k_vec_dict, symmetry="cubic")
54025243

@@ -5506,68 +5347,16 @@ def match_vector_pattern(k_vec, k_vec_dict, symmetry=None):
55065347
except ValueError:
55075348
break
55085349

5509-
5510-
# The following chunk of code is for converting the k-vector from the
5511-
# conventional setting to the primitive setting.
55125350
phase_sel = G2frame.kvecSearch['phase']
55135351
_, Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
55145352
Phase = Phases[phase_sel]
55155353

5516-
lat_type = Phase["General"]["SGData"]["SGLatt"]
5517-
lat_sym = Phase["General"]["SGData"]["SGSys"]
5518-
if lat_sym == "trigonal":
5519-
brav_sym = "hR"
5520-
else:
5521-
brav_sym = lat_sym[0] + lat_type
5522-
5523-
Trans = np.eye(3)
5524-
Uvec = np.zeros(3)
5525-
Vvec = np.zeros(3)
5526-
5527-
newPhase = copy.deepcopy(Phase)
5528-
newPhase['ranId'] = ran.randint(0, sys.maxsize)
5529-
newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
5530-
newPhase, _ = G2lat.TransformPhase(Phase, newPhase, Trans,
5531-
Uvec, Vvec, False)
5532-
atoms_pointer = newPhase['General']['AtomPtrs']
5533-
5534-
atom_coords = list()
5535-
atom_types = list()
5536-
for atom in newPhase["Atoms"]:
5537-
coord_tmp = atom[atoms_pointer[0]:atoms_pointer[0] + 3]
5538-
atom_coords.append(coord_tmp)
5539-
type_tmp = atom[atoms_pointer[1]]
5540-
atom_types.append(type_tmp)
5541-
5542-
atom_ids = kvs.unique_id_gen(atom_types)
5543-
5544-
cell_params = newPhase["General"]["Cell"][1:7]
5545-
lat_vectors = kvs.lat_params_to_vec(cell_params)
5546-
5547-
hkl_refls = list()
5548-
for i in range(6):
5549-
for j in range(6):
5550-
for k in range(6):
5551-
hkl_refls.append([i, j, k])
5552-
5553-
k_search = kvs.kVector(
5554-
brav_sym, lat_vectors,
5555-
atom_coords, atom_ids,hkl_refls,
5556-
[0.], 0.
5557-
)
5558-
kpoint = k_search.hklConvToPrim(kpoint)
5559-
kpoint_frac = (
5560-
Fraction(kpoint[0]).limit_denominator(10),
5561-
Fraction(kpoint[1]).limit_denominator(10),
5562-
Fraction(kpoint[2]).limit_denominator(10)
5563-
)
5564-
55655354
data['input'] = 'kvector'
55665355
data['irrepcount'] = '1'
55675356

55685357
kvec_dict = grab_all_kvecs(out2)
5569-
5570-
data_update = setup_kvec_input(kpoint_frac, kvec_dict, symmetry=lat_sym)
5358+
lat_sym = Phase['General']['SGData']
5359+
data_update = setup_kvec_input(kpoint, kvec_dict, symmetry=lat_sym)
55715360
for key, value in data_update.items():
55725361
data[key] = value
55735362

0 commit comments

Comments
 (0)