Changeset 1090 for trunk/exports


Ignore:
Timestamp:
Oct 8, 2013 2:33:36 PM (8 years ago)
Author:
vondreele
Message:

complete 1st version of export PDB

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/exports/G2export_PDB.py

    r1089 r1090  
    1010'''Code to demonstrate how export routines are created
    1111'''
     12import numpy as np
    1213import os.path
    1314import GSASIIpath
     
    4041        '''Export as a PDB file
    4142        '''
     43       
     44        def PDBheader():
     45            self.Write("HEADER phase "+str(phasenam)+" from "+str(self.G2frame.GSASprojectfile))
     46            self.Write("TITLE")
     47            self.Write("COMPND")
     48            self.Write("SOURCE")
     49            self.Write("KEYWDS")
     50            self.Write("EXPDTA    X-RAY POWDER DIFFRACTION")
     51            self.Write("REVDAT")
     52            self.Write("JRNL")
     53            self.Write("REMARK   1")
     54            self.Write("REMARK   2")                                                                     
     55            self.Write("REMARK   2 RESOLUTION. 2.66 ANGSTROMS.")                                         
     56            self.Write("REMARK   2 POWDER DIFFRACTION MINIMUM D-SPACING.")
     57           
     58        def PDBremark250():                               
     59            self.Write('REMARK 250')                                                                     
     60            self.Write('REMARK 250 REFINEMENT.')                                                         
     61            self.Write('REMARK 250   PROGRAM     : GSAS-II')                                                 
     62            self.Write('REMARK 250   AUTHORS     : TOBY & VON DREELE')
     63            self.Write('REMARK 250   REFRENCE    : J. APPL. CRYST. 46, 544-549(2013)')                             
     64            self.Write('REMARK 250')
     65            self.Write('REMARK 250  DATA USED IN REFINEMENT')                                             
     66            self.Write('REMARK 250   RESOLUTION RANGE HIGH (ANGSTROMS) :  x.xx')                         
     67            self.Write('REMARK 250   RESOLUTION RANGE LOW  (ANGSTROMS) : xx.xx')                         
     68            self.Write('REMARK 250   POWDER DIFFRACTION DATA.')                                           
     69            self.Write('REMARK 250')
     70            self.Write('REMARK 250  FIT TO DATA USED IN REFINEMENT')                                     
     71            self.Write('REMARK 250   NUMBER OF POWDER PATTERNS         :     x')                         
     72            self.Write('REMARK 250   PROFILE R VALUES              (%) :  x.xx')                         
     73            self.Write('REMARK 250   WEIGHTED PROFILE R VALUES     (%) :  x.xx')                         
     74            self.Write('REMARK 250   F**2 R VALUES                 (%) : xx.xx')                         
     75            self.Write('REMARK 250   NUMBERS OF POWDER PATTERN POINTS  :  xxxx')                         
     76            self.Write('REMARK 250   NUMBERS OF REFLECTIONS            :  xxxx')                         
     77            self.Write('REMARK 250   TOTAL NUMBER OF POWDER POINTS     :  xxxx')                         
     78            self.Write('REMARK 250')
     79            self.Write('REMARK 250  NUMBER OF NON-HYDROGEN ATOMS USED IN REFINEMENT.')                   
     80            self.Write('REMARK 250   PROTEIN ATOMS       :      xxxx')                                             
     81            self.Write('REMARK 250   NUCLEIC ACID ATOMS  :      xxxx')                                             
     82            self.Write('REMARK 250   HETEROGEN ATOMS     :      xxxx')                                             
     83            self.Write('REMARK 250   SOLVENT ATOMS       :      xxxx')                                             
     84            self.Write('REMARK 250')
     85            self.Write('REMARK 250  MODEL REFINEMENT.')                                                   
     86            self.Write('REMARK 250   NUMBER OF LEAST-SQUARES PARAMETERS :  xxxx')                         
     87            self.Write('REMARK 250   NUMBER OF RESTRAINTS               :  xxxx')                         
     88            self.Write('REMARK 250')
     89            self.Write('REMARK 250  RMS DEVIATIONS FROM RESTRAINT TARGET VALUES. NUMBER.')               
     90            self.Write('REMARK 250   BOND ANGLES                      (DEG) : x.xx   xxx')               
     91#            self.Write('REMARK 250   ANTI-BUMPING DISTANCE RESTRAINTS   (A) :x.xxx   xxx')               
     92#            self.Write('REMARK 250   HYDROGEN BOND DISTANCE RESTRAINTS  (A) :x.xxx   xxx')               
     93            self.Write('REMARK 250   INTERATOMIC DISTANCES              (A) :x.xxx   xxx')               
     94            self.Write('REMARK 250   DISTANCES FROM RESTRAINT PLANES    (A) :x.xxx   xxx')               
     95            self.Write('REMARK 250   TORSION PSEUDOPOTENTIAL RESTRAINTS (E) : x.xx   xxx')               
     96            self.Write('REMARK 250   TORSION ANGLE RESTRAINTS           (E) : x.xx   xxx')               
     97            self.Write('REMARK 250')
     98            self.Write('REMARK 200')                                                                     
     99            self.Write('DBREF')
     100           
     101        def PDBseqres(seqList):
     102            chains = seqList.keys()
     103            chains.sort()
     104            nSeq = 0
     105            for chain in chains:
     106                nres = len(seqList[chain])
     107                nrec = (nres-1)/13+1
     108                iB = 0
     109                for irec in range(nrec):
     110                    iF = min(iB+13,nres)
     111                    text = 'SEQRES {:3d}{:2s}{:5d}  '+(iF-iB)*'{:4s}'
     112                    self.Write(text.format(irec+1,chain,nres,*seqList[chain][iB:iF]))
     113                    nSeq += 1
     114                    iB += 13
     115            return nSeq
     116           
    42117        # the export process starts here
    43118        # load all of the tree into a set of dicts
     
    54129                print 'not macromolecular phase'
    55130                return
    56             print phasedict.keys()
    57             print General['SGData'].keys()           
    58131            i = self.Phases[phasenam]['pId']
    59132            if len(self.phasenam) > 1: # if more than one filename is included, add a phase #
     
    63136                fil = self.filename
    64137            fp = self.OpenFile(fil)
    65             self.Write("HEADER phase "+str(phasenam)+" from "+str(self.G2frame.GSASprojectfile))
    66             # get cell parameters
    67             pfx = str(phasedict['pId'])+'::'
    68             A,sigA = G2stIO.cellFill(pfx,phasedict['General']['SGData'],self.parmDict,self.sigDict)
    69             line = "CRYST1 {:8.3f} {:8.3f} {:8.3f} {:6.2f} {:6.2f} {:6.2f} ".format(*G2lat.A2cell(A))
     138            Atoms = phasedict['Atoms']
     139            cx,ct,cs,cia = General['AtomPtrs']
     140            seqList = {}
     141            AA3letter = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
     142                'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
     143            seq = 0
     144            for atom in Atoms:
     145                if atom[ct-3] in AA3letter and int(atom[ct-4]) != seq:
     146                    if atom[ct-2] not in seqList:
     147                        seqList[atom[ct-2]] = []
     148                    seqList[atom[ct-2]].append(atom[ct-3])
     149                    seq = int(atom[ct-4])
     150            PDBheader()
     151            PDBremark250()
     152            nSeq = PDBseqres(seqList)
     153           
     154            # get cell parameters
     155            Cell = General['Cell'][1:7]
     156            line = "CRYST1 {:8.3f} {:8.3f} {:8.3f} {:6.2f} {:6.2f} {:6.2f} ".format(*Cell)
    70157            line += General['SGData']['SpGrp'].ljust(13)
    71158            line += '%2d'%(len(General['SGData']['SGOps'])*len(General['SGData']['SGCen']))
     
    74161            self.Write('ORIGX2      0.000000  1.000000  0.000000        0.00000')
    75162            self.Write('ORIGX3      0.000000  0.000000  1.000000        0.00000')
    76             G = G2lat.A2Gmat(A)[0]
    77             A,B = G2lat.Gmat2AB(G)
     163            A,B = G2lat.cell2AB(Cell)
    78164            self.Write('SCALE1     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[0]))
    79165            self.Write('SCALE2     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[1]))
    80166            self.Write('SCALE3     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[2]))
    81             Atoms = phasedict['Atoms']
    82            
    83 
     167            iatom = 1
     168            nHet = 0
     169            nTer = 0
     170            for atom in Atoms:
     171                if atom[cia] == 'I':
     172                    Biso = atom[cia+1]*8.*np.pi**2
     173                xyz = np.inner(A,np.array(atom[cx:cx+3]))
     174                if atom[ct-3] in AA3letter:
     175                    fmt = 'ATOM{:7d}  {:4s}{:3s}{:2s}{:4s}    '+3*'{:8.3f}'+2*'{:6.2f}'+'{:s}'
     176                    self.Write(fmt.format(iatom,atom[ct-1],atom[ct-3],    \
     177                        atom[ct-2],atom[ct-4],xyz[0],xyz[1],xyz[2],atom[cx+3],Biso,atom[ct].rjust(12)))
     178                    if atom[ct-1] == 'OXT':
     179                        iatom += 1
     180                        self.Write('TER {:7d}  {:4s}{:3s}'.format(iatom,atom[ct-1],atom[ct-3]))
     181                        nTer += 1
     182                else:
     183                    nHet += 1
     184                    fmt = 'HETATM{:5d} {:5s}{:3s}{:2s}{:4s}    '+3*'{:8.3f}'+2*'{:6.2f}'+'{:s}'
     185                    self.Write(fmt.format(iatom,atom[ct-1].ljust(5),atom[ct-3],    \
     186                        atom[ct-2],atom[ct-4],xyz[0],xyz[1],xyz[2],atom[cx+3],Biso,atom[ct].rjust(12)))
     187                iatom += 1
     188           
     189            vals = [3,0,nHet,0,0,0,6,len(Atoms),nTer,0,nSeq]
     190            fmt = 'MASTER'+11*'{:5d}'
     191            self.Write(fmt.format(*vals))
     192            self.Write('END')
    84193            self.CloseFile()
    85194            print('Phase '+str(phasenam)+' written to file '+str(fil))
     195
     196   
Note: See TracChangeset for help on using the changeset viewer.