source: trunk/imports/G2phase_CIF.py @ 4635

Last change on this file since 4635 was 4635, checked in by vondreele, 13 months ago

change all 'Ur' to 'r' in open commands
double all '\' in G2scriptable
close dangling open files

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