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