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