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