source: trunk/imports/G2phase_CIF.py @ 3736

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

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