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