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