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