source: trunk/imports/G2phase_CIF.py @ 5065

Last change on this file since 5065 was 5065, checked in by vondreele, 8 months ago

remove check on Histograms & Phases in beginning of UpdateConstraints? - had prevented viewing constraints when no phases
Always enable SHOWISO - same reason
remove check on length of varyList in LoadParmDict? - prevented export of cif when there was no refinement
checks on failed ISODISTORT runs - now shows resulting html page with error message
fixes to RMCProfile GUI startup
fix distortion mode plotting - now moves bonds & polyhedra
cif importer looks for space group number if symbol not interpretable - i.e. full HM symbol
remove a check on k==0 in line 993 of cif importer; k is always zero in ISODISTORT cifs when there in no preset mode distortion value
fix name bug in line 1008 & remove breakpoint in cif importer
ISODISTORT makes html page & displays it in case of ISODISTORT error

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 60.4 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-11-04 13:09:15 +0000 (Thu, 04 Nov 2021) $
4# $Author: vondreele $
5# $Revision: 5065 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 5065 2021-11-04 13:09:15Z 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 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: 5065 $")
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('_','')
237                    self.Phase['General']['Type'] = 'nuclear'
238#process space group symbol
239            E,SGData = G2spc.SpcGroup(SpGrp)
240            if E and SpGrp:
241                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
242                if SpGrpNorm:
243                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
244            if E:   #try lookup from number  - found full symbol?               
245                SpGrpNorm = G2spc.spgbyNum[int(blk.get('_symmetry_Int_Tables_number'))]
246                if SpGrpNorm:
247                    E,SGData = G2spc.SpcGroup(SpGrpNorm) 
248            # nope, try the space group "out of the Box"
249            if E:
250                self.warnings += 'ERROR in space group symbol '+SpGrp
251                self.warnings += '\nThe space group has been set to "P 1". '
252                self.warnings += "Change this in phase's General tab."
253                self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
254                self.warnings += G2spc.SGErrors(E)
255                SGData = G2obj.P1SGData # P 1
256            self.Phase['General']['SGData'] = SGData
257            # save symmetry operators, if specified (use to check for origin 1 in GSASIIdataGUI OnImportPhase)
258            ops = blk.get("_symmetry_equiv_pos_as_xyz")   # try older name 1st
259            if ops:
260                self.SymOps['xyz'] = ops
261            elif blk.get("_space_group_symop_operation_xyz"):
262                self.SymOps['xyz'] = blk.get("_space_group_symop_operation_xyz")
263            else:
264                if 'xyz' in self.SymOps: del self.SymOps['xyz']
265            if Super:
266                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
267                if E:
268                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
269                    self.warnings += '\n'+E
270                    SuperSg = SuperSg[:SuperSg.index(')')+1]
271                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
272                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
273                self.Phase['General']['SSGData'] = SSGData
274                if magnetic:
275                    self.MPhase['General']['SGData'] = SGData
276                    self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
277                    self.MPhase['General']['SSGData'] = SSGData
278
279            if magnetic:    #replace std operaors with those from cif file - probably not the same!
280                SGData['SGFixed'] = True
281                SGData['SGOps'] = []
282                SGData['SGCen'] = []
283                if Super:
284                    SSGData['SSGOps'] = []
285                    SSGData['SSGCen'] = []
286                    try:
287                        sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id')
288                        sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id')
289                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1]
290                        centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1]                   
291                    except KeyError:        #old mag cif names
292                        sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id')
293                        sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id')
294                        opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1]
295                        centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1]
296                    spnflp = []
297                    for op in sgoploop:
298                        M,T,S = G2spc.MagSSText2MTS(op[opid])
299                        SSGData['SSGOps'].append([np.array(M,dtype=float),T])
300                        SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]])
301                        spnflp.append(S)
302                    censpn = []
303                    for cent in sgcenloop:
304                        M,C,S = G2spc.MagSSText2MTS(cent[centid])
305                        SSGData['SSGCen'].append(C)
306                        SGData['SGCen'].append(C[:3])
307                        censpn += list(np.array(spnflp)*S)
308                    self.MPhase['General']['SSGData'] = SSGData
309                else:
310                    try:
311                        sgoploop = blk.GetLoop('_space_group_symop_magn_operation.id')
312                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
313                        try:
314                            sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
315                            centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]
316                        except KeyError:
317                            sgcenloop = None
318                    except KeyError:                                         
319                        try:
320                            sgoploop = blk.GetLoop('_space_group_symop_magn.id')
321                            sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
322                            opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
323                            centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]                   
324                        except KeyError:        #old mag cif names
325                            sgoploop = blk.GetLoop('_space_group_symop.magn_id')
326                            sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
327                            opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
328                            centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
329                    spnflp = []
330                    for op in sgoploop:
331                        try:
332                            M,T,S = G2spc.MagText2MTS(op[opid])
333                            SGData['SGOps'].append([np.array(M,dtype=float),T])
334                            spnflp.append(S)
335                        except KeyError:
336                            self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II'
337                            return False
338                    censpn = []
339                    if sgcenloop:
340                        for cent in sgcenloop:
341                            M,C,S = G2spc.MagText2MTS(cent[centid])
342                            SGData['SGCen'].append(C)
343                            censpn += list(np.array(spnflp)*S)
344                    else:
345                            M,C,S = G2spc.MagText2MTS('x,y,z,+1')
346                            SGData['SGCen'].append(C)
347                            censpn += list(np.array(spnflp)*S)                       
348                self.MPhase['General']['SGData'] = SGData
349                self.MPhase['General']['SGData']['SpnFlp'] = censpn
350                G2spc.GenMagOps(SGData)         #set magMom
351                self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
352                if "1'" in MSpGrp:
353                    SGData['SGGray'] = True
354                MagPtGp = blk.get('_space_group.magn_point_group')
355                if not MagPtGp:
356                    MagPtGp = blk.get('_space_group.magn_point_group_name')
357                if not MagPtGp:
358                    MagPtGp = blk.get('_space_group_magn.point_group_name')
359                self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
360
361            # cell parameters
362            cell = []
363            for lbl in cellitems:
364                cell.append(cif.get_number_with_esd(blk[lbl])[0])
365            Volume = G2lat.calc_V(G2lat.cell2A(cell))
366            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
367            if magnetic:
368                self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]               
369            if Super:
370                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
371                waveDict = dict(waveloop.items())
372                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
373                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
374                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
375
376            # read in atoms
377            self.errors = 'Error during reading of atoms'
378            atomlbllist = [] # table to look up atom IDs
379            atomloop = blk.GetLoop('_atom_site_label')
380            atomkeys = [i.lower() for i in atomloop.keys()]
381            if not blk.get('_atom_site_type_symbol'):
382                isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
383            if magnetic:
384                try:
385                    magmoment = '_atom_site_moment.label'
386                    magatomloop = blk.GetLoop(magmoment)
387                    magatomkeys = [i.lower() for i in magatomloop.keys()]
388                    magatomlabels = blk.get(magmoment)
389                    G2MagDict = {'_atom_site_moment.label': 0,
390                                 '_atom_site_moment.crystalaxis_x':7,
391                                 '_atom_site_moment.crystalaxis_y':8,
392                                 '_atom_site_moment.crystalaxis_z':9}
393                except KeyError:
394                    magmoment = '_atom_site_moment_label'
395                    magatomloop = blk.GetLoop(magmoment)
396                    magatomkeys = [i.lower() for i in magatomloop.keys()]
397                    magatomlabels = blk.get(magmoment)
398                    G2MagDict = {'_atom_site_moment_label': 0,
399                                 '_atom_site_moment_crystalaxis_x':7,
400                                 '_atom_site_moment_crystalaxis_y':8,
401                                 '_atom_site_moment_crystalaxis_z':9}
402                   
403            if blk.get('_atom_site_aniso_label'):
404                anisoloop = blk.GetLoop('_atom_site_aniso_label')
405                anisokeys = [i.lower() for i in anisoloop.keys()]
406                anisolabels = blk.get('_atom_site_aniso_label')
407            else:
408                anisoloop = None
409                anisokeys = []
410                anisolabels = []
411            if Super:
412                occFloop = None
413                occCloop = None
414                occFdict = {}
415                occCdict = {}
416                displSloop = None
417                displFloop = None
418                MagFloop = None
419                displSdict = {}
420                displFdict = {}
421                MagFdict = {}
422                UijFloop = None
423                UijFdict = {}
424                #occupancy modulation
425                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
426                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
427                    occFdict = dict(occFloop.items())
428                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
429                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
430                    occCdict = dict(occCloop.items())
431                #position modulation
432                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
433                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
434                    displFdict = dict(displFloop.items())                           
435                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
436                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
437                    displSdict = dict(displSloop.items())
438                #U modulation
439                if blk.get('_atom_site_U_Fourier_atom_site_label'):
440                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
441                    UijFdict = dict(UijFloop.items())
442                #Mag moment modulation
443                if blk.get('_atom_site_moment_Fourier_atom_site_label'):
444                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
445                    MagFdict = dict(MagFloop.items())
446                    Mnames =  ['_atom_site_moment_fourier_atom_site_label',
447                               '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
448                               '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
449                elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
450                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
451                    MagFdict = dict(MagFloop.items())                           
452                    Mnames =  ['_atom_site_moment_fourier.atom_site_label',
453                               '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
454                               '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
455            self.Phase['Atoms'] = []
456            if magnetic:
457                self.MPhase['Atoms'] = []
458            G2AtomDict = {  '_atom_site_type_symbol' : 1,
459                            '_atom_site_label' : 0,
460                            '_atom_site_fract_x' : 3,
461                            '_atom_site_fract_y' : 4,
462                            '_atom_site_fract_z' : 5,
463                            '_atom_site_occupancy' : 6,
464                            '_atom_site_aniso_u_11' : 11,
465                            '_atom_site_aniso_u_22' : 12,
466                            '_atom_site_aniso_u_33' : 13,
467                            '_atom_site_aniso_u_12' : 14,
468                            '_atom_site_aniso_u_13' : 15,
469                            '_atom_site_aniso_u_23' : 16, }
470
471            ranIdlookup = {}
472            for aitem in atomloop:
473                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
474                for val,key in zip(aitem,atomkeys):
475                    col = G2AtomDict.get(key,-1)
476                    if col >= 3:
477                        atomlist[col] = cif.get_number_with_esd(val)[0]
478                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
479                    elif col is not None and col != -1:
480                        atomlist[col] = val
481                    elif key in ('_atom_site_thermal_displace_type',
482                               '_atom_site_adp_type'):   #Iso or Aniso?
483                        if val.lower() == 'uani':
484                            atomlist[9] = 'A'
485                    elif key == '_atom_site_u_iso_or_equiv':
486                        uisoval = cif.get_number_with_esd(val)[0]
487                        if uisoval is not None: 
488                            atomlist[10] = uisoval
489                    elif key == '_atom_site_b_iso_or_equiv':
490                        uisoval = cif.get_number_with_esd(val)[0]
491                        if uisoval is not None: 
492                            atomlist[10] = uisoval/(8*np.pi**2)
493                if not atomlist[1] and atomlist[0]:
494                    typ = atomlist[0].rstrip('0123456789-+')
495                    if G2elem.CheckElement(typ):
496                        atomlist[1] = typ
497                    if not atomlist[1]:
498                        atomlist[1] = 'Xe'
499                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
500                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
501                    atomlist[9] = 'A'
502                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
503                        col = G2AtomDict.get(key)
504                        if col:
505                            atomlist[col] = cif.get_number_with_esd(val)[0]
506                if None in atomlist[11:17]:
507                    atomlist[9] = 'I'
508                    atomlist[11:17] =  [0.,0.,0.,0.,0.,0.]
509                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
510                atomlist[1] = G2elem.FixValence(atomlist[1])
511                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
512                self.Phase['Atoms'].append(atomlist)
513                ranIdlookup[atomlist[0]] = atomlist[-1]
514                if atomlist[0] in atomlbllist:
515                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
516                else:
517                    atomlbllist.append(atomlist[0])
518
519                if magnetic and atomlist[0] in magatomlabels:
520                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
521                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
522                        mcol = G2MagDict.get(mkey,-1)
523                        if mcol > 0:
524                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
525                    self.MPhase['Atoms'].append(matomlist)
526                if Super:
527                    Sfrac = np.zeros((4,2))
528                    Sadp = np.zeros((4,12))
529                    Spos = np.zeros((4,6))
530                    Smag = np.zeros((4,6))
531                    nim = -1
532                    waveType = 'Fourier'
533                    if  occCdict:
534                        for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']):
535                            if item == atomlist[0]:
536                                waveType = 'Crenel'
537                                val = occCdict['_atom_site_occ_special_func_crenel_c'][i]
538                                Sfrac[0][0] = cif.get_number_with_esd(val)[0]
539                                val = occCdict['_atom_site_occ_special_func_crenel_w'][i]
540                                Sfrac[0][1] = cif.get_number_with_esd(val)[0]
541                                nim = 1
542                   
543                    if nim >= 0:
544                        Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]]
545                    else:
546                        Sfrac = []
547                    nim = -1
548                    if displFdict:
549                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
550                            if item == atomlist[0]:
551                                waveType = 'Fourier'                               
552                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
553                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
554                                if im != nim:
555                                    nim = im
556                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
557                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
558                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
559                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
560                    if nim >= 0:
561                        Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]]
562                    else:
563                        Spos = []
564                    nim = -1
565                    if UijFdict:
566                        nim = -1
567                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
568                            if item == atomlist[0]:
569                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
570                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
571                                if im != nim:
572                                    nim = im
573                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
574                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
575                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
576                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
577                    if nim >= 0:
578                        Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]]
579                    else:
580                        Sadp = []
581                    nim = -1
582                    if MagFdict:
583                        nim = -1
584                        for i,item in enumerate(MagFdict[Mnames[0]]):
585                            if item == atomlist[0]:
586                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
587                                im = int(MagFdict[Mnames[2]][i])
588                                if im != nim:
589                                    nim = im
590                                val = MagFdict[Mnames[3]][i]
591                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
592                                val = MagFdict[Mnames[4]][i]
593                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
594                    if nim >= 0:
595                        Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]]
596                    else:
597                        Smag = []
598                    SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
599                    if magnetic and atomlist[0] in magatomlabels:
600                        matomlist.append(SSdict)
601                    atomlist.append(SSdict)
602            if len(atomlbllist) != len(self.Phase['Atoms']):
603                isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
604            for lbl in phasenamefields: # get a name for the phase
605                try:
606                    name = blk.get(lbl)
607                except:
608                    continue
609                if name is None:
610                    continue
611                name = name.strip()
612                if name == '?' or name == '.':
613                    continue
614                else:
615                    break
616            else: # no name found, use block name for lack of a better choice; for isodistort use filename
617                name = blknm
618            if 'isodistort' in name:
619                self.Phase['General']['Name'] = os.path.split(filename)[1].split('.')[0]
620            else:
621                self.Phase['General']['Name'] = name.strip()
622            self.Phase['General']['Super'] = Super
623            self.Phase = copy.deepcopy(self.Phase)      #clean copy
624            if magnetic:
625                self.MPhase = copy.deepcopy(self.MPhase)    #clean copy
626                self.MPhase['General']['Type'] = 'magnetic'               
627                self.MPhase['General']['Name'] = name.strip()+' mag'
628                self.MPhase['General']['Super'] = Super
629                if Super:
630                    self.MPhase['General']['Modulated'] = True
631                    self.MPhase['General']['SuperVec'] = SuperVec
632                    self.MPhase['General']['SuperSg'] = SuperSg
633                    if self.MPhase['General']['SGData']['SGGray']:
634                        self.MPhase['General']['SuperSg'] += 's'
635                if 'mcif' not in filename:
636                    self.Phase = copy.deepcopy(self.MPhase)
637                    del self.MPhase
638            else:
639                self.MPhase = None
640            if Super:
641                if self.Phase['General']['SGData']['SGGray']:
642                    SGData = self.Phase['General']['SGData']
643                    SGData['SGGray'] = False
644                    ncen = len(SGData['SGCen'])
645                    SGData['SGCen'] = SGData['SGCen'][:ncen//2]
646                    self.Phase['General']['SGData'].update(SGData)
647                self.Phase['General']['Modulated'] = True
648                self.Phase['General']['SuperVec'] = SuperVec
649                self.Phase['General']['SuperSg'] = SuperSg
650                if not self.Phase['General']['SGData']['SGFixed']:
651                    self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
652            if self.ISODISTORT_test(blk):
653                if isodistort_warnings:
654                    self.warnings += isodistort_warnings
655                else:
656                    self.errors = "Error while processing ISODISTORT constraints"
657                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup,filename)
658                    self.errors = ""
659            returnstat = True
660        return returnstat
661       
662    def ISODISTORT_test(self,blk):
663        '''Test if there is any ISODISTORT information in CIF
664
665        At present only _iso_displacivemode... and _iso_occupancymode... are
666        tested.
667        '''
668        for i in ('_iso_displacivemode_label',
669                  '_iso_occupancymode_label'):
670            if blk.get(i): return True
671        return False
672       
673    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup,filename):
674        '''Process ISODISTORT items to create constraints etc.
675        Constraints are generated from information extracted from
676        loops beginning with _iso_ and are placed into
677        self.Constraints, which contains a list of
678        :ref:`constraints tree items <Constraint_definitions_table>`
679        and one dict.
680        The dict contains help text for each generated ISODISTORT variable
681
682        At present only _iso_displacivemode... and _iso_occupancymode... are
683        processed. Not yet processed: _iso_magneticmode...,
684        _iso_rotationalmode... & _iso_strainmode...
685        '''
686        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
687        'Maps ISODISTORT parm names to GSAS-II names'
688        # used for all types of modes
689        self.Constraints = []
690        explaination = {}
691        G2obj.AddPhase2Index(self,filename)   # put phase info into Var index
692        #----------------------------------------------------------------------
693        # read in the ISODISTORT displacement modes
694        #----------------------------------------------------------------------
695        if blk.get('_iso_displacivemode_label'):
696            modelist = []
697            shortmodelist = []
698            idlist = []
699            for id,lbl in zip(
700                blk.get('_iso_displacivemode_ID'),
701                blk.get('_iso_displacivemode_label')):
702                idlist.append(int(id))
703                modelist.append(lbl)
704                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
705            # just in case the items are not ordered increasing by id, sort them here
706            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
707            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
708            # read in the coordinate offset variables names and map them to G2 names/objects
709            coordVarLbl = []
710            G2varObj = []
711            idlist = []
712            error = False
713            for id,lbl in zip(
714                blk.get('_iso_deltacoordinate_ID'),
715                blk.get('_iso_deltacoordinate_label') ):
716                idlist.append(int(id))
717                coordVarLbl.append(lbl)
718                if '_' in lbl:
719                    albl = lbl[:lbl.rfind('_')]
720                    vlbl = lbl[lbl.rfind('_')+1:]
721                else:
722                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
723                    error = True
724                    continue
725                if albl not in atomlbllist:
726                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
727                    error = True
728                    continue
729                # else:
730                #     anum = atomlbllist.index(albl)
731                var = varLookup.get(vlbl)
732                if not var:
733                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
734                    error = True
735                    continue
736                G2varObj.append(G2obj.G2VarObj(
737                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
738                    ))
739            if error:
740                raise Exception("Error decoding variable labels")
741            # just in case the items are not ordered increasing by id, sort them here
742            coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
743            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
744
745            if len(G2varObj) != len(modelist):
746                print ("non-square input")
747                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
748
749            error = False
750            ParentCoordinates = {}
751            for lbl,exp in zip(
752                blk.get('_iso_coordinate_label'),
753                blk.get('_iso_coordinate_formula') ):
754                if '_' in lbl:
755                    albl = lbl[:lbl.rfind('_')]
756                    vlbl = lbl[lbl.rfind('_')+1:]
757                else:
758                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
759                    error = True
760                    continue
761                if vlbl not in 'xyz' or len(vlbl) != 1:
762                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
763                    error = True
764                    continue
765                i = 'xyz'.index(vlbl)
766                if not ParentCoordinates.get(albl):
767                    ParentCoordinates[albl] = [None,None,None]
768                if '+' in exp:
769                    val = exp.split('+')[0].strip()
770                    val = G2p3.FormulaEval(val)
771                    if val is None:
772                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
773                        error = True
774                        continue
775                    ParentCoordinates[albl][i] = val
776                else:
777                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
778            if error:
779                print (self.warnings)
780                raise Exception("Error decoding variable labels")
781            # get mapping of modes to atomic coordinate displacements
782            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
783            for row,col,val in zip(
784                blk.get('_iso_displacivemodematrix_row'),
785                blk.get('_iso_displacivemodematrix_col'),
786                blk.get('_iso_displacivemodematrix_value'),):
787                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
788            # Invert to get mapping of atom displacements to modes
789            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
790            # create the constraints
791            modeVarList = []
792            modeDispl = []
793            for i,row in enumerate(displacivemodeInvmatrix):
794                constraint = []
795                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
796                    if k == 0: continue
797                    constraint.append([k,G2varObj[j]])
798                modeVar = G2obj.G2VarObj(
799                    (self.Phase['ranId'],None,shortmodelist[i],None))
800                modeVarList.append(modeVar)               
801                constraint += [modeVar,False,'f']
802                self.Constraints.append(constraint)
803                modeDispl.append(0.0)
804            # normilization constants
805            normlist = []
806            idlist = []
807            for id,exp in zip(
808                blk.get('_iso_displacivemodenorm_ID'),
809                blk.get('_iso_displacivemodenorm_value'),
810                ):
811                idlist.append(int(id))
812                normlist.append(float(exp))
813            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
814            #----------------------------------------------------------------------
815            # save the ISODISTORT info for "mode analysis"
816            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
817            self.Phase['ISODISTORT'].update({
818                # coordinate items
819                'IsoVarList' : coordVarLbl,
820                'G2VarList' : G2varObj,
821                'ParentStructure' : ParentCoordinates,
822                # mode items
823                'IsoModeList' : modelist,
824                'G2ModeList' : modeVarList,
825                'NormList' : normlist,
826                # transform matrices
827                'Var2ModeMatrix' : displacivemodeInvmatrix,
828                'Mode2VarMatrix' : displacivemodematrix,
829                'modeDispl' : modeDispl
830                })
831            # make explaination dictionary
832            for mode,shortmode in zip(modelist,shortmodelist):
833                explaination[str(self.Phase['ranId'])+shortmode] = (
834                    "ISODISTORT full name "+str(mode))
835        #----------------------------------------------------------------------
836        # now read in the ISODISTORT occupancy modes
837        #----------------------------------------------------------------------
838        if blk.get('_iso_occupancymode_label'):
839            modelist = []
840            shortmodelist = []
841            idlist = []
842            for id,lbl in zip(
843                blk.get('_iso_occupancymode_ID'),
844                blk.get('_iso_occupancymode_label')):
845                idlist.append(int(id))
846                modelist.append(lbl)
847                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
848            # just in case the items are not ordered increasing by id, sort them here
849            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
850            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
851            # read in the coordinate offset variables names and map them to G2 names/objects
852            occVarLbl = []
853            G2varObj = []
854            idlist = []
855            error = False
856            for id,lbl in zip(
857                blk.get('_iso_deltaoccupancy_ID'),
858                blk.get('_iso_deltaoccupancy_label') ):
859                idlist.append(int(id))
860                occVarLbl.append(lbl)
861                if '_' in lbl:
862                    albl = lbl[:lbl.rfind('_')]
863                    vlbl = lbl[lbl.rfind('_')+1:]
864                else:
865                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
866                    error = True
867                    continue
868                if albl not in atomlbllist:
869                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
870                    error = True
871                    continue
872                # else:
873                #     anum = atomlbllist.index(albl)
874                var = varLookup.get(vlbl)
875                if not var:
876                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
877                    error = True
878                    continue
879                G2varObj.append(G2obj.G2VarObj(
880                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
881                    ))
882            if error:
883                raise Exception("Error decoding variable labels")
884            # just in case the items are not ordered increasing by id, sort them here
885            occVarLbl = [i for i,j in sorted(zip(occVarLbl,idlist),key=lambda k:k[1])]
886            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
887
888            if len(G2varObj) != len(modelist):
889                print ("non-square input")
890                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
891
892            error = False
893            ParentOcc = {}
894            for lbl,exp in zip(
895                blk.get('_iso_occupancy_label'),
896                blk.get('_iso_occupancy_formula') ):
897                if '_' in lbl:
898                    albl = lbl[:lbl.rfind('_')]
899                    vlbl = lbl[lbl.rfind('_')+1:]
900                else:
901                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
902                    error = True
903                    continue
904                if vlbl != 'occ':
905                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
906                    error = True
907                    continue
908                if '+' in exp:
909                    val = exp.split('+')[0].strip()
910                    val = G2p3.FormulaEval(val)
911                    if val is None:
912                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
913                        error = True
914                        continue
915                    ParentOcc[albl] = val
916            if error:
917                raise Exception("Error decoding occupancy labels")
918            # get mapping of modes to atomic coordinate displacements
919            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
920            for row,col,val in zip(
921                blk.get('_iso_occupancymodematrix_row'),
922                blk.get('_iso_occupancymodematrix_col'),
923                blk.get('_iso_occupancymodematrix_value'),):
924                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
925            # Invert to get mapping of atom displacements to modes
926            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
927            # create the constraints
928            modeVarList = []
929            for i,row in enumerate(occupancymodeInvmatrix):
930                constraint = []
931                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
932                    if k == 0: continue
933                    constraint.append([k,G2varObj[j]])
934                modeVar = G2obj.G2VarObj(
935                    (self.Phase['ranId'],None,shortmodelist[i],None))
936                modeVarList.append(modeVar)               
937                constraint += [modeVar,False,'f']
938                self.Constraints.append(constraint)
939            # normilization constants
940            normlist = []
941            idlist = []
942            for id,exp in zip(
943                blk.get('_iso_occupancymodenorm_ID'),
944                blk.get('_iso_occupancymodenorm_value'),
945                ):
946                idlist.append(int(id))
947                normlist.append(float(exp))
948            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
949            #----------------------------------------------------------------------
950            # save the ISODISTORT info for "mode analysis"
951            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
952            self.Phase['ISODISTORT'].update({
953                # coordinate items
954                'OccVarList' : occVarLbl,
955                'G2OccVarList' : G2varObj,
956                'BaseOcc' : ParentOcc,
957                # mode items
958                'OccModeList' : modelist,
959                'G2OccModeList' : modeVarList,
960                'OccNormList' : normlist,
961                # transform matrices
962                'Var2OccMatrix' : occupancymodeInvmatrix,
963                'Occ2VarMatrix' : occupancymodematrix,
964                })
965            # make explaination dictionary
966            for mode,shortmode in zip(modelist,shortmodelist):
967                explaination[str(self.Phase['ranId'])+shortmode] = (
968                    "ISODISTORT full name "+str(mode))
969        #----------------------------------------------------------------------
970        # done with read
971        #----------------------------------------------------------------------
972        if explaination: self.Constraints.append(explaination)
973
974        if not debug: return
975
976        # debug: show displacive mode var to mode relations
977        if 'IsoVarList' in self.Phase['ISODISTORT']:
978            # coordinate items
979            #coordVarLbl = self.Phase['ISODISTORT']['IsoVarList']
980            G2varObj = self.Phase['ISODISTORT']['G2VarList']
981            #ParentCoordinates = self.Phase['ISODISTORT']['ParentStructure']
982            # mode items
983            modelist = self.Phase['ISODISTORT']['IsoModeList']
984            modeVarList = self.Phase['ISODISTORT']['G2ModeList']
985            normlist = self.Phase['ISODISTORT']['NormList']
986            # transform matrices
987            #displacivemodeInvmatrix = self.Phase['ISODISTORT']['Var2ModeMatrix']
988            #displacivemodematrix = self.Phase['ISODISTORT']['Mode2VarMatrix']
989       
990            print( 70*'=')
991            print('\nVar2ModeMatrix' ,'IsoVarList' )
992            def fmtConstr(i,head,l,var,k):
993                if len(head) + len(l) > 65:
994                    print(head+l)
995                    head = 20*' '
996                    l = ''
997#                if k == 0: return
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            for i,row in enumerate(displacivemodeInvmatrix):
1009                head = str(i) + ': ' + str(modeVarList[i]) + ' = '
1010                line = ''
1011                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
1012                    head,line = fmtConstr(j,head,line,G2varObj[j].name,k)
1013                print(head+line)
1014
1015            # Get the ISODISTORT offset values
1016            coordVarDelta = {}
1017            for lbl,val in zip(
1018                blk.get('_iso_deltacoordinate_label'),
1019                blk.get('_iso_deltacoordinate_value'),):
1020                coordVarDelta[lbl] = float(val)
1021            modeVarDelta = {}
1022            for lbl,val in zip(
1023                blk.get('_iso_displacivemode_label'),
1024                blk.get('_iso_displacivemode_value'),):
1025                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
1026
1027            print( 70*'=')
1028            print('\nInverse relations using Var2ModeMatrix, NormList, IsoVarList')
1029            # compute the mode values from the reported coordinate deltas
1030            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
1031                line = ''
1032                head = 'a'+str(i)+': '+str(modeVarList[i])+' = ('
1033                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
1034                    head,line = fmtConstr(j,head,line,lbl,k)
1035                print(head+line+') / '+('%.3f'%n))
1036            print('\nCalculation checks\n')
1037            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
1038                #l = ''
1039                sl = ''
1040                s = 0.
1041                for lbl,k in zip(coordVarLbl,row):
1042                    if k == 0: continue
1043                    #if l: l += ' + '
1044                    #l += lbl+' * '+str(k)
1045                    if sl: sl += ' + '
1046                    sl += str(coordVarDelta[lbl])+' * '+str(k)
1047                    s += coordVarDelta[lbl] * k
1048                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1049                print(' ?= ',modeVarDelta[modelist[i]])
1050                print()
1051
1052            print( 70*'=')
1053            print('\nDirect relations using Mode2VarMatrix, NormList, IsoVarList')
1054            # compute the coordinate displacements from the reported mode values
1055            DeltaCoords = {}
1056            for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
1057                l = ''
1058                s = 0.0
1059                at,d = lbl.rsplit('_',1)
1060                if at not in DeltaCoords:
1061                    DeltaCoords[at] = [0,0,0]
1062                for j,(k,n) in enumerate(zip(row,normlist)):
1063                    if k == 0: continue
1064                    if l: l += ' + '
1065                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1066                    s += n * modeVarDelta[modelist[j]] * k
1067                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1068                if d == 'dx':
1069                    DeltaCoords[at][0] = s
1070                elif d == 'dy':
1071                    DeltaCoords[at][1] = s
1072                elif d == 'dz':
1073                    DeltaCoords[at][2] = s
1074                else:
1075                    print('unexpected',d)
1076
1077            print( 70*'=')
1078            print('\nCoordinate checks')
1079            print('\nComputed')
1080            for at in sorted(DeltaCoords):
1081                s = at
1082                for i in range(3):
1083                    s += '  '
1084                    s += str(ParentCoordinates[at][i]+DeltaCoords[at][i])
1085                print(s)
1086
1087            # determine the coordinate delta values from deviations from the parent structure
1088            print('\nFrom CIF')
1089            for atmline in self.Phase['Atoms']:
1090                lbl = atmline[0]
1091                x,y,z = atmline[3:6]
1092                print( lbl,x,y,z)
1093
1094            print( 70*'=')
1095            print('\nGenerated constraints')
1096            for i in self.Constraints:
1097                if type(i) is dict:
1098                    print('\nconstraint help dict')
1099                    for key in i:
1100                        print('\t',key,':',i[key])
1101                elif i[-1] == 'f':
1102                    print('\n\t',i[-3],' =')
1103                    for m,j in i[:-3]:
1104                        print('\t\t+',m,' * ',j,'   ',repr(j))
1105                else:
1106                    print('  unexpected: ',repr(i))
1107            print("\nG2name ==> ISODISTORT full name",
1108                      ".Phase['ISODISTORT']['IsoModeList']",
1109                      ".Phase['ISODISTORT']['G2ModeList']")
1110            for mode,G2mode in zip(modelist,modeVarList):
1111                print("  ?::"+str(G2mode),' ==>', mode)
1112
1113        #======================================================================
1114        # debug: show occupancy mode var to mode relations
1115        if 'G2OccVarList' in self.Phase['ISODISTORT']:
1116            # coordinate items
1117            #occVarLbl = self.Phase['ISODISTORT']['OccVarList']
1118            G2varObj = self.Phase['ISODISTORT']['G2OccVarList']
1119            #ParentOcc = self.Phase['ISODISTORT']['BaseOcc']
1120            # mode items
1121            modelist = self.Phase['ISODISTORT']['OccModeList']
1122            modeVarList = self.Phase['ISODISTORT']['G2OccModeList']
1123            normlist = self.Phase['ISODISTORT']['OccNormList']
1124            # transform matrices
1125            #occupancymodeInvmatrix = self.Phase['ISODISTORT']['Var2OccMatrix']
1126            #occupancymodematrix = self.Phase['ISODISTORT']['Occ2VarMatrix']
1127           
1128            print( 70*'=')
1129            print('\nVar2OccMatrix' ,'OccVarList' )
1130            for i,row in enumerate(occupancymodeInvmatrix):
1131                l = ''
1132                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
1133                    if k == 0: continue
1134                    if l: l += ' + '
1135                    #l += lbl+' * '+str(k)
1136                    l += str(G2varObj[j])+' * '+str(k)
1137                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
1138
1139            # Get the ISODISTORT offset values
1140            occVarDelta = {}
1141            for lbl,val in zip(
1142                blk.get('_iso_deltaoccupancy_label'),
1143                blk.get('_iso_deltaoccupancy_value'),):
1144                occVarDelta[lbl] = float(val)
1145            modeVarDelta = {}
1146            for lbl,val in zip(
1147                blk.get('_iso_occupancymode_label'),
1148                blk.get('_iso_occupancymode_value'),):
1149                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
1150
1151            print( 70*'=')
1152            print('\nInverse relations using Var2OccModeMatrix, OccNormList, OccVarList')
1153            # compute the mode values from the reported coordinate deltas
1154            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1155                l = ''
1156                for lbl,k in zip(occVarLbl,row):
1157                    if k == 0: continue
1158                    if l: l += ' + '
1159                    l += lbl+' * '+str(k)
1160                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
1161            print('\nCalculation checks\n')
1162            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1163                #l = ''
1164                sl = ''
1165                s = 0.
1166                for lbl,k in zip(occVarLbl,row):
1167                    if k == 0: continue
1168                    #if l: l += ' + '
1169                    #l += lbl+' * '+str(k)
1170                    if sl: sl += ' + '
1171                    sl += str(occVarDelta[lbl])+' * '+str(k)
1172                    s += occVarDelta[lbl] * k
1173                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1174                print(' ?= ',modeVarDelta[modelist[i]])
1175                print()
1176
1177            print( 70*'=')
1178            print('\nDirect relations using Occ2VarMatrix, OccNormList, OccVarList')
1179            # compute the coordinate displacements from the reported mode values
1180            Occ = {}
1181            for i,lbl,row in zip(range(len(occVarLbl)),occVarLbl,occupancymodematrix):
1182                l = ''
1183                s = 0.0
1184                for j,(k,n) in enumerate(zip(row,normlist)):
1185                    if k == 0: continue
1186                    if l: l += ' + '
1187                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1188                    s += n * modeVarDelta[modelist[j]] * k
1189                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1190                j = lbl.split('_')[0] 
1191                Occ[j] = ParentOcc[j]+s
1192               
1193            # determine the coordinate delta values from deviations from the parent structure
1194            print('\nOccupancy from CIF vs computed')
1195            for atmline in self.Phase['Atoms']:
1196                lbl = atmline[0]
1197                if lbl in Occ: print( lbl,atmline[6],Occ[lbl])
1198
1199            print( 70*'=')
1200            print('\nGenerated constraints')
1201            for i in self.Constraints:
1202                if type(i) is dict:
1203                    print('\nconstraint help dict')
1204                    for key in i:
1205                        print('\t',key,':',i[key])
1206                elif i[-1] == 'f':
1207                    print('\n\t',i[-3],' =')
1208                    for m,j in i[:-3]:
1209                        print('\t\t+',m,' * ',j,'   ',repr(j))
1210                else:
1211                    print('  unexpected: ',repr(i))
1212            print("\nG2name ==> ISODISTORT full name",
1213                      ".Phase['ISODISTORT']['OccModeList']",
1214                      ".Phase['ISODISTORT']['G2OccModeList']")
1215            for mode,G2mode in zip(modelist,modeVarList):
1216                print("  ?::"+str(G2mode),' ==>', mode)
1217               
1218def ISODISTORT_shortLbl(lbl,shortmodelist):
1219    '''Shorten model labels and remove special characters
1220    '''
1221    regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)\[.*?\](.*)',lbl)
1222    # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)[...]BBBBB
1223    # where SSSSS is the parent spacegroup, [x,y,z] is a location, etc.
1224    if regexp:
1225        # this extracts the AAAAA and BBBBB parts of the string
1226        lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1227    else:
1228        # try form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
1229        regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
1230        if regexp:
1231            lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1232    lbl = lbl.replace('order','o')
1233    lbl = lbl.replace('+','_')
1234    lbl = lbl.replace('-','_')
1235    G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
1236               
Note: See TracBrowser for help on using the repository browser.