← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/moltemplate/examples/all_atom_examples/force_field_OPLSAA/ethylene+benzene/moltemplate_files/oplsaa_lt_generator/oplsaa_moltemplate.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/moltemplate/examples/all_atom_examples/force_field_OPLSAA/ethylene+benzene/moltemplate_files/oplsaa_lt_generator/oplsaa_moltemplate.py
 
1
#! /usr/bin/env python
 
2
#
 
3
# The purpose of this script is to create a moltemplate lt file for the oplsaa.
 
4
# forcefield.  This will assist researchers in building complex simulations using
 
5
# this OPLS-UA and the OPLS-AA forcefields.
 
6
 
 
7
__author__="Jason Lambert"
 
8
# (some additional corrections by Miguel Gonzalez and Andrew Jewett)
 
9
__version__="0.18"
 
10
 
 
11
import sys
 
12
import os
 
13
from operator import itemgetter
 
14
 
 
15
g_program_name    = __file__.split('/')[-1]
 
16
 
 
17
#  To do that, first make a copy of the \"oplsaa.prm\" file
 
18
#  (which can be downloaded from the TINKER web site).
 
19
#  The lines in this file beginning with the word \"atoms\" should
 
20
#  define the atoms which you plan to put in your simulation. All other 
 
21
#  lines beginning with the word \"atoms\" should be deleted.
 
22
#  (Leave the other sections of this file alone.)
 
23
#""")
 
24
 
 
25
 
 
26
 
 
27
#input data from file containing opls aa force field parameters.
 
28
try:
 
29
   f=open(sys.argv[1],"r")
 
30
except:
 
31
   sys.stderr.write("Error: \n"
 
32
   "    You need to specify a file name as an input argument:\n"
 
33
   "    python oplsaa_moltemplate.py <forcefield file name>\n"
 
34
   "    (or the file name is specified incorrectly)\n")
 
35
   sys.exit()
 
36
 
 
37
 
 
38
sys.stderr.write(g_program_name+", version "+__version__+"\n"
 
39
                 "Reading parameter file...\n")
 
40
 
 
41
 
 
42
#output lt file
 
43
g=open("oplsaa.lt","w")
 
44
 
 
45
 
 
46
 
 
47
lines = f.readlines()
 
48
 
 
49
 
 
50
 
 
51
# Ignore/Comment out lines before the "##  Atom Type Definitions  ##" section.
 
52
 
 
53
for i in range(0, len(lines)):
 
54
   if (lines[i].find("##  Atom Type Definitions  ##") != -1):
 
55
      break
 
56
   else:
 
57
      lines[i] = '# ' + lines[i]
 
58
 
 
59
 
 
60
# As of late 2014, there appear to be 906 atom types, but we don't assume this.
 
61
# First try to infer out how many atom types there were in the original 
 
62
# oplsaa.prm file, or at least find an upper bound on the atom-type numbers.
 
63
# (Keep track of the maximum value of the first column in the "atom" section.)
 
64
max_atomType = 0
 
65
num_atomTypes = 0
 
66
for line in lines:
 
67
   # skip over text after a # comment character
 
68
   ic = line.find('#') 
 
69
   if ic != -1:
 
70
      line = (line[:ic]).strip()
 
71
   else:
 
72
      line = line.strip()
 
73
   # now look for lines beginning with the word "atom"
 
74
   tokens = line.split()
 
75
   if ((len(tokens)>2) and (tokens[0] == "atom")):
 
76
       num_atomTypes += 1
 
77
       if (int(tokens[1]) > max_atomType):
 
78
          max_atomType = int(tokens[1])
 
79
 
 
80
if num_atomTypes > 25:
 
81
   sys.stderr.write("\n"
 
82
      "(Note: If your computer freezes while running "+g_program_name+",\n"
 
83
      "       it could be because you forgot to edit the .prm file.\n"
 
84
      "       The original \"oplsaa.prm\" file distributed with TINKER has over 900 atom\n"
 
85
      "       types.  If you run "+g_program_name+" on this file, it may freeze or\n"
 
86
      "       crash.  Instead, run "+g_program_name+" on a SUBSET of the OPLS atoms\n"
 
87
      "       relevant to your problem.  To do that, delete the lines from the .prm\n"
 
88
      "       file beginning with \"atom\" which you do not need.)\n\n")
 
