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