source: trunk/imports/G2phase_CIF.py @ 3338

Last change on this file since 3338 was 3338, checked in by vondreele, 4 years ago

fix SpnFlp? issues for mcif imports - had messed up mag structure factors
make correct Mag moments for crystal axes aligned components - both structure factors & derivatives
checked with very oblique sarcopside.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 45.2 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-04-06 20:35:44 +0000 (Fri, 06 Apr 2018) $
4# $Author: vondreele $
5# $Revision: 3338 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3338 2018-04-06 20:35:44Z 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 GSASIIIO as G2IO
27import GSASIIobj as G2obj
28import GSASIIspc as G2spc
29import GSASIIElem as G2elem
30import GSASIIlattice as G2lat
31import GSASIIpy3 as G2p3
32import GSASIIpath
33GSASIIpath.SetVersionNumber("$Revision: 3338 $")
34import CifFile as cif # PyCifRW from James Hester
35
36class CIFPhaseReader(G2obj.ImportPhase):
37    'Implements a phase importer from a possibly multi-block CIF file'
38    def __init__(self):
39        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
40            extensionlist=('.CIF','.cif','.mcif'),
41            strictExtension=False,
42            formatName = 'CIF',
43            longFormatName = 'Crystallographic Information File import'
44            )
45       
46    def ContentsValidator(self, filename):
47        fp = open(filename,'r')
48        return self.CIFValidator(fp)
49        fp.close()
50
51    def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused):
52        self.isodistort_warnings = ''
53        self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
54        # make sure the ranId is really unique!
55        while self.Phase['ranId'] in usedRanIdList:
56            self.Phase['ranId'] = ran.randint(0,sys.maxsize)
57        self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
58        # make sure the ranId is really unique!
59        while self.MPhase['ranId'] in usedRanIdList:
60            self.MPhase['ranId'] = ran.randint(0,sys.maxsize)
61        returnstat = False
62        cellitems = (
63            '_cell_length_a','_cell_length_b','_cell_length_c',
64            '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
65#        cellwaveitems = (
66#            '_cell_wave_vector_seq_id',
67#            '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z')
68        reqitems = (
69             '_atom_site_fract_x',
70             '_atom_site_fract_y',
71             '_atom_site_fract_z',
72            )
73        phasenamefields = (
74            '_chemical_name_common',
75            '_pd_phase_name',
76            '_chemical_formula_sum'
77            )
78        try:
79            cf = G2obj.ReadCIF(filename)
80        except cif.StarError as msg:
81            msg  = 'Unreadable cif file\n'+str(msg)
82            self.errors = msg
83            self.warnings += msg
84            return False
85        # scan blocks for structural info
86        self.errors = 'Error during scan of blocks for datasets'
87        str_blklist = []
88        for blk in cf.keys():
89            for r in reqitems+cellitems:
90                if r not in cf[blk].keys():
91                    break
92            else:
93                str_blklist.append(blk)
94        if not str_blklist:
95            selblk = None # no block to choose
96        elif len(str_blklist) == 1: # only one choice
97            selblk = 0
98        else:                       # choose from options
99            choice = []
100            for blknm in str_blklist:
101                choice.append('')
102                # accumumlate some info about this phase
103                choice[-1] += blknm + ': '
104                for i in phasenamefields: # get a name for the phase
105                    name = cf[blknm].get(i,'phase name').strip()
106                    if name is None or name == '?' or name == '.':
107                        continue
108                    else:
109                        choice[-1] += name.strip()[:20] + ', '
110                        break
111                na = len(cf[blknm].get("_atom_site_fract_x"))
112                if na == 1:
113                    choice[-1] += '1 atom'
114                else:
115                    choice[-1] += ('%d' % na) + ' atoms'
116                choice[-1] += ', cell: '
117                fmt = "%.2f,"
118                for i,key in enumerate(cellitems):
119                    if i == 3: fmt = "%.f,"
120                    if i == 5: fmt = "%.f"
121                    choice[-1] += fmt % cif.get_number_with_esd(
122                        cf[blknm].get(key))[0]
123                sg = cf[blknm].get("_symmetry_space_group_name_H-M",'')
124                if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'')
125                if not sg: sg = cf[blknm].get("_space_group_ssg_name",'')
126                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'')
127                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'')
128                #how about checking for super/magnetic ones as well? - reject 'X'?
129                sg = sg.replace('_','')
130                if sg: choice[-1] += ', (' + sg.strip() + ')'
131            selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame,
132                title= 'Select a phase from one the CIF data_ blocks below',size=(600,100))
133        self.errors = 'Error during reading of selected block'
134#process selected phase
135        if selblk is None:
136            returnstat = False # no block selected or available
137        else:   #do space group symbol & phase type first
138            blknm = str_blklist[selblk]
139            blk = cf[str_blklist[selblk]]
140            E = True
141            Super = False
142            magnetic = False
143            moddim = int(blk.get("_cell_modulation_dimension",'0'))
144            if moddim:      #incommensurate
145                if moddim > 1:
146                    msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
147                    self.errors = msg
148                    self.warnings += '\n'+msg
149                    return False
150                if blk.get('_cell_subsystems_number'):
151                    msg = 'Composite super structures not allowed in GSAS-II'
152                    self.errors = msg
153                    self.warnings += '\n'+msg
154                    return False
155                sspgrp = blk.get("_space_group_ssg_name",'')
156                if not sspgrp:          #might be incommensurate magnetic
157                    MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'')
158                    if not MSSpGrp:
159                        MSSpGrp = blk.get("_space_group.magn_ssg_name",'') 
160                    if not MSSpGrp:
161                        msg = 'No incommensurate space group name was found in the CIF.'
162                        self.errors = msg
163                        self.warnings += '\n'+msg
164                        return False                                                           
165                    if 'X' in MSSpGrp:
166                        msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II'
167                        self.warnings += '\n'+msg
168                        self.errors = msg
169                        return False
170                    magnetic = True
171                if 'X' in sspgrp:
172                    msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II'
173                    self.warnings += '\n'+msg
174                    self.errors = msg
175                    return False
176                Super = True
177                if magnetic:
178                    sspgrp = MSSpGrp.split('(')
179                    sspgrp[1] = "("+sspgrp[1]
180                    SpGrp = G2spc.StandardizeSpcName(sspgrp[0])
181                    if "1'" in SpGrp: sspgrp[1] = sspgrp[1][:-1]  #take off extra 's'; gets put back later
182                    MSpGrp = sspgrp[0]
183                    self.MPhase['General']['Type'] = 'magnetic'
184                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
185                else:
186                    sspgrp = sspgrp.split('(')
187                    sspgrp[1] = "("+sspgrp[1]
188                    SpGrp = sspgrp[0]
189                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
190                    self.Phase['General']['Type'] = 'nuclear'
191                if not SpGrp:
192                    print (sspgrp)
193                    self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II'
194                    return False
195                SuperSg = sspgrp[1].replace('\\','')
196                SuperVec = [[0,0,.1],False,4]
197            else:   #not incommensurate
198                SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
199                if not SpGrp:
200                    SpGrp = blk.get("_space_group_name_H-M_alt",'')
201                if not SpGrp:   #try magnetic           
202                    MSpGrp = blk.get("_space_group.magn_name_BNS",'')
203                    if not MSpGrp:
204                        MSpGrp = blk.get("_space_group_magn.name_BNS",'')
205                        if not MSpGrp:
206                            msg = 'No recognizable space group name was found in the CIF.'
207                            self.errors = msg
208                            self.warnings += '\n'+msg
209                            return False
210                    SpGrp = blk.get('_parent_space_group.name_H-M_alt')
211                    if not SpGrp:
212                        SpGrp = blk.get('_parent_space_group.name_H-M')
213#                    SpGrp = MSpGrp.replace("'",'')
214                    SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','')   #get rid of screw '_'
215                    if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:]
216                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
217                    magnetic = True
218                    self.MPhase['General']['Type'] = 'magnetic'
219                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
220                    if not SpGrp:
221                        print (MSpGrp)
222                        self.warnings += 'No space group name was found in the CIF.'
223                        return False
224                else:
225                    SpGrp = SpGrp.replace('_','')
226                    self.Phase['General']['Type'] = 'nuclear'
227#process space group symbol
228            E,SGData = G2spc.SpcGroup(SpGrp)
229            if E and SpGrp:
230                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
231                if SpGrpNorm:
232                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
233            # nope, try the space group "out of the Box"
234            if E:
235                self.warnings += 'ERROR in space group symbol '+SpGrp
236                self.warnings += '\nThe space group has been set to "P 1". '
237                self.warnings += "Change this in phase's General tab."
238                self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
239                self.warnings += G2spc.SGErrors(E)
240                SGData = G2obj.P1SGData # P 1
241            self.Phase['General']['SGData'] = SGData
242            if Super:
243                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
244                if E:
245                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
246                    self.warnings += '\n'+E
247                    SuperSg = SuperSg[:SuperSg.index(')')+1]
248                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
249                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
250                self.Phase['General']['SSGData'] = SSGData
251                if magnetic:
252                    self.MPhase['General']['SGData'] = SGData
253                    self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
254                    self.MPhase['General']['SSGData'] = SSGData
255
256            if magnetic:    #replace std operaors with those from cif file - probably not the same!
257                SGData['SGFixed'] = True
258                SGData['SGOps'] = []
259                SGData['SGCen'] = []
260                if Super:
261                    SSGData['SSGOps'] = []
262                    SSGData['SSGCen'] = []
263                    try:
264                        sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id')
265                        sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id')
266                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1]
267                        centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1]                   
268                    except KeyError:        #old mag cif names
269                        sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id')
270                        sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id')
271                        opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1]
272                        centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1]
273                    spnflp = []
274                    for op in sgoploop:
275                        M,T,S = G2spc.MagSSText2MTS(op[opid])
276                        SSGData['SSGOps'].append([np.array(M,dtype=float),T])
277                        SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]])
278                        spnflp.append(S)
279                    censpn = []
280                    for cent in sgcenloop:
281                        M,C,S = G2spc.MagSSText2MTS(cent[centid])
282                        SSGData['SSGCen'].append(C)
283                        SGData['SGCen'].append(C[:3])
284                        censpn += list(np.array(spnflp)*S)
285                    self.MPhase['General']['SSGData'] = SSGData
286                else:   
287                    try:
288                        sgoploop = blk.GetLoop('_space_group_symop_magn.id')
289                        sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
290                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
291                        centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]                   
292                    except KeyError:        #old mag cif names
293                        sgoploop = blk.GetLoop('_space_group_symop.magn_id')
294                        sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
295                        opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
296                        centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
297                    spnflp = []
298                    for op in sgoploop:
299                        try:
300                            M,T,S = G2spc.MagText2MTS(op[opid])
301                            SGData['SGOps'].append([np.array(M,dtype=float),T])
302                            spnflp.append(S)
303                        except KeyError:
304                            self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II'
305                            return False
306                    censpn = []
307                    for cent in sgcenloop:
308                        M,C,S = G2spc.MagText2MTS(cent[centid])
309                        SGData['SGCen'].append(C)
310                        censpn += list(np.array(spnflp)*S)
311                self.MPhase['General']['SGData'] = SGData
312                self.MPhase['General']['SGData']['SpnFlp'] = censpn
313                G2spc.GenMagOps(SGData)         #set magMom
314                self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
315                MagPtGp = blk.get('_space_group.magn_point_group')
316                if not MagPtGp:
317                    MagPtGp = blk.get('_space_group.magn_point_group_name')
318                if not MagPtGp:
319                    MagPtGp = blk.get('_space_group_magn.point_group_name')
320                self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
321
322            # cell parameters
323            cell = []
324            for lbl in cellitems:
325                cell.append(cif.get_number_with_esd(blk[lbl])[0])
326            Volume = G2lat.calc_V(G2lat.cell2A(cell))
327            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
328            if magnetic:
329                self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]               
330            if Super:
331                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
332                waveDict = dict(waveloop.items())
333                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
334                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
335                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
336
337            # read in atoms
338            self.errors = 'Error during reading of atoms'
339            atomlbllist = [] # table to look up atom IDs
340            atomloop = blk.GetLoop('_atom_site_label')
341            atomkeys = [i.lower() for i in atomloop.keys()]
342            if not blk.get('_atom_site_type_symbol'):
343                self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
344            if magnetic:
345                try:
346                    magmoment = '_atom_site_moment.label'
347                    magatomloop = blk.GetLoop(magmoment)
348                    magatomkeys = [i.lower() for i in magatomloop.keys()]
349                    magatomlabels = blk.get(magmoment)
350                    G2MagDict = {'_atom_site_moment.label': 0,
351                                 '_atom_site_moment.crystalaxis_x':7,
352                                 '_atom_site_moment.crystalaxis_y':8,
353                                 '_atom_site_moment.crystalaxis_z':9}
354                except KeyError:
355                    magmoment = '_atom_site_moment_label'
356                    magatomloop = blk.GetLoop(magmoment)
357                    magatomkeys = [i.lower() for i in magatomloop.keys()]
358                    magatomlabels = blk.get(magmoment)
359                    G2MagDict = {'_atom_site_moment_label': 0,
360                                 '_atom_site_moment_crystalaxis_x':7,
361                                 '_atom_site_moment_crystalaxis_y':8,
362                                 '_atom_site_moment_crystalaxis_z':9}
363                   
364            if blk.get('_atom_site_aniso_label'):
365                anisoloop = blk.GetLoop('_atom_site_aniso_label')
366                anisokeys = [i.lower() for i in anisoloop.keys()]
367                anisolabels = blk.get('_atom_site_aniso_label')
368            else:
369                anisoloop = None
370                anisokeys = []
371                anisolabels = []
372            if Super:
373                occFloop = None
374                occCloop = None
375                occFdict = {}
376                occCdict = {}
377                displSloop = None
378                displFloop = None
379                MagFloop = None
380                displSdict = {}
381                displFdict = {}
382                MagFdict = {}
383                UijFloop = None
384                UijFdict = {}
385                #occupancy modulation
386                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
387                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
388                    occFdict = dict(occFloop.items())
389                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
390                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
391                    occCdict = dict(occCloop.items())
392                #position modulation
393                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
394                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
395                    displFdict = dict(displFloop.items())                           
396                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
397                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
398                    displSdict = dict(displSloop.items())
399                #U modulation
400                if blk.get('_atom_site_U_Fourier_atom_site_label'):
401                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
402                    UijFdict = dict(UijFloop.items())
403                #Mag moment modulation
404                if blk.get('_atom_site_moment_Fourier_atom_site_label'):
405                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
406                    MagFdict = dict(MagFloop.items())
407                    Mnames =  ['_atom_site_moment_fourier_atom_site_label',
408                               '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
409                               '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
410                elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
411                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
412                    MagFdict = dict(MagFloop.items())                           
413                    Mnames =  ['_atom_site_moment_fourier.atom_site_label',
414                               '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
415                               '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
416            self.Phase['Atoms'] = []
417            if magnetic:
418                self.MPhase['Atoms'] = []
419            G2AtomDict = {  '_atom_site_type_symbol' : 1,
420                            '_atom_site_label' : 0,
421                            '_atom_site_fract_x' : 3,
422                            '_atom_site_fract_y' : 4,
423                            '_atom_site_fract_z' : 5,
424                            '_atom_site_occupancy' : 6,
425                            '_atom_site_aniso_u_11' : 11,
426                            '_atom_site_aniso_u_22' : 12,
427                            '_atom_site_aniso_u_33' : 13,
428                            '_atom_site_aniso_u_12' : 14,
429                            '_atom_site_aniso_u_13' : 15,
430                            '_atom_site_aniso_u_23' : 16, }
431
432            ranIdlookup = {}
433            for aitem in atomloop:
434                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
435                for val,key in zip(aitem,atomkeys):
436                    col = G2AtomDict.get(key,-1)
437                    if col >= 3:
438                        atomlist[col] = cif.get_number_with_esd(val)[0]
439                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
440                    elif col is not None and col != -1:
441                        atomlist[col] = val
442                    elif key in ('_atom_site_thermal_displace_type',
443                               '_atom_site_adp_type'):   #Iso or Aniso?
444                        if val.lower() == 'uani':
445                            atomlist[9] = 'A'
446                    elif key == '_atom_site_u_iso_or_equiv':
447                        uisoval = cif.get_number_with_esd(val)[0]
448                        if uisoval is not None: 
449                            atomlist[10] = uisoval
450                if not atomlist[1] and atomlist[0]:
451                    typ = atomlist[0].rstrip('0123456789-+')
452                    if G2elem.CheckElement(typ):
453                        atomlist[1] = typ
454                    if not atomlist[1]:
455                        atomlist[1] = 'Xe'
456                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
457                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
458                    atomlist[9] = 'A'
459                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
460                        col = G2AtomDict.get(key)
461                        if col:
462                            atomlist[col] = cif.get_number_with_esd(val)[0]
463                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
464                atomlist[1] = G2elem.FixValence(atomlist[1])
465                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
466                self.Phase['Atoms'].append(atomlist)
467                ranIdlookup[atomlist[0]] = atomlist[-1]
468                if atomlist[0] in atomlbllist:
469                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
470                else:
471                    atomlbllist.append(atomlist[0])
472
473                if magnetic and atomlist[0] in magatomlabels:
474                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
475                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
476                        mcol = G2MagDict.get(mkey,-1)
477                        if mcol:
478                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
479                    self.MPhase['Atoms'].append(matomlist)
480                if Super:
481                    Sfrac = np.zeros((4,2))
482                    Sadp = np.zeros((4,12))
483                    Spos = np.zeros((4,6))
484                    Smag = np.zeros((4,6))
485                    nim = -1
486                    waveType = 'Fourier'
487                    if  occCdict:
488                        for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']):
489                            if item == atomlist[0]:
490                                waveType = 'Crenel'
491                                val = occCdict['_atom_site_occ_special_func_crenel_c'][i]
492                                Sfrac[0][0] = cif.get_number_with_esd(val)[0]
493                                val = occCdict['_atom_site_occ_special_func_crenel_w'][i]
494                                Sfrac[0][1] = cif.get_number_with_esd(val)[0]
495                                nim = 1
496                   
497                    if nim >= 0:
498                        Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]]
499                    else:
500                        Sfrac = []
501                    nim = -1
502                    if displFdict:
503                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
504                            if item == atomlist[0]:
505                                waveType = 'Fourier'                               
506                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
507                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
508                                if im != nim:
509                                    nim = im
510                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
511                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
512                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
513                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
514                    if nim >= 0:
515                        Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]]
516                    else:
517                        Spos = []
518                    nim = -1
519                    if UijFdict:
520                        nim = -1
521                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
522                            if item == atomlist[0]:
523                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
524                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
525                                if im != nim:
526                                    nim = im
527                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
528                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
529                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
530                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
531                    if nim >= 0:
532                        Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]]
533                    else:
534                        Sadp = []
535                    nim = -1
536                    if MagFdict:
537                        nim = -1
538                        for i,item in enumerate(MagFdict[Mnames[0]]):
539                            if item == atomlist[0]:
540                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
541                                im = int(MagFdict[Mnames[2]][i])
542                                if im != nim:
543                                    nim = im
544                                val = MagFdict[Mnames[3]][i]
545                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
546                                val = MagFdict[Mnames[4]][i]
547                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
548                    if nim >= 0:
549                        Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]]
550                    else:
551                        Smag = []
552                    SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
553                    if magnetic and atomlist[0] in magatomlabels:
554                        matomlist.append(SSdict)
555                    atomlist.append(SSdict)
556            if len(atomlbllist) != len(self.Phase['Atoms']):
557                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
558            for lbl in phasenamefields: # get a name for the phase
559                name = blk.get(lbl)
560                if name is None:
561                    continue
562                name = name.strip()
563                if name == '?' or name == '.':
564                    continue
565                else:
566                    break
567            else: # no name found, use block name for lack of a better choice
568                name = blknm
569            self.Phase['General']['Name'] = name.strip()[:20]
570            self.Phase['General']['Super'] = Super
571            if magnetic:
572                self.MPhase['General']['Type'] = 'magnetic'               
573                self.MPhase['General']['Name'] = name.strip()[:20]+' mag'
574                self.MPhase['General']['Super'] = Super
575                if Super:
576                    if self.MPhase['General']['SGData']['SGGray']:
577                        SuperSg += 's'
578                    self.MPhase['General']['Modulated'] = True
579                    self.MPhase['General']['SuperVec'] = SuperVec
580                    self.MPhase['General']['SuperSg'] = SuperSg
581            else:
582                self.MPhase = None
583            if Super:
584                self.Phase['General']['Modulated'] = True
585                self.Phase['General']['SuperVec'] = SuperVec
586                self.Phase['General']['SuperSg'] = SuperSg
587                if not self.Phase['General']['SGData']['SGFixed']:
588                    self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
589            if not self.isodistort_warnings:
590                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
591                    self.errors = "Error while processing ISODISTORT constraints"
592                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
593            else:
594                self.warnings += self.isodistort_warnings
595            returnstat = True
596        return returnstat
597
598    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
599        'Process ISODISTORT items to create constraints etc.'
600        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
601        'Maps ISODISTORT parm names to GSAS-II names'
602        #----------------------------------------------------------------------
603        # read in the ISODISTORT displacement modes
604        #----------------------------------------------------------------------
605        self.Constraints = []
606        explaination = {}
607        if blk.get('_iso_displacivemode_label'):
608            modelist = []
609            shortmodelist = []
610            for lbl in blk.get('_iso_displacivemode_label'):
611                modelist.append(lbl)
612                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
613                # where SSSSS is the parent spacegroup, [x,y,z] is a location
614                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
615                # this extracts the AAAAA and BBBBB parts of the string
616                if regexp:
617                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
618                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
619            # read in the coordinate offset variables names and map them to G2 names/objects
620            coordVarLbl = []
621            G2varLbl = []
622            G2varObj = []
623            error = False
624            for lbl in blk.get('_iso_deltacoordinate_label'):
625                coordVarLbl.append(lbl)
626                if '_' in lbl:
627                    albl = lbl[:lbl.rfind('_')]
628                    vlbl = lbl[lbl.rfind('_')+1:]
629                else:
630                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
631                    error = True
632                    continue
633                if albl not in atomlbllist:
634                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
635                    error = True
636                    continue
637                else:
638                    anum = atomlbllist.index(albl)
639                var = varLookup.get(vlbl)
640                if not var:
641                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
642                    error = True
643                    continue
644                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
645                G2varObj.append(G2obj.G2VarObj(
646                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
647                    ))
648            if error:
649                raise Exception("Error decoding variable labels")
650
651            if len(G2varObj) != len(modelist):
652                print ("non-square input")
653                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
654
655            error = False
656            ParentCoordinates = {}
657            for lbl,exp in zip(
658                blk.get('_iso_coordinate_label'),
659                blk.get('_iso_coordinate_formula'),
660                ):
661                if '_' in lbl:
662                    albl = lbl[:lbl.rfind('_')]
663                    vlbl = lbl[lbl.rfind('_')+1:]
664                else:
665                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
666                    error = True
667                    continue
668                if vlbl not in 'xyz' or len(vlbl) != 1:
669                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
670                    error = True
671                    continue
672                i = 'xyz'.index(vlbl)
673                if not ParentCoordinates.get(albl):
674                    ParentCoordinates[albl] = [None,None,None]
675                if '+' in exp:
676                    val = exp.split('+')[0].strip()
677                    val = G2p3.FormulaEval(val)
678                    if val is None:
679                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
680                        error = True
681                        continue
682                    ParentCoordinates[albl][i] = val
683                else:
684                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
685            if error:
686                print (self.warnings)
687                raise Exception("Error decoding variable labels")
688            # get mapping of modes to atomic coordinate displacements
689            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
690            for row,col,val in zip(
691                blk.get('_iso_displacivemodematrix_row'),
692                blk.get('_iso_displacivemodematrix_col'),
693                blk.get('_iso_displacivemodematrix_value'),):
694                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
695            # Invert to get mapping of atom displacements to modes
696            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
697            # create the constraints
698            for i,row in enumerate(displacivemodeInvmatrix):
699                constraint = []
700                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
701                    if k == 0: continue
702                    constraint.append([k,G2varObj[j]])
703                constraint += [shortmodelist[i],False,'f']
704                self.Constraints.append(constraint)
705            #----------------------------------------------------------------------
706            # save the ISODISTORT info for "mode analysis"
707            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
708            self.Phase['ISODISTORT'].update({
709                'IsoModeList' : modelist,
710                'G2ModeList' : shortmodelist,
711                'IsoVarList' : coordVarLbl,
712                'G2VarList' : G2varObj,
713                'ParentStructure' : ParentCoordinates,
714                'Var2ModeMatrix' : displacivemodeInvmatrix,
715                'Mode2VarMatrix' : displacivemodematrix,
716                })
717            # make explaination dictionary
718            for mode,shortmode in zip(modelist,shortmodelist):
719                explaination[shortmode] = "ISODISTORT full name "+str(mode)
720        #----------------------------------------------------------------------
721        # now read in the ISODISTORT occupancy modes
722        #----------------------------------------------------------------------
723        if blk.get('_iso_occupancymode_label'):
724            modelist = []
725            shortmodelist = []
726            for lbl in blk.get('_iso_occupancymode_label'):
727                modelist.append(lbl)
728                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
729                # where SSSSS is the parent spacegroup, [x,y,z] is a location
730                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
731                # this extracts the AAAAA and BBBBB parts of the string
732                if regexp:
733                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
734                lbl = lbl.replace('order','o')
735                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
736            # read in the coordinate offset variables names and map them to G2 names/objects
737            occVarLbl = []
738            G2varLbl = []
739            G2varObj = []
740            error = False
741            for lbl in blk.get('_iso_deltaoccupancy_label'):
742                occVarLbl.append(lbl)
743                if '_' in lbl:
744                    albl = lbl[:lbl.rfind('_')]
745                    vlbl = lbl[lbl.rfind('_')+1:]
746                else:
747                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
748                    error = True
749                    continue
750                if albl not in atomlbllist:
751                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
752                    error = True
753                    continue
754                else:
755                    anum = atomlbllist.index(albl)
756                var = varLookup.get(vlbl)
757                if not var:
758                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
759                    error = True
760                    continue
761                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
762                G2varObj.append(G2obj.G2VarObj(
763                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
764                    ))
765            if error:
766                raise Exception("Error decoding variable labels")
767
768            if len(G2varObj) != len(modelist):
769                print ("non-square input")
770                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
771
772            error = False
773            ParentCoordinates = {}
774            for lbl,exp in zip(
775                blk.get('_iso_occupancy_label'),
776                blk.get('_iso_occupancy_formula'),
777                ):
778                if '_' in lbl:
779                    albl = lbl[:lbl.rfind('_')]
780                    vlbl = lbl[lbl.rfind('_')+1:]
781                else:
782                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
783                    error = True
784                    continue
785                if vlbl != 'occ':
786                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
787                    error = True
788                    continue
789                if '+' in exp:
790                    val = exp.split('+')[0].strip()
791                    val = G2p3.FormulaEval(val)
792                    if val is None:
793                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
794                        error = True
795                        continue
796                    ParentCoordinates[albl] = val
797            if error:
798                raise Exception("Error decoding occupancy labels")
799            # get mapping of modes to atomic coordinate displacements
800            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
801            for row,col,val in zip(
802                blk.get('_iso_occupancymodematrix_row'),
803                blk.get('_iso_occupancymodematrix_col'),
804                blk.get('_iso_occupancymodematrix_value'),):
805                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
806            # Invert to get mapping of atom displacements to modes
807            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
808            # create the constraints
809            for i,row in enumerate(occupancymodeInvmatrix):
810                constraint = []
811                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
812                    if k == 0: continue
813                    constraint.append([k,G2varObj[j]])
814                constraint += [shortmodelist[i],False,'f']
815                self.Constraints.append(constraint)
816            #----------------------------------------------------------------------
817            # save the ISODISTORT info for "mode analysis"
818            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
819            self.Phase['ISODISTORT'].update({
820                'OccModeList' : modelist,
821                'G2OccModeList' : shortmodelist,
822                'OccVarList' : occVarLbl,
823                'G2OccVarList' : G2varObj,
824                'BaseOcc' : ParentCoordinates,
825                'Var2OccMatrix' : occupancymodeInvmatrix,
826                'Occ2VarMatrix' : occupancymodematrix,
827                })
828            # make explaination dictionary
829            for mode,shortmode in zip(modelist,shortmodelist):
830                explaination[shortmode] = "ISODISTORT full name "+str(mode)
831        #----------------------------------------------------------------------
832        # done with read
833        #----------------------------------------------------------------------
834        if explaination: self.Constraints.append(explaination)
835
836        # # debug: show the mode var to mode relations
837        # for i,row in enumerate(displacivemodeInvmatrix):
838        #     l = ''
839        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
840        #         if k == 0: continue
841        #         if l: l += ' + '
842        #         #l += lbl+' * '+str(k)
843        #         l += G2varLbl[j]+' * '+str(k)
844        #     print str(i) + ': '+shortmodelist[i]+' = '+l
845        # print 70*'='
846
847        # # debug: Get the ISODISTORT offset values
848        # coordVarDelta = {}
849        # for lbl,val in zip(
850        #     blk.get('_iso_deltacoordinate_label'),
851        #     blk.get('_iso_deltacoordinate_value'),):
852        #     coordVarDelta[lbl] = float(val)
853        # modeVarDelta = {}
854        # for lbl,val in zip(
855        #     blk.get('_iso_displacivemode_label'),
856        #     blk.get('_iso_displacivemode_value'),):
857        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
858
859        # print 70*'='
860        # # compute the mode values from the reported coordinate deltas
861        # for i,row in enumerate(displacivemodeInvmatrix):
862        #     l = ''
863        #     sl = ''
864        #     s = 0.
865        #     for lbl,k in zip(coordVarLbl,row):
866        #         if k == 0: continue
867        #         if l: l += ' + '
868        #         l += lbl+' * '+str(k)
869        #         if sl: sl += ' + '
870        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
871        #         s += coordVarDelta[lbl] * k
872        #     print 'a'+str(i)+' = '+l
873        #     print '\t= '+sl
874        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
875        #     print
876
877        # print 70*'='
878        # # compute the coordinate displacements from the reported mode values
879        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
880        #     l = ''
881        #     sl = ''
882        #     s = 0.0
883        #     for j,k in enumerate(row):
884        #         if k == 0: continue
885        #         if l: l += ' + '
886        #         l += 'a'+str(j+1)+' * '+str(k)
887        #         if sl: sl += ' + '
888        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
889        #         s += modeVarDelta[modelist[j]] * k
890        #     print lbl+' = '+l
891        #     print '\t= '+sl
892        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
893        #     print
894
895        # determine the coordinate delta values from deviations from the parent structure
896        # for atmline in self.Phase['Atoms']:
897        #     lbl = atmline[0]
898        #     x,y,z = atmline[3:6]
899        #     if lbl not in ParentCoordinates:
900        #         print lbl,x,y,z
901        #         continue
902        #     px,py,pz = ParentCoordinates[lbl]
903        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.