psfft.py

Tuesday, 08 December 2009 00:00
#!/bin/env python
# -*- coding: utf-8 -*-
#
#--------------------------------------------------------------------------------
#psfft.py v1.2, Copyright Bjoern Olausson
#--------------------------------------------------------------------------------
#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 2 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.
#
#To view the license visit
#http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
#or write to
#Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
#--------------------------------------------------------------------------------
#--------------------------------------------------------------------------------
#
#This program can convert your current PSF/CRD files from all27FF to all36FF without invoking CHARMM.
#Furthermore it can adapt the MASS ID in the PSF file to changes in the RTF topologie file. Use "-s" to
#skipt the tranlsation for the new naming scheme in all36FF if you only want to adjust your MASS IDs
#
##NOTE: It does not check for patches!!!!
import re, time, sys
 
from optparse import OptionParser, OptionGroup
usage='''NOTE, DOES NOT TAKE PATCHES INTO ACCOUNT!'''
opa = OptionParser(usage=usage)
 
opa.add_option("-p", action="store", dest="psf", metavar="FILE",
	help="Old PSF file. [REQUIRED]")
 
opa.add_option("-c", action="store", dest="crd", metavar="FILE",
	help="Old CRD (human readable) file. [OPTIONAL]")
 
opa.add_option("-r", action="append", dest="rtf", metavar="FILE",
	help="Path to top_all36_lipid.rtf file. Default ./top_all36_lipid.rtf [REQUIRED]")
 
opa.add_option("-a", action="append", dest="sub", metavar="FILE",
	help="Supplementary rtf or stream file containing topology information. Can be used more then once. [OPTIONAL]")
 
opa.add_option("-o", action="store", dest="out", metavar="FILE", default="converted-to-c36-with-pssft", type="str",
	help="Name of the new PSF and/or CRD file. Default: converted-to-c36-with-pssft(.crd/.psf)")
 
opa.add_option("-s", action="store_true", dest="skip", default=False, metavar="FILE",
	help="Skip conversion from c27ff to c36ff. Usefull for merged FFs where the MASS ID changed. [OPTIONAL]")
 
# instruct optparse to parse the program's command line
(options, args) = opa.parse_args()
 
if len(sys.argv) == 1:
    opa.print_help()
    sys.exit()
print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
print "!!!!!!!! NOTE, DOES NOT TAKE PATCHES INTO ACCOUNT !!!!!!!!"
print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
# Checking for required otions
if not options.rtf:
	parser.error("-r must be set")
elif not options.psf:
	parser.error("-p must be set")
elif not options.out:
	parser.error("-o must be set")
 
PSF = options.psf
CRD = options.crd
RTF = options.rtf
OUT = options.out
SUP = options.sub
SKIP = options.skip
 
if SUP :
	RTF.extend(SUP)
 
#Regular expression to find DOPC and POPC residues
re_PC = re.compile("[A-Z][A-Z]PC")
 
name_dict = {
	"P1": "P",
	"O2" : "O11",  
	"O1" : "O12",
	"O3" : "O13",
	"O4" : "O14",
	# choline fragment; temp names
	"C11" : "C12", 
	"C12" : "C13",
	# choline; C14 is unchanged
	"C13" : "C14",
	"C14" : "C15",
	"C15" : "C11",
	# methylene protons
	"H11" : "H12A",
	"H12" : "H12B",
	"H51" : "H11A",
	"H52" : "H11B",
	# methyl protons
	"H21" : "H13A",
	"H22" : "H13B",
	"H23" : "H13C",
	"H31" : "H14A",
	"H32" : "H14B",
	"H33" : "H14C",
	"H41" : "H15A",
	"H42" : "H15B",
	"H43" : "H15C"
}
 
def linefilter(lines, start, stop, YIELD_MATCHING=None) :
	"""Filter lines based on regular expressions
	    YIELD_MATCHING can be set to:
	    First - to also yield the first matching line.
	    Last - to also yield the last matching linefilter
	    Both - to yield first and last matching linefilter
	    None - to not yield the matching lines"""
	re_START_MATCH = re.compile(start)
	re_STOP_MATCH = re.compile(stop)
 
	for line in lines :
		if re_START_MATCH.search(line) :
			start_line = line
			if YIELD_MATCHING == "First" or YIELD_MATCHING == "Both":
				yield line
			break
 
	for line in lines:
		if re_STOP_MATCH.search(line) and line != start_line:
			if YIELD_MATCHING == "Last" or YIELD_MATCHING == "Both":
				yield line
			break
		yield line
 
