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