← Back to branch summary

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

« back to all changes in this revision

Viewing changes to tools/amber2lmp/amber2lammps.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/amber2lmp/amber2lammps.py
1
 
#! /usr/bin/python
2
 
 
3
 
#
4
 
# This is amber2lammps, a program written by Keir E. Novik to convert
5
 
# Amber files to Lammps files.
6
 
#
7
 
# Copyright 1999, 2000 Keir E. Novik; all rights reserved.
8
 
# 
9
 
# Modified by Vikas Varshney, U Akron, 5 July 2005, as described in README
10
 
# Bug Fixed :Third argument in Dihedral Coeffs section is an integer - Ketan S Khare September 26, 2011
11
 
# Modified by Vikas Varshney, Oct 8, 2013 to include additional flags (Atomic_Number, Coulombic and van der Waals 1-4 factors which are included in newer vesions of .top and .crd files in amber12. 
12
 
 
13
 
#============================================================
14
 
 
15
 
def Pop(S, I=-1):
16
 
    'Pop item I from list'
17
 
    X = S[I]
18
 
    del S[I]
19
 
    return X
20
 
 
21
 
#============================================================
22
 
 
23
 
class Lammps:
24
 
 
25
 
    #--------------------------------------------------------
26
 
 
27
 
    def Dump(self):
28
 
        'Write out contents of self (intended for debugging)'
29
 
        Name_list = self.__dict__.keys()
30
 
        Name_list.sort()
31
 
        for Name in Name_list:
32
 
            print Name + ':', self.__dict__[Name]
33
 
 
34
 
    #--------------------------------------------------------
35
 
 
36
 
    def Write_data(self, Basename, Item_list):
37
 
        'Write the Lammps data to file (used by Write_Lammps)'
38
 
        import os, sys
39
 
 
40
 
        Filename = 'data.' + Basename
41
 
 
42
 
        Dir_list = os.listdir('.')
43
 
        i = 1
44
 
        while Filename in Dir_list:
45
 
            Filename = 'data' + `i` + '.' + Basename
46
 
            i = i +1
47
 
        del i
48
 
 
49
 
        print 'Writing', Filename + '...',
50
 
        sys.stdout.flush()
51
 
 
52
 
        try:
53
 
            F = open(Filename, 'w')
54
 
        except IOError, Detail:
55
 
            print '(error:', Detail[1] + '!)'
56
 
            return
57
 
 
58
 
        try:
59
 
            F.writelines(Item_list)
60
 
        except IOError, Detail:
61
 
            print '(error:', Detail[1] + '!)'
62
 
            F.close()
63
 
            return
64
 
 
65
 
        F.close()
66
 
        print 'done.'
67
 
 
68
 
    #--------------------------------------------------------
69
 
 
70
 
    def Write_Lammps(self, Basename):
71
 
        'Write the Lammps data file, ignoring blank sections'
72
 
        import string
73
 
        L = []
74
 
 
75
 
        L.append('LAMMPS data file for ' + self.name + '\n\n')
76
 
 
77
 
        L.append(`self.atoms` + ' atoms\n')
78
 
        L.append(`self.bonds` + ' bonds\n')
79
 
        L.append(`self.angles` + ' angles\n')
80
 
        L.append(`self.dihedrals` + ' dihedrals\n')
81
 
        L.append(`self.impropers` + ' impropers\n\n')
82
 
 
83
 
        L.append(`self.atom_types` + ' atom types\n')
84
 
        if self.bonds > 0:
85
 
            L.append(`self.bond_types` + ' bond types\n')
86
 
        if self.angles > 0:
87
 
            L.append(`self.angle_types` + ' angle types\n')
88
 
        if self.dihedrals > 0:
89
 
            L.append(`self.dihedral_types` + ' dihedral types\n')
90
 
        L.append('\n')
91
 
        
92
 
        L.append(`self.xlo` + ' ' + `self.xhi` + ' xlo xhi\n')
93
 
        L.append(`self.ylo` + ' ' + `self.yhi` + ' ylo yhi\n')
94
 
        L.append(`self.zlo` + ' ' + `self.zhi` + ' zlo zhi\n\n')
95
 
 
96
 
        if self.atom_types != 0:
97
 
            L.append('Masses\n\n')
98
 
            for i in range(self.atom_types):
99
 
                L.append(`i+1` + ' ' + `self.Masses[i]` + '\n')
100
 
            L.append('\n')
101
 
 
102
 
            L.append('Pair Coeffs\n\n')
103
 
            for i in range(self.atom_types):
104
 
                L.append(`i+1`)
105
 
                for j in range(len(self.Nonbond_Coeffs[0])):
106
 
                    L.append(' ' + `self.Nonbond_Coeffs[i][j]`)
107
 
                L.append('\n')
108
 
            L.append('\n')
109
 
 
110
 
        if self.bonds != 0 and self.bond_types != 0:
111
 
            L.append('Bond Coeffs\n\n')
112
 
            for i in range(self.bond_types):
113
 
                L.append(`i+1`)
114
 
                for j in range(len(self.Bond_Coeffs[0])):
115
 
                    L.append(' ' + `self.Bond_Coeffs[i][j]`)
116
 
                L.append('\n')
117
 
            L.append('\n')
118
 
 
119
 
        if self.angles != 0 and self.angle_types != 0:
120
 
            L.append('Angle Coeffs\n\n')
121
 
            for i in range(self.angle_types):
122
 
                L.append(`i+1`)
123
 
                for j in range(len(self.Angle_Coeffs[0])):
124
 
                    L.append(' ' + `self.Angle_Coeffs[i][j]`)
125
 
                L.append('\n')
126
 
            L.append('\n')
127
 
 
128
 
        if self.dihedrals != 0 and self.dihedral_types != 0:
129
 
            L.append('Dihedral Coeffs\n\n')
130
 
            for i in range(self.dihedral_types):
131
 
                L.append(`i+1`)
132
 
                for j in range(len(self.Dihedral_Coeffs[0])):
133
 
                    L.append(' ' + `self.Dihedral_Coeffs[i][j]`)
134
 
                L.append('\n')
135
 
            L.append('\n')
136
 
 
137
 
        if self.atoms != 0:
138
 
            L.append('Atoms\n\n')
139
 
            for i in range(self.atoms):
140
 
                L.append(`i+1`)
141
 
                for j in range(len(self.Atoms[0])):
142
 
                    L.append(' ' + `self.Atoms[i][j]`)
143
 
                L.append('\n')
144
 
            L.append('\n')
145
 
 
146
 
        if self.bonds != 0 and self.bond_types != 0:
147
 
            L.append('Bonds\n\n')
148
 
            for i in range(self.bonds):
149
 
                L.append(`i+1`)
150
 
                for j in range(len(self.Bonds[0])):
151
 
                    L.append(' ' + `self.Bonds[i][j]`)
152
 
                L.append('\n')
153
 
            L.append('\n')
154
 
 
155
 
        if self.angles != 0 and self.angle_types != 0:
156
 
            L.append('Angles\n\n')
157
 
            for i in range(self.angles):
158
 
                L.append(`i+1`)
159
 
                for j in range(len(self.Angles[0])):
160
 
                    L.append(' ' + `self.Angles[i][j]`)
161
 
                L.append('\n')
162
 
            L.append('\n')
163
 
 
164
 
        if self.dihedrals != 0 and self.dihedral_types != 0:
165
 
            L.append('Dihedrals\n\n')
166
 
            for i in range(self.dihedrals):
167
 
                L.append(`i+1`)
168
 
                for j in range(len(self.Dihedrals[0])):
169
 
                    L.append(' ' + `self.Dihedrals[i][j]`)
170
 
                L.append('\n')
171
 
            L.append('\n')
172
 
 
173
 
        self.Write_data(Basename, L)
174
 
 
175
 
#============================================================
176
 
 
177
 
class Amber:
178
 
    def __init__(self):
179
 
        'Initialise the Amber class'
180
 
        self.CRD_is_read = 0
181
 
        self.TOP_is_read = 0
182
 
 
183
 
    #--------------------------------------------------------
184
 
 
185
 
    def Dump(self):
186
 
        'Write out contents of self (intended for debugging)'
187
 
        Name_list = self.__dict__.keys()
188
 
        Name_list.sort()
189
 
        for Name in Name_list:
190
 
            print Name + ':', self.__dict__[Name]
191
 
 
192
 
    #--------------------------------------------------------
193
 
 
194
 
    def Coerce_to_Lammps(self):
195
 
        'Return the Amber data converted to Lammps format'
196
 
 
197
 
        import math
198
 
 
199
 
        if self.CRD_is_read and self.TOP_is_read:
200
 
            l = Lammps()
201
 
            print 'Converting...',
202
 
 
203
 
            l.name = self.ITITL
204
 
            l.atoms = self.NATOM
205
 
            l.bonds = self.NBONH + self.MBONA
206
 
            l.angles = self.NTHETH + self.MTHETA
207
 
            l.dihedrals = self.NPHIH + self.MPHIA
208
 
            l.impropers = 0
209
 
            l.atom_types = self.NTYPES
210
 
            l.bond_types = self.NUMBND
211
 
            l.angle_types = self.NUMANG
212
 
            l.dihedral_types = self.NPTRA
213
 
 
214
 
            Shift = 0
215
 
            if self.__dict__.has_key('BOX'):
216
 
                l.xlo = 0.0
217
 
                l.xhi = self.BOX[0]
218
 
                l.ylo = 0.0
219
 
                l.yhi = self.BOX[1]
220
 
                l.zlo = 0.0
221
 
                l.zhi = self.BOX[2]
222
 
                if (l.xlo > min(self.X)) or (l.xhi < max(self.X)) or \
223
 
                   (l.ylo > min(self.Y)) or (l.yhi < max(self.Y)) or \
224
 
                   (l.zlo > min(self.Z)) or (l.zhi < max(self.Z)):
225
 
                    #  Vikas Modification: Disabling Shifting. This means I am intend to send exact coordinates of each atom and let LAMMPS 
226
 
                    #  take care of imaging into periodic image cells. If one wants to shift all atoms in the periodic box, 
227
 
                    #  please uncomment the below 2 lines.
228
 
                    print '(warning: Currently not shifting the atoms to the periodic box)' 
229
 
                    #Shift = 1
230
 
            else:
231
 
                print '(warning: Guessing at periodic box!)',
232
 
                l.xlo = min(self.X)
233
 
                l.xhi = max(self.X)
234
 
                l.ylo = min(self.Y)
235
 
                l.yhi = max(self.Y)
236
 
                l.zlo = min(self.Z)
237
 
                l.zhi = max(self.Z)
238
 
 
239
 
            # This doesn't check duplicate values
240
 
            l.Masses = []
241
 
            for i in range(l.atom_types):
242
 
                l.Masses.append(0)
243
 
            for i in range(self.NATOM):
244
 
                l.Masses[self.IAC[i] - 1] = self.AMASS[i]
245
 
 
246
 
            l.Nonbond_Coeffs = []
247
 
            for i in range(self.NTYPES):
248
 
                l.Nonbond_Coeffs.append([0,0])
249
 
            for i in range(self.NTYPES):
250
 
                j = self.ICO[i * (self.NTYPES + 1)] - 1
251
 
                if self.CN1[j] == 0.0:
252
 
                    l.Nonbond_Coeffs[i][0] = 0.0
253
 
                else:
254
 
                    l.Nonbond_Coeffs[i][0] = \
255
 
                                0.25 * (self.CN2[j])**2 / self.CN1[j]
256
 
                if self.CN2[j] == 0.0:
257
 
                    l.Nonbond_Coeffs[i][1] = 0.0
