A freeze group can now be allowed to move rigidly in some dimension by using "freezed...
[gromacs/rigid-bodies.git] / scripts / make_gromos_nb.pl
blob19f9f415aa6d10a91b1f81a80f45a1da919c92db
1 #!/usr/bin/perl
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
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";
21 # start arrays with 1 instead of 0
22 $[=1;
24 $i=1;
26 # read in all entries (43 in ifp43a1.dat). This will BREAK for other input!
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];
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;
45 $i++;
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
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";
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
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";}
75 print sprintf("%5s 0.000 0.000 A %10.8g %10.8g\n",
76 $name[$j],$c6[$j]*$c6[$j],$c12*$c12);
79 # Now make the LJ matrix. It's ok to do some double work in shell scripts,
80 # trust me.
82 print "\n";
83 print "[ nonbond_params ]\n";
84 print "; i j func c6 c12\n";
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";
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";
113 print sprintf("%8s %8s 1 %10.8g %10.8g\n", $name[$j], $name[$k],
114 $c6[$j]*$c6[$k], $c12_j*$c12_k);
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";
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]);