2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
44 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/utility/fatalerror.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/utility/futil.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/random/random.h"
62 int gmx_nmens(int argc
, char *argv
[])
64 const char *desc
[] = {
65 "[THISMODULE] generates an ensemble around an average structure",
66 "in a subspace that is defined by a set of normal modes (eigenvectors).",
67 "The eigenvectors are assumed to be mass-weighted.",
68 "The position along each eigenvector is randomly taken from a Gaussian",
69 "distribution with variance kT/eigenvalue.[PAR]",
70 "By default the starting eigenvector is set to 7, since the first six",
71 "normal modes are the translational and rotational degrees of freedom."
73 static int nstruct
= 100, first
= 7, last
= -1, seed
= -1;
74 static real temp
= 300.0;
76 { "-temp", FALSE
, etREAL
, {&temp
},
77 "Temperature in Kelvin" },
78 { "-seed", FALSE
, etINT
, {&seed
},
79 "Random seed, -1 generates a seed from time and pid" },
80 { "-num", FALSE
, etINT
, {&nstruct
},
81 "Number of structures to generate" },
82 { "-first", FALSE
, etINT
, {&first
},
83 "First eigenvector to use (-1 is select)" },
84 { "-last", FALSE
, etINT
, {&last
},
85 "Last eigenvector to use (-1 is till the last)" }
94 rvec
*xtop
, *xref
, *xav
, *xout1
, *xout2
;
95 gmx_bool bDMR
, bDMA
, bFit
;
96 int nvec
, *eignr
= NULL
;
99 real
*eigval
, totmass
, *invsqrtm
, t
, disp
;
101 char *grpname
, title
[STRLEN
];
102 const char *indexfile
;
104 int nout
, *iout
, noutvec
, *outvec
;
106 real rfac
, invfr
, rhalf
, jr
;
111 const unsigned long im
= 0xffff;
112 const unsigned long ia
= 1093;
113 const unsigned long ic
= 18257;
117 { efTRN
, "-v", "eigenvec", ffREAD
},
118 { efXVG
, "-e", "eigenval", ffREAD
},
119 { efTPS
, NULL
, NULL
, ffREAD
},
120 { efNDX
, NULL
, NULL
, ffOPTRD
},
121 { efTRO
, "-o", "ensemble", ffWRITE
}
123 #define NFILE asize(fnm)
125 if (!parse_common_args(&argc
, argv
, PCA_BE_NICE
,
126 NFILE
, fnm
, NPA
, pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
131 indexfile
= ftp2fn_null(efNDX
, NFILE
, fnm
);
133 read_eigenvectors(opt2fn("-v", NFILE
, fnm
), &natoms
, &bFit
,
134 &xref
, &bDMR
, &xav
, &bDMA
, &nvec
, &eignr
, &eigvec
, &eigval
);
136 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), title
, &top
, &ePBC
, &xtop
, NULL
, box
, bDMA
);
139 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms
);
140 get_index(atoms
, indexfile
, 1, &i
, &index
, &grpname
);
143 gmx_fatal(FARGS
, "you selected a group with %d elements instead of %d",
148 snew(invsqrtm
, natoms
);
151 for (i
= 0; (i
< natoms
); i
++)
153 invsqrtm
[i
] = gmx_invsqrt(atoms
->atom
[index
[i
]].m
);
158 for (i
= 0; (i
< natoms
); i
++)
170 /* make an index from first to last */
173 for (i
= 0; i
< nout
; i
++)
180 printf("Select eigenvectors for output, end your selection with 0\n");
186 srenew(iout
, nout
+1);
187 if (1 != scanf("%d", &iout
[nout
]))
189 gmx_fatal(FARGS
, "Error reading user input");
193 while (iout
[nout
] >= 0);
197 /* make an index of the eigenvectors which are present */
200 for (i
= 0; i
< nout
; i
++)
203 while ((j
< nvec
) && (eignr
[j
] != iout
[i
]))
207 if ((j
< nvec
) && (eignr
[j
] == iout
[i
]))
210 iout
[noutvec
] = iout
[i
];
215 fprintf(stderr
, "%d eigenvectors selected for output\n", noutvec
);
219 seed
= (int)gmx_rng_make_seed();
220 rng
= gmx_rng_init(seed
);
224 rng
= gmx_rng_init(seed
);
226 fprintf(stderr
, "Using seed %d and a temperature of %g K\n", seed
, temp
);
229 snew(xout2
, atoms
->nr
);
230 out
= open_trx(ftp2fn(efTRO
, NFILE
, fnm
), "w");
231 jran
= (unsigned long)((real
)im
*gmx_rng_uniform_real(rng
));
232 gmx_rng_destroy(rng
);
233 for (s
= 0; s
< nstruct
; s
++)
235 for (i
= 0; i
< natoms
; i
++)
237 copy_rvec(xav
[i
], xout1
[i
]);
239 for (j
= 0; j
< noutvec
; j
++)
242 /* (r-0.5) n times: var_n = n * var_1 = n/12
243 n=4: var_n = 1/3, so multiply with 3 */
245 rfac
= sqrt(3.0 * BOLTZ
*temp
/eigval
[iout
[j
]]);
247 rfac
= rfac
/(real
)im
;
249 jran
= (jran
*ia
+ic
) & im
;
251 jran
= (jran
*ia
+ic
) & im
;
253 jran
= (jran
*ia
+ic
) & im
;
255 jran
= (jran
*ia
+ic
) & im
;
257 disp
= rfac
* jr
- rhalf
;
259 for (i
= 0; i
< natoms
; i
++)
261 for (d
= 0; d
< DIM
; d
++)
263 xout1
[i
][d
] += disp
*eigvec
[v
][i
][d
]*invsqrtm
[i
];
267 for (i
= 0; i
< natoms
; i
++)
269 copy_rvec(xout1
[i
], xout2
[index
[i
]]);
272 write_trx(out
, natoms
, index
, atoms
, 0, t
, box
, xout2
, NULL
, NULL
);
273 fprintf(stderr
, "\rGenerated %d structures", s
+1);
275 fprintf(stderr
, "\n");