def glob_rtf(FILES):
	start_time = time.time()
	print "\nLoading topologie..."
	# Create a dictionary for MASS
	#{"ATOMTYPE" : ["MASSID", "ATOMTYPE", "MASS"]}
	M_DICT = {}
	# Create a dictionary for the Residues
	#{"RESID" : A_DICT}
	R_DICT = {}
	# Creat a dictionary for the Atoms in the Residue
	#{"ATOMNAME" : ["ATOMTYPE", "MASSID", "ATOMCHG", "MASS"]}
	A_DICT = {}
	counter = 0
	for f in FILES:
		FILE = open(f, 'r')
		for line in FILE:
			counter += 1
			# MASS Section
			if re.match( '^MASS [0-9]+ .*', line ):
				line_list = line.split()
				# MASS 1    HL    1.008000 H ! polar H (equivalent to protein H)
				M_LIST = [ line_list[1], line_list[2], line_list[3], line_list[4] ]
				M_DICT[M_LIST[1]] = M_LIST
			# Create a dictionary for the components
			if re.match( '^RESI [A-Z]+[0-9]*.*', line ):
				# RESI LPPC        0.00 ! deoxylysophosphatidylcholine
				RESIDUE = line.split()[1]
				R_DICT[RESIDUE] = ""
				A_DICT = {}
			if re.match('ATOM [A-Z]+[0-9]*.*', line):
				line_list = line.split()
				## ATOM N    NTL   -0.60 !
				## Get the Atom Type and fetch the MASS-ID and mass from M_DICT
				ATOMNAME = line_list[1]
				ATOMTYPE = line_list[2]
				ATOMCHG = line_list[3]
				## MASSID, ATOMTYPE, MASS, ATOM
				ATOM_INFO = M_DICT[ATOMTYPE]
				MASSID = ATOM_INFO[0]
				MASS = ATOM_INFO[2]
				R_LIST = [ ATOMNAME, ATOMTYPE, MASSID, ATOMCHG, MASS ]
				A_DICT[ATOMNAME] = R_LIST
				R_DICT[RESIDUE] = A_DICT
		FILE.close()
	stop_time = time.time()
	print "Done loading topologie (%.5fs)" %(stop_time - start_time)
	return R_DICT
 
#RTF_DICT = {'POPC' : {H13X': ['HAL2', '6', '0.09', '1.008000'], 'H13Y': ['HAL2', '6', '0.09', '1.008000']...}, 'POPS' : {...}, ...}
RTF_DICT = glob_rtf(RTF)
 
