scripts/make_gromos_rtp.pl

   1 #!/usr/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
  13 open (TYPES, "gromos.atp") || die "Can't open gromos.atp: $!";
  14 while (<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
  24 print "[ bondedtypes ]\n";
  25 print "; bonds  angles  dihedrals  impropers\n";
  26 print "    2       2          1          2\n\n";
  27
  28 $first = 1;
  29
  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;
  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