src/gromacs/gmxana/gmx_nmens.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,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/pdbio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/gmxana/eigio.h"
  48 #include "gromacs/gmxana/gmx_ana.h"
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/random/threefry.h"
  53 #include "gromacs/random/uniformintdistribution.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/arraysize.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
  63 int gmx_nmens(int argc, char* argv[])
  64 {
  65     const char* desc[] = {
  66         "[THISMODULE] generates an ensemble around an average structure",
  67         "in a subspace that is defined by a set of normal modes (eigenvectors).",
  68         "The eigenvectors are assumed to be mass-weighted.",
  69         "The position along each eigenvector is randomly taken from a Gaussian",
  70         "distribution with variance kT/eigenvalue.[PAR]",
  71         "By default the starting eigenvector is set to 7, since the first six",
  72         "normal modes are the translational and rotational degrees of freedom."
  73     };
  74     static int  nstruct = 100, first = 7, last = -1, seed = 0;
  75     static real temp = 300.0;
  76     t_pargs     pa[] = {
  77         { "-temp", FALSE, etREAL, { &temp }, "Temperature in Kelvin" },
  78         { "-seed", FALSE, etINT, { &seed }, "Random seed (0 means generate)" },
  79         { "-num", FALSE, etINT, { &nstruct }, "Number of structures to generate" },
  80         { "-first", FALSE, etINT, { &first }, "First eigenvector to use (-1 is select)" },
  81         { "-last", FALSE, etINT, { &last }, "Last eigenvector to use (-1 is till the last)" }
  82     };
  83 #define NPA asize(pa)
  84
  85     t_trxstatus*        out;
  86     t_topology          top;
  87     PbcType             pbcType;
  88     t_atoms*            atoms;
  89     rvec *              xtop, *xref, *xav, *xout1, *xout2;
  90     gmx_bool            bDMR, bDMA, bFit;
  91     int                 nvec, *eignr = nullptr;
  92     rvec**              eigvec = nullptr;
  93     matrix              box;
  94     real *              eigval, *invsqrtm, t, disp;
  95     int                 natoms;
  96     char*               grpname;
  97     const char*         indexfile;
  98     int                 i, j, d, s, v;
  99     int                 nout, *iout, noutvec, *outvec;
 100     int*                index;
 101     real                rfac, rhalf, jr;
 102     gmx_output_env_t*   oenv;
 103     int                 jran;
 104     const unsigned long im = 0xffff;
 105     const unsigned long ia = 1093;
 106     const unsigned long ic = 18257;
 107
 108
 109     t_filenm fnm[] = { { efTRN, "-v", "eigenvec", ffREAD },
 110                        { efXVG, "-e", "eigenval", ffREAD },
 111                        { efTPS, nullptr, nullptr, ffREAD },
 112                        { efNDX, nullptr, nullptr, ffOPTRD },
 113                        { efTRO, "-o", "ensemble", ffWRITE } };
 114 #define NFILE asize(fnm)
 115
 116     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
 117     {
 118         return 0;
 119     }
 120
 121     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
 122
 123     read_eigenvectors(opt2fn("-v", NFILE, fnm), &natoms, &bFit, &xref, &bDMR, &xav, &bDMA, &nvec,
 124                       &eignr, &eigvec, &eigval);
 125
 126     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, box, bDMA);
 127     atoms = &top.atoms;
 128
 129     printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
 130     get_index(atoms, indexfile, 1, &i, &index, &grpname);
 131     if (i != natoms)
 132     {
 133         gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
 134     }
 135     printf("\n");
 136
 137     snew(invsqrtm, natoms);
 138     if (bDMA)
 139     {
 140         for (i = 0; (i < natoms); i++)
 141         {
 142             invsqrtm[i] = gmx::invsqrt(atoms->atom[index[i]].m);
 143         }
 144     }
 145     else
 146     {
 147         for (i = 0; (i < natoms); i++)
 148         {
 149             invsqrtm[i] = 1.0;
 150         }
 151     }
 152
 153     if (last == -1)
 154     {
 155         last = natoms * DIM;
 156     }
 157     if (first > -1)
 158     {
 159         /* make an index from first to last */
 160         nout = last - first + 1;
 161         snew(iout, nout);
 162         for (i = 0; i < nout; i++)
 163         {
 164             iout[i] = first - 1 + i;
 165         }
 166     }
 167     else
 168     {
 169         printf("Select eigenvectors for output, end your selection with 0\n");
 170         nout = -1;
 171         iout = nullptr;
 172         do
 173         {
 174             nout++;
 175             srenew(iout, nout + 1);
 176             if (1 != scanf("%d", &iout[nout]))
 177             {
 178                 gmx_fatal(FARGS, "Error reading user input");
 179             }
 180             iout[nout]--;
 181         } while (iout[nout] >= 0);
 182         printf("\n");
 183     }
 184
 185     /* make an index of the eigenvectors which are present */
 186     snew(outvec, nout);
 187     noutvec = 0;
 188     for (i = 0; i < nout; i++)
 189     {
 190         j = 0;
 191         while ((j < nvec) && (eignr[j] != iout[i]))
 192         {
 193             j++;
 194         }
 195         if ((j < nvec) && (eignr[j] == iout[i]))
 196         {
 197             outvec[noutvec] = j;
 198             iout[noutvec]   = iout[i];
 199             noutvec++;
 200         }
 201     }
 202
 203     fprintf(stderr, "%d eigenvectors selected for output\n", noutvec);
 204
 205
 206     if (seed == 0)
 207     {
 208         // Make do with 32 bits for now to avoid changing user input to hex
 209         seed = static_cast<int>(gmx::makeRandomSeed());
 210     }
 211
 212     gmx::DefaultRandomEngine rng(seed);
 213
 214     fprintf(stderr, "Using random seed %d and a temperature of %g K.\n", seed, temp);
 215
 216     gmx::UniformIntDistribution<int> dist(0, im - 1);
 217     jran = dist(rng);
 218
 219     snew(xout1, natoms);
 220     snew(xout2, atoms->nr);
 221     out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
 222
 223     for (s = 0; s < nstruct; s++)
 224     {
 225         for (i = 0; i < natoms; i++)
 226         {
 227             copy_rvec(xav[i], xout1[i]);
 228         }
 229         for (j = 0; j < noutvec; j++)
 230         {
 231             v = outvec[j];
 232             /* (r-0.5) n times:  var_n = n * var_1 = n/12
 233                n=4:  var_n = 1/3, so multiply with 3 */
 234
 235             rfac  = std::sqrt(3.0 * BOLTZ * temp / eigval[iout[j]]);
 236             rhalf = 2.0 * rfac;
 237             rfac  = rfac / im;
 238
 239             jran = (jran * ia + ic) & im;
 240             jr   = jran;
 241             jran = (jran * ia + ic) & im;
 242             jr += jran;
 243             jran = (jran * ia + ic) & im;
 244             jr += jran;
 245             jran = (jran * ia + ic) & im;
 246             jr += jran;
 247             disp = rfac * jr - rhalf;
 248
 249             for (i = 0; i < natoms; i++)
 250             {
 251                 for (d = 0; d < DIM; d++)
 252                 {
 253                     xout1[i][d] += disp * eigvec[v][i][d] * invsqrtm[i];
 254                 }
 255             }
 256         }
 257         for (i = 0; i < natoms; i++)
 258         {
 259             copy_rvec(xout1[i], xout2[index[i]]);
 260         }
 261         t = s + 1;
 262         write_trx(out, natoms, index, atoms, 0, t, box, xout2, nullptr, nullptr);
 263         fprintf(stderr, "\rGenerated %d structures", s + 1);
 264         fflush(stderr);
 265     }
 266     fprintf(stderr, "\n");
 267     close_trx(out);
 268
 269     return 0;
 270 }