def line_formating_psf(line, ext):
	"""#NOTE: CHARMM Element source/io/psfres.src 1.1
	#-----------------------------------------------------------------------#
	#NOTE: if (qextfmt) then
	#NOTE :c23456789     +123456789+123456789+123456789+123456789+123456789+123456789+
	#NOTE: heading='PSF EXT'
	#NOTE: fmt00='(/I10,A)'
	#NOTE: fmt01='  (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8)'
	#NOTE: fmt01a='(I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8,2G14.6)'
	#NOTE: fmt02='  (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8)'
	#NOTE: fmt02a='(I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8,2G14.6)'
	#NOTE: fmt03='(8I10)'
	#NOTE: fmt04='(9I10)'
	#NOTE: fmt05='(/2I10,A)'
	#NOTE: fmt06='(2I10,3X,L1,3G14.6)'
	#NOTE: else
	#NOTE: heading='PSF'
	#NOTE: fmt00='(/I8,A)'
	#NOTE: fmt01='(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8)'
	#NOTE: fmt01a='(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8,2G14.6)'
	#NOTE: fmt02='(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8)'
	#NOTE: fmt02a='(I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8,2G14.6)'
	#NOTE: fmt03='(8I8)'
	#NOTE: fmt04='(9I8)'
	#NOTE: fmt05='(/2I8,A)'
	#NOTE: fmt06='(2I8,3X,L1,3G14.6)'
	#NOTE: endif
	#-----------------------------------------------------------------------#
 
	#-----------------------------------------------------------------------#
	#NOTE X1 = 1 space
	#NOTE I10 = int of len 10
	#NOTE A8 = ASCII of len 8
	#NOTE F10.5 = float of len 10 with a precision of 5
	#NOTE 2G10.5 = sames as F10.5 but E notation (but what is does the 2G mean?)
	#-----------------------------------------------------------------------#"""
	L = line.split()
	lmass = len(L[7])
 
	# truncate or enlarge the MASS value to exactly 8 chars
	if lmass > 7:
		L[7] = L[7][0:7]
	if lmass < 7:
		nzeros = 7 - lmass
		L[7] = L[7] + '0' * nzeros
	Llengths = len(L)
	if ext:
		# Extended file format
		#           I10    1X    A8      1X    A8    1X    A8  1X    A8   1X    I4     1X1X   2G14.6   1X   I8   1X   2G11.6   1X F10     1X     G16.6
		if Llengths == 11:
			line = '%(ATOMNO)10i %(SEGMENT)-8s %(RESID)-8i %(RES)-8s %(TYPE)-8s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G %(UNK2) 10f %(UNK3) 16.6E\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8]), 'UNK2': float(L[9]), 'UNK3': float(L[10])}
		elif Llengths == 10:
			line = '%(ATOMNO)10i %(SEGMENT)-8s %(RESID)-8i %(RES)-8s %(TYPE)-8s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G %(UNK2) 10f\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8]), 'UNK2': float(L[9])}
		elif Llengths == 9:
			line = '%(ATOMNO)10i %(SEGMENT)-8s %(RESID)-8i %(RES)-8s %(TYPE)-8s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8])}
		elif Llengths == 8:
			line = '%(ATOMNO)10i %(SEGMENT)-8s %(RESID)-8i %(RES)-8s %(TYPE)-8s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7])}
		else:
			raise Exception("Unknown number of elements per line")
	else:
		#normal file format
		#           I8    1X    A4      1X    A4    1X    A4  1X    A4   1X    I4     1X1X   2G14.6   1X   I8   1X   2G11.6   1X F10     1X     G16.6
		if Llengths == 11:
			line = '%(ATOMNO)8i %(SEGMENT)-4s %(RESID)-4i %(RES)-4s %(TYPE)-4s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G %(UNK2) 10f %(UNK3) 16.6E\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8]), 'UNK2': float(L[9]), 'UNK3': float(L[10])}
		elif Llengths == 10:
			line = '%(ATOMNO)8i %(SEGMENT)-4s %(RESID)-4i %(RES)-4s %(TYPE)-4s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G %(UNK2) 10f\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8]), 'UNK2': float(L[9])}
		elif Llengths == 9:
			line = '%(ATOMNO)8i %(SEGMENT)-4s %(RESID)-4i %(RES)-4s %(TYPE)-4s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s %(UNK1)11.6G\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7]), 'UNK1': int(L[8])}
		elif Llengths == 8:
			line = '%(ATOMNO)8i %(SEGMENT)-4s %(RESID)-4i %(RES)-4s %(TYPE)-4s %(MASSID) 4i  %(CHG) -14.6G  %(MASS)7s\n' \
				    %{'ATOMNO': int(L[0]), 'SEGMENT': str(L[1]), 'RESID': int(L[2]), 'RES': str(L[3]), 'TYPE': str(L[4]), 'MASSID': int(L[5]), 'CHG': float(L[6]), 'MASS': str(L[7])}
		else:
			raise Exception("Unknown number of elements per line")
	return line
 
def get_new_parameters(line):
	LINE_LIST = line.split()
	ANAME = LINE_LIST[4]
	RESIDUE = LINE_LIST[3]
	try:
		# ['N', 'NTL', '40', '-0.60', '14.007000']
		NEW_PARA = RTF_DICT[RESIDUE][ANAME]
		NEW_MASSID = NEW_PARA[2]
		NEW_CHARGE = NEW_PARA[3]
		NEW_MASS = NEW_PARA[4]
		NEW_LINE_LIST = LINE_LIST[:5]
		NEW_LINE_LIST.extend([NEW_MASSID, NEW_CHARGE, NEW_MASS])
		NEW_LINE_LIST.extend(LINE_LIST[8:])
		new_parameter = ' '.join(NEW_LINE_LIST)
	except KeyError:
		new_parameter = ' '.join(LINE_LIST)
	return new_parameter
 
 
