source: trunk/imports/G2phase.py @ 4415

Last change on this file since 4415 was 4415, checked in by vondreele, 3 years ago

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