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