src/kernel/g_x2top.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.3.3
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2008, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "maths.h"
  40 #include "macros.h"
  41 #include "copyrite.h"
  42 #include "bondf.h"
  43 #include "gmxfio.h"
  44 #include "string2.h"
  45 #include "smalloc.h"
  46 #include "strdb.h"
  47 #include "sysstuff.h"
  48 #include "confio.h"
  49 #include "physics.h"
  50 #include "statutil.h"
  51 #include "vec.h"
  52 #include "random.h"
  53 #include "3dview.h"
  54 #include "txtdump.h"
  55 #include "readinp.h"
  56 #include "names.h"
  57 #include "toppush.h"
  58 #include "pdb2top.h"
  59 #include "gen_ad.h"
  60 #include "gpp_nextnb.h"
  61 #include "vec.h"
  62 #include "g_x2top.h"
  63 #include "atomprop.h"
  64
  65 char atp[7] = "HCNOSX";
  66 #define NATP (asize(atp)-1)
  67
  68 real blen[NATP][NATP] = {
  69   {  0.00,  0.108, 0.105, 0.10, 0.10, 0.10 },
  70   {  0.108, 0.15,  0.14,  0.14, 0.16, 0.14 },
  71   {  0.105, 0.14,  0.14,  0.14, 0.16, 0.14 },
  72   {  0.10,  0.14,  0.14,  0.14, 0.17, 0.14 },
  73   {  0.10,  0.16,  0.16,  0.17, 0.20, 0.17 },
  74   {  0.10,  0.14,  0.14,  0.14, 0.17, 0.17 }
  75 };
  76
  77 #define MARGIN_FAC 1.1
  78
  79 static real set_x_blen(real scale)
  80 {
  81   real maxlen;
  82   int  i,j;
  83
  84   for(i=0; i<NATP-1; i++) {
  85     blen[NATP-1][i] *= scale;
  86     blen[i][NATP-1] *= scale;
  87   }
  88   blen[NATP-1][NATP-1] *= scale;
  89
  90   maxlen = 0;
  91   for(i=0; i<NATP; i++)
  92     for(j=0; j<NATP; j++)
  93       if (blen[i][j] > maxlen)
  94         maxlen = blen[i][j];
  95
  96   return maxlen*MARGIN_FAC;
  97 }
  98
  99 static gmx_bool is_bond(int nnm,t_nm2type nmt[],char *ai,char *aj,real blen)
 100 {
 101   int i,j;
 102
 103   for(i=0; (i<nnm); i++) {
 104     for(j=0; (j<nmt[i].nbonds); j++) {
 105       if ((((gmx_strncasecmp(ai,nmt[i].elem,1) == 0) &&
 106            (gmx_strncasecmp(aj,nmt[i].bond[j],1) == 0)) ||
 107           ((gmx_strncasecmp(ai,nmt[i].bond[j],1) == 0) &&
 108            (gmx_strncasecmp(aj,nmt[i].elem,1) == 0))) &&
 109           (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
 110         return TRUE;
 111     }
 112   }
 113   return FALSE;
 114 }
 115
 116 static int get_atype(char *nm)
 117 {
 118   int i,aai=NATP-1;
 119
 120   for(i=0; (i<NATP-1); i++) {
 121     if (nm[0] == atp[i]) {
 122       aai=i;
 123       break;
 124     }
 125   }
 126   return aai;
 127 }
 128
 129 void mk_bonds(int nnm,t_nm2type nmt[],
 130               t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
 131               gmx_bool bPBC,matrix box,gmx_atomprop_t aps)
 132 {
 133   t_param b;
 134   int     i,j;
 135   t_pbc   pbc;
 136   rvec    dx;
 137   real    dx2;
 138
 139   for(i=0; (i<MAXATOMLIST); i++)
 140     b.a[i] = -1;
 141   for(i=0; (i<MAXFORCEPARAM); i++)
 142     b.c[i] = 0.0;
 143
 144   if (bPBC)
 145     set_pbc(&pbc,-1,box);
 146   for(i=0; (i<atoms->nr); i++) {
 147     if ((i % 10) == 0)
 148       fprintf(stderr,"\ratom %d",i);
 149     for(j=i+1; (j<atoms->nr); j++) {
 150       if (bPBC)
 151         pbc_dx(&pbc,x[i],x[j],dx);
 152       else
 153         rvec_sub(x[i],x[j],dx);
 154
 155       dx2 = iprod(dx,dx);
 156       if (is_bond(nnm,nmt,*atoms->atomname[i],*atoms->atomname[j],
 157                   sqrt(dx2))) {
 158         b.AI = i;
 159         b.AJ = j;
 160         b.C0 = sqrt(dx2);
 161         add_param_to_list (bond,&b);
 162         nbond[i]++;
 163         nbond[j]++;
 164         if (debug)
 165           fprintf(debug,"Bonding atoms %s-%d and %s-%d\n",
 166                   *atoms->atomname[i],i+1,*atoms->atomname[j],j+1);
 167       }
 168     }
 169   }
 170   fprintf(stderr,"\ratom %d\n",i);
 171 }
 172
 173 int *set_cgnr(t_atoms *atoms,gmx_bool bUsePDBcharge,real *qtot,real *mtot)
 174 {
 175   int    i,n=1;
 176   int    *cgnr;
 177   double qt=0,mt=0;
 178
 179   *qtot = *mtot = 0;
 180   snew(cgnr,atoms->nr);
 181   for(i=0; (i<atoms->nr); i++) {
 182     if (atoms->pdbinfo && bUsePDBcharge)
 183       atoms->atom[i].q = atoms->pdbinfo[i].bfac;
 184     qt += atoms->atom[i].q;
 185     *qtot += atoms->atom[i].q;
 186     *mtot += atoms->atom[i].m;
 187     cgnr[i] = n;
 188     if (is_int(qt)) {
 189       n++;
 190       qt=0;
 191     }
 192   }
 193   return cgnr;
 194 }
 195
 196 gpp_atomtype_t set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
 197                              int *nbonds,int nnm,t_nm2type nm2t[])
 198 {
 199   gpp_atomtype_t atype;
 200   int nresolved;
 201
 202   atype = init_atomtype();
 203   snew(atoms->atomtype,atoms->nr);
 204   nresolved = nm2type(nnm,nm2t,tab,atoms,atype,nbonds,bonds);
 205   if (nresolved != atoms->nr)
 206     gmx_fatal(FARGS,"Could only find a forcefield type for %d out of %d atoms",
 207               nresolved,atoms->nr);
 208
 209   fprintf(stderr,"There are %d different atom types in your sample\n",
 210           get_atomtype_ntypes(atype));
 211
 212   return atype;
 213 }
 214
 215 void lo_set_force_const(t_params *plist,real c[],int nrfp,gmx_bool bRound,
 216                         gmx_bool bDih,gmx_bool bParam)
 217 {
 218   int    i,j;
 219   double cc;
 220   char   buf[32];
 221
 222   for(i=0; (i<plist->nr); i++) {
 223     if (!bParam)
 224       for(j=0; j<nrfp; j++)
 225         c[j] = NOTSET;
 226     else {
 227       if (bRound) {
 228         sprintf(buf,"%.2e",plist->param[i].c[0]);
 229         sscanf(buf,"%lf",&cc);
 230         c[0] = cc;
 231       }
 232       else
 233         c[0] = plist->param[i].c[0];
 234       if (bDih) {
 235         c[0] *= c[2];
 236         c[0] = ((int)(c[0] + 3600)) % 360;
 237         if (c[0] > 180)
 238           c[0] -= 360;
 239         /* To put the minimum at the current angle rather than the maximum */
 240         c[0] += 180;
 241       }
 242     }
 243     for(j=0; (j<nrfp); j++) {
 244       plist->param[i].c[j]      = c[j];
 245       plist->param[i].c[nrfp+j] = c[j];
 246     }
 247     set_p_string(&(plist->param[i]),"");
 248   }
 249 }
 250
 251 void set_force_const(t_params plist[],real kb,real kt,real kp,gmx_bool bRound,
 252                      gmx_bool bParam)
 253 {
 254   int i;
 255   real c[MAXFORCEPARAM];
 256
 257   c[0] = 0;
 258   c[1] = kb;
 259   lo_set_force_const(&plist[F_BONDS],c,2,bRound,FALSE,bParam);
 260   c[1] = kt;
 261   lo_set_force_const(&plist[F_ANGLES],c,2,bRound,FALSE,bParam);
 262   c[1] = kp;
 263   c[2] = 3;
 264   lo_set_force_const(&plist[F_PDIHS],c,3,bRound,TRUE,bParam);
 265 }
 266
 267 void calc_angles_dihs(t_params *ang,t_params *dih,rvec x[],gmx_bool bPBC,
 268                       matrix box)
 269 {
 270   int    i,ai,aj,ak,al,t1,t2,t3;
 271   rvec   r_ij,r_kj,r_kl,m,n;
 272   real   sign,th,costh,ph;
 273   t_pbc  pbc;
 274
 275   if (bPBC)
 276     set_pbc(&pbc,epbcXYZ,box);
 277   if (debug)
 278     pr_rvecs(debug,0,"X2TOP",box,DIM);
 279   for(i=0; (i<ang->nr); i++) {
 280     ai = ang->param[i].AI;
 281     aj = ang->param[i].AJ;
 282     ak = ang->param[i].AK;
 283     th = RAD2DEG*bond_angle(x[ai],x[aj],x[ak],bPBC ? &pbc : NULL,
 284                             r_ij,r_kj,&costh,&t1,&t2);
 285     if (debug)
 286       fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
 287               ai,aj,ak,norm(r_ij),norm(r_kj),th);
 288     ang->param[i].C0 = th;
 289   }
 290   for(i=0; (i<dih->nr); i++) {
 291     ai = dih->param[i].AI;
 292     aj = dih->param[i].AJ;
 293     ak = dih->param[i].AK;
 294     al = dih->param[i].AL;
 295     ph = RAD2DEG*dih_angle(x[ai],x[aj],x[ak],x[al],bPBC ? & pbc : NULL,
 296                            r_ij,r_kj,r_kl,m,n,&sign,&t1,&t2,&t3);
 297     if (debug)
 298       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",
 299               ai,aj,ak,al,norm(r_ij),norm(r_kj),norm(r_kl),ph);
 300     dih->param[i].C0 = ph;
 301   }
 302 }
 303
 304 static void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[])
 305 {
 306   int i;
 307
 308   for(i=0; (i<atoms->nr); i++) {
 309     fprintf(fp,"Atom %5s has %1d bonds\n",*atoms->atomname[i],nbonds[i]);
 310   }
 311 }
 312
 313 static void print_pl(FILE *fp,t_params plist[],int ftp,const char *name,
 314                      char ***atomname)
 315 {
 316   int i,j,nral,nrfp;
 317
 318   if (plist[ftp].nr > 0) {
 319     fprintf(fp,"\n");
 320     fprintf(fp,"[ %s ]\n",name);
 321     nral = interaction_function[ftp].nratoms;
 322     nrfp = interaction_function[ftp].nrfpA;
 323     for(i=0; (i<plist[ftp].nr); i++) {
 324       for(j=0; (j<nral); j++)
 325         fprintf(fp,"  %5s",*atomname[plist[ftp].param[i].a[j]]);
 326       for(j=0; (j<nrfp); j++) {
 327         if (plist[ftp].param[i].c[j] != NOTSET)
 328           fprintf(fp,"  %10.3e",plist[ftp].param[i].c[j]);
 329       }
 330       fprintf(fp,"\n");
 331     }
 332   }
 333 }
 334
 335 static void print_rtp(const char *filenm,const char *title,t_atoms *atoms,
 336                       t_params plist[],gpp_atomtype_t atype,int cgnr[],
 337                       int nbts,int bts[])
 338 {
 339   FILE *fp;
 340   int i,tp;
 341   char *tpnm;
 342
 343   fp = gmx_fio_fopen(filenm,"w");
 344   fprintf(fp,"; %s\n",title);
 345   fprintf(fp,"\n");
 346   fprintf(fp,"[ %s ]\n",*atoms->resinfo[0].name);
 347   fprintf(fp,"\n");
 348   fprintf(fp,"[ atoms ]\n");
 349   for(i=0; (i<atoms->nr); i++) {
 350     tp = atoms->atom[i].type;
 351     if ((tpnm = get_atomtype_name(tp,atype)) == NULL)
 352       gmx_fatal(FARGS,"tp = %d, i = %d in print_rtp",tp,i);
 353     fprintf(fp,"%-8s  %12s  %8.4f  %5d\n",
 354             *atoms->atomname[i],tpnm,
 355             atoms->atom[i].q,cgnr[i]);
 356   }
 357   print_pl(fp,plist,F_BONDS,"bonds",atoms->atomname);
 358   print_pl(fp,plist,F_ANGLES,"angles",atoms->atomname);
 359   print_pl(fp,plist,F_PDIHS,"dihedrals",atoms->atomname);
 360   print_pl(fp,plist,F_IDIHS,"impropers",atoms->atomname);
 361
 362   gmx_fio_fclose(fp);
 363 }
 364
 365 int main(int argc, char *argv[])
 366 {
 367   const char *desc[] = {
 368     "[TT]g_x2top[tt] generates a primitive topology from a coordinate file.",
 369     "The program assumes all hydrogens are present when defining",
 370     "the hybridization from the atom name and the number of bonds.",
 371     "The program can also make an [TT].rtp[tt] entry, which you can then add",
 372     "to the [TT].rtp[tt] database.[PAR]",
 373     "When [TT]-param[tt] is set, equilibrium distances and angles",
 374     "and force constants will be printed in the topology for all",
 375     "interactions. The equilibrium distances and angles are taken",
 376     "from the input coordinates, the force constant are set with",
 377     "command line options.",
 378     "The force fields somewhat supported currently are:[PAR]",
 379     "G53a5  GROMOS96 53a5 Forcefield (official distribution)[PAR]",
 380     "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
 381     "The corresponding data files can be found in the library directory",
 382     "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
 383     "information about file formats. By default, the force field selection",
 384     "is interactive, but you can use the [TT]-ff[tt] option to specify",
 385     "one of the short names above on the command line instead. In that",
 386     "case [TT]g_x2top[tt] just looks for the corresponding file.[PAR]",
 387   };
 388   const char *bugs[] = {
 389     "The atom type selection is primitive. Virtually no chemical knowledge is used",
 390     "Periodic boundary conditions screw up the bonding",
 391     "No improper dihedrals are generated",
 392     "The atoms to atomtype translation table is incomplete ([TT]atomname2type.n2t[tt] file in the data directory). Please extend it and send the results back to the GROMACS crew."
 393   };
 394   FILE       *fp;
 395   t_params   plist[F_NRE];
 396   t_excls    *excls;
 397   t_atoms    *atoms;       /* list with all atoms */
 398   gpp_atomtype_t atype;
 399   t_nextnb   nnb;
 400   t_nm2type  *nm2t;
 401   t_mols     mymol;
 402   gmx_atomprop_t aps;
 403   int        nnm;
 404   char       title[STRLEN],forcefield[32],ffdir[STRLEN];
 405   rvec       *x;        /* coordinates? */
 406   int        *nbonds,*cgnr;
 407   int        bts[] = { 1,1,1,2 };
 408   matrix     box;          /* box length matrix */
 409   int        natoms;       /* number of atoms in one molecule  */
 410   int        nres;         /* number of molecules? */
 411   int        i,j,k,l,m,ndih;
 412   int        epbc;
 413   gmx_bool       bRTP,bTOP,bOPLS;
 414   t_symtab   symtab;
 415   real       cutoff,qtot,mtot;
 416   char       n2t[STRLEN];
 417   output_env_t oenv;
 418
 419   t_filenm fnm[] = {
 420     { efSTX, "-f", "conf", ffREAD  },
 421     { efTOP, "-o", "out",  ffOPTWR },
 422     { efRTP, "-r", "out",  ffOPTWR }
 423   };
 424 #define NFILE asize(fnm)
 425   static real scale = 1.1, kb = 4e5,kt = 400,kp = 5;
 426   static int  nexcl = 3;
 427   static gmx_bool bRemoveDih = FALSE;
 428   static gmx_bool bParam = TRUE, bH14 = TRUE,bAllDih = FALSE,bRound = TRUE;
 429   static gmx_bool bPairs = TRUE, bPBC = TRUE;
 430   static gmx_bool bUsePDBcharge = FALSE,bVerbose=FALSE;
 431   static const char *molnm = "ICE";
 432   static const char *ff = "oplsaa";
 433   t_pargs pa[] = {
 434     { "-ff",     FALSE, etSTR, {&ff},
 435       "Force field for your simulation. Type \"select\" for interactive selection." },
 436     { "-v",      FALSE, etBOOL, {&bVerbose},
 437       "Generate verbose output in the top file." },
 438     { "-nexcl", FALSE, etINT,  {&nexcl},
 439       "Number of exclusions" },
 440     { "-H14",    FALSE, etBOOL, {&bH14},
 441       "Use 3rd neighbour interactions for hydrogen atoms" },
 442     { "-alldih", FALSE, etBOOL, {&bAllDih},
 443       "Generate all proper dihedrals" },
 444     { "-remdih", FALSE, etBOOL, {&bRemoveDih},
 445       "Remove dihedrals on the same bond as an improper" },
 446     { "-pairs",  FALSE, etBOOL, {&bPairs},
 447       "Output 1-4 interactions (pairs) in topology file" },
 448     { "-name",   FALSE, etSTR,  {&molnm},
 449       "Name of your molecule" },
 450     { "-pbc",    FALSE, etBOOL, {&bPBC},
 451       "Use periodic boundary conditions." },
 452     { "-pdbq",  FALSE, etBOOL, {&bUsePDBcharge},
 453       "Use the B-factor supplied in a [TT].pdb[tt] file for the atomic charges" },
 454     { "-param", FALSE, etBOOL, {&bParam},
 455       "Print parameters in the output" },
 456     { "-round",  FALSE, etBOOL, {&bRound},
 457       "Round off measured values" },
 458     { "-kb",    FALSE, etREAL, {&kb},
 459       "Bonded force constant (kJ/mol/nm^2)" },
 460     { "-kt",    FALSE, etREAL, {&kt},
 461       "Angle force constant (kJ/mol/rad^2)" },
 462     { "-kp",    FALSE, etREAL, {&kp},
 463       "Dihedral angle force constant (kJ/mol/rad^2)" }
 464   };
 465
 466   CopyRight(stdout,argv[0]);
 467
 468   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 469                     asize(desc),desc,asize(bugs),bugs,&oenv);
 470   bRTP = opt2bSet("-r",NFILE,fnm);
 471   bTOP = opt2bSet("-o",NFILE,fnm);
 472
 473   if (!bRTP && !bTOP)
 474     gmx_fatal(FARGS,"Specify at least one output file");
 475
 476   aps = gmx_atomprop_init();
 477
 478   /* Force field selection, interactive or direct */
 479   choose_ff(strcmp(ff,"select") == 0 ? NULL : ff,
 480             forcefield,sizeof(forcefield),
 481             ffdir,sizeof(ffdir));
 482
 483   bOPLS = (strcmp(forcefield,"oplsaa") == 0);
 484
 485
 486   mymol.name = strdup(molnm);
 487   mymol.nr   = 1;
 488
 489   /* Init parameter lists */
 490   init_plist(plist);
 491
 492   /* Read coordinates */
 493   get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
 494   snew(atoms,1);
 495
 496   /* make space for all the atoms */
 497   init_t_atoms(atoms,natoms,TRUE);
 498   snew(x,natoms);
 499
 500   read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,NULL,&epbc,box);
 501
 502   sprintf(n2t,"%s",ffdir);
 503   nm2t = rd_nm2type(n2t,&nnm);
 504   if (nnm == 0)
 505     gmx_fatal(FARGS,"No or incorrect atomname2type.n2t file found (looking for %s)",
 506               n2t);
 507   else
 508     printf("There are %d name to type translations in file %s\n",nnm,n2t);
 509   if (debug)
 510     dump_nm2type(debug,nnm,nm2t);
 511   printf("Generating bonds from distances...\n");
 512   snew(nbonds,atoms->nr);
 513   mk_bonds(nnm,nm2t,atoms,x,&(plist[F_BONDS]),nbonds,forcefield,
 514            bPBC,box,aps);
 515
 516   open_symtab(&symtab);
 517   atype = set_atom_type(&symtab,atoms,&(plist[F_BONDS]),nbonds,nnm,nm2t);
 518
 519   /* Make Angles and Dihedrals */
 520   snew(excls,atoms->nr);
 521   printf("Generating angles and dihedrals from bonds...\n");
 522   init_nnb(&nnb,atoms->nr,4);
 523   gen_nnb(&nnb,plist);
 524   print_nnb(&nnb,"NNB");
 525   gen_pad(&nnb,atoms,bH14,nexcl,plist,excls,NULL,bAllDih,bRemoveDih,TRUE);
 526   done_nnb(&nnb);
 527
 528   if (!bPairs)
 529     plist[F_LJ14].nr = 0;
 530   fprintf(stderr,
 531           "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
 532           "          %4d pairs,     %4d bonds and  %4d atoms\n",
 533           plist[F_PDIHS].nr,
 534           bOPLS ? "Ryckaert-Bellemans" : "proper",
 535           plist[F_IDIHS].nr, plist[F_ANGLES].nr,
 536           plist[F_LJ14].nr, plist[F_BONDS].nr,atoms->nr);
 537
 538   calc_angles_dihs(&plist[F_ANGLES],&plist[F_PDIHS],x,bPBC,box);
 539
 540   set_force_const(plist,kb,kt,kp,bRound,bParam);
 541
 542   cgnr = set_cgnr(atoms,bUsePDBcharge,&qtot,&mtot);
 543   printf("Total charge is %g, total mass is %g\n",qtot,mtot);
 544   if (bOPLS) {
 545     bts[2] = 3;
 546     bts[3] = 1;
 547   }
 548
 549   if (bTOP) {
 550     fp = ftp2FILE(efTOP,NFILE,fnm,"w");
 551     print_top_header(fp,ftp2fn(efTOP,NFILE,fnm),
 552                      "Generated by x2top",TRUE,ffdir,1.0);
 553
 554     write_top(fp,NULL,mymol.name,atoms,FALSE,bts,plist,excls,atype,
 555               cgnr,nexcl);
 556     print_top_mols(fp,mymol.name,ffdir,NULL,0,NULL,1,&mymol);
 557
 558     ffclose(fp);
 559   }
 560   if (bRTP)
 561     print_rtp(ftp2fn(efRTP,NFILE,fnm),"Generated by x2top",
 562               atoms,plist,atype,cgnr,asize(bts),bts);
 563
 564   if (debug) {
 565     dump_hybridization(debug,atoms,nbonds);
 566   }
 567   close_symtab(&symtab);
 568
 569   printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",Program());
 570   printf("         Please verify atomtypes and charges by comparison to other\n");
 571   printf("         topologies.\n");
 572
 573   thanx(stderr);
 574
 575   return 0;
 576 }