258
 
                else:
259
 
                    l.Nonbond_Coeffs[i][1] = \
260
 
                                (self.CN1[j] / self.CN2[j])**(1.0/6.0)
261
 
 
262
 
            l.Bond_Coeffs = []
263
 
            for i in range(self.NUMBND):
264
 
                l.Bond_Coeffs.append([0,0])
265
 
            for i in range(self.NUMBND):
266
 
                l.Bond_Coeffs[i][0] = self.RK[i]
267
 
                l.Bond_Coeffs[i][1] = self.REQ[i]
268
 
 
269
 
            l.Angle_Coeffs = []
270
 
            for i in range(self.NUMANG):
271
 
                l.Angle_Coeffs.append([0,0])
272
 
            for i in range(self.NUMANG):
273
 
                l.Angle_Coeffs[i][0] = self.TK[i]
274
 
                l.Angle_Coeffs[i][1] = (180/math.pi) * self.TEQ[i]
275
 
 
276
 
            l.Dihedral_Coeffs = []
277
 
            for i in range(self.NPTRA):
278
 
                l.Dihedral_Coeffs.append([0,0,0])
279
 
            for i in range(self.NPTRA):
280
 
                l.Dihedral_Coeffs[i][0] = self.PK[i]
281
 
                if self.PHASE[i] == 0:
282
 
                    l.Dihedral_Coeffs[i][1] = 1
283
 
                else:
284
 
                    l.Dihedral_Coeffs[i][1] = -1
285
 
                l.Dihedral_Coeffs[i][2] = int(self.PN[i])
286
 
 
287
 
            l.Atoms = []
288
 
            for i in range(self.NATOM):
289
 
                x = self.X[i]
290
 
                y = self.Y[i]
291
 
                z = self.Z[i]
292
 
                if Shift:
293
 
                    while x < l.xlo:
294
 
                        x = x + self.BOX[0]
295
 
                    while x > l.xhi:
296
 
                        x = x - self.BOX[0]
297
 
                    while y < l.ylo:
298
 
                        y = y + self.BOX[1]
299
 
                    while y > l.yhi:
300
 
                        y = y - self.BOX[1]
301
 
                    while z < l.zlo:
302
 
                        z = z + self.BOX[2]
303
 
                    while z > l.zhi:
304
 
                        z = z - self.BOX[2]
305
 
                l.Atoms.append([0, self.IAC[i], self.CHRG[i]/18.2223, \
306
 
                                x, y, z])
307
 
 
308
 
            l.Bonds = []
309
 
            for i in range(l.bonds):
310
 
                l.Bonds.append([0,0,0])
311
 
            for i in range(self.NBONH):
312
 
                l.Bonds[i][0] = self.ICBH[i]
313
 
                l.Bonds[i][1] = abs(self.IBH[i])/3 + 1
314
 
                l.Bonds[i][2] = abs(self.JBH[i])/3 + 1
315
 
            for i in range(self.NBONA):
316
 
                l.Bonds[self.NBONH + i][0] = self.ICB[i]
317
 
                l.Bonds[self.NBONH + i][1] = abs(self.IB[i])/3 + 1
318
 
                l.Bonds[self.NBONH + i][2] = abs(self.JB[i])/3 + 1
319
 
 
320
 
            l.Angles = []
321
 
            for i in range(l.angles):
322
 
                l.Angles.append([0,0,0,0])
323
 
            for i in range(self.NTHETH):
324
 
                l.Angles[i][0] = self.ICTH[i]
325
 
                l.Angles[i][1] = abs(self.ITH[i])/3 + 1
326
 
                l.Angles[i][2] = abs(self.JTH[i])/3 + 1
327
 
                l.Angles[i][3] = abs(self.KTH[i])/3 + 1
328
 
            for i in range(self.NTHETA):
329
 
                l.Angles[self.NTHETH + i][0] = self.ICT[i]
330
 
                l.Angles[self.NTHETH + i][1] = abs(self.IT[i])/3 + 1
331
 
                l.Angles[self.NTHETH + i][2] = abs(self.JT[i])/3 + 1
332
 
                l.Angles[self.NTHETH + i][3] = abs(self.KT[i])/3 + 1
333
 
 
334
 
            l.Dihedrals = []
335
 
            for i in range(l.dihedrals):
336
 
                l.Dihedrals.append([0,0,0,0,0])
337
 
            for i in range(self.NPHIH):
338
 
                l.Dihedrals[i][0] = self.ICPH[i]
339
 
                l.Dihedrals[i][1] = abs(self.IPH[i])/3 + 1
340
 
                l.Dihedrals[i][2] = abs(self.JPH[i])/3 + 1
341
 
                l.Dihedrals[i][3] = abs(self.KPH[i])/3 + 1
342
 
                l.Dihedrals[i][4] = abs(self.LPH[i])/3 + 1
343
 
            for i in range(self.NPHIA):
344
 
                l.Dihedrals[self.NPHIH + i][0] = self.ICP[i]
345
 
                l.Dihedrals[self.NPHIH + i][1] = abs(self.IP[i])/3 + 1
346
 
                l.Dihedrals[self.NPHIH + i][2] = abs(self.JP[i])/3 + 1
347
 
                l.Dihedrals[self.NPHIH + i][3] = abs(self.KP[i])/3 + 1
348
 
                l.Dihedrals[self.NPHIH + i][4] = abs(self.LP[i])/3 + 1
349
 
 
350
 
            print 'done.'
351
 
            return l
352
 
        else:
353
 
            print '(Error: Not all the Amber data has been read!)'
354
 
    
355
 
    #--------------------------------------------------------
356
 
 
357
 
    def Read_data(self, Filename):
358
 
        'Read the filename, returning a list of strings'
359
 
 
360
 
        import string, sys
361
 
 
362
 
        print 'Reading', Filename + '...',
363
 
        sys.stdout.flush()
364
 
 
365
 
        try:
366
 
            F = open(Filename)
367
 
        except IOError, Detail:
368
 
            print '(error:', Detail[1] + '!)'
369
 
            return
370
 
 
371
 
        try:
372
 
            Lines = F.readlines()
373
 
        except IOError, Detail:
374
 
            print '(error:', Detail[1] + '!)'
375
 
            F.close()
376
 
            return
377
 
 
378
 
        F.close()
379
 
 
380
 
        # If the first line is empty, use the Basename
381
 
        if Filename[-4:] == '.crd':
382
 
            if string.split(Lines[0]) == []: # This line corresponds to TITLE name in CRD file
383
 
                Basename = Filename[:string.find(Filename, '.')]
384
 
                Item_list = [Basename]
385
 
                print 'Warning: Title not present... Assigning Basename as Title'
386
 
            else:
387
 
                Item_list = []
388
 
        else:
389
 
            if string.split(Lines[3]) == []: # This line corresponds to TITLE name in TOPOLOGY file
390
 
                Basename = Filename[:string.find(Filename, '.')]
391
 
                Item_list = [Basename]
392
 
                print 'Warning: Title not present... Assigning Basename as Title'
393
 
            else:
394
 
                Item_list = []
395
 
 
396
 
        for Line in Lines:
397
 
            if Line[0]!='%': #Vikas' Modification: This condition ignores all the lines starting with % in the topology file.
398
 
                Item_list.extend(string.split(Line))
399
 
 
400
 
        return Item_list
401
 
 
402
 
    #--------------------------------------------------------
403
 
 
404
 
    def Read_CRD(self, Basename):
405
 
        'Read the Amber coordinate/restart (.crd) file'
406
 
 
407
 
        # The optional velocities and periodic box size are not yet parsed.
408
 
 
409
 
        Item_list = self.Read_data(Basename + '.crd')
410
 
 
411
 
        if Item_list == None:
412
 
            return
413
 
        elif len(Item_list) < 2:
414
 
            print '(error: File too short!)'
415
 
            return
416
 
 
417
 
        # Parse the data
418
 
        if self.__dict__.has_key('ITITL'):
419
 
            if Pop(Item_list,0) != self.ITITL:
420
 
                print '(warning: ITITL differs!)',
421
 
        else:
422
 
            self.ITITL = Pop(Item_list,0)
423
 
        print self.ITITL #Vikas Modification : Priting the Title
424
 
 
425
 
        if self.__dict__.has_key('NATOM'):
426
 
            if eval(Pop(Item_list,0)) != self.NATOM:
427
 
                print '(error: NATOM differs!)'
428
 
                return
429
 
        else:
430
 
            self.NATOM = eval(Pop(Item_list,0))
431
 
        print self.NATOM # Vikas' Modification: Printing number of atoms just to make sure that the program is reading the correct value. 
432
 
 
433
 
        #if len(Item_list) == 1 + 3 * self.NATOM:
434
 
        # Vikas' Modification: I changed the condition.
435
 
        if (len(Item_list)%3) != 0:
436
 
            self.TIME = eval(Pop(Item_list,0))
437
 
        else:
438
 
            self.TIME = 0
439
 
        print self.TIME # Vikas' Modification : Printing simulation time, just to make sure that the program is readint the correct value.
440
 
        if len(Item_list) < 3 * self.NATOM:
441
 
            print '(error: File too short!)'
442
 
            return
443
 
 
444
 
        self.X = []
445
 
        self.Y = []
446
 
        self.Z = []
447
 
        for i in range(self.NATOM):
448
 
            self.X.append(eval(Pop(Item_list,0)))
449
 
            self.Y.append(eval(Pop(Item_list,0)))
450
 
            self.Z.append(eval(Pop(Item_list,0)))
451
 
 
452
 
        if (self.NATOM == 1) and len(Item_list):
453
 
            print '(warning: Ambiguity!)',
454
 
 
455
 
        if len(Item_list) >= 3 * self.NATOM:
456
 
                self.VX = []
457
 
                self.VY = []
458
 
                self.VZ = []
459
 
                for i in range(self.NATOM):
460
 
                    self.VX.append(eval(Pop(Item_list,0)))
461
 
                    self.VY.append(eval(Pop(Item_list,0)))
462
 
                    self.VZ.append(eval(Pop(Item_list,0)))
463
 
 
464
 
        if len(Item_list) >= 3:
465
 
            self.BOX = []
466
 
            for i in range(3):
467
 
                self.BOX.append(eval(Pop(Item_list,0)))
468
 
 
469
 
        if len(Item_list):
470
 
            print '(warning: File too large!)',
471
 
        
472
 
        print 'done.'
473
 
        self.CRD_is_read = 1
474
 
 
475
 
    #--------------------------------------------------------
476
 
 
477
 
    def Read_TOP(self, Basename):
478
 
        'Read the Amber parameter/topology (.top) file'
479
 
        Item_list = self.Read_data(Basename + '.top')
480
 
 
481
 
        if Item_list == None:
482
 
            return
483
 
        elif len(Item_list) < 31:
484
 
            print '(error: File too short!)'
485
 
            return
486
 
 
487
 
       # Parse the data
488
 
        if self.__dict__.has_key('ITITL'):
489
 
            if Pop(Item_list,0) != self.ITITL:
490
 
                print '(warning: ITITL differs!)'
491
 
        else:
492
 
            self.ITITL = Pop(Item_list,0)
493
 
        print self.ITITL # Printing Self Title
494
 
 
495
 
        if self.__dict__.has_key('NATOM'):
496
 
            if eval(Pop(Item_list,0)) != self.NATOM:
497
 
                print '(error: NATOM differs!)'
498
 
                return
499
 
        else:
500
 
            self.NATOM = eval(Pop(Item_list,0))
501
 
        print self.NATOM # Printing total number of atoms just to make sure that thing are going right
