| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | ########### SVN repository information ################### |
|---|
| 3 | # $Date: 2023-05-11 23:08:12 +0000 (Thu, 11 May 2023) $ |
|---|
| 4 | # $Author: toby $ |
|---|
| 5 | # $Revision: 5577 $ |
|---|
| 6 | # $URL: trunk/imports/G2phase_CIF.py $ |
|---|
| 7 | # $Id: G2phase_CIF.py 5577 2023-05-11 23:08:12Z toby $ |
|---|
| 8 | ########### SVN repository information ################### |
|---|
| 9 | ''' |
|---|
| 10 | ''' |
|---|
| 11 | # Routines to import Phase information from CIF files |
|---|
| 12 | from __future__ import division, print_function |
|---|
| 13 | import sys |
|---|
| 14 | import random as ran |
|---|
| 15 | import numpy as np |
|---|
| 16 | import re |
|---|
| 17 | import copy |
|---|
| 18 | import os.path |
|---|
| 19 | import GSASIIobj as G2obj |
|---|
| 20 | import GSASIIspc as G2spc |
|---|
| 21 | import GSASIIElem as G2elem |
|---|
| 22 | import GSASIIlattice as G2lat |
|---|
| 23 | import GSASIIpath |
|---|
| 24 | import GSASIIfiles as G2fil |
|---|
| 25 | try: |
|---|
| 26 | import GSASIIctrlGUI as G2G |
|---|
| 27 | except ImportError: |
|---|
| 28 | pass |
|---|
| 29 | GSASIIpath.SetVersionNumber("$Revision: 5577 $") |
|---|
| 30 | import CifFile as cif # PyCifRW from James Hester |
|---|
| 31 | debug = GSASIIpath.GetConfigValue('debug') |
|---|
| 32 | #debug = False |
|---|
| 33 | |
|---|
| 34 | class CIFPhaseReader(G2obj.ImportPhase): |
|---|
| 35 | 'Implements a phase importer from a possibly multi-block CIF file' |
|---|
| 36 | def __init__(self): |
|---|
| 37 | super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__ |
|---|
| 38 | extensionlist=('.CIF','.cif','.mcif'), |
|---|
| 39 | strictExtension=False, |
|---|
| 40 | formatName = 'CIF', |
|---|
| 41 | longFormatName = 'Crystallographic Information File import' |
|---|
| 42 | ) |
|---|
| 43 | |
|---|
| 44 | def ContentsValidator(self, filename): |
|---|
| 45 | fp = open(filename,'r') |
|---|
| 46 | ok = self.CIFValidator(fp) |
|---|
| 47 | fp.close() |
|---|
| 48 | return ok |
|---|
| 49 | |
|---|
| 50 | def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused): |
|---|
| 51 | isodistort_warnings = '' # errors that would prevent an isodistort analysis |
|---|
| 52 | self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 53 | # make sure the ranId is really unique! |
|---|
| 54 | while self.Phase['ranId'] in usedRanIdList: |
|---|
| 55 | self.Phase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 56 | self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 57 | # make sure the ranId is really unique! |
|---|
| 58 | while self.MPhase['ranId'] in usedRanIdList: |
|---|
| 59 | self.MPhase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 60 | returnstat = False |
|---|
| 61 | cellitems = ( |
|---|
| 62 | '_cell_length_a','_cell_length_b','_cell_length_c', |
|---|
| 63 | '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',) |
|---|
| 64 | # cellwaveitems = ( |
|---|
| 65 | # '_cell_wave_vector_seq_id', |
|---|
| 66 | # '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z') |
|---|
| 67 | reqitems = ( |
|---|
| 68 | '_atom_site_fract_x', |
|---|
| 69 | '_atom_site_fract_y', |
|---|
| 70 | '_atom_site_fract_z', |
|---|
| 71 | ) |
|---|
| 72 | phasenamefields = ( |
|---|
| 73 | '_chemical_name_common', |
|---|
| 74 | '_pd_phase_name', |
|---|
| 75 | '_chemical_formula_sum' |
|---|
| 76 | ) |
|---|
| 77 | try: |
|---|
| 78 | cf = G2obj.ReadCIF(filename) |
|---|
| 79 | except cif.StarError as msg: |
|---|
| 80 | msg = 'Unreadable cif file\n'+str(msg) |
|---|
| 81 | self.errors = msg |
|---|
| 82 | self.warnings += msg |
|---|
| 83 | return False |
|---|
| 84 | # scan blocks for structural info |
|---|
| 85 | self.errors = 'Error during scan of blocks for datasets' |
|---|
| 86 | str_blklist = [] |
|---|
| 87 | for blk in cf.keys(): |
|---|
| 88 | for r in reqitems+cellitems: |
|---|
| 89 | if r not in cf[blk].keys(): |
|---|
| 90 | break |
|---|
| 91 | else: |
|---|
| 92 | str_blklist.append(blk) |
|---|
| 93 | if not str_blklist: |
|---|
| 94 | selblk = None # no block to choose |
|---|
| 95 | elif len(str_blklist) == 1: # only one choice |
|---|
| 96 | selblk = 0 |
|---|
| 97 | else: # choose from options |
|---|
| 98 | choice = [] |
|---|
| 99 | for blknm in str_blklist: |
|---|
| 100 | choice.append('') |
|---|
| 101 | # accumumlate some info about this phase |
|---|
| 102 | choice[-1] += blknm + ': ' |
|---|
| 103 | for i in phasenamefields: # get a name for the phase |
|---|
| 104 | name = cf[blknm].get(i,'phase name').strip() |
|---|
| 105 | if name is None or name == '?' or name == '.': |
|---|
| 106 | continue |
|---|
| 107 | else: |
|---|
| 108 | choice[-1] += name.strip() + ', ' |
|---|
| 109 | break |
|---|
| 110 | na = len(cf[blknm].get("_atom_site_fract_x")) |
|---|
| 111 | if na == 1: |
|---|
| 112 | choice[-1] += '1 atom' |
|---|
| 113 | else: |
|---|
| 114 | choice[-1] += ('%d' % na) + ' atoms' |
|---|
| 115 | choice[-1] += ', cell: ' |
|---|
| 116 | fmt = "%.2f," |
|---|
| 117 | for i,key in enumerate(cellitems): |
|---|
| 118 | if i == 3: fmt = "%.f," |
|---|
| 119 | if i == 5: fmt = "%.f" |
|---|
| 120 | choice[-1] += fmt % cif.get_number_with_esd( |
|---|
| 121 | cf[blknm].get(key))[0] |
|---|
| 122 | sg = cf[blknm].get("_symmetry_space_group_name_H-M",'') |
|---|
| 123 | if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'') |
|---|
| 124 | if not sg: sg = cf[blknm].get("_space_group_ssg_name",'') |
|---|
| 125 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 126 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'') |
|---|
| 127 | #how about checking for super/magnetic ones as well? - reject 'X'? |
|---|
| 128 | sg = sg.replace('_','') |
|---|
| 129 | if sg: choice[-1] += ', (' + sg.strip() + ')' |
|---|
| 130 | selblk = G2G.PhaseSelector(choice,ParentFrame=ParentFrame, |
|---|
| 131 | title= 'Select a phase from one the CIF data_ blocks below',size=(600,100)) |
|---|
| 132 | self.errors = 'Error during reading of selected block' |
|---|
| 133 | #process selected phase |
|---|
| 134 | if selblk is None: |
|---|
| 135 | returnstat = False # no block selected or available |
|---|
| 136 | else: #do space group symbol & phase type first |
|---|
| 137 | blknm = str_blklist[selblk] |
|---|
| 138 | blk = cf[str_blklist[selblk]] |
|---|
| 139 | E = True |
|---|
| 140 | Super = False |
|---|
| 141 | magnetic = False |
|---|
| 142 | moddim = int(blk.get("_cell_modulation_dimension",'0')) |
|---|
| 143 | if moddim: #incommensurate |
|---|
| 144 | if moddim > 1: |
|---|
| 145 | msg = 'more than 3+1 super symmetry is not allowed in GSAS-II' |
|---|
| 146 | self.errors = msg |
|---|
| 147 | self.warnings += '\n'+msg |
|---|
| 148 | return False |
|---|
| 149 | if blk.get('_cell_subsystems_number'): |
|---|
| 150 | msg = 'Composite super structures not allowed in GSAS-II' |
|---|
| 151 | self.errors = msg |
|---|
| 152 | self.warnings += '\n'+msg |
|---|
| 153 | return False |
|---|
| 154 | sspgrp = blk.get("_space_group_ssg_name",'') |
|---|
| 155 | if not sspgrp: #might be incommensurate magnetic |
|---|
| 156 | MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 157 | if not MSSpGrp: |
|---|
| 158 | MSSpGrp = blk.get("_space_group.magn_ssg_name",'') |
|---|
| 159 | if not MSSpGrp: |
|---|
| 160 | msg = 'No incommensurate space group name was found in the CIF.' |
|---|
| 161 | self.errors = msg |
|---|
| 162 | self.warnings += '\n'+msg |
|---|
| 163 | return False |
|---|
| 164 | if 'X' in MSSpGrp: |
|---|
| 165 | msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II' |
|---|
| 166 | self.warnings += '\n'+msg |
|---|
| 167 | self.errors = msg |
|---|
| 168 | return False |
|---|
| 169 | magnetic = True |
|---|
| 170 | if 'X' in sspgrp: |
|---|
| 171 | msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II' |
|---|
| 172 | self.warnings += '\n'+msg |
|---|
| 173 | self.errors = msg |
|---|
| 174 | return False |
|---|
| 175 | Super = True |
|---|
| 176 | if magnetic: |
|---|
| 177 | sspgrp = MSSpGrp.split('(') |
|---|
| 178 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 179 | SpGrp = G2spc.StandardizeSpcName(sspgrp[0]) |
|---|
| 180 | if "1'" in SpGrp: sspgrp[1] = sspgrp[1][:-1] #take off extra 's'; gets put back later |
|---|
| 181 | MSpGrp = sspgrp[0] |
|---|
| 182 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 183 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 184 | else: |
|---|
| 185 | sspgrp = sspgrp.split('(') |
|---|
| 186 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 187 | SpGrp = sspgrp[0] |
|---|
| 188 | if "1'" in SpGrp: #nonmagnetics can't be gray |
|---|
| 189 | SpGrp = SpGrp.replace("1'",'') |
|---|
| 190 | sspgrp[1] = sspgrp[1][:-1] #take off extra 's' |
|---|
| 191 | # SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 192 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 193 | if not SpGrp: |
|---|
| 194 | print (sspgrp) |
|---|
| 195 | self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II' |
|---|
| 196 | return False |
|---|
| 197 | SuperSg = sspgrp[1].replace('\\','').replace(',','') |
|---|
| 198 | SuperVec = [[0,0,.1],False,4] |
|---|
| 199 | else: #not incommensurate |
|---|
| 200 | SpGrp = blk.get("_symmetry_space_group_name_H-M",'') |
|---|
| 201 | if not SpGrp: |
|---|
| 202 | SpGrp = blk.get("_space_group_name_H-M_alt",'') |
|---|
| 203 | if not SpGrp: #try magnetic |
|---|
| 204 | MSpGrp = blk.get("_space_group.magn_name_BNS",'') |
|---|
| 205 | if not MSpGrp: |
|---|
| 206 | MSpGrp = blk.get("_space_group_magn.name_BNS",'') |
|---|
| 207 | if not MSpGrp: |
|---|
| 208 | msg = 'No recognizable space group name was found in the CIF.' |
|---|
| 209 | self.errors = msg |
|---|
| 210 | self.warnings += '\n'+msg |
|---|
| 211 | return False |
|---|
| 212 | SpGrp = blk.get('_parent_space_group.name_H-M_alt') |
|---|
| 213 | if not SpGrp: |
|---|
| 214 | SpGrp = blk.get('_parent_space_group.name_H-M') |
|---|
| 215 | # SpGrp = MSpGrp.replace("'",'') |
|---|
| 216 | SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','') #get rid of screw '_' |
|---|
| 217 | if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:] |
|---|
| 218 | SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 219 | magnetic = True |
|---|
| 220 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 221 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 222 | if not SpGrp: |
|---|
| 223 | print (MSpGrp) |
|---|
| 224 | self.warnings += 'No space group name was found in the CIF.' |
|---|
| 225 | return False |
|---|
| 226 | else: |
|---|
| 227 | SpGrp = SpGrp.replace('_','').split('(')[0] |
|---|
| 228 | SpGrp = G2spc.fullHM2shortHM(SpGrp) |
|---|
| 229 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 230 | #process space group symbol |
|---|
| 231 | E,SGData = G2spc.SpcGroup(SpGrp) |
|---|
| 232 | if E and SpGrp: |
|---|
| 233 | SpGrpNorm = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 234 | if SpGrpNorm: |
|---|
| 235 | E,SGData = G2spc.SpcGroup(SpGrpNorm) |
|---|
| 236 | # if E: #try lookup from number - found full symbol? |
|---|
| 237 | # SpGrpNorm = G2spc.spgbyNum[int(blk.get('_symmetry_Int_Tables_number'))] |
|---|
| 238 | # if SpGrpNorm: |
|---|
| 239 | # E,SGData = G2spc.SpcGroup(SpGrpNorm) |
|---|
| 240 | # nope, try the space group "out of the Box" |
|---|
| 241 | if E: |
|---|
| 242 | self.warnings += 'ERROR in space group symbol '+SpGrp |
|---|
| 243 | self.warnings += '\nThe space group has been set to "P 1". ' |
|---|
| 244 | self.warnings += "Change this in phase's General tab." |
|---|
| 245 | self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: ' |
|---|
| 246 | self.warnings += G2spc.SGErrors(E) |
|---|
| 247 | SGData = G2obj.P1SGData # P 1 |
|---|
| 248 | self.Phase['General']['SGData'] = SGData |
|---|
| 249 | # save symmetry operators, if specified (use to check for origin 1 in GSASIIdataGUI OnImportPhase) |
|---|
| 250 | ops = blk.get("_symmetry_equiv_pos_as_xyz") # try older name 1st |
|---|
| 251 | if ops: |
|---|
| 252 | self.SymOps['xyz'] = ops |
|---|
| 253 | elif blk.get("_space_group_symop_operation_xyz"): |
|---|
| 254 | self.SymOps['xyz'] = blk.get("_space_group_symop_operation_xyz") |
|---|
| 255 | else: |
|---|
| 256 | if 'xyz' in self.SymOps: del self.SymOps['xyz'] |
|---|
| 257 | if Super: |
|---|
| 258 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 259 | if E: |
|---|
| 260 | self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg |
|---|
| 261 | self.warnings += '\n'+E |
|---|
| 262 | SuperSg = SuperSg[:SuperSg.index(')')+1] |
|---|
| 263 | self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg |
|---|
| 264 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 265 | self.Phase['General']['SSGData'] = SSGData |
|---|
| 266 | if magnetic: |
|---|
| 267 | self.MPhase['General']['SGData'] = SGData |
|---|
| 268 | self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','') |
|---|
| 269 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 270 | |
|---|
| 271 | if magnetic: #replace std operaors with those from cif file - probably not the same! |
|---|
| 272 | SGData['SGFixed'] = True |
|---|
| 273 | SGData['SGOps'] = [] |
|---|
| 274 | SGData['SGCen'] = [] |
|---|
| 275 | if Super: |
|---|
| 276 | SSGData['SSGOps'] = [] |
|---|
| 277 | SSGData['SSGCen'] = [] |
|---|
| 278 | try: |
|---|
| 279 | sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id') |
|---|
| 280 | sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id') |
|---|
| 281 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1] |
|---|
| 282 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1] |
|---|
| 283 | except KeyError: #old mag cif names |
|---|
| 284 | sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id') |
|---|
| 285 | sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id') |
|---|
| 286 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1] |
|---|
| 287 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1] |
|---|
| 288 | spnflp = [] |
|---|
| 289 | for op in sgoploop: |
|---|
| 290 | M,T,S = G2spc.MagSSText2MTS(op[opid]) |
|---|
| 291 | SSGData['SSGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 292 | SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]]) |
|---|
| 293 | spnflp.append(S) |
|---|
| 294 | censpn = [] |
|---|
| 295 | for cent in sgcenloop: |
|---|
| 296 | M,C,S = G2spc.MagSSText2MTS(cent[centid]) |
|---|
| 297 | SSGData['SSGCen'].append(C) |
|---|
| 298 | SGData['SGCen'].append(C[:3]) |
|---|
| 299 | censpn += list(np.array(spnflp)*S) |
|---|
| 300 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 301 | else: |
|---|
| 302 | try: |
|---|
| 303 | sgoploop = blk.GetLoop('_space_group_symop_magn_operation.id') |
|---|
| 304 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1] |
|---|
| 305 | try: |
|---|
| 306 | sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id') |
|---|
| 307 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1] |
|---|
| 308 | except KeyError: |
|---|
| 309 | sgcenloop = None |
|---|
| 310 | except KeyError: |
|---|
| 311 | try: |
|---|
| 312 | sgoploop = blk.GetLoop('_space_group_symop_magn.id') |
|---|
| 313 | sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id') |
|---|
| 314 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1] |
|---|
| 315 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1] |
|---|
| 316 | except KeyError: #old mag cif names |
|---|
| 317 | sgoploop = blk.GetLoop('_space_group_symop.magn_id') |
|---|
| 318 | sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id') |
|---|
| 319 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1] |
|---|
| 320 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1] |
|---|
| 321 | spnflp = [] |
|---|
| 322 | for op in sgoploop: |
|---|
| 323 | try: |
|---|
| 324 | M,T,S = G2spc.MagText2MTS(op[opid]) |
|---|
| 325 | SGData['SGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 326 | spnflp.append(S) |
|---|
| 327 | except KeyError: |
|---|
| 328 | self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II' |
|---|
| 329 | return False |
|---|
| 330 | censpn = [] |
|---|
| 331 | if sgcenloop: |
|---|
| 332 | for cent in sgcenloop: |
|---|
| 333 | M,C,S = G2spc.MagText2MTS(cent[centid]) |
|---|
| 334 | SGData['SGCen'].append(C) |
|---|
| 335 | censpn += list(np.array(spnflp)*S) |
|---|
| 336 | else: |
|---|
| 337 | M,C,S = G2spc.MagText2MTS('x,y,z,+1') |
|---|
| 338 | SGData['SGCen'].append(C) |
|---|
| 339 | censpn += list(np.array(spnflp)*S) |
|---|
| 340 | self.MPhase['General']['SGData'] = SGData |
|---|
| 341 | self.MPhase['General']['SGData']['SpnFlp'] = censpn |
|---|
| 342 | G2spc.GenMagOps(SGData) #set magMom |
|---|
| 343 | self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp |
|---|
| 344 | if "1'" in MSpGrp: |
|---|
| 345 | SGData['SGGray'] = True |
|---|
| 346 | MagPtGp = blk.get('_space_group.magn_point_group') |
|---|
| 347 | if not MagPtGp: |
|---|
| 348 | MagPtGp = blk.get('_space_group.magn_point_group_name') |
|---|
| 349 | if not MagPtGp: |
|---|
| 350 | MagPtGp = blk.get('_space_group_magn.point_group_name') |
|---|
| 351 | self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp |
|---|
| 352 | |
|---|
| 353 | # cell parameters |
|---|
| 354 | cell = [] |
|---|
| 355 | for lbl in cellitems: |
|---|
| 356 | cell.append(cif.get_number_with_esd(blk[lbl])[0]) |
|---|
| 357 | Volume = G2lat.calc_V(G2lat.cell2A(cell)) |
|---|
| 358 | self.Phase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 359 | if magnetic: |
|---|
| 360 | self.MPhase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 361 | if Super: |
|---|
| 362 | waveloop = blk.GetLoop('_cell_wave_vector_seq_id') |
|---|
| 363 | waveDict = dict(waveloop.items()) |
|---|
| 364 | SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0], |
|---|
| 365 | cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0], |
|---|
| 366 | cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4] |
|---|
| 367 | |
|---|
| 368 | # read in atoms |
|---|
| 369 | self.errors = 'Error during reading of atoms' |
|---|
| 370 | atomlbllist = [] # table to look up atom IDs |
|---|
| 371 | atomloop = blk.GetLoop('_atom_site_label') |
|---|
| 372 | atomkeys = [i.lower() for i in atomloop.keys()] |
|---|
| 373 | if not blk.get('_atom_site_type_symbol'): |
|---|
| 374 | isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed' |
|---|
| 375 | if magnetic: |
|---|
| 376 | try: |
|---|
| 377 | magmoment = '_atom_site_moment.label' |
|---|
| 378 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 379 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 380 | magatomlabels = blk.get(magmoment) |
|---|
| 381 | G2MagDict = {'_atom_site_moment.label': 0, |
|---|
| 382 | '_atom_site_moment.crystalaxis_x':7, |
|---|
| 383 | '_atom_site_moment.crystalaxis_y':8, |
|---|
| 384 | '_atom_site_moment.crystalaxis_z':9} |
|---|
| 385 | except KeyError: |
|---|
| 386 | magmoment = '_atom_site_moment_label' |
|---|
| 387 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 388 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 389 | magatomlabels = blk.get(magmoment) |
|---|
| 390 | G2MagDict = {'_atom_site_moment_label': 0, |
|---|
| 391 | '_atom_site_moment_crystalaxis_x':7, |
|---|
| 392 | '_atom_site_moment_crystalaxis_y':8, |
|---|
| 393 | '_atom_site_moment_crystalaxis_z':9} |
|---|
| 394 | |
|---|
| 395 | if blk.get('_atom_site_aniso_label'): |
|---|
| 396 | anisoloop = blk.GetLoop('_atom_site_aniso_label') |
|---|
| 397 | anisokeys = [i.lower() for i in anisoloop.keys()] |
|---|
| 398 | anisolabels = blk.get('_atom_site_aniso_label') |
|---|
| 399 | else: |
|---|
| 400 | anisoloop = None |
|---|
| 401 | anisokeys = [] |
|---|
| 402 | anisolabels = [] |
|---|
| 403 | if Super: |
|---|
| 404 | # occFloop = None |
|---|
| 405 | occCloop = None |
|---|
| 406 | # occFdict = {} |
|---|
| 407 | occCdict = {} |
|---|
| 408 | # displSloop = None |
|---|
| 409 | displFloop = None |
|---|
| 410 | MagFloop = None |
|---|
| 411 | # displSdict = {} |
|---|
| 412 | displFdict = {} |
|---|
| 413 | MagFdict = {} |
|---|
| 414 | UijFloop = None |
|---|
| 415 | UijFdict = {} |
|---|
| 416 | #occupancy modulation |
|---|
| 417 | # if blk.get('_atom_site_occ_Fourier_atom_site_label'): |
|---|
| 418 | # occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label') |
|---|
| 419 | # occFdict = dict(occFloop.items()) |
|---|
| 420 | if blk.get('_atom_site_occ_special_func_atom_site_label'): #Crenel (i.e. Block Wave) occ |
|---|
| 421 | occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label') |
|---|
| 422 | occCdict = dict(occCloop.items()) |
|---|
| 423 | #position modulation |
|---|
| 424 | if blk.get('_atom_site_displace_Fourier_atom_site_label'): |
|---|
| 425 | displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label') |
|---|
| 426 | displFdict = dict(displFloop.items()) |
|---|
| 427 | # if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth |
|---|
| 428 | # displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label') |
|---|
| 429 | # displSdict = dict(displSloop.items()) |
|---|
| 430 | #U modulation |
|---|
| 431 | if blk.get('_atom_site_U_Fourier_atom_site_label'): |
|---|
| 432 | UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label') |
|---|
| 433 | UijFdict = dict(UijFloop.items()) |
|---|
| 434 | #Mag moment modulation |
|---|
| 435 | if blk.get('_atom_site_moment_Fourier_atom_site_label'): |
|---|
| 436 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label') |
|---|
| 437 | MagFdict = dict(MagFloop.items()) |
|---|
| 438 | Mnames = ['_atom_site_moment_fourier_atom_site_label', |
|---|
| 439 | '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id', |
|---|
| 440 | '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos'] |
|---|
| 441 | elif blk.get('_atom_site_moment_Fourier.atom_site_label'): |
|---|
| 442 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label') |
|---|
| 443 | MagFdict = dict(MagFloop.items()) |
|---|
| 444 | Mnames = ['_atom_site_moment_fourier.atom_site_label', |
|---|
| 445 | '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id', |
|---|
| 446 | '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos'] |
|---|
| 447 | self.Phase['Atoms'] = [] |
|---|
| 448 | if magnetic: |
|---|
| 449 | self.MPhase['Atoms'] = [] |
|---|
| 450 | G2AtomDict = { '_atom_site_type_symbol' : 1, |
|---|
| 451 | '_atom_site_label' : 0, |
|---|
| 452 | '_atom_site_fract_x' : 3, |
|---|
| 453 | '_atom_site_fract_y' : 4, |
|---|
| 454 | '_atom_site_fract_z' : 5, |
|---|
| 455 | '_atom_site_occupancy' : 6, |
|---|
| 456 | '_atom_site_aniso_u_11' : 11, |
|---|
| 457 | '_atom_site_aniso_u_22' : 12, |
|---|
| 458 | '_atom_site_aniso_u_33' : 13, |
|---|
| 459 | '_atom_site_aniso_u_12' : 14, |
|---|
| 460 | '_atom_site_aniso_u_13' : 15, |
|---|
| 461 | '_atom_site_aniso_u_23' : 16, } |
|---|
| 462 | |
|---|
| 463 | ranIdlookup = {} |
|---|
| 464 | for aitem in atomloop: |
|---|
| 465 | atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,] |
|---|
| 466 | for val,key in zip(aitem,atomkeys): |
|---|
| 467 | col = G2AtomDict.get(key,-1) |
|---|
| 468 | if col >= 3: |
|---|
| 469 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 470 | if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag |
|---|
| 471 | elif col is not None and col != -1: |
|---|
| 472 | atomlist[col] = val |
|---|
| 473 | elif key in ('_atom_site_thermal_displace_type', |
|---|
| 474 | '_atom_site_adp_type'): #Iso or Aniso? |
|---|
| 475 | if val.lower() == 'uani': |
|---|
| 476 | atomlist[9] = 'A' |
|---|
| 477 | elif key == '_atom_site_u_iso_or_equiv': |
|---|
| 478 | uisoval = cif.get_number_with_esd(val)[0] |
|---|
| 479 | if uisoval is not None: |
|---|
| 480 | atomlist[10] = uisoval |
|---|
| 481 | elif key == '_atom_site_b_iso_or_equiv': |
|---|
| 482 | uisoval = cif.get_number_with_esd(val)[0] |
|---|
| 483 | if uisoval is not None: |
|---|
| 484 | atomlist[10] = uisoval/(8*np.pi**2) |
|---|
| 485 | if not atomlist[1] and atomlist[0]: |
|---|
| 486 | typ = atomlist[0].rstrip('0123456789-+') |
|---|
| 487 | if G2elem.CheckElement(typ): |
|---|
| 488 | atomlist[1] = typ |
|---|
| 489 | if not atomlist[1]: |
|---|
| 490 | atomlist[1] = 'Xe' |
|---|
| 491 | self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n' |
|---|
| 492 | if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop? |
|---|
| 493 | atomlist[9] = 'A' |
|---|
| 494 | for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys): |
|---|
| 495 | col = G2AtomDict.get(key) |
|---|
| 496 | if col: |
|---|
| 497 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 498 | if None in atomlist[11:17]: |
|---|
| 499 | atomlist[9] = 'I' |
|---|
| 500 | atomlist[11:17] = [0.,0.,0.,0.,0.,0.] |
|---|
| 501 | atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2] |
|---|
| 502 | atomlist[1] = G2elem.FixValence(atomlist[1]) |
|---|
| 503 | atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id |
|---|
| 504 | self.Phase['Atoms'].append(atomlist) |
|---|
| 505 | ranIdlookup[atomlist[0]] = atomlist[-1] |
|---|
| 506 | if atomlist[0] in atomlbllist: |
|---|
| 507 | self.warnings += ' ERROR: repeated atom label: '+atomlist[0] |
|---|
| 508 | else: |
|---|
| 509 | atomlbllist.append(atomlist[0]) |
|---|
| 510 | |
|---|
| 511 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 512 | matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:] |
|---|
| 513 | for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys): |
|---|
| 514 | mcol = G2MagDict.get(mkey,-1) |
|---|
| 515 | if mcol > 0: |
|---|
| 516 | matomlist[mcol] = cif.get_number_with_esd(mval)[0] |
|---|
| 517 | self.MPhase['Atoms'].append(matomlist) |
|---|
| 518 | if Super: |
|---|
| 519 | Sfrac = np.zeros((4,2)) |
|---|
| 520 | Sadp = np.zeros((4,12)) |
|---|
| 521 | Spos = np.zeros((4,6)) |
|---|
| 522 | Smag = np.zeros((4,6)) |
|---|
| 523 | nim = -1 |
|---|
| 524 | waveType = 'Fourier' |
|---|
| 525 | if occCdict: |
|---|
| 526 | for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']): |
|---|
| 527 | if item == atomlist[0]: |
|---|
| 528 | waveType = 'Crenel' |
|---|
| 529 | val = occCdict['_atom_site_occ_special_func_crenel_c'][i] |
|---|
| 530 | Sfrac[0][0] = cif.get_number_with_esd(val)[0] |
|---|
| 531 | val = occCdict['_atom_site_occ_special_func_crenel_w'][i] |
|---|
| 532 | Sfrac[0][1] = cif.get_number_with_esd(val)[0] |
|---|
| 533 | nim = 1 |
|---|
| 534 | |
|---|
| 535 | if nim >= 0: |
|---|
| 536 | Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]] |
|---|
| 537 | else: |
|---|
| 538 | Sfrac = [] |
|---|
| 539 | nim = -1 |
|---|
| 540 | if displFdict: |
|---|
| 541 | for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']): |
|---|
| 542 | if item == atomlist[0]: |
|---|
| 543 | waveType = 'Fourier' |
|---|
| 544 | ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i]) |
|---|
| 545 | im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i]) |
|---|
| 546 | if im != nim: |
|---|
| 547 | nim = im |
|---|
| 548 | val = displFdict['_atom_site_displace_fourier_param_sin'][i] |
|---|
| 549 | Spos[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 550 | val = displFdict['_atom_site_displace_fourier_param_cos'][i] |
|---|
| 551 | Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 552 | if nim >= 0: |
|---|
| 553 | Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]] |
|---|
| 554 | else: |
|---|
| 555 | Spos = [] |
|---|
| 556 | nim = -1 |
|---|
| 557 | if UijFdict: |
|---|
| 558 | nim = -1 |
|---|
| 559 | for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']): |
|---|
| 560 | if item == atomlist[0]: |
|---|
| 561 | ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i]) |
|---|
| 562 | im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i]) |
|---|
| 563 | if im != nim: |
|---|
| 564 | nim = im |
|---|
| 565 | val = UijFdict['_atom_site_u_fourier_param_sin'][i] |
|---|
| 566 | Sadp[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 567 | val = UijFdict['_atom_site_u_fourier_param_cos'][i] |
|---|
| 568 | Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0] |
|---|
| 569 | if nim >= 0: |
|---|
| 570 | Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]] |
|---|
| 571 | else: |
|---|
| 572 | Sadp = [] |
|---|
| 573 | nim = -1 |
|---|
| 574 | if MagFdict: |
|---|
| 575 | nim = -1 |
|---|
| 576 | for i,item in enumerate(MagFdict[Mnames[0]]): |
|---|
| 577 | if item == atomlist[0]: |
|---|
| 578 | ix = ['x','y','z'].index(MagFdict[Mnames[1]][i]) |
|---|
| 579 | im = int(MagFdict[Mnames[2]][i]) |
|---|
| 580 | if im != nim: |
|---|
| 581 | nim = im |
|---|
| 582 | val = MagFdict[Mnames[3]][i] |
|---|
| 583 | Smag[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 584 | val = MagFdict[Mnames[4]][i] |
|---|
| 585 | Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 586 | if nim >= 0: |
|---|
| 587 | Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]] |
|---|
| 588 | else: |
|---|
| 589 | Smag = [] |
|---|
| 590 | SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}} |
|---|
| 591 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 592 | matomlist.append(SSdict) |
|---|
| 593 | atomlist.append(SSdict) |
|---|
| 594 | if len(atomlbllist) != len(self.Phase['Atoms']): |
|---|
| 595 | isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode' |
|---|
| 596 | for lbl in phasenamefields: # get a name for the phase |
|---|
| 597 | try: |
|---|
| 598 | name = blk.get(lbl) |
|---|
| 599 | except: |
|---|
| 600 | continue |
|---|
| 601 | if name is None: |
|---|
| 602 | continue |
|---|
| 603 | name = name.strip() |
|---|
| 604 | if name == '?' or name == '.': |
|---|
| 605 | continue |
|---|
| 606 | else: |
|---|
| 607 | break |
|---|
| 608 | else: # no name found, use block name for lack of a better choice; for isodistort use filename |
|---|
| 609 | name = blknm |
|---|
| 610 | if 'isodistort' in name: |
|---|
| 611 | self.Phase['General']['Name'] = os.path.split(filename)[1].split('.')[0] |
|---|
| 612 | else: |
|---|
| 613 | self.Phase['General']['Name'] = name.strip() |
|---|
| 614 | self.Phase['General']['Super'] = Super |
|---|
| 615 | self.Phase = copy.deepcopy(self.Phase) #clean copy |
|---|
| 616 | if magnetic: |
|---|
| 617 | self.MPhase = copy.deepcopy(self.MPhase) #clean copy |
|---|
| 618 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 619 | self.MPhase['General']['Name'] = name.strip()+' mag' |
|---|
| 620 | self.MPhase['General']['Super'] = Super |
|---|
| 621 | if Super: |
|---|
| 622 | self.MPhase['General']['Modulated'] = True |
|---|
| 623 | self.MPhase['General']['SuperVec'] = SuperVec |
|---|
| 624 | self.MPhase['General']['SuperSg'] = SuperSg |
|---|
| 625 | if self.MPhase['General']['SGData']['SGGray']: |
|---|
| 626 | self.MPhase['General']['SuperSg'] += 's' |
|---|
| 627 | if 'mcif' not in filename: |
|---|
| 628 | self.Phase = copy.deepcopy(self.MPhase) |
|---|
| 629 | del self.MPhase |
|---|
| 630 | else: |
|---|
| 631 | self.MPhase = None |
|---|
| 632 | if Super: |
|---|
| 633 | if self.Phase['General']['SGData']['SGGray']: |
|---|
| 634 | SGData = self.Phase['General']['SGData'] |
|---|
| 635 | SGData['SGGray'] = False |
|---|
| 636 | ncen = len(SGData['SGCen']) |
|---|
| 637 | SGData['SGCen'] = SGData['SGCen'][:ncen//2] |
|---|
| 638 | self.Phase['General']['SGData'].update(SGData) |
|---|
| 639 | self.Phase['General']['Modulated'] = True |
|---|
| 640 | self.Phase['General']['SuperVec'] = SuperVec |
|---|
| 641 | self.Phase['General']['SuperSg'] = SuperSg |
|---|
| 642 | if not self.Phase['General']['SGData']['SGFixed']: |
|---|
| 643 | self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1] |
|---|
| 644 | if self.ISODISTORT_test(blk): |
|---|
| 645 | if isodistort_warnings: |
|---|
| 646 | self.warnings += isodistort_warnings |
|---|
| 647 | else: |
|---|
| 648 | self.errors = "Error while processing ISODISTORT constraints" |
|---|
| 649 | self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup,filename) |
|---|
| 650 | self.errors = "" |
|---|
| 651 | returnstat = True |
|---|
| 652 | return returnstat |
|---|
| 653 | |
|---|
| 654 | def ISODISTORT_test(self,blk): |
|---|
| 655 | '''Test if there is any ISODISTORT information in CIF |
|---|
| 656 | |
|---|
| 657 | At present only _iso_displacivemode... and _iso_occupancymode... are |
|---|
| 658 | tested. |
|---|
| 659 | ''' |
|---|
| 660 | for i in ('_iso_displacivemode_label', |
|---|
| 661 | '_iso_occupancymode_label'): |
|---|
| 662 | if blk.get(i): return True |
|---|
| 663 | return False |
|---|
| 664 | |
|---|
| 665 | def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup,filename): |
|---|
| 666 | '''Process ISODISTORT items to create constraints etc. |
|---|
| 667 | Constraints are generated from information extracted from |
|---|
| 668 | loops beginning with _iso_ and are placed into |
|---|
| 669 | self.Constraints, which contains a list of |
|---|
| 670 | :ref:`constraints tree items <Constraint_definitions_table>` |
|---|
| 671 | and one dict. |
|---|
| 672 | The dict contains help text for each generated ISODISTORT variable |
|---|
| 673 | |
|---|
| 674 | At present only _iso_displacivemode... and _iso_occupancymode... are |
|---|
| 675 | processed. Not yet processed: _iso_magneticmode..., |
|---|
| 676 | _iso_rotationalmode... & _iso_strainmode... |
|---|
| 677 | ''' |
|---|
| 678 | varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'} |
|---|
| 679 | 'Maps ISODISTORT parm names to GSAS-II names' |
|---|
| 680 | # used for all types of modes |
|---|
| 681 | self.Constraints = [] |
|---|
| 682 | explaination = {} |
|---|
| 683 | G2obj.AddPhase2Index(self,filename) # put phase info into Var index |
|---|
| 684 | #---------------------------------------------------------------------- |
|---|
| 685 | # read in the ISODISTORT displacement modes |
|---|
| 686 | #---------------------------------------------------------------------- |
|---|
| 687 | if blk.get('_iso_displacivemode_label'): |
|---|
| 688 | modelist = [] |
|---|
| 689 | shortmodelist = [] |
|---|
| 690 | modedispl = [] |
|---|
| 691 | idlist = [] |
|---|
| 692 | for id,lbl,val in zip( |
|---|
| 693 | blk.get('_iso_displacivemode_ID'), |
|---|
| 694 | blk.get('_iso_displacivemode_label'), |
|---|
| 695 | blk.get('_iso_displacivemode_value')): |
|---|
| 696 | idlist.append(int(id)) |
|---|
| 697 | modelist.append(lbl) |
|---|
| 698 | modedispl.append(float(val)) |
|---|
| 699 | ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique |
|---|
| 700 | # just in case the items are not ordered increasing by id, sort them here |
|---|
| 701 | modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])] |
|---|
| 702 | shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])] |
|---|
| 703 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 704 | coordVarLbl = [] |
|---|
| 705 | G2varObj = [] |
|---|
| 706 | G2coordOffset = [] |
|---|
| 707 | idlist = [] |
|---|
| 708 | error = False |
|---|
| 709 | ParentCoordinates = {} |
|---|
| 710 | for lbl,exp in zip( |
|---|
| 711 | blk.get('_iso_coordinate_label'), |
|---|
| 712 | blk.get('_iso_coordinate_formula') ): |
|---|
| 713 | if '_' in lbl: |
|---|
| 714 | albl = lbl[:lbl.rfind('_')] |
|---|
| 715 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 716 | else: |
|---|
| 717 | self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl |
|---|
| 718 | error = True |
|---|
| 719 | continue |
|---|
| 720 | if vlbl not in 'xyz' or len(vlbl) != 1: |
|---|
| 721 | self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl |
|---|
| 722 | error = True |
|---|
| 723 | continue |
|---|
| 724 | i = 'xyz'.index(vlbl) |
|---|
| 725 | if not ParentCoordinates.get(albl): |
|---|
| 726 | ParentCoordinates[albl] = [None,None,None] |
|---|
| 727 | if '+' in exp: |
|---|
| 728 | val = exp.split('+')[0].strip() |
|---|
| 729 | val = G2fil.FormulaEval(val) |
|---|
| 730 | elif '-' in exp: |
|---|
| 731 | val = exp.split('-')[0].strip() |
|---|
| 732 | val = G2fil.FormulaEval(val) |
|---|
| 733 | else: |
|---|
| 734 | val = G2fil.FormulaEval(exp) |
|---|
| 735 | if val is None: |
|---|
| 736 | self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl |
|---|
| 737 | error = True |
|---|
| 738 | continue |
|---|
| 739 | else: |
|---|
| 740 | ParentCoordinates[albl][i] = val |
|---|
| 741 | if error: |
|---|
| 742 | print (self.warnings) |
|---|
| 743 | raise Exception("Error decoding variable labels") |
|---|
| 744 | |
|---|
| 745 | error = False |
|---|
| 746 | coordOffset = {xyz:i for i,xyz in enumerate(('dx','dy','dz'))} |
|---|
| 747 | for id,lbl,val in zip( |
|---|
| 748 | blk.get('_iso_deltacoordinate_ID'), |
|---|
| 749 | blk.get('_iso_deltacoordinate_label'), |
|---|
| 750 | blk.get('_iso_deltacoordinate_value') ): |
|---|
| 751 | idlist.append(int(id)) |
|---|
| 752 | coordVarLbl.append(lbl) |
|---|
| 753 | if '_' in lbl: |
|---|
| 754 | albl = lbl[:lbl.rfind('_')] |
|---|
| 755 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 756 | else: |
|---|
| 757 | self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl |
|---|
| 758 | error = True |
|---|
| 759 | continue |
|---|
| 760 | if albl not in atomlbllist: |
|---|
| 761 | self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl |
|---|
| 762 | error = True |
|---|
| 763 | continue |
|---|
| 764 | # else: |
|---|
| 765 | # anum = atomlbllist.index(albl) |
|---|
| 766 | var = varLookup.get(vlbl) |
|---|
| 767 | if not var: |
|---|
| 768 | self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl |
|---|
| 769 | error = True |
|---|
| 770 | continue |
|---|
| 771 | G2coordOffset.append(ParentCoordinates[albl][coordOffset[vlbl]] - float(val)) |
|---|
| 772 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 773 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 774 | )) |
|---|
| 775 | if error: |
|---|
| 776 | raise Exception("Error decoding variable labels") |
|---|
| 777 | # just in case the items are not ordered increasing by id, sort them here |
|---|
| 778 | coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])] |
|---|
| 779 | G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])] |
|---|
| 780 | |
|---|
| 781 | if len(G2varObj) != len(modelist): |
|---|
| 782 | print ("non-square input") |
|---|
| 783 | raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate") |
|---|
| 784 | |
|---|
| 785 | # normalization constants |
|---|
| 786 | normlist = [] |
|---|
| 787 | idlist = [] |
|---|
| 788 | for id,exp in zip( |
|---|
| 789 | blk.get('_iso_displacivemodenorm_ID'), |
|---|
| 790 | blk.get('_iso_displacivemodenorm_value'), |
|---|
| 791 | ): |
|---|
| 792 | idlist.append(int(id)) |
|---|
| 793 | normlist.append(float(exp)) |
|---|
| 794 | normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])] |
|---|
| 795 | # get mapping of modes to atomic coordinate displacements |
|---|
| 796 | displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 797 | for row,col,val in zip( |
|---|
| 798 | blk.get('_iso_displacivemodematrix_row'), |
|---|
| 799 | blk.get('_iso_displacivemodematrix_col'), |
|---|
| 800 | blk.get('_iso_displacivemodematrix_value'),): |
|---|
| 801 | displacivemodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 802 | # Invert to get mapping of atom displacements to modes |
|---|
| 803 | Var2ModeMatrix = np.linalg.inv(displacivemodematrix) |
|---|
| 804 | # create the constraints |
|---|
| 805 | modeVarList = [] |
|---|
| 806 | for i,(row,norm) in enumerate(zip(Var2ModeMatrix,normlist)): |
|---|
| 807 | constraint = [] |
|---|
| 808 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 809 | if k == 0: continue |
|---|
| 810 | constraint.append([k/norm,G2varObj[j]]) |
|---|
| 811 | modeVar = G2obj.G2VarObj( |
|---|
| 812 | (self.Phase['ranId'],None,shortmodelist[i],None)) |
|---|
| 813 | modeVarList.append(modeVar) |
|---|
| 814 | constraint += [modeVar,False,'f'] |
|---|
| 815 | self.Constraints.append(constraint) |
|---|
| 816 | #---------------------------------------------------------------------- |
|---|
| 817 | # save the ISODISTORT info for "mode analysis" |
|---|
| 818 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 819 | self.Phase['ISODISTORT'].update({ |
|---|
| 820 | # coordinate items |
|---|
| 821 | 'IsoVarList' : coordVarLbl, |
|---|
| 822 | 'G2VarList' : G2varObj, |
|---|
| 823 | 'G2coordOffset' : G2coordOffset, |
|---|
| 824 | 'G2parentCoords' : [ParentCoordinates[item] for item in ParentCoordinates], #Assumes python 3.7 dict ordering! |
|---|
| 825 | # mode items |
|---|
| 826 | 'IsoModeList' : modelist, |
|---|
| 827 | 'G2ModeList' : modeVarList, |
|---|
| 828 | 'NormList' : normlist, |
|---|
| 829 | 'modeDispl' : modedispl, |
|---|
| 830 | 'ISOmodeDispl' : copy.deepcopy(modedispl), |
|---|
| 831 | # transform matrices |
|---|
| 832 | 'Var2ModeMatrix' : Var2ModeMatrix, |
|---|
| 833 | 'Mode2VarMatrix' : displacivemodematrix, |
|---|
| 834 | }) |
|---|
| 835 | # make explaination dictionary |
|---|
| 836 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 837 | modeVar = G2obj.G2VarObj( |
|---|
| 838 | (self.Phase['ranId'],None,shortmode,None)) |
|---|
| 839 | explaination[modeVar] = ("Full ISODISTORT name for " + |
|---|
| 840 | shortmode + " is " + str(mode)) |
|---|
| 841 | #---------------------------------------------------------------------- |
|---|
| 842 | # now read in the ISODISTORT occupancy modes |
|---|
| 843 | #---------------------------------------------------------------------- |
|---|
| 844 | if blk.get('_iso_occupancymode_label'): |
|---|
| 845 | modelist = [] |
|---|
| 846 | shortmodelist = [] |
|---|
| 847 | idlist = [] |
|---|
| 848 | for id,lbl in zip( |
|---|
| 849 | blk.get('_iso_occupancymode_ID'), |
|---|
| 850 | blk.get('_iso_occupancymode_label')): |
|---|
| 851 | idlist.append(int(id)) |
|---|
| 852 | modelist.append(lbl) |
|---|
| 853 | ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique |
|---|
| 854 | # just in case the items are not ordered increasing by id, sort them here |
|---|
| 855 | modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])] |
|---|
| 856 | shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])] |
|---|
| 857 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 858 | occVarLbl = [] |
|---|
| 859 | G2varObj = [] |
|---|
| 860 | idlist = [] |
|---|
| 861 | error = False |
|---|
| 862 | for id,lbl in zip( |
|---|
| 863 | blk.get('_iso_deltaoccupancy_ID'), |
|---|
| 864 | blk.get('_iso_deltaoccupancy_label') ): |
|---|
| 865 | idlist.append(int(id)) |
|---|
| 866 | occVarLbl.append(lbl) |
|---|
| 867 | if '_' in lbl: |
|---|
| 868 | albl = lbl[:lbl.rfind('_')] |
|---|
| 869 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 870 | else: |
|---|
| 871 | self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl |
|---|
| 872 | error = True |
|---|
| 873 | continue |
|---|
| 874 | if albl not in atomlbllist: |
|---|
| 875 | self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl |
|---|
| 876 | error = True |
|---|
| 877 | continue |
|---|
| 878 | # else: |
|---|
| 879 | # anum = atomlbllist.index(albl) |
|---|
| 880 | var = varLookup.get(vlbl) |
|---|
| 881 | if not var: |
|---|
| 882 | self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl |
|---|
| 883 | error = True |
|---|
| 884 | continue |
|---|
| 885 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 886 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 887 | )) |
|---|
| 888 | if error: |
|---|
| 889 | raise Exception("Error decoding variable labels") |
|---|
| 890 | # just in case the items are not ordered increasing by id, sort them here |
|---|
| 891 | occVarLbl = [i for i,j in sorted(zip(occVarLbl,idlist),key=lambda k:k[1])] |
|---|
| 892 | G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])] |
|---|
| 893 | |
|---|
| 894 | if len(G2varObj) != len(modelist): |
|---|
| 895 | print ("non-square input") |
|---|
| 896 | raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy") |
|---|
| 897 | |
|---|
| 898 | error = False |
|---|
| 899 | ParentOcc = {} |
|---|
| 900 | for lbl,exp in zip( |
|---|
| 901 | blk.get('_iso_occupancy_label'), |
|---|
| 902 | blk.get('_iso_occupancy_formula') ): |
|---|
| 903 | if '_' in lbl: |
|---|
| 904 | albl = lbl[:lbl.rfind('_')] |
|---|
| 905 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 906 | else: |
|---|
| 907 | self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl |
|---|
| 908 | error = True |
|---|
| 909 | continue |
|---|
| 910 | if vlbl != 'occ': |
|---|
| 911 | self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl |
|---|
| 912 | error = True |
|---|
| 913 | continue |
|---|
| 914 | if '+' in exp: |
|---|
| 915 | val = exp.split('+')[0].strip() |
|---|
| 916 | val = G2fil.FormulaEval(val) |
|---|
| 917 | if val is None: |
|---|
| 918 | self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl |
|---|
| 919 | error = True |
|---|
| 920 | continue |
|---|
| 921 | ParentOcc[albl] = val |
|---|
| 922 | if error: |
|---|
| 923 | raise Exception("Error decoding occupancy labels") |
|---|
| 924 | # get mapping of modes to atomic coordinate displacements |
|---|
| 925 | occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 926 | for row,col,val in zip( |
|---|
| 927 | blk.get('_iso_occupancymodematrix_row'), |
|---|
| 928 | blk.get('_iso_occupancymodematrix_col'), |
|---|
| 929 | blk.get('_iso_occupancymodematrix_value'),): |
|---|
| 930 | occupancymodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 931 | # Invert to get mapping of atom displacements to modes |
|---|
| 932 | occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix) |
|---|
| 933 | # create the constraints |
|---|
| 934 | modeVarList = [] |
|---|
| 935 | for i,row in enumerate(occupancymodeInvmatrix): |
|---|
| 936 | constraint = [] |
|---|
| 937 | for j,(lbl,k) in enumerate(zip(occVarLbl,row)): |
|---|
| 938 | if k == 0: continue |
|---|
| 939 | constraint.append([k,G2varObj[j]]) |
|---|
| 940 | modeVar = G2obj.G2VarObj( |
|---|
| 941 | (self.Phase['ranId'],None,shortmodelist[i],None)) |
|---|
| 942 | modeVarList.append(modeVar) |
|---|
| 943 | constraint += [modeVar,False,'f'] |
|---|
| 944 | self.Constraints.append(constraint) |
|---|
| 945 | # normilization constants |
|---|
| 946 | normlist = [] |
|---|
| 947 | idlist = [] |
|---|
| 948 | for id,exp in zip( |
|---|
| 949 | blk.get('_iso_occupancymodenorm_ID'), |
|---|
| 950 | blk.get('_iso_occupancymodenorm_value'), |
|---|
| 951 | ): |
|---|
| 952 | idlist.append(int(id)) |
|---|
| 953 | normlist.append(float(exp)) |
|---|
| 954 | normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])] |
|---|
| 955 | #---------------------------------------------------------------------- |
|---|
| 956 | # save the ISODISTORT info for "mode analysis" |
|---|
| 957 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 958 | self.Phase['ISODISTORT'].update({ |
|---|
| 959 | # coordinate items |
|---|
| 960 | 'OccVarList' : occVarLbl, |
|---|
| 961 | 'G2OccVarList' : G2varObj, |
|---|
| 962 | 'BaseOcc' : ParentOcc, |
|---|
| 963 | # mode items |
|---|
| 964 | 'OccModeList' : modelist, |
|---|
| 965 | 'G2OccModeList' : modeVarList, |
|---|
| 966 | 'OccNormList' : normlist, |
|---|
| 967 | # transform matrices |
|---|
| 968 | 'Var2OccMatrix' : occupancymodeInvmatrix, |
|---|
| 969 | 'Occ2VarMatrix' : occupancymodematrix, |
|---|
| 970 | }) |
|---|
| 971 | # make explaination dictionary |
|---|
| 972 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 973 | modeVar = G2obj.G2VarObj( |
|---|
| 974 | (self.Phase['ranId'],None,shortmode,None)) |
|---|
| 975 | explaination[modeVar] = ("Full ISODISTORT name for " + |
|---|
| 976 | shortmode + " is " + str(mode)) |
|---|
| 977 | if explaination: self.Constraints.append(explaination) |
|---|
| 978 | #---------------------------------------------------------------------- |
|---|
| 979 | # done with read |
|---|
| 980 | #---------------------------------------------------------------------- |
|---|
| 981 | |
|---|
| 982 | def fmtEqn(i,head,l,var,k): |
|---|
| 983 | 'format a section of a row of variables and multipliers' |
|---|
| 984 | if np.isclose(k,0): return head,l |
|---|
| 985 | if len(head) + len(l) > 65: |
|---|
| 986 | print(head+l) |
|---|
| 987 | head = 20*' ' |
|---|
| 988 | l = '' |
|---|
| 989 | if k < 0 and i > 0: |
|---|
| 990 | l += ' - ' |
|---|
| 991 | k = -k |
|---|
| 992 | elif i > 0: |
|---|
| 993 | l += ' + ' |
|---|
| 994 | if k == 1: |
|---|
| 995 | l += '%s ' % str(var) |
|---|
| 996 | else: |
|---|
| 997 | l += '%.3f * %s' % (k,str(var)) |
|---|
| 998 | return head,l |
|---|
| 999 | |
|---|
| 1000 | # debug: show displacive mode var to mode relations |
|---|
| 1001 | if debug and 'IsoVarList' in self.Phase['ISODISTORT']: |
|---|
| 1002 | print('\n' + 70*'=') |
|---|
| 1003 | print('ISO modes from Iso coordinate vars (using Var2ModeMatrix, IsoVarList, G2VarList & G2ModeList)' ) |
|---|
| 1004 | for i,row in enumerate(self.Phase['ISODISTORT']['Var2ModeMatrix']): |
|---|
| 1005 | norm = self.Phase['ISODISTORT']['NormList'][i] |
|---|
| 1006 | head = ' ' + str(self.Phase['ISODISTORT']['G2ModeList'][i]) + ' = (' |
|---|
| 1007 | line = '' |
|---|
| 1008 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 1009 | var = self.Phase['ISODISTORT']['IsoVarList'][j] |
|---|
| 1010 | head,line = fmtEqn(j,head,line,var,k) |
|---|
| 1011 | print(head+line+') / {:.3g}'.format(norm)) |
|---|
| 1012 | head = ' = ' |
|---|
| 1013 | line = '' |
|---|
| 1014 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 1015 | var = self.Phase['ISODISTORT']['IsoVarList'][j] |
|---|
| 1016 | head,line = fmtEqn(j,head,line,var,k/norm) |
|---|
| 1017 | print(head+line) |
|---|
| 1018 | print('\nConstraints') |
|---|
| 1019 | for c in self.Constraints: |
|---|
| 1020 | if type(c) is dict: continue |
|---|
| 1021 | if c[-1] != 'f': continue |
|---|
| 1022 | line = '' |
|---|
| 1023 | head = ' ' + str(c[-3]) + ' = ' |
|---|
| 1024 | for j,(k,var) in enumerate(c[:-3]): |
|---|
| 1025 | head,line = fmtEqn(j,head,line,var,k) |
|---|
| 1026 | print(head+line) |
|---|
| 1027 | |
|---|
| 1028 | # Get the ISODISTORT offset values |
|---|
| 1029 | coordVarDelta = {} |
|---|
| 1030 | for lbl,val in zip( |
|---|
| 1031 | blk.get('_iso_deltacoordinate_label'), |
|---|
| 1032 | blk.get('_iso_deltacoordinate_value'),): |
|---|
| 1033 | coordVarDelta[lbl] = float(val) |
|---|
| 1034 | modeVarDelta = {} |
|---|
| 1035 | for lbl,val in zip( |
|---|
| 1036 | blk.get('_iso_displacivemode_label'), |
|---|
| 1037 | blk.get('_iso_displacivemode_value'),): |
|---|
| 1038 | modeVarDelta[lbl] = cif.get_number_with_esd(val)[0] |
|---|
| 1039 | |
|---|
| 1040 | print('\n' + 70*'=') |
|---|
| 1041 | print('Confirming inverse mode relations computed from displacement values', |
|---|
| 1042 | '\nusing Var2ModeMatrix, NormList, IsoVarList') |
|---|
| 1043 | # compute the mode values from the reported coordinate deltas |
|---|
| 1044 | for i,(row,n) in enumerate(zip(self.Phase['ISODISTORT']['Var2ModeMatrix'], |
|---|
| 1045 | self.Phase['ISODISTORT']['NormList'])): |
|---|
| 1046 | line = '' |
|---|
| 1047 | print(str(self.Phase['ISODISTORT']['IsoModeList'][i])+' = ') |
|---|
| 1048 | head = ' = (' |
|---|
| 1049 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 1050 | head,line = fmtEqn(j,head,line,lbl,k) |
|---|
| 1051 | print(head+line+') / '+('%.3f'%n)) |
|---|
| 1052 | line = '' |
|---|
| 1053 | head = ' = (' |
|---|
| 1054 | vsum = 0. |
|---|
| 1055 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 1056 | val = "{:3g}".format(coordVarDelta[lbl]) |
|---|
| 1057 | head,line = fmtEqn(j,head,line,val,k) |
|---|
| 1058 | vsum += coordVarDelta[lbl] * k |
|---|
| 1059 | print(head+line+') / '+('%.3f'%n)) |
|---|
| 1060 | fileval = modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][i]] |
|---|
| 1061 | print("{} = {:4g} (value read from CIF = {:4g})\n".format( |
|---|
| 1062 | self.Phase['ISODISTORT']['IsoModeList'][i], vsum, fileval)) |
|---|
| 1063 | |
|---|
| 1064 | print( 70*'=') |
|---|
| 1065 | print('Direct displacement relations computed from ISO modes in CIF', |
|---|
| 1066 | '\nusing Mode2VarMatrix, NormList, IsoModeList, IsoVarList',) |
|---|
| 1067 | # compute the coordinate displacements from the reported mode values |
|---|
| 1068 | for lbl,row in zip(self.Phase['ISODISTORT']['IsoVarList'], |
|---|
| 1069 | self.Phase['ISODISTORT']['Mode2VarMatrix']): |
|---|
| 1070 | l = '' |
|---|
| 1071 | s = 0.0 |
|---|
| 1072 | head = lbl+' =' |
|---|
| 1073 | for j,(k,n) in enumerate(zip(row,self.Phase['ISODISTORT']['NormList'])): |
|---|
| 1074 | if k == 0: continue |
|---|
| 1075 | if len(l) > 65: |
|---|
| 1076 | print(head,l) |
|---|
| 1077 | head = 20*' ' |
|---|
| 1078 | l = '' |
|---|
| 1079 | l1 = '' |
|---|
| 1080 | k1 = k |
|---|
| 1081 | if j > 0 and k < 0: |
|---|
| 1082 | k1 = -k |
|---|
| 1083 | l1 = ' - ' |
|---|
| 1084 | elif j > 0: |
|---|
| 1085 | l1 += ' + ' |
|---|
| 1086 | l += '{:} {:3g} * {:4g} * {:}'.format( |
|---|
| 1087 | l1, k1, n, self.Phase['ISODISTORT']['IsoModeList'][j]) |
|---|
| 1088 | |
|---|
| 1089 | s += n * modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][j]] * k |
|---|
| 1090 | print(head,l) |
|---|
| 1091 | print(lbl,'=',s) |
|---|
| 1092 | print(lbl,'==>',str(self.Phase['ISODISTORT']['G2VarList'][i]),'\n') |
|---|
| 1093 | DeltaCoords = {} |
|---|
| 1094 | for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix): |
|---|
| 1095 | s = 0.0 |
|---|
| 1096 | for j,(k,n) in enumerate(zip(row,self.Phase['ISODISTORT']['NormList'])): |
|---|
| 1097 | s += k * n * modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][j]] |
|---|
| 1098 | at,d = lbl.rsplit('_',1) |
|---|
| 1099 | if at not in DeltaCoords: |
|---|
| 1100 | DeltaCoords[at] = [0,0,0] |
|---|
| 1101 | if d == 'dx': |
|---|
| 1102 | DeltaCoords[at][0] = s |
|---|
| 1103 | elif d == 'dy': |
|---|
| 1104 | DeltaCoords[at][1] = s |
|---|
| 1105 | elif d == 'dz': |
|---|
| 1106 | DeltaCoords[at][2] = s |
|---|
| 1107 | #else: |
|---|
| 1108 | # print('unexpected',d) |
|---|
| 1109 | |
|---|
| 1110 | print( 70*'=') |
|---|
| 1111 | print('Coordinate checks') |
|---|
| 1112 | print("\nxyz's Computed from ISO mode values, as above") |
|---|
| 1113 | for at in sorted(DeltaCoords): |
|---|
| 1114 | s = at |
|---|
| 1115 | for i in range(3): |
|---|
| 1116 | s += ' ' |
|---|
| 1117 | s += str(ParentCoordinates[at][i]+DeltaCoords[at][i]) |
|---|
| 1118 | print(s) |
|---|
| 1119 | |
|---|
| 1120 | # determine the coordinate delta values from deviations from the parent structure |
|---|
| 1121 | print("\nxyz Values read directly from CIF") |
|---|
| 1122 | for atmline in self.Phase['Atoms']: |
|---|
| 1123 | lbl = atmline[0] |
|---|
| 1124 | x,y,z = atmline[3:6] |
|---|
| 1125 | print( lbl,x,y,z) |
|---|
| 1126 | |
|---|
| 1127 | print('\n' + 70*'=') |
|---|
| 1128 | print("G2 short name ==> ISODISTORT full name", |
|---|
| 1129 | " (from IsoModeList and G2ModeList)") |
|---|
| 1130 | for mode,G2mode in zip(self.Phase['ISODISTORT']['IsoModeList'], |
|---|
| 1131 | self.Phase['ISODISTORT']['G2ModeList']): |
|---|
| 1132 | print('{} ==> {}'.format(str(G2mode), mode)) |
|---|
| 1133 | print('\nConstraint help dict info') |
|---|
| 1134 | for i in self.Constraints: |
|---|
| 1135 | if type(i) is dict: |
|---|
| 1136 | for key in i: |
|---|
| 1137 | print('\t',key,':',i[key]) |
|---|
| 1138 | print( 70*'=') |
|---|
| 1139 | |
|---|
| 1140 | #====================================================================== |
|---|
| 1141 | # debug: show occupancy mode var to mode relations |
|---|
| 1142 | if debug and 'G2OccVarList' in self.Phase['ISODISTORT']: |
|---|
| 1143 | # coordinate items |
|---|
| 1144 | #occVarLbl = self.Phase['ISODISTORT']['OccVarList'] |
|---|
| 1145 | G2varObj = self.Phase['ISODISTORT']['G2OccVarList'] |
|---|
| 1146 | #ParentOcc = self.Phase['ISODISTORT']['BaseOcc'] |
|---|
| 1147 | # mode items |
|---|
| 1148 | modelist = self.Phase['ISODISTORT']['OccModeList'] |
|---|
| 1149 | modeVarList = self.Phase['ISODISTORT']['G2OccModeList'] |
|---|
| 1150 | normlist = self.Phase['ISODISTORT']['OccNormList'] |
|---|
| 1151 | # transform matrices |
|---|
| 1152 | #occupancymodeInvmatrix = self.Phase['ISODISTORT']['Var2OccMatrix'] |
|---|
| 1153 | #occupancymodematrix = self.Phase['ISODISTORT']['Occ2VarMatrix'] |
|---|
| 1154 | |
|---|
| 1155 | print( 70*'=') |
|---|
| 1156 | print('\nVar2OccMatrix' ,'OccVarList' ) |
|---|
| 1157 | for i,row in enumerate(occupancymodeInvmatrix): |
|---|
| 1158 | l = '' |
|---|
| 1159 | for j,(lbl,k) in enumerate(zip(occVarLbl,row)): |
|---|
| 1160 | if k == 0: continue |
|---|
| 1161 | if l: l += ' + ' |
|---|
| 1162 | #l += lbl+' * '+str(k) |
|---|
| 1163 | l += str(G2varObj[j])+' * '+str(k) |
|---|
| 1164 | print( str(i) + ': '+str(modeVarList[i])+' = '+l) |
|---|
| 1165 | |
|---|
| 1166 | # Get the ISODISTORT offset values |
|---|
| 1167 | occVarDelta = {} |
|---|
| 1168 | for lbl,val in zip( |
|---|
| 1169 | blk.get('_iso_deltaoccupancy_label'), |
|---|
| 1170 | blk.get('_iso_deltaoccupancy_value'),): |
|---|
| 1171 | occVarDelta[lbl] = float(val) |
|---|
| 1172 | modeVarDelta = {} |
|---|
| 1173 | for lbl,val in zip( |
|---|
| 1174 | blk.get('_iso_occupancymode_label'), |
|---|
| 1175 | blk.get('_iso_occupancymode_value'),): |
|---|
| 1176 | modeVarDelta[lbl] = cif.get_number_with_esd(val)[0] |
|---|
| 1177 | |
|---|
| 1178 | print( 70*'=') |
|---|
| 1179 | print('\nInverse relations using Var2OccModeMatrix, OccNormList, OccVarList') |
|---|
| 1180 | # compute the mode values from the reported coordinate deltas |
|---|
| 1181 | for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)): |
|---|
| 1182 | l = '' |
|---|
| 1183 | for lbl,k in zip(occVarLbl,row): |
|---|
| 1184 | if k == 0: continue |
|---|
| 1185 | if l: l += ' + ' |
|---|
| 1186 | l += lbl+' * '+str(k) |
|---|
| 1187 | print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n)) |
|---|
| 1188 | print('\nCalculation checks\n') |
|---|
| 1189 | for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)): |
|---|
| 1190 | #l = '' |
|---|
| 1191 | sl = '' |
|---|
| 1192 | s = 0. |
|---|
| 1193 | for lbl,k in zip(occVarLbl,row): |
|---|
| 1194 | if k == 0: continue |
|---|
| 1195 | #if l: l += ' + ' |
|---|
| 1196 | #l += lbl+' * '+str(k) |
|---|
| 1197 | if sl: sl += ' + ' |
|---|
| 1198 | sl += str(occVarDelta[lbl])+' * '+str(k) |
|---|
| 1199 | s += occVarDelta[lbl] * k |
|---|
| 1200 | print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n) |
|---|
| 1201 | print(' ?= ',modeVarDelta[modelist[i]]) |
|---|
| 1202 | print() |
|---|
| 1203 | |
|---|
| 1204 | print( 70*'=') |
|---|
| 1205 | print('\nDirect relations using Occ2VarMatrix, OccNormList, OccVarList') |
|---|
| 1206 | # compute the coordinate displacements from the reported mode values |
|---|
| 1207 | Occ = {} |
|---|
| 1208 | for i,lbl,row in zip(range(len(occVarLbl)),occVarLbl,occupancymodematrix): |
|---|
| 1209 | l = '' |
|---|
| 1210 | s = 0.0 |
|---|
| 1211 | for j,(k,n) in enumerate(zip(row,normlist)): |
|---|
| 1212 | if k == 0: continue |
|---|
| 1213 | if l: l += ' + ' |
|---|
| 1214 | l += str(n)+' * '+str(modeVarList[j])+' * '+str(k) |
|---|
| 1215 | s += n * modeVarDelta[modelist[j]] * k |
|---|
| 1216 | print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n') |
|---|
| 1217 | j = lbl.split('_')[0] |
|---|
| 1218 | Occ[j] = ParentOcc[j]+s |
|---|
| 1219 | |
|---|
| 1220 | # determine the coordinate delta values from deviations from the parent structure |
|---|
| 1221 | print('\nOccupancy from CIF vs computed') |
|---|
| 1222 | for atmline in self.Phase['Atoms']: |
|---|
| 1223 | lbl = atmline[0] |
|---|
| 1224 | if lbl in Occ: print( lbl,atmline[6],Occ[lbl]) |
|---|
| 1225 | |
|---|
| 1226 | print( 70*'=') |
|---|
| 1227 | print('\nGenerated constraints') |
|---|
| 1228 | for i in self.Constraints: |
|---|
| 1229 | if type(i) is dict: |
|---|
| 1230 | print('\nconstraint help dict') |
|---|
| 1231 | for key in i: |
|---|
| 1232 | print('\t',key,':',i[key]) |
|---|
| 1233 | elif i[-1] == 'f': |
|---|
| 1234 | print('\n\t',i[-3],' =') |
|---|
| 1235 | for m,j in i[:-3]: |
|---|
| 1236 | print('\t\t+',m,' * ',j,' ',repr(j)) |
|---|
| 1237 | else: |
|---|
| 1238 | print(' unexpected: ',repr(i)) |
|---|
| 1239 | print("\nG2name ==> ISODISTORT full name", |
|---|
| 1240 | ".Phase['ISODISTORT']['OccModeList']", |
|---|
| 1241 | ".Phase['ISODISTORT']['G2OccModeList']") |
|---|
| 1242 | for mode,G2mode in zip(modelist,modeVarList): |
|---|
| 1243 | print(" ?::"+str(G2mode),' ==>', mode) |
|---|
| 1244 | |
|---|
| 1245 | def ISODISTORT_shortLbl(lbl,shortmodelist): |
|---|
| 1246 | '''Shorten model labels and remove special characters |
|---|
| 1247 | ''' |
|---|
| 1248 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)\[.*?\](.*)',lbl) |
|---|
| 1249 | # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)[...]BBBBB |
|---|
| 1250 | # where SSSSS is the parent spacegroup, [x,y,z] is a location, etc. |
|---|
| 1251 | if regexp: |
|---|
| 1252 | # this extracts the AAAAA and BBBBB parts of the string |
|---|
| 1253 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 1254 | else: |
|---|
| 1255 | # try form SSSSS[x,y,z]AAAA(a,b,...)BBBBB |
|---|
| 1256 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl) |
|---|
| 1257 | if regexp: |
|---|
| 1258 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 1259 | lbl = lbl.replace('order','o') |
|---|
| 1260 | lbl = lbl.replace('+','_') |
|---|
| 1261 | lbl = lbl.replace('-','_') |
|---|
| 1262 | G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list |
|---|
| 1263 | |
|---|
