source: trunk/imports/G2phase_CIF.py

Last change on this file was 5170, checked in by toby, 3 months ago

deal with ISODISP _iso_coordinate_formula with - rather than +

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