502
 
        self.NTYPES = eval(Pop(Item_list,0))
503
 
        self.NBONH  = eval(Pop(Item_list,0))
504
 
        self.MBONA  = eval(Pop(Item_list,0))
505
 
        self.NTHETH = eval(Pop(Item_list,0))
506
 
        self.MTHETA = eval(Pop(Item_list,0))
507
 
        self.NPHIH  = eval(Pop(Item_list,0))
508
 
        self.MPHIA  = eval(Pop(Item_list,0))
509
 
        self.NHPARM = eval(Pop(Item_list,0))
510
 
        self.NPARM  = eval(Pop(Item_list,0))
511
 
        self.NEXT   = eval(Pop(Item_list,0))
512
 
        self.NRES   = eval(Pop(Item_list,0))
513
 
        self.NBONA  = eval(Pop(Item_list,0))
514
 
        self.NTHETA = eval(Pop(Item_list,0))
515
 
        self.NPHIA  = eval(Pop(Item_list,0))
516
 
        self.NUMBND = eval(Pop(Item_list,0))
517
 
        self.NUMANG = eval(Pop(Item_list,0))
518
 
        self.NPTRA  = eval(Pop(Item_list,0))
519
 
        self.NATYP  = eval(Pop(Item_list,0))
520
 
        self.NPHB   = eval(Pop(Item_list,0))
521
 
        self.IFPERT = eval(Pop(Item_list,0))
522
 
        self.NBPER  = eval(Pop(Item_list,0))
523
 
        self.NGPER  = eval(Pop(Item_list,0))
524
 
        self.NDPER  = eval(Pop(Item_list,0))
525
 
        self.MBPER  = eval(Pop(Item_list,0))
526
 
        self.MGPER  = eval(Pop(Item_list,0))
527
 
        self.MDPER  = eval(Pop(Item_list,0))
528
 
        self.IFBOX  = eval(Pop(Item_list,0))
529
 
        self.NMXRS  = eval(Pop(Item_list,0))
530
 
        self.IFCAP  = eval(Pop(Item_list,0))
531
 
 
532
 
        #....................................................
533
 
 
534
 
        if len(Item_list) < 5 * self.NATOM + self.NTYPES**2 + \
535
 
           2*(self.NRES + self.NUMBND + self.NUMANG) + \
536
 
           3*self.NPTRA + self.NATYP:
537
 
            print '(error: File too short!)'
538
 
            return -1
539
 
 
540
 
        self.IGRAPH = []
541
 
        Pop(Item_list,0)
542
 
 
543
 
          # A little kludge is needed here, since the IGRAPH strings are
544
 
        # not separated by spaces if 4 characters in length.
545
 
        for i in range(self.NATOM):
546
 
            if len(Item_list[0]) > 4:
547
 
                Item_list.insert(1, Item_list[0][4:])
548
 
                Item_list.insert(1, Item_list[0][0:4])
549
 
                del Item_list[0]
550
 
            self.IGRAPH.append(Pop(Item_list,0))
551
 
        
552
 
        # Vikas' Modification : In the following section, I am printing out each quantity which is currently being read from the topology file.
553
 
        print 'Reading Charges...'
554
 
        self.CHRG = []
555
 
        for i in range(self.NATOM):
556
 
            self.CHRG.append(eval(Pop(Item_list,0)))
557
 
 
558
 
        print 'Reading Atomic Number...'
559
 
        self.ANUMBER = []
560
 
        for i in range(self.NATOM):
561
 
            self.ANUMBER.append(eval(Pop(Item_list,0)))
562
 
 
563
 
        print 'Reading Atomic Masses...'
564
 
        self.AMASS = []
565
 
        for i in range(self.NATOM):
566
 
            self.AMASS.append(eval(Pop(Item_list,0)))
567
 
 
568
 
        print 'Reading Atom Types...'
569
 
        self.IAC = []
570
 
        for i in range(self.NATOM):
571
 
            self.IAC.append(eval(Pop(Item_list,0)))
572
 
 
573
 
        print 'Reading Excluded Atoms...'
574
 
        self.NUMEX = []
575
 
        for i in range(self.NATOM):
576
 
            self.NUMEX.append(eval(Pop(Item_list,0)))
577
 
 
578
 
        print 'Reading Non-bonded Parameter Index...'
579
 
        self.ICO = []
580
 
        for i in range(self.NTYPES**2):
581
 
            self.ICO.append(eval(Pop(Item_list,0)))
582
 
 
583
 
        print 'Reading Residue Labels...'
584
 
        self.LABRES = []
585
 
        for i in range(self.NRES):
586
 
            self.LABRES.append(Pop(Item_list,0))
587
 
 
588
 
        print 'Reading Residues Starting Pointers...'
589
 
        self.IPRES = []
590
 
        for i in range(self.NRES):
591
 
            self.IPRES.append(eval(Pop(Item_list,0)))
592
 
 
593
 
        print 'Reading Bond Force Constants...'
594
 
        self.RK = []
595
 
        for i in range(self.NUMBND):
596
 
            self.RK.append(eval(Pop(Item_list,0)))
597
 
 
598
 
        print 'Reading Equilibrium Bond Values...'
599
 
        self.REQ = []
600
 
        for i in range(self.NUMBND):
601
 
            self.REQ.append(eval(Pop(Item_list,0)))
602
 
 
603
 
        print 'Reading Angle Force Constants...' 
604
 
        self.TK = []
605
 
        for i in range(self.NUMANG):
606
 
            self.TK.append(eval(Pop(Item_list,0)))
607
 
 
608
 
        print 'Reading Equilibrium Angle Values...'
609
 
        self.TEQ = []
610
 
        for i in range(self.NUMANG):
611
 
            self.TEQ.append(eval(Pop(Item_list,0)))
612
 
 
613
 
        print 'Reading Dihedral Force Constants...'
614
 
        self.PK = []
615
 
        for i in range(self.NPTRA):
616
 
            self.PK.append(eval(Pop(Item_list,0)))
617
 
        
618
 
        print 'Reading Dihedral Periodicity...' 
619
 
        self.PN = []
620
 
        for i in range(self.NPTRA):
621
 
            self.PN.append(eval(Pop(Item_list,0)))
622
 
        
623
 
        print 'Reading Dihedral Phase...'
624
 
        self.PHASE = []
625
 
        for i in range(self.NPTRA):
626
 
            self.PHASE.append(eval(Pop(Item_list,0)))
627
 
        
628
 
        print 'Reading 1-4 Electrostatic Scaling Factor...'
629
 
        self.SCEEFAC = []
630
 
        for i in range(self.NPTRA):
631
 
            self.SCEEFAC.append(eval(Pop(Item_list,0)))
632
 
        
633
 
        print 'Reading 1-4 Van der Waals Scaling Factor...'
634
 
        self.SCNBFAC = []
635
 
        for i in range(self.NPTRA):
636
 
            self.SCNBFAC.append(eval(Pop(Item_list,0)))
637
 
        
638
 
        print 'Reading Solty...' #I think this is currently not used in AMBER. Check it out, though
639
 
        self.SOLTY = []
640
 
        for i in range(self.NATYP):
641
 
            self.SOLTY.append(eval(Pop(Item_list,0)))
642
 
 
643
 
        #....................................................
644
 
 
645
 
        if len(Item_list) < 2 * self.NTYPES * (self.NTYPES + 1) / 2:
646
 
            print '(error: File too short!)'
647
 
            return -1
648
 
 
649
 
        print 'Reading LJ A Coefficient...'
650
 
        self.CN1 = []
651
 
        for i in range(self.NTYPES * (self.NTYPES + 1) / 2):
652
 
            self.CN1.append(eval(Pop(Item_list,0)))
653
 
 
654
 
        print 'Reading LJ B Coefficient...'
655
 
        self.CN2 = []
656
 
        for i in range(self.NTYPES * (self.NTYPES + 1) / 2):
657
 
            self.CN2.append(eval(Pop(Item_list,0)))
658
 
 
659
 
        #....................................................
660
 
 
661
 
        if len(Item_list) < 3 * (self.NBONH + self.NBONA) + \
662
 
           4 * (self.NTHETH + self.NTHETA) + 5 * (self.NPHIH + self.NPHIA):
663
 
            print '(error: File too short!)'
664
 
            return -1
665
 
 
666
 
        print 'Reading Bonds which include hydrogen...'
667
 
        self.IBH = []
668
 
        self.JBH = []
669
 
        self.ICBH = []
670
 
        for i in range(self.NBONH):
671
 
            self.IBH.append(eval(Pop(Item_list,0)))
672
 
            self.JBH.append(eval(Pop(Item_list,0)))
673
 
            self.ICBH.append(eval(Pop(Item_list,0)))
674
 
 
675
 
        print 'Reading Bonds which dont include hydrogen...'
676
 
        self.IB = []
677
 
        self.JB = []
678
 
        self.ICB = []
679
 
        for i in range(self.NBONA):
680
 
            self.IB.append(eval(Pop(Item_list,0)))
681
 
            self.JB.append(eval(Pop(Item_list,0)))
682
 
            self.ICB.append(eval(Pop(Item_list,0)))
683
 
 
684
 
        print 'Reading Angles which include hydrogen...' 
685
 
        self.ITH = []
686
 
        self.JTH = []
687
 
        self.KTH = []
688
 
        self.ICTH = []
689
 
        for i in range(self.NTHETH):
690
 
            self.ITH.append(eval(Pop(Item_list,0)))
691
 
            self.JTH.append(eval(Pop(Item_list,0)))
692
 
            self.KTH.append(eval(Pop(Item_list,0)))
693
 
            self.ICTH.append(eval(Pop(Item_list,0)))
694
 
 
695
 
        print 'Reading Angles which dont include hydrogen...'
696
 
        self.IT = []
697
 
        self.JT = []
698
 
        self.KT = []
699
 
        self.ICT = []
700
 
        for i in range(self.NTHETA):
701
 
            self.IT.append(eval(Pop(Item_list,0)))
702
 
            self.JT.append(eval(Pop(Item_list,0)))
703
 
            self.KT.append(eval(Pop(Item_list,0)))
704
 
            self.ICT.append(eval(Pop(Item_list,0)))
705
 
 
706
 
        print 'Reading Dihedrals which include hydrogen...'
707
 
        self.IPH = []
708
 
        self.JPH = []
709
 
        self.KPH = []
710
 
        self.LPH = []
711
 
        self.ICPH = []
712
 
        for i in range(self.NPHIH):
713
 
            self.IPH.append(eval(Pop(Item_list,0)))
714
 
            self.JPH.append(eval(Pop(Item_list,0)))
715
 
            self.KPH.append(eval(Pop(Item_list,0)))
716
 
            self.LPH.append(eval(Pop(Item_list,0)))
717
 
            self.ICPH.append(eval(Pop(Item_list,0)))
718
 
 
719
 
        print 'Reading Dihedrals which dont include hydrogen...'
720
 
        self.IP = []
721
 
        self.JP = []
722
 
        self.KP = []
723
 
        self.LP = []
724
 
        self.ICP = []
725
 
        for i in range(self.NPHIA):
726
 
            self.IP.append(eval(Pop(Item_list,0)))
727
 
            self.JP.append(eval(Pop(Item_list,0)))
728
 
            self.KP.append(eval(Pop(Item_list,0)))
729
 
            self.LP.append(eval(Pop(Item_list,0)))
730
 
            self.ICP.append(eval(Pop(Item_list,0)))
731
 
 
732
 
        #....................................................
733
 
 
734
 
        if len(Item_list) < self.NEXT + 3 * self.NPHB + 4 * self.NATOM:
735
 
            print '(error: File too short!)'
