'''
    SASSIE: Copyright (C) 2011 Joseph E. Curtis, Ph.D. 

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
'''
import os
import sys
import string
import locale
import struct
import numpy
import time
import dcdio

#	SASIO
#
#	12/5/2009	--	initial coding			:	jc
#	12/10/2009	--	doc strings 			:	jc
#	01/01/2011	--	added dcdio wrappers		:	jc
#
#LC	 1         2         3         4         5         6         7
#LC4567890123456789012345678901234567890123456789012345678901234567890123456789
#								       *      **
'''
	Sasio is the main module that contains the base classes that 
	read and write atomic information from and to the hard disk,
	and (eventually) deal with logging of input parameters and
	project runs.

	The methods in class Files are used to read and write data to
	the Charmm/Xplor binary data format (DCD) and textual protein
	data bank (PDB) format.  

	See the following sites for the DCD format:

	http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html
	http://www.lrz-muenchen.de/~heller/ego/manual/node93.html

	The methods (read_pdb, write_pdb) are used to read and write data
	using the PDB file format as described at the following
	web-site:

	http://deposit.rcsb.org/adit/docs/pdb_atom_format.html

	These classes are accessed by the SasAtm class found in
	the sasmol module.

'''

class Files(object):

	def __init__(self,filename,flag):
		pass

	def open_dcd_read(self,filename):
		'''
		This method opens a file to read in the Charmm/Xplor data format.
		by calling a pre-compiled C module (dcdio).
		'''

		filepointer=dcdio.open_dcd_read(filename)

		num_fixed=0
		result = 1

		readheaderresult,nnatoms,nset,istart,nsavc,delta,namnf,reverseEndian,charmm=dcdio.read_dcdheader(filepointer)
		if(readheaderresult!=0):
			print 'failed to read header'
			print 'readheaderresult = ',readheaderresult	
	
		print 'result = ',readheaderresult
	
		dcdfile = [filepointer,nnatoms,nset,reverseEndian,charmm]
	
		return dcdfile

	def open_dcd_write(self,filename):
		'''
		This method opens a file and writes the headerfile in the Charmm/Xplor data format.
		by calling a pre-compiled C module (dcdio).

		This function will OVERWRITE a file with the same name without prompting.
		'''

		filepointer=dcdio.open_dcd_write(filename)

		self.write_dcd_header(filepointer,1)

		return filepointer

	def write_dcd_header(self,filepointer,nset):
		'''
		This method writes the headerfile in the Charmm/Xplor data format.
		by calling a pre-compiled C module (dcdio).
		'''
		filename=" "
		natoms = self._coor[0,:,0].shape[0]
		istart = 0 ; nsavc = 1 ; delta = 1.0
		headerresult=dcdio.write_dcdheader(filepointer,filename,natoms,nset,istart,nsavc,delta)

		return 

	def write_dcd_step(self,filepointer,frame,step):
		'''
		This method writes a single step in the Charmm/Xplor data format.
		by calling a pre-compiled C module (dcdio).
		'''

		tx=self._coor[frame,:,0].astype(numpy.float32)	
		ty=self._coor[frame,:,1].astype(numpy.float32)	
		tz=self._coor[frame,:,2].astype(numpy.float32)	

		stepresult=dcdio.write_dcdstep(filepointer,tx,ty,tz,step)

		return

	def close_dcd_write(self,filepointer):
		'''
		This method closes a dcd file.
		'''
		result = dcdio.close_dcd_write(filepointer)
		print 'result = ',result
	
		return
	
	def close_dcd_read(self,filepointer):
		'''
		This method closes a dcd file.
		'''
		result = dcdio.close_dcd_read(filepointer)
		print 'result = ',result
	
		return

	def write_dcd(self,filename):
		'''
		This method writes data in the Charmm/Xplor data format.
		by calling a pre-compiled C module (dcdio).

		'''
		
		outfile=dcdio.open_dcd_write(filename)

		nset = self._coor[:,0,0].shape[0]
		natoms = self._coor[0,:,0].shape[0]
		istart = 0 ; nsavc = 1 ; delta = 1.0
		
		print 'writing nset = ',nset,' to headerfile'		
		print 'writing nset = ',nset,' to headerfile'		
		print 'writing nset = ',nset,' to headerfile'		
		print 'writing nset = ',nset,' to headerfile'		

		headerresult=dcdio.write_dcdheader(outfile,filename,natoms,nset,istart,nsavc,delta)

		for frame in xrange(nset):
			print ".",
			sys.stdout.flush()

			tx=self._coor[frame,:,0].astype(numpy.float32)	
			ty=self._coor[frame,:,1].astype(numpy.float32)	
			tz=self._coor[frame,:,2].astype(numpy.float32)	

			stepresult=dcdio.write_dcdstep(outfile,tx,ty,tz,frame+1)

		result = dcdio.close_dcd_write(outfile)

	def read_single_dcd_step(self,filename,frame):
		'''
		This method reads a single dcd step in the Charmm/Xplor data format.
		'''

		infile=dcdio.open_dcd_read(filename)
		num_fixed=0
		result = 1

		readheaderresult,nnatoms,nset,istart,nsavc,delta,namnf,reverseEndian,charmm=dcdio.read_dcdheader(infile)
		if(readheaderresult!=0):
			print 'failed to read header'
			print 'readheaderresult = ',readheaderresult	

		coor=numpy.zeros((1,nnatoms,3),numpy.float)	
	
		tx=numpy.zeros(nnatoms,dtype=numpy.float32)
		ty=numpy.zeros(nnatoms,dtype=numpy.float32)
		tz=numpy.zeros(nnatoms,dtype=numpy.float32)

		for i in xrange(frame+1):	
			result=dcdio.read_dcdstep(infile,tx,ty,tz,num_fixed,i,reverseEndian,charmm)

		coor[0,:,0]=tx.astype(numpy.float) ; coor[0,:,1]=ty.astype(numpy.float) ; coor[0,:,2]=tz.astype(numpy.float)
		
		result = dcdio.close_dcd_read(infile)
		self._coor=numpy.array(coor)
	
		if(result!=0):
			print 'failed to read coordinates'	
			print 'result = ',result

		return
	
	def read_dcd_step(self,dcdfile,frame):
		'''
		This method reads a single dcd step in the Charmm/Xplor data format.
		'''
		num_fixed=0

		filepointer = dcdfile[0]
		nnatoms = dcdfile[1]
		reverseEndian = dcdfile[3]
		charmm = dcdfile[4]

		tx=numpy.zeros(nnatoms,dtype=numpy.float32)
		ty=numpy.zeros(nnatoms,dtype=numpy.float32)
		tz=numpy.zeros(nnatoms,dtype=numpy.float32)

		result=dcdio.read_dcdstep(filepointer,tx,ty,tz,num_fixed,frame,reverseEndian,charmm)

		self._coor[0,:,0]=tx.astype(numpy.float) ; self._coor[0,:,1]=ty.astype(numpy.float) ; self._coor[0,:,2]=tz.astype(numpy.float)
	
		print '.',

		return

	def read_dcd(self,filename):

		'''
		This method reads data in the Charmm/Xplor data format.
		'''

		infile=dcdio.open_dcd_read(filename)

		nnatoms=0 ; nset=0 ; istart=0 ; nsavc=0 ; delta=0.0
		namnf=0 ; freeindexes=[] ; reverseEndian=0 ; charmm=0

		readheaderresult,nnatoms,nset,istart,nsavc,delta,namnf,reverseEndian,charmm=dcdio.read_dcdheader(infile)
		coor=numpy.zeros((nset,nnatoms,3),numpy.float)	
	
		num_fixed=0 
		result=1

		sum=0.0
		for i in xrange(nset):
			print '.',
 			sys.stdout.flush()
 			read_start_time=time.time()
			tx=numpy.zeros(nnatoms,dtype=numpy.float32)
			ty=numpy.zeros(nnatoms,dtype=numpy.float32)
			tz=numpy.zeros(nnatoms,dtype=numpy.float32)
		
			result=dcdio.read_dcdstep(infile,tx,ty,tz,num_fixed,i,reverseEndian,charmm)
 			read_end_time=time.time()
	
			sum+=read_end_time-read_start_time

			coor[i,:,0]=tx.astype(numpy.float) ; coor[i,:,1]=ty.astype(numpy.float) ; coor[i,:,2]=tz.astype(numpy.float)
	
		result = dcdio.close_dcd_read(infile)
		self._coor=numpy.array(coor)

	def element_filter(self):
		'''	
		This function filters the PDB file to sraighten
		out various things (add to this as you go along).
		'''

		'''
		Filter element list... 
		'''
		unique_elements = []
		error = []
		single_atom_names = ['H','F','B','D','C','N','O','S','P','I','K','U','V','W','Y']
		for i in range(len(self._name)):
			name = self._name[i].upper()
			resname = self._resname[i].upper()
			if (self._element[i]=='' or self._element[i]==' ' or self._element[i]=='  '):
				natoms = len(self._name[i])
				if(natoms == 1):
					if(name in single_atom_names):
						self._element[i]=self._name[i][0]
					else:
						message = 'atom name'+name+' not found in periodic table or charmm list'
						message += ' stopping now: name = '+name+'\n'
						error.append(message)
				elif(natoms > 1):
					error,self._element[i] = self.get_element(name,resname)

			if(self._element[i] not in unique_elements): unique_elements.append(self._element[i])
