1#!/usr/local/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
28while(<>) {
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
57print "[ atomtypes ]\n";
58print ";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
63for ($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
82print "\n";
83print "[ nonbond_params ]\n";
84print "; i    j func          c6           c12\n";
85
86for ($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
119print "\n";
120print "[ pairtypes ]\n";
121print "; i    j func          c6           c12\n";
122
123for ($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
137