736
 
            return -1
737
 
 
738
 
        print 'Reading Excluded Atom List...'
739
 
        self.NATEX = []
740
 
        for i in range(self.NEXT):
741
 
            self.NATEX.append(eval(Pop(Item_list,0)))
742
 
 
743
 
        print 'Reading H-Bond A Coefficient, corresponding to r**12 term for all possible types...'
744
 
        self.ASOL = []
745
 
        for i in range(self.NPHB):
746
 
            self.ASOL.append(eval(Pop(Item_list,0)))
747
 
 
748
 
        print 'Reading H-Bond B Coefficient, corresponding to r**10 term for all possible types...'
749
 
        self.BSOL = []
750
 
        for i in range(self.NPHB):
751
 
            self.BSOL.append(eval(Pop(Item_list,0)))
752
 
 
753
 
        print 'Reading H-Bond Cut...' # I think it is not being used nowadays
754
 
        self.HBCUT = []
755
 
        for i in range(self.NPHB):
756
 
            self.HBCUT.append(eval(Pop(Item_list,0)))
757
 
 
758
 
        print 'Reading Amber Atom Types for each atom...'
759
 
        self.ISYMBL = []
760
 
        for i in range(self.NATOM):
761
 
            self.ISYMBL.append(Pop(Item_list,0))
762
 
 
763
 
        print 'Reading Tree Chain Classification...'
764
 
        self.ITREE = []
765
 
        for i in range(self.NATOM):
766
 
            self.ITREE.append(Pop(Item_list,0))
767
 
 
768
 
        print 'Reading Join Array: Tree joining information' # Currently unused in Sander, an AMBER module
769
 
        self.JOIN = []
770
 
        for i in range(self.NATOM):
771
 
            self.JOIN.append(eval(Pop(Item_list,0)))
772
 
 
773
 
        print 'Reading IRotate...' # Currently unused in Sander and Gibbs
774
 
        self.IROTAT = []
775
 
        for i in range(self.NATOM):
776
 
            self.IROTAT.append(eval(Pop(Item_list,0)))
777
 
 
778
 
        #....................................................
779
 
 
780
 
        if self.IFBOX > 0:
781
 
            if len(Item_list) < 3:
782
 
                print '(error: File too short!)'
783
 
                return -1
784
 
 
785
 
            print 'Reading final residue which is part of solute...'
786
 
            self.IPTRES = eval(Pop(Item_list,0))
787
 
            print 'Reading total number of molecules...'
788
 
            self.NSPM = eval(Pop(Item_list,0))
789
 
            print 'Reading first solvent moleule index...'
790
 
            self.NSPSOL = eval(Pop(Item_list,0))
791
 
 
792
 
            if len(Item_list) < self.NSPM + 4:
793
 
                print '(error: File too short!)'
794
 
                return -1
795
 
 
796
 
            print 'Reading atom per molecule...' 
797
 
            self.NSP = []
798
 
            for i in range(self.NSPM):
799
 
                self.NSP.append(eval(Pop(Item_list,0)))
800
 
 
801
 
            self.BETA = eval(Pop(Item_list,0))
802
 
 
803
 
            print 'Reading Box Dimensions...' 
804
 
            if self.__dict__.has_key('BOX'):
805
 
                BOX = []
806
 
                for i in range(3):
807
 
                    BOX.append(eval(Pop(Item_list,0)))
808
 
                for i in range(3):
809
 
                    if BOX[i] != self.BOX[i]:
810
 
                        print '(warning: BOX differs!)',
811
 
                        break
812
 
                del BOX
813
 
            else:
814
 
                self.BOX = []
815
 
                for i in range(3):
816
 
                    self.BOX.append(eval(Pop(Item_list,0)))
817
 
 
818
 
        #....................................................
819
 
 
820
 
        if self.IFCAP > 0:
821
 
            if len(Item_list) < 5:
822
 
                print '(error: File too short!)'
823
 
                return -1
824
 
            print 'Reading ICAP variables::: For details, refer to online AMBER format manual'
825
 
            self.NATCAP = eval(Pop(Item_list,0))
826
 
            self.CUTCAP = eval(Pop(Item_list,0))
827
 
            self.XCAP = eval(Pop(Item_list,0))
828
 
            self.YCAP = eval(Pop(Item_list,0))
829
 
            self.ZCAP = eval(Pop(Item_list,0))
830
 
 
831
 
        #....................................................
832
 
 
833
 
        if self.IFPERT > 0:
834
 
            if len(Item_list) < 4 * self.NBPER + 5 * self.NGPER + \
835
 
               6 * self.NDPER + self.NRES + 6 * self.NATOM:
836
 
                print '(error: File too short!)'
837
 
                return -1
838
 
 
839
 
            print 'Reading perturb variables, 1. Bond, 2. Angles, 3. Dihedrals, etc etc.::: For details, refer to online AMBER format manual' 
840
 
            self.IBPER = []
841
 
            self.JBPER = []
842
 
            for i in range(self.NBPER):
843
 
                self.IBPER.append(eval(Pop(Item_list,0)))
844
 
                self.JBPER.append(eval(Pop(Item_list,0)))
845
 
 
846
 
            self.ICBPER = []
847
 
            for i in range(2 * self.NBPER):
848
 
                self.ICBPER.append(eval(Pop(Item_list,0)))
849
 
 
850
 
            self.ITPER = []
851
 
            self.JTPER = []
852
 
            self.KTPER = []
853
 
            for i in range(self.NGPER):
854
 
                self.ITPER.append(eval(Pop(Item_list,0)))
855
 
                self.JTPER.append(eval(Pop(Item_list,0)))
856
 
                self.KTPER.append(eval(Pop(Item_list,0)))
857
 
 
858
 
            self.ICTPER = []
859
 
            for i in range(2 * self.NGPER):
860
 
                self.ICTPER.append(eval(Pop(Item_list,0)))
861
 
 
862
 
            self.IPPER = []
863
 
            self.JPPER = []
864
 
            self.KPPER = []
865
 
            self.LPPER = []
866
 
            for i in range(self.NDPER):
867
 
                self.IPPER.append(eval(Pop(Item_list,0)))
868
 
                self.JPPER.append(eval(Pop(Item_list,0)))
869
 
                self.KPPER.append(eval(Pop(Item_list,0)))
870
 
                self.LPPER.append(eval(Pop(Item_list,0)))
871
 
 
872
 
            self.ICPPER = []
873
 
            for i in range(2 * self.NDPER):
874
 
                self.ICPPER.append(eval(Pop(Item_list,0)))
875
 
 
876
 
            LABRES = []
877
 
            for i in range(self.NRES):
878
 
                LABRES.append(Pop(Item_list,0))
879
 
            for i in range(self.NRES):
880
 
                if LABRES[i] != self.LABRES[i]:
881
 
                    print '(warning: BOX differs!)',
882
 
                    break
883
 
 
884
 
            self.IGRPER = []
885
 
            for i in range(self.NATOM):
886
 
                self.IGRPER.append(eval(Pop(Item_list,0)))
887
 
 
888
 
            self.ISMPER = []
889
 
            for i in range(self.NATOM):
890
 
                self.ISMPER.append(eval(Pop(Item_list,0)))
891
 
 
892
 
            self.ALMPER = []
893
 
            for i in range(self.NATOM):
894
 
                self.ALMPER.append(eval(Pop(Item_list,0)))
895
 
 
896
 
            self.IAPER = []
897
 
            for i in range(self.NATOM):
898
 
                self.IAPER.append(eval(Pop(Item_list,0)))
899
 
 
900
 
            self.IACPER = []
901
 
            for i in range(self.NATOM):
902
 
                self.IACPER.append(eval(Pop(Item_list,0)))
903
 
 
904
 
            self.CGPER = []
905
 
            for i in range(self.NATOM):
906
 
                self.CGPER.append(eval(Pop(Item_list,0)))
907
 
 
908
 
        #....................................................
909
 
 
910
 
        self.IPOL = 0
911
 
        if self.IPOL == 1:
912
 
            if len(Item_list) < self.NATOM:
913
 
                print '(error: File too short!)'
914
 
                return -1
915
 
            print 'Reading Polarizability Data. For details, refer to online AMBER format manual'
916
 
            self.ATPOL = []
917
 
            for i in range(self.NATOM):
918
 
                self.ATPOL.append(eval(Pop(Item_list,0)))
919
 
 
920
 
            if self.IFPERT == 1:
921
 
                if len(Item_list) < self.NATOM:
922
 
                    print '(error: File too short!)'
923
 
                    return -1
924
 
                self.ATPOL1 = []
925
 
                for i in range(self.NATOM):
926
 
                    self.ATPOL1.append(eval(Pop(Item_list,0)))
927
 
 
928
 
        #....................................................
929
 
 
930
 
        if len(Item_list):
931
 
            print '(warning: File too large!)',
932
 
 
933
 
        print 'done.'
934
 
        self.TOP_is_read = 1
935
 
 
936
 
#============================================================
937
 
 
938
 
def Find_Amber_files():
939
 
    'Look for sets of Amber files to process'
940
 
    '''If not passed anything on the command line, look for pairs of
941
 
    Amber files (.crd and .top) in the current directory.  For
942
 
    each set if there is no corresponding Lammps file (data.), or it is
943
 
    older than any of the Amber files, add its basename to a list of
944
 
    strings.  This list is returned by the function'''
945
 
 
946
 
    # Date and existence checks not yet implemented
947
 
 
948
 
    import os, sys
949
 
 
950
 
    Basename_list = []
951
 
 
952
 
    # Extract basenames from command line
953
 
    for Name in sys.argv[1:]:
954
 
        if Name[-4:] == '.crd':
955
 
            Basename_list.append(Name[:-4])
956
 
        else:
957
 
            if Name[-4:] == '.top':
958
 
                Basename_list.append(Name[:-4])
959
 
            else:
960
 
                Basename_list.append(Name)
961
 
 
962
 
    # Remove duplicate basenames
963
 
    for Basename in Basename_list[:]:
964
 
        while Basename_list.count(Basename) > 1:
965
 
            Basename_list.remove(Basename)
966
 
 
967
 
    if Basename_list == []:
968
 
        print 'Looking for Amber files...',
969
 
        Dir_list = os.listdir('.')
970
 
        Dir_list.sort()
971
 
        for File in Dir_list:
972
 
            if File[-4:] == '.top':
973
 
                Basename = File[:-4]
974
 
                if (Basename + '.crd') in Dir_list:
975
 
                    Basename_list.append(Basename)
976
 
        if Basename_list != []:
977
 
            print 'found',
978
 
            for i in range(len(Basename_list)-1):
979
 
                print Basename_list[i] + ',',
980
 
            print Basename_list[-1] + '\n'
981
 
    
982
 
    if Basename_list == []:
983
 
        print 'none.\n'
984
 
 
985
 
    return Basename_list
986
 
    
987
 
#============================================================
988
 
 
989
 
def Convert_Amber_files():
990
 
    'Handle the whole conversion process'
991
 
    print
992
 
    print 'Welcome to amber2lammps, a program to convert Amber files to Lammps format!'
993
 
    print
994
 
    Basename_list = Find_Amber_files()
995
 
    for Basename in Basename_list:
996
 
        a = Amber()
997
 
        a.Read_CRD(Basename)
998
 
        if a.CRD_is_read:
999
 
            a.Read_TOP(Basename)
1000
 
            if a.TOP_is_read:
1001
 
                l = a.Coerce_to_Lammps()
1002
 
                l.Write_Lammps(Basename)
1003
 
                del l
1004
 
        del a
1005
 
        print
1006
 
 
1007
 
