src/programs/alexandria/gentop.cpp

   1 /*
   2  * This source file is part of the Alexandria Chemistry Toolkit.
   3  *
   4  * Copyright (C) 2014-2020
   5  *
   6  * Developers:
   7  *             Mohammad Mehdi Ghahremanpour,
   8  *             Paul J. van Maaren,
   9  *             David van der Spoel (Project leader)
  10  *
  11  * This program is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU General Public License
  13  * as published by the Free Software Foundation; either version 2
  14  * of the License, or (at your option) any later version.
  15  *
  16  * This program is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  * GNU General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU General Public License
  22  * along with this program; if not, write to the Free Software
  23  * Foundation, Inc., 51 Franklin Street, Fifth Floor,
  24  * Boston, MA  02110-1301, USA.
  25  */
  26
  27 /*! \internal \brief
  28  * Implements part of the alexandria program.
  29  * \author Mohammad Mehdi Ghahremanpour <mohammad.ghahremanpour@icm.uu.se>
  30  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  31  */
  32
  33 #include <ctype.h>
  34 #include <stdlib.h>
  35
  36 #include "gromacs/commandline/filenm.h"
  37 #include "gromacs/commandline/pargs.h"
  38 #include "gromacs/gmxlib/network.h"
  39 #include "gromacs/mdtypes/commrec.h"
  40 #include "gromacs/mdtypes/enerdata.h"
  41 #include "gromacs/mdtypes/inputrec.h"
  42 #include "gromacs/topology/atomprop.h"
  43 #include "gromacs/utility/arraysize.h"
  44 #include "gromacs/utility/cstringutil.h"
  45 #include "gromacs/utility/futil.h"
  46 #include "gromacs/utility/stringutil.h"
  47 #include "gromacs/utility/smalloc.h"
  48
  49 #include "alex_modules.h"
  50 #include "babel_io.h"
  51 #include "fill_inputrec.h"
  52 #include "molprop_util.h"
  53 #include "molprop_xml.h"
  54 #include "mymol.h"
  55 #include "poldata_xml.h"
  56
  57 static void print_errors(const char                     *fn,
  58                          const std::vector<std::string> &errors,
  59                          immStatus                       imm)
  60 {
  61     FILE  *fp   = gmx_ffopen(fn, "w");
  62     fprintf(fp, "The detected error is \"%s\"\n", alexandria::immsg(imm));
  63     for (auto error : errors)
  64     {
  65         fprintf(fp, "Error: %s\n", error.c_str());
  66     }
  67     fclose(fp);
  68 }
  69
  70 int alex_gentop(int argc, char *argv[])
  71 {
  72     static const char               *desc[] = {
  73         "gentop generates a topology from molecular coordinates",
  74         "either from a file, from a database, or from a gaussian log file.",
  75         "The program assumes all hydrogens are present when defining",
  76         "the hybridization from the atom name and the number of bonds.[PAR]",
  77         "If the [TT]-oi[tt] option is set an [TT]itp[tt] file will be generated",
  78         "instead of a [TT]top[tt] file.",
  79         "When [TT]-param[tt] is set, equilibrium distances and angles",
  80         "and force constants will be printed in the topology for all",
  81         "interactions. The equilibrium distances and angles are taken",
  82         "from the input coordinates, the force constant are set with",
  83         "command line options.",
  84         "With the [TT]-db molecule[tt] option a file is extracted from the",
  85         "database from one of the specified QM calculations (given with [TT]-lot[tt]).",
  86         "An alternative to the system-wide database [TT]molprops.dat[tt]",
  87         "can be passed along using the [TT]-mpdb[tt] flag.[PAR]",
  88         "If the flag [TT]-qgen[tt] is given, charges will be generated using the",
  89         "specified algorithm. Without the flag the charges from the QM calculation",
  90         "will be used.",
  91         "The supported force field for this tool are Alexandria with two",
  92         "different flavors (specify with the -ff flag): [BR]",
  93         "ACM-g  : Alexandria Charge Model with Gaussian charges, and[BR]",
  94         "ACM-s  : Alexandria Charge Model with Slater charges, and[BR]",
  95         "ACM-pg : Alexandria Charge Model with Polarizable Gaussian charges.[PAR]",
  96         "ACM-ps : Alexandria Charge Model with Polarizable Slater charges.[PAR]",
  97         "A few other choices are available for historical reasons and for",
  98         "development:[PAR]",
  99         "ESP-p  : ElectroStatic Potential with point charges[PAR]",
 100         "ESP-pp : ESP with with polarizable point charges[PAR]",
 101         "ESP-pg : ESP with with polarizable Gaussian-distributed charges[PAR]",
 102         "ESP-ps : ESP with with polarizable Slater-distributed charges[PAR]",
 103         "Rappe  : Rappe and Goddard (J Phys Chem 95 (1991) 3358)[PAR]",
 104         "Yang   : Yang & Sharp (J Chem Theory Comput 2 (2006), 1152)[PAR]",
 105         "Bultinck: Bultinck et al. (J Phys Chem A 106 (2002) 7887)[PAR]",
 106         "The corresponding data files can be found in the library directory",
 107         "in subdirectory alexandria.ff. Check chapter 5 of the manual for more",
 108         "information about file formats.[PAR]"
 109     };
 110     const char                      *bugs[] = {
 111         "No force constants for impropers are generated"
 112     };
 113     gmx_output_env_t                *oenv;
 114     gmx_atomprop_t                   aps;
 115     std::vector<alexandria::MolProp> mps;
 116     alexandria::MolPropIterator      mpi;
 117     alexandria::MyMol                mymol;
 118     immStatus                        imm;
 119
 120     t_filenm                         fnm[] = {
 121         { efTOP, "-p",        "out",           ffOPTWR },
 122         { efITP, "-oi",       "out",           ffOPTWR },
 123         { efSTO, "-c",        "out",           ffWRITE },
 124         { efNDX, "-n",        "renum",         ffOPTWR },
 125         { efDAT, "-q",        "qout",          ffOPTWR },
 126         { efDAT, "-mpdb",     "molprops",      ffOPTRD },
 127         { efDAT, "-d",        "gentop",        ffOPTRD },
 128         { efXVG, "-table",    "table",         ffOPTRD },
 129         { efCUB, "-pot",      "potential",     ffOPTWR },
 130         { efCUB, "-ref",      "refpot",        ffOPTRD },
 131         { efCUB, "-diff",     "diffpot",       ffOPTWR },
 132         { efCUB, "-rho",      "density",       ffOPTWR },
 133         { efXVG, "-diffhist", "diffpot",       ffOPTWR },
 134         { efXVG, "-his",      "pot-histo",     ffOPTWR },
 135         { efXVG, "-pc",       "pot-comp",      ffOPTWR },
 136         { efPDB, "-pdbdiff",  "pdbdiff",       ffOPTWR },
 137         { efXVG, "-plotESP",  "ESPcorr",       ffOPTWR },
 138         { efLOG, "-g",        "gentop_errors", ffWRITE }
 139     };
 140
 141     const  int                       NFILE          = asize(fnm);
 142
 143     static int                       maxpot         = 100;
 144     static int                       nsymm          = 0;
 145     static int                       qcycle         = 1000;
 146     static int                       nexcl          = 2;
 147     static real                      qtol           = 1e-6;
 148     static real                      qtot           = 0;
 149     static real                      hfac           = 0;
 150     static real                      watoms         = 0;
 151     static real                      spacing        = 0.1;
 152     static real                      efield         = 0;
 153     static char                     *molnm          = (char *)"";
 154     static char                     *iupac          = (char *)"";
 155     static char                     *dbname         = (char *)"";
 156     static char                     *symm_string    = (char *)"";
 157     static char                     *conf           = (char *)"minimum";
 158     static char                     *jobtype        = (char *)"unknown";
 159     static char                     *filename       = (char *)"";
 160     static gmx_bool                  bQsym          = false;
 161     static gmx_bool                  bITP           = false;
 162     static gmx_bool                  bPairs         = false;
 163     static gmx_bool                  bUsePDBcharge  = false;
 164     static gmx_bool                  bGenVSites     = false;
 165     static gmx_bool                  bDihedral      = false;
 166     static gmx_bool                  bCUBE          = false;
 167     static gmx_bool                  bH14           = true;
 168     static gmx_bool                  bVerbose       = true;
 169     static gmx_bool                  addHydrogens   = false;
 170
 171     static const char               *ff[]           = {nullptr, "ACM-g", "ACM-pg", "ACM-s", "ACM-ps", "ESP-p", "ESP-pp", "ESP-pg", "ESP-ps", "Yang", "Bultinck", "Rappe", nullptr};
 172     static const char               *cgopt[]        = {nullptr, "Atom", "Group", "Neutral", nullptr};
 173     static const char               *lot            = nullptr;
 174
 175     t_pargs                          pa[]     =
 176     {
 177         { "-f",      FALSE, etSTR,  {&filename},
 178            "Input file name to be turned into GROMACS input" },
 179         { "-v",      FALSE, etBOOL, {&bVerbose},
 180           "Generate verbose output in the top file and on terminal." },
 181         { "-cube",   FALSE, etBOOL, {&bCUBE},
 182           "Generate cube." },
 183         { "-db",     FALSE, etSTR,  {&dbname},
 184           "Read a molecule from the database rather than from a file" },
 185         { "-lot",    FALSE, etSTR,  {&lot},
 186           "Use this method and level of theory when selecting coordinates and charges" },
 187         { "-dih",    FALSE, etBOOL, {&bDihedral},
 188           "Add dihedrals to the topology" },
 189         { "-H14",    FALSE, etBOOL, {&bH14},
 190           "HIDDENUse 3rd neighbour interactions for hydrogen atoms" },
 191         { "-pairs",  FALSE, etBOOL, {&bPairs},
 192           "HIDDENOutput 1-4 interactions (pairs) in topology file. Check consistency of your option with the [TT]-nexcl[tt] flag." },
 193         { "-name",   FALSE, etSTR,  {&molnm},
 194           "Name of your molecule" },
 195         { "-iupac",   FALSE, etSTR,  {&iupac},
 196           "IUPAC Name of your molecule" },
 197         { "-conf",  FALSE, etSTR, {&conf},
 198           "Conformation of the molecule" },
 199         { "-maxpot", FALSE, etINT, {&maxpot},
 200           "Fraction of potential points to read from the gaussian file (percent). If 100 all points are registered, else a selection of points evenly spread over the range of values is taken" },
 201         { "-nsymm", FALSE, etINT, {&nsymm},
 202           "Symmetry number of the molecule can be supplied here if you know there is an error in the input file" },
 203         { "-genvsites", FALSE, etBOOL, {&bGenVSites},
 204           "Generate virtual sites. Check and double check." },
 205         { "-pdbq",  FALSE, etBOOL, {&bUsePDBcharge},
 206           "HIDDENUse the B-factor supplied in a pdb file for the atomic charges" },
 207         { "-spacing", FALSE, etREAL, {&spacing},
 208           "Spacing of grid points for computing the potential (not used when a reference file is read)." },
 209         { "-watoms", FALSE, etREAL, {&watoms},
 210           "Weight for the atoms when fitting the charges to the electrostatic potential. The potential on atoms is usually two orders of magnitude larger than on other points (and negative). For point charges or single smeared charges use 0. For point+smeared charges 1 is recommended." },
 211         { "-ff",     FALSE, etENUM, {ff},
 212           "Force field model. Note that only ACM-xx will yield complete topologies but see help text ([TT]-h[tt])." },
 213         { "-qtol",   FALSE, etREAL, {&qtol},
 214           "Tolerance for assigning charge generation algorithm" },
 215         { "-qtot",   FALSE, etREAL, {&qtot},
 216           "Total charge of the molecule. If the input file is a Gaussian log file, qtot will be taken from the log file." },
 217         { "-addh",   FALSE, etBOOL, {&addHydrogens},
 218           "Add hydrogen atoms to the compound - useful for PDB files." },
 219         { "-qcycle", FALSE, etINT, {&qcycle},
 220           "Max number of tries for optimizing the charges. The trial with lowest chi2 will be used for generating a topology. Will be turned off if randzeta is No." },
 221         { "-hfac",    FALSE, etREAL, {&hfac},
 222           "HIDDENFudge factor for AXx algorithms that modulates J00 for hydrogen atoms by multiplying it by (1 + hfac*qH). This hack is originally due to Rappe & Goddard." },
 223         { "-qsymm",  FALSE, etBOOL, {&bQsym},
 224           "Symmetrize the charges on symmetric groups, e.g. CH3, NH2." },
 225         { "-symm",   FALSE, etSTR, {&symm_string},
 226           "Use the order given here for symmetrizing, e.g. when specifying [TT]-symm '0 1 0'[tt] for a water molecule (H-O-H) the hydrogens will have obtain the same charge. For simple groups, like methyl (or water) this is done automatically, but higher symmetry is not detected by the program. The numbers should correspond to atom numbers minus 1, and point to either the atom itself or to a previous atom." },
 227         { "-cgsort", FALSE, etSTR, {cgopt},
 228           "HIDDENOption for assembling charge groups: based on Atom (default, does not change the atom order), Group (e.g. CH3 groups are kept together), or Neutral sections (try to find groups that together are neutral). If the order of atoms is changed an index file is written in order to facilitate changing the order in old files." },
 229         { "-nexcl",    FALSE, etINT, {&nexcl},
 230           "HIDDENNumber of exclusion" },
 231         { "-jobtype",  FALSE, etSTR, {&jobtype},
 232           "The job type used in the Gaussian calculation: Opt, Polar, SP, and etc." },
 233         { "-efield",  FALSE, etREAL, {&efield},
 234           "The magnitude of the external electeric field to calculate polarizability tensor." },
 235     };
 236
 237     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 238                            asize(desc), desc, asize(bugs), bugs, &oenv))
 239     {
 240         return 0;
 241     }
 242
 243     alexandria::Poldata       pd;
 244     t_inputrec               *inputrec = new t_inputrec();
 245     t_commrec                *cr       = init_commrec();
 246     const char               *tabfn    = opt2fn_null("-table", NFILE, fnm);
 247     eChargeGroup              ecg      = (eChargeGroup) get_option(cgopt);
 248     gmx::MDLogger             mdlog {};
 249     ChargeModel               iModel;
 250     std::string               method, basis;
 251     splitLot(lot, &method, &basis);
 252
 253     /* Check the options */
 254     bITP = opt2bSet("-oi", NFILE, fnm);
 255     if ((qtol < 0) || (qtol > 1))
 256     {
 257         gmx_fatal(FARGS, "Charge tolerance should be between 0 and 1 (not %g)", qtol);
 258     }
 259     // Check whether there is something to read
 260     if (strlen(dbname) == 0 && strlen(filename) == 0)
 261     {
 262         gmx_fatal(FARGS, "Specify either the -db or the -f option. No output without input");
 263     }
 264     const char *gentop_fnm = opt2fn_null("-d", NFILE, fnm);
 265     if (opt2parg_bSet("-ff", asize(pa), pa) && nullptr == gentop_fnm)
 266     {
 267         iModel     = name2eemtype(ff[0]);
 268         gentop_fnm = gmx::formatString("%s_2020.dat", ff[0]).c_str();
 269     }
 270     if (nullptr == gentop_fnm)
 271     {
 272         fprintf(stderr, "Please specify either a force field file name or use the -ff flag");
 273         return 0;
 274     }
 275
 276     /* Read standard atom properties */
 277     aps = gmx_atomprop_init();
 278     try
 279     {
 280         alexandria::readPoldata(gentop_fnm, pd, aps);
 281     }
 282     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 283     iModel = pd.getChargeModel();
 284     printf("Using force field file %s and charge distribution model %s\n",
 285            gentop_fnm, getEemtypeName(iModel));
 286     if (bVerbose)
 287     {
 288         printf("Reading force field information. There are %d atomtypes.\n",
 289                static_cast<int>(pd.getNatypes()));
 290     }
 291
 292     if (pd.getNexcl() != nexcl && opt2parg_bSet("-nexcl", asize(pa), pa))
 293     {
 294         fprintf(stderr, "WARNING: Changing exclusion number from %d in force field file\n", pd.getNexcl());
 295         fprintf(stderr, "         to %d (command line), Please check your output carefully.\n", nexcl);
 296         pd.setNexcl(nexcl);
 297     }
 298     if (strlen(dbname) > 0)
 299     {
 300         const char *molpropDatabase = opt2fn_null("-mpdb", NFILE, fnm);
 301         if (!molpropDatabase || strlen(molpropDatabase) == 0)
 302         {
 303             gmx_fatal(FARGS, "Empty database file name");
 304         }
 305         if (bVerbose)
 306         {
 307             printf("Looking up %s in molecule database %s.\n",
 308                    dbname, molpropDatabase);
 309         }
 310         MolPropRead(molpropDatabase, &mps);
 311         for (mpi = mps.begin(); (mpi < mps.end()); mpi++)
 312         {
 313             if (strcasecmp(dbname, mpi->getMolname().c_str()) == 0)
 314             {
 315                 break;
 316             }
 317         }
 318         if (mpi == mps.end())
 319         {
 320             gmx_fatal(FARGS, "Molecule %s not found in database", dbname);
 321         }
 322         mymol.Merge(mpi);
 323     }
 324     else
 325     {
 326         std::string fnm = filename;
 327         auto fns = gmx::splitString(fnm);
 328         if (strlen(molnm) == 0)
 329         {
 330             molnm = (char *)"MOL";
 331         }
 332         if (fns.size() > 0)
 333         {
 334             for (auto &i : fns)
 335             {
 336                 alexandria::MolProp  mp;
 337                 if (readBabel(i.c_str(),
 338                               &mp,
 339                               molnm,
 340                               iupac,
 341                               conf,
 342                               basis.c_str(),
 343                               maxpot,
 344                               nsymm,
 345                               jobtype,
 346                               qtot,
 347                               addHydrogens))
 348                 {
 349                     mps.push_back(mp);
 350                 }
 351             }
 352         }
 353         else
 354         {
 355             gmx_fatal(FARGS, "No input file has been specified.");
 356         }
 357         for (auto mpi = mps.begin(); mpi < mps.end(); mpi++)
 358         {
 359             mymol.Merge(mpi);
 360         }
 361     }
 362     mymol.SetForceField(ff[0]);
 363     fill_inputrec(inputrec);
 364     mymol.setInputrec(inputrec);
 365     std::string mylot;
 366     imm = mymol.GenerateTopology(aps,
 367                                  &pd,
 368                                  method,
 369                                  basis,
 370                                  &mylot,
 371                                  bGenVSites,
 372                                  bPairs,
 373                                  bDihedral,
 374                                  false,
 375                                  tabfn);
 376
 377     if (immOK == imm)
 378     {
 379         maxpot = 100; //Use 100 percent of the ESP read from Gaussian file.
 380         imm    = mymol.GenerateCharges(&pd,
 381                                        mdlog,
 382                                        aps,
 383                                        watoms,
 384                                        hfac,
 385                                        method,
 386                                        basis,
 387                                        &mylot,
 388                                        bQsym,
 389                                        symm_string,
 390                                        cr,
 391                                        tabfn,
 392                                        nullptr,
 393                                        qcycle,
 394                                        maxpot,
 395                                        qtol,
 396                                        oenv,
 397                                        opt2fn_null("-plotESP", NFILE, fnm));
 398     }
 399     if (bCUBE && immOK == imm)
 400     {
 401         fprintf(stderr, "Fix me: GenerateCube is broken\n");
 402         mymol.GenerateCube(&pd,
 403                            spacing,
 404                            opt2fn_null("-ref",      NFILE, fnm),
 405                            opt2fn_null("-pc",       NFILE, fnm),
 406                            opt2fn_null("-pdbdiff",  NFILE, fnm),
 407                            opt2fn_null("-pot",      NFILE, fnm),
 408                            opt2fn_null("-rho",      NFILE, fnm),
 409                            opt2fn_null("-his",      NFILE, fnm),
 410                            opt2fn_null("-diff",     NFILE, fnm),
 411                            opt2fn_null("-diffhist", NFILE, fnm),
 412                            oenv);
 413     }
 414
 415     if (immOK == imm)
 416     {
 417         imm = mymol.GenerateChargeGroups(ecg, bUsePDBcharge);
 418     }
 419
 420     if (immOK == imm && mymol.errors().size() == 0)
 421     {
 422         splitLot(mylot.c_str(), &method, &basis);
 423         mymol.PrintConformation(opt2fn("-c", NFILE, fnm));
 424         mymol.PrintTopology(bITP ? ftp2fn(efITP, NFILE, fnm) : ftp2fn(efTOP, NFILE, fnm),
 425                             bVerbose,
 426                             &pd,
 427                             aps,
 428                             cr,
 429                             efield,
 430                             method,
 431                             basis);
 432     }
 433     else
 434     {
 435         auto fn = opt2fn("-g", NFILE, fnm);
 436         fprintf(stderr, "\nFatal Error. Please check the %s file for error messages.\n", fn);
 437         print_errors(fn, mymol.errors(), imm);
 438     }
 439     return 0;
 440 }