source: trunk/imports/G2phase_CIF.py @ 4832

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

fix gaim map calc to include cylindrical sample absorption
fix phase cif importer

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