src/kernel/g_x2top.c

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