scripts/make_gromos_bon.pl

   1 #!/usr/bin/perl
   2
   3 # usage: make_gromos_bon.pl ifp41a1.dat > gromos-bon.itp
   4 # this generates gromos-bon.itp, a list of bonds, angles, impropers
   5 # and dihedrals from ifp41a1.dat.
   6
   7 while (<>) {
   8     if (/BONDTYPE/)         {$prefix = "#define gb_"; $bonds = 1; $go = 1;}
   9     if (/BONDANGLE/)        {$prefix = "#define ga_"; $bonds = 0; $angles = 1;}
  10     if (/IMPDIH/)           {$prefix = "#define gi_"; $angles= 0; $imp = 1;}
  11     if (/DIHEDRALTYPECODE/) {$prefix = "#define gd_"; $imp=0; $dih=1;}
  12     if (/SINGLE/)           {$go = 0;}
  13
  14     if (/^#/ && $go) {tr/#/;/; print $_;}
  15     if (/^\d+/ && $go) {
  16         # need to switch the order of terms for bonds and angles for use
  17         # GROMACS
  18         if ($bonds || $angles ) {
  19             ($nr,$k,$dist) = split(' ',$_);
  20             print sprintf("$prefix%-3d  %10s  %10s\n", $nr, $dist, $k);
  21         }
  22         # impropers have the wrong units for k, convert to degrees
  23         if ($imp) {
  24             ($nr,$k,$dist) = split(' ',$_);
  25             $k = $k*180*180/(3.141593*3.141593);
  26             print sprintf("$prefix%-3d  %10s  %10.5f\n", $nr, $dist, $k);
  27         }
  28         # same for dihedrals, also convert phase from cos(phi) to phi
  29         if ($dih) {
  30             ($nr, $k, $phase, $mult) = split(' ',$_);
  31             if ($phase > 0) {$phase = 0.0;}
  32             else {$phase = 180.0;}
  33             print sprintf("$prefix%-3d %8.3f %10s %10s\n",$nr, $phase,
  34                           $k, $mult);
  35         }
  36         if (/^[a-zA-Z]+/ && $go) {print ";" . $_;}
  37     }
  38 }
  39
  40 # add stuff for the termina for now. Can be removed later if we have
  41 # a working termini database
  42
  43 print "\n[ bondtypes ]\n";
  44 print "NL     H       2    gb_2\n";
  45 print "C      OM      2    gb_5\n";
  46 print "OA     H       2    gb_1\n";
  47 print "C      OA      2    gb_12\n";
  48 print "C      O       2    gb_4\n";
  49 print "S      S       2    gb_33\n";
  50 print "NR     FE      2    gb_32\n";
  51
  52 print "\n[ angletypes ]\n";
  53 print "H      NL     H     2   ga_9\n";
  54 print "H      NL     CH1   2   ga_10\n";
  55 print "CH1    C      OM    2   ga_21\n";
  56 print "OM     C      OM    2   ga_37\n";
  57 print "O      C      OA    2   ga_32\n";
  58 print "C      OA     H     2   ga_11\n";
  59 print "CH1    C      O     2   ga_29\n";
  60 print "CH1    CH2    S     2   ga_15\n";
  61 print "CH2    S      S     2   ga_5\n";
  62 print "CH2    C     OM     2   ga_21\n";
  63 print "CR1    NR    FE     2    ga_33\n";
  64 print "NR     FE    NR     2   ga_16\n";
  65
  66 print "\n[ dihedraltypes ]\n";
  67 print "S      S      1   gd_10\n";
  68 print "NR     FE     1   gd_18\n";
  69 print "CH2    S      1   gd_13\n";
  70