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