source: trunk/imports/G2phase_CIF.py @ 4410

Last change on this file since 4410 was 4410, checked in by toby, 2 years ago

do much more cleanup after phase delete; refactor some importer imports to avoid G2IO

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 58.6 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2020-04-24 16:06:35 +0000 (Fri, 24 Apr 2020) $
4# $Author: toby $
5# $Revision: 4410 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 4410 2020-04-24 16:06:35Z 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 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: 4410 $")
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                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
493                atomlist[1] = G2elem.FixValence(atomlist[1])
494                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
495                self.Phase['Atoms'].append(atomlist)
496                ranIdlookup[atomlist[0]] = atomlist[-1]
497                if atomlist[0] in atomlbllist:
498                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
499                else:
500                    atomlbllist.append(atomlist[0])
501
502                if magnetic and atomlist[0] in magatomlabels:
503                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
504                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
505                        mcol = G2MagDict.get(mkey,-1)
506                        if mcol > 0:
507                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
508                    self.MPhase['Atoms'].append(matomlist)
509                if Super:
510                    Sfrac = np.zeros((4,2))
511                    Sadp = np.zeros((4,12))
512                    Spos = np.zeros((4,6))
513                    Smag = np.zeros((4,6))
514                    nim = -1
515                    waveType = 'Fourier'
516                    if  occCdict:
517                        for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']):
518                            if item == atomlist[0]:
519                                waveType = 'Crenel'
520                                val = occCdict['_atom_site_occ_special_func_crenel_c'][i]
521                                Sfrac[0][0] = cif.get_number_with_esd(val)[0]
522                                val = occCdict['_atom_site_occ_special_func_crenel_w'][i]
523                                Sfrac[0][1] = cif.get_number_with_esd(val)[0]
524                                nim = 1
525                   
526                    if nim >= 0:
527                        Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]]
528                    else:
529                        Sfrac = []
530                    nim = -1
531                    if displFdict:
532                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
533                            if item == atomlist[0]:
534                                waveType = 'Fourier'                               
535                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
536                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
537                                if im != nim:
538                                    nim = im
539                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
540                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
541                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
542                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
543                    if nim >= 0:
544                        Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]]
545                    else:
546                        Spos = []
547                    nim = -1
548                    if UijFdict:
549                        nim = -1
550                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
551                            if item == atomlist[0]:
552                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
553                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
554                                if im != nim:
555                                    nim = im
556                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
557                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
558                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
559                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
560                    if nim >= 0:
561                        Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]]
562                    else:
563                        Sadp = []
564                    nim = -1
565                    if MagFdict:
566                        nim = -1
567                        for i,item in enumerate(MagFdict[Mnames[0]]):
568                            if item == atomlist[0]:
569                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
570                                im = int(MagFdict[Mnames[2]][i])
571                                if im != nim:
572                                    nim = im
573                                val = MagFdict[Mnames[3]][i]
574                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
575                                val = MagFdict[Mnames[4]][i]
576                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
577                    if nim >= 0:
578                        Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]]
579                    else:
580                        Smag = []
581                    SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
582                    if magnetic and atomlist[0] in magatomlabels:
583                        matomlist.append(SSdict)
584                    atomlist.append(SSdict)
585            if len(atomlbllist) != len(self.Phase['Atoms']):
586                isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
587            for lbl in phasenamefields: # get a name for the phase
588                name = blk.get(lbl)
589                if name is None:
590                    continue
591                name = name.strip()
592                if name == '?' or name == '.':
593                    continue
594                else:
595                    break
596            else: # no name found, use block name for lack of a better choice
597                name = blknm
598            self.Phase['General']['Name'] = name.strip()
599            self.Phase['General']['Super'] = Super
600            self.Phase = copy.deepcopy(self.Phase)      #clean copy
601            if magnetic:
602                self.MPhase = copy.deepcopy(self.MPhase)    #clean copy
603                self.MPhase['General']['Type'] = 'magnetic'               
604                self.MPhase['General']['Name'] = name.strip()+' mag'
605                self.MPhase['General']['Super'] = Super
606                if Super:
607                    self.MPhase['General']['Modulated'] = True
608                    self.MPhase['General']['SuperVec'] = SuperVec
609                    self.MPhase['General']['SuperSg'] = SuperSg
610                    if self.MPhase['General']['SGData']['SGGray']:
611                        self.MPhase['General']['SuperSg'] += 's'
612                if 'mcif' not in filename:
613                    self.Phase = copy.deepcopy(self.MPhase)
614                    del self.MPhase
615            else:
616                self.MPhase = None
617            if Super:
618                if self.Phase['General']['SGData']['SGGray']:
619                    SGData = self.Phase['General']['SGData']
620                    SGData['SGGray'] = False
621                    ncen = len(SGData['SGCen'])
622                    SGData['SGCen'] = SGData['SGCen'][:ncen//2]
623                    self.Phase['General']['SGData'].update(SGData)
624                self.Phase['General']['Modulated'] = True
625                self.Phase['General']['SuperVec'] = SuperVec
626                self.Phase['General']['SuperSg'] = SuperSg
627                if not self.Phase['General']['SGData']['SGFixed']:
628                    self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
629            if self.ISODISTORT_test(blk):
630                if isodistort_warnings:
631                    self.warnings += isodistort_warnings
632                else:
633                    self.errors = "Error while processing ISODISTORT constraints"
634                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
635                    self.errors = ""
636            returnstat = True
637        return returnstat
638       
639    def ISODISTORT_test(self,blk):
640        '''Test if there is any ISODISTORT information in CIF
641
642        At present only _iso_displacivemode... and _iso_occupancymode... are
643        tested.
644        '''
645        for i in ('_iso_displacivemode_label',
646                  '_iso_occupancymode_label'):
647            if blk.get(i): return True
648        return False
649       
650    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
651        '''Process ISODISTORT items to create constraints etc.
652        Constraints are generated from information extracted from
653        loops beginning with _iso_ and are placed into
654        self.Constraints, which contains a list of
655        :ref:`constraints tree items <Constraint_definitions_table>`
656        and one dict.
657        The dict contains help text for each generated ISODISTORT variable
658
659        At present only _iso_displacivemode... and _iso_occupancymode... are
660        processed. Not yet processed: _iso_magneticmode...,
661        _iso_rotationalmode... & _iso_strainmode...
662        '''
663        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
664        'Maps ISODISTORT parm names to GSAS-II names'
665        # used for all types of modes
666        self.Constraints = []
667        explaination = {}
668        #----------------------------------------------------------------------
669        # read in the ISODISTORT displacement modes
670        #----------------------------------------------------------------------
671        if blk.get('_iso_displacivemode_label'):
672            modelist = []
673            shortmodelist = []
674            idlist = []
675            for id,lbl in zip(
676                blk.get('_iso_displacivemode_ID'),
677                blk.get('_iso_displacivemode_label')):
678                idlist.append(int(id))
679                modelist.append(lbl)
680                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
681            # just in case the items are not ordered increasing by id, sort them here
682            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
683            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
684            # read in the coordinate offset variables names and map them to G2 names/objects
685            coordVarLbl = []
686            G2varObj = []
687            idlist = []
688            error = False
689            for id,lbl in zip(
690                blk.get('_iso_deltacoordinate_ID'),
691                blk.get('_iso_deltacoordinate_label') ):
692                idlist.append(int(id))
693                coordVarLbl.append(lbl)
694                if '_' in lbl:
695                    albl = lbl[:lbl.rfind('_')]
696                    vlbl = lbl[lbl.rfind('_')+1:]
697                else:
698                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
699                    error = True
700                    continue
701                if albl not in atomlbllist:
702                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
703                    error = True
704                    continue
705                else:
706                    anum = atomlbllist.index(albl)
707                var = varLookup.get(vlbl)
708                if not var:
709                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
710                    error = True
711                    continue
712                G2varObj.append(G2obj.G2VarObj(
713                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
714                    ))
715            if error:
716                raise Exception("Error decoding variable labels")
717            # just in case the items are not ordered increasing by id, sort them here
718            coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
719            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
720
721            if len(G2varObj) != len(modelist):
722                print ("non-square input")
723                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
724
725            error = False
726            ParentCoordinates = {}
727            for lbl,exp in zip(
728                blk.get('_iso_coordinate_label'),
729                blk.get('_iso_coordinate_formula') ):
730                if '_' in lbl:
731                    albl = lbl[:lbl.rfind('_')]
732                    vlbl = lbl[lbl.rfind('_')+1:]
733                else:
734                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
735                    error = True
736                    continue
737                if vlbl not in 'xyz' or len(vlbl) != 1:
738                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
739                    error = True
740                    continue
741                i = 'xyz'.index(vlbl)
742                if not ParentCoordinates.get(albl):
743                    ParentCoordinates[albl] = [None,None,None]
744                if '+' in exp:
745                    val = exp.split('+')[0].strip()
746                    val = G2p3.FormulaEval(val)
747                    if val is None:
748                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
749                        error = True
750                        continue
751                    ParentCoordinates[albl][i] = val
752                else:
753                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
754            if error:
755                print (self.warnings)
756                raise Exception("Error decoding variable labels")
757            # get mapping of modes to atomic coordinate displacements
758            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
759            for row,col,val in zip(
760                blk.get('_iso_displacivemodematrix_row'),
761                blk.get('_iso_displacivemodematrix_col'),
762                blk.get('_iso_displacivemodematrix_value'),):
763                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
764            # Invert to get mapping of atom displacements to modes
765            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
766            # create the constraints
767            modeVarList = []
768            for i,row in enumerate(displacivemodeInvmatrix):
769                constraint = []
770                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
771                    if k == 0: continue
772                    constraint.append([k,G2varObj[j]])
773                modeVar = G2obj.G2VarObj(
774                    (self.Phase['ranId'],None,shortmodelist[i],None))
775                modeVarList.append(modeVar)               
776                constraint += [modeVar,False,'f']
777                self.Constraints.append(constraint)
778            # normilization constants
779            normlist = []
780            idlist = []
781            for id,exp in zip(
782                blk.get('_iso_displacivemodenorm_ID'),
783                blk.get('_iso_displacivemodenorm_value'),
784                ):
785                idlist.append(int(id))
786                normlist.append(float(exp))
787            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
788            #----------------------------------------------------------------------
789            # save the ISODISTORT info for "mode analysis"
790            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
791            self.Phase['ISODISTORT'].update({
792                # coordinate items
793                'IsoVarList' : coordVarLbl,
794                'G2VarList' : G2varObj,
795                'ParentStructure' : ParentCoordinates,
796                # mode items
797                'IsoModeList' : modelist,
798                'G2ModeList' : modeVarList,
799                'NormList' : normlist,
800                # transform matrices
801                'Var2ModeMatrix' : displacivemodeInvmatrix,
802                'Mode2VarMatrix' : displacivemodematrix,
803                })
804            # make explaination dictionary
805            for mode,shortmode in zip(modelist,shortmodelist):
806                explaination[str(self.Phase['ranId'])+shortmode] = (
807                    "ISODISTORT full name "+str(mode))
808        #----------------------------------------------------------------------
809        # now read in the ISODISTORT occupancy modes
810        #----------------------------------------------------------------------
811        if blk.get('_iso_occupancymode_label'):
812            modelist = []
813            shortmodelist = []
814            idlist = []
815            for id,lbl in zip(
816                blk.get('_iso_occupancymode_ID'),
817                blk.get('_iso_occupancymode_label')):
818                idlist.append(int(id))
819                modelist.append(lbl)
820                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
821            # just in case the items are not ordered increasing by id, sort them here
822            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
823            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
824            # read in the coordinate offset variables names and map them to G2 names/objects
825            occVarLbl = []
826            G2varObj = []
827            idlist = []
828            error = False
829            for id,lbl in zip(
830                blk.get('_iso_deltaoccupancy_ID'),
831                blk.get('_iso_deltaoccupancy_label') ):
832                idlist.append(int(id))
833                occVarLbl.append(lbl)
834                if '_' in lbl:
835                    albl = lbl[:lbl.rfind('_')]
836                    vlbl = lbl[lbl.rfind('_')+1:]
837                else:
838                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
839                    error = True
840                    continue
841                if albl not in atomlbllist:
842                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
843                    error = True
844                    continue
845                else:
846                    anum = atomlbllist.index(albl)
847                var = varLookup.get(vlbl)
848                if not var:
849                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
850                    error = True
851                    continue
852                G2varObj.append(G2obj.G2VarObj(
853                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
854                    ))
855            if error:
856                raise Exception("Error decoding variable labels")
857            # just in case the items are not ordered increasing by id, sort them here
858            occVarLbl = [i for i,j in sorted(zip(occVarLbl,idlist),key=lambda k:k[1])]
859            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
860
861            if len(G2varObj) != len(modelist):
862                print ("non-square input")
863                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
864
865            error = False
866            ParentOcc = {}
867            for lbl,exp in zip(
868                blk.get('_iso_occupancy_label'),
869                blk.get('_iso_occupancy_formula') ):
870                if '_' in lbl:
871                    albl = lbl[:lbl.rfind('_')]
872                    vlbl = lbl[lbl.rfind('_')+1:]
873                else:
874                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
875                    error = True
876                    continue
877                if vlbl != 'occ':
878                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
879                    error = True
880                    continue
881                if '+' in exp:
882                    val = exp.split('+')[0].strip()
883                    val = G2p3.FormulaEval(val)
884                    if val is None:
885                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
886                        error = True
887                        continue
888                    ParentOcc[albl] = val
889            if error:
890                raise Exception("Error decoding occupancy labels")
891            # get mapping of modes to atomic coordinate displacements
892            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
893            for row,col,val in zip(
894                blk.get('_iso_occupancymodematrix_row'),
895                blk.get('_iso_occupancymodematrix_col'),
896                blk.get('_iso_occupancymodematrix_value'),):
897                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
898            # Invert to get mapping of atom displacements to modes
899            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
900            # create the constraints
901            modeVarList = []
902            for i,row in enumerate(occupancymodeInvmatrix):
903                constraint = []
904                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
905                    if k == 0: continue
906                    constraint.append([k,G2varObj[j]])
907                modeVar = G2obj.G2VarObj(
908                    (self.Phase['ranId'],None,shortmodelist[i],None))
909                modeVarList.append(modeVar)               
910                constraint += [modeVar,False,'f']
911                self.Constraints.append(constraint)
912            # normilization constants
913            normlist = []
914            idlist = []
915            for id,exp in zip(
916                blk.get('_iso_occupancymodenorm_ID'),
917                blk.get('_iso_occupancymodenorm_value'),
918                ):
919                idlist.append(int(id))
920                normlist.append(float(exp))
921            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
922            #----------------------------------------------------------------------
923            # save the ISODISTORT info for "mode analysis"
924            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
925            self.Phase['ISODISTORT'].update({
926                # coordinate items
927                'OccVarList' : occVarLbl,
928                'G2OccVarList' : G2varObj,
929                'BaseOcc' : ParentOcc,
930                # mode items
931                'OccModeList' : modelist,
932                'G2OccModeList' : modeVarList,
933                'OccNormList' : normlist,
934                # transform matrices
935                'Var2OccMatrix' : occupancymodeInvmatrix,
936                'Occ2VarMatrix' : occupancymodematrix,
937                })
938            # make explaination dictionary
939            for mode,shortmode in zip(modelist,shortmodelist):
940                explaination[str(self.Phase['ranId'])+shortmode] = (
941                    "ISODISTORT full name "+str(mode))
942        #----------------------------------------------------------------------
943        # done with read
944        #----------------------------------------------------------------------
945        if explaination: self.Constraints.append(explaination)
946
947        if not debug: return
948
949        # debug: show displacive mode var to mode relations
950        if 'IsoVarList' in self.Phase['ISODISTORT']:
951            # coordinate items
952            #coordVarLbl = self.Phase['ISODISTORT']['IsoVarList']
953            G2varObj = self.Phase['ISODISTORT']['G2VarList']
954            #ParentCoordinates = self.Phase['ISODISTORT']['ParentStructure']
955            # mode items
956            modelist = self.Phase['ISODISTORT']['IsoModeList']
957            modeVarList = self.Phase['ISODISTORT']['G2ModeList']
958            normlist = self.Phase['ISODISTORT']['NormList']
959            # transform matrices
960            #displacivemodeInvmatrix = self.Phase['ISODISTORT']['Var2ModeMatrix']
961            #displacivemodematrix = self.Phase['ISODISTORT']['Mode2VarMatrix']
962       
963            print( 70*'=')
964            print('\nVar2ModeMatrix' ,'IsoVarList' )
965            for i,row in enumerate(displacivemodeInvmatrix):
966                l = ''
967                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
968                    if k == 0: continue
969                    if l: l += ' + '
970                    #l += lbl+' * '+str(k)
971                    l += str(G2varObj[j])+' * '+str(k)
972                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
973
974            # Get the ISODISTORT offset values
975            coordVarDelta = {}
976            for lbl,val in zip(
977                blk.get('_iso_deltacoordinate_label'),
978                blk.get('_iso_deltacoordinate_value'),):
979                coordVarDelta[lbl] = float(val)
980            modeVarDelta = {}
981            for lbl,val in zip(
982                blk.get('_iso_displacivemode_label'),
983                blk.get('_iso_displacivemode_value'),):
984                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
985
986            print( 70*'=')
987            print('\nInverse relations using Var2ModeMatrix, NormList, IsoVarList')
988            # compute the mode values from the reported coordinate deltas
989            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
990                l = ''
991                for lbl,k in zip(coordVarLbl,row):
992                    if k == 0: continue
993                    if l: l += ' + '
994                    l += lbl+' * '+str(k)
995                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
996            print('\nCalculation checks\n')
997            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
998                #l = ''
999                sl = ''
1000                s = 0.
1001                for lbl,k in zip(coordVarLbl,row):
1002                    if k == 0: continue
1003                    #if l: l += ' + '
1004                    #l += lbl+' * '+str(k)
1005                    if sl: sl += ' + '
1006                    sl += str(coordVarDelta[lbl])+' * '+str(k)
1007                    s += coordVarDelta[lbl] * k
1008                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1009                print(' ?= ',modeVarDelta[modelist[i]])
1010                print()
1011
1012            print( 70*'=')
1013            print('\nDirect relations using Mode2VarMatrix, NormList, IsoVarList')
1014            # compute the coordinate displacements from the reported mode values
1015            DeltaCoords = {}
1016            for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
1017                l = ''
1018                s = 0.0
1019                at,d = lbl.split('_')
1020                if at not in DeltaCoords:
1021                    DeltaCoords[at] = [0,0,0]
1022                for j,(k,n) in enumerate(zip(row,normlist)):
1023                    if k == 0: continue
1024                    if l: l += ' + '
1025                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1026                    s += n * modeVarDelta[modelist[j]] * k
1027                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1028                if d == 'dx':
1029                    DeltaCoords[at][0] = s
1030                elif d == 'dy':
1031                    DeltaCoords[at][1] = s
1032                elif d == 'dz':
1033                    DeltaCoords[at][2] = s
1034                else:
1035                    print('unexpected',d)
1036
1037            print( 70*'=')
1038            print('\nCoordinate checks')
1039            print('\nComputed')
1040            for at in sorted(DeltaCoords):
1041                s = at
1042                for i in range(3):
1043                    s += '  '
1044                    s += str(ParentCoordinates[at][i]+DeltaCoords[at][i])
1045                print(s)
1046
1047            # determine the coordinate delta values from deviations from the parent structure
1048            print('\nFrom CIF')
1049            for atmline in self.Phase['Atoms']:
1050                lbl = atmline[0]
1051                x,y,z = atmline[3:6]
1052                print( lbl,x,y,z)
1053
1054            print( 70*'=')
1055            print('\nGenerated constraints')
1056            for i in self.Constraints:
1057                if type(i) is dict:
1058                    print('\nconstraint help dict')
1059                    for key in i:
1060                        print('\t',key,':',i[key])
1061                elif i[-1] == 'f':
1062                    print('\n\t',i[-3],' =')
1063                    for m,j in i[:-3]:
1064                        print('\t\t+',m,' * ',j,'   ',repr(j))
1065                else:
1066                    print('  unexpected: ',repr(i))
1067            print("\nG2name ==> ISODISTORT full name",
1068                      ".Phase['ISODISTORT']['IsoModeList']",
1069                      ".Phase['ISODISTORT']['G2ModeList']")
1070            for mode,G2mode in zip(modelist,modeVarList):
1071                print("  ?::"+str(G2mode),' ==>', mode)
1072
1073        #======================================================================
1074        # debug: show occupancy mode var to mode relations
1075        if 'G2OccVarList' in self.Phase['ISODISTORT']:
1076            # coordinate items
1077            #occVarLbl = self.Phase['ISODISTORT']['OccVarList']
1078            G2varObj = self.Phase['ISODISTORT']['G2OccVarList']
1079            #ParentOcc = self.Phase['ISODISTORT']['BaseOcc']
1080            # mode items
1081            modelist = self.Phase['ISODISTORT']['OccModeList']
1082            modeVarList = self.Phase['ISODISTORT']['G2OccModeList']
1083            normlist = self.Phase['ISODISTORT']['OccNormList']
1084            # transform matrices
1085            #occupancymodeInvmatrix = self.Phase['ISODISTORT']['Var2OccMatrix']
1086            #occupancymodematrix = self.Phase['ISODISTORT']['Occ2VarMatrix']
1087           
1088            print( 70*'=')
1089            print('\nVar2OccMatrix' ,'OccVarList' )
1090            for i,row in enumerate(occupancymodeInvmatrix):
1091                l = ''
1092                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
1093                    if k == 0: continue
1094                    if l: l += ' + '
1095                    #l += lbl+' * '+str(k)
1096                    l += str(G2varObj[j])+' * '+str(k)
1097                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
1098
1099            # Get the ISODISTORT offset values
1100            occVarDelta = {}
1101            for lbl,val in zip(
1102                blk.get('_iso_deltaoccupancy_label'),
1103                blk.get('_iso_deltaoccupancy_value'),):
1104                occVarDelta[lbl] = float(val)
1105            modeVarDelta = {}
1106            for lbl,val in zip(
1107                blk.get('_iso_occupancymode_label'),
1108                blk.get('_iso_occupancymode_value'),):
1109                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
1110
1111            print( 70*'=')
1112            print('\nInverse relations using Var2OccModeMatrix, OccNormList, OccVarList')
1113            # compute the mode values from the reported coordinate deltas
1114            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1115                l = ''
1116                for lbl,k in zip(occVarLbl,row):
1117                    if k == 0: continue
1118                    if l: l += ' + '
1119                    l += lbl+' * '+str(k)
1120                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
1121            print('\nCalculation checks\n')
1122            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1123                #l = ''
1124                sl = ''
1125                s = 0.
1126                for lbl,k in zip(occVarLbl,row):
1127                    if k == 0: continue
1128                    #if l: l += ' + '
1129                    #l += lbl+' * '+str(k)
1130                    if sl: sl += ' + '
1131                    sl += str(occVarDelta[lbl])+' * '+str(k)
1132                    s += occVarDelta[lbl] * k
1133                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1134                print(' ?= ',modeVarDelta[modelist[i]])
1135                print()
1136
1137            print( 70*'=')
1138            print('\nDirect relations using Occ2VarMatrix, OccNormList, OccVarList')
1139            # compute the coordinate displacements from the reported mode values
1140            Occ = {}
1141            for i,lbl,row in zip(range(len(occVarLbl)),occVarLbl,occupancymodematrix):
1142                l = ''
1143                s = 0.0
1144                for j,(k,n) in enumerate(zip(row,normlist)):
1145                    if k == 0: continue
1146                    if l: l += ' + '
1147                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1148                    s += n * modeVarDelta[modelist[j]] * k
1149                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1150                j = lbl.split('_')[0] 
1151                Occ[j] = ParentOcc[j]+s
1152               
1153            # determine the coordinate delta values from deviations from the parent structure
1154            print('\nOccupancy from CIF vs computed')
1155            for atmline in self.Phase['Atoms']:
1156                lbl = atmline[0]
1157                if lbl in Occ: print( lbl,atmline[6],Occ[lbl])
1158
1159            print( 70*'=')
1160            print('\nGenerated constraints')
1161            for i in self.Constraints:
1162                if type(i) is dict:
1163                    print('\nconstraint help dict')
1164                    for key in i:
1165                        print('\t',key,':',i[key])
1166                elif i[-1] == 'f':
1167                    print('\n\t',i[-3],' =')
1168                    for m,j in i[:-3]:
1169                        print('\t\t+',m,' * ',j,'   ',repr(j))
1170                else:
1171                    print('  unexpected: ',repr(i))
1172            print("\nG2name ==> ISODISTORT full name",
1173                      ".Phase['ISODISTORT']['OccModeList']",
1174                      ".Phase['ISODISTORT']['G2OccModeList']")
1175            for mode,G2mode in zip(modelist,modeVarList):
1176                print("  ?::"+str(G2mode),' ==>', mode)
1177               
1178def ISODISTORT_shortLbl(lbl,shortmodelist):
1179    '''Shorten model labels and remove special characters
1180    '''
1181    regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)\[.*?\](.*)',lbl)
1182    # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)[...]BBBBB
1183    # where SSSSS is the parent spacegroup, [x,y,z] is a location, etc.
1184    if regexp:
1185        # this extracts the AAAAA and BBBBB parts of the string
1186        lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1187    else:
1188        # try form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
1189        regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
1190        if regexp:
1191            lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1192    lbl = lbl.replace('order','o')
1193    lbl = lbl.replace('+','_')
1194    lbl = lbl.replace('-','_')
1195    G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
1196               
Note: See TracBrowser for help on using the repository browser.