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