src/gromacs/tools/convert_tpr.cpp

   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-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include "convert_tpr.h"
  41
  42 #include <cmath>
  43
  44 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  45 #include "gromacs/fileio/checkpoint.h"
  46 #include "gromacs/fileio/enxio.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trrio.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/mdtypes/inputrec.h"
  51 #include "gromacs/mdtypes/md_enums.h"
  52 #include "gromacs/mdtypes/state.h"
  53 #include "gromacs/options/basicoptions.h"
  54 #include "gromacs/options/filenameoption.h"
  55 #include "gromacs/options/ioptionscontainer.h"
  56 #include "gromacs/random/seed.h"
  57 #include "gromacs/topology/ifunc.h"
  58 #include "gromacs/topology/index.h"
  59 #include "gromacs/topology/mtop_util.h"
  60 #include "gromacs/topology/topology.h"
  61 #include "gromacs/utility/arraysize.h"
  62 #include "gromacs/utility/cstringutil.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/real.h"
  66 #include "gromacs/utility/smalloc.h"
  67 #include "gromacs/utility/stringutil.h"
  68
  69 static void rangeCheck(int numberInIndexFile, int maxAtomNumber)
  70 {
  71     if ((numberInIndexFile) >= (maxAtomNumber))
  72     {
  73         gmx_fatal(FARGS,
  74                   "Your index file contains atomnumbers (e.g. %d)\nthat are larger than the number "
  75                   "of atoms in the tpr file (%d)",
  76                   (numberInIndexFile),
  77                   (maxAtomNumber));
  78     }
  79 }
  80
  81 static std::vector<bool> bKeepIt(int gnx, int natoms, int index[])
  82 {
  83     std::vector<bool> b(natoms);
  84
  85     for (int i = 0; (i < gnx); i++)
  86     {
  87         rangeCheck(index[i], natoms);
  88         b[index[i]] = TRUE;
  89     }
  90
  91     return b;
  92 }
  93
  94 static std::vector<int> invind(int gnx, int natoms, int index[])
  95 {
  96     std::vector<int> inv(natoms);
  97
  98     for (int i = 0; (i < gnx); i++)
  99     {
 100         rangeCheck(index[i], natoms);
 101         inv[index[i]] = i;
 102     }
 103
 104     return inv;
 105 }
 106
 107 static gmx::ListOfLists<int> reduce_listoflists(gmx::ArrayRef<const int>     invindex,
 108                                                 const std::vector<bool>&     bKeep,
 109                                                 const gmx::ListOfLists<int>& src,
 110                                                 const char*                  name)
 111 {
 112     gmx::ListOfLists<int> lists;
 113
 114     std::vector<int> exclusionsForAtom;
 115     for (gmx::index i = 0; i < src.ssize(); i++)
 116     {
 117         if (bKeep[i])
 118         {
 119             exclusionsForAtom.clear();
 120             for (const int j : src[i])
 121             {
 122                 if (bKeep[j])
 123                 {
 124                     exclusionsForAtom.push_back(invindex[j]);
 125                 }
 126             }
 127             lists.pushBack(exclusionsForAtom);
 128         }
 129     }
 130
 131     fprintf(stderr,
 132             "Reduced block %8s from %6zu to %6zu index-, %6d to %6d a-entries\n",
 133             name,
 134             src.size(),
 135             lists.size(),
 136             src.numElements(),
 137             lists.numElements());
 138
 139     return lists;
 140 }
 141
 142 static void reduce_rvec(int gnx, const int index[], rvec vv[])
 143 {
 144     rvec* ptr;
 145     int   i;
 146
 147     snew(ptr, gnx);
 148     for (i = 0; (i < gnx); i++)
 149     {
 150         copy_rvec(vv[index[i]], ptr[i]);
 151     }
 152     for (i = 0; (i < gnx); i++)
 153     {
 154         copy_rvec(ptr[i], vv[i]);
 155     }
 156     sfree(ptr);
 157 }
 158
 159 static void reduce_atom(int gnx, const int index[], t_atom atom[], char*** atomname, int* nres, t_resinfo* resinfo)
 160 {
 161     t_atom*    ptr;
 162     char***    aname;
 163     t_resinfo* rinfo;
 164     int        i, nr;
 165
 166     snew(ptr, gnx);
 167     snew(aname, gnx);
 168     snew(rinfo, atom[index[gnx - 1]].resind + 1);
 169     for (i = 0; (i < gnx); i++)
 170     {
 171         ptr[i]   = atom[index[i]];
 172         aname[i] = atomname[index[i]];
 173     }
 174     nr = -1;
 175     for (i = 0; (i < gnx); i++)
 176     {
 177         atom[i]     = ptr[i];
 178         atomname[i] = aname[i];
 179         if ((i == 0) || (atom[i].resind != atom[i - 1].resind))
 180         {
 181             nr++;
 182             rinfo[nr] = resinfo[atom[i].resind];
 183         }
 184         atom[i].resind = nr;
 185     }
 186     nr++;
 187     for (i = 0; (i < nr); i++)
 188     {
 189         resinfo[i] = rinfo[i];
 190     }
 191     *nres = nr;
 192
 193     sfree(aname);
 194     sfree(ptr);
 195     sfree(rinfo);
 196 }
 197
 198 static void reduce_ilist(gmx::ArrayRef<const int> invindex,
 199                          const std::vector<bool>& bKeep,
 200                          InteractionList*         il,
 201                          int                      nratoms,
 202                          const char*              name)
 203 {
 204     if (!il->empty())
 205     {
 206         std::vector<int> newAtoms(nratoms);
 207         InteractionList  ilReduced;
 208         for (int i = 0; i < il->size(); i += nratoms + 1)
 209         {
 210             bool bB = true;
 211             for (int j = 0; j < nratoms; j++)
 212             {
 213                 bB = bB && bKeep[il->iatoms[i + 1 + j]];
 214             }
 215             if (bB)
 216             {
 217                 for (int j = 0; j < nratoms; j++)
 218                 {
 219                     newAtoms[j] = invindex[il->iatoms[i + 1 + j]];
 220                 }
 221                 ilReduced.push_back(il->iatoms[i], nratoms, newAtoms.data());
 222             }
 223         }
 224         fprintf(stderr,
 225                 "Reduced ilist %8s from %6d to %6d entries\n",
 226                 name,
 227                 il->size() / (nratoms + 1),
 228                 ilReduced.size() / (nratoms + 1));
 229
 230         *il = std::move(ilReduced);
 231     }
 232 }
 233
 234 static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[], rvec v[])
 235 {
 236     gmx_localtop_t top(mtop->ffparams);
 237     gmx_mtop_generate_local_top(*mtop, &top, false);
 238     t_atoms atoms = gmx_mtop_global_atoms(mtop);
 239
 240     const std::vector<bool> bKeep    = bKeepIt(gnx, atoms.nr, index);
 241     const std::vector<int>  invindex = invind(gnx, atoms.nr, index);
 242
 243     reduce_rvec(gnx, index, x);
 244     reduce_rvec(gnx, index, v);
 245     reduce_atom(gnx, index, atoms.atom, atoms.atomname, &(atoms.nres), atoms.resinfo);
 246
 247     for (int i = 0; (i < F_NRE); i++)
 248     {
 249         reduce_ilist(invindex,
 250                      bKeep,
 251                      &(top.idef.il[i]),
 252                      interaction_function[i].nratoms,
 253                      interaction_function[i].name);
 254     }
 255
 256     atoms.nr = gnx;
 257
 258     mtop->moltype.resize(1);
 259     mtop->moltype[0].name  = mtop->name;
 260     mtop->moltype[0].atoms = atoms;
 261     mtop->moltype[0].excls = reduce_listoflists(invindex, bKeep, top.excls, "excls");
 262     for (int i = 0; i < F_NRE; i++)
 263     {
 264         mtop->moltype[0].ilist[i] = std::move(top.idef.il[i]);
 265     }
 266
 267     mtop->molblock.resize(1);
 268     mtop->molblock[0].type = 0;
 269     mtop->molblock[0].nmol = 1;
 270
 271     mtop->natoms = atoms.nr;
 272 }
 273
 274 static void zeroq(const int index[], gmx_mtop_t* mtop)
 275 {
 276     for (gmx_moltype_t& moltype : mtop->moltype)
 277     {
 278         for (int i = 0; i < moltype.atoms.nr; i++)
 279         {
 280             moltype.atoms.atom[index[i]].q  = 0;
 281             moltype.atoms.atom[index[i]].qB = 0;
 282         }
 283     }
 284 }
 285
 286 namespace gmx
 287 {
 288
 289 namespace
 290 {
 291
 292 class ConvertTpr : public ICommandLineOptionsModule
 293 {
 294 public:
 295     ConvertTpr() {}
 296
 297     // From ICommandLineOptionsModule
 298     void init(CommandLineModuleSettings* /*settings*/) override {}
 299     void initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings) override;
 300     void optionsFinished() override;
 301     int  run() override;
 302
 303 private:
 304     //! Name of input tpr file.
 305     std::string inputTprFileName_;
 306     //! Name of input index file.
 307     std::string inputIndexFileName_;
 308     //! Name of output tpr file.
 309     std::string outputTprFileName_;
 310     //! If we have read in an index file.
 311     bool haveReadIndexFile_ = false;
 312     //! Time to extend simulation by.
 313     real extendTime_ = 0;
 314     //! If the option to extend simulation time is set.
 315     bool extendTimeIsSet_ = false;
 316     //! Final run time value.
 317     real runToMaxTime_ = 0;
 318     //! If the option to run simulation until specified time is set.
 319     bool runToMaxTimeIsSet_ = false;
 320     //! Maximum number of steps to run.
 321     int64_t maxSteps_ = 0;
 322     //! If the option to use maximumstep number is set.
 323     bool maxStepsIsSet_ = false;
 324     //! If the option to zero charge is set.
 325     bool zeroQIsSet_ = false;
 326 };
 327
 328 void ConvertTpr::initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings)
 329 {
 330     std::vector<const char*> desc = {
 331         "[THISMODULE] can edit run input files in three ways.[PAR]",
 332         "[BB]1.[bb] by modifying the number of steps in a run input file",
 333         "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
 334         "(nsteps=-1 means unlimited number of steps)[PAR]",
 335         "[BB]2.[bb] by creating a [REF].tpx[ref] file for a subset of your original",
 336         "tpx file, which is useful when you want to remove the solvent from",
 337         "your [REF].tpx[ref] file, or when you want to make e.g. a pure C[GRK]alpha[grk] ",
 338         "[REF].tpx[ref] file.",
 339         "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
 340         "this to work.",
 341         "[BB]WARNING: this [REF].tpx[ref] file is not fully functional[bb].[PAR]",
 342         "[BB]3.[bb] by setting the charges of a specified group",
 343         "to zero. This is useful when doing free energy estimates",
 344         "using the LIE (Linear Interaction Energy) method."
 345     };
 346
 347     settings->setHelpText(desc);
 348
 349     options->addOption(FileNameOption("s")
 350                                .filetype(eftTopology)
 351                                .inputFile()
 352                                .required()
 353                                .store(&inputTprFileName_)
 354                                .defaultBasename("topol")
 355                                .description("Run input file to modify"));
 356     options->addOption(FileNameOption("n")
 357                                .filetype(eftIndex)
 358                                .inputFile()
 359                                .store(&inputIndexFileName_)
 360                                .storeIsSet(&haveReadIndexFile_)
 361                                .defaultBasename("index")
 362                                .description("File containing additional index groups"));
 363     options->addOption(FileNameOption("o")
 364                                .filetype(eftTopology)
 365                                .outputFile()
 366                                .store(&outputTprFileName_)
 367                                .defaultBasename("tprout")
 368                                .description("Generated modified run input file"));
 369     options->addOption(RealOption("extend")
 370                                .store(&extendTime_)
 371                                .storeIsSet(&extendTimeIsSet_)
 372                                .timeValue()
 373                                .description("Extend runtime by this amount (ps)"));
 374     options->addOption(RealOption("until")
 375                                .store(&runToMaxTime_)
 376                                .storeIsSet(&runToMaxTimeIsSet_)
 377                                .timeValue()
 378                                .description("Extend runtime until this ending time (ps)"));
 379     options->addOption(
 380             Int64Option("nsteps").store(&maxSteps_).storeIsSet(&maxStepsIsSet_).description("Change the number of steps"));
 381     options->addOption(
 382             BooleanOption("zeroq").store(&zeroQIsSet_).description("Set the charges of a group (from the index) to zero"));
 383 }
 384
 385 void ConvertTpr::optionsFinished() {}
 386
 387 int ConvertTpr::run()
 388 {
 389     gmx_mtop_t mtop;
 390     t_atoms    atoms;
 391     t_state    state;
 392     char       buf[200], buf2[200];
 393
 394     fprintf(stderr, "Reading toplogy and stuff from %s\n", inputTprFileName_.c_str());
 395
 396     t_inputrec  irInstance;
 397     t_inputrec* ir = &irInstance;
 398     read_tpx_state(inputTprFileName_.c_str(), ir, &state, &mtop);
 399     int64_t currentMaxStep    = ir->init_step;
 400     double  currentRunTime    = ir->init_step * ir->delta_t + ir->init_t;
 401     real    currentMaxRunTime = 0.0;
 402
 403     if (maxStepsIsSet_)
 404     {
 405         fprintf(stderr, "Setting nsteps to %s\n", gmx_step_str(maxSteps_, buf));
 406         ir->nsteps = maxSteps_;
 407     }
 408     else
 409     {
 410         /* Determine total number of steps remaining */
 411         if (extendTimeIsSet_)
 412         {
 413             ir->nsteps = ir->nsteps - (currentMaxStep - ir->init_step)
 414                          + gmx::roundToInt64(extendTime_ / ir->delta_t);
 415             printf("Extending remaining runtime of by %g ps (now %s steps)\n",
 416                    extendTime_,
 417                    gmx_step_str(ir->nsteps, buf));
 418         }
 419         else if (runToMaxTimeIsSet_)
 420         {
 421             printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
 422                    gmx_step_str(ir->nsteps, buf),
 423                    gmx_step_str(currentMaxStep, buf2),
 424                    currentRunTime,
 425                    runToMaxTime_);
 426             ir->nsteps = gmx::roundToInt64((currentMaxRunTime - currentRunTime) / ir->delta_t);
 427             printf("Extending remaining runtime until %g ps (now %s steps)\n",
 428                    currentMaxRunTime,
 429                    gmx_step_str(ir->nsteps, buf));
 430         }
 431         else
 432         {
 433             ir->nsteps -= currentMaxStep - ir->init_step;
 434             /* Print message */
 435             printf("%s steps (%g ps) remaining from first run.\n",
 436                    gmx_step_str(ir->nsteps, buf),
 437                    ir->nsteps * ir->delta_t);
 438         }
 439     }
 440
 441     if (maxStepsIsSet_ || zeroQIsSet_ || (ir->nsteps > 0))
 442     {
 443         ir->init_step = currentMaxStep;
 444
 445         if (haveReadIndexFile_ || !(maxStepsIsSet_ || extendTimeIsSet_ || runToMaxTimeIsSet_))
 446         {
 447             atoms         = gmx_mtop_global_atoms(&mtop);
 448             int   gnx     = 0;
 449             int*  index   = nullptr;
 450             char* grpname = nullptr;
 451             get_index(&atoms, inputIndexFileName_.c_str(), 1, &gnx, &index, &grpname);
 452             bool bSel = false;
 453             if (!zeroQIsSet_)
 454             {
 455                 bSel = (gnx != state.natoms);
 456                 for (int i = 0; ((i < gnx) && (!bSel)); i++)
 457                 {
 458                     bSel = (i != index[i]);
 459                 }
 460             }
 461             else
 462             {
 463                 bSel = false;
 464             }
 465             if (bSel)
 466             {
 467                 fprintf(stderr,
 468                         "Will write subset %s of original tpx containing %d "
 469                         "atoms\n",
 470                         grpname,
 471                         gnx);
 472                 reduce_topology_x(gnx, index, &mtop, state.x.rvec_array(), state.v.rvec_array());
 473                 state.natoms = gnx;
 474             }
 475             else if (zeroQIsSet_)
 476             {
 477                 zeroq(index, &mtop);
 478                 fprintf(stderr, "Zero-ing charges for group %s\n", grpname);
 479             }
 480             else
 481             {
 482                 fprintf(stderr, "Will write full tpx file (no selection)\n");
 483             }
 484         }
 485
 486         double stateTime = ir->init_t + ir->init_step * ir->delta_t;
 487         sprintf(buf,
 488                 "Writing statusfile with starting step %s%s and length %s%s steps...\n",
 489                 "%10",
 490                 PRId64,
 491                 "%10",
 492                 PRId64);
 493         fprintf(stderr, buf, ir->init_step, ir->nsteps);
 494         fprintf(stderr,
 495                 "                                 time %10.3f and length %10.3f ps\n",
 496                 stateTime,
 497                 ir->nsteps * ir->delta_t);
 498         write_tpx_state(outputTprFileName_.c_str(), ir, &state, &mtop);
 499     }
 500     else
 501     {
 502         printf("You've simulated long enough. Not writing tpr file\n");
 503     }
 504
 505     return 0;
 506 }
 507
 508 } // namespace
 509
 510 const char ConvertTprInfo::name[]             = "convert-tpr";
 511 const char ConvertTprInfo::shortDescription[] = "Make a modifed run-input file";
 512 ICommandLineOptionsModulePointer ConvertTprInfo::create()
 513 {
 514     return ICommandLineOptionsModulePointer(std::make_unique<ConvertTpr>());
 515 }
 516
 517 } // namespace gmx