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 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.
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.
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.
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.
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.
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.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/topology/mtop_util.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/arraysize.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
59 int gmx_genrestr(int argc
, char* argv
[])
61 const char* desc
[] = {
62 "[THISMODULE] produces an #include file for a topology containing",
63 "a list of atom numbers and three force constants for the",
64 "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-direction based on",
65 "the contents of the [TT]-f[tt] file. A single isotropic force constant may",
66 "be given on the command line instead of three components.[PAR]",
67 "WARNING: Position restraints are interactions within molecules, therefore",
68 "they must be included within the correct [TT][ moleculetype ][tt]",
69 "block in the topology. The atom indices within the",
70 "[TT][ position_restraints ][tt] block must be within the range of the",
71 "atom indices for that molecule type. Since the atom numbers in every",
72 "moleculetype in the topology start at 1 and the numbers in the input file",
73 "for [THISMODULE] number consecutively from 1, [THISMODULE] will only",
74 "produce a useful file for the first molecule. You may wish to",
75 "edit the resulting index file to remove the lines for later atoms,",
76 "or construct a suitable index group to provide",
77 "as input to [THISMODULE].[PAR]",
78 "The [TT]-of[tt] option produces an index file that can be used for",
79 "freezing atoms. In this case, the input file must be a [REF].pdb[ref] file.[PAR]",
80 "With the [TT]-disre[tt] option, half a matrix of distance restraints",
81 "is generated instead of position restraints. With this matrix, that",
82 "one typically would apply to C[GRK]alpha[grk] atoms in a protein, one can",
83 "maintain the overall conformation of a protein without tieing it to",
84 "a specific position (as with position restraints)."
86 static rvec fc
= { 1000.0, 1000.0, 1000.0 };
87 static real freeze_level
= 0.0;
88 static real disre_dist
= 0.1;
89 static real disre_frac
= 0.0;
90 static real disre_up2
= 1.0;
91 static gmx_bool bDisre
= FALSE
;
92 static gmx_bool bConstr
= FALSE
;
93 static real cutoff
= -1.0;
96 { "-fc", FALSE
, etRVEC
, { fc
}, "Force constants (kJ/mol nm^2)" },
101 "If the [TT]-of[tt] option or this one is given an index file will be written containing "
102 "atom numbers of all atoms that have a B-factor less than the level given here" },
107 "Generate a distance restraint matrix for all the atoms in index" },
112 "Distance range around the actual distance for generating distance restraints" },
117 "Fraction of distance to be used as interval rather than a fixed distance. If the "
118 "fraction of the distance that you specify here is less than the distance given in the "
119 "previous option, that one is used instead." },
124 "Distance between upper bound for distance restraints, and the distance at which the "
125 "force becomes constant (see manual)" },
130 "Only generate distance restraints for atoms pairs within cutoff (nm)" },
135 "Generate a constraint matrix rather than distance restraints. Constraints of type 2 "
136 "will be generated that do generate exclusions." }
138 #define npargs asize(pa)
140 gmx_output_env_t
* oenv
;
146 const char * xfn
, *nfn
;
149 rvec dx
, *x
= nullptr, *v
= nullptr;
151 t_filenm fnm
[] = { { efSTX
, "-f", nullptr, ffREAD
},
152 { efNDX
, "-n", nullptr, ffOPTRD
},
153 { efITP
, "-o", "posre", ffWRITE
},
154 { efNDX
, "-of", "freeze", ffOPTWR
} };
155 #define NFILE asize(fnm)
157 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, npargs
, pa
, asize(desc
), desc
, 0, nullptr, &oenv
))
161 output_env_done(oenv
);
163 bFreeze
= opt2bSet("-of", NFILE
, fnm
) || opt2parg_bSet("-freeze", asize(pa
), pa
);
164 bDisre
= bDisre
|| opt2parg_bSet("-disre_dist", npargs
, pa
);
165 xfn
= opt2fn_null("-f", NFILE
, fnm
);
166 nfn
= opt2fn_null("-n", NFILE
, fnm
);
168 if ((nfn
== nullptr) && (xfn
== nullptr))
170 gmx_fatal(FARGS
, "no index file and no structure file supplied");
173 if ((disre_frac
< 0) || (disre_frac
>= 1))
175 gmx_fatal(FARGS
, "disre_frac should be between 0 and 1");
179 gmx_fatal(FARGS
, "disre_dist should be >= 0");
182 const char* title
= "";
183 bool haveTopology
= false;
185 int* indexGroups
= nullptr;
186 char* indexGroupNames
= nullptr;
190 fprintf(stderr
, "\nReading structure file\n");
191 readConfAndTopology(xfn
, &haveTopology
, &mtop
, nullptr, &x
, &v
, box
);
193 atoms
= gmx_mtop_global_atoms(&mtop
);
194 if (atoms
.pdbinfo
== nullptr)
196 snew(atoms
.pdbinfo
, atoms
.nr
);
203 if (!haveTopology
|| !atoms
.pdbinfo
)
205 gmx_fatal(FARGS
, "No B-factors in input file %s, use a pdb file next time.", xfn
);
208 out
= opt2FILE("-of", NFILE
, fnm
, "w");
209 fprintf(out
, "[ freeze ]\n");
210 for (i
= 0; (i
< atoms
.nr
); i
++)
212 if (atoms
.pdbinfo
[i
].bfac
<= freeze_level
)
214 fprintf(out
, "%d\n", i
+ 1);
219 else if ((bDisre
|| bConstr
) && x
)
221 printf("Select group to generate %s matrix from\n", bConstr
? "constraint" : "distance restraint");
222 get_index(&atoms
, nfn
, 1, &igrp
, &indexGroups
, &indexGroupNames
);
224 out
= ftp2FILE(efITP
, NFILE
, fnm
, "w");
227 fprintf(out
, "; constraints for %s of %s\n\n", indexGroupNames
, title
);
228 fprintf(out
, "[ constraints ]\n");
229 fprintf(out
, ";%4s %5s %1s %10s\n", "i", "j", "tp", "dist");
233 fprintf(out
, "; distance restraints for %s of %s\n\n", indexGroupNames
, title
);
234 fprintf(out
, "[ distance_restraints ]\n");
235 fprintf(out
, ";%4s %5s %1s %5s %10s %10s %10s %10s %10s\n", "i", "j", "?", "label",
236 "funct", "lo", "up1", "up2", "weight");
238 for (i
= k
= 0; i
< igrp
; i
++)
240 for (j
= i
+ 1; j
< igrp
; j
++, k
++)
242 rvec_sub(x
[indexGroups
[i
]], x
[indexGroups
[j
]], dx
);
246 fprintf(out
, "%5d %5d %1d %10g\n", indexGroups
[i
] + 1, indexGroups
[j
] + 1, 2, d
);
250 if (cutoff
< 0 || d
< cutoff
)
254 dd
= std::min(disre_dist
, disre_frac
* d
);
260 lo
= std::max(0.0_real
, d
- dd
);
262 fprintf(out
, "%5d %5d %1d %5d %10d %10g %10g %10g %10g\n", indexGroups
[i
] + 1,
263 indexGroups
[j
] + 1, 1, k
, 1, lo
, hi
, hi
+ disre_up2
, 1.0);
272 printf("Select group to position restrain\n");
273 get_index(&atoms
, nfn
, 1, &igrp
, &indexGroups
, &indexGroupNames
);
275 out
= ftp2FILE(efITP
, NFILE
, fnm
, "w");
276 fprintf(out
, "; position restraints for %s of %s\n\n", indexGroupNames
, title
);
277 fprintf(out
, "[ position_restraints ]\n");
278 fprintf(out
, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
279 for (i
= 0; i
< igrp
; i
++)
281 fprintf(out
, "%4d %4d %10g %10g %10g\n", indexGroups
[i
] + 1, 1, fc
[XX
], fc
[YY
], fc
[ZZ
]);
291 sfree(indexGroupNames
);