source: trunk/exports/G2export_PDB.py @ 3245

Last change on this file since 3245 was 3245, checked in by vondreele, 7 years ago

fix problem with incommensurate magnetic structures from CIFhklReader
make "Edit Body" as "Edit Residue Body" in RB menu
change protein validation bar plots to be vertical
change 'twin' to 'flag' in Reflection List table & add '-2' as 'free' (for proteins - not yet implemented)
fix PDB phase import to get atom type from 76:78 of ATOM record
change importer names for single crystal TOF data to be more explicit (SNS vs ISIS)
change cif structure factor importer - F2, sig(F2) & Fcalc now allowed
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 12.6 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2018-01-27 14:49:12 +0000 (Sat, 27 Jan 2018) $
5# $Author: vondreele $
6# $Revision: 3245 $
7# $URL: trunk/exports/G2export_PDB.py $
8# $Id: G2export_PDB.py 3245 2018-01-27 14:49:12Z vondreele $
9########### SVN repository information ###################
10'''
11*Module G2export_PDB: Macromolecular export*
12--------------------------------------------
13Code to export a phase into the venerated/obsolete (pick one)
14ASCII PDB format. Also defines exporter :class:`ExportPhaseCartXYZ`
15which writes atom positions in orthogonal coordinates for a phase.
16
17'''
18from __future__ import division, print_function
19import numpy as np
20import os.path
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 3245 $")
23import GSASIIIO as G2IO
24import GSASIIlattice as G2lat
25
26class ExportPhasePDB(G2IO.ExportBaseclass):
27    '''Used to create a PDB file for a phase
28
29    :param wx.Frame G2frame: reference to main GSAS-II frame
30    '''
31    def __init__(self,G2frame):
32        super(self.__class__,self).__init__( # fancy way to say <parentclass>.__init__
33            G2frame=G2frame,
34            formatName = 'PDB',
35            extension='.PDB',
36            longFormatName = 'Export phase as .PDB file'
37            )
38        self.exporttype = ['phase']
39        self.multiple = True
40
41    def Exporter(self,event=None):
42        '''Export as a PDB file
43        '''
44       
45        def PDBheader():
46            self.Write("HEADER phase "+str(phasenam)+" from "+str(self.G2frame.GSASprojectfile))
47            self.Write("TITLE")
48            self.Write("COMPND")
49            self.Write("SOURCE")
50            self.Write("KEYWDS")
51            self.Write("EXPDTA    X-RAY POWDER DIFFRACTION")
52            self.Write("REVDAT")
53            self.Write("JRNL")
54            self.Write("REMARK   1")
55            self.Write("REMARK   2")                                                                     
56            self.Write("REMARK   2 RESOLUTION. 2.66 ANGSTROMS.")                                         
57            self.Write("REMARK   2 POWDER DIFFRACTION MINIMUM D-SPACING.")
58           
59        def PDBremark250():                               
60            self.Write('REMARK 250')                                                                     
61            self.Write('REMARK 250 REFINEMENT.')                                                         
62            self.Write('REMARK 250   PROGRAM     : GSAS-II')                                                 
63            self.Write('REMARK 250   AUTHORS     : TOBY & VON DREELE')
64            self.Write('REMARK 250   REFRENCE    : J. APPL. CRYST. 46, 544-549(2013)')                             
65            self.Write('REMARK 250')
66            self.Write('REMARK 250  DATA USED IN REFINEMENT')                                             
67            self.Write('REMARK 250   RESOLUTION RANGE HIGH (ANGSTROMS) :  x.xx')                         
68            self.Write('REMARK 250   RESOLUTION RANGE LOW  (ANGSTROMS) : xx.xx')                         
69            self.Write('REMARK 250   POWDER DIFFRACTION DATA.')                                           
70            self.Write('REMARK 250')
71            self.Write('REMARK 250  FIT TO DATA USED IN REFINEMENT')                                     
72            self.Write('REMARK 250   NUMBER OF POWDER PATTERNS         :     x')                         
73            self.Write('REMARK 250   PROFILE R VALUES              (%) :  x.xx')                         
74            self.Write('REMARK 250   WEIGHTED PROFILE R VALUES     (%) :  x.xx')                         
75            self.Write('REMARK 250   F**2 R VALUES                 (%) : xx.xx')                         
76            self.Write('REMARK 250   NUMBERS OF POWDER PATTERN POINTS  :  xxxx')                         
77            self.Write('REMARK 250   NUMBERS OF REFLECTIONS            :  xxxx')                         
78            self.Write('REMARK 250   TOTAL NUMBER OF POWDER POINTS     :  xxxx')                         
79            self.Write('REMARK 250')
80            self.Write('REMARK 250  NUMBER OF NON-HYDROGEN ATOMS USED IN REFINEMENT.')                   
81            self.Write('REMARK 250   PROTEIN ATOMS       :      xxxx')                                             
82            self.Write('REMARK 250   NUCLEIC ACID ATOMS  :      xxxx')                                             
83            self.Write('REMARK 250   HETEROGEN ATOMS     :      xxxx')                                             
84            self.Write('REMARK 250   SOLVENT ATOMS       :      xxxx')                                             
85            self.Write('REMARK 250')
86            self.Write('REMARK 250  MODEL REFINEMENT.')                                                   
87            self.Write('REMARK 250   NUMBER OF LEAST-SQUARES PARAMETERS :  xxxx')                         
88            self.Write('REMARK 250   NUMBER OF RESTRAINTS               :  xxxx')                         
89            self.Write('REMARK 250')
90            self.Write('REMARK 250  RMS DEVIATIONS FROM RESTRAINT TARGET VALUES. NUMBER.')               
91            self.Write('REMARK 250   BOND ANGLES                      (DEG) : x.xx   xxx')               
92#            self.Write('REMARK 250   ANTI-BUMPING DISTANCE RESTRAINTS   (A) :x.xxx   xxx')               
93#            self.Write('REMARK 250   HYDROGEN BOND DISTANCE RESTRAINTS  (A) :x.xxx   xxx')               
94            self.Write('REMARK 250   INTERATOMIC DISTANCES              (A) :x.xxx   xxx')               
95            self.Write('REMARK 250   DISTANCES FROM RESTRAINT PLANES    (A) :x.xxx   xxx')               
96            self.Write('REMARK 250   TORSION PSEUDOPOTENTIAL RESTRAINTS (E) : x.xx   xxx')               
97            self.Write('REMARK 250   TORSION ANGLE RESTRAINTS           (E) : x.xx   xxx')               
98            self.Write('REMARK 250')
99            self.Write('REMARK 200')                                                                     
100            self.Write('DBREF')
101           
102        def PDBseqres(seqList):
103            chains = list(seqList.keys())
104            chains.sort()
105            nSeq = 0
106            for chain in chains:
107                nres = len(seqList[chain])
108                nrec = (nres-1)//13+1
109                iB = 0
110                for irec in range(nrec):
111                    iF = min(iB+13,nres)
112                    text = 'SEQRES {:3d}{:2s}{:5d}  '+(iF-iB)*'{:4s}'
113                    self.Write(text.format(irec+1,chain,nres,*seqList[chain][iB:iF]))
114                    nSeq += 1
115                    iB += 13
116            return nSeq
117           
118        # the export process starts here
119        self.InitExport(event)
120        # load all of the tree into a set of dicts
121        self.loadTree()
122        # create a dict with refined values and their uncertainties
123        self.loadParmDict()
124        if self.ExportSelect():    # set export parameters; ask for file name
125            return
126        filename = self.filename
127        for phasenam in self.phasenam:
128            phasedict = self.Phases[phasenam] # pointer to current phase info
129            General = phasedict['General']
130            if General['Type'] != 'macromolecular':
131                print ('phase '+phasenam+' not macromolecular, skipping')
132                continue
133            i = self.Phases[phasenam]['pId']
134            if len(self.phasenam) > 1: # if more than one filename is included, add a phase #
135                self.filename = os.path.splitext(filename)[1] + "_" + str(i) + self.extension
136            fp = self.OpenFile()
137            Atoms = phasedict['Atoms']
138            cx,ct,cs,cia = General['AtomPtrs']
139            seqList = {}
140            AA3letter = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
141                'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
142            seq = 0
143            for atom in Atoms:
144                if atom[ct-3] in AA3letter and int(atom[ct-4]) != seq:
145                    if atom[ct-2] not in seqList:
146                        seqList[atom[ct-2]] = []
147                    seqList[atom[ct-2]].append(atom[ct-3])
148                    seq = int(atom[ct-4])
149            PDBheader()
150            PDBremark250()
151            nSeq = PDBseqres(seqList)
152           
153            # get cell parameters
154            Cell = General['Cell'][1:7]
155            line = "CRYST1 {:8.3f} {:8.3f} {:8.3f} {:6.2f} {:6.2f} {:6.2f} ".format(*Cell)
156            line += General['SGData']['SpGrp'].ljust(13)
157            line += '%2d'%(len(General['SGData']['SGOps'])*len(General['SGData']['SGCen']))
158            self.Write(line)
159            self.Write('ORIGX1      1.000000  0.000000  0.000000        0.00000')
160            self.Write('ORIGX2      0.000000  1.000000  0.000000        0.00000')
161            self.Write('ORIGX3      0.000000  0.000000  1.000000        0.00000')
162            A,B = G2lat.cell2AB(Cell)
163            self.Write('SCALE1     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[0]))
164            self.Write('SCALE2     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[1]))
165            self.Write('SCALE3     {:9.6f} {:9.6f} {:9.6f}        0.00000'.format(*B[2]))
166            iatom = 1
167            nHet = 0
168            nTer = 0
169            fmt = '{:6s}{:5d}  {:4s}{:3s} {:1s}{:4s}    '+3*'{:8.3f}'+2*'{:6.2f}'+'{:s}'
170            for atom in Atoms:
171                if atom[cia] == 'I':    #need to deal with aniso thermals for proteins = "ANISOU" records
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                    self.Write(fmt.format('ATOM  ',iatom,atom[ct-1],atom[ct-3].strip(),    \
176                        atom[ct-2].strip(),atom[ct-4].rjust(4),xyz[0],xyz[1],xyz[2],atom[cx+3], \
177                        Biso,atom[ct].rjust(12)))
178                    if atom[ct-1] == 'OXT':
179                        iatom += 1
180                        self.Write('{:6s}{:5d}  {:4s}{:3s}'.format('TER   ',iatom,atom[ct-1],atom[ct-3].strip()))
181                        nTer += 1
182                else:
183                    nHet += 1
184                    self.Write(fmt.format('HETATM',iatom,atom[ct-1],atom[ct-3].strip(),    \
185                        atom[ct-2].strip(),atom[ct-4].rjust(4),xyz[0],xyz[1],xyz[2],atom[cx+3], \
186                        Biso,atom[ct].rjust(12)))
187                #if atim[cia] == 'a':
188                #   put in 'ANISOU' record
189                #'ANISOU    1  N   ALA A 340     4392   4159   4615    249   -189     73       N' 
190                iatom += 1
191           
192            vals = [3,0,nHet,0,0,0,6,len(Atoms),nTer,0,nSeq]
193            fmt = 'MASTER'+11*'{:5d}'
194            self.Write(fmt.format(*vals))
195            self.Write('END')
196            self.CloseFile()
197            print('Phase '+phasenam+' written to PDB file '+self.fullpath)
198
199class ExportPhaseCartXYZ(G2IO.ExportBaseclass):
200    '''Used to create a Cartesian XYZ file for a phase
201
202    :param wx.Frame G2frame: reference to main GSAS-II frame
203    '''
204    def __init__(self,G2frame):
205        super(self.__class__,self).__init__( # fancy way to say <parentclass>.__init__
206            G2frame=G2frame,
207            formatName = 'Cartesian XYZ',
208            extension='.XYZ',
209            longFormatName = 'Export phase with Cartesian coordinates as .XYZ file'
210            )
211        self.exporttype = ['phase']
212        self.multiple = True
213
214    def Exporter(self,event=None):
215        '''Export as a XYZ file
216        '''
217        # the export process starts here
218        self.InitExport(event)
219        # load all of the tree into a set of dicts
220        self.loadTree()
221        # create a dict with refined values and their uncertainties
222        self.loadParmDict()
223        if self.ExportSelect():    # set export parameters; ask for file name
224            return
225        filename = self.filename
226        for phasenam in self.phasenam:
227            phasedict = self.Phases[phasenam] # pointer to current phase info
228            General = phasedict['General']
229            i = self.Phases[phasenam]['pId']
230            Atoms = phasedict['Atoms']
231            if not len(Atoms):
232                print('**** ERROR - Phase '+phasenam+' has no atoms! ****')
233                continue
234            if len(self.phasenam) > 1: # if more than one filename is included, add a phase #
235                self.filename = os.path.splitext(filename)[1] + "_" + str(i) + self.extension
236            cx,ct,cs,cia = General['AtomPtrs']
237            Cell = General['Cell'][1:7]
238            A,B = G2lat.cell2AB(Cell)
239            fmt = '{:4s}'+3*'{:12.4f}'
240            self.Write('{:6d}'.format(len(Atoms)))
241            self.Write(' ')
242            for atom in Atoms:
243                xyz = np.inner(A,np.array(atom[cx:cx+3]))
244                self.Write(fmt.format(atom[ct],*xyz))
245            self.CloseFile()
246            print('Phase '+phasenam+' written to XYZ file '+self.fullpath)
247   
Note: See TracBrowser for help on using the repository browser.