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,2021, 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 static void print_runtime_info(t_inputrec* ir)
 287 {
 288     char buf[STEPSTRSIZE];
 289
 290     printf("  Run start step                %22s     \n", gmx_step_str(ir->init_step, buf));
 291     printf("  Run start time                %22g ps  \n", ir->init_step * ir->delta_t + ir->init_t);
 292     printf("  Step to be made during run    %22s     \n", gmx_step_str(ir->nsteps, buf));
 293     printf("  Runtime for the run           %22g ps  \n", ir->nsteps * ir->delta_t);
 294     printf("  Run end step                  %22s     \n", gmx_step_str(ir->init_step + ir->nsteps, buf));
 295     printf("  Run end time                  %22g ps  \n\n",
 296            (ir->init_step + ir->nsteps) * ir->delta_t + ir->init_t);
 297 }
 298
 299 namespace gmx
 300 {
 301
 302 namespace
 303 {
 304
 305 class ConvertTpr : public ICommandLineOptionsModule
 306 {
 307 public:
 308     ConvertTpr() {}
 309
 310     // From ICommandLineOptionsModule
 311     void init(CommandLineModuleSettings* /*settings*/) override {}
 312     void initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings) override;
 313     void optionsFinished() override;
 314     int  run() override;
 315
 316 private:
 317     //! Name of input tpr file.
 318     std::string inputTprFileName_;
 319     //! Name of input index file.
 320     std::string inputIndexFileName_;
 321     //! Name of output tpr file.
 322     std::string outputTprFileName_;
 323     //! If we have read in an index file.
 324     bool haveReadIndexFile_ = false;
 325     //! Time to extend simulation by.
 326     real extendTime_ = 0;
 327     //! If the option to extend simulation time is set.
 328     bool extendTimeIsSet_ = false;
 329     //! Final run time value.
 330     real runToMaxTime_ = 0;
 331     //! If the option to run simulation until specified time is set.
 332     bool runToMaxTimeIsSet_ = false;
 333     //! Maximum number of steps to run.
 334     int64_t maxSteps_ = 0;
 335     //! If the option to use maximumstep number is set.
 336     bool maxStepsIsSet_ = false;
 337     //! If the option to zero charge is set.
 338     bool zeroQIsSet_ = false;
 339 };
 340
 341 void ConvertTpr::initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings)
 342 {
 343     std::vector<const char*> desc = {
 344         "[THISMODULE] can edit run input files in three ways.[PAR]",
 345         "[BB]1.[bb] by modifying the number of steps in a run input file",
 346         "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
 347         "(nsteps=-1 means unlimited number of steps)[PAR]",
 348         "[BB]2.[bb] by creating a [REF].tpx[ref] file for a subset of your original",
 349         "tpx file, which is useful when you want to remove the solvent from",
 350         "your [REF].tpx[ref] file, or when you want to make e.g. a pure C[GRK]alpha[grk] ",
 351         "[REF].tpx[ref] file.",
 352         "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
 353         "this to work.",
 354         "[BB]WARNING: this [REF].tpx[ref] file is not fully functional[bb].[PAR]",
 355         "[BB]3.[bb] by setting the charges of a specified group",
 356         "to zero. This is useful when doing free energy estimates",
 357         "using the LIE (Linear Interaction Energy) method."
 358     };
 359
 360     settings->setHelpText(desc);
 361
 362     options->addOption(FileNameOption("s")
 363                                .filetype(OptionFileType::Topology)
 364                                .inputFile()
 365                                .required()
 366                                .store(&inputTprFileName_)
 367                                .defaultBasename("topol")
 368                                .description("Run input file to modify"));
 369     options->addOption(FileNameOption("n")
 370                                .filetype(OptionFileType::Index)
 371                                .inputFile()
 372                                .store(&inputIndexFileName_)
 373                                .storeIsSet(&haveReadIndexFile_)
 374                                .defaultBasename("index")
 375                                .description("File containing additional index groups"));
 376     options->addOption(FileNameOption("o")
 377                                .filetype(OptionFileType::Topology)
 378                                .outputFile()
 379                                .store(&outputTprFileName_)
 380                                .defaultBasename("tprout")
 381                                .description("Generated modified run input file"));
 382     options->addOption(RealOption("extend")
 383                                .store(&extendTime_)
 384                                .storeIsSet(&extendTimeIsSet_)
 385                                .timeValue()
 386                                .description("Extend runtime by this amount (ps)"));
 387     options->addOption(RealOption("until")
 388                                .store(&runToMaxTime_)
 389                                .storeIsSet(&runToMaxTimeIsSet_)
 390                                .timeValue()
 391                                .description("Extend runtime until this ending time (ps)"));
 392     options->addOption(Int64Option("nsteps")
 393                                .store(&maxSteps_)
 394                                .storeIsSet(&maxStepsIsSet_)
 395                                .description("Change the number of steps remaining to be made"));
 396     options->addOption(
 397             BooleanOption("zeroq").store(&zeroQIsSet_).description("Set the charges of a group (from the index) to zero"));
 398 }
 399
 400 void ConvertTpr::optionsFinished() {}
 401
 402 int ConvertTpr::run()
 403 {
 404     gmx_mtop_t mtop;
 405     t_atoms    atoms;
 406     t_state    state;
 407     char       buf[STEPSTRSIZE];
 408
 409     if (static_cast<int>(extendTimeIsSet_) + static_cast<int>(maxStepsIsSet_) + static_cast<int>(runToMaxTimeIsSet_)
 410         > 1)
 411     {
 412         printf("Multiple runtime modification operations cannot be done in a single call.\n");
 413         return 1;
 414     }
 415
 416     if ((extendTimeIsSet_ || maxStepsIsSet_ || runToMaxTimeIsSet_) && (zeroQIsSet_ || haveReadIndexFile_))
 417     {
 418         printf("Cannot do both runtime modification and charge-zeroing/index group extraction in a "
 419                "single call.\n");
 420         return 1;
 421     }
 422
 423     if (zeroQIsSet_ && !haveReadIndexFile_)
 424     {
 425         printf("Charge zeroing need an index file.\n");
 426         return 1;
 427     }
 428
 429
 430     t_inputrec  irInstance;
 431     t_inputrec* ir = &irInstance;
 432     read_tpx_state(inputTprFileName_.c_str(), ir, &state, &mtop);
 433
 434
 435     if (extendTimeIsSet_ || maxStepsIsSet_ || runToMaxTimeIsSet_)
 436     {
 437         // Doing runtime modification
 438
 439         const double  inputTimeAtStartOfRun = ir->init_step * ir->delta_t + ir->init_t;
 440         const int64_t inputStepAtEndOfRun   = ir->init_step + ir->nsteps;
 441         const double  inputTimeAtEndOfRun   = inputStepAtEndOfRun * ir->delta_t + ir->init_t;
 442
 443         printf("Input file:\n");
 444         print_runtime_info(ir);
 445
 446         if (maxStepsIsSet_)
 447         {
 448             fprintf(stderr, "Setting nsteps to %s\n", gmx_step_str(maxSteps_, buf));
 449             ir->nsteps = maxSteps_;
 450         }
 451         else if (extendTimeIsSet_)
 452         {
 453             printf("Extending remaining runtime by %g ps\n", extendTime_);
 454             ir->nsteps += gmx::roundToInt64(extendTime_ / ir->delta_t);
 455         }
 456         else if (runToMaxTimeIsSet_)
 457         {
 458             if (runToMaxTime_ <= inputTimeAtStartOfRun)
 459             {
 460                 printf("The requested run end time is at/before the run start time.\n");
 461                 return 1;
 462             }
 463             if (runToMaxTime_ < inputTimeAtEndOfRun)
 464             {
 465                 printf("The requested run end time is before the original run end time.\n");
 466                 printf("Reducing remaining runtime to %g ps\n", runToMaxTime_);
 467             }
 468             else
 469             {
 470                 printf("Extending remaining runtime to %g ps\n", runToMaxTime_);
 471             }
 472             ir->nsteps = gmx::roundToInt64((runToMaxTime_ - inputTimeAtStartOfRun) / ir->delta_t);
 473         }
 474
 475         printf("\nOutput file:\n");
 476         print_runtime_info(ir);
 477     }
 478     else
 479     {
 480         // If zeroQIsSet_, then we are doing charge zero-ing; otherwise index group extraction
 481         // In both cases an index filename has been provided
 482
 483         atoms         = gmx_mtop_global_atoms(mtop);
 484         int   gnx     = 0;
 485         int*  index   = nullptr;
 486         char* grpname = nullptr;
 487         get_index(&atoms, inputIndexFileName_.c_str(), 1, &gnx, &index, &grpname);
 488         bool bSel = false;
 489         if (!zeroQIsSet_)
 490         {
 491             bSel = (gnx != state.natoms);
 492             for (int i = 0; ((i < gnx) && (!bSel)); i++)
 493             {
 494                 bSel = (i != index[i]);
 495             }
 496         }
 497         else
 498         {
 499             bSel = false;
 500         }
 501         if (bSel)
 502         {
 503             fprintf(stderr,
 504                     "Will write subset %s of original tpx containing %d "
 505                     "atoms\n",
 506                     grpname,
 507                     gnx);
 508             reduce_topology_x(gnx, index, &mtop, state.x.rvec_array(), state.v.rvec_array());
 509             state.natoms = gnx;
 510         }
 511         else if (zeroQIsSet_)
 512         {
 513             zeroq(index, &mtop);
 514             fprintf(stderr, "Zero-ing charges for group %s\n", grpname);
 515         }
 516         else
 517         {
 518             fprintf(stderr, "Will write full tpx file (no selection)\n");
 519         }
 520     }
 521
 522     // Writing output tpx regardless of the operation performed
 523     write_tpx_state(outputTprFileName_.c_str(), ir, &state, mtop);
 524
 525     return 0;
 526 }
 527
 528 } // namespace
 529
 530 const char ConvertTprInfo::name[]             = "convert-tpr";
 531 const char ConvertTprInfo::shortDescription[] = "Make a modifed run-input file";
 532 ICommandLineOptionsModulePointer ConvertTprInfo::create()
 533 {
 534     return ICommandLineOptionsModulePointer(std::make_unique<ConvertTpr>());
 535 }
 536
 537 } // namespace gmx