89
 
 
90
#temporary storage file
 
91
h=open("temp.txt","w+")
 
92
atom_lookup={} #this dictionary contains all the atom ffid's as a key and the number of atoms with that key
 
93
#atom=[[10000,10000] for i in range(906)]  <- don't assume there are 906 atoms
 
94
atom=[[-10000,-10000] for i in range(0,max_atomType+1)]
 
95
#charge_by_type={} # lookup charge by atom type
 
96
#vdw_by_type={}    # lookup epsilon & sigma paramters by atom type
 
97
charge_by_type=[0.0 for i in range(0,max_atomType+1)]    # lookup charge by atom
 
98
vdw_by_type=[(0.0,0.0) for i in range(0,max_atomType+1)] # lookup epsilon & sigma
 
99
 
 
100
 
 
101
 
 
102
#atom is declared this way so for sorting purposes.  
 
103
#atom contains the following data upon allocation
 
104
#atom[][0]=atom_id( Important for partial charges and non_bonded interactions)
 
105
#atom[][1]=atom_ffid( Important for stretches, bending, torsions and impropers)
 
106
#atom[][2]=atom_mass
 
107
#atom[][3]=partial charge
 
108
#atom[][4]=non_bonding sigma
 
109
#atom[][5]=non_bonding epsilon
 
110
#atom[][6]=atom comment
 
111
bond=[]
 
112
#bond contains the following data
 
113
#bond[0]=atom 1 ffid
 
114
#bond[1]=atom 2 ffid
 
115
#bond[2]=bond spring constant(OPLS-aa compatible)
 
116
#bond[3]=equilibrium bond distance(Angstrom)
 
117
angle=[]
 
118
#angle contains the following data
 
119
#angle[0]=atom 1 ffid
 
120
#angle[1]=atom 2 ffid
 
121
#angle[2]=atom 3 ffid
 
122
#angle[3]=spring constant
 
123
#angle[4]=equilibrium angle (degrees)
 
124
dihedral=[]
 
125
#dihedral contains the following data
 
126
#dihedral[0]=atom 1 ffid
 
127
#dihedral[1]=atom 2 ffid
 
128
#dihedral[2]=atom 3 ffid
 
129
#dihedral[3]=atom 4 ffid
 
130
#dihedral[4]=v1
 
131
#dihedral[5]=v2
 
132
#dihedral[6]=v3
 
133
#dihedral[7]=v4
 
134
improper=[]
 
135
#improper[0]=atom 1 ffid
 
136
#improper[1]=atom 2 ffid(central atom)
 
137
#improper[2]=atom 3 ffid
 
138
#improper[3]=atom 4 ffid
 
139
#improper[4]=spring coefficient
 
140
#improper[5]=equilibrium angle
 
141
 
 
142
 
 
143
#This section gets all the parameters from the force field file
 
144
for line in lines:
 
145
 
 
146
   # skip over text after a # comment character
 
147
   ic = line.find('#') 
 
148
   if ic != -1:
 
149
      line = (line[:ic]).strip()
 
150
   else:
 
151
      line = line.strip()
 
152
 
 
153
   if line.find("atom") == 0:
 
154
      line=line.split()
 
155
      atom[int(line[1])-1]=[int(line[1]),int(line[2]),float(line[-2]),
 
156
                            0.0,0.0,0.0," ".join(line[3:-2])]    
 
157
   elif line.find("vdw") == 0:
 
158
      line=line.split()
 
159
      #vdw_temp.append([float(line[1]),float(line[2]),float(line[3])])
 
160
      if (int(line[1]) <= max_atomType):
 
161
         vdw_by_type[int(line[1])] = (float(line[2]),float(line[3]))
 
162
   elif line.find("bond") == 0:
 
163
      line=line.split()
 
164
      bond.append([int(line[1]),int(line[2]),float(line[3]),float(line[4])])
 
165
   elif line.find("angle") == 0:
 
166
      line=line.split()
 
167
      angle.append([int(line[1]),int(line[2]),int(line[3]),
 
168
      float(line[4]),float(line[5])])
 
169
   elif line.find("torsion") == 0:
 
170
      line=line.split()
 
171
      dihedral.append([int(line[1]),int(line[2]),int(line[3]),int(line[4]),
 
172
      float(line[5]),float(line[8]), float(line[11]), 0.0])
 
