1#!/usr/local/bin/perl
2
3# usage: make_gromos_rtp.pl mtb43a1.dat > gromos.rtp
4# this script tries to make a residue topology database for GROMACS from
5# the GROMOS version of this file. It needs ffG96A.atp to fill in the
6# atom types for the atom integer codes. It converts until it finds the string
7# "#RNMES", which indicates the start of the solvent part
8# of mtb43(a,b)1.dat. Solvents are treated differently in GROMACS
9
10# author: Peter Tieleman, 12 March 1999
11
12# read in the atomtypes
13open (TYPES, "gromos.atp") || die "Can't open gromos.atp: $!";
14while (<TYPES>) {
15    $type++;
16    ($gromostype[$type],$rest) = split(' ',$_);
17}
18
19$/ = "# RNME";
20# split input on residue name. The last line of each block of residue data
21# has RNME in it. The first line has the name of the residue.
22
23# write header
24print "[ bondedtypes ]\n";
25print "; bonds  angles  dihedrals  impropers\n";
26print "    2       2          1          2\n\n";
27
28$first = 1;
29
30# loop over the actual residues
31while (<>) {
32    if ($first) {$first = 0; next;} #skip over the first nonsense residue
33    # initialise for a new residue:
34    $read_atoms = 0; $j = 0;
35    $read_bonds = 0; $k = 0;
36    $read_angles = 0; $l = 0;
37    $read_dihedrals = 0; $m = 0;
38    $read_impropers = 0; $n = 0;
39
40    # now $_ has all data for a single residue in it.
41    # split it up in lines
42    @residue_data = split('\n',$_);
43    $residue_name = @residue_data[1];
44
45    # loop over the lines of a residue
46    for ($i=0;$i<@residue_data;$i++) {
47	# do we have to skip a line?
48	if ($skip == 1) {$skip = 0; next;}
49
50	# look for the headers of each section, and move to the start
51	# of the relevant data.
52
53	# if you find ATOM ANM, start reading atoms
54	if ($residue_data[$i] =~ /\#ATOM ANM/) {
55	    $read_atoms = 1;
56	}
57	# if you find NB, read bonds
58	if ($residue_data[$i] =~ /\#\s+NB/) {
59	    $read_atoms = 0; $read_bonds = 1; $skip = 1; #skip next line
60	}
61	# if you find NBA, read angles
62	if ($residue_data[$i] =~ /NBA/) {
63	    $read_bonds = 0; $read_angles = 1; $skip = 1; #skip next line
64	}
65	# if you find NIDA, read impropers
66	if ($residue_data[$i] =~ /NIDA/) {
67	    $read_angles =0; $read_impropers = 1; $skip = 1; #skip next line
68	}
69	# if you find NDA, read dihedrals
70	if ($residue_data[$i] =~ /NDA/) {
71	    $read_impropers= 0; $read_dihedrals=1; $skip = 1; #skip next line
72	}
73
74	# stop parsing alltogether if we find #RNMES
75	if ($residue_data[$i] =~ /\#RNMES/) { last;}
76        # skip lines that start with #
77	if ( $residue_data[$i] =~ /^\#/) { next; }
78        # also skip lines with END on it, and with MTBUILDBLSOLUTE
79	if ($residue_data[$i] =~ /END/) { next; }
80        if ($residue_data[$i] =~ /BUILD/) {next;}
81
82        if (!($read_atoms || $read_bonds || $read_angles || $read_dihedrals
83	      || $read_impropers)) { next; } # we're not reading anything yet
84
85	if ($read_atoms) {
86	    ($atomnr[$j],$atomname[$j],$atomtype[$j],$mass,$charge[$j],
87	     $chargegroup[$j],$dummy) = split(' ',$residue_data[$i]);
88	    # some lines only have exclusions, check if there is a name
89	    if ($atomname[$j] !~ /[A-Z]+/) {$j--;}
90	    $j++;
91	}
92	if ($read_bonds) {
93	    ($bondi[$k],$bondj[$k],$bondtype[$k]) =
94		split(' ',$residue_data[$i]);
95	    $k++;
96	}
97	if ($read_angles) {
98	    ($anglei[$l],$anglej[$l],$anglek[$l],$angletype[$l]) =
99		split(' ',$residue_data[$i]);
100	    $l++;
101	}
102
103	if ($read_impropers) {
104	    ($imp_i[$m],$imp_j[$m],$imp_k[$m],$imp_l[$m],$imp_type[$m]) =
105		split(' ',$residue_data[$i]);
106	    $m++;
107	}
108
109        if ($read_dihedrals) {
110	    ($dih_i[$n],$dih_j[$n],$dih_k[$n],$dih_l[$n],$dih_type[$n]) =
111		split(' ',$residue_data[$i]);
112	    $n++;
113	}
114    }
115    $natoms = $j;
116
117    #  print out this residue to the GROMACS file and go to the next residue
118    print "[ $residue_name ]\n";
119    print " [ atoms ]\n";
120    $chargegroup = 0;
121    for ($t=0;$t<$j;$t++) {
122	print sprintf("%5s %5s %8.3f %5s\n", $atomname[$t],
123		      $gromostype[$atomtype[$t]], $charge[$t], $chargegroup);
124	$chargegroup += $chargegroup[$t];
125    }
126
127    print " [ bonds ]\n";
128    for ($t=0;$t<$k;$t++) {
129	$bond ="gb_$bondtype[$t]";
130	# convert the nrs from GROMOS to atomnames
131	if ($bondi[$t] < 0 || $bondj[$t] < 0) {
132	    print STDERR "Skipping specbond $bondi[$t] $bondj[$t]\n";
133	    next;
134	}
135	if ($bondi[$t] == $natoms+1) {$ati = "+N";}
136	else {$ati = $atomname[$bondi[$t]-1];}
137	if ($bondj[$t] == $natoms+1) {$atj = "+N";}
138	else {$atj = $atomname[$bondj[$t]-1];}
139	# write them out.
140	print sprintf("%5s %5s    %-5s\n", $ati, $atj, $bond);
141    }
142
143    print " [ angles ]\n";
144    print ";   ai    aj    ak  gromos type\n";
145    for ($t=0;$t<$l;$t++) {
146
147	# convert numbers to names.
148	# i may be in the previous residue
149	if ($anglei[$t] == -1) {$ati = "-C";}
150	else {$ati = $atomname[$anglei[$t]-1];}
151	# j is always in this residue
152	$atj = $atomname[$anglej[$t]-1];
153	# k may be in the next residue
154	if ($anglek[$t] == $natoms+1) {$atk = "+N";}
155	else {$atk = $atomname[$anglek[$t]-1];}
156
157	# print them out
158	$angle = "ga_$angletype[$t]";
159	print sprintf("%5s %5s %5s    %-5s\n", $ati, $atj, $atk, $angle);
160    }
161
162    print " [ impropers ]\n";
163    print ";   ai    aj    ak    al  gromos type\n";
164    for ($t=0;$t<$m;$t++) {
165
166	# convert numbers to names.
167	#  i may be in the next residue
168	if ($imp_i[$t] == $natoms+1) {$ati = "+N";}
169	else {$ati = $atomname[$imp_i[$t]-1];}
170	# j may be in the next or in the previous residue
171	if ($imp_j[$t] == -1) {$atj = "-C";}
172	elsif ($imp_j[$t] == $natoms+1) {$atj = "+N";}
173	else {$atj = $atomname[$imp_j[$t]-1];}
174	# k may be in the next residue or in the previous residue
175	if ($imp_k[$t] == -1) {$atk = "-C";}
176	elsif ($imp_k[$t] == $natoms+1) {$atk = "+N";}
177	else {$atk = $atomname[$imp_k[$t]-1];}
178	# l may be in the next residue
179	if ($imp_l[$t] == $natoms+1) {$atl = "+N";}
180	else {$atl = $atomname[$imp_l[$t]-1];}
181
182	$improper = "gi_$imp_type[$t]";
183	print sprintf("%5s %5s %5s %5s    %-5s\n", $ati, $atj, $atk, $atl,
184		      $improper);
185    }
186
187    print " [ dihedrals ]\n";
188    print ";   ai    aj    ak    al  gromos type\n";
189    for ($t=0;$t<$n;$t++) {
190	# convert numbers to names.
191	#  i may be in the previous residue
192	if ($dih_i[$t] == -1) {$ati = "-C";}
193	elsif ($dih_i[$t] == -2) {$ati = "-CA";}
194	else {$ati = $atomname[$dih_i[$t]-1];}
195	# j may be in the previous residue
196	if ($dih_j[$t] == -1) {$atj = "-C";}
197	else {$atj = $atomname[$dih_j[$t]-1];}
198	# k must be in this residue
199	$atk = $atomname[$dih_k[$t]-1];
200	# l may be in the next residue
201	if ($dih_l[$t] == $natoms+1) {$atl = "+N";}
202	else {$atl = $atomname[$dih_l[$t]-1];}
203
204	$dihedral = "gd_$dih_type[$t]";
205	print sprintf("%5s %5s %5s %5s    %-5s\n", $ati, $atj, $atk, $atl,
206		      $dihedral);
207    }
208    if ($residue_name eq $last) {last;}
209}
210
211
212
213
214
215
216
217
218