#============================================================
1008
 
 
1009
 
Convert_Amber_files()
 
1
#! /usr/bin/evn python2
 
2
 
 
3
# This is amber2lammps, a program written by Keir E. Novik to convert
 
4
# Amber files to Lammps files.
 
5
#
 
6
# Copyright 1999, 2000 Keir E. Novik; all rights reserved.
 
7
# 
 
8
# Modified by Vikas Varshney, U Akron, 5 July 2005, as described in README
 
9
# Bug Fixed :Third argument in Dihedral Coeffs section is an integer - Ketan S Khare September 26, 2011
 
10
# Modified by Vikas Varshney, Oct 8, 2013 to include additional flags (Atomic_Number, Coulombic and van der Waals 1-4 factors which are included in newer vesions of .top and .crd files in amber12. 
 
11
 
 
12
#============================================================
 
13
 
 
14
def Pop(S, I=-1):
 
15
    'Pop item I from list'
 
16
    X = S[I]
 
17
    del S[I]
 
18
    return X
 
19
 
 
20
#============================================================
 
21
 
 
22
class Lammps:
 
23
 
 
24
    #--------------------------------------------------------
 
25
 
 
26
    def Dump(self):
 
27
        'Write out contents of self (intended for debugging)'
 
28
        Name_list = self.__dict__.keys()
 
29
        Name_list.sort()
 
30
        for Name in Name_list:
 
31
            print Name + ':', self.__dict__[Name]
 
32
 
 
33
    #--------------------------------------------------------
 
34
 
 
35
    def Write_data(self, Basename, Item_list):
 
36
        'Write the Lammps data to file (used by Write_Lammps)'
 
37
        import os, sys
 
38
 
 
39
        Filename = 'data.' + Basename
 
40
 
 
41
        Dir_list = os.listdir('.')
 
42
        i = 1
 
43
        while Filename in Dir_list:
 
44
            Filename = 'data' + `i` + '.' + Basename
 
45
            i = i +1
 
46
        del i
 
47
 
 
48
        print 'Writing', Filename + '...',
 
49
        sys.stdout.flush()
 
50
 
 
51
        try:
 
52
            F = open(Filename, 'w')
 
53
        except IOError, Detail:
 
54
            print '(error:', Detail[1] + '!)'
 
55
            return
 
56
 
 
57
        try:
 
58
            F.writelines(Item_list)
 
59
        except IOError, Detail:
 
60
            print '(error:', Detail[1] + '!)'
 
61
            F.close()
 
62
            return
 
63
 
 
64
        F.close()
 
65
        print 'done.'
 
66
 
 
67
    #--------------------------------------------------------
 
68
 
 
69
    def Write_Lammps(self, Basename):
 
70
        'Write the Lammps data file, ignoring blank sections'
 
71
        import string
 
72
        L = []
 
73
 
 
74
        L.append('LAMMPS data file for ' + self.name + '\n\n')
 
75
 
 
76
        L.append(`self.atoms` + ' atoms\n')
 
77
        L.append(`self.bonds` + ' bonds\n')
 
78
        L.append(`self.angles` + ' angles\n')
 
79
        L.append(`self.dihedrals` + ' dihedrals\n')
 
80
        L.append(`self.impropers` + ' impropers\n\n')
 
81
 
 
82
        L.append(`self.atom_types` + ' atom types\n')
 
83
        if self.bonds > 0:
 
84
            L.append(`self.bond_types` + ' bond types\n')
 
85
        if self.angles > 0:
 
86
            L.append(`self.angle_types` + ' angle types\n')
 
87
        if self.dihedrals > 0:
 
88
            L.append(`self.dihedral_types` + ' dihedral types\n')
 
89
        L.append('\n')
 
90
        
 
91
        L.append(`self.xlo` + ' ' + `self.xhi` + ' xlo xhi\n')
 
92
        L.append(`self.ylo` + ' ' + `self.yhi` + ' ylo yhi\n')
 
93
        L.append(`self.zlo` + ' ' + `self.zhi` + ' zlo zhi\n\n')
 
94
 
 
95
        if self.atom_types != 0:
 
96
            L.append('Masses\n\n')
 
97
            for i in range(self.atom_types):
 
98
                L.append(`i+1` + ' ' + `self.Masses[i]` + '\n')
 
99
            L.append('\n')
 
100
 
 
101
            L.append('Pair Coeffs\n\n')
 
102
            for i in range(self.atom_types):
 
103
                L.append(`i+1`)
 
104
                for j in range(len(self.Nonbond_Coeffs[0])):
 
105
                    L.append(' ' + `self.Nonbond_Coeffs[i][j]`)
 
106
                L.append('\n')
 
107
            L.append('\n')
 
108
 
 
109
        if self.bonds != 0 and self.bond_types != 0:
 
110
            L.append('Bond Coeffs\n\n')
 
111
            for i in range(self.bond_types):
 
112
                L.append(`i+1`)
 
113
                for j in range(len(self.Bond_Coeffs[0])):
 
114
                    L.append(' ' + `self.Bond_Coeffs[i][j]`)
 
115
                L.append('\n')
 
116
            L.append('\n')
 
117
 
 
118
        if self.angles != 0 and self.angle_types != 0:
 
119
            L.append('Angle Coeffs\n\n')
 
120
            for i in range(self.angle_types):
 
121
                L.append(`i+1`)
 
122
                for j in range(len(self.Angle_Coeffs[0])):
 
123
                    L.append(' ' + `self.Angle_Coeffs[i][j]`)
 
124
                L.append('\n')
 
125
            L.append('\n')
 
126
 
 
127
        if self.dihedrals != 0 and self.dihedral_types != 0:
 
128
            L.append('Dihedral Coeffs\n\n')
 
129
            for i in range(self.dihedral_types):
 
130
                L.append(`i+1`)
 
131
                for j in range(len(self.Dihedral_Coeffs[0])):
 
132
                    L.append(' ' + `self.Dihedral_Coeffs[i][j]`)
 
133
                L.append('\n')
 
134
            L.append('\n')
 
135
 
 
136
        if self.atoms != 0:
 
137
            L.append('Atoms\n\n')
 
138
            for i in range(self.atoms):
 
139
                L.append(`i+1`)
 
140
                for j in range(len(self.Atoms[0])):
 
141
                    L.append(' ' + `self.Atoms[i][j]`)
 
142
                L.append('\n')
 
143
            L.append('\n')
 
144
 
 
145
        if self.bonds != 0 and self.bond_types != 0:
 
146
            L.append('Bonds\n\n')
 
147
            for i in range(self.bonds):
 
148
                L.append(`i+1`)
 
149
                for j in range(len(self.Bonds[0])):
 
150
                    L.append(' ' + `self.Bonds[i][j]`)
 
151
                L.append('\n')
 
152
            L.append('\n')
 
153
 
 
154
        if self.angles != 0 and self.angle_types != 0:
 
155
            L.append('Angles\n\n')
 
156
            for i in range(self.angles):
 
157
                L.append(`i+1`)
 
158
                for j in range(len(self.Angles[0])):
 
159
                    L.append(' ' + `self.Angles[i][j]`)
 
160
                L.append('\n')
 
161
            L.append('\n')
 
162
 
 
163
        if self.dihedrals != 0 and self.dihedral_types != 0:
 
164
            L.append('Dihedrals\n\n')
 
165
            for i in range(self.dihedrals):
 
166
                L.append(`i+1`)
 
167
                for j in range(len(self.Dihedrals[0])):
 
168
                    L.append(' ' + `self.Dihedrals[i][j]`)
 
169
                L.append('\n')
 
170
            L.append('\n')
 
171
 
 
172
        self.Write_data(Basename, L)
 
173
 
 
174
#============================================================
 
175
 
 
176
class Amber:
 
177
    def __init__(self):
 
178
        'Initialise the Amber class'
 
179
        self.CRD_is_read = 0
 
180
        self.TOP_is_read = 0
 
181
 
 
182
    #--------------------------------------------------------
 
183
 
 
184
    def Dump(self):
 
185
        'Write out contents of self (intended for debugging)'
 
186
        Name_list = self.__dict__.keys()
 
187
        Name_list.sort()
 
188
        for Name in Name_list:
 
189
            print Name + ':', self.__dict__[Name]
 
190
 
 
191
    #--------------------------------------------------------
 
192
 
 
193
    def Coerce_to_Lammps(self):
 
194
        'Return the Amber data converted to Lammps format'
 
195
 
 
196
        import math
 
197
 
 
198
        if self.CRD_is_read and self.TOP_is_read:
 
199
            l = Lammps()
 
200
            print 'Converting...',
 
201
 
 
202
            l.name = self.ITITL
 
203
            l.atoms = self.NATOM
 
204
            l.bonds = self.NBONH + self.MBONA
 
205
            l.angles = self.NTHETH + self.MTHETA
 
206
            l.dihedrals = self.NPHIH + self.MPHIA
 
207
            l.impropers = 0
 
208
            l.atom_types = self.NTYPES
 
209
            l.bond_types = self.NUMBND
 
210
            l.angle_types = self.NUMANG
 
211
            l.dihedral_types = self.NPTRA
 
212
 
 
213
            Shift = 0
 
214
            if self.__dict__.has_key('BOX'):
 
215
                l.xlo = 0.0
 
216
                l.xhi = self.BOX[0]
 
217
                l.ylo = 0.0
 
218
                l.yhi = self.BOX[1]
 
219
                l.zlo = 0.0
 
220
                l.zhi = self.BOX[2]
 
221
                if (l.xlo > min(self.X)) or (l.xhi < max(self.X)) or \
 
222
                   (l.ylo > min(self.Y)) or (l.yhi < max(self.Y)) or \
 
223
                   (l.zlo > min(self.Z)) or (l.zhi < max(self.Z)):
 
224
                    #  Vikas Modification: Disabling Shifting. This means I am intend to send exact coordinates of each atom and let LAMMPS 
 
225
                    #  take care of imaging into periodic image cells. If one wants to shift all atoms in the periodic box, 
 
226
                    #  please uncomment the below 2 lines.
 
227
                    print '(warning: Currently not shifting the atoms to the periodic box)' 
 
228
                    #Shift = 1
 
229
            else:
 
230
                print '(warning: Guessing at periodic box!)',
 
231
                l.xlo = min(self.X)
 
232
                l.xhi = max(self.X)
 
233
                l.ylo = min(self.Y)
 
234
                l.yhi = max(self.Y)
 
235
                l.zlo = min(self.Z)
 
236
                l.zhi = max(self.Z)
 
237
 
 
238
            # This doesn't check duplicate values
 
239
            l.Masses = []
 
240
            for i in range(l.atom_types):
 
241
                l.Masses.append(0)
 
242
            for i in range(self.NATOM):
 
243
                l.Masses[self.IAC[i] - 1] = self.AMASS[i]
 
244
 
 
245
            l.Nonbond_Coeffs = []
 
246
            for i in range(self.NTYPES):
 
247
                l.Nonbond_Coeffs.append([0,0])
 
248
            for i in range(self.NTYPES):
 
249
                j = self.ICO[i * (self.NTYPES + 1)] - 1
 
250
                if self.CN1[j] == 0.0:
 
251
                    l.Nonbond_Coeffs[i][0] = 0.0
 
252
                else:
 
253
                    l.Nonbond_Coeffs[i][0] = \
 
254
                                0.25 * (self.CN2[j])**2 / self.CN1[j]
 
255
                if self.CN2[j] == 0.0:
 
256
                    l.Nonbond_Coeffs[i][1] = 0.0
 
257
                else:
 
258
                    l.Nonbond_Coeffs[i][1] = \
 