173
   elif line.find("charge") == 0:
 
174
      line=line.split()
 
175
      #charge_temp.append([int(line[1]),float(line[2])])
 
176
      if (int(line[1]) <= max_atomType):
 
177
         charge_by_type[int(line[1])] = float(line[2])
 
178
   elif line.find("imptors") == 0:
 
179
      line=line.split()
 
180
      improper.append([int(line[1]), int(line[2]), 
 
181
      int(line[3]), int(line[4]), float(line[5]), float(line[6])])
 
182
 
 
183
if len(atom) > 600:
 
184
   sys.stderr.write("WARNING: The number of atom types in your file exceeds 600\n"
 
185
                    "         (You were supposed to edit out the atoms you don't need.\n"
 
186
                    "          Not doing this may crash your computer.)\n"
 
187
                    "\n"
 
188
                    "         Proceed? (Y/N): ")
 
189
   reply = sys.stdin.readline()
 
190
   if find(reply.strip().lower(), 'y') != 0:
 
191
      exit(0)
 
192
 
 
193
#adding the charge and Lennard Jones parameters to
 
194
#to each atom type.
 
195
#----------------------------------------------#
 
196
for i in range(0,len(atom)):
 
197
   atom_type_num = atom[i][0]
 
198
   #q = charge_by_type.get(atomTypeNum)
 
199
   #if q:
 
200
   #   atom[i][3] = q
 
201
   if atom_type_num != -10000:
 
202
      q = charge_by_type[atom_type_num]
 
203
      atom[i][3] = q
 
204
 
 
205
for i in range(0,len(atom)):
 
206
   atom_type_num = atom[i][0]
 
207
   #vdw_params = vdw_by_type.get(atomTypeNum)
 
208
   #if vdw_params:
 
209
   #   atom[i][4] = vdw_params[0]
 
210
   #   atom[i][5] = vdw_params[1]
 
211
   if atom_type_num != -10000:
 
212
      vdw_params = vdw_by_type[atom_type_num]
 
213
      atom[i][4] = vdw_params[0]
 
214
      atom[i][5] = vdw_params[1]
 
215
   
 
216
del(charge_by_type)
 
217
del(vdw_by_type)
 
218
 
 
219
#----------------------------------------------------------#
 
220
#begin writing content to lt file
 
221
g.write("OPLSAA {\n\n" )
 
222
 
 
223
#write out the atom masses
 
224
#----------------------------------------------------------#
 
225
g.write("  write_once(\"Data Masses\"){\n")#checked with gaff
 
226
for i,x in enumerate(atom):
 
227
   if x[0] != -10000:
 
228
     g.write("    @atom:{} {} #{} partial charge={}\n".format(
 
229
     x[0],x[2],x[6],x[3]))
 
230
g.write("  } #(end of atom masses)\n\n")
 
231
#----------------------------------------------------------#
 
232
 
 
233
 
 
234
#write out the pair coefficients
 
235
#----------------------------------------------------------#
 
236
g.write("  write_once(\"In Settings\"){\n")#checked with gaff
 
237
for i,x in enumerate(atom):
 
238
  if x[0] != -10000:
 
239
   g.write("    pair_coeff @atom:{0} @atom:{0} lj/cut/coul/long {1} {2}\n".format(x[0],x[5],x[4]))
 
240
g.write("  } #(end of pair coeffs)\n\n")
 
241
 
 
242
g.write("  write_once(\"In Charges\"){\n")#checked with gaff
 
243
for i,x in enumerate(atom):
 
244
  if x[0] != -10000:
 
245
   g.write("    set type @atom:{0} charge {1}\n".format(x[0],x[3]))
 
246
g.write("  } #(end of atom charges)\n\n")
 
247
 
 
248
#-----------------------------------------------------------#
 
249
 
 
250
# This part of the code creates a lookup dictionary
 
251
# that allows you to find every type of atom by its 
 
252
# force field id. force field id is the id number
 
253
# relevant to bonds, angles, dihedrals, and impropers.
 
254
# This greatly increases the speed of angle, bond, dihedral
 
255
# and improper assignment.
 
256
#------------------------------------------------------------#
 
257
atom=sorted(atom,key=itemgetter(1))
 
258
atom_ffid=0
 
259
for x in atom:
 
