source: trunk/imports/G2phase_CIF.py @ 3215

Last change on this file since 3215 was 3215, checked in by vondreele, 5 years ago

mag incommensurate stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 42.4 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-01-07 23:02:06 +0000 (Sun, 07 Jan 2018) $
4# $Author: vondreele $
5# $Revision: 3215 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 3215 2018-01-07 23:02:06Z 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: 3215 $")
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                    MSpGrp = sspgrp[0]
182                    self.MPhase['General']['Type'] = 'magnetic'
183                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
184                else:
185                    sspgrp = sspgrp.split('(')
186                    sspgrp[1] = "("+sspgrp[1]
187                    SpGrp = sspgrp[0]
188                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
189                    self.Phase['General']['Type'] = 'nuclear'
190                if not SpGrp:
191                    print (sspgrp)
192                    self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II'
193                    return False
194                SuperSg = sspgrp[1].replace('\\','')
195                SuperVec = [[0,0,.1],False,4]
196            else:   #not incommensurate
197                SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
198                if not SpGrp:
199                    SpGrp = blk.get("_space_group_name_H-M_alt",'')
200                if not SpGrp:   #try magnetic           
201                    MSpGrp = blk.get("_space_group.magn_name_BNS",'')
202                    if not MSpGrp:
203                        MSpGrp = blk.get("_space_group_magn.name_BNS",'')
204                        if not MSpGrp:
205                            msg = 'No recognizable space group name was found in the CIF.'
206                            self.errors = msg
207                            self.warnings += '\n'+msg
208                            return False
209                    SpGrp = blk.get('_parent_space_group.name_H-M_alt')
210                    if not SpGrp:
211                        SpGrp = blk.get('_parent_space_group.name_H-M')
212#                    SpGrp = MSpGrp.replace("'",'')
213                    SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','')   #get rid of screw '_'
214                    if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:]
215                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
216                    magnetic = True
217                    self.MPhase['General']['Type'] = 'magnetic'
218                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
219                    if not SpGrp:
220                        print (MSpGrp)
221                        self.warnings += 'No space group name was found in the CIF.'
222                        return False
223                else:
224                    SpGrp = SpGrp.replace('_','')
225                    self.Phase['General']['Type'] = 'nuclear'
226#process space group symbol
227            E,SGData = G2spc.SpcGroup(SpGrp)
228            if E and SpGrp:
229                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
230                if SpGrpNorm:
231                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
232            # nope, try the space group "out of the Box"
233            if E:
234                self.warnings += 'ERROR in space group symbol '+SpGrp
235                self.warnings += '\nThe space group has been set to "P 1". '
236                self.warnings += "Change this in phase's General tab."
237                self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
238                self.warnings += G2spc.SGErrors(E)
239                SGData = G2obj.P1SGData # P 1
240            self.Phase['General']['SGData'] = SGData
241
242            if magnetic and not Super:
243                SGData['SGFixed'] = True
244                try:
245                    sgoploop = blk.GetLoop('_space_group_symop_magn.id')
246                    sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
247                    opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
248                    centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]                   
249                except KeyError:        #old mag cif names
250                    sgoploop = blk.GetLoop('_space_group_symop.magn_id')
251                    sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
252                    opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
253                    centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
254                SGData['SGOps'] = []
255                SGData['SGCen'] = []
256                spnflp = []
257                for op in sgoploop:
258                    M,T,S = G2spc.MagText2MTS(op[opid])
259                    SGData['SGOps'].append([np.array(M,dtype=float),T])
260                    spnflp.append(S)
261                censpn = []
262                for cent in sgcenloop:
263                    M,C,S = G2spc.MagText2MTS(cent[centid])
264                    SGData['SGCen'].append(C)
265                    censpn += list(np.array(spnflp)*S)
266                self.MPhase['General']['SGData'] = SGData
267                self.MPhase['General']['SGData']['SpnFlp'] = censpn
268                self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
269                MagPtGp = blk.get('_space_group.magn_point_group')
270                if not MagPtGp:
271                    MagPtGp = blk.get('_space_group_magn.point_group_name')
272                self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
273#                GenSym,GenFlg = G2spc.GetGenSym(SGData)
274#                self.MPhase['General']['SGData']['GenSym'] = GenSym
275#                self.MPhase['General']['SGData']['GenFlg'] = GenFlg
276
277            if Super:
278                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
279                if E:
280                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
281                    self.warnings += '\n'+E
282                    SuperSg = SuperSg[:SuperSg.index(')')+1]
283                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
284                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
285                self.Phase['General']['SSGData'] = SSGData
286                if magnetic:
287                    self.MPhase['General']['SGData'] = SGData
288#                    self.MPhase['General']['SGData']['SpnFlp'] = censpn
289                    self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
290                    self.MPhase['General']['SGData']['MagPtGp'] = blk.get('_space_group.magn_point_group')
291                    self.MPhase['General']['SSGData'] = SSGData
292    #                GenSym,GenFlg = G2spc.GetGenSym(SGData)
293    #                self.MPhase['General']['SGData']['GenSym'] = GenSym
294    #                self.MPhase['General']['SGData']['GenFlg'] = GenFlg
295                   
296
297            # cell parameters
298            cell = []
299            for lbl in cellitems:
300                cell.append(cif.get_number_with_esd(blk[lbl])[0])
301            Volume = G2lat.calc_V(G2lat.cell2A(cell))
302            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
303            if magnetic:
304                self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]               
305            if Super:
306                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
307                waveDict = dict(waveloop.items())
308                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
309                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
310                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
311
312            # read in atoms
313            self.errors = 'Error during reading of atoms'
314            atomlbllist = [] # table to look up atom IDs
315            atomloop = blk.GetLoop('_atom_site_label')
316            atomkeys = [i.lower() for i in atomloop.keys()]
317            if not blk.get('_atom_site_type_symbol'):
318                self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
319            if magnetic:
320                try:
321                    magmoment = '_atom_site_moment.label'
322                    magatomloop = blk.GetLoop(magmoment)
323                    magatomkeys = [i.lower() for i in magatomloop.keys()]
324                    magatomlabels = blk.get(magmoment)
325                    G2MagDict = {'_atom_site_moment.label': 0,
326                                 '_atom_site_moment.crystalaxis_x':7,
327                                 '_atom_site_moment.crystalaxis_y':8,
328                                 '_atom_site_moment.crystalaxis_z':9}
329                except KeyError:
330                    magmoment = '_atom_site_moment_label'
331                    magatomloop = blk.GetLoop(magmoment)
332                    magatomkeys = [i.lower() for i in magatomloop.keys()]
333                    magatomlabels = blk.get(magmoment)
334                    G2MagDict = {'_atom_site_moment_label': 0,
335                                 '_atom_site_moment_crystalaxis_x':7,
336                                 '_atom_site_moment_crystalaxis_y':8,
337                                 '_atom_site_moment_crystalaxis_z':9}
338                   
339            if blk.get('_atom_site_aniso_label'):
340                anisoloop = blk.GetLoop('_atom_site_aniso_label')
341                anisokeys = [i.lower() for i in anisoloop.keys()]
342                anisolabels = blk.get('_atom_site_aniso_label')
343            else:
344                anisoloop = None
345                anisokeys = []
346                anisolabels = []
347            if Super:
348                occFloop = None
349                occCloop = None
350                occFdict = {}
351                occCdict = {}
352                displSloop = None
353                displFloop = None
354                MagFloop = None
355                displSdict = {}
356                displFdict = {}
357                MagFdict = {}
358                UijFloop = None
359                UijFdict = {}
360                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
361                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
362                    occFdict = dict(occFloop.items())
363                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
364                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
365                    occCdict = dict(occCloop.items())
366                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
367                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
368                    displFdict = dict(displFloop.items())                           
369                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
370                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
371                    displSdict = dict(displSloop.items())
372                if blk.get('_atom_site_U_Fourier_atom_site_label'):
373                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
374                    UijFdict = dict(UijFloop.items())
375                if blk.get('_atom_site_moment_Fourier_atom_site_label'):
376                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
377                    MagFdict = dict(MagFloop.items())
378                    Mnames =  ['_atom_site_moment_fourier_atom_site_label',
379                               '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
380                               '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
381                elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
382                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
383                    MagFdict = dict(MagFloop.items())                           
384                    Mnames =  ['_atom_site_moment_fourier.atom_site_label',
385                               '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
386                               '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
387            self.Phase['Atoms'] = []
388            if magnetic:
389                self.MPhase['Atoms'] = []
390            G2AtomDict = {  '_atom_site_type_symbol' : 1,
391                            '_atom_site_label' : 0,
392                            '_atom_site_fract_x' : 3,
393                            '_atom_site_fract_y' : 4,
394                            '_atom_site_fract_z' : 5,
395                            '_atom_site_occupancy' : 6,
396                            '_atom_site_aniso_u_11' : 11,
397                            '_atom_site_aniso_u_22' : 12,
398                            '_atom_site_aniso_u_33' : 13,
399                            '_atom_site_aniso_u_12' : 14,
400                            '_atom_site_aniso_u_13' : 15,
401                            '_atom_site_aniso_u_23' : 16, }
402
403            ranIdlookup = {}
404            for aitem in atomloop:
405                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
406                for val,key in zip(aitem,atomkeys):
407                    col = G2AtomDict.get(key,-1)
408                    if col >= 3:
409                        atomlist[col] = cif.get_number_with_esd(val)[0]
410                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
411                    elif col is not None and col != -1:
412                        atomlist[col] = val
413                    elif key in ('_atom_site_thermal_displace_type',
414                               '_atom_site_adp_type'):   #Iso or Aniso?
415                        if val.lower() == 'uani':
416                            atomlist[9] = 'A'
417                    elif key == '_atom_site_u_iso_or_equiv':
418                        uisoval = cif.get_number_with_esd(val)[0]
419                        if uisoval is not None: 
420                            atomlist[10] = uisoval
421                if not atomlist[1] and atomlist[0]:
422                    typ = atomlist[0].rstrip('0123456789-+')
423                    if G2elem.CheckElement(typ):
424                        atomlist[1] = typ
425                    if not atomlist[1]:
426                        atomlist[1] = 'Xe'
427                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
428                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
429                    atomlist[9] = 'A'
430                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
431                        col = G2AtomDict.get(key)
432                        if col:
433                            atomlist[col] = cif.get_number_with_esd(val)[0]
434                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
435                atomlist[1] = G2elem.FixValence(atomlist[1])
436                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
437                self.Phase['Atoms'].append(atomlist)
438                ranIdlookup[atomlist[0]] = atomlist[-1]
439                if atomlist[0] in atomlbllist:
440                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
441                else:
442                    atomlbllist.append(atomlist[0])
443
444                if magnetic and atomlist[0] in magatomlabels:
445                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
446                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
447                        mcol = G2MagDict.get(mkey,-1)
448                        if mcol:
449                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
450                    self.MPhase['Atoms'].append(matomlist)
451                if Super:
452                    Sfrac = []
453                    Sadp = np.zeros((4,12))
454                    Spos = np.zeros((4,6))
455                    Smag = np.zeros((4,6))
456                    nim = -1
457                    waveType = 'Fourier'                               
458                    if displFdict:
459                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
460                            if item == atomlist[0]:
461                                waveType = 'Fourier'                               
462                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
463                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
464                                if im != nim:
465                                    nim = im
466                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
467                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
468                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
469                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
470                    if nim >= 0:
471                        Spos = [[spos,False] for spos in Spos[:nim]]
472                    else:
473                        Spos = []
474                    if UijFdict:
475                        nim = -1
476                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
477                            if item == atomlist[0]:
478                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
479                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
480                                if im != nim:
481                                    nim = im
482                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
483                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
484                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
485                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
486                    if nim >= 0:
487                        Sadp = [[sadp,False] for sadp in Sadp[:nim]]
488                    else:
489                        Sadp = []
490                    if MagFdict:
491                        nim = -1
492                        for i,item in enumerate(MagFdict[Mnames[0]]):
493                            if item == atomlist[0]:
494                                waveType = 'Fourier'                               
495                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
496                                im = int(MagFdict[Mnames[2]][i])
497                                if im != nim:
498                                    nim = im
499                                val = MagFdict[Mnames[3]][i]
500                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
501                                val = MagFdict[Mnames[4]][i]
502                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
503                        if nim >= 0:
504                            Smag = [[smag,False] for smag in Smag[:nim]]
505                        else:
506                            Smag = []
507                    SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
508                    if magnetic and atomlist[0] in magatomlabels:
509                        matomlist.append(SSdict)
510                    atomlist.append(SSdict)
511            if len(atomlbllist) != len(self.Phase['Atoms']):
512                self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
513            for lbl in phasenamefields: # get a name for the phase
514                name = blk.get(lbl)
515                if name is None:
516                    continue
517                name = name.strip()
518                if name == '?' or name == '.':
519                    continue
520                else:
521                    break
522            else: # no name found, use block name for lack of a better choice
523                name = blknm
524            self.Phase['General']['Name'] = name.strip()[:20]
525            self.Phase['General']['Super'] = Super
526            if magnetic:
527                self.MPhase['General']['Type'] = 'magnetic'               
528                self.MPhase['General']['Name'] = name.strip()[:20]+' mag'
529                self.MPhase['General']['Super'] = Super
530                if Super:
531                    self.MPhase['General']['Modulated'] = True
532                    self.MPhase['General']['SuperVec'] = SuperVec
533                    self.MPhase['General']['SuperSg'] = SuperSg
534                    self.MPhase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
535            else:
536                self.MPhase = None
537            if Super:
538                self.Phase['General']['Modulated'] = True
539                self.Phase['General']['SuperVec'] = SuperVec
540                self.Phase['General']['SuperSg'] = SuperSg
541                self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
542            if not self.isodistort_warnings:
543                if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
544                    self.errors = "Error while processing ISODISTORT constraints"
545                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
546            else:
547                self.warnings += self.isodistort_warnings
548            returnstat = True
549        return returnstat
550
551    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
552        'Process ISODISTORT items to create constraints etc.'
553        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
554        'Maps ISODISTORT parm names to GSAS-II names'
555        #----------------------------------------------------------------------
556        # read in the ISODISTORT displacement modes
557        #----------------------------------------------------------------------
558        self.Constraints = []
559        explaination = {}
560        if blk.get('_iso_displacivemode_label'):
561            modelist = []
562            shortmodelist = []
563            for lbl in blk.get('_iso_displacivemode_label'):
564                modelist.append(lbl)
565                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
566                # where SSSSS is the parent spacegroup, [x,y,z] is a location
567                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
568                # this extracts the AAAAA and BBBBB parts of the string
569                if regexp:
570                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
571                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
572            # read in the coordinate offset variables names and map them to G2 names/objects
573            coordVarLbl = []
574            G2varLbl = []
575            G2varObj = []
576            error = False
577            for lbl in blk.get('_iso_deltacoordinate_label'):
578                coordVarLbl.append(lbl)
579                if '_' in lbl:
580                    albl = lbl[:lbl.rfind('_')]
581                    vlbl = lbl[lbl.rfind('_')+1:]
582                else:
583                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
584                    error = True
585                    continue
586                if albl not in atomlbllist:
587                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
588                    error = True
589                    continue
590                else:
591                    anum = atomlbllist.index(albl)
592                var = varLookup.get(vlbl)
593                if not var:
594                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
595                    error = True
596                    continue
597                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
598                G2varObj.append(G2obj.G2VarObj(
599                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
600                    ))
601            if error:
602                raise Exception("Error decoding variable labels")
603
604            if len(G2varObj) != len(modelist):
605                print ("non-square input")
606                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
607
608            error = False
609            ParentCoordinates = {}
610            for lbl,exp in zip(
611                blk.get('_iso_coordinate_label'),
612                blk.get('_iso_coordinate_formula'),
613                ):
614                if '_' in lbl:
615                    albl = lbl[:lbl.rfind('_')]
616                    vlbl = lbl[lbl.rfind('_')+1:]
617                else:
618                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
619                    error = True
620                    continue
621                if vlbl not in 'xyz' or len(vlbl) != 1:
622                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
623                    error = True
624                    continue
625                i = 'xyz'.index(vlbl)
626                if not ParentCoordinates.get(albl):
627                    ParentCoordinates[albl] = [None,None,None]
628                if '+' in exp:
629                    val = exp.split('+')[0].strip()
630                    val = G2p3.FormulaEval(val)
631                    if val is None:
632                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
633                        error = True
634                        continue
635                    ParentCoordinates[albl][i] = val
636                else:
637                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
638            if error:
639                print (self.warnings)
640                raise Exception("Error decoding variable labels")
641            # get mapping of modes to atomic coordinate displacements
642            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
643            for row,col,val in zip(
644                blk.get('_iso_displacivemodematrix_row'),
645                blk.get('_iso_displacivemodematrix_col'),
646                blk.get('_iso_displacivemodematrix_value'),):
647                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
648            # Invert to get mapping of atom displacements to modes
649            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
650            # create the constraints
651            for i,row in enumerate(displacivemodeInvmatrix):
652                constraint = []
653                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
654                    if k == 0: continue
655                    constraint.append([k,G2varObj[j]])
656                constraint += [shortmodelist[i],False,'f']
657                self.Constraints.append(constraint)
658            #----------------------------------------------------------------------
659            # save the ISODISTORT info for "mode analysis"
660            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
661            self.Phase['ISODISTORT'].update({
662                'IsoModeList' : modelist,
663                'G2ModeList' : shortmodelist,
664                'IsoVarList' : coordVarLbl,
665                'G2VarList' : G2varObj,
666                'ParentStructure' : ParentCoordinates,
667                'Var2ModeMatrix' : displacivemodeInvmatrix,
668                'Mode2VarMatrix' : displacivemodematrix,
669                })
670            # make explaination dictionary
671            for mode,shortmode in zip(modelist,shortmodelist):
672                explaination[shortmode] = "ISODISTORT full name "+str(mode)
673        #----------------------------------------------------------------------
674        # now read in the ISODISTORT occupancy modes
675        #----------------------------------------------------------------------
676        if blk.get('_iso_occupancymode_label'):
677            modelist = []
678            shortmodelist = []
679            for lbl in blk.get('_iso_occupancymode_label'):
680                modelist.append(lbl)
681                # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
682                # where SSSSS is the parent spacegroup, [x,y,z] is a location
683                regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
684                # this extracts the AAAAA and BBBBB parts of the string
685                if regexp:
686                    lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
687                lbl = lbl.replace('order','o')
688                G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
689            # read in the coordinate offset variables names and map them to G2 names/objects
690            occVarLbl = []
691            G2varLbl = []
692            G2varObj = []
693            error = False
694            for lbl in blk.get('_iso_deltaoccupancy_label'):
695                occVarLbl.append(lbl)
696                if '_' in lbl:
697                    albl = lbl[:lbl.rfind('_')]
698                    vlbl = lbl[lbl.rfind('_')+1:]
699                else:
700                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
701                    error = True
702                    continue
703                if albl not in atomlbllist:
704                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
705                    error = True
706                    continue
707                else:
708                    anum = atomlbllist.index(albl)
709                var = varLookup.get(vlbl)
710                if not var:
711                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
712                    error = True
713                    continue
714                G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
715                G2varObj.append(G2obj.G2VarObj(
716                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
717                    ))
718            if error:
719                raise Exception("Error decoding variable labels")
720
721            if len(G2varObj) != len(modelist):
722                print ("non-square input")
723                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
724
725            error = False
726            ParentCoordinates = {}
727            for lbl,exp in zip(
728                blk.get('_iso_occupancy_label'),
729                blk.get('_iso_occupancy_formula'),
730                ):
731                if '_' in lbl:
732                    albl = lbl[:lbl.rfind('_')]
733                    vlbl = lbl[lbl.rfind('_')+1:]
734                else:
735                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
736                    error = True
737                    continue
738                if vlbl != 'occ':
739                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
740                    error = True
741                    continue
742                if '+' in exp:
743                    val = exp.split('+')[0].strip()
744                    val = G2p3.FormulaEval(val)
745                    if val is None:
746                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
747                        error = True
748                        continue
749                    ParentCoordinates[albl] = val
750            if error:
751                raise Exception("Error decoding occupancy labels")
752            # get mapping of modes to atomic coordinate displacements
753            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
754            for row,col,val in zip(
755                blk.get('_iso_occupancymodematrix_row'),
756                blk.get('_iso_occupancymodematrix_col'),
757                blk.get('_iso_occupancymodematrix_value'),):
758                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
759            # Invert to get mapping of atom displacements to modes
760            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
761            # create the constraints
762            for i,row in enumerate(occupancymodeInvmatrix):
763                constraint = []
764                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
765                    if k == 0: continue
766                    constraint.append([k,G2varObj[j]])
767                constraint += [shortmodelist[i],False,'f']
768                self.Constraints.append(constraint)
769            #----------------------------------------------------------------------
770            # save the ISODISTORT info for "mode analysis"
771            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
772            self.Phase['ISODISTORT'].update({
773                'OccModeList' : modelist,
774                'G2OccModeList' : shortmodelist,
775                'OccVarList' : occVarLbl,
776                'G2OccVarList' : G2varObj,
777                'BaseOcc' : ParentCoordinates,
778                'Var2OccMatrix' : occupancymodeInvmatrix,
779                'Occ2VarMatrix' : occupancymodematrix,
780                })
781            # make explaination dictionary
782            for mode,shortmode in zip(modelist,shortmodelist):
783                explaination[shortmode] = "ISODISTORT full name "+str(mode)
784        #----------------------------------------------------------------------
785        # done with read
786        #----------------------------------------------------------------------
787        if explaination: self.Constraints.append(explaination)
788
789        # # debug: show the mode var to mode relations
790        # for i,row in enumerate(displacivemodeInvmatrix):
791        #     l = ''
792        #     for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
793        #         if k == 0: continue
794        #         if l: l += ' + '
795        #         #l += lbl+' * '+str(k)
796        #         l += G2varLbl[j]+' * '+str(k)
797        #     print str(i) + ': '+shortmodelist[i]+' = '+l
798        # print 70*'='
799
800        # # debug: Get the ISODISTORT offset values
801        # coordVarDelta = {}
802        # for lbl,val in zip(
803        #     blk.get('_iso_deltacoordinate_label'),
804        #     blk.get('_iso_deltacoordinate_value'),):
805        #     coordVarDelta[lbl] = float(val)
806        # modeVarDelta = {}
807        # for lbl,val in zip(
808        #     blk.get('_iso_displacivemode_label'),
809        #     blk.get('_iso_displacivemode_value'),):
810        #     modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
811
812        # print 70*'='
813        # # compute the mode values from the reported coordinate deltas
814        # for i,row in enumerate(displacivemodeInvmatrix):
815        #     l = ''
816        #     sl = ''
817        #     s = 0.
818        #     for lbl,k in zip(coordVarLbl,row):
819        #         if k == 0: continue
820        #         if l: l += ' + '
821        #         l += lbl+' * '+str(k)
822        #         if sl: sl += ' + '
823        #         sl += str(coordVarDelta[lbl])+' * '+str(k)
824        #         s += coordVarDelta[lbl] * k
825        #     print 'a'+str(i)+' = '+l
826        #     print '\t= '+sl
827        #     print  modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
828        #     print
829
830        # print 70*'='
831        # # compute the coordinate displacements from the reported mode values
832        # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
833        #     l = ''
834        #     sl = ''
835        #     s = 0.0
836        #     for j,k in enumerate(row):
837        #         if k == 0: continue
838        #         if l: l += ' + '
839        #         l += 'a'+str(j+1)+' * '+str(k)
840        #         if sl: sl += ' + '
841        #         sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
842        #         s += modeVarDelta[modelist[j]] * k
843        #     print lbl+' = '+l
844        #     print '\t= '+sl
845        #     print lbl,G2varLbl[i],coordVarDelta[lbl],s
846        #     print
847
848        # determine the coordinate delta values from deviations from the parent structure
849        # for atmline in self.Phase['Atoms']:
850        #     lbl = atmline[0]
851        #     x,y,z = atmline[3:6]
852        #     if lbl not in ParentCoordinates:
853        #         print lbl,x,y,z
854        #         continue
855        #     px,py,pz = ParentCoordinates[lbl]
856        #     print lbl,x,y,z,x-px,y-py,z-pz
Note: See TracBrowser for help on using the repository browser.