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,2015,2017,2018,2019,2020, 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.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/topology/index.h"
50 #include "gromacs/topology/mtop_util.h"
51 #include "gromacs/topology/topology.h"
52 #include "gromacs/utility/arraysize.h"
53 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
58 int gmx_genrestr(int argc
, char* argv
[])
60 const char* desc
[] = {
61 "[THISMODULE] produces an #include file for a topology containing",
62 "a list of atom numbers and three force constants for the",
63 "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-direction based on",
64 "the contents of the [TT]-f[tt] file. A single isotropic force constant may",
65 "be given on the command line instead of three components.[PAR]",
66 "WARNING: Position restraints are interactions within molecules, therefore",
67 "they must be included within the correct [TT][ moleculetype ][tt]",
68 "block in the topology. The atom indices within the",
69 "[TT][ position_restraints ][tt] block must be within the range of the",
70 "atom indices for that molecule type. Since the atom numbers in every",
71 "moleculetype in the topology start at 1 and the numbers in the input file",
72 "for [THISMODULE] number consecutively from 1, [THISMODULE] will only",
73 "produce a useful file for the first molecule. You may wish to",
74 "edit the resulting index file to remove the lines for later atoms,",
75 "or construct a suitable index group to provide",
76 "as input to [THISMODULE].[PAR]",
77 "The [TT]-of[tt] option produces an index file that can be used for",
78 "freezing atoms. In this case, the input file must be a [REF].pdb[ref] file.[PAR]",
79 "With the [TT]-disre[tt] option, half a matrix of distance restraints",
80 "is generated instead of position restraints. With this matrix, that",
81 "one typically would apply to C[GRK]alpha[grk] atoms in a protein, one can",
82 "maintain the overall conformation of a protein without tieing it to",
83 "a specific position (as with position restraints)."
85 static rvec fc
= { 1000.0, 1000.0, 1000.0 };
86 static real freeze_level
= 0.0;
87 static real disre_dist
= 0.1;
88 static real disre_frac
= 0.0;
89 static real disre_up2
= 1.0;
90 static gmx_bool bDisre
= FALSE
;
91 static gmx_bool bConstr
= FALSE
;
92 static real cutoff
= -1.0;
95 { "-fc", FALSE
, etRVEC
, { fc
}, "Force constants (kJ/mol nm^2)" },
100 "If the [TT]-of[tt] option or this one is given an index file will be written containing "
101 "atom numbers of all atoms that have a B-factor less than the level given here" },
106 "Generate a distance restraint matrix for all the atoms in index" },
111 "Distance range around the actual distance for generating distance restraints" },
116 "Fraction of distance to be used as interval rather than a fixed distance. If the "
117 "fraction of the distance that you specify here is less than the distance given in the "
118 "previous option, that one is used instead." },
123 "Distance between upper bound for distance restraints, and the distance at which the "
124 "force becomes constant (see manual)" },
129 "Only generate distance restraints for atoms pairs within cutoff (nm)" },
134 "Generate a constraint matrix rather than distance restraints. Constraints of type 2 "
135 "will be generated that do generate exclusions." }
137 #define npargs asize(pa)
139 gmx_output_env_t
* oenv
;
145 const char * xfn
, *nfn
;
148 rvec dx
, *x
= nullptr, *v
= nullptr;
150 t_filenm fnm
[] = { { efSTX
, "-f", nullptr, ffREAD
},
151 { efNDX
, "-n", nullptr, ffOPTRD
},
152 { efITP
, "-o", "posre", ffWRITE
},
153 { efNDX
, "-of", "freeze", ffOPTWR
} };
154 #define NFILE asize(fnm)
156 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, npargs
, pa
, asize(desc
), desc
, 0, nullptr, &oenv
))
160 output_env_done(oenv
);
162 bFreeze
= opt2bSet("-of", NFILE
, fnm
) || opt2parg_bSet("-freeze", asize(pa
), pa
);
163 bDisre
= bDisre
|| opt2parg_bSet("-disre_dist", npargs
, pa
);
164 xfn
= opt2fn_null("-f", NFILE
, fnm
);
165 nfn
= opt2fn_null("-n", NFILE
, fnm
);
167 if ((nfn
== nullptr) && (xfn
== nullptr))
169 gmx_fatal(FARGS
, "no index file and no structure file supplied");
172 if ((disre_frac
< 0) || (disre_frac
>= 1))
174 gmx_fatal(FARGS
, "disre_frac should be between 0 and 1");
178 gmx_fatal(FARGS
, "disre_dist should be >= 0");
181 const char* title
= "";
182 bool haveTopology
= false;
184 int* indexGroups
= nullptr;
185 char* indexGroupNames
= nullptr;
189 fprintf(stderr
, "\nReading structure file\n");
190 readConfAndTopology(xfn
, &haveTopology
, &mtop
, nullptr, &x
, &v
, box
);
192 atoms
= gmx_mtop_global_atoms(&mtop
);
193 if (atoms
.pdbinfo
== nullptr)
195 snew(atoms
.pdbinfo
, atoms
.nr
);
202 if (!haveTopology
|| !atoms
.pdbinfo
)
204 gmx_fatal(FARGS
, "No B-factors in input file %s, use a pdb file next time.", xfn
);
207 out
= opt2FILE("-of", NFILE
, fnm
, "w");
208 fprintf(out
, "[ freeze ]\n");
209 for (i
= 0; (i
< atoms
.nr
); i
++)
211 if (atoms
.pdbinfo
[i
].bfac
<= freeze_level
)
213 fprintf(out
, "%d\n", i
+ 1);
218 else if ((bDisre
|| bConstr
) && x
)
220 printf("Select group to generate %s matrix from\n", bConstr
? "constraint" : "distance restraint");
221 get_index(&atoms
, nfn
, 1, &igrp
, &indexGroups
, &indexGroupNames
);
223 out
= ftp2FILE(efITP
, NFILE
, fnm
, "w");
226 fprintf(out
, "; constraints for %s of %s\n\n", indexGroupNames
, title
);
227 fprintf(out
, "[ constraints ]\n");
228 fprintf(out
, ";%4s %5s %1s %10s\n", "i", "j", "tp", "dist");
232 fprintf(out
, "; distance restraints for %s of %s\n\n", indexGroupNames
, title
);
233 fprintf(out
, "[ distance_restraints ]\n");
234 fprintf(out
, ";%4s %5s %1s %5s %10s %10s %10s %10s %10s\n", "i", "j", "?", "label",
235 "funct", "lo", "up1", "up2", "weight");
237 for (i
= k
= 0; i
< igrp
; i
++)
239 for (j
= i
+ 1; j
< igrp
; j
++, k
++)
241 rvec_sub(x
[indexGroups
[i
]], x
[indexGroups
[j
]], dx
);
245 fprintf(out
, "%5d %5d %1d %10g\n", indexGroups
[i
] + 1, indexGroups
[j
] + 1, 2, d
);
249 if (cutoff
< 0 || d
< cutoff
)
253 dd
= std::min(disre_dist
, disre_frac
* d
);
259 lo
= std::max(0.0_real
, d
- dd
);
261 fprintf(out
, "%5d %5d %1d %5d %10d %10g %10g %10g %10g\n", indexGroups
[i
] + 1,
262 indexGroups
[j
] + 1, 1, k
, 1, lo
, hi
, hi
+ disre_up2
, 1.0);
271 printf("Select group to position restrain\n");
272 get_index(&atoms
, nfn
, 1, &igrp
, &indexGroups
, &indexGroupNames
);
274 out
= ftp2FILE(efITP
, NFILE
, fnm
, "w");
275 fprintf(out
, "; position restraints for %s of %s\n\n", indexGroupNames
, title
);
276 fprintf(out
, "[ position_restraints ]\n");
277 fprintf(out
, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
278 for (i
= 0; i
< igrp
; i
++)
280 fprintf(out
, "%4d %4d %10g %10g %10g\n", indexGroups
[i
] + 1, 1, fc
[XX
], fc
[YY
], fc
[ZZ
]);
290 sfree(indexGroupNames
);