On 2012-03-23 15:31, Ricardo O. S. Soares wrote:
Hi GMX-users,
Klauda et al (J. Phys. Chem. B, 2012, 116(1), pp 203–210) recently
provided Cholesterol parameters for Charmm FF.
Does anyone know/have a protocol or script to convert the .str file to a
valid .itp file for Charmm within GROMACS?
I understand that Dr. Spoel and colleagues (*J. Am. Chem. Soc.*, 2012,
DOI: 10.1021/ja211929h) are providing cholesterol .itp for OPLS-AA.
Could that maybe be a different starting point?
Thanks for eventual reply,
Here's a script that has been used for this purpose. Please use with
care and *check your output*.
Ricardo.
////
---
Ricardo O. S. Soares , PhD Student.
Group of Biological Physics - Department of Physics & Chemistry
Faculty of Pharmaceutical Sciences at Ribeirão Preto - University of São
Paulo.
Av.do Café, S/N - ZIP:14040-903 - Ribeirão Preto, São Paulo, Brazil.
Phone: +55 16 36024840.
.
.
--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
sp...@xray.bmc.uu.se http://folding.bmc.uu.se
#!/usr/bin/python
"""
Script for parsing charmm force field to gromacs format
inparameters:
command line parameters:
1 charmm topology file
2 corresponding charmm parameter file
3 opt foldername, default cgenff.ff
outfiles:
1 foldername/atomtypes.atp
2 foldername/forcefield.itp
3 foldername/forcefield.doc
4 foldername/aminoacids.rtp
5 foldername/ffbonded.itp
6 foldername/ffnonbonded.itp
7 foldername/forcefield.r2b
8 opt foldername/lipids.rtp (if '!lipid section' statement in CHARMM top file)
9 opt foldername/cmap.itp (if genCMAP = True)
"""
import sys
import math
import re
import os
#------------------
# System parameters
#------------------
# infiles
parFile = open(sys.argv[2], 'r')
topFile = open(sys.argv[1], 'r')
# test to see if there's a name specified
if len(sys.argv)>3:
ffName = sys.argv[3]
print('Creating '+ffName+' files...')
# if not, use cgenff
else:
ffName = 'cgenff-2b7.ff'
print('Creating '+ffName+' files...')
# conversion constant between kcal and kJ
kcal2kJ = 4.184
#-------------------------
# User specific parameters
#-------------------------
# specification of character used for comments in charmm ff file
comment = '!'
# create folder for output files
os.mkdir(ffName)
os.chdir(ffName)
# outfiles
nbFile = open('ffnonbonded.itp', 'w') # nonbonded file
bonFile = open('ffbonded.itp', 'w') # bonded file
atpFile = open('atomtypes.atp', 'w') # atom type file
itpFile = open('forcefield.itp', 'w') # ffcharmm**.itp file
docFile = open('forcefield.doc', 'w') # ffcharmm**.itp file
aartpFile = open('aminoacids.rtp', 'w') # aminoacids rtp file
for line in topFile:
if line.startswith('CMAP'):
print "\n NOTE: This force field seems to support CMAP so trying to port it!\n"
genCMAP = True
cmapFile = open('cmap.itp', 'w') # cmap itp file
break
else:
genCMAP = False
topFile.close()
topFile = open('../'+sys.argv[1], 'r')
# set the func parameter for bonds, angles and proper/improper dihedrals
# for further information see section 5.3.2 in gromacs documentation:
# ftp://ftp.gromacs.org/pub/manual/manual-3.3.pdf
funcForBonds = '1'
funcForAngles = '5' # Urey-Bradley angle type
funcForDihedrals = '9' # special type for treating multiple entries (modification in source code)
funcForImpropers = '2'
funcFor14 = '1' # function
# particle type
ptype = 'A'
# dictionary for atom numbers, used for the nb file
element2atomNumber= {}
element2atomNumber['H']='1'
element2atomNumber['HE']='2'
element2atomNumber['C']='6'
element2atomNumber['N']='7'
element2atomNumber['O']='8'
element2atomNumber['F']='9'
element2atomNumber['NE']='10'
element2atomNumber['NA']='11'
element2atomNumber['MG']='12'
element2atomNumber['AL']='13'
element2atomNumber['P']='15'
element2atomNumber['S']='16'
element2atomNumber['CL']='17'
element2atomNumber['K']='19'
element2atomNumber['CA']='20'
element2atomNumber['Fe']='26'
element2atomNumber['ZN']='30'
element2atomNumber['BR']='35'
element2atomNumber['I']='53'
element2atomNumber['CS']='55'
#-----------------------------------------------------------------------
# parsing the charmm top file and writing to gromacs .atp and .rtp files
#-----------------------------------------------------------------------
# position flags
mass = False
postMass = False
type2element = {}
element2mass = {}
type2charge = {}
firstBond = True
firstImpr = True
presFlag = False
lipidFlag = False
lipidFlagCounter = 0
# group counter
groupCounter = 0
# initiation of rtp file, defaults etc
aartpFile.write('[ bondedtypes ] \n')
aartpFile.write('; Col 1: Type of bond \n')
aartpFile.write('; Col 2: Type of angles \n')
aartpFile.write('; Col 3: Type of proper dihedrals \n')
aartpFile.write('; Col 4: Type of improper dihedrals \n')
aartpFile.write('; Col 5: Generate all dihedrals if 1, only heavy atoms of 0. \n')
aartpFile.write('; Col 6: Number of excluded neighbors for nonbonded interactions \n')
aartpFile.write('; Col 7: Generate 1,4 interactions between pairs of hydrogens if 1 \n')
aartpFile.write('; Col 8: Remove propers over the same bond as an improper if it is 1 \n')
aartpFile.write('; bonds angles dihedrals impropers all_dihedrals nrexcl HH14 RemoveDih \n')
aartpFile.write(' 1 5 9 2 1 3 1 0 \n')
# parse line by line
for line in topFile:
# parse masses and write to gromacs atp file
# match with MASS
if line.startswith('MASS'):
mass = True
# match with something after MASS
if line.startswith('DECL') or line.startswith('DEFA') or line.startswith('AUTO'):
postMass = True
if line.startswith('!lipid section'):
lipidFlagCounter += 1
# match for lipid section in CHARMM top file and create lipids rtp file
if lipidFlagCounter == 2:
lipidFlag = True
lipidrtpFile = open('lipids.rtp', 'w') # lipids rtp file
# initiation of rtp file, defaults etc
lipidrtpFile.write('[ bondedtypes ] \n')
lipidrtpFile.write('; Col 1: Type of bond \n')
lipidrtpFile.write('; Col 2: Type of angles \n')
lipidrtpFile.write('; Col 3: Type of proper dihedrals \n')
lipidrtpFile.write('; Col 4: Type of improper dihedrals \n')
lipidrtpFile.write('; Col 5: Generate all dihedrals if 1, only heavy atoms of 0. \n')
lipidrtpFile.write('; Col 6: Number of excluded neighbors for nonbonded interactions \n')
lipidrtpFile.write('; Col 7: Generate 1,4 interactions between pairs of hydrogens if 1 \n')
lipidrtpFile.write('; Col 8: Remove propers over the same bond as an improper if it is 1 \n')
lipidrtpFile.write('; bonds angles dihedrals impropers all_dihedrals nrexcl HH14 RemoveDih \n')
lipidrtpFile.write(' 1 5 9 2 1 3 1 0 \n')
# not empty line
line = line.strip()
if len(line) > 1:
# mass part
#-----------
if mass and not postMass:
segments = line.split()
# ignore comments
if line[0] != comment:
type = segments[2]
mass = segments[3]
if type[0] == 'H':
element = segments[4]
element = 'H'
elif type[0] == 'B' and type[1] == 'R':
element = 'BR'
elif type[0] == 'C' and type[1] == 'L':
element = 'CL'
elif type[0] == 'A' and type[1] == 'L':
element = 'AL'
elif type[0] == 'C':
element = 'C'
elif type[0] == 'N':
element = 'N'
elif type[0] == 'O':
element = 'O'
elif type[0] == 'S':
element = 'S'
elif type[0] == 'F':
element = 'F'
elif type[0] == 'P':
element = 'P'
elif type[0] == 'I':
element = 'I'
atomComment = segments[6:]
# construct a string from the rest of the elements in the list
string = ''
for word in atomComment:
string = string + word + ' '
# write to the .atp file
atpFile.write(type+'\t'+
mass+' '+
';'+'\t'+
string+'\n')
# set dictionaries:
type2element[type]=element
element2mass[element]=mass
# parse the topologies of the charmm top file and write to gromacs rtp file
#--------------------------------------------------------------------------
# save the comments if available
comments = line.split(' ! ')
# delete comments in the end of the lines
segments = line.split()
line = ''
for seg in segments:
if seg == comment:
break
else:
line = line+str(seg)+'\t'
# new residues to parse
if line.startswith('RESI'):
# reset flags and group counter because this is a new residue
presFlag = False # reset
firstBond = True # reset
firstImpr = True # reset
groupCounter = 0 # reset group counter since new residue
# read and print name of residue
name = line.split()[1]
if not lipidFlag:
aartpFile.write('\n; '+comments[1])
aartpFile.write('\n[ '+name+' ]'+'\n')
aartpFile.write(' [ atoms ]\n') # write header
else:
lipidrtpFile.write('\n[ '+name+' ]'+'\n')
lipidrtpFile.write(' [ atoms ]\n') # write header
# set flag for pres entries, except for ACE and CT2 pres residues
if line.startswith('PRES'):
# find ACE and CT2 PRES residues, set presFlag = False for adding to rtp file
name = line.split()[1]
if name == 'ACE' or name == 'CT2':
presFlag = True # reset
firstBond = True # reset
firstImpr = True # reset
groupCounter = 0 # reset group counter since new residue
# read and print name of residue
if not presFlag:
if not lipidFlag:
aartpFile.write('\n[ '+name+' ]'+'\n')
aartpFile.write(' [ atoms ]\n') # write header
else:
lipidrtpFile.write('\n[ '+name+' ]'+'\n')
lipidrtpFile.write(' [ atoms ]\n') # write header
else:
presFlag = True
# discard if entry is a pres
if not presFlag:
if line.startswith('GROUP'):
#groupCounter += 1
None
if line.startswith('ATOM'):
segments = line.split()
name = segments[1]
type = segments[2]
charge = segments[3]
try:
type2charge[type]=charge
except KeyError:
None
if not lipidFlag:
aartpFile.write('\t'+name+'\t'+type+'\t'+charge+'\t'+str(groupCounter)+'\n')
else:
lipidrtpFile.write('\t'+name+'\t'+type+'\t'+charge+'\t'+str(groupCounter)+'\n')
groupCounter += 1
if line.startswith('BOND') or line.startswith('DOUBLE'):
# several bond statements...
if firstBond:
# write header
if not lipidFlag:
aartpFile.write(' [ bonds ]\n')
else:
lipidrtpFile.write(' [ bonds ]\n')
segments = line.split()
bondNumber = len(segments[1:])/2
for i in range(bondNumber):
atom1 = segments[i+i+1]
atom2 = segments[i+i+2]
if not lipidFlag:
aartpFile.write('\t'+atom1+'\t'+atom2+'\n')
else:
lipidrtpFile.write('\t'+atom1+'\t'+atom2+'\n')
firstBond = False
else:
segments = line.split()
bondNumber = len(segments[1:])/2
for i in range(bondNumber):
atom1 = segments[i+i+1]
atom2 = segments[i+i+2]
if not lipidFlag:
aartpFile.write('\t'+atom1+'\t'+atom2+'\n')
else:
lipidrtpFile.write('\t'+atom1+'\t'+atom2+'\n')
if line.startswith('IMPR'):
if firstImpr:
# write header
if not lipidFlag:
aartpFile.write(' [ impropers ]\n')
else:
lipidrtpFile.write(' [ impropers ]\n')
segments = line.split()
imprNumber = len(segments[1:])/4
for i in range(imprNumber):
atom1 = segments[4*i+1]
atom2 = segments[4*i+2]
atom3 = segments[4*i+3]
atom4 = segments[4*i+4]
if not lipidFlag:
aartpFile.write('\t'+atom1+'\t'+atom2+'\t'+atom3+'\t'+atom4+'\n')
else:
lipidrtpFile.write('\t'+atom1+'\t'+atom2+'\t'+atom3+'\t'+atom4+'\n')
firstImpr = False
else:
segments = line.split()
imprNumber = len(segments[1:])/4
for i in range(imprNumber):
atom1 = segments[4*i+1]
atom2 = segments[4*i+2]
atom3 = segments[4*i+3]
atom4 = segments[4*i+4]
if not lipidFlag:
aartpFile.write('\t'+atom1+'\t'+atom2+'\t'+atom3+'\t'+atom4+'\n')
else:
lipidrtpFile.write('\t'+atom1+'\t'+atom2+'\t'+atom3+'\t'+atom4+'\n')
# parse cmap foursomes
if line.startswith('CMAP'):
segments = line.split()
if genCMAP:
if not lipidFlag:
aartpFile.write(' [ cmap ]\n')
aartpFile.write('\t'+segments[1]+'\t'+segments[2]+'\t'+segments[3]+'\t'+segments[4]+'\t'+segments[8]+'\n')
else:
lipidrtpFile.write(' [ cmap ]\n')
lipidrtpFile.write('\t'+segments[1]+'\t'+segments[2]+'\t'+segments[3]+'\t'+segments[4]+'\t'+segments[8]+'\n')
#--------------------------------------------------------------------
# parsing the charmm par file and writing to gromacs bon and nb files
#--------------------------------------------------------------------
# position flags
bonds = False
angles = False
dihedrals = False
impropers = False
cmap = False
nonbonded = False
hbond = False
# lists for saving 1-4 and LJ params
paramList = []
LJlist = []
# cmap variables
if genCMAP:
cmapType = 1
cmapData = []
cmapParamCounter = 0
# matrix (2D list) for saving cmaps
cmapValues = []
# initiation of temporary storage for dihedrals and impropers containing wildcards
dihWilds = []
impWilds = []
### parse charmm parFile
# parse line by line
for line in parFile:
# match with BONDS
if line.startswith('BONDS'):
bonds = True
# print header
bonFile.write('[ bondtypes ]'+ '\n')
bonFile.write('; i'+'\t'+
'j'+'\t'+
'func'+'\t'+
'b0'+'\t'+
'kb'+'\n')
# match with ANGLES
if line.startswith('ANGLES'):
angles = True
bonds = False
# print header
bonFile.write('\n'+'[ angletypes ]'+'\n')
bonFile.write('; i'+'\t'+
'j'+'\t'+
'k'+'\t'+
'func'+'\t'+
'th0'+'\t'+
'cth'+'\t'+
'ub0'+'\t'+
'cub'+'\n')
# match with (proper) DIHEDRALS
if line.startswith('DIHEDRALS'):
dihedrals = True
angles = False
# print header
bonFile.write('\n'+'[ dihedraltypes ]'+'\n')
bonFile.write('; i\tj\tk\tl\t'+
'func'+'\t'+
'phi0'+'\t'+
'cp'+'\t'
'mult'+'\n')
# match with IMPROPER (dihedrals)
if line.startswith('IMPROPER'):
## print the dihedrals containing wildcards before!
for wilds in dihWilds:
cp = wilds[4]
# conversion to kJ
cp = str(float(cp)*kcal2kJ) # not a factor 2!!!!
bonFile.write(wilds[0]+'\t'+wilds[1]+'\t'+wilds[2]+'\t'+wilds[3]+'\t'+funcForDihedrals+'\t'+wilds[6]+'\t'+cp+'\t'+wilds[5]+'\n')
# continue with the impropers...
impropers = True
dihedrals = False
# print header
bonFile.write('\n'+'[ dihedraltypes ]'+'\n')
bonFile.write('; i'+'\t'+
'j'+'\t'+
'k'+'\t'+
'l'+'\t'+
'func'+'\t'+
'q0'+'\t'+
'cq'+'\n')
# match with CMAP
if line.startswith('CMAP'):
cmap = True
impropers = False
if genCMAP:
cmapFile.write('[ cmaptypes ]')
# match with NONBONDED
if line.startswith('NONBONDED'):
## print the impropers containing wildcards
for wilds in impWilds:
cq = wilds[4]
# converstion to kJ
cq = str(float(cq)*2*kcal2kJ) # Need a factor 2 here too of course!!!
bonFile.write(wilds[0]+'\t'+wilds[1]+'\t'+wilds[2]+'\t'+wilds[3]+'\t'+funcForImpropers+'\t'+wilds[6]+'\t'+cq+'\n')
# continue with the nonbonded...
nonbonded = True
impropers = False
# write last CMAP
if genCMAP:
cmapParamCounter = 0
for i in range(len(cmapData)):
cmapParamCounter = cmapParamCounter + 1
if cmapParamCounter < 11:
cmapFile.write(str(cmapData[i]))
if not cmapParamCounter == 10:
cmapFile.write(' ')
else:
cmapParamCounter = 0
cmapFile.write('\\'+'\n'+str(cmapData[i])+' ')
cmapParamCounter = 1
cmap = False
# print header
nbFile.write('[ atomtypes ]'+'\n')
nbFile.write(';name'+'\t'+
'at.num'+'\t'+
'mass'+'\t'+
'charge'+'\t'+
'ptype'+'\t'+
'sigma'+'\t'+
'epsilon'+'\n')
# match with NBFIX
if line.startswith('NBFIX'):
nonbonded = False
# match with HBOND => after NONBONDED -> 1-4 parameters
if line.startswith('HBOND'):
hbond = True
# write the header for the pairwise 1-4 parameters
nbFile.write('\n[ pairtypes ]'+'\n')
nbFile.write('; i'+'\t'+
'j'+'\t'+
'func'+'\t'+
'sigma1-4'+'\t'+
'epsilon1-4 ; THESE ARE 1-4 INTERACTIONS\n')
###
### Write bonFile
###
# not empty line
line = line.strip()
if len(line) > 0:
# bonds part
#-----------
if bonds and not angles:
segments = line.split()
# ignore comments
if line[0] != comment and segments[0] != 'BONDS':
typei = segments[0]
typej = segments[1]
Kb = segments[2]
# converstion from kcal/mole/A**2 -> kJ/mole/nm**2 incl factor 2 (see definitions)
Kb = str(float(Kb)*2*kcal2kJ*10*10)
b0 = segments[3]
# conversion from A -> nm
b0 = str(float(b0)/10)
bonFile.write(typei+'\t'+
typej+'\t'+
funcForBonds+'\t'+
b0+'\t'+
Kb+'\n')
# angles part, using Urey-Bradley type on all angles (type 5)
#------------------------------------------------------------
if angles and not dihedrals:
segments = line.split()
# ignore comments and the first ANGLES line
if line[0] != comment and segments[0] != 'ANGLES':
typei = segments[0]
typej = segments[1]
typek = segments[2]
th0 = segments [4]
cth = segments [3]
cth = str(float(cth)*2*kcal2kJ) # -> kJ/mol and an factor 2 (see definitions)
# check for Urey-Bradley parameters
if len(segments)>6:
try:
Kub = float(segments[5])*10*10
Kub = Kub*2*kcal2kJ
S0 = float(segments[6])
S0 = S0/10
ubFlag = True
except ValueError:
ubFlag = False
else:
ubFlag = False
if not ubFlag:
Kub = 0.0
S0 = 0.0
## add comments also
#lineComment = segments[6:]
# construct a string from the rest of the elements in the list
#string = ''
#for word in lineComment:
# string = string + word + ' '
bonFile.write(typei+'\t'+
typej+'\t'+
typek+'\t'+
funcForAngles+'\t'+
th0+'\t'+
cth+'\t'+
str(S0)+'\t'+
str(Kub)+'\n') #' ;'+
#string+'\n')
# dihedrals part
#---------------
if dihedrals and not impropers:
segments = line.split()
# ignore comments and the first DIHEDRALS line
if line[0] != comment and segments[0] != 'DIHEDRALS':
typei = segments[0]
typej = segments[1]
typek = segments[2]
typel = segments[3]
# look for wildcards in positions 1 and 4
if typei == 'X' and typel == 'X':
dihWilds.append(segments) # save them in a list
else:
phi0 = segments[6]
cp = segments[4]
# conversion to kJ
cp = str(float(cp)*kcal2kJ)
mult = segments[5]
bonFile.write(typei+'\t'+
typej+'\t'+
typek+'\t'+
typel+'\t'+
funcForDihedrals+'\t'+
phi0+'\t'+
cp+'\t'+
mult+'\n')
# impropers part
#---------------
if impropers and not nonbonded and not cmap:
segments = line.split()
# ignore comments and the first IMPROPERS line
if line[0] != comment and segments[0] != 'IMPROPERS':
typei = segments[0]
typej = segments[1]
typek = segments[2]
typel = segments[3]
# look for wildcards in positions 2 and 3
if typej == 'X' and typek == 'X':
impWilds.append(segments) # save them in a list
else: # no wildcard - write to bon file
q0 = segments [6]
cq = segments [4]
# converstion to kJ
cq = str(float(cq)*2*kcal2kJ) # factor 2 from definition difference
bonFile.write(typei+'\t'+
typej+'\t'+
typek+'\t'+
typel+'\t'+
funcForImpropers+'\t'+
q0+'\t'+
cq+'\n')
# cmap part (new part for the new version of gromacs supporting charmm cmaps)
#----------------------------------------------------------------------------
if genCMAP:
if cmap and not nonbonded:
segments = line.split()
# discard lines starting with a comment
if line[0]!= comment and segments[0]!='CMAP':
# find start of a surface
try: # if float cmap parameters start
segments[0] = float(segments[0])
for i in range(len(segments)):
cmapData.append(float(segments[i])*kcal2kJ)
#cmapFile.write(line+'\n')
except: # if not float it's a string and defines the cmap atom types
# write the parameters
cmapParamCounter = 0
for i in range(len(cmapData)):
cmapParamCounter = cmapParamCounter + 1
if cmapParamCounter < 11:
cmapFile.write(str(cmapData[i]))
if not cmapParamCounter == 10:
cmapFile.write(' ')
else:
cmapParamCounter = 0
cmapFile.write('\\'+'\n'+str(cmapData[i])+' ')
cmapParamCounter = 1
# write the next CMAP header
cmapData = []
cmapFile.write('\n\n'+segments[0]+' '+segments[1]+' '+segments[2]+' '+segments[3]+' '+segments[7]+' '+str(cmapType)+' '+segments[-1]+' '+segments[-1]+'\\'+'\n')
# IF COMMENTS !-180, !-164, etc WANTED
# else:
# # find phi = 0
# try:
# phi = re.search('^!\s*(0)', line).group(1)
# except:
# # find the other phi values
# try:
# phi = re.search('^!\s*(-*[0-9]{2,3})', line).group(1)
# cmapFile.write('!'+phi+'\n')
# except:
# None
###
### Write nbFile
###
# nonbonded part: LJ parameters
#------------------------------
charge = '0.000' # charge for nb definitions
if nonbonded and not hbond:
segments = line.split()
# ignore comments and the first NONBONDED lines
if line[0]!=comment and segments[0]!='NONBONDED' and segments[0]!='cutnb' and segments[0]!='CUTNB':
type = segments[0]
epsilon = segments[2]
eps = str(abs(float(epsilon)*kcal2kJ)) # ->kJ and positive
RminHalf = segments[3]
sigma = str(2*float(RminHalf)/(10.0*2.0**(1.0/6.0))) # -> nm, double distance and rmin2sigma factor
LJlist.append([type, eps, sigma])
# test length to avoid IndexError
if len(segments)> 6:
try: # if possible, convert element 5 to float
segments[5] = float(segments[5])
except:
None
# is segment 5 and 6 floats => there's 1-4 defined
if not isinstance(segments[5],str): # not string?
# read charmm epsilon
epsilon14 = segments[5]
# conversion to gromacs units
eps14 = str(abs(float(epsilon14)*kcal2kJ))
# read charmm Rmin*1/2
Rmin14Half = segments[6]
# conversion to gromacs units
sigma14 = str(2*float(Rmin14Half)/(10.0*2.0**(1.0/6.0)))
# add to list
paramList.append([type, eps14, sigma14])
# test if partial charge is defined
try:
charge = type2charge[type]
noChargeComment = ''
except KeyError:
noChargeComment = '; partial charge def not found'
# add special types of TIP3p model, in both with and without HEAVY_H
noChargeComment = ''
if type == 'HT':
nbFile.write('#ifdef HEAVY_H'+'\n')
nbFile.write(type+'\t'+
element2atomNumber[type2element[type]]+'\t'+
str(4*float(element2mass[type2element[type]]))+'\t'+
charge+'\t'+
ptype+'\t'+
sigma+'\t'+
eps+' '+noChargeComment+'; CHARMM TIP3p H\n')
nbFile.write('#else \n')
nbFile.write(type+'\t'+
element2atomNumber[type2element[type]]+'\t'+
element2mass[type2element[type]]+'\t'+
charge+'\t'+
ptype+'\t'+
sigma+'\t'+
eps+' '+noChargeComment+'\n')
nbFile.write('#endif'+'\n')
else:
if type2element.has_key(type):
print type, " ->", type2element
nbFile.write(type+'\t'+
element2atomNumber[type2element[type]]+'\t'+
element2mass[type2element[type]]+'\t'+
charge+'\t'+
ptype+'\t'+
sigma+'\t'+
eps +
' ' +
# noChargeComment +
'\n')
if type == 'OT':
nbFile.write('#ifdef HEAVY_H'+'\n')
nbFile.write(type+'\t'+
element2atomNumber[type2element[type]]+'\t'+
'9.951400'+'\t'+
charge+'\t'+
ptype+'\t'+
sigma+'\t'+
eps+' '+noChargeComment+'; CHARMM TIP3p O\n')
nbFile.write('#endif'+'\n')
# nonbonded part: 1-4 parameters, where available
# NOTE: not in the for line in par loop since the 1-4 params are saved already
#------------------------------------------------------------------------------
# loop through all atom types with 1-4 parameters specified
for i in range(len(paramList)): # outer loop
# look for other types with 1-4 parameters
for j in range(i, len(paramList)): # inner loop
# use combination rule 2 for epsilon: epsij=sqrt(epsi*epsj)
combEps14 = math.sqrt(float(paramList[i][1])*float(paramList[j][1]))
# use combination rule 2 for sigma
combSigma14 = (float(paramList[i][2])+float(paramList[j][2]))/2.0
# write pairs to [ pairtypes ] section
nbFile.write(paramList[i][0]+'\t'+
paramList[j][0]+'\t'+
funcFor14+'\t'+
str(combSigma14)+'\t'+
str(combEps14)+'\n')
# generate 1-4 params with non-1-4 specified types also
for k in range(len(LJlist)):
not14 = True
# check if atom is present in 1-4 list
for l in range(len(paramList)):
if LJlist[k][0] == paramList[l][0]:
not14 = False
# if not calculate pair 1-4 interactions according to combination rules
if not14:
# use combination rule 2 for epsilon: epsij=sqrt(epsi*epsj)
combEps14 = math.sqrt(float(paramList[i][1])*float(LJlist[k][1]))
# use combination rule 2 for sigma
combSigma14 = (float(paramList[i][2])+float(LJlist[k][2]))/2.0
# write pairs to [ pairtypes ] section
nbFile.write(paramList[i][0]+'\t'+
LJlist[k][0]+'\t'+
funcFor14+'\t'+
str(combSigma14)+'\t'+
str(combEps14)+'\n')
#-----------------------------------------
# construction of the forcefield.itp file
#-----------------------------------------
# '#' is for the c-preprocessor
itpFile.write('*******************************************************************************\n')
itpFile.write('* CHARMM port writted by *\n')
itpFile.write('* Par Bjelkmar, Per Larsson, Michel Cuendet, *\n')
itpFile.write('* Berk Hess and Erik Lindahl. *\n')
itpFile.write('* Correspondance: *\n')
itpFile.write('* bjelk...@cbr.su.se or lind...@cbr.su.se *\n')
itpFile.write('*******************************************************************************\n\n\n')
itpFile.write('#define _FF_CHARMM\n')
itpFile.write('[ defaults ]\n')
itpFile.write('; nbfunc\tcomb-rule\tgen-pairs\tfudgeLJ\tfudgeQQ\n')
itpFile.write('1\t2\tyes\t1.0\t1.0\n\n')
itpFile.write('#include "ffnonbonded.itp"\n')
itpFile.write('#include "ffbonded.itp"\n')
itpFile.write('#include "gb.itp"\n')
if genCMAP:
itpFile.write('#include "cmap.itp"\n')
itpFile.write('; Nucleic acids nonbonded and bonded parameters"\n')
itpFile.write('#include "ffnanonbonded.itp"\n')
itpFile.write('#include "ffnabonded.itp"\n')
#-----------------------------------------
# construction of the forcefield.doc file
#-----------------------------------------
# '#' is for the c-preprocessor
docFile.write('CHARMM27 all-atom force field (with CMAP) - version 2.0\n\n')
docFile.write('*******************************************************************************\n')
docFile.write('* CHARMM port writted by *\n')
docFile.write('* Par Bjelkmar, Per Larsson, Michel Cuendet, *\n')
docFile.write('* Berk Hess and Erik Lindahl. *\n')
docFile.write('* Correspondance: *\n')
docFile.write('* bjelk...@cbr.su.se or lind...@cbr.su.se *\n')
docFile.write('*******************************************************************************\n\n\n')
docFile.write('Parameters derived from c32b1 version of CHARMM\n')
docFile.write('NOTE: Atom-based charge groups\n\n')
docFile.write('References:\n\n')
docFile.write('Proteins\n\n')
docFile.write('MacKerell, Jr., A. D., Feig, M., Brooks, C.L., III, Extending the\n')
docFile.write('treatment of backbone energetics in protein force fields: limitations\n')
docFile.write('of gas-phase quantum mechanics in reproducing protein conformational\n')
docFile.write('distributions in molecular dynamics simulations, Journal of\n')
docFile.write('Computational Chemistry, 25: 1400-1415, 2004.\n\n')
docFile.write('and \n\n')
docFile.write('MacKerell, Jr., A. D., et al. All-atom\n')
docFile.write('empirical potential for molecular modeling and dynamics Studies of\n')
docFile.write('proteins. Journal of Physical Chemistry B, 1998, 102, 3586-3616.\n\n')
docFile.write('Lipids\n\n')
docFile.write('Feller, S. and MacKerell, Jr., A.D. An Improved Empirical Potential\n')
docFile.write('Energy Function for Molecular Simulations of Phospholipids, Journal\n')
docFile.write('of Physical Chemistry B, 2000, 104: 7510-7515.\n\n')
docFile.write('Nucleic Acids\n\n')
docFile.write('Foloppe, N. and MacKerell, Jr., A.D. "All-Atom Empirical Force Field for\n')
docFile.write('Nucleic Acids: 2) Parameter Optimization Based on Small Molecule and\n')
docFile.write('Condensed Phase Macromolecular Target Data. 2000, 21: 86-104.\n\n')
docFile.write('and \n\n')
docFile.write('MacKerell, Jr., A.D. and Banavali, N. "All-Atom Empirical Force Field for\n')
docFile.write('Nucleic Acids: 2) Application to Molecular Dynamics Simulations of DNA\n')
docFile.write('and RNA in Solution. 2000, 21: 105-120.\n\n')
docFile.write('\n\n')
docFile.write('If using there parameters for research please cite:\n')
docFile.write('Bjelkmar, P., Larsson, P., Cuendet, M. A, Bess, B., Lindahl, E.\n')
docFile.write('Implementation of the CHARMM force field in GROMACS: Analysis of protein\n')
docFile.write('stability effects from correction maps, virtual interaction sites, and\n')
docFile.write('water models., Journal of Chemical Theory and Computation, 6: 459-466, 2010.\n')
# Close all files
nbFile.close()
bonFile.close()
atpFile.close()
itpFile.close()
docFile.close()
aartpFile.close()
if genCMAP:
cmapFile.close()
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
