scripts/make_gromos_nb.pl

   1 #!/usr/bin/perl
   2
   3 # usage: make_gromos_nb.pl gromos_atoms > gromos_nb.itp
   4 # this script generates a GROMOS96 nonbonded forcefield file for GROMACS,
   5 # with as input the file gromos_atoms, which contains records of the form
   6 #:        1       O       0.04756 0.8611E-3       1.125E-3        0.0
   7 #: #CS6 CS12 parameters LJ14PAIR
   8 #:       0.04756 0.8611E-3
   9 #:  1   1   2   2   2   2   2   2   2   2   1   1   1   1   1   1   1   1   1   1
  10 #:  2   2   2   2   2   2   2   1   1   1   1   1   2   2   1   1   1   1   1   1
  11 #:  1   1   1
  12 #: #---
  13 # taken directly from ifp43a1.dat. For a description of what the numbers
  14 # mean, see the GROMOS96 manual and the fie ifp43a1.dat
  15
  16 # set the input separator to #--- so we read one atom entry at a time
  17 # make sure the records are actually separated by #---\n and not by say #--\n!!
  18 # don't put a separator at the end of the file, only between records
  19 $/= "#---\n";
  20
  21 # start arrays with 1 instead of 0
  22 $[=1;
  23
  24 $i=1;
  25
  26 # read in all entries (43 in ifp43a1.dat). This will BREAK for other input!
  27
  28 while(<>) {
  29     @lines = split('\n',$_);
  30     ($number[$i],$name[$i],$c6[$i],$c12_1[$i],$c12_2[$i],$c12_3[$i]) =
  31         split(' ',$lines[1]);
  32     ($c6_14[$i],$c12_14[$i]) = split(' ',$lines[3]);
  33     $combination[$i] = $lines[4] . $lines[5] . $lines[6];
  34
  35     # one type is called P,SI, the same LJ parameters for both P and SI
  36     # treat P,SI different: create a 44th type for Si, just a copy of Si,P
  37     # and rename P,SI to P
  38     if ($name[$i] =~ /P,SI/) {
  39         $number[44] = 44; $name[44] = "SI"; $c6[44] = $c6[$i];
  40         $c12_1[44] = $c12_1[$i]; $c12_2[44] = $c12_2[$i];
  41         $c12_3[44] = $c12_3[$i]; $combination[44] = $combination[$i];
  42         $name[$i] = "P";
  43         $P = $i;
  44     }
  45     $i++;
  46 }
  47
  48 # now $number[$i] has the atom nr, $name[$i] the name, $c6 the C6 LJ parameter
  49 # $c12_1[$i], $c12_2[$i], $c12_3[$i] the three C12 parameters and
  50 # $combination[$i] the matrix with 43 elements that tells which C12 parameter
  51 # you need in combination with each of the 44 atomtypes. $i runs from 1 to 44,
  52 # one for each atom type. This goes nicely wrong because of the Stupid SI,P
  53 # entry so we have to give an extra element 44 with the same value as that of
  54 # P
  55
  56 # start printing to the output file: header for gromos-nb.itp
  57 print "[ atomtypes ]\n";
  58 print ";name        mass      charge   ptype            c6           c12\n";
  59
  60 # print out the atomtypes, plus the C6 and C12 for interactions with themselves
  61 # the masses and charges are set to 0, the masses should come from gromos.atp
  62
  63 for ($j=1;$j<=44;$j++) {
  64 # lookup the C12 with itself
  65     @c12types = split(' ',$combination[$j]);
  66     $c12types[44] = $c12types[$P]; # SI has the same as P
  67     if ($c12types[$j] == 1)
  68     {$c12 = $c12_1[$j];}
  69     elsif ($c12types[$j] == 2)
  70     {$c12 = $c12_2[$j];}
  71     elsif ($c12types[$j] == 3)
  72     {$c12 = $c12_3[$j];}
  73     else {die "Error [atomtypes]: c12 type is not 1,2,3:j=$j,c12=$c12\n";}
  74
  75     print sprintf("%5s     0.000      0.000     A  %10.8g  %10.8g\n",
  76                   $name[$j],$c6[$j]*$c6[$j],$c12*$c12);
  77 }
  78
  79 # Now make the LJ matrix. It's ok to do some double work in shell scripts,
  80 # trust me.
  81
  82 print "\n";
  83 print "[ nonbond_params ]\n";
  84 print "; i    j func          c6           c12\n";
  85
  86 for ($j=1;$j<=44;$j++) {
  87     for ($k=1;$k<$j;$k++) {
  88 # lookup the C12 of j for k
  89         @c12types_j = split(' ',$combination[$j]);
  90         $c12types_j[44] = $c12types_j[$P]; # SI has the same as P
  91         if ($c12types_j[$k] == 1)
  92         {$c12_j = $c12_1[$j];}
  93         elsif ($c12types_j[$k] == 2)
  94         {$c12_j = $c12_2[$j];}
  95         elsif ($c12types_j[$k] == 3)
  96         {$c12_j = $c12_3[$j];}
  97         else {
  98             die "Error [nonbond-params] j=$j,k=$k: c12 type is not one of 1,2,3\n";
  99         }
 100 # lookup the C12 of k for j
 101         @c12types_k = split(' ',$combination[$k]);
 102         $c12types_k[44] = $c12types_k[$P]; # SI has the same as P
 103         if ($c12types_k[$j] == 1)
 104         {$c12_k = $c12_1[$k];}
 105         elsif ($c12types_k[$j] == 2)
 106         {$c12_k = $c12_2[$k];}
 107         elsif ($c12types_k[$j] == 3)
 108         {$c12_k = $c12_3[$k];}
 109         else {
 110             die "Error [nonbond-params] j=$j,k=$k: c12 type is not one of 1,2,3\n";
 111         }
 112
 113         print sprintf("%8s %8s  1  %10.8g  %10.8g\n", $name[$j], $name[$k],
 114                       $c6[$j]*$c6[$k], $c12_j*$c12_k);
 115     }
 116 }
 117
 118 # Now do the same for the 1-4 interactions
 119 print "\n";
 120 print "[ pairtypes ]\n";
 121 print "; i    j func          c6           c12\n";
 122
 123 for ($j=1;$j<=44;$j++) {
 124     for ($k=1;$k<=$j;$k++) {
 125         print sprintf("%8s %8s  1  %10.8g  %10.8g\n", $name[$j], $name[$k],
 126                       $c6_14[$j]*$c6_14[$k], $c12_14[$j]*$c12_14[$k]);
 127
 128     }
 129 }
 130
 131
 132
 133
 134
 135
 136