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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "maths.h"
  34 #include "macros.h"
  35 #include "copyrite.h"
  36 #include "bondf.h"
  37 #include "string2.h"
  38 #include "smalloc.h"
  39 #include "strdb.h"
  40 #include "sysstuff.h"
  41 #include "confio.h"
  42 #include "physics.h"
  43 #include "statutil.h"
  44 #include "vec.h"
  45 #include "random.h"
  46 #include "3dview.h"
  47 #include "txtdump.h"
  48 #include "readinp.h"
  49 #include "names.h"
  50 #include "toppush.h"
  51 #include "pdb2top.h"
  52 #include "gen_ad.h"
  53 #include "topexcl.h"
  54 #include "vec.h"
  55 #include "x2top.h"
  56
  57 char atp[6] = "HCNOSX";
  58 #define NATP asize(atp)
  59
  60 real blen[NATP][NATP] = {
  61   {  0.00,  0.108, 0.105, 0.10, 0.10, 0.10 },
  62   {  0.108, 0.15,  0.14,  0.14, 0.16, 0.14 },
  63   {  0.105, 0.14,  0.14,  0.14, 0.16, 0.14 },
  64   {  0.10,  0.14,  0.14,  0.14, 0.17, 0.14 },
  65   {  0.10,  0.16,  0.16,  0.17, 0.20, 0.17 },
  66   {  0.10,  0.14,  0.14,  0.14, 0.17, 0.17 }
  67 };
  68
  69 #define MARGIN_FAC 1.1
  70
  71 static real set_x_blen(real scale)
  72 {
  73   real maxlen;
  74   int  i,j;
  75
  76   for(i=0; i<NATP-1; i++) {
  77     blen[NATP-1][i] *= scale;
  78     blen[i][NATP-1] *= scale;
  79   }
  80   blen[NATP-1][NATP-1] *= scale;
  81
  82   maxlen = 0;
  83   for(i=0; i<NATP; i++)
  84     for(j=0; j<NATP; j++)
  85       if (blen[i][j] > maxlen)
  86         maxlen = blen[i][j];
  87
  88   return maxlen*MARGIN_FAC;
  89 }
  90
  91 static int get_atype(char *nm)
  92 {
  93   int i,aai=NATP-1;
  94
  95   for(i=0; (i<NATP-1); i++) {
  96     if (nm[0] == atp[i]) {
  97       aai=i;
  98       break;
  99     }
 100   }
 101   return aai;
 102 }
 103
 104 static bool is_bond(int aai,int aaj,real len2)
 105 {
 106   bool bIsBond;
 107   real bl;
 108
 109   if (len2 == 0.0)
 110     bIsBond = FALSE;
 111   else {
 112     /* There is a bond when the deviation from ideal length is less than
 113      * 10 %
 114      */
 115     bl = blen[aai][aaj]*MARGIN_FAC;
 116
 117     bIsBond = (len2 < bl*bl);
 118   }
 119 #ifdef DEBUG
 120   if (debug)
 121     fprintf(debug,"ai: %5s  aj: %5s  len: %8.3f  bond: %s\n",
 122             ai,aj,len,BOOL(bIsBond));
 123 #endif
 124   return bIsBond;
 125 }
 126
 127 real get_amass(char *aname,int nmass,char **nm2mass)
 128 {
 129   int    i;
 130   char   nmbuf[32];
 131   double mass;
 132   real   m;
 133
 134   m = 12;
 135   for(i=0; (i<nmass); i++) {
 136     sscanf(nm2mass[i],"%s%lf",nmbuf,&mass);
 137     trim(nmbuf);
 138     if (strcmp(aname,nmbuf) == 0) {
 139       m = mass;
 140       break;
 141     }
 142   }
 143   return m;
 144 }
 145
 146 void mk_bonds(t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
 147               real cutoff,bool bPBC,matrix box)
 148 {
 149   t_param b;
 150   t_atom  *atom;
 151   char    **nm2mass=NULL,buf[128];
 152   int     i,j,aai,nmass;
 153   rvec    dx;
 154   real    dx2,c2;
 155
 156   for(i=0; (i<MAXATOMLIST); i++)
 157     b.a[i] = -1;
 158   for(i=0; (i<MAXFORCEPARAM); i++)
 159     b.c[i] = 0.0;
 160
 161   sprintf(buf,"%s.atp",ff);
 162   nmass = get_file(buf,&nm2mass);
 163   fprintf(stderr,"There are %d type to mass translations\n",nmass);
 164   atom  = atoms->atom;
 165   for(i=0; (i<atoms->nr); i++) {
 166     atom[i].type = get_atype(*atoms->atomname[i]);
 167     atom[i].m    = get_amass(*atoms->atomname[i],nmass,nm2mass);
 168   }
 169   c2 = sqr(cutoff);
 170   if (bPBC)
 171     init_pbc(box);
 172   for(i=0; (i<atoms->nr); i++) {
 173     if ((i % 10) == 0)
 174       fprintf(stderr,"\ratom %d",i);
 175     aai = atom[i].type;
 176     for(j=i+1; (j<atoms->nr); j++) {
 177       if (bPBC)
 178         pbc_dx(x[i],x[j],dx);
 179       else
 180         rvec_sub(x[i],x[j],dx);
 181
 182       dx2 = iprod(dx,dx);
 183       if ((dx2 < c2) && (is_bond(aai,atom[j].type,dx2))) {
 184         b.AI = i;
 185         b.AJ = j;
 186         b.C0 = sqrt(dx2);
 187         push_bondnow (bond,&b);
 188         nbond[i]++;
 189         nbond[j]++;
 190       }
 191     }
 192   }
 193   fprintf(stderr,"\ratom %d\n",i);
 194 }
 195
 196 int *set_cgnr(t_atoms *atoms)
 197 {
 198   int    i,n=0;
 199   int    *cgnr;
 200   double qt = 0;
 201
 202   snew(cgnr,atoms->nr);
 203   for(i=0; (i<atoms->nr); i++) {
 204     qt += atoms->atom[i].q;
 205     cgnr[i] = n;
 206     if (is_int(qt)) {
 207       n++;
 208       qt=0;
 209     }
 210   }
 211   return cgnr;
 212 }
 213
 214 t_atomtype *set_atom_type(t_atoms *atoms,int nbonds[],
 215                           int nnm,t_nm2type nm2t[])
 216 {
 217   static t_symtab symtab;
 218   t_atomtype *atype;
 219   char *type;
 220   int  i,k;
 221
 222   open_symtab(&symtab);
 223   snew(atype,1);
 224   atype->nr       = 0;
 225   atype->atom     = NULL;
 226   atype->atomname = NULL;
 227   atype->bondatomtype = NULL;
 228   k=0;
 229   for(i=0; (i<atoms->nr); i++) {
 230     if ((type = nm2type(nnm,nm2t,*atoms->atomname[i],nbonds[i])) == NULL)
 231       fatal_error(0,"No forcefield type for atom %s (%d) with %d bonds",
 232                   *atoms->atomname[i],i+1,nbonds[i]);
 233     else if (debug)
 234       fprintf(debug,"Selected atomtype %s for atom %s\n",
 235               type,*atoms->atomname[i]);
 236     for(k=0; (k<atype->nr); k++) {
 237       if (strcmp(type,*atype->atomname[k]) == 0) {
 238         atoms->atom[i].type  = k;
 239         atoms->atom[i].typeB = k;
 240         break;
 241       }
 242     }
 243     if (k==atype->nr) {
 244       /* New atomtype */
 245       atype->nr++;
 246       srenew(atype->atomname,atype->nr);
 247       srenew(atype->atom,atype->nr);
 248       srenew(atype->bondatomtype,atype->nr);
 249       atype->atomname[k]   = put_symtab(&symtab,type);
 250       atype->bondatomtype[k] = k; /* Set bond_atomtype identical to atomtype */
 251       atype->atom[k].type  = k;
 252       atoms->atom[i].type  = k;
 253       atype->atom[k].typeB = k;
 254       atoms->atom[i].typeB = k;
 255     }
 256   }
 257   /* MORE CODE */
 258
 259   close_symtab(&symtab);
 260
 261   fprintf(stderr,"There are %d different atom types in your sample\n",
 262           atype->nr);
 263
 264   return atype;
 265 }
 266
 267 void lo_set_force_const(t_params *plist,real c[],int nrfp,bool bRound,
 268                         bool bDih,bool bParam)
 269 {
 270   int    i,j;
 271   double cc;
 272   char   buf[32];
 273
 274   for(i=0; (i<plist->nr); i++) {
 275     if (!bParam)
 276       for(j=0; j<nrfp; j++)
 277         c[j] = NOTSET;
 278     else {
 279       if (bRound) {
 280         sprintf(buf,"%.2e",plist->param[i].c[0]);
 281         sscanf(buf,"%lf",&cc);
 282         c[0] = cc;
 283       }
 284       else
 285         c[0] = plist->param[i].c[0];
 286       if (bDih) {
 287         c[0] *= c[2];
 288         c[0] = ((int)(c[0] + 3600)) % 360;
 289         if (c[0] > 180)
 290           c[0] -= 360;
 291         /* To put the minimum at the current angle rather than the maximum */
 292         c[0] += 180;
 293       }
 294     }
 295     for(j=0; (j<nrfp); j++) {
 296       plist->param[i].c[j]      = c[j];
 297       plist->param[i].c[nrfp+j] = c[j];
 298     }
 299     set_p_string(&(plist->param[i]),"");
 300   }
 301 }
 302
 303 void set_force_const(t_params plist[],real kb,real kt,real kp,bool bRound,
 304                      bool bParam)
 305 {
 306   int i;
 307   real c[MAXFORCEPARAM];
 308
 309   c[0] = 0;
 310   c[1] = kb;
 311   lo_set_force_const(&plist[F_BONDS],c,2,bRound,FALSE,bParam);
 312   c[1] = kt;
 313   lo_set_force_const(&plist[F_ANGLES],c,2,bRound,FALSE,bParam);
 314   c[1] = kp;
 315   c[2] = 3;
 316   lo_set_force_const(&plist[F_PDIHS],c,3,bRound,TRUE,bParam);
 317 }
 318
 319 void calc_angles_dihs(t_params *ang,t_params *dih,rvec x[],bool bPBC,
 320                       matrix box)
 321 {
 322   int    i,ai,aj,ak,al;
 323   rvec   r_ij,r_kj,r_kl,m,n;
 324   real   sign,th,costh,ph,cosph;
 325
 326   if (!bPBC)
 327     box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1000;
 328   if (debug)
 329     pr_rvecs(debug,0,"X2TOP",box,DIM);
 330   for(i=0; (i<ang->nr); i++) {
 331     ai = ang->param[i].AI;
 332     aj = ang->param[i].AJ;
 333     ak = ang->param[i].AK;
 334     th = RAD2DEG*bond_angle(box,x[ai],x[aj],x[ak],r_ij,r_kj,&costh);
 335     if (debug)
 336       fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
 337               ai,aj,ak,norm(r_ij),norm(r_kj),th);
 338     ang->param[i].C0 = th;
 339   }
 340   for(i=0; (i<dih->nr); i++) {
 341     ai = dih->param[i].AI;
 342     aj = dih->param[i].AJ;
 343     ak = dih->param[i].AK;
 344     al = dih->param[i].AL;
 345     ph = RAD2DEG*dih_angle(box,x[ai],x[aj],x[ak],x[al],
 346                            r_ij,r_kj,r_kl,m,n,&cosph,&sign);
 347     if (debug)
 348       fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d al=%3d r_ij=%8.3f r_kj=%8.3f r_kl=%8.3f ph=%8.3f\n",
 349               ai,aj,ak,al,norm(r_ij),norm(r_kj),norm(r_kl),ph);
 350     dih->param[i].C0 = ph;
 351   }
 352 }
 353
 354 static void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[])
 355 {
 356   int i;
 357
 358   for(i=0; (i<atoms->nr); i++) {
 359     fprintf(fp,"Atom %5s has %1d bonds\n",*atoms->atomname[i],nbonds[i]);
 360   }
 361 }
 362
 363 static void print_rtp(char *filenm,char *title,char *name,t_atoms *atoms,
 364                       t_params plist[],t_atomtype *atype,int cgnr[])
 365 {
 366   FILE *fp;
 367 }
 368
 369 int main(int argc, char *argv[])
 370 {
 371   static char *desc[] = {
 372     "x2top generates a primitive topology from a coordinate file.",
 373     "The program assumes all hydrogens are present when defining",
 374     "the hybridization from the atom name and the number of bonds.",
 375     "The program can also make an rtp entry, which you can then add",
 376     "to the rtp database.[PAR]",
 377     "When [TT]-param[tt] is set, equilibrium distances and angles",
 378     "and force constants will be printed in the topology for all",
 379     "interactions. The equilibrium distances and angles are taken",
 380     "from the input coordinates, the force constant are set with",
 381     "command line options."
 382   };
 383   static char *bugs[] = {
 384     "The atom type selection is primitive. Virtually no chemical knowledge is used",
 385     "Periodic boundary conditions screw up the bonding",
 386     "No improper dihedrals are generated",
 387     "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."
 388   };
 389   FILE       *fp;
 390   t_params   plist[F_NRE];
 391   t_excls    *excls;
 392   t_atoms    *atoms;       /* list with all atoms */
 393   t_atomtype *atype;
 394   t_nextnb   nnb;
 395   t_nm2type  *nm2t;
 396   t_mols     mymol;
 397   char       ff[256];
 398   int        nnm;
 399   char       title[STRLEN];
 400   rvec       *x;        /* coordinates? */
 401   int        *nbonds,*cgnr;
 402   int        bts[] = { 1,1,1,2 };
 403   matrix     box;          /* box length matrix */
 404   int        natoms;       /* number of atoms in one molecule  */
 405   int        nres;         /* number of molecules? */
 406   int        i,j,k,l,m;
 407   bool       bRTP,bTOP;
 408   real       cutoff;
 409
 410   t_filenm fnm[] = {
 411     { efSTX, "-f", "conf", ffREAD  },
 412     { efTOP, "-o", "out",  ffOPTWR },
 413     { efRTP, "-r", "out",  ffOPTWR }
 414   };
 415 #define NFILE asize(fnm)
 416   static real scale = 1.1, kb = 4e5,kt = 400,kp = 5;
 417   static int  nexcl = 3;
 418   static bool bParam = FALSE,bH14 = FALSE,bAllDih = FALSE,bRound = TRUE;
 419   static bool bPairs = TRUE, bPBC = TRUE;
 420   static char *molnm = "ICE";
 421   t_pargs pa[] = {
 422     { "-scale", FALSE, etREAL, {&scale},
 423       "Scaling factor for bonds with unknown atom types relative to atom type O" },
 424     { "-nexcl", FALSE, etINT,  {&nexcl},
 425       "Number of exclusions" },
 426     { "-H14",    FALSE, etBOOL, {&bH14},
 427       "Use 3rd neighbour interactions for hydrogen atoms" },
 428     { "-alldih", FALSE, etBOOL, {&bAllDih},
 429       "Generate all proper dihedrals" },
 430     { "-pairs",  FALSE, etBOOL, {&bPairs},
 431       "Output 1-4 interactions (pairs) in topology file" },
 432     { "-name",   FALSE, etSTR,  {&molnm},
 433       "Name of your molecule" },
 434     { "-pbc",    FALSE, etBOOL, {&bPBC},
 435       "Use periodic boundary conditions. Please set the GMXFULLPBC environment variable as well." },
 436     { "-param", FALSE, etBOOL, {&bParam},
 437       "Print parameters in the output" },
 438     { "-round",  FALSE, etBOOL, {&bRound},
 439       "Round off measured values" },
 440     { "-kb",    FALSE, etREAL, {&kb},
 441       "Bonded force constant (kJ/mol/nm^2)" },
 442     { "-kt",    FALSE, etREAL, {&kt},
 443       "Angle force constant (kJ/mol/rad^2)" },
 444     { "-kp",    FALSE, etREAL, {&kp},
 445       "Dihedral angle force constant (kJ/mol/rad^2)" }
 446   };
 447
 448   CopyRight(stdout,argv[0]);
 449
 450   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 451                     asize(desc),desc,asize(bugs),bugs);
 452   bRTP = opt2bSet("-r",NFILE,fnm);
 453   bTOP = opt2bSet("-o",NFILE,fnm);
 454   if (!bRTP && !bTOP)
 455     fatal_error(0,"Specify at least one output file");
 456
 457   cutoff = set_x_blen(scale);
 458
 459   if (bPBC)
 460     set_gmx_full_pbc(stdout);
 461
 462   mymol.name = strdup(molnm);
 463   mymol.nr   = 1;
 464
 465   /* Init parameter lists */
 466   init_plist(plist);
 467
 468   /* Read coordinates */
 469   get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
 470   snew(atoms,1);
 471
 472   /* make space for all the atoms */
 473   init_t_atoms(atoms,natoms,FALSE);
 474   snew(x,natoms);
 475
 476   read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,NULL,box);
 477
 478   choose_ff(ff,255);
 479
 480   snew(nbonds,atoms->nr);
 481
 482   printf("Generating bonds from distances...\n");
 483   mk_bonds(atoms,x,&(plist[F_BONDS]),nbonds,ff,cutoff,bPBC,box);
 484
 485   nm2t = rd_nm2type(ff,&nnm);
 486   printf("There are %d name to type translations\n",nnm);
 487   if (debug)
 488     dump_nm2type(debug,nnm,nm2t);
 489   atype = set_atom_type(atoms,nbonds,nnm,nm2t);
 490
 491   /* Make Angles and Dihedrals */
 492   snew(excls,atoms->nr);
 493   printf("Generating angles and dihedrals from bonds...\n");
 494   init_nnb(&nnb,atoms->nr,4);
 495   gen_nnb(&nnb,plist);
 496   print_nnb(&nnb,"NNB");
 497   gen_pad(&nnb,atoms,bH14,nexcl,plist,excls,NULL,bAllDih);
 498   done_nnb(&nnb);
 499
 500   if (!bPairs)
 501     plist[F_LJ14].nr = 0;
 502   fprintf(stderr,
 503           "There are %4d dihedrals, %4d impropers, %4d angles\n"
 504           "          %4d pairs,     %4d bonds and  %4d atoms\n",
 505           plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
 506           plist[F_LJ14].nr, plist[F_BONDS].nr,atoms->nr);
 507
 508   calc_angles_dihs(&plist[F_ANGLES],&plist[F_PDIHS],x,bPBC,box);
 509
 510   set_force_const(plist,kb,kt,kp,bRound,bParam);
 511
 512   cgnr = set_cgnr(atoms);
 513
 514   if (bTOP) {
 515     fp = ftp2FILE(efTOP,NFILE,fnm,"w");
 516     print_top_header(fp,ftp2fn(efTOP,NFILE,fnm),
 517                      "Generated by x2top",TRUE, ff,1.0);
 518
 519     write_top(fp,NULL,mymol.name,atoms,bts,plist,excls,atype,cgnr,nexcl);
 520     print_top_mols(fp,mymol.name,0,NULL,1,&mymol);
 521
 522     fclose(fp);
 523   }
 524   if (bRTP)
 525     print_rtp(ftp2fn(efRTP,NFILE,fnm),"Generated by x2top",
 526               mymol.name,atoms,plist,atype,cgnr);
 527
 528   if (debug) {
 529     dump_hybridization(debug,atoms,nbonds);
 530   }
 531
 532   thanx(stderr);
 533
 534   return 0;
 535 }