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 "insert-molecules.h"
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43
  44 #include "sysstuff.h"
  45 #include "typedefs.h"
  46 #include "gromacs/utility/smalloc.h"
  47 #include "gromacs/math/utilities.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "macros.h"
  50 #include "gromacs/random/random.h"
  51 #include "gromacs/fileio/futil.h"
  52 #include "atomprop.h"
  53 #include "names.h"
  54 #include "vec.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/commandline/pargs.h"
  57 #include "gromacs/gmxlib/conformation-utilities.h"
  58 #include "addconf.h"
  59 #include "read-conformation.h"
  60 #include "pbc.h"
  61 #include "xvgr.h"
  62
  63 static gmx_bool in_box(t_pbc *pbc, rvec x)
  64 {
  65     rvec box_center, dx;
  66     int  shift;
  67
  68     /* pbc_dx_aiuc only works correctly with the rectangular box center */
  69     calc_box_center(ecenterRECT, pbc->box, box_center);
  70
  71     shift = pbc_dx_aiuc(pbc, x, box_center, dx);
  72
  73     return (shift == CENTRAL);
  74 }
  75
  76 /* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
  77  * leaks memory (May 2012). The function could be deleted as soon as the memory leaks
  78  * there are fixed.
  79  * However, when inserting a small molecule in a system containing not too many atoms,
  80  * allPairsDistOk is probably even faster than the other code.
  81  */
  82 static gmx_bool
  83 allPairsDistOk(t_atoms *atoms, rvec *x, real *exclusionDistances,
  84                int ePBC, matrix box,
  85                t_atoms *atoms_insrt, rvec *x_n, real *exclusionDistances_insrt)
  86 {
  87     int   i, j;
  88     rvec  dx;
  89     real  n2, r2;
  90     t_pbc pbc;
  91
  92     set_pbc(&pbc, ePBC, box);
  93     for (i = 0; i < atoms->nr; i++)
  94     {
  95         for (j = 0; j < atoms_insrt->nr; j++)
  96         {
  97             pbc_dx(&pbc, x[i], x_n[j], dx);
  98             n2 = norm2(dx);
  99             r2 = sqr(exclusionDistances[i]+exclusionDistances_insrt[j]);
 100             if (n2 < r2)
 101             {
 102                 return FALSE;
 103             }
 104         }
 105     }
 106     return TRUE;
 107 }
 108
 109 /* enum for random rotations of inserted solutes */
 110 enum {
 111     en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR
 112 };
 113
 114 static char *insert_mols(const char *mol_insrt, int nmol_insrt, int ntry, int seed,
 115                          t_atoms *atoms, rvec **x, real **exclusionDistances, int ePBC, matrix box,
 116                          gmx_atomprop_t aps,
 117                          real defaultDistance, real scaleFactor, real rshell,
 118                          const output_env_t oenv,
 119                          const char* posfn, const rvec deltaR, int enum_rot,
 120                          gmx_bool bCheckAllPairDist)
 121 {
 122     t_pbc            pbc;
 123     static  char    *title_insrt;
 124     t_atoms          atoms_insrt;
 125     rvec            *x_insrt, *x_n;
 126     real            *exclusionDistances_insrt;
 127     int              ePBC_insrt;
 128     matrix           box_insrt;
 129     int              i, mol, onr, ncol;
 130     real             alfa = 0., beta = 0., gamma = 0.;
 131     rvec             offset_x;
 132     int              trial;
 133     double         **rpos;
 134     gmx_rng_t        rng;
 135
 136     rng = gmx_rng_init(seed);
 137     set_pbc(&pbc, ePBC, box);
 138
 139     /* read number of atoms of insert molecules */
 140     {
 141         int natoms;
 142         get_stx_coordnum(mol_insrt, &natoms);
 143         if (natoms == 0)
 144         {
 145             gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
 146         }
 147         init_t_atoms(&atoms_insrt, natoms, FALSE);
 148     }
 149     /* allocate memory for atom coordinates of insert molecules */
 150     snew(x_insrt, atoms_insrt.nr);
 151     snew(atoms_insrt.resinfo, atoms_insrt.nr);
 152     snew(atoms_insrt.atomname, atoms_insrt.nr);
 153     snew(atoms_insrt.atom, atoms_insrt.nr);
 154     atoms_insrt.pdbinfo = NULL;
 155     snew(x_n, atoms_insrt.nr);
 156     snew(title_insrt, STRLEN);
 157
 158     /* read residue number, residue names, atomnames, coordinates etc. */
 159     fprintf(stderr, "Reading molecule configuration \n");
 160     read_stx_conf(mol_insrt, title_insrt, &atoms_insrt, x_insrt, NULL,
 161                   &ePBC_insrt, box_insrt);
 162     fprintf(stderr, "%s\nContaining %d atoms in %d residue\n",
 163             title_insrt, atoms_insrt.nr, atoms_insrt.nres);
 164     srenew(atoms_insrt.resinfo, atoms_insrt.nres);
 165
 166     /* initialise distance arrays for inserted molecules */
 167     exclusionDistances_insrt = makeExclusionDistances(&atoms_insrt, aps, defaultDistance, scaleFactor);
 168
 169     /* With -ip, take nmol_insrt from file posfn */
 170     if (posfn != NULL)
 171     {
 172         nmol_insrt = read_xvg(posfn, &rpos, &ncol);
 173         if (ncol != 3)
 174         {
 175             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n", ncol, posfn);
 176         }
 177         fprintf(stderr, "Read %d positions from file %s\n\n", nmol_insrt, posfn);
 178     }
 179
 180     srenew(atoms->resinfo, (atoms->nres+nmol_insrt*atoms_insrt.nres));
 181     srenew(atoms->atomname, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 182     srenew(atoms->atom, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 183     srenew(*x, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 184     srenew(*exclusionDistances, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 185
 186     trial = mol = 0;
 187     while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
 188     {
 189         fprintf(stderr, "\rTry %d", trial++);
 190         for (i = 0; (i < atoms_insrt.nr); i++)
 191         {
 192             copy_rvec(x_insrt[i], x_n[i]);
 193         }
 194         switch (enum_rot)
 195         {
 196             case en_rotXYZ:
 197                 alfa  = 2*M_PI * gmx_rng_uniform_real(rng);
 198                 beta  = 2*M_PI * gmx_rng_uniform_real(rng);
 199                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 200                 break;
 201             case en_rotZ:
 202                 alfa  = beta = 0.;
 203                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 204                 break;
 205             case en_rotNone:
 206                 alfa = beta = gamma = 0.;
 207                 break;
 208         }
 209         if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
 210         {
 211             rotate_conf(atoms_insrt.nr, x_n, NULL, alfa, beta, gamma);
 212         }
 213         if (posfn == NULL)
 214         {
 215             /* insert at random positions */
 216             offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
 217             offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
 218             offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
 219             make_new_box(atoms_insrt.nr, x_n, box_insrt, offset_x, TRUE);
 220             if (!in_box(&pbc, x_n[0]) || !in_box(&pbc, x_n[atoms_insrt.nr-1]))
 221             {
 222                 continue;
 223             }
 224         }
 225         else
 226         {
 227             /* Insert at positions taken from option -ip file */
 228             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
 229             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
 230             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
 231             for (i = 0; i < atoms_insrt.nr; i++)
 232             {
 233                 rvec_inc(x_n[i], offset_x);
 234             }
 235         }
 236
 237         onr = atoms->nr;
 238
 239         /* This is a (maybe) slow workaround to avoid too many calls of add_conf, which
 240          * leaks memory (status May 2012). If the momory leaks in add_conf() are fixed,
 241          * this check could be removed. Note, however, that allPairsDistOk is probably
 242          * even faster than add_conf() when inserting a small molecule into a moderately
 243          * small system.
 244          */
 245         if (bCheckAllPairDist && !allPairsDistOk(atoms, *x, *exclusionDistances, ePBC, box, &atoms_insrt, x_n, exclusionDistances_insrt))
 246         {
 247             continue;
 248         }
 249
 250         add_conf(atoms, x, NULL, exclusionDistances, FALSE, ePBC, box, TRUE,
 251                  &atoms_insrt, x_n, NULL, exclusionDistances_insrt, FALSE, rshell, 0, oenv);
 252
 253         if (atoms->nr == (atoms_insrt.nr+onr))
 254         {
 255             mol++;
 256             fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
 257         }
 258     }
 259     gmx_rng_destroy(rng);
 260     srenew(atoms->resinfo,  atoms->nres);
 261     srenew(atoms->atomname, atoms->nr);
 262     srenew(atoms->atom,     atoms->nr);
 263     srenew(*x,              atoms->nr);
 264     srenew(*exclusionDistances, atoms->nr);
 265
 266     fprintf(stderr, "\n");
 267     /* print number of molecules added */
 268     fprintf(stderr, "Added %d molecules (out of %d requested) of %s\n",
 269             mol, nmol_insrt, *atoms_insrt.resinfo[0].name);
 270
 271     sfree(x_n);
 272     sfree(exclusionDistances_insrt);
 273     done_atom(&atoms_insrt);
 274
 275     return title_insrt;
 276 }
 277
 278 int gmx_insert_molecules(int argc, char *argv[])
 279 {
 280     const char *desc[] = {
 281         "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
 282         "the [TT]-ci[tt] input file. The insertions take place either into",
 283         "vacant space in the solute conformation given with [TT]-f[tt], or",
 284         "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
 285         "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
 286         "around the solute before insertions. Any velocities present are",
 287         "discarded.[PAR]",
 288
 289         "By default, the insertion positions are random (with initial seed",
 290         "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
 291         "molecules have been inserted in the box. Molecules are not inserted",
 292         "where the distance between any existing atom and any atom of the",
 293         "inserted molecule is less than the sum based on the van der Waals",
 294         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 295         "Waals radii is read by the program, and the resulting radii scaled",
 296         "by [TT]-scale[tt]. If radii are not found in the database, those"
 297         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 298
 299         "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
 300         "before giving up. Increase [TT]-try[tt] if you have several small",
 301         "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
 302         "molecules are randomly oriented before insertion attempts.[PAR]",
 303
 304         "Alternatively, the molecules can be inserted only at positions defined in",
 305         "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
 306         "that give the displacements compared to the input molecule position",
 307         "([TT]-ci[tt]). Hence, if that file should contain the absolute",
 308         "positions, the molecule must be centered on (0,0,0) before using",
 309         "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
 310         "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
 311         "defines the maximally allowed displacements during insertial trials.",
 312         "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above).",
 313         "[PAR]",
 314     };
 315
 316     const char *bugs[] = {
 317         "Molecules must be whole in the initial configurations.",
 318         "Many repeated neighbor searchings with -ci blows up the allocated memory. "
 319         "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
 320     };
 321
 322     /* parameter data */
 323     gmx_bool       bProt, bBox;
 324     const char    *conf_prot, *confout;
 325     real          *exclusionDistances = NULL;
 326     char          *title_ins          = NULL;
 327     gmx_atomprop_t aps;
 328
 329     /* protein configuration data */
 330     char          *title = NULL;
 331     t_atoms       *atoms;
 332     rvec          *x    = NULL;
 333     int            ePBC = -1;
 334     matrix         box;
 335
 336     t_filenm       fnm[] = {
 337         { efSTX, "-f", "protein", ffOPTRD },
 338         { efSTX, "-ci", "insert",  ffREAD},
 339         { efDAT, "-ip", "positions",  ffOPTRD},
 340         { efSTO, NULL,  NULL,      ffWRITE},
 341     };
 342 #define NFILE asize(fnm)
 343
 344     static int      nmol_ins               = 0, nmol_try = 10, seed = 1997, enum_rot;
 345     static real     defaultDistance        = 0.105, scaleFactor = 0.57;
 346     static rvec     new_box                = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
 347     static gmx_bool bCheckAllPairDist      = FALSE;
 348     output_env_t    oenv;
 349     const char     *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
 350     t_pargs         pa[]              = {
 351         { "-box",    FALSE, etRVEC, {new_box},
 352           "Box size (in nm)" },
 353         { "-nmol",   FALSE, etINT, {&nmol_ins},
 354           "Number of extra molecules to insert" },
 355         { "-try",    FALSE, etINT, {&nmol_try},
 356           "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
 357         { "-seed",   FALSE, etINT, {&seed},
 358           "Random generator seed"},
 359         { "-radius",   FALSE, etREAL, {&defaultDistance},
 360           "Default van der Waals distance"},
 361         { "-scale", FALSE, etREAL, {&scaleFactor},
 362           "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." },
 363         { "-dr",    FALSE, etRVEC, {deltaR},
 364           "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
 365         { "-rot", FALSE,  etENUM, {enum_rot_string},
 366           "rotate inserted molecules randomly" },
 367         { "-allpair",    FALSE, etBOOL, {&bCheckAllPairDist},
 368           "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
 369     };
 370
 371     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 372                            asize(desc), desc, asize(bugs), bugs, &oenv))
 373     {
 374         return 0;
 375     }
 376
 377     bProt     = opt2bSet("-f", NFILE, fnm);
 378     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 379     enum_rot  = nenum(enum_rot_string);
 380
 381     /* check input */
 382     const char *insertionMoleculeFileName = opt2fn("-ci", NFILE, fnm);
 383     if (!gmx_fexist(insertionMoleculeFileName))
 384     {
 385         gmx_fatal(FARGS,
 386                   "A molecule conformation to insert is required in -ci. %s was not found!",
 387                   insertionMoleculeFileName);
 388     }
 389     if (nmol_ins <= 0 && !opt2bSet("-ip", NFILE, fnm))
 390     {
 391         gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
 392                   "or positions must be given with -ip");
 393     }
 394     if (!bProt && !bBox)
 395     {
 396         gmx_fatal(FARGS, "When no solute (-f) is specified, "
 397                   "a box size (-box) must be specified");
 398     }
 399
 400     aps = gmx_atomprop_init();
 401
 402     snew(atoms, 1);
 403     init_t_atoms(atoms, 0, FALSE);
 404     if (bProt)
 405     {
 406         /* Generate a solute configuration */
 407         conf_prot = opt2fn("-f", NFILE, fnm);
 408         title     = readConformation(conf_prot, atoms, &x, NULL,
 409                                      &ePBC, box);
 410         exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
 411         if (atoms->nr == 0)
 412         {
 413             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 414             bProt = FALSE;
 415         }
 416     }
 417     if (bBox)
 418     {
 419         ePBC = epbcXYZ;
 420         clear_mat(box);
 421         box[XX][XX] = new_box[XX];
 422         box[YY][YY] = new_box[YY];
 423         box[ZZ][ZZ] = new_box[ZZ];
 424     }
 425     if (det(box) == 0)
 426     {
 427         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 428                   "or give explicit -box command line option");
 429     }
 430
 431     /* add nmol_ins molecules of atoms_ins
 432        in random orientation at random place */
 433     title_ins = insert_mols(insertionMoleculeFileName, nmol_ins, nmol_try, seed,
 434                             atoms, &x, &exclusionDistances, ePBC, box, aps,
 435                             defaultDistance, scaleFactor, 0,
 436                             oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
 437                             bCheckAllPairDist);
 438
 439     /* write new configuration to file confout */
 440     confout = ftp2fn(efSTO, NFILE, fnm);
 441     fprintf(stderr, "Writing generated configuration to %s\n", confout);
 442     if (bProt)
 443     {
 444         write_sto_conf(confout, title, atoms, x, NULL, ePBC, box);
 445         /* print box sizes and box type to stderr */
 446         fprintf(stderr, "%s\n", title);
 447     }
 448     else
 449     {
 450         write_sto_conf(confout, title_ins, atoms, x, NULL, ePBC, box);
 451     }
 452
 453     /* print size of generated configuration */
 454     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 455             atoms->nr, atoms->nres);
 456
 457     gmx_atomprop_destroy(aps);
 458     sfree(exclusionDistances);
 459     sfree(x);
 460     done_atom(atoms);
 461     sfree(atoms);
 462
 463     return 0;
 464 }