def get_new_name(line, ANAME):
	if re_PC.search(line) and not SKIP:
		try:
			newa = name_dict[ANAME]
			la = len(ANAME)
			lan = len(newa)
			# The following is actually only relevant for CRD files, sinc I don't want to reasambe the line from scratch like
			# for the PSF. Here we kee the chars which are replaced identical so we don't destroy the line layout'
			# I the new atom name is shorter add blanks to keep the correct number of chars
			if lan < la:
				newa = newa + " " * (la - lan)
			# If the old atom name is longer, add ablank to the old name so the same ammount of chars will be replaced
			if lan > la:
				ANAME = ANAME + " " * (lan - la)
			new_names = re.sub(ANAME, newa, line)
			return new_names
		except KeyError:
			return line
	else:
		return line
 
def crd_conv(CRD, OUT):
	NEW_CRD = open(OUT+".crd", "w")
	start_time = time.time()
	print "Processing %s" %(CRD)
	re_START = re.compile('^ *[0-9]+ *[0-9]+')
	START = False
	STOP = False
	""" Must only be converted when Atom names change, like in the conversion from all27FF to all36FF.
	    #NOTE:ATOMNO RESNO RES TYPE   X          Y          Z   SEGID RESID Weighting
	    #NOTE:   I5  I5 1X A4 1X A4 F10.5      F10.5     F10.5 1X A4 1X A4  F10.5
	    #NOTE:    !    !!   !!   !         !         !         !!   !!   !         !
	    #NOTE: Example:
	    #NOTE:    1    1 GLY  N     -5.37898  10.13102   3.52667 MON1 1      2.23000
	    #NOTE:    2    1 GLY  HT1   -5.89532  10.67297   2.78945 MON1 1      0.00000
	    #NOTE:    3    1 GLY  HT2   -4.49310   9.78471   3.06925 MON1 1      0.00000
	    #NOTE:    4    1 GLY  HT3   -5.10293  10.77107   4.31571 MON1 1      0.00000"""
	START = False
	STOP = False
	for line in open(CRD, 'r'):
		if not START:
			START = re_START.match(line)
			#Check if the file is in the extended format
			if line.find('EXT') != -1:
				ext = True
			else:
				ext = False
			NEW_CRD.write(line)
		elif START:
			START = True
			#Check for an empty line: line.strip() --> False
			if line.strip():
				ANAME = str(line.split()[3])
				new_names = get_new_name(line, ANAME)
				NEW_CRD.write(new_names)
			else:
				NEW_CRD.write(line)
		else:
			NEW_CRD.write(line)
	stop_time = time.time()
	print "Done processing %s (%.5fs)" %(CRD, stop_time - start_time)
 
def psf_conv(PSF, OUT):
	NEW_PSF = open(OUT+".psf", "w")
	start_time = time.time()
	print "\nProcessing %s" %(PSF)
	#Regular expression to find DOPC and POPC residues
	re_PC = re.compile("[A-Z][A-Z]PC")
	re_START = re.compile('[0-9]* !NATOM')
	re_STOP = re.compile('[0-9]* !NBOND')
	START = False
	STOP = False
	for line in open(PSF, 'r'):
		if not START and not STOP:
			START = re_START.search(line)
			#Check if the file is in the extended format
			if line.find('EXT') != -1:
				ext = True
			else:
				ext = False
			NEW_PSF.write(line)
		elif START and not STOP:
			START = True
			STOP = re_STOP.search(line)
			#Check for an empty line: line.strip() --> False
			if not STOP and line.strip():
				ANAME = str(line.split()[4])
				new_names = get_new_name(line, ANAME)
				new_parameter = get_new_parameters(new_names)
				new_line = line_formating_psf(new_parameter, ext)
				NEW_PSF.write(new_line)
			else:
				NEW_PSF.write(line)
		else:
			NEW_PSF.write(line)
	stop_time = time.time()
	print "Done processing %s (%.5fs)\n" %(PSF, stop_time - start_time)
 
if PSF:
	psf_conv(PSF, OUT)
if CRD and not SKIP:
	crd_conv(CRD, OUT)