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, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "insert-molecules.h"
  40
  41 #include <algorithm>
  42 #include <string>
  43
  44 #include "gromacs/commandline/cmdlineoptionsmodule.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/filenm.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxlib/conformation-utilities.h"
  49 #include "gromacs/gmxpreprocess/read-conformation.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/options/basicoptions.h"
  53 #include "gromacs/options/filenameoption.h"
  54 #include "gromacs/options/ioptionscontainer.h"
  55 #include "gromacs/options/options.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/random/random.h"
  58 #include "gromacs/selection/nbsearch.h"
  59 #include "gromacs/topology/atomprop.h"
  60 #include "gromacs/topology/atoms.h"
  61 #include "gromacs/utility/cstringutil.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/smalloc.h"
  64
  65 /* enum for random rotations of inserted solutes */
  66 enum {
  67     en_rotXYZ, en_rotZ, en_rotNone
  68 };
  69
  70 static void center_molecule(int atomCount, rvec x[])
  71 {
  72     rvec center;
  73     clear_rvec(center);
  74     for (int i = 0; i < atomCount; ++i)
  75     {
  76         rvec_inc(center, x[i]);
  77     }
  78     svmul(1.0/atomCount, center, center);
  79     for (int i = 0; i < atomCount; ++i)
  80     {
  81         rvec_dec(x[i], center);
  82     }
  83 }
  84
  85 static void generate_trial_conf(int atomCount, const rvec xin[],
  86                                 const rvec offset, int enum_rot, gmx_rng_t rng,
  87                                 rvec xout[])
  88 {
  89     for (int i = 0; i < atomCount; ++i)
  90     {
  91         copy_rvec(xin[i], xout[i]);
  92     }
  93     real alfa = 0.0, beta = 0.0, gamma = 0.0;
  94     switch (enum_rot)
  95     {
  96         case en_rotXYZ:
  97             alfa  = 2*M_PI * gmx_rng_uniform_real(rng);
  98             beta  = 2*M_PI * gmx_rng_uniform_real(rng);
  99             gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 100             break;
 101         case en_rotZ:
 102             alfa  = beta = 0.;
 103             gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 104             break;
 105         case en_rotNone:
 106             alfa = beta = gamma = 0.;
 107             break;
 108     }
 109     if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
 110     {
 111         rotate_conf(atomCount, xout, NULL, alfa, beta, gamma);
 112     }
 113     for (int i = 0; i < atomCount; ++i)
 114     {
 115         rvec_inc(xout[i], offset);
 116     }
 117 }
 118
 119 static bool is_insertion_allowed(gmx::AnalysisNeighborhoodSearch *search,
 120                                  const real *exclusionDistances,
 121                                  int atomCount, const rvec *x,
 122                                  const real *exclusionDistances_insrt)
 123 {
 124     gmx::AnalysisNeighborhoodPositions  pos(x, atomCount);
 125     gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
 126     gmx::AnalysisNeighborhoodPair       pair;
 127     while (pairSearch.findNextPair(&pair))
 128     {
 129         const real r1 = exclusionDistances[pair.refIndex()];
 130         const real r2 = exclusionDistances_insrt[pair.testIndex()];
 131         if (pair.distance2() < sqr(r1 + r2))
 132         {
 133             return false;
 134         }
 135     }
 136     return true;
 137 }
 138
 139 static void merge_atoms_noalloc(t_atoms *atoms, const t_atoms *atoms_add)
 140 {
 141     int resnr = 0;
 142     if (atoms->nr > 0)
 143     {
 144         resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
 145     }
 146     int prevResInd = -1;
 147     for (int i = 0; i < atoms_add->nr; ++i)
 148     {
 149         if (atoms_add->atom[i].resind != prevResInd)
 150         {
 151             prevResInd = atoms_add->atom[i].resind;
 152             ++resnr;
 153             atoms->resinfo[atoms->nres]    = atoms_add->resinfo[prevResInd];
 154             atoms->resinfo[atoms->nres].nr = resnr;
 155             ++atoms->nres;
 156         }
 157         atoms->atom[atoms->nr]        = atoms_add->atom[i];
 158         atoms->atomname[atoms->nr]    = atoms_add->atomname[i];
 159         atoms->atom[atoms->nr].resind = atoms->nres-1;
 160         ++atoms->nr;
 161     }
 162 }
 163
 164 static void insert_mols(int nmol_insrt, int ntry, int seed,
 165                         real defaultDistance, real scaleFactor,
 166                         t_atoms *atoms, rvec **x,
 167                         const t_atoms *atoms_insrt, const rvec *x_insrt,
 168                         int ePBC, matrix box,
 169                         const std::string &posfn, const rvec deltaR, int enum_rot)
 170 {
 171     t_pbc            pbc;
 172     rvec            *x_n;
 173
 174     fprintf(stderr, "Initialising inter-atomic distances...\n");
 175     gmx_atomprop_t   aps = gmx_atomprop_init();
 176     real            *exclusionDistances
 177         = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
 178     real            *exclusionDistances_insrt
 179         = makeExclusionDistances(atoms_insrt, aps, defaultDistance, scaleFactor);
 180     gmx_atomprop_destroy(aps);
 181
 182     const real       maxInsertRadius
 183         = *std::max_element(exclusionDistances_insrt,
 184                             exclusionDistances_insrt + atoms_insrt->nr);
 185     real             maxRadius = maxInsertRadius;
 186     if (atoms->nr > 0)
 187     {
 188         const real maxExistingRadius
 189             = *std::max_element(exclusionDistances,
 190                                 exclusionDistances + atoms->nr);
 191         maxRadius = std::max(maxInsertRadius, maxExistingRadius);
 192     }
 193
 194     // TODO: Make all of this exception-safe.
 195     gmx::AnalysisNeighborhood nb;
 196     nb.setCutoff(maxInsertRadius + maxRadius);
 197
 198     gmx_rng_t        rng = gmx_rng_init(seed);
 199     set_pbc(&pbc, ePBC, box);
 200
 201     snew(x_n, atoms_insrt->nr);
 202
 203     /* With -ip, take nmol_insrt from file posfn */
 204     double         **rpos = NULL;
 205     if (!posfn.empty())
 206     {
 207         int ncol;
 208         nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
 209         if (ncol != 3)
 210         {
 211             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n",
 212                       posfn.c_str());
 213         }
 214         fprintf(stderr, "Read %d positions from file %s\n\n",
 215                 nmol_insrt, posfn.c_str());
 216     }
 217
 218     {
 219         const int finalAtomCount    = atoms->nr + nmol_insrt * atoms_insrt->nr;
 220         const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt->nres;
 221         srenew(atoms->resinfo,      finalResidueCount);
 222         srenew(atoms->atomname,     finalAtomCount);
 223         srenew(atoms->atom,         finalAtomCount);
 224         srenew(*x,                  finalAtomCount);
 225         srenew(exclusionDistances,  finalAtomCount);
 226     }
 227
 228     int mol        = 0;
 229     int trial      = 0;
 230     int firstTrial = 0;
 231     int failed     = 0;
 232     while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
 233     {
 234         rvec offset_x;
 235         if (posfn.empty())
 236         {
 237             // Insert at random positions.
 238             offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
 239             offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
 240             offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
 241         }
 242         else
 243         {
 244             // Skip a position if ntry trials were not successful.
 245             if (trial >= firstTrial + ntry)
 246             {
 247                 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n",
 248                         rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
 249                 ++mol;
 250                 ++failed;
 251             }
 252             // Insert at positions taken from option -ip file.
 253             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
 254             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
 255             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
 256         }
 257         fprintf(stderr, "\rTry %d", ++trial);
 258         generate_trial_conf(atoms_insrt->nr, x_insrt, offset_x, enum_rot, rng,
 259                             x_n);
 260         gmx::AnalysisNeighborhoodPositions pos(*x, atoms->nr);
 261         gmx::AnalysisNeighborhoodSearch    search = nb.initSearch(&pbc, pos);
 262         if (is_insertion_allowed(&search, exclusionDistances, atoms_insrt->nr,
 263                                  x_n, exclusionDistances_insrt))
 264         {
 265             const int firstIndex = atoms->nr;
 266             for (int i = 0; i < atoms_insrt->nr; ++i)
 267             {
 268                 copy_rvec(x_n[i], (*x)[firstIndex + i]);
 269                 exclusionDistances[firstIndex + i] = exclusionDistances_insrt[i];
 270             }
 271             merge_atoms_noalloc(atoms, atoms_insrt);
 272             ++mol;
 273             firstTrial = trial;
 274             fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
 275         }
 276     }
 277     gmx_rng_destroy(rng);
 278     srenew(atoms->resinfo,  atoms->nres);
 279     srenew(atoms->atomname, atoms->nr);
 280     srenew(atoms->atom,     atoms->nr);
 281     srenew(*x,              atoms->nr);
 282
 283     fprintf(stderr, "\n");
 284     /* print number of molecules added */
 285     fprintf(stderr, "Added %d molecules (out of %d requested)\n",
 286             mol - failed, nmol_insrt);
 287
 288     sfree(x_n);
 289     sfree(exclusionDistances);
 290     sfree(exclusionDistances_insrt);
 291     if (rpos != NULL)
 292     {
 293         for (int i = 0; i < DIM; ++i)
 294         {
 295             sfree(rpos[i]);
 296         }
 297         sfree(rpos);
 298     }
 299 }
 300
 301 namespace gmx
 302 {
 303
 304 namespace
 305 {
 306
 307 class InsertMolecules : public ICommandLineOptionsModule
 308 {
 309     public:
 310         InsertMolecules()
 311             : bBox_(false), nmolIns_(0), nmolTry_(10), seed_(1997),
 312               defaultDistance_(0.105), scaleFactor_(0.57), enumRot_(en_rotXYZ)
 313         {
 314             clear_rvec(newBox_);
 315             clear_rvec(deltaR_);
 316         }
 317
 318         virtual void init(CommandLineModuleSettings * /*settings*/)
 319         {
 320         }
 321
 322         virtual void initOptions(IOptionsContainer                 *options,
 323                                  ICommandLineOptionsModuleSettings *settings);
 324         virtual void optionsFinished(Options *options);
 325
 326         virtual int run();
 327
 328     private:
 329         std::string inputConfFile_;
 330         std::string insertConfFile_;
 331         std::string positionFile_;
 332         std::string outputConfFile_;
 333         rvec        newBox_;
 334         bool        bBox_;
 335         int         nmolIns_;
 336         int         nmolTry_;
 337         int         seed_;
 338         real        defaultDistance_;
 339         real        scaleFactor_;
 340         rvec        deltaR_;
 341         int         enumRot_;
 342 };
 343
 344 void InsertMolecules::initOptions(IOptionsContainer                 *options,
 345                                   ICommandLineOptionsModuleSettings *settings)
 346 {
 347     const char *const desc[] = {
 348         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 349         "the [TT]-ci[tt] input file. The insertions take place either into",
 350         "vacant space in the solute conformation given with [TT]-f[tt], or",
 351         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 352         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 353         "around the solute before insertions. Any velocities present are",
 354         "discarded.[PAR]",
 355
 356         "By default, the insertion positions are random (with initial seed",
 357         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 358         "molecules have been inserted in the box. Molecules are not inserted",
 359         "where the distance between any existing atom and any atom of the",
 360         "inserted molecule is less than the sum based on the van der Waals",
 361         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 362         "Waals radii is read by the program, and the resulting radii scaled",
 363         "by [TT]-scale[tt]. If radii are not found in the database, those"
 364         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 365
 366         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 367         "before giving up. Increase [TT]-try[tt] if you have several small",
 368         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 369         "molecules are randomly oriented before insertion attempts.[PAR]",
 370
 371         "Alternatively, the molecules can be inserted only at positions defined in",
 372         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 373         "that give the displacements compared to the input molecule position",
 374         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 375         "positions, the molecule must be centered on (0,0,0) before using",
 376         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 377         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 378         "defines the maximally allowed displacements during insertial trials.",
 379         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
 380     };
 381
 382     settings->setHelpText(desc);
 383
 384     // TODO: Replace use of legacyType.
 385     options->addOption(FileNameOption("f")
 386                            .legacyType(efSTX).inputFile()
 387                            .store(&inputConfFile_)
 388                            .defaultBasename("protein")
 389                            .description("Existing configuration to insert into"));
 390     options->addOption(FileNameOption("ci")
 391                            .legacyType(efSTX).inputFile().required()
 392                            .store(&insertConfFile_)
 393                            .defaultBasename("insert")
 394                            .description("Configuration to insert"));
 395     options->addOption(FileNameOption("ip")
 396                            .filetype(eftGenericData).inputFile()
 397                            .store(&positionFile_)
 398                            .defaultBasename("positions")
 399                            .description("Predefined insertion trial positions"));
 400     options->addOption(FileNameOption("o")
 401                            .legacyType(efSTO).outputFile().required()
 402                            .store(&outputConfFile_)
 403                            .defaultBasename("out")
 404                            .description("Output configuration after insertion"));
 405
 406     options->addOption(RealOption("box").vector()
 407                            .store(newBox_)
 408                            .description("Box size (in nm)"));
 409     options->addOption(IntegerOption("nmol")
 410                            .store(&nmolIns_)
 411                            .description("Number of extra molecules to insert"));
 412     options->addOption(IntegerOption("try")
 413                            .store(&nmolTry_)
 414                            .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
 415     options->addOption(IntegerOption("seed")
 416                            .store(&seed_)
 417                            .description("Random generator seed"));
 418     options->addOption(RealOption("radius")
 419                            .store(&defaultDistance_)
 420                            .description("Default van der Waals distance"));
 421     options->addOption(RealOption("scale")
 422                            .store(&scaleFactor_)
 423                            .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."));
 424     options->addOption(RealOption("dr").vector()
 425                            .store(deltaR_)
 426                            .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
 427     const char *const cRotationEnum[] = {"xyz", "z", "none"};
 428     options->addOption(StringOption("rot").enumValue(cRotationEnum)
 429                            .storeEnumIndex(&enumRot_)
 430                            .description("Rotate inserted molecules randomly"));
 431 }
 432
 433 void InsertMolecules::optionsFinished(Options *options)
 434 {
 435     bBox_ = options->isSet("box");
 436 }
 437
 438 int InsertMolecules::run()
 439 {
 440     const bool bProt = !inputConfFile_.empty();
 441
 442     /* check input */
 443     if (nmolIns_ <= 0 && positionFile_.empty())
 444     {
 445         gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
 446                   "or positions must be given with -ip");
 447     }
 448     if (!bProt && !bBox_)
 449     {
 450         gmx_fatal(FARGS, "When no solute (-f) is specified, "
 451                   "a box size (-box) must be specified");
 452     }
 453
 454     char    *title = NULL;
 455     t_atoms *atoms;
 456     rvec    *x = NULL;
 457     matrix   box;
 458     int      ePBC = -1;
 459     snew(atoms, 1);
 460     init_t_atoms(atoms, 0, FALSE);
 461     if (bProt)
 462     {
 463         /* Generate a solute configuration */
 464         title = readConformation(inputConfFile_.c_str(), atoms, &x, NULL,
 465                                  &ePBC, box, "solute");
 466         if (atoms->nr == 0)
 467         {
 468             fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
 469             sfree(title);
 470             title = NULL;
 471         }
 472     }
 473     if (bBox_)
 474     {
 475         ePBC = epbcXYZ;
 476         clear_mat(box);
 477         box[XX][XX] = newBox_[XX];
 478         box[YY][YY] = newBox_[YY];
 479         box[ZZ][ZZ] = newBox_[ZZ];
 480     }
 481     if (det(box) == 0)
 482     {
 483         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 484                   "or give explicit -box command line option");
 485     }
 486
 487     t_atoms *atoms_insrt;
 488     rvec    *x_insrt = NULL;
 489     snew(atoms_insrt, 1);
 490     init_t_atoms(atoms_insrt, 0, FALSE);
 491     {
 492         int         ePBC_dummy;
 493         matrix      box_dummy;
 494         char       *title_ins
 495             = readConformation(insertConfFile_.c_str(), atoms_insrt, &x_insrt,
 496                                NULL, &ePBC_dummy, box_dummy, "molecule");
 497         if (atoms_insrt->nr == 0)
 498         {
 499             gmx_fatal(FARGS, "No molecule in %s, please check your input",
 500                       insertConfFile_.c_str());
 501         }
 502         if (title == NULL)
 503         {
 504             title = title_ins;
 505         }
 506         else
 507         {
 508             sfree(title_ins);
 509         }
 510         if (positionFile_.empty())
 511         {
 512             center_molecule(atoms_insrt->nr, x_insrt);
 513         }
 514     }
 515
 516     /* add nmol_ins molecules of atoms_ins
 517        in random orientation at random place */
 518     insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
 519                 atoms, &x, atoms_insrt, x_insrt,
 520                 ePBC, box, positionFile_, deltaR_, enumRot_);
 521
 522     /* write new configuration to file confout */
 523     fprintf(stderr, "Writing generated configuration to %s\n",
 524             outputConfFile_.c_str());
 525     write_sto_conf(outputConfFile_.c_str(), title, atoms, x, NULL, ePBC, box);
 526
 527     /* print size of generated configuration */
 528     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 529             atoms->nr, atoms->nres);
 530
 531     sfree(x);
 532     sfree(x_insrt);
 533     done_atom(atoms);
 534     done_atom(atoms_insrt);
 535     sfree(atoms);
 536     sfree(atoms_insrt);
 537     sfree(title);
 538
 539     return 0;
 540 }
 541
 542 }   // namespace
 543
 544 const char InsertMoleculesInfo::name[]             = "insert-molecules";
 545 const char InsertMoleculesInfo::shortDescription[] =
 546     "Insert molecules into existing vacancies";
 547 ICommandLineOptionsModule *InsertMoleculesInfo::create()
 548 {
 549     return new InsertMolecules();
 550 }
 551
 552 } // namespace gmx