source: trunk/imports/G2phase_CIF.py @ 4659

Last change on this file since 4659 was 4659, checked in by toby, 13 months ago

extensive Origin 1->2 updates; merge SetupGeneral? into a single routine in G2Elem; read sym ops from CIF; fix sites sym/mult after cood xform

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