| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | ########### SVN repository information ################### |
|---|
| 3 | # $Date: 2018-04-04 14:47:03 +0000 (Wed, 04 Apr 2018) $ |
|---|
| 4 | # $Author: vondreele $ |
|---|
| 5 | # $Revision: 3335 $ |
|---|
| 6 | # $URL: trunk/imports/G2phase_CIF.py $ |
|---|
| 7 | # $Id: G2phase_CIF.py 3335 2018-04-04 14:47:03Z vondreele $ |
|---|
| 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 GSASIIIO as G2IO |
|---|
| 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 | GSASIIpath.SetVersionNumber("$Revision: 3335 $") |
|---|
| 34 | import CifFile as cif # PyCifRW from James Hester |
|---|
| 35 | |
|---|
| 36 | class CIFPhaseReader(G2obj.ImportPhase): |
|---|
| 37 | 'Implements a phase importer from a possibly multi-block CIF file' |
|---|
| 38 | def __init__(self): |
|---|
| 39 | super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__ |
|---|
| 40 | extensionlist=('.CIF','.cif','.mcif'), |
|---|
| 41 | strictExtension=False, |
|---|
| 42 | formatName = 'CIF', |
|---|
| 43 | longFormatName = 'Crystallographic Information File import' |
|---|
| 44 | ) |
|---|
| 45 | |
|---|
| 46 | def ContentsValidator(self, filename): |
|---|
| 47 | fp = open(filename,'r') |
|---|
| 48 | return self.CIFValidator(fp) |
|---|
| 49 | fp.close() |
|---|
| 50 | |
|---|
| 51 | def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused): |
|---|
| 52 | self.isodistort_warnings = '' |
|---|
| 53 | self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 54 | # make sure the ranId is really unique! |
|---|
| 55 | while self.Phase['ranId'] in usedRanIdList: |
|---|
| 56 | self.Phase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 57 | self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 58 | # make sure the ranId is really unique! |
|---|
| 59 | while self.MPhase['ranId'] in usedRanIdList: |
|---|
| 60 | self.MPhase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 61 | returnstat = False |
|---|
| 62 | cellitems = ( |
|---|
| 63 | '_cell_length_a','_cell_length_b','_cell_length_c', |
|---|
| 64 | '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',) |
|---|
| 65 | # cellwaveitems = ( |
|---|
| 66 | # '_cell_wave_vector_seq_id', |
|---|
| 67 | # '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z') |
|---|
| 68 | reqitems = ( |
|---|
| 69 | '_atom_site_fract_x', |
|---|
| 70 | '_atom_site_fract_y', |
|---|
| 71 | '_atom_site_fract_z', |
|---|
| 72 | ) |
|---|
| 73 | phasenamefields = ( |
|---|
| 74 | '_chemical_name_common', |
|---|
| 75 | '_pd_phase_name', |
|---|
| 76 | '_chemical_formula_sum' |
|---|
| 77 | ) |
|---|
| 78 | try: |
|---|
| 79 | cf = G2obj.ReadCIF(filename) |
|---|
| 80 | except cif.StarError as msg: |
|---|
| 81 | msg = 'Unreadable cif file\n'+str(msg) |
|---|
| 82 | self.errors = msg |
|---|
| 83 | self.warnings += msg |
|---|
| 84 | return False |
|---|
| 85 | # scan blocks for structural info |
|---|
| 86 | self.errors = 'Error during scan of blocks for datasets' |
|---|
| 87 | str_blklist = [] |
|---|
| 88 | for blk in cf.keys(): |
|---|
| 89 | for r in reqitems+cellitems: |
|---|
| 90 | if r not in cf[blk].keys(): |
|---|
| 91 | break |
|---|
| 92 | else: |
|---|
| 93 | str_blklist.append(blk) |
|---|
| 94 | if not str_blklist: |
|---|
| 95 | selblk = None # no block to choose |
|---|
| 96 | elif len(str_blklist) == 1: # only one choice |
|---|
| 97 | selblk = 0 |
|---|
| 98 | else: # choose from options |
|---|
| 99 | choice = [] |
|---|
| 100 | for blknm in str_blklist: |
|---|
| 101 | choice.append('') |
|---|
| 102 | # accumumlate some info about this phase |
|---|
| 103 | choice[-1] += blknm + ': ' |
|---|
| 104 | for i in phasenamefields: # get a name for the phase |
|---|
| 105 | name = cf[blknm].get(i,'phase name').strip() |
|---|
| 106 | if name is None or name == '?' or name == '.': |
|---|
| 107 | continue |
|---|
| 108 | else: |
|---|
| 109 | choice[-1] += name.strip()[:20] + ', ' |
|---|
| 110 | break |
|---|
| 111 | na = len(cf[blknm].get("_atom_site_fract_x")) |
|---|
| 112 | if na == 1: |
|---|
| 113 | choice[-1] += '1 atom' |
|---|
| 114 | else: |
|---|
| 115 | choice[-1] += ('%d' % na) + ' atoms' |
|---|
| 116 | choice[-1] += ', cell: ' |
|---|
| 117 | fmt = "%.2f," |
|---|
| 118 | for i,key in enumerate(cellitems): |
|---|
| 119 | if i == 3: fmt = "%.f," |
|---|
| 120 | if i == 5: fmt = "%.f" |
|---|
| 121 | choice[-1] += fmt % cif.get_number_with_esd( |
|---|
| 122 | cf[blknm].get(key))[0] |
|---|
| 123 | sg = cf[blknm].get("_symmetry_space_group_name_H-M",'') |
|---|
| 124 | if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'') |
|---|
| 125 | if not sg: sg = cf[blknm].get("_space_group_ssg_name",'') |
|---|
| 126 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 127 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'') |
|---|
| 128 | #how about checking for super/magnetic ones as well? - reject 'X'? |
|---|
| 129 | sg = sg.replace('_','') |
|---|
| 130 | if sg: choice[-1] += ', (' + sg.strip() + ')' |
|---|
| 131 | selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame, |
|---|
| 132 | title= 'Select a phase from one the CIF data_ blocks below',size=(600,100)) |
|---|
| 133 | self.errors = 'Error during reading of selected block' |
|---|
| 134 | #process selected phase |
|---|
| 135 | if selblk is None: |
|---|
| 136 | returnstat = False # no block selected or available |
|---|
| 137 | else: #do space group symbol & phase type first |
|---|
| 138 | blknm = str_blklist[selblk] |
|---|
| 139 | blk = cf[str_blklist[selblk]] |
|---|
| 140 | E = True |
|---|
| 141 | Super = False |
|---|
| 142 | magnetic = False |
|---|
| 143 | moddim = int(blk.get("_cell_modulation_dimension",'0')) |
|---|
| 144 | if moddim: #incommensurate |
|---|
| 145 | if moddim > 1: |
|---|
| 146 | msg = 'more than 3+1 super symmetry is not allowed in GSAS-II' |
|---|
| 147 | self.errors = msg |
|---|
| 148 | self.warnings += '\n'+msg |
|---|
| 149 | return False |
|---|
| 150 | if blk.get('_cell_subsystems_number'): |
|---|
| 151 | msg = 'Composite super structures not allowed in GSAS-II' |
|---|
| 152 | self.errors = msg |
|---|
| 153 | self.warnings += '\n'+msg |
|---|
| 154 | return False |
|---|
| 155 | sspgrp = blk.get("_space_group_ssg_name",'') |
|---|
| 156 | if not sspgrp: #might be incommensurate magnetic |
|---|
| 157 | MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 158 | if not MSSpGrp: |
|---|
| 159 | MSSpGrp = blk.get("_space_group.magn_ssg_name",'') |
|---|
| 160 | if not MSSpGrp: |
|---|
| 161 | msg = 'No incommensurate space group name was found in the CIF.' |
|---|
| 162 | self.errors = msg |
|---|
| 163 | self.warnings += '\n'+msg |
|---|
| 164 | return False |
|---|
| 165 | if 'X' in MSSpGrp: |
|---|
| 166 | msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II' |
|---|
| 167 | self.warnings += '\n'+msg |
|---|
| 168 | self.errors = msg |
|---|
| 169 | return False |
|---|
| 170 | magnetic = True |
|---|
| 171 | if 'X' in sspgrp: |
|---|
| 172 | msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II' |
|---|
| 173 | self.warnings += '\n'+msg |
|---|
| 174 | self.errors = msg |
|---|
| 175 | return False |
|---|
| 176 | Super = True |
|---|
| 177 | if magnetic: |
|---|
| 178 | sspgrp = MSSpGrp.split('(') |
|---|
| 179 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 180 | SpGrp = G2spc.StandardizeSpcName(sspgrp[0]) |
|---|
| 181 | if "1'" in SpGrp: sspgrp[1] = sspgrp[1][:-1] #take off extra 's'; gets put back later |
|---|
| 182 | MSpGrp = sspgrp[0] |
|---|
| 183 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 184 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 185 | else: |
|---|
| 186 | sspgrp = sspgrp.split('(') |
|---|
| 187 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 188 | SpGrp = sspgrp[0] |
|---|
| 189 | SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 190 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 191 | if not SpGrp: |
|---|
| 192 | print (sspgrp) |
|---|
| 193 | self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II' |
|---|
| 194 | return False |
|---|
| 195 | SuperSg = sspgrp[1].replace('\\','') |
|---|
| 196 | SuperVec = [[0,0,.1],False,4] |
|---|
| 197 | else: #not incommensurate |
|---|
| 198 | SpGrp = blk.get("_symmetry_space_group_name_H-M",'') |
|---|
| 199 | if not SpGrp: |
|---|
| 200 | SpGrp = blk.get("_space_group_name_H-M_alt",'') |
|---|
| 201 | if not SpGrp: #try magnetic |
|---|
| 202 | MSpGrp = blk.get("_space_group.magn_name_BNS",'') |
|---|
| 203 | if not MSpGrp: |
|---|
| 204 | MSpGrp = blk.get("_space_group_magn.name_BNS",'') |
|---|
| 205 | if not MSpGrp: |
|---|
| 206 | msg = 'No recognizable space group name was found in the CIF.' |
|---|
| 207 | self.errors = msg |
|---|
| 208 | self.warnings += '\n'+msg |
|---|
| 209 | return False |
|---|
| 210 | SpGrp = blk.get('_parent_space_group.name_H-M_alt') |
|---|
| 211 | if not SpGrp: |
|---|
| 212 | SpGrp = blk.get('_parent_space_group.name_H-M') |
|---|
| 213 | # SpGrp = MSpGrp.replace("'",'') |
|---|
| 214 | SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','') #get rid of screw '_' |
|---|
| 215 | if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:] |
|---|
| 216 | SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 217 | magnetic = True |
|---|
| 218 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 219 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 220 | if not SpGrp: |
|---|
| 221 | print (MSpGrp) |
|---|
| 222 | self.warnings += 'No space group name was found in the CIF.' |
|---|
| 223 | return False |
|---|
| 224 | else: |
|---|
| 225 | SpGrp = SpGrp.replace('_','') |
|---|
| 226 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 227 | #process space group symbol |
|---|
| 228 | E,SGData = G2spc.SpcGroup(SpGrp) |
|---|
| 229 | if E and SpGrp: |
|---|
| 230 | SpGrpNorm = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 231 | if SpGrpNorm: |
|---|
| 232 | E,SGData = G2spc.SpcGroup(SpGrpNorm) |
|---|
| 233 | # nope, try the space group "out of the Box" |
|---|
| 234 | if E: |
|---|
| 235 | self.warnings += 'ERROR in space group symbol '+SpGrp |
|---|
| 236 | self.warnings += '\nThe space group has been set to "P 1". ' |
|---|
| 237 | self.warnings += "Change this in phase's General tab." |
|---|
| 238 | self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: ' |
|---|
| 239 | self.warnings += G2spc.SGErrors(E) |
|---|
| 240 | SGData = G2obj.P1SGData # P 1 |
|---|
| 241 | self.Phase['General']['SGData'] = SGData |
|---|
| 242 | if Super: |
|---|
| 243 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 244 | if E: |
|---|
| 245 | self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg |
|---|
| 246 | self.warnings += '\n'+E |
|---|
| 247 | SuperSg = SuperSg[:SuperSg.index(')')+1] |
|---|
| 248 | self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg |
|---|
| 249 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 250 | self.Phase['General']['SSGData'] = SSGData |
|---|
| 251 | if magnetic: |
|---|
| 252 | self.MPhase['General']['SGData'] = SGData |
|---|
| 253 | self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','') |
|---|
| 254 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 255 | |
|---|
| 256 | if magnetic: #replace std operaors with those from cif file - probably not the same! |
|---|
| 257 | SGData['SGFixed'] = True |
|---|
| 258 | SGData['SGOps'] = [] |
|---|
| 259 | SGData['SGCen'] = [] |
|---|
| 260 | if Super: |
|---|
| 261 | SSGData['SSGOps'] = [] |
|---|
| 262 | SSGData['SSGCen'] = [] |
|---|
| 263 | try: |
|---|
| 264 | sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id') |
|---|
| 265 | sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id') |
|---|
| 266 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1] |
|---|
| 267 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1] |
|---|
| 268 | except KeyError: #old mag cif names |
|---|
| 269 | sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id') |
|---|
| 270 | sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id') |
|---|
| 271 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1] |
|---|
| 272 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1] |
|---|
| 273 | spnflp = [] |
|---|
| 274 | for op in sgoploop: |
|---|
| 275 | M,T,S = G2spc.MagSSText2MTS(op[opid]) |
|---|
| 276 | SSGData['SSGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 277 | SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]]) |
|---|
| 278 | spnflp.append(S) |
|---|
| 279 | censpn = [] |
|---|
| 280 | for cent in sgcenloop: |
|---|
| 281 | M,C,S = G2spc.MagSSText2MTS(cent[centid]) |
|---|
| 282 | SSGData['SSGCen'].append(C) |
|---|
| 283 | SGData['SGCen'].append(C[:3]) |
|---|
| 284 | censpn += list(np.array(spnflp)*S) |
|---|
| 285 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 286 | else: |
|---|
| 287 | try: |
|---|
| 288 | sgoploop = blk.GetLoop('_space_group_symop_magn.id') |
|---|
| 289 | sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id') |
|---|
| 290 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1] |
|---|
| 291 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1] |
|---|
| 292 | except KeyError: #old mag cif names |
|---|
| 293 | sgoploop = blk.GetLoop('_space_group_symop.magn_id') |
|---|
| 294 | sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id') |
|---|
| 295 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1] |
|---|
| 296 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1] |
|---|
| 297 | spnflp = [] |
|---|
| 298 | for op in sgoploop: |
|---|
| 299 | try: |
|---|
| 300 | M,T,S = G2spc.MagText2MTS(op[opid]) |
|---|
| 301 | SGData['SGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 302 | spnflp.append(S) |
|---|
| 303 | except KeyError: |
|---|
| 304 | self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II' |
|---|
| 305 | return False |
|---|
| 306 | censpn = [] |
|---|
| 307 | for cent in sgcenloop: |
|---|
| 308 | M,C,S = G2spc.MagText2MTS(cent[centid]) |
|---|
| 309 | SGData['SGCen'].append(C) |
|---|
| 310 | censpn += list(np.array(spnflp)*S) |
|---|
| 311 | self.MPhase['General']['SGData'] = SGData |
|---|
| 312 | self.MPhase['General']['SGData']['SpnFlp'] = censpn |
|---|
| 313 | self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp |
|---|
| 314 | MagPtGp = blk.get('_space_group.magn_point_group') |
|---|
| 315 | if not MagPtGp: |
|---|
| 316 | MagPtGp = blk.get('_space_group.magn_point_group_name') |
|---|
| 317 | if not MagPtGp: |
|---|
| 318 | MagPtGp = blk.get('_space_group_magn.point_group_name') |
|---|
| 319 | self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp |
|---|
| 320 | |
|---|
| 321 | # cell parameters |
|---|
| 322 | cell = [] |
|---|
| 323 | for lbl in cellitems: |
|---|
| 324 | cell.append(cif.get_number_with_esd(blk[lbl])[0]) |
|---|
| 325 | Volume = G2lat.calc_V(G2lat.cell2A(cell)) |
|---|
| 326 | self.Phase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 327 | if magnetic: |
|---|
| 328 | self.MPhase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 329 | if Super: |
|---|
| 330 | waveloop = blk.GetLoop('_cell_wave_vector_seq_id') |
|---|
| 331 | waveDict = dict(waveloop.items()) |
|---|
| 332 | SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0], |
|---|
| 333 | cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0], |
|---|
| 334 | cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4] |
|---|
| 335 | |
|---|
| 336 | # read in atoms |
|---|
| 337 | self.errors = 'Error during reading of atoms' |
|---|
| 338 | atomlbllist = [] # table to look up atom IDs |
|---|
| 339 | atomloop = blk.GetLoop('_atom_site_label') |
|---|
| 340 | atomkeys = [i.lower() for i in atomloop.keys()] |
|---|
| 341 | if not blk.get('_atom_site_type_symbol'): |
|---|
| 342 | self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed' |
|---|
| 343 | if magnetic: |
|---|
| 344 | try: |
|---|
| 345 | magmoment = '_atom_site_moment.label' |
|---|
| 346 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 347 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 348 | magatomlabels = blk.get(magmoment) |
|---|
| 349 | G2MagDict = {'_atom_site_moment.label': 0, |
|---|
| 350 | '_atom_site_moment.crystalaxis_x':7, |
|---|
| 351 | '_atom_site_moment.crystalaxis_y':8, |
|---|
| 352 | '_atom_site_moment.crystalaxis_z':9} |
|---|
| 353 | except KeyError: |
|---|
| 354 | magmoment = '_atom_site_moment_label' |
|---|
| 355 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 356 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 357 | magatomlabels = blk.get(magmoment) |
|---|
| 358 | G2MagDict = {'_atom_site_moment_label': 0, |
|---|
| 359 | '_atom_site_moment_crystalaxis_x':7, |
|---|
| 360 | '_atom_site_moment_crystalaxis_y':8, |
|---|
| 361 | '_atom_site_moment_crystalaxis_z':9} |
|---|
| 362 | |
|---|
| 363 | if blk.get('_atom_site_aniso_label'): |
|---|
| 364 | anisoloop = blk.GetLoop('_atom_site_aniso_label') |
|---|
| 365 | anisokeys = [i.lower() for i in anisoloop.keys()] |
|---|
| 366 | anisolabels = blk.get('_atom_site_aniso_label') |
|---|
| 367 | else: |
|---|
| 368 | anisoloop = None |
|---|
| 369 | anisokeys = [] |
|---|
| 370 | anisolabels = [] |
|---|
| 371 | if Super: |
|---|
| 372 | occFloop = None |
|---|
| 373 | occCloop = None |
|---|
| 374 | occFdict = {} |
|---|
| 375 | occCdict = {} |
|---|
| 376 | displSloop = None |
|---|
| 377 | displFloop = None |
|---|
| 378 | MagFloop = None |
|---|
| 379 | displSdict = {} |
|---|
| 380 | displFdict = {} |
|---|
| 381 | MagFdict = {} |
|---|
| 382 | UijFloop = None |
|---|
| 383 | UijFdict = {} |
|---|
| 384 | #occupancy modulation |
|---|
| 385 | if blk.get('_atom_site_occ_Fourier_atom_site_label'): |
|---|
| 386 | occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label') |
|---|
| 387 | occFdict = dict(occFloop.items()) |
|---|
| 388 | if blk.get('_atom_site_occ_special_func_atom_site_label'): #Crenel (i.e. Block Wave) occ |
|---|
| 389 | occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label') |
|---|
| 390 | occCdict = dict(occCloop.items()) |
|---|
| 391 | #position modulation |
|---|
| 392 | if blk.get('_atom_site_displace_Fourier_atom_site_label'): |
|---|
| 393 | displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label') |
|---|
| 394 | displFdict = dict(displFloop.items()) |
|---|
| 395 | if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth |
|---|
| 396 | displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label') |
|---|
| 397 | displSdict = dict(displSloop.items()) |
|---|
| 398 | #U modulation |
|---|
| 399 | if blk.get('_atom_site_U_Fourier_atom_site_label'): |
|---|
| 400 | UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label') |
|---|
| 401 | UijFdict = dict(UijFloop.items()) |
|---|
| 402 | #Mag moment modulation |
|---|
| 403 | if blk.get('_atom_site_moment_Fourier_atom_site_label'): |
|---|
| 404 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label') |
|---|
| 405 | MagFdict = dict(MagFloop.items()) |
|---|
| 406 | Mnames = ['_atom_site_moment_fourier_atom_site_label', |
|---|
| 407 | '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id', |
|---|
| 408 | '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos'] |
|---|
| 409 | elif blk.get('_atom_site_moment_Fourier.atom_site_label'): |
|---|
| 410 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label') |
|---|
| 411 | MagFdict = dict(MagFloop.items()) |
|---|
| 412 | Mnames = ['_atom_site_moment_fourier.atom_site_label', |
|---|
| 413 | '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id', |
|---|
| 414 | '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos'] |
|---|
| 415 | self.Phase['Atoms'] = [] |
|---|
| 416 | if magnetic: |
|---|
| 417 | self.MPhase['Atoms'] = [] |
|---|
| 418 | G2AtomDict = { '_atom_site_type_symbol' : 1, |
|---|
| 419 | '_atom_site_label' : 0, |
|---|
| 420 | '_atom_site_fract_x' : 3, |
|---|
| 421 | '_atom_site_fract_y' : 4, |
|---|
| 422 | '_atom_site_fract_z' : 5, |
|---|
| 423 | '_atom_site_occupancy' : 6, |
|---|
| 424 | '_atom_site_aniso_u_11' : 11, |
|---|
| 425 | '_atom_site_aniso_u_22' : 12, |
|---|
| 426 | '_atom_site_aniso_u_33' : 13, |
|---|
| 427 | '_atom_site_aniso_u_12' : 14, |
|---|
| 428 | '_atom_site_aniso_u_13' : 15, |
|---|
| 429 | '_atom_site_aniso_u_23' : 16, } |
|---|
| 430 | |
|---|
| 431 | ranIdlookup = {} |
|---|
| 432 | for aitem in atomloop: |
|---|
| 433 | atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,] |
|---|
| 434 | for val,key in zip(aitem,atomkeys): |
|---|
| 435 | col = G2AtomDict.get(key,-1) |
|---|
| 436 | if col >= 3: |
|---|
| 437 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 438 | if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag |
|---|
| 439 | elif col is not None and col != -1: |
|---|
| 440 | atomlist[col] = val |
|---|
| 441 | elif key in ('_atom_site_thermal_displace_type', |
|---|
| 442 | '_atom_site_adp_type'): #Iso or Aniso? |
|---|
| 443 | if val.lower() == 'uani': |
|---|
| 444 | atomlist[9] = 'A' |
|---|
| 445 | elif key == '_atom_site_u_iso_or_equiv': |
|---|
| 446 | uisoval = cif.get_number_with_esd(val)[0] |
|---|
| 447 | if uisoval is not None: |
|---|
| 448 | atomlist[10] = uisoval |
|---|
| 449 | if not atomlist[1] and atomlist[0]: |
|---|
| 450 | typ = atomlist[0].rstrip('0123456789-+') |
|---|
| 451 | if G2elem.CheckElement(typ): |
|---|
| 452 | atomlist[1] = typ |
|---|
| 453 | if not atomlist[1]: |
|---|
| 454 | atomlist[1] = 'Xe' |
|---|
| 455 | self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n' |
|---|
| 456 | if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop? |
|---|
| 457 | atomlist[9] = 'A' |
|---|
| 458 | for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys): |
|---|
| 459 | col = G2AtomDict.get(key) |
|---|
| 460 | if col: |
|---|
| 461 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 462 | atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2] |
|---|
| 463 | atomlist[1] = G2elem.FixValence(atomlist[1]) |
|---|
| 464 | atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id |
|---|
| 465 | self.Phase['Atoms'].append(atomlist) |
|---|
| 466 | ranIdlookup[atomlist[0]] = atomlist[-1] |
|---|
| 467 | if atomlist[0] in atomlbllist: |
|---|
| 468 | self.warnings += ' ERROR: repeated atom label: '+atomlist[0] |
|---|
| 469 | else: |
|---|
| 470 | atomlbllist.append(atomlist[0]) |
|---|
| 471 | |
|---|
| 472 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 473 | matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:] |
|---|
| 474 | for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys): |
|---|
| 475 | mcol = G2MagDict.get(mkey,-1) |
|---|
| 476 | if mcol: |
|---|
| 477 | matomlist[mcol] = cif.get_number_with_esd(mval)[0] |
|---|
| 478 | self.MPhase['Atoms'].append(matomlist) |
|---|
| 479 | if Super: |
|---|
| 480 | Sfrac = np.zeros((4,2)) |
|---|
| 481 | Sadp = np.zeros((4,12)) |
|---|
| 482 | Spos = np.zeros((4,6)) |
|---|
| 483 | Smag = np.zeros((4,6)) |
|---|
| 484 | nim = -1 |
|---|
| 485 | waveType = 'Fourier' |
|---|
| 486 | if occCdict: |
|---|
| 487 | for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']): |
|---|
| 488 | if item == atomlist[0]: |
|---|
| 489 | waveType = 'Crenel' |
|---|
| 490 | val = occCdict['_atom_site_occ_special_func_crenel_c'][i] |
|---|
| 491 | Sfrac[0][0] = cif.get_number_with_esd(val)[0] |
|---|
| 492 | val = occCdict['_atom_site_occ_special_func_crenel_w'][i] |
|---|
| 493 | Sfrac[0][1] = cif.get_number_with_esd(val)[0] |
|---|
| 494 | nim = 1 |
|---|
| 495 | |
|---|
| 496 | if nim >= 0: |
|---|
| 497 | Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]] |
|---|
| 498 | else: |
|---|
| 499 | Sfrac = [] |
|---|
| 500 | nim = -1 |
|---|
| 501 | if displFdict: |
|---|
| 502 | for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']): |
|---|
| 503 | if item == atomlist[0]: |
|---|
| 504 | waveType = 'Fourier' |
|---|
| 505 | ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i]) |
|---|
| 506 | im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i]) |
|---|
| 507 | if im != nim: |
|---|
| 508 | nim = im |
|---|
| 509 | val = displFdict['_atom_site_displace_fourier_param_sin'][i] |
|---|
| 510 | Spos[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 511 | val = displFdict['_atom_site_displace_fourier_param_cos'][i] |
|---|
| 512 | Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 513 | if nim >= 0: |
|---|
| 514 | Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]] |
|---|
| 515 | else: |
|---|
| 516 | Spos = [] |
|---|
| 517 | nim = -1 |
|---|
| 518 | if UijFdict: |
|---|
| 519 | nim = -1 |
|---|
| 520 | for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']): |
|---|
| 521 | if item == atomlist[0]: |
|---|
| 522 | ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i]) |
|---|
| 523 | im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i]) |
|---|
| 524 | if im != nim: |
|---|
| 525 | nim = im |
|---|
| 526 | val = UijFdict['_atom_site_u_fourier_param_sin'][i] |
|---|
| 527 | Sadp[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 528 | val = UijFdict['_atom_site_u_fourier_param_cos'][i] |
|---|
| 529 | Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0] |
|---|
| 530 | if nim >= 0: |
|---|
| 531 | Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]] |
|---|
| 532 | else: |
|---|
| 533 | Sadp = [] |
|---|
| 534 | nim = -1 |
|---|
| 535 | if MagFdict: |
|---|
| 536 | nim = -1 |
|---|
| 537 | for i,item in enumerate(MagFdict[Mnames[0]]): |
|---|
| 538 | if item == atomlist[0]: |
|---|
| 539 | ix = ['x','y','z'].index(MagFdict[Mnames[1]][i]) |
|---|
| 540 | im = int(MagFdict[Mnames[2]][i]) |
|---|
| 541 | if im != nim: |
|---|
| 542 | nim = im |
|---|
| 543 | val = MagFdict[Mnames[3]][i] |
|---|
| 544 | Smag[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 545 | val = MagFdict[Mnames[4]][i] |
|---|
| 546 | Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 547 | if nim >= 0: |
|---|
| 548 | Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]] |
|---|
| 549 | else: |
|---|
| 550 | Smag = [] |
|---|
| 551 | SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}} |
|---|
| 552 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 553 | matomlist.append(SSdict) |
|---|
| 554 | atomlist.append(SSdict) |
|---|
| 555 | if len(atomlbllist) != len(self.Phase['Atoms']): |
|---|
| 556 | self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode' |
|---|
| 557 | for lbl in phasenamefields: # get a name for the phase |
|---|
| 558 | name = blk.get(lbl) |
|---|
| 559 | if name is None: |
|---|
| 560 | continue |
|---|
| 561 | name = name.strip() |
|---|
| 562 | if name == '?' or name == '.': |
|---|
| 563 | continue |
|---|
| 564 | else: |
|---|
| 565 | break |
|---|
| 566 | else: # no name found, use block name for lack of a better choice |
|---|
| 567 | name = blknm |
|---|
| 568 | self.Phase['General']['Name'] = name.strip()[:20] |
|---|
| 569 | self.Phase['General']['Super'] = Super |
|---|
| 570 | if magnetic: |
|---|
| 571 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 572 | self.MPhase['General']['Name'] = name.strip()[:20]+' mag' |
|---|
| 573 | self.MPhase['General']['Super'] = Super |
|---|
| 574 | if Super: |
|---|
| 575 | if self.MPhase['General']['SGData']['SGGray']: |
|---|
| 576 | SuperSg += 's' |
|---|
| 577 | self.MPhase['General']['Modulated'] = True |
|---|
| 578 | self.MPhase['General']['SuperVec'] = SuperVec |
|---|
| 579 | self.MPhase['General']['SuperSg'] = SuperSg |
|---|
| 580 | else: |
|---|
| 581 | self.MPhase = None |
|---|
| 582 | if Super: |
|---|
| 583 | self.Phase['General']['Modulated'] = True |
|---|
| 584 | self.Phase['General']['SuperVec'] = SuperVec |
|---|
| 585 | self.Phase['General']['SuperSg'] = SuperSg |
|---|
| 586 | if not self.Phase['General']['SGData']['SGFixed']: |
|---|
| 587 | self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1] |
|---|
| 588 | if not self.isodistort_warnings: |
|---|
| 589 | if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'): |
|---|
| 590 | self.errors = "Error while processing ISODISTORT constraints" |
|---|
| 591 | self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup) |
|---|
| 592 | else: |
|---|
| 593 | self.warnings += self.isodistort_warnings |
|---|
| 594 | returnstat = True |
|---|
| 595 | return returnstat |
|---|
| 596 | |
|---|
| 597 | def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup): |
|---|
| 598 | 'Process ISODISTORT items to create constraints etc.' |
|---|
| 599 | varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'} |
|---|
| 600 | 'Maps ISODISTORT parm names to GSAS-II names' |
|---|
| 601 | #---------------------------------------------------------------------- |
|---|
| 602 | # read in the ISODISTORT displacement modes |
|---|
| 603 | #---------------------------------------------------------------------- |
|---|
| 604 | self.Constraints = [] |
|---|
| 605 | explaination = {} |
|---|
| 606 | if blk.get('_iso_displacivemode_label'): |
|---|
| 607 | modelist = [] |
|---|
| 608 | shortmodelist = [] |
|---|
| 609 | for lbl in blk.get('_iso_displacivemode_label'): |
|---|
| 610 | modelist.append(lbl) |
|---|
| 611 | # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB |
|---|
| 612 | # where SSSSS is the parent spacegroup, [x,y,z] is a location |
|---|
| 613 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl) |
|---|
| 614 | # this extracts the AAAAA and BBBBB parts of the string |
|---|
| 615 | if regexp: |
|---|
| 616 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 617 | G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list |
|---|
| 618 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 619 | coordVarLbl = [] |
|---|
| 620 | G2varLbl = [] |
|---|
| 621 | G2varObj = [] |
|---|
| 622 | error = False |
|---|
| 623 | for lbl in blk.get('_iso_deltacoordinate_label'): |
|---|
| 624 | coordVarLbl.append(lbl) |
|---|
| 625 | if '_' in lbl: |
|---|
| 626 | albl = lbl[:lbl.rfind('_')] |
|---|
| 627 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 628 | else: |
|---|
| 629 | self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl |
|---|
| 630 | error = True |
|---|
| 631 | continue |
|---|
| 632 | if albl not in atomlbllist: |
|---|
| 633 | self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl |
|---|
| 634 | error = True |
|---|
| 635 | continue |
|---|
| 636 | else: |
|---|
| 637 | anum = atomlbllist.index(albl) |
|---|
| 638 | var = varLookup.get(vlbl) |
|---|
| 639 | if not var: |
|---|
| 640 | self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl |
|---|
| 641 | error = True |
|---|
| 642 | continue |
|---|
| 643 | G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID |
|---|
| 644 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 645 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 646 | )) |
|---|
| 647 | if error: |
|---|
| 648 | raise Exception("Error decoding variable labels") |
|---|
| 649 | |
|---|
| 650 | if len(G2varObj) != len(modelist): |
|---|
| 651 | print ("non-square input") |
|---|
| 652 | raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate") |
|---|
| 653 | |
|---|
| 654 | error = False |
|---|
| 655 | ParentCoordinates = {} |
|---|
| 656 | for lbl,exp in zip( |
|---|
| 657 | blk.get('_iso_coordinate_label'), |
|---|
| 658 | blk.get('_iso_coordinate_formula'), |
|---|
| 659 | ): |
|---|
| 660 | if '_' in lbl: |
|---|
| 661 | albl = lbl[:lbl.rfind('_')] |
|---|
| 662 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 663 | else: |
|---|
| 664 | self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl |
|---|
| 665 | error = True |
|---|
| 666 | continue |
|---|
| 667 | if vlbl not in 'xyz' or len(vlbl) != 1: |
|---|
| 668 | self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl |
|---|
| 669 | error = True |
|---|
| 670 | continue |
|---|
| 671 | i = 'xyz'.index(vlbl) |
|---|
| 672 | if not ParentCoordinates.get(albl): |
|---|
| 673 | ParentCoordinates[albl] = [None,None,None] |
|---|
| 674 | if '+' in exp: |
|---|
| 675 | val = exp.split('+')[0].strip() |
|---|
| 676 | val = G2p3.FormulaEval(val) |
|---|
| 677 | if val is None: |
|---|
| 678 | self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl |
|---|
| 679 | error = True |
|---|
| 680 | continue |
|---|
| 681 | ParentCoordinates[albl][i] = val |
|---|
| 682 | else: |
|---|
| 683 | ParentCoordinates[albl][i] = G2p3.FormulaEval(exp) |
|---|
| 684 | if error: |
|---|
| 685 | print (self.warnings) |
|---|
| 686 | raise Exception("Error decoding variable labels") |
|---|
| 687 | # get mapping of modes to atomic coordinate displacements |
|---|
| 688 | displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 689 | for row,col,val in zip( |
|---|
| 690 | blk.get('_iso_displacivemodematrix_row'), |
|---|
| 691 | blk.get('_iso_displacivemodematrix_col'), |
|---|
| 692 | blk.get('_iso_displacivemodematrix_value'),): |
|---|
| 693 | displacivemodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 694 | # Invert to get mapping of atom displacements to modes |
|---|
| 695 | displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix) |
|---|
| 696 | # create the constraints |
|---|
| 697 | for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 698 | constraint = [] |
|---|
| 699 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 700 | if k == 0: continue |
|---|
| 701 | constraint.append([k,G2varObj[j]]) |
|---|
| 702 | constraint += [shortmodelist[i],False,'f'] |
|---|
| 703 | self.Constraints.append(constraint) |
|---|
| 704 | #---------------------------------------------------------------------- |
|---|
| 705 | # save the ISODISTORT info for "mode analysis" |
|---|
| 706 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 707 | self.Phase['ISODISTORT'].update({ |
|---|
| 708 | 'IsoModeList' : modelist, |
|---|
| 709 | 'G2ModeList' : shortmodelist, |
|---|
| 710 | 'IsoVarList' : coordVarLbl, |
|---|
| 711 | 'G2VarList' : G2varObj, |
|---|
| 712 | 'ParentStructure' : ParentCoordinates, |
|---|
| 713 | 'Var2ModeMatrix' : displacivemodeInvmatrix, |
|---|
| 714 | 'Mode2VarMatrix' : displacivemodematrix, |
|---|
| 715 | }) |
|---|
| 716 | # make explaination dictionary |
|---|
| 717 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 718 | explaination[shortmode] = "ISODISTORT full name "+str(mode) |
|---|
| 719 | #---------------------------------------------------------------------- |
|---|
| 720 | # now read in the ISODISTORT occupancy modes |
|---|
| 721 | #---------------------------------------------------------------------- |
|---|
| 722 | if blk.get('_iso_occupancymode_label'): |
|---|
| 723 | modelist = [] |
|---|
| 724 | shortmodelist = [] |
|---|
| 725 | for lbl in blk.get('_iso_occupancymode_label'): |
|---|
| 726 | modelist.append(lbl) |
|---|
| 727 | # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB |
|---|
| 728 | # where SSSSS is the parent spacegroup, [x,y,z] is a location |
|---|
| 729 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl) |
|---|
| 730 | # this extracts the AAAAA and BBBBB parts of the string |
|---|
| 731 | if regexp: |
|---|
| 732 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 733 | lbl = lbl.replace('order','o') |
|---|
| 734 | G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list |
|---|
| 735 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 736 | occVarLbl = [] |
|---|
| 737 | G2varLbl = [] |
|---|
| 738 | G2varObj = [] |
|---|
| 739 | error = False |
|---|
| 740 | for lbl in blk.get('_iso_deltaoccupancy_label'): |
|---|
| 741 | occVarLbl.append(lbl) |
|---|
| 742 | if '_' in lbl: |
|---|
| 743 | albl = lbl[:lbl.rfind('_')] |
|---|
| 744 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 745 | else: |
|---|
| 746 | self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl |
|---|
| 747 | error = True |
|---|
| 748 | continue |
|---|
| 749 | if albl not in atomlbllist: |
|---|
| 750 | self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl |
|---|
| 751 | error = True |
|---|
| 752 | continue |
|---|
| 753 | else: |
|---|
| 754 | anum = atomlbllist.index(albl) |
|---|
| 755 | var = varLookup.get(vlbl) |
|---|
| 756 | if not var: |
|---|
| 757 | self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl |
|---|
| 758 | error = True |
|---|
| 759 | continue |
|---|
| 760 | G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID |
|---|
| 761 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 762 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 763 | )) |
|---|
| 764 | if error: |
|---|
| 765 | raise Exception("Error decoding variable labels") |
|---|
| 766 | |
|---|
| 767 | if len(G2varObj) != len(modelist): |
|---|
| 768 | print ("non-square input") |
|---|
| 769 | raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy") |
|---|
| 770 | |
|---|
| 771 | error = False |
|---|
| 772 | ParentCoordinates = {} |
|---|
| 773 | for lbl,exp in zip( |
|---|
| 774 | blk.get('_iso_occupancy_label'), |
|---|
| 775 | blk.get('_iso_occupancy_formula'), |
|---|
| 776 | ): |
|---|
| 777 | if '_' in lbl: |
|---|
| 778 | albl = lbl[:lbl.rfind('_')] |
|---|
| 779 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 780 | else: |
|---|
| 781 | self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl |
|---|
| 782 | error = True |
|---|
| 783 | continue |
|---|
| 784 | if vlbl != 'occ': |
|---|
| 785 | self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl |
|---|
| 786 | error = True |
|---|
| 787 | continue |
|---|
| 788 | if '+' in exp: |
|---|
| 789 | val = exp.split('+')[0].strip() |
|---|
| 790 | val = G2p3.FormulaEval(val) |
|---|
| 791 | if val is None: |
|---|
| 792 | self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl |
|---|
| 793 | error = True |
|---|
| 794 | continue |
|---|
| 795 | ParentCoordinates[albl] = val |
|---|
| 796 | if error: |
|---|
| 797 | raise Exception("Error decoding occupancy labels") |
|---|
| 798 | # get mapping of modes to atomic coordinate displacements |
|---|
| 799 | occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 800 | for row,col,val in zip( |
|---|
| 801 | blk.get('_iso_occupancymodematrix_row'), |
|---|
| 802 | blk.get('_iso_occupancymodematrix_col'), |
|---|
| 803 | blk.get('_iso_occupancymodematrix_value'),): |
|---|
| 804 | occupancymodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 805 | # Invert to get mapping of atom displacements to modes |
|---|
| 806 | occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix) |
|---|
| 807 | # create the constraints |
|---|
| 808 | for i,row in enumerate(occupancymodeInvmatrix): |
|---|
| 809 | constraint = [] |
|---|
| 810 | for j,(lbl,k) in enumerate(zip(occVarLbl,row)): |
|---|
| 811 | if k == 0: continue |
|---|
| 812 | constraint.append([k,G2varObj[j]]) |
|---|
| 813 | constraint += [shortmodelist[i],False,'f'] |
|---|
| 814 | self.Constraints.append(constraint) |
|---|
| 815 | #---------------------------------------------------------------------- |
|---|
| 816 | # save the ISODISTORT info for "mode analysis" |
|---|
| 817 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 818 | self.Phase['ISODISTORT'].update({ |
|---|
| 819 | 'OccModeList' : modelist, |
|---|
| 820 | 'G2OccModeList' : shortmodelist, |
|---|
| 821 | 'OccVarList' : occVarLbl, |
|---|
| 822 | 'G2OccVarList' : G2varObj, |
|---|
| 823 | 'BaseOcc' : ParentCoordinates, |
|---|
| 824 | 'Var2OccMatrix' : occupancymodeInvmatrix, |
|---|
| 825 | 'Occ2VarMatrix' : occupancymodematrix, |
|---|
| 826 | }) |
|---|
| 827 | # make explaination dictionary |
|---|
| 828 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 829 | explaination[shortmode] = "ISODISTORT full name "+str(mode) |
|---|
| 830 | #---------------------------------------------------------------------- |
|---|
| 831 | # done with read |
|---|
| 832 | #---------------------------------------------------------------------- |
|---|
| 833 | if explaination: self.Constraints.append(explaination) |
|---|
| 834 | |
|---|
| 835 | # # debug: show the mode var to mode relations |
|---|
| 836 | # for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 837 | # l = '' |
|---|
| 838 | # for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 839 | # if k == 0: continue |
|---|
| 840 | # if l: l += ' + ' |
|---|
| 841 | # #l += lbl+' * '+str(k) |
|---|
| 842 | # l += G2varLbl[j]+' * '+str(k) |
|---|
| 843 | # print str(i) + ': '+shortmodelist[i]+' = '+l |
|---|
| 844 | # print 70*'=' |
|---|
| 845 | |
|---|
| 846 | # # debug: Get the ISODISTORT offset values |
|---|
| 847 | # coordVarDelta = {} |
|---|
| 848 | # for lbl,val in zip( |
|---|
| 849 | # blk.get('_iso_deltacoordinate_label'), |
|---|
| 850 | # blk.get('_iso_deltacoordinate_value'),): |
|---|
| 851 | # coordVarDelta[lbl] = float(val) |
|---|
| 852 | # modeVarDelta = {} |
|---|
| 853 | # for lbl,val in zip( |
|---|
| 854 | # blk.get('_iso_displacivemode_label'), |
|---|
| 855 | # blk.get('_iso_displacivemode_value'),): |
|---|
| 856 | # modeVarDelta[lbl] = cif.get_number_with_esd(val)[0] |
|---|
| 857 | |
|---|
| 858 | # print 70*'=' |
|---|
| 859 | # # compute the mode values from the reported coordinate deltas |
|---|
| 860 | # for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 861 | # l = '' |
|---|
| 862 | # sl = '' |
|---|
| 863 | # s = 0. |
|---|
| 864 | # for lbl,k in zip(coordVarLbl,row): |
|---|
| 865 | # if k == 0: continue |
|---|
| 866 | # if l: l += ' + ' |
|---|
| 867 | # l += lbl+' * '+str(k) |
|---|
| 868 | # if sl: sl += ' + ' |
|---|
| 869 | # sl += str(coordVarDelta[lbl])+' * '+str(k) |
|---|
| 870 | # s += coordVarDelta[lbl] * k |
|---|
| 871 | # print 'a'+str(i)+' = '+l |
|---|
| 872 | # print '\t= '+sl |
|---|
| 873 | # print modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s |
|---|
| 874 | # print |
|---|
| 875 | |
|---|
| 876 | # print 70*'=' |
|---|
| 877 | # # compute the coordinate displacements from the reported mode values |
|---|
| 878 | # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix): |
|---|
| 879 | # l = '' |
|---|
| 880 | # sl = '' |
|---|
| 881 | # s = 0.0 |
|---|
| 882 | # for j,k in enumerate(row): |
|---|
| 883 | # if k == 0: continue |
|---|
| 884 | # if l: l += ' + ' |
|---|
| 885 | # l += 'a'+str(j+1)+' * '+str(k) |
|---|
| 886 | # if sl: sl += ' + ' |
|---|
| 887 | # sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k) |
|---|
| 888 | # s += modeVarDelta[modelist[j]] * k |
|---|
| 889 | # print lbl+' = '+l |
|---|
| 890 | # print '\t= '+sl |
|---|
| 891 | # print lbl,G2varLbl[i],coordVarDelta[lbl],s |
|---|
| 892 | # print |
|---|
| 893 | |
|---|
| 894 | # determine the coordinate delta values from deviations from the parent structure |
|---|
| 895 | # for atmline in self.Phase['Atoms']: |
|---|
| 896 | # lbl = atmline[0] |
|---|
| 897 | # x,y,z = atmline[3:6] |
|---|
| 898 | # if lbl not in ParentCoordinates: |
|---|
| 899 | # print lbl,x,y,z |
|---|
| 900 | # continue |
|---|
| 901 | # px,py,pz = ParentCoordinates[lbl] |
|---|
| 902 | # print lbl,x,y,z,x-px,y-py,z-pz |
|---|
