1#!/usr/local/bin/perl
2
3# usage: make_gromos_bon.pl ifp41a1.dat > gromos-bon.itp
4# this generates gromos-bon.itp, a list of bonds, angles, impropers
5# and dihedrals from ifp41a1.dat.
6
7while (<>) {
8    if (/BONDTYPE/)         {$prefix = "#define gb_"; $bonds = 1; $go = 1;}
9    if (/BONDANGLE/)        {$prefix = "#define ga_"; $bonds = 0; $angles = 1;}
10    if (/IMPDIH/)           {$prefix = "#define gi_"; $angles= 0; $imp = 1;}
11    if (/DIHEDRALTYPECODE/) {$prefix = "#define gd_"; $imp=0; $dih=1;}
12    if (/SINGLE/)           {$go = 0;}
13
14    if (/^#/ && $go) {tr/#/;/; print $_;}
15    if (/^\d+/ && $go) {
16	# need to switch the order of terms for bonds and angles for use
17	# GROMACS
18	if ($bonds || $angles ) {
19	    ($nr,$k,$dist) = split(' ',$_);
20	    print sprintf("$prefix%-3d  %10s  %10s\n", $nr, $dist, $k);
21	}
22	# impropers have the wrong units for k, convert to degrees
23	if ($imp) {
24	    ($nr,$k,$dist) = split(' ',$_);
25	    $k = $k*180*180/(3.141593*3.141593);
26	    print sprintf("$prefix%-3d  %10s  %10.5f\n", $nr, $dist, $k);
27	}
28	# same for dihedrals, also convert phase from cos(phi) to phi
29	if ($dih) {
30	    ($nr, $k, $phase, $mult) = split(' ',$_);
31	    if ($phase > 0) {$phase = 0.0;}
32	    else {$phase = 180.0;}
33	    print sprintf("$prefix%-3d %8.3f %10s %10s\n",$nr, $phase,
34			  $k, $mult);
35	}
36	if (/^[a-zA-Z]+/ && $go) {print ";" . $_;}
37    }
38}
39
40# add stuff for the termina for now. Can be removed later if we have
41# a working termini database
42
43print "\n[ bondtypes ]\n";
44print "NL     H       2    gb_2\n";
45print "C      OM      2    gb_5\n";
46print "OA     H       2    gb_1\n";
47print "C      OA      2    gb_12\n";
48print "C      O       2    gb_4\n";
49print "S      S       2    gb_33\n";
50print "NR     FE      2    gb_32\n";
51
52print "\n[ angletypes ]\n";
53print "H      NL     H     2   ga_9\n";
54print "H      NL     CH1   2   ga_10\n";
55print "CH1    C      OM    2   ga_21\n";
56print "OM     C      OM    2   ga_37\n";
57print "O      C      OA    2   ga_32\n";
58print "C      OA     H     2   ga_11\n";
59print "CH1    C      O     2   ga_29\n";
60print "CH1    CH2    S     2   ga_15\n";
61print "CH2    S      S     2   ga_5\n";
62print "CH2    C     OM     2   ga_21\n";
63print "CR1    NR    FE     2	ga_33\n";
64print "NR     FE    NR     2   ga_16\n";
65
66print "\n[ dihedraltypes ]\n";
67print "S      S      1   gd_10\n";
68print "NR     FE     1   gd_18\n";
69print "CH2    S      1   gd_13\n";
70
71