259
                                (self.CN1[j] / self.CN2[j])**(1.0/6.0)
 
260
 
 
261
            l.Bond_Coeffs = []
 
262
            for i in range(self.NUMBND):
 
263
                l.Bond_Coeffs.append([0,0])
 
264
            for i in range(self.NUMBND):
 
265
                l.Bond_Coeffs[i][0] = self.RK[i]
 
266
                l.Bond_Coeffs[i][1] = self.REQ[i]
 
267
 
 
268
            l.Angle_Coeffs = []
 
269
            for i in range(self.NUMANG):
 
270
                l.Angle_Coeffs.append([0,0])
 
271
            for i in range(self.NUMANG):
 
272
                l.Angle_Coeffs[i][0] = self.TK[i]
 
273
                l.Angle_Coeffs[i][1] = (180/math.pi) * self.TEQ[i]
 
274
 
 
275
            l.Dihedral_Coeffs = []
 
276
            for i in range(self.NPTRA):
 
277
                l.Dihedral_Coeffs.append([0,0,0])
 
278
            for i in range(self.NPTRA):
 
279
                l.Dihedral_Coeffs[i][0] = self.PK[i]
 
280
                if self.PHASE[i] == 0:
 
281
                    l.Dihedral_Coeffs[i][1] = 1
 
282
                else:
 
283
                    l.Dihedral_Coeffs[i][1] = -1
 
284
                l.Dihedral_Coeffs[i][2] = int(self.PN[i])
 
285
 
 
286
            l.Atoms = []
 
287
            for i in range(self.NATOM):
 
288
                x = self.X[i]
 
289
                y = self.Y[i]
 
290
                z = self.Z[i]
 
291
                if Shift:
 
292
                    while x < l.xlo:
 
293
                        x = x + self.BOX[0]
 
294
                    while x > l.xhi:
 
295
                        x = x - self.BOX[0]
 
296
                    while y < l.ylo:
 
297
                        y = y + self.BOX[1]
 
298
                    while y > l.yhi:
 
299
                        y = y - self.BOX[1]
 
300
                    while z < l.zlo:
 
301
                        z = z + self.BOX[2]
 
302
                    while z > l.zhi:
 
303
                        z = z - self.BOX[2]
 
304
                l.Atoms.append([0, self.IAC[i], self.CHRG[i]/18.2223, \
 
305
                                x, y, z])
 
306
 
 
307
            l.Bonds = []
 
308
            for i in range(l.bonds):
 
309
                l.Bonds.append([0,0,0])
 
310
            for i in range(self.NBONH):
 
311
                l.Bonds[i][0] = self.ICBH[i]
 
312
                l.Bonds[i][1] = abs(self.IBH[i])/3 + 1
 
313
                l.Bonds[i][2] = abs(self.JBH[i])/3 + 1
 
314
            for i in range(self.NBONA):
 
315
                l.Bonds[self.NBONH + i][0] = self.ICB[i]
 
316
                l.Bonds[self.NBONH + i][1] = abs(self.IB[i])/3 + 1
 
317
                l.Bonds[self.NBONH + i][2] = abs(self.JB[i])/3 + 1
 
318
 
 
319
            l.Angles = []
 
320
            for i in range(l.angles):
 
321
                l.Angles.append([0,0,0,0])
 
322
            for i in range(self.NTHETH):
 
323
                l.Angles[i][0] = self.ICTH[i]
 
324
                l.Angles[i][1] = abs(self.ITH[i])/3 + 1
 
325
                l.Angles[i][2] = abs(self.JTH[i])/3 + 1
 
326
                l.Angles[i][3] = abs(self.KTH[i])/3 + 1
 
327
            for i in range(self.NTHETA):
 
328
                l.Angles[self.NTHETH + i][0] = self.ICT[i]
 
329
                l.Angles[self.NTHETH + i][1] = abs(self.IT[i])/3 + 1
 
330
                l.Angles[self.NTHETH + i][2] = abs(self.JT[i])/3 + 1
 
331
                l.Angles[self.NTHETH + i][3] = abs(self.KT[i])/3 + 1
 
332
 
 
333
            l.Dihedrals = []
 
334
            for i in range(l.dihedrals):
 
335
                l.Dihedrals.append([0,0,0,0,0])
 
336
            for i in range(self.NPHIH):
 
337
                l.Dihedrals[i][0] = self.ICPH[i]
 
338
                l.Dihedrals[i][1] = abs(self.IPH[i])/3 + 1
 
339
                l.Dihedrals[i][2] = abs(self.JPH[i])/3 + 1
 
340
                l.Dihedrals[i][3] = abs(self.KPH[i])/3 + 1
 
341
                l.Dihedrals[i][4] = abs(self.LPH[i])/3 + 1
 
342
            for i in range(self.NPHIA):
 
343
                l.Dihedrals[self.NPHIH + i][0] = self.ICP[i]
 
344
                l.Dihedrals[self.NPHIH + i][1] = abs(self.IP[i])/3 + 1
 
345
                l.Dihedrals[self.NPHIH + i][2] = abs(self.JP[i])/3 + 1
 
346
                l.Dihedrals[self.NPHIH + i][3] = abs(self.KP[i])/3 + 1
 
347
                l.Dihedrals[self.NPHIH + i][4] = abs(self.LP[i])/3 + 1
 
348
 
 
349
            print 'done.'
 
350
            return l
 
351
        else:
 
352
            print '(Error: Not all the Amber data has been read!)'
 
353
    
 
354
    #--------------------------------------------------------
 
355
 
 
356
    def Read_data(self, Filename):
 
357
        'Read the filename, returning a list of strings'
 
358
 
 
359
        import string, sys
 
360
 
 
361
        print 'Reading', Filename + '...',
 
362
        sys.stdout.flush()
 
363
 
 
364
        try:
 
365
            F = open(Filename)
 
366
        except IOError, Detail:
 
367
            print '(error:', Detail[1] + '!)'
 
368
            return
 
369
 
 
370
        try:
 
371
            Lines = F.readlines()
 
372
        except IOError, Detail:
 
373
            print '(error:', Detail[1] + '!)'
 
374
            F.close()
 
375
            return
 
376
 
 
377
        F.close()
 
378
 
 
379
        # If the first line is empty, use the Basename
 
380
        if Filename[-4:] == '.crd':
 
381
            if string.split(Lines[0]) == []: # This line corresponds to TITLE name in CRD file
 
382
                Basename = Filename[:string.find(Filename, '.')]
 
383
                Item_list = [Basename]
 
384
                print 'Warning: Title not present... Assigning Basename as Title'
 
385
            else:
 
386
                Item_list = []
 
387
        else:
 
388
            if string.split(Lines[3]) == []: # This line corresponds to TITLE name in TOPOLOGY file
 
389
                Basename = Filename[:string.find(Filename, '.')]
 
390
                Item_list = [Basename]
 
391
                print 'Warning: Title not present... Assigning Basename as Title'
 
392
            else:
 
393
                Item_list = []
 
394
 
 
395
        for Line in Lines:
 
396
            if Line[0]!='%': #Vikas' Modification: This condition ignores all the lines starting with % in the topology file.
 
397
                Item_list.extend(string.split(Line))
 
398
 
 
399
        return Item_list
 
400
 
 
401
    #--------------------------------------------------------
 
402
 
 
403
    def Read_CRD(self, Basename):
 
404
        'Read the Amber coordinate/restart (.crd) file'
 
405
 
 
406
        # The optional velocities and periodic box size are not yet parsed.
 
407
 
 
408
        Item_list = self.Read_data(Basename + '.crd')
 
409
 
 
410
        if Item_list == None:
 
411
            return
 
412
        elif len(Item_list) < 2:
 
413
            print '(error: File too short!)'
 
414
            return
 
415
 
 
416
        # Parse the data
 
417
        if self.__dict__.has_key('ITITL'):
 
418
            if Pop(Item_list,0) != self.ITITL:
 
419
                print '(warning: ITITL differs!)',
 
420
        else:
 
421
            self.ITITL = Pop(Item_list,0)
 
422
        print self.ITITL #Vikas Modification : Priting the Title
 
423
 
 
424
        if self.__dict__.has_key('NATOM'):
 
425
            if eval(Pop(Item_list,0)) != self.NATOM:
 
426
                print '(error: NATOM differs!)'
 
427
                return
 
428
        else:
 
429
            self.NATOM = eval(Pop(Item_list,0))
 
430
        print self.NATOM # Vikas' Modification: Printing number of atoms just to make sure that the program is reading the correct value. 
 
431
 
 
432
        #if len(Item_list) == 1 + 3 * self.NATOM:
 
433
        # Vikas' Modification: I changed the condition.
 
434
        if (len(Item_list)%3) != 0:
 
435
            self.TIME = eval(Pop(Item_list,0))
 
436
        else:
 
437
            self.TIME = 0
 
438
        print self.TIME # Vikas' Modification : Printing simulation time, just to make sure that the program is readint the correct value.
 
439
        if len(Item_list) < 3 * self.NATOM:
 
440
            print '(error: File too short!)'
 
441
            return
 
442
 
 
443
        self.X = []
 
444
        self.Y = []
 
445
        self.Z = []
 
446
        for i in range(self.NATOM):
 
447
            self.X.append(eval(Pop(Item_list,0)))
 
448
            self.Y.append(eval(Pop(Item_list,0)))
 
449
            self.Z.append(eval(Pop(Item_list,0)))
 
450
 
 
451
        if (self.NATOM == 1) and len(Item_list):
 
452
            print '(warning: Ambiguity!)',
 
453
 
 
454
        if len(Item_list) >= 3 * self.NATOM:
 
455
                self.VX = []
 
456
                self.VY = []
 
457
                self.VZ = []
 
458
                for i in range(self.NATOM):
 
459
                    self.VX.append(eval(Pop(Item_list,0)))
 
460
                    self.VY.append(eval(Pop(Item_list,0)))
 
461
                    self.VZ.append(eval(Pop(Item_list,0)))
 
462
 
 
463
        if len(Item_list) >= 3:
 
464
            self.BOX = []
 
465
            for i in range(3):
 
466
                self.BOX.append(eval(Pop(Item_list,0)))
 
467
 
 
468
        if len(Item_list):
 
469
            print '(warning: File too large!)',
 
470
        
 
471
        print 'done.'
 
472
        self.CRD_is_read = 1
 
473
 
 
474
    #--------------------------------------------------------
 
475
 
 
476
    def Read_TOP(self, Basename):
 
477
        'Read the Amber parameter/topology (.top) file'
 
478
        Item_list = self.Read_data(Basename + '.top')
 
479
 
 
480
        if Item_list == None:
 
481
            return
 
482
        elif len(Item_list) < 31:
 
483
            print '(error: File too short!)'
 
484
            return
 
485
 
 
486
       # Parse the data
 
487
        if self.__dict__.has_key('ITITL'):
 
488
            if Pop(Item_list,0) != self.ITITL:
 
489
                print '(warning: ITITL differs!)'
 
490
        else:
 
491
            self.ITITL = Pop(Item_list,0)
 
492
        print self.ITITL # Printing Self Title
 
493
 
 
494
        if self.__dict__.has_key('NATOM'):
 
495
            if eval(Pop(Item_list,0)) != self.NATOM:
 
496
                print '(error: NATOM differs!)'
 
497
                return
 
498
        else:
 
499
            self.NATOM = eval(Pop(Item_list,0))
 
500
        print self.NATOM # Printing total number of atoms just to make sure that thing are going right
 
501
        self.NTYPES = eval(Pop(Item_list,0))
 
