source: trunk/imports/G2phase_CIF.py @ 4476

Last change on this file since 4476 was 4476, checked in by vondreele, 2 years ago

make sure if Uij = None that it is zero, type is 'I' and that after a reimport atoms that Draw Atoms is updated

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