260
        if x[1]==atom_ffid:
 
261
                atom_lookup[x[1]].append(x[0])
 
262
        elif x[1]>atom_ffid:
 
263
           atom_lookup[x[1]]=[x[0]]
 
264
           atom_ffid=x[1]
 
265
atom_lookup[0]=["*"]
 
266
#-------------------------------------------------------------#
 
267
 
 
268
 
 
269
#writing out the bond coefficients and bond parameters#
 
270
#-------------------------------------------------------------#
 
271
g.write("  write_once(\"In Settings\") {\n")
 
272
index1=0
 
273
for x in bond:
 
274
  for y in atom_lookup.get(x[0],[]):
 
275
    for z in atom_lookup.get(x[1],[]):
 
276
      #g.write("    bond_coeff @bond:{}-{} harmonic {} {}\n".format(y,z,x[2]/2,x[3]))
 
277
      # Miguel Gonzales corrected this line to:
 
278
      g.write("    bond_coeff @bond:{}-{} harmonic {} {}\n".format(y,z,x[2],x[3]))
 
279
      h.write("    @bond:{0}-{1} @atom:{0} @atom:{1}\n".format(y,z))
 
280
g.write("  } #(end of bond_coeffs)\n\n")
 
281
h.seek(0,0)      
 
282
g.write("  write_once(\"Data Bonds By Type\") {\n")
 
283
for line in h.readlines():
 
284
  g.write(line)
 
285
g.write("  } #(end of bonds by type)\n\n")
 
286
del(bond)
 
287
h.close()
 
288
#-----------------------------------------------------------#
 
289
h=open("temp.txt","w+")
 
290
 
 
291
#writing out angle coefficients and angles by type.---------#
 
292
#-----------------------------------------------------------#  
 
293
g.write("  write_once(\"Data Angles By Type\"){\n") 
 
294
for x in angle:
 
295
  for y in atom_lookup.get(x[0],[]):
 
296
    for z in atom_lookup.get(x[1],[]):
 
297
      for u in atom_lookup.get(x[2],[]):
 
298
         #print(y,z,u,x)
 
299
         #h.write("    angle_coeff @angle:{}-{}-{} harmonic {} {}\n".format(y,z,u,x[3]/2.0,x[4]))
 
300
         # Miguel Gonzales corrected this line:
 
301
         h.write("    angle_coeff @angle:{}-{}-{} harmonic {} {}\n".format(y,z,u,x[3],x[4]))
 
302
         g.write("    @angle:{0}-{1}-{2} @atom:{0} @atom:{1} @atom:{2}\n".format(y,z,u))
 
303
 
 
304
 
 
305
g.write("  } #(end of angles by type)\n\n") 
 
306
h.seek(0,0)
 
307
g.write("  write_once(\"In Settings\" ){\n")
 
308
for line in h.readlines():
 
309
        g.write(line)
 
310
g.write("  } #(end of angle_coeffs)\n\n")                    
 
311
del(angle)           
 
312
h.close()
 
313
#----------------------------------------------------------#
 
314
 
 
315
#writing dihedrals by type and dihedral coefficients-------#
 
316
h=h=open("temp.txt","w+")
 
317
g.write("  write_once(\"Data Dihedrals By Type\") {\n")
 
318
#print(atom_lookup)
 
319
for x in dihedral:
 
320
  for y in atom_lookup.get(x[0],[]):
 
321
    for z in atom_lookup.get(x[1],[]):
 
322
      for u in atom_lookup.get(x[2],[]): 
 
323
        for v in atom_lookup.get(x[3],[]):
 
324
          if x[0]!=0 and x[3]!=0:
 
325
            g.write("    @dihedral:{0}-{1}-{2}-{3} @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
 
326
            y,z,u,v))
 
327
            h.write("    dihedral_coeff @dihedral:{}-{}-{}-{} opls {} {} {} {}\n".format(
 
328
            y,z,u,v,x[4],x[5],x[6],x[7]))
 
329
          elif x[0]==0 and x[3]!=0:
 
330
            g.write("    @dihedral:0-{1}-{2}-{3} @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
 
331
            y,z,u,v))
 
332
            h.write("    dihedral_coeff @dihedral:0-{}-{}-{} opls {} {} {} {}\n".format(
 
333
            z,u,v,x[4],x[5],x[6],x[7]))
 