#
###	OPEN	Error exception handling stub
#

		if(len(error) > 0):
			print error
			sys.exit()

		return unique_elements

	def get_element(self,name,resname):
		'''
		Get elment from charmm27 atom name list plus extras from periodic table
		and resolve name conflicts	
		'''

		error = []
		element_name = ''
	
		hydrogen,carbon,nitrogen,oxygen,sulfur,phosphorus,other = self.charmm_names()	
	
		conflict_atoms = ['CD','CE', 'HE','HG','NE','ND','NB','PB','PA']
	
		if(name in conflict_atoms and (name == resname)):
				element_name = name	
		elif(name in hydrogen): element_name = 'H'
		elif(name in carbon): element_name = 'C'
		elif(name in nitrogen): element_name = 'N'
		elif(name in oxygen): element_name = 'O'
		elif(name in sulfur): element_name = 'S'
		elif(name in phosphorus): element_name = 'P'
		elif(name in other): element_name = name
		elif(name == '1H'): element_name = name
		elif(name == '2H'): element_name = "D"
		elif(name == "SOD"): element_name = 'NA'
		elif(name == "POT"): element_name = 'K'
		elif(name == "CAL"): element_name = 'CA'
		elif(name == "CLA"): element_name = 'CL'
		elif(name == "CES"): element_name = 'CS'
		elif(name == "F2'"): element_name = 'F'
		else:
			message = '\nELEMENT NAME NOT IN CHARMM DATABASE\n'
			message += ' stopping now: name = '+name+'\n'
			error.append(message)
		
		return error,element_name

	def write_pdb(self,filename,frame,flag):
		'''
		This method writes the PDB file
		'''
		debug=0
		result=1
	
		if(flag=='w' or flag=='W'):
			infile=open(filename,'w')
		elif(flag=='a' or flag=='A'):				
			infile=open(filename,'a')
	
		if(debug==1):
			infile.write("          1         2         3         4         5         6         7         \n")
			infile.write("01234567890123456789012345678901234567890123456789012345678901234567890123456789\n")

		for i in range(len(self._atom)):
			infile.write("%-6s%5s %-4s%1s%-4s%1s%4s%1s   %8.3f%8.3f%8.3f%6s%6s      %-4s%2s%2s\n" % (self._atom[i],self._index[i],self._name[i],self._loc[i],self._resname[i],self._chain[i],self._resid[i],self._rescode[i],self._coor[frame,i,0],self._coor[frame,i,1],self._coor[frame,i,2],self._occupancy[i],self._beta[i],self._segname[i],self._element[i],self._charge[i]))
		infile.write("END\n")
		infile.close()
	
		return result

	def get_resnames(self):
		'''
		This method holds names of residues to use to set moltype.  Based on Charmm 27 naming.
		'''

		protein_resnames=['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','HSD','HSE','HSP','ILE','LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL']
	
		dna_resnames=['NUSA','NUSG','NUSC','NUSU']

		rna_resnames=['RNUS','RNUA','RUUG','RNUC']

		nucleic_resnames = ['GUA','ADE','CYT','THY','URA','G', 'A', 'C', 'T', 'U']
		
		water_resnames=['TIP3','SPCE','TIP','SPC','TIP4','TP3M'] 
			
		return protein_resnames,dna_resnames,rna_resnames,nucleic_resnames,water_resnames

	def initialize_children(self):

		'''
			This initializes the following "children" with their
			masks already defined to the "parent" molecule

			names() 		: names_mask()
			resnames()		: resnames_mask()
			resids()		: resids_mask()
			chains()		: chains_mask()
			segnames()		: segnames_mask()
			occupancies()	: occupancies_mask()
			betas()		: betas_mask()
			elements()		: elements_mask()
			
			The objects on the left contain the unique values and
			the objects on the right contain the masks that have
			the indices to extract the information for each unique
			value from the parent molecule.

		'''

		number_of_names = self.number_of_names() ; unique_names = self.names()
		number_of_resnames = self.number_of_resnames() ; unique_resnames = self.resnames()
		number_of_resids = self.number_of_resids() ; unique_resids = self.resids()
		number_of_chains = self.number_of_chains() ; unique_chains = self.chains()
		number_of_occupancies = self.number_of_occupancies() ; unique_occupancies = self.occupancies()
		number_of_betas = self.number_of_betas() ; unique_betas = self.betas()
		number_of_elements = self.number_of_elements() ; unique_elements = self.elements()
		number_of_segnames = self.number_of_segnames() ; unique_segnames = self.segnames()

		natoms = self.natoms()

		name = self.name() ; resname = self.resname() ; resid = self.resid() ; chain = self.chain()
		occupancy = self.occupancy() ; beta = self.beta() ; element = self.element() ; segname = self.segname()

		names = [] ; resnames = [] ; resids = [] ; chains = [] 
		occupancies = [] ; betas = [] ; elements = [] ; segments = []

		names_mask = numpy.zeros((number_of_names,natoms),numpy.int)
		resnames_mask = numpy.zeros((number_of_resnames,natoms),numpy.int)
		resids_mask = numpy.zeros((number_of_resids,natoms),numpy.int)
		chains_mask = numpy.zeros((number_of_chains,natoms),numpy.int)
		occupancies_mask = numpy.zeros((number_of_occupancies,natoms),numpy.int)
		betas_mask = numpy.zeros((number_of_betas,natoms),numpy.int)
		elements_mask = numpy.zeros((number_of_elements,natoms),numpy.int)
		segnames_mask = numpy.zeros((number_of_segnames,natoms),numpy.int)

		nresid = 0

		for i in xrange(natoms):

			names_mask[unique_names.index(name[i])][i] = 1
			resnames_mask[unique_resnames.index(resname[i])][i] = 1
			resids_mask[unique_resids.index(resid[i])][i] = 1
			chains_mask[unique_chains.index(chain[i])][i] = 1
			occupancies_mask[unique_occupancies.index(occupancy[i])][i] = 1
			betas_mask[unique_betas.index(beta[i])][i] = 1
			elements_mask[unique_elements.index(element[i])][i] = 1
			segnames_mask[unique_segnames.index(segname[i])][i] = 1

		self._names_mask = names_mask
		self._resnames_mask = resnames_mask
		self._resids_mask = resids_mask
		self._chains_mask = chains_mask
		self._occupancies_mask = occupancies_mask
		self._betas_mask = betas_mask
		self._elements_mask = elements_mask
		self._segnames_mask = segnames_mask

		return

	def read_pdb(self,filename):
		'''
		This method reads a PDB file.  
		'''
		debug=0
		result=1

		protein_resnames,dna_resnames,rna_resnames,nucleic_resnames,water_resnames = self.get_resnames()


	### LONG	Need to properly import header information
	
	# see: http://deposit.rcsb.org/adit/docs/pdb_atom_format.html

		infile=open(filename,'r').readlines()
		print 'reading filename: ',filename
		
		atom=[] ; index=[] ; name=[] ; loc=[] ; resname=[] ; chain=[] ; resid=[] ; rescode=[]
		x=[] ; y=[] ; z=[]
		occupancy=[] ; beta=[] ; segname=[] ; element=[] ; charge=[] ; moltype=[]
	
	# first check to see how many frames are in the file

		num_model = 0
		num_endmdl = 0
		num_end = 0

		num_frames = 1

		count_index = 0

		num_counts_this_model = 0 # number of atoms encompassed by "MODEL" and "ENDMDL" lines
		num_counts_this_end = 0 # number of atoms before the first "END" line or between two consecutive "END" lines
		num_counts_per_model = [] # array of number of atoms encompassed by "MODEL" and "ENDMDL" lines
		num_counts_per_end = [] # array of number of atoms before the first "END" line or between two consecutive "END" lines
		modelON = False # Flag to indicate that a "MODEL" frame is being read
		for i in range(len(infile)):
			lin=infile[i]
			lins=string.split(infile[i])
			record_name = string.strip(lin[0:6])
