amberparm2molby.pl

   1 #!/usr/bin/perl
   2 #  amberparm2molby: read amber .dat file and convert to a Molby parameter file
   3 #  input (stdin): a .dat file
   4 #  output (stdout): a Molby parm file
   5 #  Written as amberparm2namd.pl by Toshi Nagata, Dec 2 2004
   6 #  Updated for Molby, Nov 13 2009
   7 #  Revised Jul 2 2011
   8
   9 $title = <>;
  10 print "! ", $title;
  11
  12 $blanks = ' ' x 100;
  13
  14 sub format_com {
  15         my ($com) = @_;
  16         $com =~ s/(^\s*)|(\s*$)//g;
  17         $com =~ s/^\!\s*//;
  18         if ($com eq "!") {
  19                 $com = "";
  20         } elsif ($com ne "") {
  21                 $com = " ! " . $com;
  22         }
  23         return $com;
  24 }
  25
  26 while (<>) {
  27         chomp;
  28         last if /^\s*$/;
  29         ($aname, $molw, $pol, $com) = unpack("A2 x1 A10 A10 A77", $_ . $blanks);
  30         push @atoms, $aname;
  31         $molws{$aname} = $molw;
  32         $comments{$aname} = format_com($com);
  33 }
  34
  35 #  Hydrophilic atoms
  36 $line = <>;
  37
  38 #  Bonds
  39 while (<>) {
  40         chomp;
  41         last if /^\s*$/;
  42         ($a1, $a2, $fc, $len, $com) = unpack("A2 x1 A2 A10 A10 A75", $_ . $blanks);
  43         printf "bond %-2s   %-2s   %7.1f %5.3f%s\n", $a1, $a2, $fc, $len, format_com($com);
  44 }
  45 print "\n";
  46
  47 #  Angles
  48 while (<>) {
  49         chomp;
  50         last if /^\s*$/;
  51         ($a1, $a2, $a3, $fc, $ang, $com) = unpack("A2 x1 A2 x1 A2 A10 A10 A72", $_ . $blanks);
  52         printf "angle %-2s   %-2s   %-2s   %7.1f %5.2f%s\n", $a1, $a2, $a3, $fc, $ang, format_com($com);
  53 }
  54 print "\n";
  55
  56 #  Dihedrals
  57 while (<>) {
  58         chomp;
  59         last if /^\s*$/;
  60         ($a1, $a2, $a3, $a4, $div, $fc, $ang, $per, $com) = unpack("A2 x1 A2 x1 A2 x1 A2 A4 A15 A15 A15 A40", $_ . $blanks);
  61         if ($div > 0) {
  62                 $fc /= $div;
  63         }
  64         next if $per < 0;
  65         printf "dihe %-2s   %-2s   %-2s   %-2s %7.2f %7d %7.2f%s\n", $a1, $a2, $a3, $a4, $fc, $per, $ang, format_com($com);
  66 }
  67 print "\n";
  68
  69 #  Impropers
  70 while (<>) {
  71         chomp;
  72         last if /^\s*$/;
  73         ($a1, $a2, $a3, $a4, $div, $fc, $ang, $per, $com) = unpack("A2 x1 A2 x1 A2 x1 A2 A4 A15 A15 A15 A40", $_ . $blanks);
  74         printf "impr %-2s   %-2s   %-2s   %-2s %7.2f %7d %7.2f%s\n", $a1, $a2, $a3, $a4, $fc, $per, $ang, format_com($com);
  75 }
  76 print "\n";
  77
  78 $line = <>;  #  H-bond param
  79 $line = <>;  #  blank
  80
  81 #  Atom equivalences
  82 while (<>) {
  83         chomp;
  84         last if /^\s*$/;
  85         $label = "";
  86         foreach $i (unpack("A2 x2" x 20, $_ . $blanks)) {
  87                 if ($label eq "") {
  88                         $label = $i;
  89                 } elsif ($i !~ /^\s*$/) {
  90                         push @{$equil{$label}}, $i;
  91                 }
  92         }
  93 }
  94
  95 while (<>) {
  96         chomp;
  97         last if /^\s*$/;
  98         ($label, $kind) = unpack("A4 x6 A2", $_. $blanks);
  99         if ($label eq "END") {
 100                 last;
 101         }
 102         if ($kind ne "RE") {
 103                 print STDERR "The vdW parameter other than 'RE' format ($kind) is not supported yet.\n";
 104                 print STDERR "\$_='$_'\n\$label='$label'\n\$kind='$kind'\n";
 105                 exit(1);
 106         }
 107         while (<>) {
 108                 chomp;
 109                 last if /^\s*$/;
 110                 ($sym, $r, $edep) = unpack("x2 A2 x6 A10 A10", $_ . $blanks);
 111                 $sigma = 2 * $r * 0.890899;  #  sigma = 2R * (1/2)**(1/6)
 112                 printf "vdw %-2s %7.4f %7.4f %7.4f %7.4f %d %7.3f%s\n", $sym, $edep, $r, $edep, $r, 0, $molws{$sym}, $comments{$sym};
 113                 foreach $sym2 (@{$equil{$sym}}) {
 114                         printf "vdw %-2s %7.4f %7.4f %7.4f %7.4f %d %7.3f%s\n", $sym2, $edep, $r, $edep, $r, 0, $molws{$sym2}, $comments{$sym2};
 115                 }
 116         }
 117 }
 118 print "\n";