502
        self.NBONH  = eval(Pop(Item_list,0))
 
503
        self.MBONA  = eval(Pop(Item_list,0))
 
504
        self.NTHETH = eval(Pop(Item_list,0))
 
505
        self.MTHETA = eval(Pop(Item_list,0))
 
506
        self.NPHIH  = eval(Pop(Item_list,0))
 
507
        self.MPHIA  = eval(Pop(Item_list,0))
 
508
        self.NHPARM = eval(Pop(Item_list,0))
 
509
        self.NPARM  = eval(Pop(Item_list,0))
 
510
        self.NEXT   = eval(Pop(Item_list,0))
 
511
        self.NRES   = eval(Pop(Item_list,0))
 
512
        self.NBONA  = eval(Pop(Item_list,0))
 
513
        self.NTHETA = eval(Pop(Item_list,0))
 
514
        self.NPHIA  = eval(Pop(Item_list,0))
 
515
        self.NUMBND = eval(Pop(Item_list,0))
 
516
        self.NUMANG = eval(Pop(Item_list,0))
 
517
        self.NPTRA  = eval(Pop(Item_list,0))
 
518
        self.NATYP  = eval(Pop(Item_list,0))
 
519
        self.NPHB   = eval(Pop(Item_list,0))
 
520
        self.IFPERT = eval(Pop(Item_list,0))
 
521
        self.NBPER  = eval(Pop(Item_list,0))
 
522
        self.NGPER  = eval(Pop(Item_list,0))
 
523
        self.NDPER  = eval(Pop(Item_list,0))
 
524
        self.MBPER  = eval(Pop(Item_list,0))
 
525
        self.MGPER  = eval(Pop(Item_list,0))
 
526
        self.MDPER  = eval(Pop(Item_list,0))
 
527
        self.IFBOX  = eval(Pop(Item_list,0))
 
528
        self.NMXRS  = eval(Pop(Item_list,0))
 
529
        self.IFCAP  = eval(Pop(Item_list,0))
 
530
 
 
531
        #....................................................
 
532
 
 
533
        if len(Item_list) < 5 * self.NATOM + self.NTYPES**2 + \
 
534
           2*(self.NRES + self.NUMBND + self.NUMANG) + \
 
535
           3*self.NPTRA + self.NATYP:
 
536
            print '(error: File too short!)'
 
537
            return -1
 
538
 
 
539
        self.IGRAPH = []
 
540
        Pop(Item_list,0)
 
541
 
 
542
          # A little kludge is needed here, since the IGRAPH strings are
 
543
        # not separated by spaces if 4 characters in length.
 
544
        for i in range(self.NATOM):
 
545
            if len(Item_list[0]) > 4:
 
546
                Item_list.insert(1, Item_list[0][4:])
 
547
                Item_list.insert(1, Item_list[0][0:4])
 
548
                del Item_list[0]
 
549
            self.IGRAPH.append(Pop(Item_list,0))
 
550
        
 
551
        # Vikas' Modification : In the following section, I am printing out each quantity which is currently being read from the topology file.
 
552
        print 'Reading Charges...'
 
553
        self.CHRG = []
 
554
        for i in range(self.NATOM):
 
555
            self.CHRG.append(eval(Pop(Item_list,0)))
 
556
 
 
557
        print 'Reading Atomic Number...'
 
558
        self.ANUMBER = []
 
559
        for i in range(self.NATOM):
 
560
            self.ANUMBER.append(eval(Pop(Item_list,0)))
 
561
 
 
562
        print 'Reading Atomic Masses...'
 
563
        self.AMASS = []
 
564
        for i in range(self.NATOM):
 
565
            self.AMASS.append(eval(Pop(Item_list,0)))
 
566
 
 
567
        print 'Reading Atom Types...'
 
568
        self.IAC = []
 
569
        for i in range(self.NATOM):
 
570
            self.IAC.append(eval(Pop(Item_list,0)))
 
571
 
 
572
        print 'Reading Excluded Atoms...'
 
573
        self.NUMEX = []
 
574
        for i in range(self.NATOM):
 
575
            self.NUMEX.append(eval(Pop(Item_list,0)))
 
576
 
 
577
        print 'Reading Non-bonded Parameter Index...'
 
578
        self.ICO = []
 
579
        for i in range(self.NTYPES**2):
 
580
            self.ICO.append(eval(Pop(Item_list,0)))
 
581
 
 
582
        print 'Reading Residue Labels...'
 
583
        self.LABRES = []
 
584
        for i in range(self.NRES):
 
585
            self.LABRES.append(Pop(Item_list,0))
 
586
 
 
587
        print 'Reading Residues Starting Pointers...'
 
588
        self.IPRES = []
 
589
        for i in range(self.NRES):
 
590
            self.IPRES.append(eval(Pop(Item_list,0)))
 
591
 
 
592
        print 'Reading Bond Force Constants...'
 
593
        self.RK = []
 
594
        for i in range(self.NUMBND):
 
595
            self.RK.append(eval(Pop(Item_list,0)))
 
596
 
 
597
        print 'Reading Equilibrium Bond Values...'
 
598
        self.REQ = []
 
599
        for i in range(self.NUMBND):
 
600
            self.REQ.append(eval(Pop(Item_list,0)))
 
601
 
 
602
        print 'Reading Angle Force Constants...' 
 
603
        self.TK = []
 
604
        for i in range(self.NUMANG):
 
605
            self.TK.append(eval(Pop(Item_list,0)))
 
606
 
 
607
        print 'Reading Equilibrium Angle Values...'
 
608
        self.TEQ = []
 
609
        for i in range(self.NUMANG):
 
610
            self.TEQ.append(eval(Pop(Item_list,0)))
 
611
 
 
612
        print 'Reading Dihedral Force Constants...'
 
613
        self.PK = []
 
614
        for i in range(self.NPTRA):
 
615
            self.PK.append(eval(Pop(Item_list,0)))
 
616
        
 
617
        print 'Reading Dihedral Periodicity...' 
 
618
        self.PN = []
 
619
        for i in range(self.NPTRA):
 
620
            self.PN.append(eval(Pop(Item_list,0)))
 
621
        
 
622
        print 'Reading Dihedral Phase...'
 
623
        self.PHASE = []
 
624
        for i in range(self.NPTRA):
 
625
            self.PHASE.append(eval(Pop(Item_list,0)))
 
626
        
 
627
        print 'Reading 1-4 Electrostatic Scaling Factor...'
 
628
        self.SCEEFAC = []
 
629
        for i in range(self.NPTRA):
 
630
            self.SCEEFAC.append(eval(Pop(Item_list,0)))
 
631
        
 
632
        print 'Reading 1-4 Van der Waals Scaling Factor...'
 
633
        self.SCNBFAC = []
 
634
        for i in range(self.NPTRA):
 
635
            self.SCNBFAC.append(eval(Pop(Item_list,0)))
 
636
        
 
637
        print 'Reading Solty...' #I think this is currently not used in AMBER. Check it out, though
 
638
        self.SOLTY = []
 
639
        for i in range(self.NATYP):
 
640
            self.SOLTY.append(eval(Pop(Item_list,0)))
 
641
 
 
642
        #....................................................
 
643
 
 
644
        if len(Item_list) < 2 * self.NTYPES * (self.NTYPES + 1) / 2:
 
645
            print '(error: File too short!)'
 
646
            return -1
 
647
 
 
648
        print 'Reading LJ A Coefficient...'
 
649
        self.CN1 = []
 
650
        for i in range(self.NTYPES * (self.NTYPES + 1) / 2):
 
651
            self.CN1.append(eval(Pop(Item_list,0)))
 
652
 
 
653
        print 'Reading LJ B Coefficient...'
 
654
        self.CN2 = []
 
655
        for i in range(self.NTYPES * (self.NTYPES + 1) / 2):
 
656
            self.CN2.append(eval(Pop(Item_list,0)))
 
657
 
 
658
        #....................................................
 
659
 
 
660
        if len(Item_list) < 3 * (self.NBONH + self.NBONA) + \
 
661
           4 * (self.NTHETH + self.NTHETA) + 5 * (self.NPHIH + self.NPHIA):
 
662
            print '(error: File too short!)'
 
663
            return -1
 
664
 
 
665
        print 'Reading Bonds which include hydrogen...'
 
666
        self.IBH = []
 
667
        self.JBH = []
 
668
        self.ICBH = []
 
669
        for i in range(self.NBONH):
 
670
            self.IBH.append(eval(Pop(Item_list,0)))
 
671
            self.JBH.append(eval(Pop(Item_list,0)))
 
672
            self.ICBH.append(eval(Pop(Item_list,0)))
 
673
 
 
674
        print 'Reading Bonds which dont include hydrogen...'
 
675
        self.IB = []
 
676
        self.JB = []
 
677
        self.ICB = []
 
678
        for i in range(self.NBONA):
 
679
            self.IB.append(eval(Pop(Item_list,0)))
 
680
            self.JB.append(eval(Pop(Item_list,0)))
 
681
            self.ICB.append(eval(Pop(Item_list,0)))
 
682
 
 
683
        print 'Reading Angles which include hydrogen...' 
 
684
        self.ITH = []
 
685
        self.JTH = []
 
686
        self.KTH = []
 
687
        self.ICTH = []
 
688
        for i in range(self.NTHETH):
 
689
            self.ITH.append(eval(Pop(Item_list,0)))
 
690
            self.JTH.append(eval(Pop(Item_list,0)))
 
691
            self.KTH.append(eval(Pop(Item_list,0)))
 
692
            self.ICTH.append(eval(Pop(Item_list,0)))
 
693
 
 
694
        print 'Reading Angles which dont include hydrogen...'
 
695
        self.IT = []
 
696
        self.JT = []
 
697
        self.KT = []
 
698
        self.ICT = []
 
699
        for i in range(self.NTHETA):
 
700
            self.IT.append(eval(Pop(Item_list,0)))
 
701
            self.JT.append(eval(Pop(Item_list,0)))
 
702
            self.KT.append(eval(Pop(Item_list,0)))
 
703
            self.ICT.append(eval(Pop(Item_list,0)))
 
704
 
 
705
        print 'Reading Dihedrals which include hydrogen...'
 
706
        self.IPH = []
 
707
        self.JPH = []
 
708
        self.KPH = []
 
709
        self.LPH = []
 
710
        self.ICPH = []
 
711
        for i in range(self.NPHIH):
 
712
            self.IPH.append(eval(Pop(Item_list,0)))
 
713
            self.JPH.append(eval(Pop(Item_list,0)))
 
714
            self.KPH.append(eval(Pop(Item_list,0)))
 
715
            self.LPH.append(eval(Pop(Item_list,0)))
 
716
            self.ICPH.append(eval(Pop(Item_list,0)))
 
717
 
 
718
        print 'Reading Dihedrals which dont include hydrogen...'
 
719
        self.IP = []
 
720
        self.JP = []
 
721
        self.KP = []
 
722
        self.LP = []
 
723
        self.ICP = []
 
724
        for i in range(self.NPHIA):
 
725
            self.IP.append(eval(Pop(Item_list,0)))
 
726
            self.JP.append(eval(Pop(Item_list,0)))
 
727
            self.KP.append(eval(Pop(Item_list,0)))
 
728
            self.LP.append(eval(Pop(Item_list,0)))
 
729
            self.ICP.append(eval(Pop(Item_list,0)))
 
730
 
 
731
        #....................................................
 
732
 
 
733
        if len(Item_list) < self.NEXT + 3 * self.NPHB + 4 * self.NATOM:
 
734
            print '(error: File too short!)'
 
735
            return -1
 
736
 
 
737
        print 'Reading Excluded Atom List...'
 
