src/kernel/x2top.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_x2top_c = "$Id$";
  30
  31 #include "maths.h"
  32 #include "macros.h"
  33 #include "copyrite.h"
  34 #include "bondf.h"
  35 #include "string2.h"
  36 #include "smalloc.h"
  37 #include "strdb.h"
  38 #include "sysstuff.h"
  39 #include "confio.h"
  40 #include "physics.h"
  41 #include "statutil.h"
  42 #include "vec.h"
  43 #include "random.h"
  44 #include "3dview.h"
  45 #include "txtdump.h"
  46 #include "readinp.h"
  47 #include "names.h"
  48 #include "toppush.h"
  49 #include "pdb2top.h"
  50 #include "gen_ad.h"
  51 #include "topexcl.h"
  52 #include "vec.h"
  53 #include "x2top.h"
  54
  55 int get_atype(char *nm)
  56 {
  57   char atp[6] = "HCNOSX";
  58 #define NATP asize(atp)
  59   int i,aai=NATP-1;
  60
  61   for(i=0; (i<NATP-1); i++) {
  62     if (nm[0] == atp[i]) {
  63       aai=i;
  64       break;
  65     }
  66   }
  67   return aai;
  68 }
  69
  70 bool is_bond(int aai,int aaj,real len2)
  71 {
  72   real blen[6][6] = {
  73     {  0.00,  0.108, 0.105, 0.10, 0.10, 0.10 },
  74     {  0.108, 0.15,  0.14,  0.14, 0.16, 0.15 },
  75     {  0.105, 0.14,  0.14,  0.14, 0.16, 0.15 },
  76     {  0.10,  0.14,  0.14,  0.14, 0.17, 0.16 },
  77     {  0.10,  0.16,  0.16,  0.17, 0.20, 0.20 },
  78     {  0.10,  0.15,  0.15,  0.16, 0.20, 0.20 }
  79   };
  80   bool bIsBond;
  81   real bl;
  82
  83   if (len2 == 0.0)
  84     bIsBond = FALSE;
  85   else {
  86     /* There is a bond when the deviation from ideal length is less than
  87      * 10 %
  88      */
  89     bl = blen[aai][aaj]*1.1;
  90     bIsBond = (len2 < bl*bl);
  91   }
  92 #ifdef DEBUG
  93   if (debug)
  94     fprintf(debug,"ai: %5s  aj: %5s  len: %8.3f  bond: %s\n",
  95             ai,aj,len,BOOL(bIsBond));
  96 #endif
  97   return bIsBond;
  98 }
  99
 100 real get_amass(char *aname,int nmass,char **nm2mass)
 101 {
 102   int    i;
 103   char   nmbuf[32];
 104   double mass;
 105   real   m;
 106
 107   m = 12;
 108   for(i=0; (i<nmass); i++) {
 109     sscanf(nm2mass[i],"%s%lf",nmbuf,&mass);
 110     trim(nmbuf);
 111     if (strcmp(aname,nmbuf) == 0) {
 112       m = mass;
 113       break;
 114     }
 115   }
 116   return m;
 117 }
 118
 119 void mk_bonds(t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff)
 120 {
 121   t_param b;
 122   t_atom  *atom;
 123   char    **nm2mass=NULL,buf[128];
 124   int     i,j,aai,nmass;
 125   rvec    dx;
 126   real    dx2;
 127
 128   for(i=0; (i<MAXATOMLIST); i++)
 129     b.a[i] = -1;
 130   for(i=0; (i<MAXFORCEPARAM); i++)
 131     b.c[i] = 0.0;
 132
 133   sprintf(buf,"%s.atp",ff);
 134   nmass = get_file(buf,&nm2mass);
 135   fprintf(stderr,"There are %d type to mass translations\n",nmass);
 136   atom  = atoms->atom;
 137   for(i=0; (i<atoms->nr); i++) {
 138     atom[i].type = get_atype(*atoms->atomname[i]);
 139     atom[i].m    = get_amass(*atoms->atomname[i],nmass,nm2mass);
 140   }
 141   for(i=0; (i<atoms->nr); i++) {
 142     if ((i % 10) == 0)
 143       fprintf(stderr,"\ratom %d",i);
 144     aai = atom[i].type;
 145     for(j=i+1; (j<atoms->nr); j++) {
 146       rvec_sub(x[i],x[j],dx);
 147
 148       if ((dx[XX] < 0.25) && (dx[YY] < 0.25) && (dx[ZZ] < 0.25)) {
 149         dx2 = iprod(dx,dx);
 150
 151         if (is_bond(aai,atom[j].type,dx2)) {
 152           b.AI = i;
 153           b.AJ = j;
 154           b.C0 = sqrt(dx2);
 155           push_bondnow (bond,&b);
 156           nbond[i]++;
 157           nbond[j]++;
 158         }
 159       }
 160     }
 161   }
 162   fprintf(stderr,"\n");
 163 }
 164
 165 int *set_cgnr(t_atoms *atoms)
 166 {
 167   int    i,n=0;
 168   int    *cgnr;
 169   double qt = 0;
 170
 171   snew(cgnr,atoms->nr);
 172   for(i=0; (i<atoms->nr); i++) {
 173     qt += atoms->atom[i].q;
 174     cgnr[i] = n;
 175     if (is_int(qt)) {
 176       n++;
 177       qt=0;
 178     }
 179   }
 180   return cgnr;
 181 }
 182
 183 t_atomtype *set_atom_type(t_atoms *atoms,int nbonds[],
 184                           int nnm,t_nm2type nm2t[])
 185 {
 186   static t_symtab symtab;
 187   t_atomtype *atype;
 188   char *type;
 189   int  i,k;
 190
 191   open_symtab(&symtab);
 192   snew(atype,1);
 193   atype->nr       = 0;
 194   atype->atom     = NULL;
 195   atype->atomname = NULL;
 196   k=0;
 197   for(i=0; (i<atoms->nr); i++) {
 198     if ((type = nm2type(nnm,nm2t,*atoms->atomname[i],nbonds[i])) == NULL)
 199       fatal_error(0,"No forcefield type for atom %s (%d) with %d bonds",
 200                   *atoms->atomname[i],i+1,nbonds[i]);
 201     else if (debug)
 202       fprintf(debug,"Selected atomtype %s for atom %s\n",
 203               type,*atoms->atomname[i]);
 204     for(k=0; (k<atype->nr); k++) {
 205       if (strcmp(type,*atype->atomname[k]) == 0) {
 206         atoms->atom[i].type  = k;
 207         atoms->atom[i].typeB = k;
 208         break;
 209       }
 210     }
 211     if (k==atype->nr) {
 212       /* New atomtype */
 213       atype->nr++;
 214       srenew(atype->atomname,atype->nr);
 215       srenew(atype->atom,atype->nr);
 216       atype->atomname[k]   = put_symtab(&symtab,type);
 217       atype->atom[k].type  = k;
 218       atoms->atom[i].type  = k;
 219       atype->atom[k].typeB = k;
 220       atoms->atom[i].typeB = k;
 221     }
 222   }
 223   /* MORE CODE */
 224
 225   close_symtab(&symtab);
 226
 227   fprintf(stderr,"There are %d different atom types in your sample\n",
 228           atype->nr);
 229
 230   return atype;
 231 }
 232
 233 void lo_set_force_const(t_params *plist,real c[],int nrfp,bool bRound)
 234 {
 235   int    i,j;
 236   double cc;
 237   char   buf[32];
 238
 239   for(i=0; (i<plist->nr); i++) {
 240     if (bRound) {
 241       sprintf(buf,"%.2e",plist->param[i].c[0]);
 242       sscanf(buf,"%lf",&cc);
 243       c[0] = cc;
 244     }
 245     else
 246       c[0] = plist->param[i].c[0];
 247
 248     for(j=0; (j<nrfp); j++) {
 249       plist->param[i].c[j]      = c[j];
 250       plist->param[i].c[nrfp+j] = c[j];
 251     }
 252     plist->param[i].s = strdup("");
 253   }
 254 }
 255
 256 void set_force_const(t_params plist[],real kb,real kt,real kp,bool bRound)
 257 {
 258   int i;
 259   real c[MAXFORCEPARAM];
 260
 261   c[0] = 0;
 262   c[1] = kb;
 263   lo_set_force_const(&plist[F_BONDS],c,2,bRound);
 264   c[1] = kt;
 265   lo_set_force_const(&plist[F_ANGLES],c,2,bRound);
 266   c[1] = kp;
 267   c[2] = 3;
 268   lo_set_force_const(&plist[F_PDIHS],c,3,bRound);
 269 }
 270
 271 void calc_angles_dihs(t_params *ang,t_params *dih,rvec x[])
 272 {
 273   int    i,ai,aj,ak,al;
 274   rvec   r_ij,r_kj,r_kl,m,n;
 275   real   sign,th,costh,ph,cosph;
 276   matrix box;
 277
 278   clear_mat(box);
 279   box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1000;
 280   for(i=0; (i<ang->nr); i++) {
 281     ai = ang->param[i].AI;
 282     aj = ang->param[i].AJ;
 283     ak = ang->param[i].AK;
 284     th = RAD2DEG*bond_angle(box,x[ai],x[aj],x[ak],r_ij,r_kj,&costh);
 285     ang->param[i].C0 = th;
 286   }
 287   for(i=0; (i<dih->nr); i++) {
 288     ai = dih->param[i].AI;
 289     aj = dih->param[i].AJ;
 290     ak = dih->param[i].AK;
 291     al = dih->param[i].AL;
 292     ph = RAD2DEG*dih_angle(box,x[ai],x[aj],x[ak],x[al],
 293                            r_ij,r_kj,r_kl,m,n,&cosph,&sign);
 294     dih->param[i].C0 = ph;
 295   }
 296 }
 297
 298 void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[])
 299 {
 300   int i;
 301
 302   for(i=0; (i<atoms->nr); i++) {
 303     fprintf(fp,"Atom %5s has %1d bonds\n",*atoms->atomname[i],nbonds[i]);
 304   }
 305 }
 306
 307 int main(int argc, char *argv[])
 308 {
 309   static char *desc[] = {
 310     "x2top generates a primitive topology from a coordinate file.",
 311     "The program assumes all hydrogens are present when defining",
 312     "the hybridization from the atom name and the number of bonds."
 313   };
 314   static char *bugs[] = {
 315     "The atom type selection is primitive. Virtually no chemical knowledge is used",
 316     "Periodic boundary conditions screw up the bonding",
 317     "No improper dihedrals are generated",
 318     "The atoms to atomtype translation table is incomplete (ffG43a1.n2t file in the $GMXLIB directory). Please extend it and send the results back to the GROMACS crew."
 319   };
 320   FILE       *fp;
 321   t_params   plist[F_NRE];
 322   t_atoms    *atoms;       /* list with all atoms */
 323   t_atomtype *atype;
 324   t_block    excl;
 325   t_nextnb   nnb;
 326   t_nm2type  *nm2t;
 327   t_mols     mymol;
 328   char       *ff;
 329   int        nnm;
 330   char       title[STRLEN];
 331   rvec       *x;        /* coordinates? */
 332   int        *nbonds,*cgnr;
 333   int        bts[] = { 1,1,1,2 };
 334   matrix     box;          /* box length matrix */
 335   int        natoms;       /* number of atoms in one molecule  */
 336   int        nres;         /* number of molecules? */
 337   int        i,j,k,l,m;
 338
 339   t_filenm fnm[] = {
 340     { efSTX, "-f", "conf", ffREAD  },
 341     { efTOP, "-o", "out",  ffWRITE }
 342   };
 343 #define NFILE asize(fnm)
 344   static real kb    = 4e5,kt = 400,kp = 5;
 345   static int  nexcl = 3;
 346   static bool bH14  = FALSE,bAllDih = FALSE,bRound = TRUE,bPairs = TRUE;
 347   static char *molnm = "ICE";
 348   t_pargs pa[] = {
 349     { "-kb",    FALSE, etREAL, {&kb},
 350       "Bonded force constant (kJ/mol/nm^2)" },
 351     { "-kt",    FALSE, etREAL, {&kt},
 352       "Angle force constant (kJ/mol/rad^2)" },
 353     { "-kp",    FALSE, etREAL, {&kp},
 354       "Dihedral angle force constant (kJ/mol/rad^2)" },
 355     { "-nexcl", FALSE, etINT,  {&nexcl},
 356       "Number of exclusions" },
 357     { "-H14",    FALSE, etBOOL, {&bH14},
 358       "Use 3rd neighbour interactions for hydrogen atoms" },
 359     { "-alldih", FALSE, etBOOL, {&bAllDih},
 360       "Generate all proper dihedrals" },
 361     { "-round",  FALSE, etBOOL, {&bRound},
 362       "Round off measured values" },
 363     { "-pairs",  FALSE, etBOOL, {&bPairs},
 364       "Output 1-4 interactions (pairs) in topology file" },
 365     { "-name",   FALSE, etSTR,  {&molnm},
 366       "Name of your molecule" }
 367   };
 368
 369   CopyRight(stdout,argv[0]);
 370
 371   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
 372                     asize(desc),desc,asize(bugs),bugs);
 373
 374   mymol.name = strdup(molnm);
 375   mymol.nr   = 1;
 376
 377   /* Init lookup table for invsqrt */
 378   init_lookup_table(stdout);
 379
 380   /* Init parameter lists */
 381   init_plist(plist);
 382
 383   /* Read coordinates */
 384   get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
 385   snew(atoms,1);
 386
 387   /* make space for all the atoms */
 388   init_t_atoms(atoms,natoms,FALSE);
 389   snew(x,natoms);
 390
 391   read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,NULL,box);
 392
 393   ff = choose_ff();
 394
 395   snew(nbonds,atoms->nr);
 396
 397   printf("Generating bonds from distances...\n");
 398   mk_bonds(atoms,x,&(plist[F_BONDS]),nbonds,ff);
 399
 400   nm2t = rd_nm2type(ff,&nnm);
 401   printf("There are %d name to type translations\n",nnm);
 402   if (debug)
 403     dump_nm2type(debug,nnm,nm2t);
 404   atype = set_atom_type(atoms,nbonds,nnm,nm2t);
 405
 406   /* Make Angles and Dihedrals */
 407   printf("Generating angles and dihedrals from bonds...\n");
 408   init_nnb(&nnb,atoms->nr,4);
 409   gen_nnb(&nnb,plist);
 410   print_nnb(&nnb,"NNB");
 411   gen_pad(&nnb,atoms,bH14,plist,NULL,bAllDih);
 412   done_nnb(&nnb);
 413
 414   if (!bPairs)
 415     plist[F_LJ14].nr = 0;
 416   fprintf(stderr,
 417           "There are %4d dihedrals, %4d impropers, %4d angles\n"
 418           "          %4d pairs,     %4d bonds and  %4d atoms\n",
 419           plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
 420           plist[F_LJ14].nr, plist[F_BONDS].nr,atoms->nr);
 421
 422   init_block(&excl);
 423
 424   calc_angles_dihs(&plist[F_ANGLES],&plist[F_PDIHS],x);
 425
 426   set_force_const(plist,kb,kt,kp,bRound);
 427
 428   cgnr = set_cgnr(atoms);
 429
 430   fp = ftp2FILE(efTOP,NFILE,fnm,"w");
 431   print_top_header(fp,"Generated by x2top",TRUE, ff,1.0);
 432
 433   write_top(fp,NULL,mymol.name,atoms,bts,plist,&excl,atype,cgnr,nexcl);
 434   print_top_mols(fp,mymol.name,0,NULL,1,&mymol);
 435
 436   fclose(fp);
 437
 438   if (debug) {
 439     dump_hybridization(debug,atoms,nbonds);
 440   }
 441
 442   thanx(stdout);
 443
 444   return 0;
 445 }