source: trunk/exports/G2export_PDB.py @ 1090

Last change on this file since 1090 was 1090, checked in by vondreele, 9 years ago

complete 1st version of export PDB

File size: 10.5 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: $
5# $Author: $
6# $Revision: -1 $
7# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/exports/G2export_CIF.py $
8# $Id: G2export_CIF.py -1   $
9########### SVN repository information ###################
10'''Code to demonstrate how export routines are created
11'''
12import numpy as np
13import os.path
14import GSASIIpath
15GSASIIpath.SetVersionNumber("$Revision: -1 $")
16import GSASIIIO as G2IO
17#import GSASIIgrid as G2gd
18import GSASIIstrIO as G2stIO
19#import GSASIImath as G2mth
20import GSASIIlattice as G2lat
21import GSASIIspc as G2spc
22#import GSASIIphsGUI as G2pg
23#import GSASIIstrMain as G2stMn
24
25class ExportPhasePDB(G2IO.ExportBaseclass):
26    '''Used to create a PDB file for a phase
27
28    :param wx.Frame G2frame: reference to main GSAS-II frame
29    '''
30    def __init__(self,G2frame):
31        super(self.__class__,self).__init__( # fancy way to say <parentclass>.__init__
32            G2frame=G2frame,
33            formatName = 'PDB',
34            extension='.PDB',
35            longFormatName = 'Export phase as .PDB file'
36            )
37        self.exporttype = ['phase']
38        self.multiple = True
39
40    def Exporter(self,event=None):
41        '''Export as a PDB file
42        '''
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           
117        # the export process starts here
118        # load all of the tree into a set of dicts
119        self.loadTree()
120        # create a dict with refined values and their uncertainties
121        self.loadParmDict()
122        if self.SetupExport(event,                         # set export parameters
123                            AskFile=True
124                            ): return 
125        for phasenam in self.phasenam:
126            phasedict = self.Phases[phasenam] # pointer to current phase info
127            General = phasedict['General']
128            if General['Type'] != 'macromolecular':
129                print 'not macromolecular phase'
130                return
131            i = self.Phases[phasenam]['pId']
132            if len(self.phasenam) > 1: # if more than one filename is included, add a phase #
133                nam,ext = os.path.splitext(self.filename)
134                fil = nam+"_"+str(i)+ext
135            else:
136                fil = self.filename
137            fp = self.OpenFile(fil)
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)
157            line += General['SGData']['SpGrp'].ljust(13)
158            line += '%2d'%(len(General['SGData']['SGOps'])*len(General['SGData']['SGCen']))
159            self.Write(line)
160            self.Write('ORIGX1      1.000000  0.000000  0.000000        0.00000')
161            self.Write('ORIGX2      0.000000  1.000000  0.000000        0.00000')
162            self.Write('ORIGX3      0.000000  0.000000  1.000000        0.00000')
163            A,B = G2lat.cell2AB(Cell)
164            self.Write('SCALE1     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[0]))
165            self.Write('SCALE2     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[1]))
166            self.Write('SCALE3     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[2]))
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')
193            self.CloseFile()
194            print('Phase '+str(phasenam)+' written to file '+str(fil))
195
196   
Note: See TracBrowser for help on using the repository browser.