738
        self.NATEX = []
 
739
        for i in range(self.NEXT):
 
740
            self.NATEX.append(eval(Pop(Item_list,0)))
 
741
 
 
742
        print 'Reading H-Bond A Coefficient, corresponding to r**12 term for all possible types...'
 
743
        self.ASOL = []
 
744
        for i in range(self.NPHB):
 
745
            self.ASOL.append(eval(Pop(Item_list,0)))
 
746
 
 
747
        print 'Reading H-Bond B Coefficient, corresponding to r**10 term for all possible types...'
 
748
        self.BSOL = []
 
749
        for i in range(self.NPHB):
 
750
            self.BSOL.append(eval(Pop(Item_list,0)))
 
751
 
 
752
        print 'Reading H-Bond Cut...' # I think it is not being used nowadays
 
753
        self.HBCUT = []
 
754
        for i in range(self.NPHB):
 
755
            self.HBCUT.append(eval(Pop(Item_list,0)))
 
756
 
 
757
        print 'Reading Amber Atom Types for each atom...'
 
758
        self.ISYMBL = []
 
759
        for i in range(self.NATOM):
 
760
            self.ISYMBL.append(Pop(Item_list,0))
 
761
 
 
762
        print 'Reading Tree Chain Classification...'
 
763
        self.ITREE = []
 
764
        for i in range(self.NATOM):
 
765
            self.ITREE.append(Pop(Item_list,0))
 
766
 
 
767
        print 'Reading Join Array: Tree joining information' # Currently unused in Sander, an AMBER module
 
768
        self.JOIN = []
 
769
        for i in range(self.NATOM):
 
770
            self.JOIN.append(eval(Pop(Item_list,0)))
 
771
 
 
772
        print 'Reading IRotate...' # Currently unused in Sander and Gibbs
 
773
        self.IROTAT = []
 
774
        for i in range(self.NATOM):
 
775
            self.IROTAT.append(eval(Pop(Item_list,0)))
 
776
 
 
777
        #....................................................
 
778
 
 
779
        if self.IFBOX > 0:
 
780
            if len(Item_list) < 3:
 
781
                print '(error: File too short!)'
 
782
                return -1
 
783
 
 
784
            print 'Reading final residue which is part of solute...'
 
785
            self.IPTRES = eval(Pop(Item_list,0))
 
786
            print 'Reading total number of molecules...'
 
787
            self.NSPM = eval(Pop(Item_list,0))
 
788
            print 'Reading first solvent moleule index...'
 
789
            self.NSPSOL = eval(Pop(Item_list,0))
 
790
 
 
791
            if len(Item_list) < self.NSPM + 4:
 
792
                print '(error: File too short!)'
 
793
                return -1
 
794
 
 
795
            print 'Reading atom per molecule...' 
 
796
            self.NSP = []
 
797
            for i in range(self.NSPM):
 
798
                self.NSP.append(eval(Pop(Item_list,0)))
 
799
 
 
800
            self.BETA = eval(Pop(Item_list,0))
 
801
 
 
802
            print 'Reading Box Dimensions...' 
 
803
            if self.__dict__.has_key('BOX'):
 
804
                BOX = []
 
805
                for i in range(3):
 
806
                    BOX.append(eval(Pop(Item_list,0)))
 
807
                for i in range(3):
 
808
                    if BOX[i] != self.BOX[i]:
 
809
                        print '(warning: BOX differs!)',
 
810
                        break
 
811
                del BOX
 
812
            else:
 
813
                self.BOX = []
 
814
                for i in range(3):
 
815
                    self.BOX.append(eval(Pop(Item_list,0)))
 
816
 
 
817
        #....................................................
 
818
 
 
819
        if self.IFCAP > 0:
 
820
            if len(Item_list) < 5:
 
821
                print '(error: File too short!)'
 
822
                return -1
 
823
            print 'Reading ICAP variables::: For details, refer to online AMBER format manual'
 
824
            self.NATCAP = eval(Pop(Item_list,0))
 
825
            self.CUTCAP = eval(Pop(Item_list,0))
 
826
            self.XCAP = eval(Pop(Item_list,0))
 
827
            self.YCAP = eval(Pop(Item_list,0))
 
828
            self.ZCAP = eval(Pop(Item_list,0))
 
829
 
 
830
        #....................................................
 
831
 
 
832
        if self.IFPERT > 0:
 
833
            if len(Item_list) < 4 * self.NBPER + 5 * self.NGPER + \
 
834
               6 * self.NDPER + self.NRES + 6 * self.NATOM:
 
835
                print '(error: File too short!)'
 
836
                return -1
 
837
 
 
838
            print 'Reading perturb variables, 1. Bond, 2. Angles, 3. Dihedrals, etc etc.::: For details, refer to online AMBER format manual' 
 
839
            self.IBPER = []
 
840
            self.JBPER = []
 
841
            for i in range(self.NBPER):
 
842
                self.IBPER.append(eval(Pop(Item_list,0)))
 
843
                self.JBPER.append(eval(Pop(Item_list,0)))
 
844
 
 
845
            self.ICBPER = []
 
846
            for i in range(2 * self.NBPER):
 
847
                self.ICBPER.append(eval(Pop(Item_list,0)))
 
848
 
 
849
            self.ITPER = []
 
850
            self.JTPER = []
 
851
            self.KTPER = []
 
852
            for i in range(self.NGPER):
 
853
                self.ITPER.append(eval(Pop(Item_list,0)))
 
854
                self.JTPER.append(eval(Pop(Item_list,0)))
 
855
                self.KTPER.append(eval(Pop(Item_list,0)))
 
856
 
 
857
            self.ICTPER = []
 
858
            for i in range(2 * self.NGPER):
 
859
                self.ICTPER.append(eval(Pop(Item_list,0)))
 
860
 
 
861
            self.IPPER = []
 
862
            self.JPPER = []
 
863
            self.KPPER = []
 
864
            self.LPPER = []
 
865
            for i in range(self.NDPER):
 
866
                self.IPPER.append(eval(Pop(Item_list,0)))
 
867
                self.JPPER.append(eval(Pop(Item_list,0)))
 
868
                self.KPPER.append(eval(Pop(Item_list,0)))
 
869
                self.LPPER.append(eval(Pop(Item_list,0)))
 
870
 
 
871
            self.ICPPER = []
 
872
            for i in range(2 * self.NDPER):
 
873
                self.ICPPER.append(eval(Pop(Item_list,0)))
 
874
 
 
875
            LABRES = []
 
876
            for i in range(self.NRES):
 
877
                LABRES.append(Pop(Item_list,0))
 
878
            for i in range(self.NRES):
 
879
                if LABRES[i] != self.LABRES[i]:
 
880
                    print '(warning: BOX differs!)',
 
881
                    break
 
882
 
 
883
            self.IGRPER = []
 
884
            for i in range(self.NATOM):
 
885
                self.IGRPER.append(eval(Pop(Item_list,0)))
 
886
 
 
887
            self.ISMPER = []
 
888
            for i in range(self.NATOM):
 
889
                self.ISMPER.append(eval(Pop(Item_list,0)))
 
890
 
 
891
            self.ALMPER = []
 
892
            for i in range(self.NATOM):
 
893
                self.ALMPER.append(eval(Pop(Item_list,0)))
 
894
 
 
895
            self.IAPER = []
 
896
            for i in range(self.NATOM):
 
897
                self.IAPER.append(eval(Pop(Item_list,0)))
 
898
 
 
899
            self.IACPER = []
 
900
            for i in range(self.NATOM):
 
901
                self.IACPER.append(eval(Pop(Item_list,0)))
 
902
 
 
903
            self.CGPER = []
 
904
            for i in range(self.NATOM):
 
905
                self.CGPER.append(eval(Pop(Item_list,0)))
 
906
 
 
907
        #....................................................
 
908
 
 
909
        self.IPOL = 0
 
910
        if self.IPOL == 1:
 
911
            if len(Item_list) < self.NATOM:
 
912
                print '(error: File too short!)'
 
913
                return -1
 
914
            print 'Reading Polarizability Data. For details, refer to online AMBER format manual'
 
915
            self.ATPOL = []
 
916
            for i in range(self.NATOM):
 
917
                self.ATPOL.append(eval(Pop(Item_list,0)))
 
918
 
 
919
            if self.IFPERT == 1:
 
920
                if len(Item_list) < self.NATOM:
 
921
                    print '(error: File too short!)'
 
922
                    return -1
 
923
                self.ATPOL1 = []
 
924
                for i in range(self.NATOM):
 
925
                    self.ATPOL1.append(eval(Pop(Item_list,0)))
 
926
 
 
927
        #....................................................
 
928
 
 
929
        if len(Item_list):
 
930
            print '(warning: File too large!)',
 
931
 
 
932
        print 'done.'
 
933
        self.TOP_is_read = 1
 
934
 
 
935
#============================================================
 
936
 
 
937
def Find_Amber_files():
 
938
    'Look for sets of Amber files to process'
 
939
    '''If not passed anything on the command line, look for pairs of
 
940
    Amber files (.crd and .top) in the current directory.  For
 
941
    each set if there is no corresponding Lammps file (data.), or it is
 
942
    older than any of the Amber files, add its basename to a list of
 
943
    strings.  This list is returned by the function'''
 
944
 
 
945
    # Date and existence checks not yet implemented
 
946
 
 
947
    import os, sys
 
948
 
 
949
    Basename_list = []
 
950
 
 
951
    # Extract basenames from command line
 
952
    for Name in sys.argv[1:]:
 
953
        if Name[-4:] == '.crd':
 
954
            Basename_list.append(Name[:-4])
 
955
        else:
 
956
            if Name[-4:] == '.top':
 
957
                Basename_list.append(Name[:-4])
 
958
            else:
 
959
                Basename_list.append(Name)
 
960
 
 
961
    # Remove duplicate basenames
 
962
    for Basename in Basename_list[:]:
 
963
        while Basename_list.count(Basename) > 1:
 
964
            Basename_list.remove(Basename)
 
965
 
 
966
    if Basename_list == []:
 
967
        print 'Looking for Amber files...',
 
968
        Dir_list = os.listdir('.')
 
969
        Dir_list.sort()
 
970
        for File in Dir_list:
 
971
            if File[-4:] == '.top':
 
972
                Basename = File[:-4]
 
973
                if (Basename + '.crd') in Dir_list:
 
974
                    Basename_list.append(Basename)
 
975
        if Basename_list != []:
 
976
            print 'found',
 
977
            for i in range(len(Basename_list)-1):
 
978
                print Basename_list[i] + ',',
 
979
            print Basename_list[-1] + '\n'
 
980
    
 
981
    if Basename_list == []:
 
982
        print 'none.\n'
 
983
 
 
984
    return Basename_list
 
985
    
 
986
#============================================================
 
987
 
 
988
def Convert_Amber_files():
 
989
    'Handle the whole conversion process'
 
990
    print
 
991
    print 'Welcome to amber2lammps, a program to convert Amber files to Lammps format!'
 
992
    print
 
993
    Basename_list = Find_Amber_files()
 
994
    for Basename in Basename_list:
 
995
        a = Amber()
 
996
        a.Read_CRD(Basename)
 
997
        if a.CRD_is_read:
 
998
            a.Read_TOP(Basename)
 
999
            if a.TOP_is_read:
 
1000
                l = a.Coerce_to_Lammps()
 
1001
                l.Write_Lammps(Basename)
 
1002
                del l
 
1003
        del a
 
1004
        print
 
1005
 
 
1006
#============================================================
 
1007
 
 
1008
Convert_Amber_files()