334
          elif x[0]==0 and x[3]==0:
 
335
            g.write("    @dihedral:0-{1}-{2}-0 @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(
 
336
            y,z,u,v))
 
337
            #h.write("    dihedral_coeff @dihedral:0-{}-{}-0 harmonic {} {} {} {}\n".format(
 
338
            h.write("    dihedral_coeff @dihedral:0-{}-{}-0 opls {} {} {} {}\n".format(
 
339
            z,u,x[4],x[5],x[6],x[7]))       
 
340
 
 
341
del(dihedral)
 
342
g.write("  } #(end of Dihedrals by type)\n\n")               
 
343
h.seek(0,0)   
 
344
g.write("  write_once(\"In Settings\") {\n") 
 
345
for line in h.readlines():
 
346
   g.write(line)
 
347
g.write("  } #(end of dihedral_coeffs)\n\n")
 
348
h.close()
 
349
#-----------------------------------------------------------------------#
 
350
 
 
351
#----writing out improper coefficients and impropers by type------------# 
 
352
h=open("temp.txt","w+")
 
353
g.write("  write_once(\"Data Impropers By Type (opls_imp.py)\") {\n")
 
354
for x in improper:
 
355
  for y in atom_lookup.get(x[0],[]):
 
356
    for z in atom_lookup.get(x[1],[]):
 
357
      for u in atom_lookup.get(x[2],[]): 
 
358
        for v in atom_lookup.get(x[3],[]):
 
359
         # Notation: let I,J,K,L denote the atom types ("biotypes") 
 
360
         #  listed in the order they appear in the "oplsaa.prm" file.
 
361
         # (I think J and L are represented by "u" and "v" in the code here.)
 
362
         # It looks like the "oplsaa.prm" file distributed with tinker
 
363
         # treats the third atom ("K") as the central atom.
 
364
         # After checking the code, it appears that the improper angle is
 
365
         # calculated as the angle between the I,J,K and the J,K,L planes
 
366
         if x[0]==0 and x[1]==0 and x[3]==0:
 
367
            g.write("    @improper:0-0-{2}-0 @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(y,z,u,v))
 
368
            h.write("    improper_coeff @improper:0-0-{2}-0 harmonic {4} {5} \n".format(y,z,u,v,x[4]/2,180))
 
369
         else:
 
370
            g.write("    @improper:0-0-{2}-{3} @atom:{0} @atom:{1} @atom:{2} @atom:{3}\n".format(y,z,u,v))
 
371
            h.write("    improper_coeff @improper:0-0-{2}-{3} harmonic {4} {5} \n".format(y,z,u,v,x[4]/2,180))
 
372
 
 
373
 
 
374
g.write("  } #(end of impropers by type)\n\n") 
 
375
h.seek(0,0)   
 
376
g.write("  write_once(\"In Settings\") {\n") 
 
377
for line in h.readlines():
 
378
   g.write(line)
 
379
g.write("  } #(end of improp_coeffs)\n\n")
 
380
#-----------------------------------------------------------------------#
 
381
 
 
382
#This section writes out the input parameters required for an opls-aa simulation
 
383
# lammps.
 
384
g.write("  write_once(\"In Init\") {\n")
 
385
g.write("    units real\n")
 
386
g.write("    atom_style full\n")
 
387
g.write("    bond_style hybrid harmonic\n")
 
388
g.write("    angle_style hybrid harmonic\n")
 
389
g.write("    dihedral_style hybrid opls\n")
 
390
g.write("    improper_style hybrid harmonic\n")
 
391
#g.write("    pair_style hybrid lj/cut/coul/cut 10.0 10.0\n")
 
392
g.write("    pair_style hybrid lj/cut/coul/long 10.0 10.0\n")
 
393
g.write("    pair_modify mix geometric\n")
 
394
g.write("    special_bonds lj/coul 0.0 0.0 0.5\n")    
 
395
g.write("    kspace_style pppm 0.0001\n")
 
396
g.write("  } #end of init parameters\n\n")
 
397
g.write("} # OPLSAA\n")
 
398
f.close()
 
399
g.close()
 
400
h.close()
 
401
os.remove("temp.txt")
 
402
 
 
403
 
 
404
sys.stderr.write("...finished.\n")
 
405
 
 
406
 
 
407
 
 
408