src/contrib/scripts/amber2gmxff

   1 #!/bin/awk -f
   2
   3 # reads Amber forcefield definition file(s) and converts to Gromacs
   4 # Version 1.1
   5 # Copyright (c) 2002
   6 # Anton Feenstra - Vrije Universiteit Amsterdam - The Netherlands
   7
   8 function abs(r) {
   9   return sqrt(r^2);
  10 }
  11
  12 # unit conversion functions:
  13
  14 # Angstrom (1e-10 m) to nanometers (1e-9 m)
  15 function A2nm(a) {
  16   return a/10;
  17 }
  18
  19 # Angstrom^-1 to nanometers^-1
  20 function pA2pnm(a) {
  21   return a*10;
  22 }
  23
  24 # kCalories to kJoules
  25 function kcal2kj(e) {
  26   return e * 4.1868;
  27 }
  28
  29 # kCal mol^-1 A^-n to kJoule mol^-1 nm^-n
  30 # 'n' is usually 2 (for bonds) or 6 or 12 (for LJ)
  31 function kcalpmolpA2kjpmolpnm(k, n) {
  32   kk=kcal2kj(k);
  33   for(ki=0; ki<n; ki++)
  34     kk=pA2pnm(kk);
  35   return kk;
  36 }
  37
  38 function warning(s) {
  39   printf("WARNING: %s\n", s);
  40   nwarn++;
  41 }
  42
  43 function print_warn() {
  44   if (nwarn)
  45     printf("\nThere were %d warnings\n\n", nwarn);
  46 }
  47
  48 function fatal_error(s) {
  49   printf("FATAL ERROR: %s\n", s);
  50   exit -1;
  51 }
  52
  53 BEGIN {
  54   if (ARGC < 2) {
  55     print("Usage:");
  56     print("amber2gmxff [ffname=<name>] [debug=1] parm##.dat");
  57     print("");
  58     print("Reads Amber forcefield definition file (parm##.dat) and writes");
  59     print("Gromacs forcefield files:");
  60     print("ff<name>.atp, ff<name>.itp, ff<name>nb.itp, ff<name>bon.itp");
  61     print("Default for <name> is 'amber'.");
  62     print("");
  63     print("Use 'debug=1' for extremely verbose output.");
  64     print("");
  65     xxit=1;
  66     exit;
  67   }
  68   if ( !ffname ) ffname="amber";
  69   ffbon = "ff" ffname "bon.itp";
  70   ffnb  = "ff" ffname "nb.itp";
  71   ffmain= "ff" ffname ".itp";
  72   ffatp = "ff" ffname ".atp";
  73   ffpdihs = ffname "-dihedrals.txt";
  74   ffidihs = ffname "-impropers.txt";
  75   printf("Sending output to: %s, %s and %s\n", ffmain, ffnb, ffbon);
  76   printf("; Amber forcefield converted to Gromacs\n") > ffnb;
  77   printf("; from file %s\n", ARGV[1]) > ffnb;
  78   printf("; \n") > ffnb;
  79   printf("; Amber forcefield converted to Gromacs\n") > ffbon;
  80   printf("; from file %s\n", ARGV[1]) > ffbon;
  81   printf("\n") > ffbon;
  82   type=0;
  83   nat=0; # atoms
  84   nb=0;  # bonds
  85   na=0;  # angles
  86   np=0;  # propers
  87   ni=0;  # impropers
  88   nhb=0; # H-bonds (10-12)
  89   neq=0; # equivalent 6-12
  90   nlj=0; # LJ (6-12)
  91   section[tp_TITL=0]="title";
  92   section[tp_ATOM=1]="atom";
  93   section[tp_HYDR=2]="hydro";
  94   section[tp_BOND=3]="bond";
  95   section[tp_ANGL=4]="angle";
  96   section[tp_PDIH=5]="pdih";
  97   section[tp_IDIH=6]="idih";
  98   section[tp_HB  =7]="HB";
  99   section[tp_EQ  =8]="eq";
 100   section[tp_LJ  =9]="LJ";
 101
 102 }
 103
 104 # keep track of what we are reading, sections are separated by empty line
 105 NF==0 {
 106   if (type==tp_HB) {
 107     if (bAllZero) {
 108       printf("; NOTE: all H-bond (10-12) parameters are Zero\n") > ffnb;
 109     } else {
 110       printf("; WARNIGN: Not all H-bond (10-12) parameters are Zero\n") > ffnb;
 111       warning("Not all H-bond (10-12) parameters are Zero\n");
 112     }
 113   }
 114   type++;
 115   printf("Start reading section: %d = %s\n", type, section[type]);
 116   next;
 117 }
 118
 119 type==tp_TITL { # title
 120   title = $0;
 121   printf("; %s\n", title) > ffnb;
 122   printf("; \n") > ffnb;
 123   printf("; %s\n", title) > ffbon;
 124   printf("; \n") > ffbon;
 125   type++;
 126   next;
 127 }
 128
 129 type==tp_ATOM { # atoms
 130   if (debug) printf("; %s\n", $0);
 131   anm = substr($0, 1, 2)  ; #  1- 2
 132   anms[nat]  = anm;
 133   amass[anm] = substr($0, 4,10)+0; #  4-13
 134   atpol[anm] = substr($0,15,10)+0; # 15-24
 135   atcom[anm] = substr($0,38,length); # rest
 136   nat++;
 137   if (debug) printf("; '%s' %10g %10g '%s'\n",
 138                     anm, amass[anm], atpol[anm], atcom[anm]);
 139 }
 140
 141 type==tp_HYDR { # hydrophylic atoms
 142   if (debug) printf("; %s\n", $0);
 143   hydrophylics = $0;
 144   type++;
 145   next;
 146 }
 147
 148 type>=tp_BOND && type<=tp_IDIH { # bonded interaction
 149   ai = substr($0, 1, 2)  ; #  1- 2
 150   aj = substr($0, 4, 2)  ; #  4- 5
 151   ak = substr($0, 7, 2)  ; #  7- 8
 152   al = substr($0,10, 2)  ; # 10-11
 153 }
 154
 155 type==tp_BOND { # bonds
 156   if (nb==0) {
 157     printf("\n") > ffbon;
 158     printf("[ bondtypes ]\n") > ffbon;
 159     printf(";%4s %5s %4s %8s %8s\n", "ai", "aj", "ft", "b0", "kb") > ffbon;
 160     printf(";%4s %5s %4s %8s %8s\n","","","", "(nm)", "(kj/mol/nm2)") > ffbon;
 161     ft=1;
 162   }
 163   if (debug) printf("; %s\n", $0);
 164   ibt[nb] = ai;
 165   jbt[nb] = aj;
 166   rk[nb]  = kcalpmolpA2kjpmolpnm(substr($0, 6,10), 2); #  6-15
 167   req[nb] = A2nm(substr($0,16,10)); # 16-25
 168   bcom[nb]= substr($0,29,length); # rest
 169   printf("%5s %5s %4d %8g %8g ; %s\n",
 170          ibt[nb], jbt[nb], ft, req[nb], rk[nb], bcom[nb]) > ffbon;
 171   nb++;
 172 }
 173
 174 type==tp_ANGL { # angles
 175   if (na==0) {
 176     printf("\n") > ffbon;
 177     printf("[ angletypes ]\n") > ffbon;
 178     printf(";%4s %5s %5s %4s %8s %8s\n",
 179            "ai", "aj", "ak", "ft", "th0", "cth") > ffbon;
 180     printf(";%4s %5s %5s %4s %8s %8s\n",
 181            "", "", "", "", "(degr)", "(kj/mol/rad2)") > ffbon;
 182     ft=1;
 183   }
 184   if (debug) printf("; %s\n", $0);
 185   itt[na] = ai;
 186   jtt[na] = aj;
 187   ktt[na] = ak;
 188   tk[na]  = kcal2kj(substr($0, 9,10)); #  9-18
 189   teq[na] = substr($0,19,10)+0; # 19-28
 190   acom[na]= substr($0,33,length); # rest
 191   printf("%5s %5s %5s %4d %8g %8g ; %s\n",
 192          itt[na],jtt[na],ktt[na], ft, teq[na], tk[na], acom[na]) > ffbon;
 193   na++;
 194 }
 195
 196 type==tp_PDIH || type==tp_IDIH { # proper/improper dihedrals
 197   if (type==tp_PDIH && np==0) {
 198     ft=1;
 199     bFirstPdihDefine=1;
 200     printf("\n") > ffbon;
 201     printf("[ dihedraltypes ]\n") > ffbon;
 202     printf(";%4s %5s %5s %8s %8s %5s\n",
 203            "aj","ak", "ft","phi0","cp","mult") > ffbon;
 204     printf(";%4s %5s %5s %8s %8s %5s\n",
 205            "","", "","(degr)","(kJ/mol/rad)","") > ffbon;
 206   }
 207   if (type==tp_IDIH && ni==0) {
 208     ft=1; # we'll have to 'abuse' propers to correspond to Amber definition
 209     printf("\n") > ffbon;
 210     printf(";%4s %5s %5s %5s %5s %8s %8s %5s\n",
 211            "ai","aj","ak","al", "ft","q0","cq","mult") > ffbon;
 212     printf("; WARNING: using Gromacs propers to define Amber impropers\n") > ffbon;
 213     printf("; defining improper parameters for improper types\n") > ffbon;
 214     warning("using Gromacs propers to define Amber impropers\n");
 215   }
 216   if (debug) printf("; %s\n", $0);
 217   idivf   = substr($0,12, 4)+0; # 12-15
 218   pk      = kcal2kj(substr($0,16,15)); # 16-30
 219   phase   = substr($0,31,15)+0; # 31-45
 220   pn      = substr($0,46,15)+0; # 46-50
 221   com     = substr($0,61,length); # rest
 222   if (type==tp_PDIH) {
 223     ipt[np] = ai;
 224     jpt[np] = aj;
 225     kpt[np] = ak;
 226     lpt[np] = al;
 227     p_idivf[np] = idivf;
 228     p_pk[np]    = pk   ;
 229     p_phase[np] = phase;
 230     p_pn[np]    = pn   ;
 231     p_com[np]   = com  ;
 232     np++;
 233     if ( ai=="X " && al=="X " ) {
 234       # this is the type of dihedral GROMACS can manage:
 235       printf("%5s %5s %5d %8g %8g %5d ; %s\n",
 236              aj,ak, ft, phase, pk/idivf, pn, com) > ffbon;
 237       if (pn<0) {
 238         printf(" ; WARNING: multiple dihedral types not supperted") > ffbon;
 239         warning("multiple dihedral types not supperted");
 240       }
 241     } else {
 242       # here we need a hack, since we cannot hava a dihedral type
 243       # depend on all four atom types in GROMACS...
 244       if (bFirstPdihDefine) {
 245         bFirstPdihDefine=0;
 246         printf("\n") > ffbon;
 247         printf("; defining dihedral parameters for four-atom dependent\n") > ffbon;
 248         printf("; dihedral types, which Gromacs does not handle by default\n") > ffbon;
 249       }
 250       print ai, aj, ak, al > ffpdihs;
 251       gsub(" ","_",ai); gsub(" ","_",aj); gsub(" ","_",ak); gsub(" ","_",al);
 252       if (pn<0 && ninc=="")
 253         ninc=1;
 254       printf("#define ad_%2s_%2s_%2s_%2s%1s%1s %8g %8g %5d %8g %8g %5d %s%s\n",
 255              ai,aj,ak,al, ninc?"_":" ", ninc,
 256              phase, pk/idivf, abs(pn), phase, pk/idivf, abs(pn),
 257              length(com)?"; ":"", com) > ffbon;
 258       if (pn<0)
 259         ninc++;
 260       else
 261         ninc="";
 262     }
 263   }
 264   if (type==tp_IDIH) {
 265     iit[ni] = ai;
 266     jit[ni] = aj;
 267     kit[ni] = ak;
 268     lit[ni] = al;
 269     i_idivf[ni] = idivf;
 270     i_pk[ni]    = pk   ;
 271     i_phase[ni] = phase;
 272     i_pn[ni]    = pn   ;
 273     i_com[ni]   = com  ;
 274     ni++;
 275     # here we need a hack, since we cannot hava a dihedral type
 276     # depend on all four atom types in GROMACS...
 277     print ai, aj, ak, al > ffidihs;
 278     gsub(" ","_",ai); gsub(" ","_",aj); gsub(" ","_",ak); gsub(" ","_",al);
 279     printf("#define ai_%2s_%2s_%2s_%2s %8g %8g %5d %8g %8g %5d %s%s\n",
 280            ai,aj,ak,al, pk,phase,pn,pk,phase,pn,
 281            length(com)?"; ":"",com) > ffbon;
 282   }
 283 }
 284
 285 type==tp_HB { # H-bond 10-12 parameters
 286   if (nhb==0) {
 287     printf("\n") > ffnb;
 288     printf("; H-bond params not implemented\n") > ffnb;
 289     warning("H-bond params not implemented\n");
 290     bAllZero=1;
 291   }
 292   if (debug) printf("; %s\n", $0);
 293   a=substr($0, 11, 10)+0; # 11-20
 294   b=substr($0, 21, 10)+0; # 21-30
 295   if (a!=0 || b!=0) {
 296     printf("; WARNING: non-zero H-bond (10-12) parameters: %s %s %g %g\n",
 297            ai, aj, a, b) > ffnb;
 298     warning("non-zero H-bond (10-12) parameters:", ai, aj, a, b);
 299     bAllZero=0;
 300   }
 301   nhb++;
 302 }
 303
 304 type==tp_EQ { # equivalent atom symbols 6-12 parameters
 305   if (neq==0) {
 306     printf("\n") > ffnb;
 307     printf("; LJ 6-12 equivalent atom symbols\n") > ffnb;
 308
 309   }
 310   if (debug) printf("; %s\n", $0);
 311   anm=substr($0, 1, 2);
 312   for(i=1; i<length; i+=4) {
 313     eqat[neq]++;
 314     eqat[neq, int(i/4)]=substr($0, i, 2);
 315   }
 316   printf("; ") > ffnb;
 317   for(i=0; i<eqat[neq]; i++) {
 318     printf("%3s ", eqat[neq, i]) > ffnb;
 319   }
 320   printf("(%d)\n", eqat[neq]) > ffnb;
 321   neq++;
 322 }
 323
 324 type==tp_LJ { # 6-12 potential parameters
 325   if ( !kindnb ) {
 326     label = $1;
 327     kindnb= $2;
 328     next;
 329   }
 330   if (nlj==0) {
 331     printf("\n") > ffnb;
 332     printf("; Define LJ 6-12 parameter types:\n") > ffnb;
 333   }
 334   if (debug) printf("; %s\n", $0);
 335   ljt[nlj] = substr($0, 3, 2)  ; #  3- 4
 336   p1       = substr($0,11,20)+0; # 11-20
 337   p2       = substr($0,21,30)+0; # 21-30
 338   if (kindnb=="AC") {
 339     p  = kcalpmolpA2kjpmolpnm(p2, 6);
 340     p2 = kcalpmolpA2kjpmolpnm(p1, 12);
 341     p1 = p;
 342   }
 343   if (kindnb=="RE") {
 344     p1 = A2nm(p1);
 345     p2 = kcal2kj(p2);
 346   }
 347   p3 = substr($0,31,40)+0; # 31-40
 348   if (debug) printf("; '%s' %10g %10g %10g\n",
 349                     ljt[nlj], p1, p2, p3);
 350   printf("#define alj_%s %10g %10g\n", ljt[nlj], p1,p2) > ffnb;
 351   nlj++;
 352 }
 353
 354 END {
 355   if (xxit) exit;
 356
 357   # write forcefield main ff<name>.itp file
 358   # we do that here so we have the title...
 359   printf("#define _FF_%s\n", toupper(ffname)) > ffmain;
 360   printf("; Amber forcefield converted to Gromacs\n") > ffmain;
 361   printf("; from file %s\n", ARGV[1]) > ffmain;
 362   printf("; \n") > ffmain;
 363   printf("; %s\n", title) > ffmain;
 364   printf("\n") > ffmain;
 365   printf("[ defaults ]\n") > ffmain;
 366   if (kindnb=="AC") cr=1;
 367   if (kindnb=="RE") cr=2;
 368   printf("; %8s %10s %10s %10s %10s\n",
 369          "nbfunc", "comb-rule", "gen-pairs", "fudgeLJ", "fudgeQQ") > ffmain;
 370   printf("%10d %10d %10s %10.3g %10.3g\n", 1, cr, "no", 1.0, 1.0) > ffmain;
 371   printf("\n") > ffmain;
 372   printf("#include \"%s\"\n", ffnb) > ffmain;
 373   printf("#include \"%s\"\n", ffbon) > ffmain;
 374
 375   # expand LJ parameters (using 'equivalent atoms')
 376   for(i=0; i<neq; i++) {
 377     eqnm = eqat[i, 0];
 378     if (debug) print i, eqnm, eqat[i];
 379     for(j=1; j<eqat[i]; j++)
 380       atlj[eqat[i, j]] = eqnm;
 381   }
 382   for(i=0; i<nat; i++) {
 383     anm=anms[i];
 384     if (!atlj[anm]) atlj[anm]=anm;
 385     if (debug) print i, anm, atlj[anm];
 386   }
 387   printf("\n") > ffnb;
 388   printf("[ atomtypes ]\n") > ffnb;
 389   if (kindnb=="AC")
 390     printf(";%4s%10s %4s %2s %-10s\n",
 391            "name", "mass", "q", "tp", "c6/c12") > ffnb;
 392   if (kindnb=="RE")
 393     printf(";%4s%10s %4s %2s %-10s\n",
 394            "name", "mass", "q", "tp", "sigma/epsilon") > ffnb;
 395   q=0;
 396   ptype="A";
 397   for(i=0; i<nat; i++) {
 398     anm=anms[i];
 399     printf("%4s %10.6g %4.2g %2s alj_%-3s ; %s\n",
 400            anm, amass[anm], q, ptype, atlj[anm], atcom[anm]) > ffnb;
 401     # here we write .atp
 402     printf("%4s %10.6g ; %s\n", anm, amass[anm], atcom[anm]) > ffatp;
 403   }
 404   printf("\n");
 405
 406   printf("Found:\n");
 407   printf("# atoms:           %s\n", nat);
 408   printf("# bonds:           %s\n", nb);
 409   printf("# angles:          %s\n", na);
 410   printf("# propers:         %s\n", np);
 411   printf("# impropers:       %s\n", ni);
 412   printf("# H-bonds (10-12): %s\n", nhb);
 413   printf("# equivalent 6-12: %s\n", neq);
 414   printf("# LJ (6-12):       %s\n", nlj);
 415
 416   printf("\n") > ffnb;
 417   printf("; Found:\n") > ffnb;
 418   printf("; # atoms:           %s\n", nat) > ffnb;
 419   printf("; # H-bonds (10-12): %s\n", nhb) > ffnb;
 420   printf("; # equivalent 6-12: %s\n", neq) > ffnb;
 421   printf("; # LJ (6-12):       %s\n", nlj) > ffnb;
 422   printf("\n") > ffbon;
 423   printf("; Found:\n") > ffbon;
 424   printf("; # bonds:           %s\n", nb)  > ffbon;
 425   printf("; # angles:          %s\n", na)  > ffbon;
 426   printf("; # propers:         %s\n", np)  > ffbon;
 427   printf("; # impropers:       %s\n", ni)  > ffbon;
 428
 429   print_warn();
 430 }
 431
 432 #last line