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