src/gromacs/gmxpreprocess/insert_molecules.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 "insert_molecules.h"
  41
  42 #include <algorithm>
  43 #include <memory>
  44 #include <set>
  45 #include <string>
  46 #include <vector>
  47
  48 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/fileio/filetypes.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gromacs/gmxlib/conformation_utilities.h"
  53 #include "gromacs/gmxpreprocess/makeexclusiondistances.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/math/utilities.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/options/basicoptions.h"
  58 #include "gromacs/options/filenameoption.h"
  59 #include "gromacs/options/ioptionscontainer.h"
  60 #include "gromacs/pbcutil/pbc.h"
  61 #include "gromacs/random/threefry.h"
  62 #include "gromacs/random/uniformrealdistribution.h"
  63 #include "gromacs/selection/nbsearch.h"
  64 #include "gromacs/selection/selection.h"
  65 #include "gromacs/selection/selectioncollection.h"
  66 #include "gromacs/selection/selectionoption.h"
  67 #include "gromacs/selection/selectionoptionbehavior.h"
  68 #include "gromacs/topology/atomprop.h"
  69 #include "gromacs/topology/atoms.h"
  70 #include "gromacs/topology/atomsbuilder.h"
  71 #include "gromacs/topology/mtop_util.h"
  72 #include "gromacs/topology/symtab.h"
  73 #include "gromacs/topology/topology.h"
  74 #include "gromacs/trajectory/trajectoryframe.h"
  75 #include "gromacs/utility/cstringutil.h"
  76 #include "gromacs/utility/exceptions.h"
  77 #include "gromacs/utility/fatalerror.h"
  78 #include "gromacs/utility/smalloc.h"
  79
  80 using gmx::RVec;
  81
  82 /* enum for random rotations of inserted solutes */
  83 enum RotationType
  84 {
  85     en_rotXYZ,
  86     en_rotZ,
  87     en_rotNone
  88 };
  89 const char* const cRotationEnum[] = { "xyz", "z", "none" };
  90
  91 static void center_molecule(gmx::ArrayRef<RVec> x)
  92 {
  93     rvec center = { 0 };
  94     for (auto& xi : x)
  95     {
  96         rvec_inc(center, xi);
  97     }
  98     svmul(1.0 / x.size(), center, center);
  99     for (auto& xi : x)
 100     {
 101         rvec_dec(xi, center);
 102     }
 103 }
 104
 105 static void generate_trial_conf(gmx::ArrayRef<RVec>       xin,
 106                                 const rvec                offset,
 107                                 RotationType              enum_rot,
 108                                 gmx::DefaultRandomEngine* rng,
 109                                 std::vector<RVec>*        xout)
 110 {
 111     gmx::UniformRealDistribution<real> dist(0, 2.0 * M_PI);
 112     xout->assign(xin.begin(), xin.end());
 113
 114     real alfa = 0.0, beta = 0.0, gamma = 0.0;
 115     switch (enum_rot)
 116     {
 117         case en_rotXYZ:
 118             alfa  = dist(*rng);
 119             beta  = dist(*rng);
 120             gamma = dist(*rng);
 121             break;
 122         case en_rotZ:
 123             alfa = beta = 0.;
 124             gamma       = dist(*rng);
 125             break;
 126         case en_rotNone: alfa = beta = gamma = 0.; break;
 127     }
 128     if (enum_rot == en_rotXYZ || enum_rot == en_rotZ)
 129     {
 130         rotate_conf(xout->size(), as_rvec_array(xout->data()), nullptr, alfa, beta, gamma);
 131     }
 132     for (size_t i = 0; i < xout->size(); ++i)
 133     {
 134         rvec_inc((*xout)[i], offset);
 135     }
 136 }
 137
 138 static bool isInsertionAllowed(gmx::AnalysisNeighborhoodSearch* search,
 139                                const std::vector<real>&         exclusionDistances,
 140                                const std::vector<RVec>&         x,
 141                                const std::vector<real>&         exclusionDistances_insrt,
 142                                const t_atoms&                   atoms,
 143                                const std::set<int>&             removableAtoms,
 144                                gmx::AtomsRemover*               remover)
 145 {
 146     gmx::AnalysisNeighborhoodPositions  pos(x);
 147     gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
 148     gmx::AnalysisNeighborhoodPair       pair;
 149     while (pairSearch.findNextPair(&pair))
 150     {
 151         const real r1 = exclusionDistances[pair.refIndex()];
 152         const real r2 = exclusionDistances_insrt[pair.testIndex()];
 153         if (pair.distance2() < gmx::square(r1 + r2))
 154         {
 155             if (removableAtoms.count(pair.refIndex()) == 0)
 156             {
 157                 return false;
 158             }
 159             // TODO: If molecule information is available, this should ideally
 160             // use it to remove whole molecules.
 161             remover->markResidue(atoms, pair.refIndex(), true);
 162         }
 163     }
 164     return true;
 165 }
 166
 167 static void insert_mols(int                  nmol_insrt,
 168                         int                  ntry,
 169                         int                  seed,
 170                         real                 defaultDistance,
 171                         real                 scaleFactor,
 172                         t_atoms*             atoms,
 173                         t_symtab*            symtab,
 174                         std::vector<RVec>*   x,
 175                         const std::set<int>& removableAtoms,
 176                         const t_atoms&       atoms_insrt,
 177                         gmx::ArrayRef<RVec>  x_insrt,
 178                         int                  ePBC,
 179                         matrix               box,
 180                         const std::string&   posfn,
 181                         const rvec           deltaR,
 182                         RotationType         enum_rot)
 183 {
 184     fprintf(stderr, "Initialising inter-atomic distances...\n");
 185     AtomProperties aps;
 186     std::vector<real> exclusionDistances(makeExclusionDistances(atoms, &aps, defaultDistance, scaleFactor));
 187     const std::vector<real> exclusionDistances_insrt(
 188             makeExclusionDistances(&atoms_insrt, &aps, defaultDistance, scaleFactor));
 189
 190     const real maxInsertRadius =
 191             *std::max_element(exclusionDistances_insrt.begin(), exclusionDistances_insrt.end());
 192     real maxRadius = maxInsertRadius;
 193     if (!exclusionDistances.empty())
 194     {
 195         const real maxExistingRadius =
 196                 *std::max_element(exclusionDistances.begin(), exclusionDistances.end());
 197         maxRadius = std::max(maxInsertRadius, maxExistingRadius);
 198     }
 199
 200     // TODO: Make all of this exception-safe.
 201     gmx::AnalysisNeighborhood nb;
 202     nb.setCutoff(maxInsertRadius + maxRadius);
 203
 204
 205     if (seed == 0)
 206     {
 207         seed = static_cast<int>(gmx::makeRandomSeed());
 208     }
 209     fprintf(stderr, "Using random seed %d\n", seed);
 210
 211     gmx::DefaultRandomEngine rng(seed);
 212
 213     t_pbc pbc;
 214     set_pbc(&pbc, ePBC, box);
 215
 216     /* With -ip, take nmol_insrt from file posfn */
 217     double**   rpos              = nullptr;
 218     const bool insertAtPositions = !posfn.empty();
 219     if (insertAtPositions)
 220     {
 221         int ncol;
 222         nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
 223         if (ncol != 3)
 224         {
 225             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n", posfn.c_str());
 226         }
 227         fprintf(stderr, "Read %d positions from file %s\n\n", nmol_insrt, posfn.c_str());
 228     }
 229
 230     gmx::AtomsBuilder builder(atoms, symtab);
 231     gmx::AtomsRemover remover(*atoms);
 232     {
 233         const int finalAtomCount    = atoms->nr + nmol_insrt * atoms_insrt.nr;
 234         const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt.nres;
 235         builder.reserve(finalAtomCount, finalResidueCount);
 236         x->reserve(finalAtomCount);
 237         exclusionDistances.reserve(finalAtomCount);
 238     }
 239
 240     std::vector<RVec> x_n(x_insrt.size());
 241
 242     int                                mol        = 0;
 243     int                                trial      = 0;
 244     int                                firstTrial = 0;
 245     int                                failed     = 0;
 246     gmx::UniformRealDistribution<real> dist;
 247
 248     while (mol < nmol_insrt && trial < ntry * nmol_insrt)
 249     {
 250         rvec offset_x;
 251         if (!insertAtPositions)
 252         {
 253             // Insert at random positions.
 254             offset_x[XX] = box[XX][XX] * dist(rng);
 255             offset_x[YY] = box[YY][YY] * dist(rng);
 256             offset_x[ZZ] = box[ZZ][ZZ] * dist(rng);
 257         }
 258         else
 259         {
 260             // Skip a position if ntry trials were not successful.
 261             if (trial >= firstTrial + ntry)
 262             {
 263                 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n", rpos[XX][mol],
 264                         rpos[YY][mol], rpos[ZZ][mol]);
 265                 ++mol;
 266                 ++failed;
 267                 firstTrial = trial;
 268                 continue;
 269             }
 270             // Insert at positions taken from option -ip file.
 271             offset_x[XX] = rpos[XX][mol] + deltaR[XX] * (2 * dist(rng) - 1);
 272             offset_x[YY] = rpos[YY][mol] + deltaR[YY] * (2 * dist(rng) - 1);
 273             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ] * (2 * dist(rng) - 1);
 274         }
 275         fprintf(stderr, "\rTry %d", ++trial);
 276         fflush(stderr);
 277
 278         generate_trial_conf(x_insrt, offset_x, enum_rot, &rng, &x_n);
 279         gmx::AnalysisNeighborhoodPositions pos(*x);
 280         gmx::AnalysisNeighborhoodSearch    search = nb.initSearch(&pbc, pos);
 281         if (isInsertionAllowed(&search, exclusionDistances, x_n, exclusionDistances_insrt, *atoms,
 282                                removableAtoms, &remover))
 283         {
 284             x->insert(x->end(), x_n.begin(), x_n.end());
 285             exclusionDistances.insert(exclusionDistances.end(), exclusionDistances_insrt.begin(),
 286                                       exclusionDistances_insrt.end());
 287             builder.mergeAtoms(atoms_insrt);
 288             ++mol;
 289             firstTrial = trial;
 290             fprintf(stderr, " success (now %d atoms)!\n", builder.currentAtomCount());
 291         }
 292     }
 293
 294     fprintf(stderr, "\n");
 295     /* print number of molecules added */
 296     fprintf(stderr, "Added %d molecules (out of %d requested)\n", mol - failed, nmol_insrt);
 297
 298     const int originalAtomCount    = atoms->nr;
 299     const int originalResidueCount = atoms->nres;
 300     remover.refreshAtomCount(*atoms);
 301     remover.removeMarkedElements(x);
 302     remover.removeMarkedAtoms(atoms);
 303     if (atoms->nr < originalAtomCount)
 304     {
 305         fprintf(stderr, "Replaced %d residues (%d atoms)\n", originalResidueCount - atoms->nres,
 306                 originalAtomCount - atoms->nr);
 307     }
 308
 309     if (rpos != nullptr)
 310     {
 311         for (int i = 0; i < DIM; ++i)
 312         {
 313             sfree(rpos[i]);
 314         }
 315         sfree(rpos);
 316     }
 317 }
 318
 319 namespace gmx
 320 {
 321
 322 namespace
 323 {
 324
 325 class InsertMolecules : public ICommandLineOptionsModule, public ITopologyProvider
 326 {
 327 public:
 328     InsertMolecules() :
 329         bBox_(false),
 330         nmolIns_(0),
 331         nmolTry_(10),
 332         seed_(0),
 333         defaultDistance_(0.105),
 334         scaleFactor_(0.57),
 335         enumRot_(en_rotXYZ)
 336     {
 337         clear_rvec(newBox_);
 338         clear_rvec(deltaR_);
 339         clear_mat(box_);
 340     }
 341
 342     // From ITopologyProvider
 343     gmx_mtop_t* getTopology(bool /*required*/) override { return &top_; }
 344     int         getAtomCount() override { return 0; }
 345
 346     // From ICommandLineOptionsModule
 347     void init(CommandLineModuleSettings* /*settings*/) override {}
 348     void initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings) override;
 349     void optionsFinished() override;
 350     int  run() override;
 351
 352 private:
 353     SelectionCollection selections_;
 354
 355     std::string  inputConfFile_;
 356     std::string  insertConfFile_;
 357     std::string  positionFile_;
 358     std::string  outputConfFile_;
 359     rvec         newBox_;
 360     bool         bBox_;
 361     int          nmolIns_;
 362     int          nmolTry_;
 363     int          seed_;
 364     real         defaultDistance_;
 365     real         scaleFactor_;
 366     rvec         deltaR_;
 367     RotationType enumRot_;
 368     Selection    replaceSel_;
 369
 370     gmx_mtop_t        top_;
 371     std::vector<RVec> x_;
 372     matrix            box_;
 373     int               ePBC_;
 374 };
 375
 376 void InsertMolecules::initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings)
 377 {
 378     const char* const desc[] = {
 379         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 380         "the [TT]-ci[tt] input file. The insertions take place either into",
 381         "vacant space in the solute conformation given with [TT]-f[tt], or",
 382         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 383         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 384         "around the solute before insertions. Any velocities present are",
 385         "discarded.",
 386         "",
 387         "It is possible to also insert into a solvated configuration and",
 388         "replace solvent atoms with the inserted atoms. To do this, use",
 389         "[TT]-replace[tt] to specify a selection that identifies the atoms",
 390         "that can be replaced. The tool assumes that all molecules in this",
 391         "selection consist of single residues: each residue from this",
 392         "selection that overlaps with the inserted molecules will be removed",
 393         "instead of preventing insertion.",
 394         "",
 395         "By default, the insertion positions are random (with initial seed",
 396         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 397         "molecules have been inserted in the box. Molecules are not inserted",
 398         "where the distance between any existing atom and any atom of the",
 399         "inserted molecule is less than the sum based on the van der Waals",
 400         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 401         "Waals radii is read by the program, and the resulting radii scaled",
 402         "by [TT]-scale[tt]. If radii are not found in the database, those",
 403         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
 404         "Note that the usefulness of those radii depends on the atom names,",
 405         "and thus varies widely with force field.",
 406         "",
 407         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 408         "before giving up. Increase [TT]-try[tt] if you have several small",
 409         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 410         "molecules are randomly oriented before insertion attempts.",
 411         "",
 412         "Alternatively, the molecules can be inserted only at positions defined in",
 413         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 414         "that give the displacements compared to the input molecule position",
 415         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 416         "positions, the molecule must be centered on (0,0,0) before using",
 417         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 418         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 419         "defines the maximally allowed displacements during insertial trials.",
 420         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
 421     };
 422
 423     settings->setHelpText(desc);
 424
 425     std::shared_ptr<SelectionOptionBehavior> selectionOptionBehavior(
 426             new SelectionOptionBehavior(&selections_, this));
 427     settings->addOptionsBehavior(selectionOptionBehavior);
 428
 429     // TODO: Replace use of legacyType.
 430     options->addOption(FileNameOption("f")
 431                                .legacyType(efSTX)
 432                                .inputFile()
 433                                .store(&inputConfFile_)
 434                                .defaultBasename("protein")
 435                                .description("Existing configuration to insert into"));
 436     options->addOption(FileNameOption("ci")
 437                                .legacyType(efSTX)
 438                                .inputFile()
 439                                .required()
 440                                .store(&insertConfFile_)
 441                                .defaultBasename("insert")
 442                                .description("Configuration to insert"));
 443     options->addOption(FileNameOption("ip")
 444                                .filetype(eftGenericData)
 445                                .inputFile()
 446                                .store(&positionFile_)
 447                                .defaultBasename("positions")
 448                                .description("Predefined insertion trial positions"));
 449     options->addOption(FileNameOption("o")
 450                                .legacyType(efSTO)
 451                                .outputFile()
 452                                .required()
 453                                .store(&outputConfFile_)
 454                                .defaultBasename("out")
 455                                .description("Output configuration after insertion"));
 456
 457     options->addOption(
 458             SelectionOption("replace").onlyAtoms().store(&replaceSel_).description("Atoms that can be removed if overlapping"));
 459     selectionOptionBehavior->initOptions(options);
 460
 461     options->addOption(RealOption("box").vector().store(newBox_).storeIsSet(&bBox_).description(
 462             "Box size (in nm)"));
 463     options->addOption(IntegerOption("nmol").store(&nmolIns_).description(
 464             "Number of extra molecules to insert"));
 465     options->addOption(IntegerOption("try").store(&nmolTry_).description(
 466             "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
 467     options->addOption(IntegerOption("seed").store(&seed_).description(
 468             "Random generator seed (0 means generate)"));
 469     options->addOption(
 470             RealOption("radius").store(&defaultDistance_).description("Default van der Waals distance"));
 471     options->addOption(
 472             RealOption("scale").store(&scaleFactor_).description("Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water."));
 473     options->addOption(RealOption("dr").vector().store(deltaR_).description(
 474             "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
 475     options->addOption(EnumOption<RotationType>("rot")
 476                                .enumValue(cRotationEnum)
 477                                .store(&enumRot_)
 478                                .description("Rotate inserted molecules randomly"));
 479 }
 480
 481 void InsertMolecules::optionsFinished()
 482 {
 483     if (nmolIns_ <= 0 && positionFile_.empty())
 484     {
 485         GMX_THROW(
 486                 InconsistentInputError("Either -nmol must be larger than 0, "
 487                                        "or positions must be given with -ip."));
 488     }
 489     if (inputConfFile_.empty() && !bBox_)
 490     {
 491         GMX_THROW(
 492                 InconsistentInputError("When no solute (-f) is specified, "
 493                                        "a box size (-box) must be specified."));
 494     }
 495     if (replaceSel_.isValid() && inputConfFile_.empty())
 496     {
 497         GMX_THROW(
 498                 InconsistentInputError("Replacement (-replace) only makes sense "
 499                                        "together with an existing configuration (-f)."));
 500     }
 501
 502     if (!inputConfFile_.empty())
 503     {
 504         bool  bTprFileWasRead;
 505         rvec* temporaryX = nullptr;
 506         fprintf(stderr, "Reading solute configuration\n");
 507         readConfAndTopology(inputConfFile_.c_str(), &bTprFileWasRead, &top_, &ePBC_, &temporaryX,
 508                             nullptr, box_);
 509         x_.assign(temporaryX, temporaryX + top_.natoms);
 510         sfree(temporaryX);
 511         if (top_.natoms == 0)
 512         {
 513             fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
 514         }
 515     }
 516 }
 517
 518 int InsertMolecules::run()
 519 {
 520     std::set<int> removableAtoms;
 521     if (replaceSel_.isValid())
 522     {
 523         t_pbc pbc;
 524         set_pbc(&pbc, ePBC_, box_);
 525         t_trxframe* fr;
 526         snew(fr, 1);
 527         fr->natoms = x_.size();
 528         fr->bX     = TRUE;
 529         fr->x      = as_rvec_array(x_.data());
 530         selections_.evaluate(fr, &pbc);
 531         sfree(fr);
 532         removableAtoms.insert(replaceSel_.atomIndices().begin(), replaceSel_.atomIndices().end());
 533         // TODO: It could be nice to check that removableAtoms contains full
 534         // residues, since we anyways remove whole residues instead of
 535         // individual atoms.
 536     }
 537
 538     int ePBCForOutput = ePBC_;
 539     if (bBox_)
 540     {
 541         ePBCForOutput = epbcXYZ;
 542         clear_mat(box_);
 543         box_[XX][XX] = newBox_[XX];
 544         box_[YY][YY] = newBox_[YY];
 545         box_[ZZ][ZZ] = newBox_[ZZ];
 546     }
 547     if (det(box_) == 0)
 548     {
 549         gmx_fatal(FARGS,
 550                   "Undefined solute box.\nCreate one with gmx editconf "
 551                   "or give explicit -box command line option");
 552     }
 553
 554     gmx_mtop_t        topInserted;
 555     t_atoms           atomsInserted;
 556     std::vector<RVec> xInserted;
 557     {
 558         bool   bTprFileWasRead;
 559         int    ePBC_dummy;
 560         matrix box_dummy;
 561         rvec*  temporaryX;
 562         readConfAndTopology(insertConfFile_.c_str(), &bTprFileWasRead, &topInserted, &ePBC_dummy,
 563                             &temporaryX, nullptr, box_dummy);
 564         xInserted.assign(temporaryX, temporaryX + topInserted.natoms);
 565         sfree(temporaryX);
 566         atomsInserted = gmx_mtop_global_atoms(&topInserted);
 567         if (atomsInserted.nr == 0)
 568         {
 569             gmx_fatal(FARGS, "No molecule in %s, please check your input", insertConfFile_.c_str());
 570         }
 571         if (top_.name == nullptr)
 572         {
 573             top_.name = topInserted.name;
 574         }
 575         if (positionFile_.empty())
 576         {
 577             center_molecule(xInserted);
 578         }
 579     }
 580
 581     t_atoms atoms = gmx_mtop_global_atoms(&top_);
 582
 583     /* add nmol_ins molecules of atoms_ins
 584        in random orientation at random place */
 585     insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_, &atoms, &top_.symtab,
 586                 &x_, removableAtoms, atomsInserted, xInserted, ePBCForOutput, box_, positionFile_,
 587                 deltaR_, enumRot_);
 588
 589     /* write new configuration to file confout */
 590     fprintf(stderr, "Writing generated configuration to %s\n", outputConfFile_.c_str());
 591     write_sto_conf(outputConfFile_.c_str(), *top_.name, &atoms, as_rvec_array(x_.data()), nullptr,
 592                    ePBCForOutput, box_);
 593
 594     /* print size of generated configuration */
 595     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n", atoms.nr, atoms.nres);
 596
 597     done_atom(&atoms);
 598     done_atom(&atomsInserted);
 599
 600     return 0;
 601 }
 602
 603 } // namespace
 604
 605
 606 const char* InsertMoleculesInfo::name()
 607 {
 608     static const char* name = "insert-molecules";
 609     return name;
 610 }
 611
 612 const char* InsertMoleculesInfo::shortDescription()
 613 {
 614     static const char* shortDescription = "Insert molecules into existing vacancies";
 615     return shortDescription;
 616 }
 617
 618 ICommandLineOptionsModulePointer InsertMoleculesInfo::create()
 619 {
 620     return ICommandLineOptionsModulePointer(new InsertMolecules());
 621 }
 622
 623 } // namespace gmx