source: trunk/imports/G2phase.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 2 years ago

multiple small changes to allow docs build under 3.x

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 32.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2019-12-20 05:04:03 +0000 (Fri, 20 Dec 2019) $
4# $Author: toby $
5# $Revision: 4213 $
6# $URL: trunk/imports/G2phase.py $
7# $Id: G2phase.py 4213 2019-12-20 05:04:03Z toby $
8########### SVN repository information ###################
9#
10'''
11*Module G2phase: PDB, .EXP & JANA m40,m50*
12-------------------------------------------
13
14A set of short routines to read in phases using routines that were
15previously implemented in GSAS-II: PDB, GSAS .EXP and JANA m40-m50 file formats
16
17'''
18
19from __future__ import division, print_function
20import sys
21import os.path
22import math
23import random as ran
24import numpy as np
25try:
26    import wx
27except ImportError:
28    wx = None
29import GSASIIobj as G2obj
30import GSASIIspc as G2spc
31import GSASIIlattice as G2lat
32import GSASIIpath
33GSASIIpath.SetVersionNumber("$Revision: 4213 $")
34try:  # fails on doc build
35    R2pisq = 1./(2.*np.pi**2)
36except TypeError:
37    pass
38
39class PDB_ReaderClass(G2obj.ImportPhase):
40    'Routine to import Phase information from a PDB file'
41    def __init__(self):
42        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
43            extensionlist=('.pdb','.ent','.PDB','.ENT'),
44            strictExtension=True,
45            formatName = 'PDB',
46            longFormatName = 'Original Protein Data Bank (.pdb file) import'
47            )
48    def ContentsValidator(self, filename):
49        '''Taking a stab a validating a PDB file
50        (look for cell & at least one atom)
51        '''
52        fp = open(filename,'r')
53#        for i,l in enumerate(fp):
54#            if l.startswith('CRYST1'):
55#                break
56#        else:
57#            self.errors = 'no CRYST1 record found'
58#            fp.close()
59#            return False
60        for i,l in enumerate(fp):
61            if l.startswith('ATOM'):
62                fp.close()
63                return True
64        self.errors = 'no ATOM records found after CRYST1 record'
65        fp.close()
66        return False
67
68    def Reader(self,filename, ParentFrame=None, **unused):
69        'Read a PDF file using :meth:`ReadPDBPhase`'
70        self.Phase = self.ReadPDBPhase(filename, ParentFrame)
71        return True
72       
73    def ReadPDBPhase(self,filename,parent=None):
74        '''Read a phase from a PDB file.
75        '''
76        EightPiSq = 8.*math.pi**2
77        self.errors = 'Error opening file'
78        file = open(filename, 'Ur')
79        Phase = {}
80        Title = os.path.basename(filename)
81       
82        Compnd = ''
83        Atoms = []
84        A = np.zeros(shape=(3,3))
85        S = file.readline()
86        line = 1
87        SGData = None
88        cell = None
89        Dummy = True
90        while S:
91            self.errors = 'Error reading at line '+str(line)
92            Atom = []
93            if 'TITLE' in S[:5]:
94                Title = S[10:72].strip()
95            elif 'COMPND    ' in S[:10]:
96                Compnd = S[10:72].strip()
97            elif 'CRYST' in S[:5]:
98                Dummy = False
99                abc = S[7:34].split()
100                angles = S[34:55].split()
101                cell=[float(abc[0]),float(abc[1]),float(abc[2]),
102                    float(angles[0]),float(angles[1]),float(angles[2])]
103                Volume = G2lat.calc_V(G2lat.cell2A(cell))
104                AA,AB = G2lat.cell2AB(cell)
105                SpGrp = S[55:65]
106                E,SGData = G2spc.SpcGroup(SpGrp)
107                # space group processing failed, try to look up name in table
108                if E:
109                    SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
110                    if SpGrpNorm:
111                        E,SGData = G2spc.SpcGroup(SpGrpNorm)
112                while E:
113                    dlg = wx.TextEntryDialog(parent,
114                        SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
115                        'ERROR in space group symbol','',style=wx.OK)
116                    if dlg.ShowModal() == wx.ID_OK:
117                        SpGrp = dlg.GetValue()
118                        E,SGData = G2spc.SpcGroup(SpGrp)
119                    else:
120                        SGData = G2obj.P1SGData # P 1
121                        self.warnings += '\nThe space group was not interpreted and has been set to "P 1".'
122                        self.warnings += "Change this in phase's General tab."           
123                    dlg.Destroy()
124#                SGlines = G2spc.SGPrint(SGData)
125#                for l in SGlines: print (l)
126            elif 'SCALE' in S[:5]:
127                V = S[10:41].split()
128                A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
129            elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
130                if not SGData:
131                    self.warnings += '\nThe space group was not read before atoms and has been set to "P 1". '
132                    self.warnings += "Change this in phase's General tab."
133                    SGData = G2obj.P1SGData # P 1
134                    cell = [20.0,20.0,20.0,90.,90.,90.]
135                    Volume = G2lat.calc_V(G2lat.cell2A(cell))
136                    AA,AB = G2lat.cell2AB(cell)                   
137                XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
138                XYZ = np.inner(AB,XYZ)
139                XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
140                SytSym,Mult = G2spc.SytSym(XYZ,SGData)[:2]
141                Uiso = float(S[61:67])/EightPiSq
142                Type = S[76:78].lower()
143                if Dummy and S[12:17].strip() == 'CA':
144                    Type = 'C'
145                Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
146                    S[12:17].strip(),Type.strip().capitalize(),'',XYZ[0],XYZ[1],XYZ[2],
147                    float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
148                if S[16] in [' ','A','B']:
149                    Atom[3] = Atom[3][:3]
150                    Atom.append(ran.randint(0,sys.maxsize))
151                    Atoms.append(Atom)
152            elif 'ANISOU' in S[:6]:
153                Uij = S[30:72].split()
154                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
155                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
156                Atoms[-1] = Atoms[-1][:14]+Uij
157                Atoms[-1][12] = 'A'
158                Atoms[-1].append(ran.randint(0,sys.maxsize))
159            S = file.readline()
160            line += 1
161        file.close()
162        self.errors = 'Error after read complete'
163        if Title:
164            PhaseName = Title
165        elif Compnd:
166            PhaseName = Compnd
167        else:
168            PhaseName = 'None'
169        if not SGData:
170            raise self.ImportException("No space group (CRYST entry) found")
171        if not cell:
172            raise self.ImportException("No cell (CRYST entry) found")
173        Phase = G2obj.SetNewPhase(Name=PhaseName,SGData=SGData,cell=cell+[Volume,])
174        Phase['General']['Type'] = 'macromolecular'
175        Phase['General']['AtomPtrs'] = [6,4,10,12]
176        Phase['Atoms'] = Atoms
177        return Phase
178
179class EXP_ReaderClass(G2obj.ImportPhase):
180    'Routine to import Phase information from GSAS .EXP files'
181    def __init__(self):
182        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
183            extensionlist=('.EXP','.exp'),
184            strictExtension=True,
185            formatName = 'GSAS .EXP',
186            longFormatName = 'GSAS Experiment (.EXP file) import'
187            )
188       
189    def ContentsValidator(self, filename):
190        'Look for a VERSION tag in 1st line' 
191        fp = open(filename,'r')
192        if fp.read(13) == '     VERSION ':
193            fp.close()
194            return True
195        self.errors = 'File does not begin with VERSION tag'
196        fp.close()
197        return False
198
199    def Reader(self,filename,ParentFrame=None,usedRanIdList=[],**unused):
200        'Read a phase from a GSAS .EXP file using :meth:`ReadEXPPhase`'
201        self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
202        while self.Phase['ranId'] in usedRanIdList:
203            self.Phase['ranId'] = ran.randint(0,sys.maxsize)
204        # make sure the ranId is really unique!
205        self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
206        while self.MPhase['ranId'] in usedRanIdList:
207            self.MPhase['ranId'] = ran.randint(0,sys.maxsize)
208        fp = open(filename,'r')
209        self.ReadEXPPhase(ParentFrame, fp)
210        fp.close()
211        return True
212
213    def ReadEXPPhase(self, G2frame,filepointer):
214        '''Read a phase from a GSAS .EXP file.
215        '''
216        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
217        textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
218            'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
219            'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
220        shNcof = 0
221        S = 1
222        NPhas = []
223        Expr = [{},{},{},{},{},{},{},{},{}] # GSAS can have at most 9 phases
224        for line,S in enumerate(filepointer):
225            self.errors = 'Error reading at line '+str(line+1)
226            if 'EXPR NPHAS' in S[:12]:
227                NPhas = S[12:-1].split()
228            if 'CRS' in S[:3]:
229                N = int(S[3:4])-1
230                Expr[N][S[:12]] = S[12:-1]
231        PNames = []
232        if not NPhas:
233            raise self.ImportException("No EXPR NPHAS record found")
234        self.errors = 'Error interpreting file'
235        for n,N in enumerate(NPhas):
236            if N != '0':
237                result = n
238                key = 'CRS'+str(n+1)+'    PNAM'
239                PNames.append(Expr[n][key])
240        if len(PNames) == 0:
241            raise self.ImportException("No phases found")           
242        elif len(PNames) > 1:
243            dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
244            try:
245                if dlg.ShowModal() == wx.ID_OK:
246                    result = dlg.GetSelection() # I think this breaks is there are skipped phases. Cant this happen?
247            finally:
248                dlg.Destroy()       
249        EXPphase = Expr[result]
250        keyList = list(EXPphase.keys())
251        keyList.sort()
252        SGData = {}
253        MPtype = ''
254        if NPhas[result] == '1':
255            Ptype = 'nuclear'
256        elif NPhas[result] =='2':
257            Ptype = 'nuclear'
258            MPtype = 'magnetic'
259            MagDmin = 1.0
260        elif NPhas[result] =='3':
261            Ptype = 'magnetic'
262            MagDmin = 1.0
263        elif NPhas[result] == '4':
264            Ptype = 'macromolecular'
265        elif NPhas[result] == '10':
266            Ptype = 'Pawley'
267        else:
268            raise self.ImportException("Phase type not recognized") 
269           
270        for key in keyList:
271            if 'PNAM' in key:
272               PhaseName = EXPphase[key].strip()
273            elif 'ABC   ' in key:
274                abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
275            elif 'ANGLES' in key:
276                angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
277            elif 'SG SYM' in key:
278                SpGrp = EXPphase[key][:15].strip()
279                E,SGData = G2spc.SpcGroup(SpGrp)
280                if E:
281                    SGData = G2obj.P1SGData # P 1 -- unlikely to need this!
282                    self.warnings += '\nThe GSAS space group was not interpreted(!) and has been set to "P 1".'
283                    self.warnings += "Change this in phase's General tab."                       
284            elif 'SPNFLP' in key:
285                SpnFlp = np.array([int(float(s)) for s in EXPphase[key].split()])
286                SpnFlp = np.where(SpnFlp==0,1,SpnFlp)
287                SpnFlp = [1,]+list(SpnFlp)
288                if SGData['SpGrp'][0] in ['A','B','C','I','R','F']:
289                    SpnFlp = list(SpnFlp)+[1,1,1,1]
290            elif 'MXDSTR' in key:
291                MagDmin = float(EXPphase[key][:10])               
292            elif 'OD    ' in key:
293                SHdata = EXPphase[key].split() # may not have all 9 values
294                SHvals = 9*[0]
295                for i in range(9):
296                    try:
297                        float(SHdata[i])
298                        SHvals[i] = SHdata[i]
299                    except:
300                        pass
301                textureData['Order'] = int(SHvals[0])
302                textureData['Model'] = shModels[int(SHvals[2])]
303                textureData['Sample omega'] = [False,float(SHvals[6])]
304                textureData['Sample chi'] = [False,float(SHvals[7])]
305                textureData['Sample phi'] = [False,float(SHvals[8])]
306                shNcof = int(SHvals[1])
307        Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
308
309        Atoms = []
310        MAtoms = []
311        Bmat = G2lat.cell2AB(abc+angles)[1]
312        if Ptype == 'macromolecular':
313            for key in keyList:
314                if 'AT' in key[6:8]:
315                    S = EXPphase[key]
316                    Atom = [S[56:60].strip(),S[50:54].strip().upper(),S[54:56],
317                        S[46:51].strip(),S[:8].strip().capitalize(),'',
318                        float(S[16:24]),float(S[24:32]),float(S[32:40]),
319                        float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
320                    XYZ = Atom[6:9]
321                    Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)[:2]
322                    Atom.append(ran.randint(0,sys.maxsize))
323                    Atoms.append(Atom)
324        else:
325            for key in keyList:
326                if 'AT' in key:
327                    if key[11:] == 'A':
328                        S = EXPphase[key]
329                    elif key[11:] == 'B':
330                        S1 = EXPphase[key]
331                        Atom = [S[50:58].strip(),S[:10].strip().capitalize(),'',
332                            float(S[10:20]),float(S[20:30]),float(S[30:40]),
333                            float(S[40:50]),'',int(S[60:62]),S1[62:63]]
334                            #float(S[40:50]),'',int(S[60:62]),S1[130:131]]
335                        if Atom[9] == 'I':
336                            Atom += [float(S1[0:10]),0.,0.,0.,0.,0.,0.]
337                        elif Atom[9] == 'A':
338                            Atom += [0.0,
339                                float(S1[ 0:10]),float(S1[10:20]),
340                                float(S1[20:30]),float(S1[30:40]),
341                                float(S1[40:50]),float(S1[50:60])]
342                        else:
343                            print('Error in line with key: '+key)
344                            Atom += [0.,0.,0.,0.,0.,0.,0.]
345                        XYZ = Atom[3:6]
346                        Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)[:2]
347                        Atom.append(ran.randint(0,sys.maxsize))
348                        Atoms.append(Atom)
349                    elif key[11:] == 'M' and key[6:8] == 'AT':
350                        S = EXPphase[key]
351                        mom = np.array([float(S[:10]),float(S[10:20]),float(S[20:30])])
352                        mag = np.sqrt(np.sum(mom**2))
353                        mom = np.inner(Bmat,mom)*mag
354                        MAtoms.append(Atom)
355                        MAtoms[-1] = Atom[:7]+list(mom)+Atom[7:]
356                       
357        if shNcof:
358            shCoef = {}
359            nRec = [i+1 for i in range((shNcof-1)//6+1)]
360            for irec in nRec:
361                ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
362                indx = EXPphase[ODkey].split()
363                ODkey = ODkey[:-1]+'B'
364                vals = EXPphase[ODkey].split()
365                for i,val in enumerate(vals):
366                    key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
367                    shCoef[key] = float(val)
368            textureData['SH Coeff'] = [False,shCoef]
369           
370        if not SGData:
371            raise self.ImportException("No space group found in phase")
372        if not abc:
373            raise self.ImportException("No cell lengths found in phase")
374        if not angles:
375            raise self.ImportException("No cell angles found in phase")
376        if not Atoms:
377            raise self.ImportException("No atoms found in phase")
378           
379        self.Phase['General'].update({'Type':Ptype,'Name':PhaseName,'Cell':[False,]+abc+angles+[Volume,],'SGData':SGData})
380        if MPtype == 'magnetic':
381            self.MPhase['General'].update({'Type':'magnetic','Name':PhaseName+' mag','Cell':[False,]+abc+angles+[Volume,],'SGData':SGData})
382        else:
383            self.MPhase = None
384           
385        if Ptype =='macromolecular':
386            self.Phase['General']['AtomPtrs'] = [6,4,10,12]
387            self.Phase['Atoms'] = Atoms
388        elif Ptype == 'magnetic':
389            self.Phase['General']['AtomPtrs'] = [3,1,10,12]
390            self.Phase['General']['SGData']['SGSpin'] = SpnFlp
391            self.Phase['General']['MagDmin'] = MagDmin
392            self.Phase['Atoms'] = MAtoms           
393        else:   #nuclear
394            self.Phase['General']['AtomPtrs'] = [3,1,7,9]   
395            self.Phase['General']['SH Texture'] = textureData
396            self.Phase['Atoms'] = Atoms
397        if MPtype =='magnetic':
398            self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
399            self.MPhase['General']['SGData']['SGSpin'] = SpnFlp
400            self.MPhase['General']['MagDmin'] = MagDmin
401            self.MPhase['Atoms'] = MAtoms
402
403class JANA_ReaderClass(G2obj.ImportPhase):
404    'Routine to import Phase information from a JANA2006 file'
405    def __init__(self):
406        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
407            extensionlist=('.m50','.M50'),
408            strictExtension=True,
409            formatName = 'JANA m50',
410            longFormatName = 'JANA2006 phase import'
411            )
412    def ContentsValidator(self, filename):
413        '''Taking a stab a validating a .m50 file
414        (look for cell & at least one atom)
415        '''
416        fp = open(filename,'r')
417        for i,l in enumerate(fp):
418            if l.startswith('cell'):
419                break
420        else:
421            self.errors = 'no cell record found'
422            fp.close()
423            return False
424        for i,l in enumerate(fp):
425            if l.startswith('spgroup'):
426                fp.close()
427                return True
428        self.errors = 'no spgroup record found after cell record'
429        fp.close()
430        return False
431       
432    def Reader(self,filename, ParentFrame=None, **unused):
433        'Read a m50 file using :meth:`ReadJANAPhase`'
434        self.Phase = self.ReadJANAPhase(filename, ParentFrame)
435        return True
436       
437    def ReadJANAPhase(self,filename,parent=None):
438        '''Read a phase from a JANA2006 m50 & m40 files.
439        '''
440        self.errors = 'Error opening file'
441        fp = open(filename, 'Ur') #contains only cell & spcgroup
442        Phase = {}
443        Title = os.path.basename(filename)
444        Type = 'nuclear'
445        Atoms = []
446        Atypes = []
447        SuperVec = [[0,0,.1],False,4]
448        S = fp.readline()
449        line = 1
450        SGData = None
451        SuperSg = ''
452        cell = None
453        nqi = 0
454        version = '2000'
455        while S:
456            self.errors = 'Error reading at line '+str(line)
457            if 'title' in S and S != 'title\n':
458                Title = S.split()[1]
459            elif 'Jana2006' in S:
460                self.warnings += '\nJana2006 file detected'
461                version = '2006'
462            elif 'cell' in S[:4]:
463                cell = S[5:].split()
464                cell=[float(cell[0]),float(cell[1]),float(cell[2]),
465                    float(cell[3]),float(cell[4]),float(cell[5])]
466                Volume = G2lat.calc_V(G2lat.cell2A(cell))
467                G,g = G2lat.cell2Gmat(cell)
468                ast = np.sqrt(np.diag(G))
469                Mast = np.multiply.outer(ast,ast)   
470               
471            elif 'spgroup' in S:
472                if 'X' in S:
473                    raise self.ImportException("Ad hoc Supersymmetry centering "+S+" not allowed in GSAS-II")           
474                SpGrp = S.split()[1]
475                SuperSg = ''
476                if '(' in SpGrp:    #supercell symmetry - split in 2
477                    SuperStr = SpGrp.split('(')
478                    SpGrp = SuperStr[0]
479                    SuperSg = '('+SuperStr[1]
480                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
481                E,SGData = G2spc.SpcGroup(SpGrpNorm)
482                # space group processing failed, try to look up name in table
483                while E:
484                    print (G2spc.SGErrors(E))
485                    dlg = wx.TextEntryDialog(parent,
486                        SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
487                        'ERROR in space group symbol','',style=wx.OK)
488                    if dlg.ShowModal() == wx.ID_OK:
489                        SpGrp = dlg.GetValue()
490                        E,SGData = G2spc.SpcGroup(SpGrp)
491                    else:
492                        SGData = G2obj.P1SGData # P 1
493                        self.warnings += '\nThe space group was not interpreted and has been set to "P 1".'
494                        self.warnings += "Change this in phase's General tab."           
495                    dlg.Destroy()
496                G2spc.SGPrint(SGData) #silent check of space group symbol
497            elif 'qi' in S[:2]:
498                if nqi:
499                    raise self.ImportException("Supersymmetry too high; GSAS-II limited to (3+1) supersymmetry")           
500                vec = S.split()[1:]
501                SuperVec = [[float(vec[i]) for i in range(3)],False,4]
502                nqi += 1
503            elif 'atom' in S[:4]:
504                Atypes.append(S.split()[1])
505            S = fp.readline()
506            line += 1
507        fp.close()
508        #read atoms from m40 file
509        if not SGData:
510            self.warnings += '\nThe space group was not read before atoms and has been set to "P 1". '
511            self.warnings += "Change this in phase's General tab."
512            SGData = G2obj.P1SGData # P 1
513        waveTypes = ['Fourier','Sawtooth','ZigZag',]
514        filename2 = os.path.splitext(filename)[0]+'.m40'
515        file2 = open(filename2,'Ur')
516        S = file2.readline()
517        line = 1
518        self.errors = 'Error reading at line '+str(line)
519        nAtoms = int(S.split()[0])
520        for i in range(4):
521            S = file2.readline()           
522        for i in range(nAtoms):
523            S1 = file2.readline()
524            S1N = S1.split()[-3:]   # no. occ, no. pos waves, no. ADP waves
525            S1N = [int(i) for i in S1N]
526            S1T = list(S1[60:63])
527            waveType = waveTypes[int(S1T[1])]
528            Spos = []
529            Sadp = []
530            Sfrac = []
531            Smag = []
532            XYZ = [float(S1[27:36]),float(S1[36:45]),float(S1[45:54])]
533            SytSym,Mult = G2spc.SytSym(XYZ,SGData)[:2]
534            aType = Atypes[int(S1[9:11])-1]
535            Name = S1[:8].strip()
536            if S1[11:15].strip() == '1':
537                S2 = file2.readline()
538                Uiso = S2[:9]
539                if version == '2000':
540                    Uiso = R2pisq*float(Uiso)/4.      #Biso -> Uiso
541                Uij = [0,0,0,0,0,0]
542                IA = 'I'
543            elif S1[11:15].strip() == '2':
544                S2 = file2.readline()
545                IA = 'A'
546                Uiso = 0.
547                Uij = [float(S2[:9]),float(S2[9:18]),float(S2[18:27]),
548                    float(S2[27:36]),float(S2[36:45]),float(S2[45:54])] #Uij in Jana2006!
549                if version == '2000':
550                    Uij = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(Uij)/Mast) #these things are betaij in Jana2000! need to convert to Uij
551            for i in range(S1N[0]):
552                if not i:
553                    FS = file2.readline()
554                    Sfrac.append(FS[:9])    #'O' or 'delta' = 'length' for crenel
555                    if int(S1T[0]):  #"", "Legendre" or "Xharm" in 18:27 for "crenel"!
556                        waveType = 'Crenel/Fourier' #all waves 'Fourier' no other choice
557                Sfrac.append(file2.readline()[:18]) #if not crenel = Osin & Ocos
558                # else Osin & Ocos except last one is X40 = 'Center'
559            for i in range(S1N[1]): 
560                Spos.append(file2.readline()[:54])
561            for i in range(S1N[2]):
562                Sadp.append(file2.readline()[:54]+file2.readline())
563            if sum(S1N):    #if any waves: skip mystery line?
564                file2.readline()
565            for i,it in enumerate(Sfrac):
566                print (i,it)
567                if not i:
568                    if 'Crenel' in waveType:
569                        vals = [float(it),float(Sfrac[-1][:9])]
570                    else:
571                        vals = [float(it),]
572                else:
573                    vals = [float(it[:9]),float(it[9:18])]
574                if 'Crenel' in waveType and i == len(Sfrac)-1:
575                    del Sfrac[-1]
576                    break               
577                Sfrac[i] = [vals,False]
578                print (Sfrac[i])
579            for i,it in enumerate(Spos):
580                if waveType in ['Sawtooth',] and not i:
581                    vals = [float(it[:9]),float(it[9:18]),float(it[18:27]),float(it[27:36])]
582                else:
583                    vals = [float(it[:9]),float(it[9:18]),float(it[18:27]),float(it[27:36]),float(it[36:45]),float(it[45:54])]
584                Spos[i] = [vals,False]
585            for i,it in enumerate(Sadp):
586                vals = [float(it[:9]),float(it[9:18]),float(it[18:27]),float(it[27:36]),float(it[36:45]),float(it[45:54]),
587                    float(it[54:63]),float(it[63:72]),float(it[72:81]),float(it[81:90]),float(it[90:99]),float(it[99:108])]
588                #these are betaij modulations in Jana2000! need to convert to Uij modulations
589                if version == '2000':               
590                    vals[:6] = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(vals[:6])/Mast)    #convert sin bij to Uij
591                    vals[6:] = R2pisq*G2lat.UijtoU6(G2lat.U6toUij(vals[6:])/Mast)    #convert cos bij to Uij
592                Sadp[i] = [vals,False]
593            Atom = [Name,aType,'',XYZ[0],XYZ[1],XYZ[2],1.0,SytSym,Mult,IA,Uiso]
594            Atom += Uij
595            Atom.append(ran.randint(0,sys.maxsize))
596            Atom.append({'SS1':{'Sfrac':[waveType,]+Sfrac,'Spos':[waveType,]+Spos,'Sadp':['Fourier',]+Sadp,'Smag':['Fourier',]+Smag}})    #SS2 is for (3+2), etc.
597            Atoms.append(Atom)
598        file2.close()
599        self.errors = 'Error after read complete'
600        if not SGData:
601            raise self.ImportException("No space group (spcgroup entry) found")
602        if not cell:
603            raise self.ImportException("No cell found")
604        Phase = G2obj.SetNewPhase(Name=Title,SGData=SGData,cell=cell+[Volume,])
605        Phase['General']['Type'] = Type
606        Phase['General']['Modulated'] = True
607        Phase['General']['Super'] = nqi
608        Phase['General']['SuperVec'] = SuperVec
609        Phase['General']['SuperSg'] = SuperSg
610        if SuperSg:
611            Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
612        Phase['General']['AtomPtrs'] = [3,1,7,9]
613        Phase['Atoms'] = Atoms
614        return Phase
615   
616class PDF_ReaderClass(G2obj.ImportPhase):
617    'Routine to import Phase information from ICDD PDF Card files'
618    def __init__(self):
619        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
620            extensionlist=('.str',),
621            strictExtension=True,
622            formatName = 'ICDD .str',
623            longFormatName = 'ICDD PDF Card (.str file) import'
624            )
625       
626    def ContentsValidator(self, filename):
627        'Look for a str tag in 1st line' 
628        fp = open(filename,'r')
629        if fp.read(3) == 'str':
630            fp.close()
631            return True
632        self.errors = 'File does not begin with str tag'
633        fp.close()
634        return False
635
636    def Reader(self,filename, ParentFrame=None, **unused):
637        'Read phase from a ICDD .str file using :meth:`ReadPDFPhase`'
638        fp = open(filename,'r')
639        self.Phase = self.ReadPDFPhase(ParentFrame, fp)
640        fp.close()
641        return True
642
643    def ReadPDFPhase(self, G2frame,fp):
644        '''Read a phase from a ICDD .str file.
645        '''
646        EightPiSq = 8.*math.pi**2
647        self.errors = 'Error opening file'
648        Phase = {}
649        Atoms = []
650        S = fp.readline()
651        line = 1
652        SGData = None
653        cell = []
654        cellkey = []
655        while S:
656            if 'space_group' in S:
657                break
658            S = fp.readline()                   
659        while S:
660            self.errors = 'Error reading at line '+str(line)
661            if 'phase_name' in S:
662                Title = S.split('"')[1]
663            elif 'Space group (HMS)' in S:
664                SpGrp = S.split()[-1]
665                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
666                E,SGData = G2spc.SpcGroup(SpGrpNorm)
667                # space group processing failed, try to look up name in table
668                while E:
669                    print (G2spc.SGErrors(E))
670                    dlg = wx.TextEntryDialog(G2frame,
671                        SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
672                        'ERROR in space group symbol','',style=wx.OK)
673                    if dlg.ShowModal() == wx.ID_OK:
674                        SpGrp = dlg.GetValue()
675                        E,SGData = G2spc.SpcGroup(SpGrp)
676                    else:
677                        SGData = G2obj.P1SGData # P 1
678                        self.warnings += '\nThe space group was not interpreted and has been set to "P 1".'
679                        self.warnings += "Change this in phase's General tab."           
680                    dlg.Destroy()
681                G2spc.SGPrint(SGData) #silent check of space group symbol
682            elif 'a a_' in S[:7]:
683                data = S.split()
684                cell.append(float(data[2]))
685                cellkey.append(data[1])
686            elif 'b b_' in S[:7]:
687                data = S.split()
688                cell.append(float(data[2]))
689                cellkey.append(data[1])
690            elif 'b =' in S[:6]:
691                data = S.split('=')
692                indx = cellkey.index(data[1].split(';')[0])
693                cell.append(cell[indx])
694            elif 'c c_' in S[:7]:
695                data = S.split()
696                cell.append(float(data[2]))
697            elif 'c =' in S[:6]:
698                data = S.split('=')
699                indx = cellkey.index(data[1].split(';')[0])
700                cell.append(cell[indx])
701            elif 'al' in S[:5]:
702                cell.append(float(S.split()[1]))
703            elif 'be' in S[:5]:
704                cell.append(float(S.split()[1]))
705            elif 'ga' in S[:5]:
706                cell.append(float(S.split()[1]))
707                Volume = G2lat.calc_V(G2lat.cell2A(cell))
708                break
709            S = fp.readline()
710        S = fp.readline()
711        while S:
712            if '/*' in S[:5]:
713                break
714            if 'site' in S[:7]:
715                atom = []
716                xyzkey = []
717                data = S.split()
718                atom.append(data[1])    #name
719                pos = data.index('occ')+1
720                atom.append(data[pos])  #type
721                atom.append('')         #refine
722                for xid in ['x =','y =','z =']:
723                    if xid in S:
724                        xpos = S.index(xid)+3
725                        xend = xpos+S[xpos:].index(';')
726                        if S[xpos:xend] in xyzkey:
727                            atom.append(atom[3+xyzkey.index(S[xpos:xend])])
728                        else:
729                            atom.append(eval(S[xpos:xend]+'.'))
730                    else:
731                        xpos = data.index(xid[0])+2
732                        xyzkey.append(data[xpos-1][1:])
733                        atom.append(float(data[xpos]))
734                atom.append(float(data[pos+2]))
735                SytSym,Mult = G2spc.SytSym(np.array(atom[3:6]),SGData)[:2]
736                atom.append(SytSym)
737                atom.append(Mult)
738                if 'beq' in S:
739                    atom.append('I')
740                    upos = data.index('beq')
741                    atom.append(float(data[upos+2])/EightPiSq)
742                    atom += [0.,0.,0.,0.,0.,0.,]
743                elif 'ADPs' in S:
744                    upos = data.index('ADPs')
745                    atom.append('A')
746                    atom.append(0.0)
747                    for uid in ['Bani11','Bani22','Bani33','Bani12','Bani13','Bani23']:
748                        upos = data.index(uid)+1
749                        atom.append(float(data[upos])/EightPiSq)
750                else:
751                    atom.append('I')
752                    atom += [0.02,0.,0.,0.,0.,0.,0.,]                   
753                atom.append(ran.randint(0,sys.maxsize))
754                Atoms.append(atom)
755            S = fp.readline()               
756        fp.close()
757        self.errors = 'Error after read complete'
758        if not SGData:
759            raise self.ImportException("No space group (spcgroup entry) found")
760        if not cell:
761            raise self.ImportException("No cell found")
762        Phase = G2obj.SetNewPhase(Name=Title,SGData=SGData,cell=cell+[Volume,])
763        Phase['General']['Type'] = 'nuclear'
764        Phase['General']['AtomPtrs'] = [3,1,7,9]
765        Phase['Atoms'] = Atoms
766        return Phase
Note: See TracBrowser for help on using the repository browser.