src/gromacs/gmxpreprocess/genconf.c

   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 "genconf.h"
  40
  41 #include "gromacs/math/utilities.h"
  42 #include "gromacs/legacyheaders/macros.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/legacyheaders/txtdump.h"
  45 #include "gromacs/legacyheaders/readinp.h"
  46 #include "gromacs/legacyheaders/names.h"
  47 #include "sortwater.h"
  48 #include "gromacs/fileio/trxio.h"
  49
  50 #include "gromacs/commandline/pargs.h"
  51 #include "gromacs/math/3dtransforms.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/random/random.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 static void rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[],
  58                      gmx_rng_t rng, rvec max_rot)
  59 {
  60     mat4 mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
  61     rvec xcm;
  62     real phi;
  63     int  i, m;
  64
  65     clear_rvec(xcm);
  66     for (i = 0; (i < natoms); i++)
  67     {
  68         for (m = 0; (m < DIM); m++)
  69         {
  70             xcm[m] += x[i][m]/natoms; /* get center of mass of one molecule  */
  71         }
  72     }
  73     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  74
  75     /* move c.o.ma to origin */
  76     gmx_mat4_init_translation(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1);
  77     for (m = 0; (m < DIM); m++)
  78     {
  79         phi = M_PI*max_rot[m]*(2*gmx_rng_uniform_real(rng) - 1)/180;
  80         gmx_mat4_init_rotation(m, phi, mr[m]);
  81     }
  82     gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  83
  84     /* For gmx_mat4_mmul() we need to multiply in the opposite order
  85      * compared to normal mathematical notation.
  86      */
  87     gmx_mat4_mmul(mtemp1, mt1, mr[XX]);
  88     gmx_mat4_mmul(mtemp2, mr[YY], mr[ZZ]);
  89     gmx_mat4_mmul(mtemp3, mtemp1, mtemp2);
  90     gmx_mat4_mmul(mxtot, mtemp3, mt2);
  91     gmx_mat4_mmul(mvtot, mr[XX], mtemp2);
  92
  93     for (i = 0; (i < natoms); i++)
  94     {
  95         gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
  96         gmx_mat4_transform_point(mvtot, v[i], vrot[i]);
  97     }
  98 }
  99
 100 static void move_x(int natoms, rvec x[], matrix box)
 101 {
 102     int  i, m;
 103     rvec xcm;
 104
 105     clear_rvec(xcm);
 106     for (i = 0; (i < natoms); i++)
 107     {
 108         for (m = 0; (m < DIM); m++)
 109         {
 110             xcm[m] += x[i][m];
 111         }
 112     }
 113     for (m = 0; (m < DIM); m++)
 114     {
 115         xcm[m] = 0.5*box[m][m]-xcm[m]/natoms;
 116     }
 117     for (i = 0; (i < natoms); i++)
 118     {
 119         for (m = 0; (m < DIM); m++)
 120         {
 121             x[i][m] += xcm[m];
 122         }
 123     }
 124 }
 125
 126 int gmx_genconf(int argc, char *argv[])
 127 {
 128     const char     *desc[] = {
 129         "[THISMODULE] multiplies a given coordinate file by simply stacking them",
 130         "on top of each other, like a small child playing with wooden blocks.",
 131         "The program makes a grid of [IT]user-defined[it]",
 132         "proportions ([TT]-nbox[tt]), ",
 133         "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 134         "When option [TT]-rot[tt] is used the program does not check for overlap",
 135         "between molecules on grid points. It is recommended to make the box in",
 136         "the input file at least as big as the coordinates + ",
 137         "van der Waals radius.[PAR]",
 138         "If the optional trajectory file is given, conformations are not",
 139         "generated, but read from this file and translated appropriately to",
 140         "build the grid."
 141
 142     };
 143     const char     *bugs[] = {
 144         "The program should allow for random displacement of lattice points."
 145     };
 146
 147     int             vol;
 148     t_atoms        *atoms;      /* list with all atoms */
 149     char            title[STRLEN];
 150     rvec           *x, *xx, *v; /* coordinates? */
 151     real            t;
 152     vec4           *xrot, *vrot;
 153     int             ePBC;
 154     matrix          box, boxx; /* box length matrix */
 155     rvec            shift;
 156     int             natoms;    /* number of atoms in one molecule  */
 157     int             nres;      /* number of molecules? */
 158     int             i, j, k, l, m, ndx, nrdx, nx, ny, nz;
 159     t_trxstatus    *status;
 160     gmx_bool        bTRX;
 161     output_env_t    oenv;
 162     gmx_rng_t       rng;
 163
 164     t_filenm        fnm[] = {
 165         { efSTX, "-f", "conf", ffREAD  },
 166         { efSTO, "-o", "out",  ffWRITE },
 167         { efTRX, "-trj", NULL,  ffOPTRD }
 168     };
 169 #define NFILE asize(fnm)
 170     static rvec     nrbox    = {1, 1, 1};
 171     static int      seed     = 0;    /* seed for random number generator */
 172     static int      nmolat   = 3;
 173     static int      nblock   = 1;
 174     static gmx_bool bShuffle = FALSE;
 175     static gmx_bool bSort    = FALSE;
 176     static gmx_bool bRandom  = FALSE;           /* False: no random rotations */
 177     static gmx_bool bRenum   = TRUE;            /* renumber residues */
 178     static rvec     dist     = {0, 0, 0};       /* space added between molecules ? */
 179     static rvec     max_rot  = {180, 180, 180}; /* maximum rotation */
 180     t_pargs         pa[]     = {
 181         { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 182         { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 183         { "-seed",   FALSE, etINT,  {&seed},
 184           "Random generator seed, if 0 generated from the time" },
 185         { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 186         { "-shuffle", FALSE, etBOOL, {&bShuffle}, "Random shuffling of molecules" },
 187         { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 188         { "-block",  FALSE, etINT,  {&nblock},
 189           "Divide the box in blocks on this number of cpus" },
 190         { "-nmolat", FALSE, etINT,  {&nmolat},
 191           "Number of atoms per molecule, assumed to start from 0. "
 192           "If you set this wrong, it will screw up your system!" },
 193         { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 194         { "-renumber", FALSE, etBOOL, {&bRenum},  "Renumber residues" }
 195     };
 196
 197     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 198                            asize(desc), desc, asize(bugs), bugs, &oenv))
 199     {
 200         return 0;
 201     }
 202
 203     if (seed == 0)
 204     {
 205         rng = gmx_rng_init(gmx_rng_make_seed());
 206     }
 207     else
 208     {
 209         rng = gmx_rng_init(seed);
 210     }
 211
 212     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 213     nx   = (int)(nrbox[XX]+0.5);
 214     ny   = (int)(nrbox[YY]+0.5);
 215     nz   = (int)(nrbox[ZZ]+0.5);
 216
 217     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 218     {
 219         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 220     }
 221     if ((nmolat <= 0) && bShuffle)
 222     {
 223         gmx_fatal(FARGS, "Can not shuffle if the molecules only have %d atoms",
 224                   nmolat);
 225     }
 226
 227     vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */
 228
 229     get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
 230     snew(atoms, 1);
 231     /* make space for all the atoms */
 232     init_t_atoms(atoms, natoms*vol, FALSE);
 233     snew(x, natoms*vol);           /* get space for coordinates of all atoms */
 234     snew(xrot, natoms);            /* get space for rotation matrix? */
 235     snew(v, natoms*vol);           /* velocities. not really needed? */
 236     snew(vrot, natoms);
 237     /* set atoms->nr to the number in one box *
 238      * to avoid complaints in read_stx_conf   *
 239      */
 240     atoms->nr = natoms;
 241     read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, v, &ePBC, box);
 242
 243     nres = atoms->nres;            /* nr of residues in one element? */
 244
 245     if (bTRX)
 246     {
 247         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 248         {
 249             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 250         }
 251     }
 252     else
 253     {
 254         snew(xx, natoms);
 255         for (i = 0; i < natoms; i++)
 256         {
 257             copy_rvec(x[i], xx[i]);
 258         }
 259     }
 260
 261
 262     for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */
 263     {
 264         shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);
 265
 266         for (j = 0; (j < ny); j++)
 267         {
 268             shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 269
 270             for (i = 0; (i < nx); i++)
 271             {
 272                 shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 273
 274                 ndx  = (i*ny*nz+j*nz+k)*natoms;
 275                 nrdx = (i*ny*nz+j*nz+k)*nres;
 276
 277                 /* Random rotation on input coords */
 278                 if (bRandom)
 279                 {
 280                     rand_rot(natoms, xx, v, xrot, vrot, rng, max_rot);
 281                 }
 282
 283                 for (l = 0; (l < natoms); l++)
 284                 {
 285                     for (m = 0; (m < DIM); m++)
 286                     {
 287                         if (bRandom)
 288                         {
 289                             x[ndx+l][m] = xrot[l][m];
 290                             v[ndx+l][m] = vrot[l][m];
 291                         }
 292                         else
 293                         {
 294                             x[ndx+l][m] = xx[l][m];
 295                             v[ndx+l][m] = v[l][m];
 296                         }
 297                     }
 298                     if (ePBC == epbcSCREW && i % 2 == 1)
 299                     {
 300                         /* Rotate around x axis */
 301                         for (m = YY; m <= ZZ; m++)
 302                         {
 303                             x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 304                             v[ndx+l][m] = -v[ndx+l][m];
 305                         }
 306                     }
 307                     for (m = 0; (m < DIM); m++)
 308                     {
 309                         x[ndx+l][m] += shift[m];
 310                     }
 311                     atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 312                     atoms->atomname[ndx+l]    = atoms->atomname[l];
 313                 }
 314
 315                 for (l = 0; (l < nres); l++)
 316                 {
 317                     atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 318                     if (bRenum)
 319                     {
 320                         atoms->resinfo[nrdx+l].nr += nrdx;
 321                     }
 322                 }
 323                 if (bTRX)
 324                 {
 325                     if (!read_next_x(oenv, status, &t, xx, boxx) &&
 326                         ((i+1)*(j+1)*(k+1) < vol))
 327                     {
 328                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 329                     }
 330                 }
 331             }
 332         }
 333     }
 334     if (bTRX)
 335     {
 336         close_trj(status);
 337     }
 338
 339     /* make box bigger */
 340     for (m = 0; (m < DIM); m++)
 341     {
 342         box[m][m] += dist[m];
 343     }
 344     svmul(nx, box[XX], box[XX]);
 345     svmul(ny, box[YY], box[YY]);
 346     svmul(nz, box[ZZ], box[ZZ]);
 347     if (ePBC == epbcSCREW && nx % 2 == 0)
 348     {
 349         /* With an even number of boxes in x we can forgot about the screw */
 350         ePBC = epbcXYZ;
 351     }
 352
 353     /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 354
 355     atoms->nr   *= vol;
 356     atoms->nres *= vol;
 357
 358     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 359     if (bRenum)
 360     {
 361         for (i = 0; i < atoms->nres; i++)
 362         {
 363             atoms->resinfo[i].nr = i+1;
 364         }
 365     }
 366
 367
 368     if (bShuffle)
 369     {
 370         randwater(0, atoms->nr/nmolat, nmolat, x, v, rng);
 371     }
 372     else if (bSort)
 373     {
 374         sortwater(0, atoms->nr/nmolat, nmolat, x, v);
 375     }
 376     else if (opt2parg_bSet("-block", asize(pa), pa))
 377     {
 378         mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
 379     }
 380     gmx_rng_destroy(rng);
 381
 382     write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);
 383
 384     return 0;
 385 }