#
### 	OPEN 	need to re-factor the exception statements to a uniform reporting mechanism
#
			if(lins[0]=='MODEL'):
				if modelON:
					raise Exception, 'Encountered two consecutive MODEL lines' 
				if (num_counts_this_model != 0):
					raise Exception, 'There should not be atoms after ENDMDL and before MODEL lines'
				modelON = True
			elif(lins[0]=='ENDMDL'):
				if not modelON:
					raise Exception, 'Encountered two consecutive ENDMDL lines'
				modelON = False
				num_counts_per_model.append(num_counts_this_model)
				num_counts_this_model = 0
			elif(lins[0]=='END'):
				num_counts_per_end.append(num_counts_this_end)
				num_counts_this_end = 0
			#
			if((record_name == 'ATOM' or record_name == 'HETATM')):
				count_index += 1
				num_counts_this_model += 1
				num_counts_this_end += 1
		#
		#
		if ( (len(num_counts_per_end)==0) and (len(num_counts_per_model)!=0) ):
			raise Exception, 'According to Protein Data Bank Contents Guide, END line must appear in each coor entry'
		if (len(num_counts_per_model)!=0 and (len(num_counts_per_end)>1 or sum(num_counts_per_model)!=sum(num_counts_per_end))):
			print num_counts_per_model,num_counts_per_end
			raise Exception, 'Only one terminating END line is allowed for pdb entries with multiple MODEL'
		#
		if (len(num_counts_per_model)>0):
			num_frames = len(num_counts_per_model)
			num_atoms = num_counts_per_model[0]
			if not all(x == num_atoms for x in num_counts_per_model):
				raise Exception, 'number of atoms per frame is not equal'
		elif (len(num_counts_per_end)>0):
			num_frames = len(num_counts_per_end)
			num_atoms = num_counts_per_end[0]
			if not all(x == num_atoms for x in num_counts_per_end):
				raise Exception, 'number of atoms per frame is not equal'
      		elif ( (len(num_counts_per_model)==0) and (len(num_counts_per_end)==0) ):
         		num_frames = 1
         		num_atoms = num_counts_this_model
      		else:
         		raise Exception, 'unexpected error!'
         

		print 'num_atoms = ',num_atoms

		print '>>> found ',num_frames,' model(s) or frame(s)'			

		this_frame = 1		

		true_index = 0
	
		coor=numpy.zeros((num_frames,num_atoms,3),numpy.float)	


		unique_names = [] ; unique_resnames = [] ; unique_resids = [] ; unique_chains = [] 
		unique_occupancies = [] ; unique_betas = [] ; unique_segnames = []

		for i in range(len(infile)):
			lin=infile[i]
			lins=string.split(infile[i])

			record_name = string.strip(lin[0:6])

			if((record_name == 'ATOM' or record_name == 'HETATM') and this_frame == 1):
				true_index += 1
				atom.append(string.strip(lin[0:6]))		#	1-6		record name	
				#index.append(lin[6:11])				#	7-11		atom serial number
				index.append(str(true_index))   	        #   set index so that > 99,999 atoms can be read and counted
				this_name = string.strip(lin[12:16])		#	13-16		atom name
				name.append(string.strip(lin[12:16]))		#	13-16		atom name
				#loc.append(lin[17])				#	17		alternate location indicator
				loc.append(' ')				#	17		alternate location indicator
				this_resname = string.strip(lin[17:21])	#	18-20		residue name
				resname.append(string.strip(lin[17:21]))	#	18-20		residue name
				this_chain = lin[21]				#	22		chain identifier
				chain.append(lin[21])				#	22		chain identifier
				this_resid = locale.atoi(lin[22:26])			#	23-26		residue sequence number
				resid.append(lin[22:26])			#	23-26		residue sequence number
				rescode.append(lin[26])				#	27		code for insertion of residues
				x.append(lin[30:38])				#	31-38		Real(8.3) X: angstroms	
				y.append(lin[38:46])				#	39-46		Real(8.3) Y: angstroms	
				z.append(lin[46:54])				#	47-54		Real(8.3) Z: angstroms	

				try:	
					occupancy.append(string.strip(lin[54:60]))      #	55-60		occupancy
					this_occupancy =string.strip(lin[54:60])      #	55-60		occupancy
					if(occupancy[-1] == ''):
						occupancy[-1] = "  0.00"
						this_occupancy[-1] = "  0.00"
				except:
					occupancy.append("  0.00")
					this_occupancy = "  0.00"
				try:
					beta.append(string.strip(lin[60:66]))		#	61-66		temperature factor
					this_beta = string.strip(lin[60:66])		#	61-66		temperature factor
					if(beta[-1] == ''):
						beta[-1] = "  0.00"
				except:
					beta.append("  0.00")
					this_beta = "  0.00"
				try:
					segname.append(string.strip(lin[72:76]))	#	73-76		segment identifier
					this_segname = string.strip(lin[72:76])	#	73-76		segment identifier
					if(segname[-1] == ''):
						segname[-1] = "DUM0"
						this_segname = "DUM0"
				except:
					this_segname = "DUM0"
					segname.append("DUM0")
				try:
					element.append(string.strip(lin[76:78]))	#	77-78		element symbol
					if(element[-1] == ''):
						element[-1] = "  "
				except:
					element.append("  ")
				try:
					charge.append(string.strip(lin[78:80]))		#	79-80		charge on the atom
					if(charge[-1] == ''):
						charge[-1] = "  "
				except:
					charge.append("  ")

				if(this_name not in unique_names): unique_names.append(this_name)
				if(this_resname not in unique_resnames): unique_resnames.append(this_resname)
				if(this_resid not in unique_resids): unique_resids.append(this_resid)
				if(this_chain not in unique_chains): unique_chains.append(this_chain)
				if(this_segname not in unique_segnames): unique_segnames.append(this_segname)
				if(this_occupancy not in unique_occupancies): unique_occupancies.append(this_occupancy)
				if(this_beta not in unique_betas): unique_betas.append(this_beta)
	
				this_resname=(string.strip(lin[17:21]))
				if this_resname in protein_resnames:
					moltype.append('protein')
				elif this_resname in dna_resnames:	
					moltype.append('dna')
				elif this_resname in rna_resnames:	
					moltype.append('rna')
				elif this_resname in nucleic_resnames:	
					moltype.append('nucleic')
				elif this_resname in water_resnames:	
					moltype.append('water')
				else:
					moltype.append('other')
					
				if(true_index == num_atoms):				
					print 'finished reading frame = ',this_frame
					index=numpy.array(index,numpy.int)
					resid=numpy.array(resid,numpy.int)
					x=numpy.array(x,numpy.float32)
					y=numpy.array(y,numpy.float32)
					z=numpy.array(z,numpy.float32)
					coor[0,:,0]=x.astype(numpy.float) ; coor[0,:,1]=y.astype(numpy.float) ; coor[0,:,2]=z.astype(numpy.float)
					true_index = 0
					x=[] ; y=[] ; z=[]
					if(this_frame == 1):
						self._atom=atom ; self._index=index ; self._name=name ; self._loc=loc ; self._resname=resname 
						self._chain=chain ; self._resid=resid ; self._rescode=rescode 
						self._occupancy=occupancy ; self._beta=beta ; self._segname=segname ; self._element=element
						self._charge=charge ; self._moltype=moltype

						self._number_of_names = len(unique_names) ; self._names = unique_names
						self._number_of_resnames = len(unique_resnames) ; self._resnames = unique_resnames
						self._number_of_resids = len(unique_resids) ; self._resids = unique_resids
						self._number_of_chains = len(unique_chains) ; self._chains = unique_chains
						self._number_of_segnames = len(unique_segnames) ; self._segnames = unique_segnames
						self._number_of_occupancies = len(unique_occupancies) ; self._occupancies = unique_occupancies
						self._number_of_betas = len(unique_betas) ; self._betas = unique_betas

					this_frame += 1
	
			elif((record_name == 'ATOM' or record_name == 'HETATM') and this_frame > 1):
				true_index += 1
				x.append(lin[30:38])				#	31-38		Real(8.3) X: angstroms	
				y.append(lin[38:46])				#	39-46		Real(8.3) Y: angstroms	
				z.append(lin[46:54])				#	47-54		Real(8.3) Z: angstroms	
		
				if(true_index == num_atoms):
					print 'finished reading frame = ',this_frame
					x=numpy.array(x,numpy.float32)
					y=numpy.array(y,numpy.float32)
					z=numpy.array(z,numpy.float32)
					coor[this_frame-1,:,0]=x ; coor[this_frame-1,:,1]=y ; coor[this_frame-1,:,2]=z
					true_index = 0
					this_frame += 1		
					x=[] ; y=[] ; z=[]

		if((this_frame - 1) != num_frames):
			print '>>> WARNING: pdb file had ',num_frames,' sasio read ',this_frame-1,' frames'
			print '>>> WARNING: pdb file had ',num_frames,' sasio read ',this_frame-1,' frames'
			print '>>> WARNING: pdb file had ',num_frames,' sasio read ',this_frame-1,' frames'

		self._coor=numpy.array(coor)
			
		unique_elements = self.element_filter()

		self._number_of_elements = len(unique_elements) ; self._elements = unique_elements

		self._natoms=len(index)

 		#start_time=time.time()
		#self.initialize_children()
 		#end_time=time.time()
		#print 'elapsed time = ',end_time - start_time

		return result