source: trunk/imports/G2phase_CIF.py @ 5062

Last change on this file since 5062 was 5062, checked in by toby, 8 months ago

Start on isodisplace import cleanup & formatting; plotting changes for conv plots

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