(re)added "-dev" suffix for the CMake hardcoded version number
[gromacs/rigid-bodies.git] / scripts / make_gromos_rtp.pl
blob0e7dc225e19d957dba7b098f2a35f9206481c297
1 #!/usr/bin/perl
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
10 # author: Peter Tieleman, 12 March 1999
12 # read in the atomtypes
13 open (TYPES, "gromos.atp") || die "Can't open gromos.atp: $!";
14 while (<TYPES>) {
15 $type++;
16 ($gromostype[$type],$rest) = split(' ',$_);
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.
23 # write header
24 print "[ bondedtypes ]\n";
25 print "; bonds angles dihedrals impropers\n";
26 print " 2 2 1 2\n\n";
28 $first = 1;
30 # loop over the actual residues
31 while (<>) {
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;
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];
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;}
50 # look for the headers of each section, and move to the start
51 # of the relevant data.
53 # if you find ATOM ANM, start reading atoms
54 if ($residue_data[$i] =~ /\#ATOM ANM/) {
55 $read_atoms = 1;
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
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
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
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
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;}
82 if (!($read_atoms || $read_bonds || $read_angles || $read_dihedrals
83 || $read_impropers)) { next; } # we're not reading anything yet
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++;
92 if ($read_bonds) {
93 ($bondi[$k],$bondj[$k],$bondtype[$k]) =
94 split(' ',$residue_data[$i]);
95 $k++;
97 if ($read_angles) {
98 ($anglei[$l],$anglej[$l],$anglek[$l],$angletype[$l]) =
99 split(' ',$residue_data[$i]);
100 $l++;
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++;
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++;
115 $natoms = $j;
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];
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;
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);
143 print " [ angles ]\n";
144 print "; ai aj ak gromos type\n";
145 for ($t=0;$t<$l;$t++) {
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];}
157 # print them out
158 $angle = "ga_$angletype[$t]";
159 print sprintf("%5s %5s %5s %-5s\n", $ati, $atj, $atk, $angle);
162 print " [ impropers ]\n";
163 print "; ai aj ak al gromos type\n";
164 for ($t=0;$t<$m;$t++) {
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];}
182 $improper = "gi_$imp_type[$t]";
183 print sprintf("%5s %5s %5s %5s %-5s\n", $ati, $atj, $atk, $atl,
184 $improper);
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];}
204 $dihedral = "gd_$dih_type[$t]";
205 print sprintf("%5s %5s %5s %5s %-5s\n", $ati, $atj, $atk, $atl,
206 $dihedral);
208 if ($residue_name eq $last) {last;}