File: | gromacs/gmxpreprocess/grompp.c |
Location: | line 1954, column 9 |
Description: | Value stored to 'max_spacing' is never read |
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 "grompp.h" |
38 | |
39 | #ifdef HAVE_CONFIG_H1 |
40 | #include <config.h> |
41 | #endif |
42 | |
43 | #include <sys/types.h> |
44 | #include <math.h> |
45 | #include <string.h> |
46 | #include <errno(*__errno_location ()).h> |
47 | #include <limits.h> |
48 | #include <assert.h> |
49 | |
50 | #include "gromacs/utility/smalloc.h" |
51 | #include "macros.h" |
52 | #include "readir.h" |
53 | #include "toputil.h" |
54 | #include "topio.h" |
55 | #include "gromacs/fileio/confio.h" |
56 | #include "readir.h" |
57 | #include "symtab.h" |
58 | #include "names.h" |
59 | #include "grompp-impl.h" |
60 | #include "gromacs/random/random.h" |
61 | #include "gromacs/gmxpreprocess/gen_maxwell_velocities.h" |
62 | #include "gromacs/math/vec.h" |
63 | #include "gromacs/utility/futil.h" |
64 | #include "gromacs/commandline/pargs.h" |
65 | #include "splitter.h" |
66 | #include "gromacs/gmxpreprocess/sortwater.h" |
67 | #include "convparm.h" |
68 | #include "gromacs/utility/fatalerror.h" |
69 | #include "warninp.h" |
70 | #include "index.h" |
71 | #include "gromacs/fileio/gmxfio.h" |
72 | #include "gromacs/fileio/trnio.h" |
73 | #include "gromacs/fileio/tpxio.h" |
74 | #include "gromacs/fileio/trxio.h" |
75 | #include "vsite_parm.h" |
76 | #include "txtdump.h" |
77 | #include "calcgrid.h" |
78 | #include "add_par.h" |
79 | #include "gromacs/fileio/enxio.h" |
80 | #include "perf_est.h" |
81 | #include "compute_io.h" |
82 | #include "gpp_atomtype.h" |
83 | #include "mtop_util.h" |
84 | #include "genborn.h" |
85 | #include "calc_verletbuf.h" |
86 | #include "tomorse.h" |
87 | #include "gromacs/imd/imd.h" |
88 | |
89 | |
90 | static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[]) |
91 | { |
92 | int i, n; |
93 | |
94 | n = 0; |
95 | /* For all the molecule types */ |
96 | for (i = 0; i < nrmols; i++) |
97 | { |
98 | n += mols[i].plist[ifunc].nr; |
99 | mols[i].plist[ifunc].nr = 0; |
100 | } |
101 | return n; |
102 | } |
103 | |
104 | static int check_atom_names(const char *fn1, const char *fn2, |
105 | gmx_mtop_t *mtop, t_atoms *at) |
106 | { |
107 | int mb, m, i, j, nmismatch; |
108 | t_atoms *tat; |
109 | #define MAXMISMATCH20 20 |
110 | |
111 | if (mtop->natoms != at->nr) |
112 | { |
113 | gmx_incons("comparing atom names")_gmx_error("incons", "comparing atom names", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 113); |
114 | } |
115 | |
116 | nmismatch = 0; |
117 | i = 0; |
118 | for (mb = 0; mb < mtop->nmolblock; mb++) |
119 | { |
120 | tat = &mtop->moltype[mtop->molblock[mb].type].atoms; |
121 | for (m = 0; m < mtop->molblock[mb].nmol; m++) |
122 | { |
123 | for (j = 0; j < tat->nr; j++) |
124 | { |
125 | if (strcmp( *(tat->atomname[j]), *(at->atomname[i]) )__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (*(tat->atomname[j])) && __builtin_constant_p (*( at->atomname[i])) && (__s1_len = strlen (*(tat-> atomname[j])), __s2_len = strlen (*(at->atomname[i])), (!( (size_t)(const void *)((*(tat->atomname[j])) + 1) - (size_t )(const void *)(*(tat->atomname[j])) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((*(at->atomname[i ])) + 1) - (size_t)(const void *)(*(at->atomname[i])) == 1 ) || __s2_len >= 4)) ? __builtin_strcmp (*(tat->atomname [j]), *(at->atomname[i])) : (__builtin_constant_p (*(tat-> atomname[j])) && ((size_t)(const void *)((*(tat->atomname [j])) + 1) - (size_t)(const void *)(*(tat->atomname[j])) == 1) && (__s1_len = strlen (*(tat->atomname[j])), __s1_len < 4) ? (__builtin_constant_p (*(at->atomname[i])) && ((size_t)(const void *)((*(at->atomname[i])) + 1) - (size_t )(const void *)(*(at->atomname[i])) == 1) ? __builtin_strcmp (*(tat->atomname[j]), *(at->atomname[i])) : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (*(at->atomname[i])); int __result = (((const unsigned char *) (const char *) (*(tat->atomname[j])))[0] - __s2[0 ]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (*(tat->atomname [j])))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ( *(tat->atomname[j])))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (*(tat->atomname[j])))[3] - __s2[3]); } } __result; }) )) : (__builtin_constant_p (*(at->atomname[i])) && ((size_t)(const void *)((*(at->atomname[i])) + 1) - (size_t )(const void *)(*(at->atomname[i])) == 1) && (__s2_len = strlen (*(at->atomname[i])), __s2_len < 4) ? (__builtin_constant_p (*(tat->atomname[j])) && ((size_t)(const void *)( (*(tat->atomname[j])) + 1) - (size_t)(const void *)(*(tat-> atomname[j])) == 1) ? __builtin_strcmp (*(tat->atomname[j] ), *(at->atomname[i])) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (*(tat-> atomname[j])); int __result = (((const unsigned char *) (const char *) (*(at->atomname[i])))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (*(at->atomname[i])))[1] - __s2[1] ); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (*(at->atomname[ i])))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (*(at ->atomname[i])))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (*(tat->atomname[j]), *(at->atomname[i]))))); }) != 0) |
126 | { |
127 | if (nmismatch < MAXMISMATCH20) |
128 | { |
129 | fprintf(stderrstderr, |
130 | "Warning: atom name %d in %s and %s does not match (%s - %s)\n", |
131 | i+1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i])); |
132 | } |
133 | else if (nmismatch == MAXMISMATCH20) |
134 | { |
135 | fprintf(stderrstderr, "(more than %d non-matching atom names)\n", MAXMISMATCH20); |
136 | } |
137 | nmismatch++; |
138 | } |
139 | i++; |
140 | } |
141 | } |
142 | } |
143 | |
144 | return nmismatch; |
145 | } |
146 | |
147 | static void check_eg_vs_cg(gmx_mtop_t *mtop) |
148 | { |
149 | int astart, mb, m, cg, j, firstj; |
150 | unsigned char firsteg, eg; |
151 | gmx_moltype_t *molt; |
152 | |
153 | /* Go through all the charge groups and make sure all their |
154 | * atoms are in the same energy group. |
155 | */ |
156 | |
157 | astart = 0; |
158 | for (mb = 0; mb < mtop->nmolblock; mb++) |
159 | { |
160 | molt = &mtop->moltype[mtop->molblock[mb].type]; |
161 | for (m = 0; m < mtop->molblock[mb].nmol; m++) |
162 | { |
163 | for (cg = 0; cg < molt->cgs.nr; cg++) |
164 | { |
165 | /* Get the energy group of the first atom in this charge group */ |
166 | firstj = astart + molt->cgs.index[cg]; |
167 | firsteg = ggrpnr(&mtop->groups, egcENER, firstj)((&mtop->groups)->grpnr[egcENER] ? (&mtop->groups )->grpnr[egcENER][firstj] : 0); |
168 | for (j = molt->cgs.index[cg]+1; j < molt->cgs.index[cg+1]; j++) |
169 | { |
170 | eg = ggrpnr(&mtop->groups, egcENER, astart+j)((&mtop->groups)->grpnr[egcENER] ? (&mtop->groups )->grpnr[egcENER][astart+j] : 0); |
171 | if (eg != firsteg) |
172 | { |
173 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 173, "atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups", |
174 | firstj+1, astart+j+1, cg+1, *molt->name); |
175 | } |
176 | } |
177 | } |
178 | astart += molt->atoms.nr; |
179 | } |
180 | } |
181 | } |
182 | |
183 | static void check_cg_sizes(const char *topfn, t_block *cgs, warninp_t wi) |
184 | { |
185 | int maxsize, cg; |
186 | char warn_buf[STRLEN4096]; |
187 | |
188 | maxsize = 0; |
189 | for (cg = 0; cg < cgs->nr; cg++) |
190 | { |
191 | maxsize = max(maxsize, cgs->index[cg+1]-cgs->index[cg])(((maxsize) > (cgs->index[cg+1]-cgs->index[cg])) ? ( maxsize) : (cgs->index[cg+1]-cgs->index[cg]) ); |
192 | } |
193 | |
194 | if (maxsize > MAX_CHARGEGROUP_SIZE32) |
195 | { |
196 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 196, "The largest charge group contains %d atoms. The maximum is %d.", maxsize, MAX_CHARGEGROUP_SIZE32); |
197 | } |
198 | else if (maxsize > 10) |
199 | { |
200 | set_warning_line(wi, topfn, -1); |
201 | sprintf(warn_buf, |
202 | "The largest charge group contains %d atoms.\n" |
203 | "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n" |
204 | "For efficiency and accuracy, charge group should consist of a few atoms.\n" |
205 | "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.", |
206 | maxsize); |
207 | warning_note(wi, warn_buf); |
208 | } |
209 | } |
210 | |
211 | static void check_bonds_timestep(gmx_mtop_t *mtop, double dt, warninp_t wi) |
212 | { |
213 | /* This check is not intended to ensure accurate integration, |
214 | * rather it is to signal mistakes in the mdp settings. |
215 | * A common mistake is to forget to turn on constraints |
216 | * for MD after energy minimization with flexible bonds. |
217 | * This check can also detect too large time steps for flexible water |
218 | * models, but such errors will often be masked by the constraints |
219 | * mdp options, which turns flexible water into water with bond constraints, |
220 | * but without an angle constraint. Unfortunately such incorrect use |
221 | * of water models can not easily be detected without checking |
222 | * for specific model names. |
223 | * |
224 | * The stability limit of leap-frog or velocity verlet is 4.44 steps |
225 | * per oscillational period. |
226 | * But accurate bonds distributions are lost far before that limit. |
227 | * To allow relatively common schemes (although not common with Gromacs) |
228 | * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints |
229 | * we set the note limit to 10. |
230 | */ |
231 | int min_steps_warn = 5; |
232 | int min_steps_note = 10; |
233 | t_iparams *ip; |
234 | int molt; |
235 | gmx_moltype_t *moltype, *w_moltype; |
236 | t_atom *atom; |
237 | t_ilist *ilist, *ilb, *ilc, *ils; |
238 | int ftype; |
239 | int i, a1, a2, w_a1, w_a2, j; |
240 | real twopi2, limit2, fc, re, m1, m2, period2, w_period2; |
241 | gmx_bool bFound, bWater, bWarn; |
242 | char warn_buf[STRLEN4096]; |
243 | |
244 | ip = mtop->ffparams.iparams; |
245 | |
246 | twopi2 = sqr(2*M_PI3.14159265358979323846); |
247 | |
248 | limit2 = sqr(min_steps_note*dt); |
249 | |
250 | w_a1 = w_a2 = -1; |
251 | w_period2 = -1.0; |
252 | |
253 | w_moltype = NULL((void*)0); |
254 | for (molt = 0; molt < mtop->nmoltype; molt++) |
255 | { |
256 | moltype = &mtop->moltype[molt]; |
257 | atom = moltype->atoms.atom; |
258 | ilist = moltype->ilist; |
259 | ilc = &ilist[F_CONSTR]; |
260 | ils = &ilist[F_SETTLE]; |
261 | for (ftype = 0; ftype < F_NRE; ftype++) |
262 | { |
263 | if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC)) |
264 | { |
265 | continue; |
266 | } |
267 | |
268 | ilb = &ilist[ftype]; |
269 | for (i = 0; i < ilb->nr; i += 3) |
270 | { |
271 | fc = ip[ilb->iatoms[i]].harmonic.krA; |
272 | re = ip[ilb->iatoms[i]].harmonic.rA; |
273 | if (ftype == F_G96BONDS) |
274 | { |
275 | /* Convert squared sqaure fc to harmonic fc */ |
276 | fc = 2*fc*re; |
277 | } |
278 | a1 = ilb->iatoms[i+1]; |
279 | a2 = ilb->iatoms[i+2]; |
280 | m1 = atom[a1].m; |
281 | m2 = atom[a2].m; |
282 | if (fc > 0 && m1 > 0 && m2 > 0) |
283 | { |
284 | period2 = twopi2*m1*m2/((m1 + m2)*fc); |
285 | } |
286 | else |
287 | { |
288 | period2 = GMX_FLOAT_MAX3.40282346E+38; |
289 | } |
290 | if (debug) |
291 | { |
292 | fprintf(debug, "fc %g m1 %g m2 %g period %g\n", |
293 | fc, m1, m2, sqrt(period2)); |
294 | } |
295 | if (period2 < limit2) |
296 | { |
297 | bFound = FALSE0; |
298 | for (j = 0; j < ilc->nr; j += 3) |
299 | { |
300 | if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) || |
301 | (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1)) |
302 | { |
303 | bFound = TRUE1; |
304 | } |
305 | } |
306 | for (j = 0; j < ils->nr; j += 4) |
307 | { |
308 | if ((a1 == ils->iatoms[j+1] || a1 == ils->iatoms[j+2] || a1 == ils->iatoms[j+3]) && |
309 | (a2 == ils->iatoms[j+1] || a2 == ils->iatoms[j+2] || a2 == ils->iatoms[j+3])) |
310 | { |
311 | bFound = TRUE1; |
312 | } |
313 | } |
314 | if (!bFound && |
315 | (w_moltype == NULL((void*)0) || period2 < w_period2)) |
316 | { |
317 | w_moltype = moltype; |
318 | w_a1 = a1; |
319 | w_a2 = a2; |
320 | w_period2 = period2; |
321 | } |
322 | } |
323 | } |
324 | } |
325 | } |
326 | |
327 | if (w_moltype != NULL((void*)0)) |
328 | { |
329 | bWarn = (w_period2 < sqr(min_steps_warn*dt)); |
330 | /* A check that would recognize most water models */ |
331 | bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' && |
332 | w_moltype->atoms.nr <= 5); |
333 | sprintf(warn_buf, "The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n" |
334 | "%s", |
335 | *w_moltype->name, |
336 | w_a1+1, *w_moltype->atoms.atomname[w_a1], |
337 | w_a2+1, *w_moltype->atoms.atomname[w_a2], |
338 | sqrt(w_period2), bWarn ? min_steps_warn : min_steps_note, dt, |
339 | bWater ? |
340 | "Maybe you asked for fexible water." : |
341 | "Maybe you forgot to change the constraints mdp option."); |
342 | if (bWarn) |
343 | { |
344 | warning(wi, warn_buf); |
345 | } |
346 | else |
347 | { |
348 | warning_note(wi, warn_buf); |
349 | } |
350 | } |
351 | } |
352 | |
353 | static void check_vel(gmx_mtop_t *mtop, rvec v[]) |
354 | { |
355 | gmx_mtop_atomloop_all_t aloop; |
356 | t_atom *atom; |
357 | int a; |
358 | |
359 | aloop = gmx_mtop_atomloop_all_init(mtop); |
360 | while (gmx_mtop_atomloop_all_next(aloop, &a, &atom)) |
361 | { |
362 | if (atom->ptype == eptShell || |
363 | atom->ptype == eptBond || |
364 | atom->ptype == eptVSite) |
365 | { |
366 | clear_rvec(v[a]); |
367 | } |
368 | } |
369 | } |
370 | |
371 | static gmx_bool nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype) |
372 | { |
373 | int nint, mb; |
374 | |
375 | nint = 0; |
376 | for (mb = 0; mb < mtop->nmolblock; mb++) |
377 | { |
378 | nint += mtop->molblock[mb].nmol*mi[mtop->molblock[mb].type].plist[ftype].nr; |
379 | } |
380 | |
381 | return nint; |
382 | } |
383 | |
384 | /* This routine reorders the molecule type array |
385 | * in the order of use in the molblocks, |
386 | * unused molecule types are deleted. |
387 | */ |
388 | static void renumber_moltypes(gmx_mtop_t *sys, |
389 | int *nmolinfo, t_molinfo **molinfo) |
390 | { |
391 | int *order, norder, i; |
392 | int mb, mi; |
393 | t_molinfo *minew; |
394 | |
395 | snew(order, *nmolinfo)(order) = save_calloc("order", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 395, (*nmolinfo), sizeof(*(order))); |
396 | norder = 0; |
397 | for (mb = 0; mb < sys->nmolblock; mb++) |
398 | { |
399 | for (i = 0; i < norder; i++) |
400 | { |
401 | if (order[i] == sys->molblock[mb].type) |
402 | { |
403 | break; |
404 | } |
405 | } |
406 | if (i == norder) |
407 | { |
408 | /* This type did not occur yet, add it */ |
409 | order[norder] = sys->molblock[mb].type; |
410 | /* Renumber the moltype in the topology */ |
411 | norder++; |
412 | } |
413 | sys->molblock[mb].type = i; |
414 | } |
415 | |
416 | /* We still need to reorder the molinfo structs */ |
417 | snew(minew, norder)(minew) = save_calloc("minew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 417, (norder), sizeof(*(minew))); |
418 | for (mi = 0; mi < *nmolinfo; mi++) |
419 | { |
420 | for (i = 0; i < norder; i++) |
421 | { |
422 | if (order[i] == mi) |
423 | { |
424 | break; |
425 | } |
426 | } |
427 | if (i == norder) |
428 | { |
429 | done_mi(&(*molinfo)[mi]); |
430 | } |
431 | else |
432 | { |
433 | minew[i] = (*molinfo)[mi]; |
434 | } |
435 | } |
436 | sfree(*molinfo)save_free("*molinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 436, (*molinfo)); |
437 | |
438 | *nmolinfo = norder; |
439 | *molinfo = minew; |
440 | } |
441 | |
442 | static void molinfo2mtop(int nmi, t_molinfo *mi, gmx_mtop_t *mtop) |
443 | { |
444 | int m; |
445 | gmx_moltype_t *molt; |
446 | |
447 | mtop->nmoltype = nmi; |
448 | snew(mtop->moltype, nmi)(mtop->moltype) = save_calloc("mtop->moltype", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 448, (nmi), sizeof(*(mtop->moltype))); |
449 | for (m = 0; m < nmi; m++) |
450 | { |
451 | molt = &mtop->moltype[m]; |
452 | molt->name = mi[m].name; |
453 | molt->atoms = mi[m].atoms; |
454 | /* ilists are copied later */ |
455 | molt->cgs = mi[m].cgs; |
456 | molt->excls = mi[m].excls; |
457 | } |
458 | } |
459 | |
460 | static void |
461 | new_status(const char *topfile, const char *topppfile, const char *confin, |
462 | t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero, |
463 | gmx_bool bGenVel, gmx_bool bVerbose, t_state *state, |
464 | gpp_atomtype_t atype, gmx_mtop_t *sys, |
465 | int *nmi, t_molinfo **mi, t_params plist[], |
466 | int *comb, double *reppow, real *fudgeQQ, |
467 | gmx_bool bMorse, |
468 | warninp_t wi) |
469 | { |
470 | t_molinfo *molinfo = NULL((void*)0); |
471 | int nmolblock; |
472 | gmx_molblock_t *molblock, *molbs; |
473 | t_atoms *confat; |
474 | int mb, i, nrmols, nmismatch; |
475 | char buf[STRLEN4096]; |
476 | gmx_bool bGB = FALSE0; |
477 | char warn_buf[STRLEN4096]; |
478 | |
479 | init_mtop(sys); |
480 | |
481 | /* Set gmx_boolean for GB */ |
482 | if (ir->implicit_solvent) |
483 | { |
484 | bGB = TRUE1; |
485 | } |
486 | |
487 | /* TOPOLOGY processing */ |
488 | sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab), |
489 | plist, comb, reppow, fudgeQQ, |
490 | atype, &nrmols, &molinfo, ir, |
491 | &nmolblock, &molblock, bGB, |
492 | wi); |
493 | |
494 | sys->nmolblock = 0; |
495 | snew(sys->molblock, nmolblock)(sys->molblock) = save_calloc("sys->molblock", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 495, (nmolblock), sizeof(*(sys->molblock))); |
496 | |
497 | sys->natoms = 0; |
498 | for (mb = 0; mb < nmolblock; mb++) |
499 | { |
500 | if (sys->nmolblock > 0 && |
501 | molblock[mb].type == sys->molblock[sys->nmolblock-1].type) |
502 | { |
503 | /* Merge consecutive blocks with the same molecule type */ |
504 | sys->molblock[sys->nmolblock-1].nmol += molblock[mb].nmol; |
505 | sys->natoms += molblock[mb].nmol*sys->molblock[sys->nmolblock-1].natoms_mol; |
506 | } |
507 | else if (molblock[mb].nmol > 0) |
508 | { |
509 | /* Add a new molblock to the topology */ |
510 | molbs = &sys->molblock[sys->nmolblock]; |
511 | *molbs = molblock[mb]; |
512 | molbs->natoms_mol = molinfo[molbs->type].atoms.nr; |
513 | molbs->nposres_xA = 0; |
514 | molbs->nposres_xB = 0; |
515 | sys->natoms += molbs->nmol*molbs->natoms_mol; |
516 | sys->nmolblock++; |
517 | } |
518 | } |
519 | if (sys->nmolblock == 0) |
520 | { |
521 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 521, "No molecules were defined in the system"); |
522 | } |
523 | |
524 | renumber_moltypes(sys, &nrmols, &molinfo); |
525 | |
526 | if (bMorse) |
527 | { |
528 | convert_harmonics(nrmols, molinfo, atype); |
529 | } |
530 | |
531 | if (ir->eDisre == edrNone) |
532 | { |
533 | i = rm_interactions(F_DISRES, nrmols, molinfo); |
534 | if (i > 0) |
535 | { |
536 | set_warning_line(wi, "unknown", -1); |
537 | sprintf(warn_buf, "disre = no, removed %d distance restraints", i); |
538 | warning_note(wi, warn_buf); |
539 | } |
540 | } |
541 | if (opts->bOrire == FALSE0) |
542 | { |
543 | i = rm_interactions(F_ORIRES, nrmols, molinfo); |
544 | if (i > 0) |
545 | { |
546 | set_warning_line(wi, "unknown", -1); |
547 | sprintf(warn_buf, "orire = no, removed %d orientation restraints", i); |
548 | warning_note(wi, warn_buf); |
549 | } |
550 | } |
551 | |
552 | /* Copy structures from msys to sys */ |
553 | molinfo2mtop(nrmols, molinfo, sys); |
554 | |
555 | gmx_mtop_finalize(sys); |
556 | |
557 | /* COORDINATE file processing */ |
558 | if (bVerbose) |
559 | { |
560 | fprintf(stderrstderr, "processing coordinates...\n"); |
561 | } |
562 | |
563 | get_stx_coordnum(confin, &state->natoms); |
564 | if (state->natoms != sys->natoms) |
565 | { |
566 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 566, "number of coordinates in coordinate file (%s, %d)\n" |
567 | " does not match topology (%s, %d)", |
568 | confin, state->natoms, topfile, sys->natoms); |
569 | } |
570 | else |
571 | { |
572 | /* make space for coordinates and velocities */ |
573 | char title[STRLEN4096]; |
574 | snew(confat, 1)(confat) = save_calloc("confat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 574, (1), sizeof(*(confat))); |
575 | init_t_atoms(confat, state->natoms, FALSE0); |
576 | init_state(state, state->natoms, 0, 0, 0, 0); |
577 | read_stx_conf(confin, title, confat, state->x, state->v, NULL((void*)0), state->box); |
578 | /* This call fixes the box shape for runs with pressure scaling */ |
579 | set_box_rel(ir, state); |
580 | |
581 | nmismatch = check_atom_names(topfile, confin, sys, confat); |
582 | free_t_atoms(confat, TRUE1); |
583 | sfree(confat)save_free("confat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 583, (confat)); |
584 | |
585 | if (nmismatch) |
586 | { |
587 | sprintf(buf, "%d non-matching atom name%s\n" |
588 | "atom names from %s will be used\n" |
589 | "atom names from %s will be ignored\n", |
590 | nmismatch, (nmismatch == 1) ? "" : "s", topfile, confin); |
591 | warning(wi, buf); |
592 | } |
593 | if (bVerbose) |
594 | { |
595 | fprintf(stderrstderr, "double-checking input for internal consistency...\n"); |
596 | } |
597 | double_check(ir, state->box, nint_ftype(sys, molinfo, F_CONSTR), wi); |
598 | } |
599 | |
600 | if (bGenVel) |
601 | { |
602 | real *mass; |
603 | gmx_mtop_atomloop_all_t aloop; |
604 | t_atom *atom; |
605 | unsigned int useed; |
606 | |
607 | snew(mass, state->natoms)(mass) = save_calloc("mass", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 607, (state->natoms), sizeof(*(mass))); |
608 | aloop = gmx_mtop_atomloop_all_init(sys); |
609 | while (gmx_mtop_atomloop_all_next(aloop, &i, &atom)) |
610 | { |
611 | mass[i] = atom->m; |
612 | } |
613 | |
614 | useed = opts->seed; |
615 | if (opts->seed == -1) |
616 | { |
617 | useed = (int)gmx_rng_make_seed(); |
618 | fprintf(stderrstderr, "Setting gen_seed to %u\n", useed); |
619 | } |
620 | maxwell_speed(opts->tempi, useed, sys, state->v); |
621 | |
622 | stop_cm(stdoutstdout, state->natoms, mass, state->x, state->v); |
623 | sfree(mass)save_free("mass", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 623, (mass)); |
624 | } |
625 | |
626 | *nmi = nrmols; |
627 | *mi = molinfo; |
628 | } |
629 | |
630 | static void copy_state(const char *slog, t_trxframe *fr, |
631 | gmx_bool bReadVel, t_state *state, |
632 | double *use_time) |
633 | { |
634 | int i; |
635 | |
636 | if (fr->not_ok & FRAME_NOT_OK((1<<0) | (1<<1))) |
637 | { |
638 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 638, "Can not start from an incomplete frame"); |
639 | } |
640 | if (!fr->bX) |
641 | { |
642 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 642, "Did not find a frame with coordinates in file %s", |
643 | slog); |
644 | } |
645 | |
646 | for (i = 0; i < state->natoms; i++) |
647 | { |
648 | copy_rvec(fr->x[i], state->x[i]); |
649 | } |
650 | if (bReadVel) |
651 | { |
652 | if (!fr->bV) |
653 | { |
654 | gmx_incons("Trajecory frame unexpectedly does not contain velocities")_gmx_error("incons", "Trajecory frame unexpectedly does not contain velocities" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 654); |
655 | } |
656 | for (i = 0; i < state->natoms; i++) |
657 | { |
658 | copy_rvec(fr->v[i], state->v[i]); |
659 | } |
660 | } |
661 | if (fr->bBox) |
662 | { |
663 | copy_mat(fr->box, state->box); |
664 | } |
665 | |
666 | *use_time = fr->time; |
667 | } |
668 | |
669 | static void cont_status(const char *slog, const char *ener, |
670 | gmx_bool bNeedVel, gmx_bool bGenVel, real fr_time, |
671 | t_inputrec *ir, t_state *state, |
672 | gmx_mtop_t *sys, |
673 | const output_env_t oenv) |
674 | /* If fr_time == -1 read the last frame available which is complete */ |
675 | { |
676 | gmx_bool bReadVel; |
677 | t_trxframe fr; |
678 | t_trxstatus *fp; |
679 | int i; |
680 | double use_time; |
681 | |
682 | bReadVel = (bNeedVel && !bGenVel); |
683 | |
684 | fprintf(stderrstderr, |
685 | "Reading Coordinates%s and Box size from old trajectory\n", |
686 | bReadVel ? ", Velocities" : ""); |
687 | if (fr_time == -1) |
688 | { |
689 | fprintf(stderrstderr, "Will read whole trajectory\n"); |
690 | } |
691 | else |
692 | { |
693 | fprintf(stderrstderr, "Will read till time %g\n", fr_time); |
694 | } |
695 | if (!bReadVel) |
696 | { |
697 | if (bGenVel) |
698 | { |
699 | fprintf(stderrstderr, "Velocities generated: " |
700 | "ignoring velocities in input trajectory\n"); |
701 | } |
702 | read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X(1<<1)); |
703 | } |
704 | else |
705 | { |
706 | read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X(1<<1) | TRX_NEED_V(1<<3)); |
707 | |
708 | if (!fr.bV) |
709 | { |
710 | fprintf(stderrstderr, |
711 | "\n" |
712 | "WARNING: Did not find a frame with velocities in file %s,\n" |
713 | " all velocities will be set to zero!\n\n", slog); |
714 | for (i = 0; i < sys->natoms; i++) |
715 | { |
716 | clear_rvec(state->v[i]); |
717 | } |
718 | close_trj(fp); |
719 | /* Search for a frame without velocities */ |
720 | bReadVel = FALSE0; |
721 | read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X(1<<1)); |
722 | } |
723 | } |
724 | |
725 | state->natoms = fr.natoms; |
726 | |
727 | if (sys->natoms != state->natoms) |
728 | { |
729 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 729, "Number of atoms in Topology " |
730 | "is not the same as in Trajectory"); |
731 | } |
732 | copy_state(slog, &fr, bReadVel, state, &use_time); |
733 | |
734 | /* Find the appropriate frame */ |
735 | while ((fr_time == -1 || fr.time < fr_time) && |
736 | read_next_frame(oenv, fp, &fr)) |
737 | { |
738 | copy_state(slog, &fr, bReadVel, state, &use_time); |
739 | } |
740 | |
741 | close_trj(fp); |
742 | |
743 | /* Set the relative box lengths for preserving the box shape. |
744 | * Note that this call can lead to differences in the last bit |
745 | * with respect to using gmx convert-tpr to create a [TT].tpx[tt] file. |
746 | */ |
747 | set_box_rel(ir, state); |
748 | |
749 | fprintf(stderrstderr, "Using frame at t = %g ps\n", use_time); |
750 | fprintf(stderrstderr, "Starting time for run is %g ps\n", ir->init_t); |
751 | |
752 | if ((ir->epc != epcNO || ir->etc == etcNOSEHOOVER) && ener) |
753 | { |
754 | get_enx_state(ener, use_time, &sys->groups, ir, state); |
755 | preserve_box_shape(ir, state->box_rel, state->boxv); |
756 | } |
757 | } |
758 | |
759 | static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB, |
760 | char *fn, |
761 | int rc_scaling, int ePBC, |
762 | rvec com, |
763 | warninp_t wi) |
764 | { |
765 | gmx_bool bFirst = TRUE1, *hadAtom; |
766 | rvec *x, *v, *xp; |
767 | dvec sum; |
768 | double totmass; |
769 | t_atoms dumat; |
770 | matrix box, invbox; |
771 | int natoms, npbcdim = 0; |
772 | char warn_buf[STRLEN4096], title[STRLEN4096]; |
773 | int a, i, ai, j, k, mb, nat_molb; |
774 | gmx_molblock_t *molb; |
775 | t_params *pr, *prfb; |
776 | t_atom *atom; |
777 | |
778 | get_stx_coordnum(fn, &natoms); |
779 | if (natoms != mtop->natoms) |
780 | { |
781 | sprintf(warn_buf, "The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.", fn, natoms, mtop->natoms, min(mtop->natoms, natoms)(((mtop->natoms) < (natoms)) ? (mtop->natoms) : (natoms ) ), fn); |
782 | warning(wi, warn_buf); |
783 | } |
784 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 784, (natoms), sizeof(*(x))); |
785 | snew(v, natoms)(v) = save_calloc("v", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 785, (natoms), sizeof(*(v))); |
786 | init_t_atoms(&dumat, natoms, FALSE0); |
787 | read_stx_conf(fn, title, &dumat, x, v, NULL((void*)0), box); |
788 | |
789 | npbcdim = ePBC2npbcdim(ePBC); |
790 | clear_rvec(com); |
791 | if (rc_scaling != erscNO) |
792 | { |
793 | copy_mat(box, invbox); |
794 | for (j = npbcdim; j < DIM3; j++) |
795 | { |
796 | clear_rvec(invbox[j]); |
797 | invbox[j][j] = 1; |
798 | } |
799 | m_inv_ur0(invbox, invbox); |
800 | } |
801 | |
802 | /* Copy the reference coordinates to mtop */ |
803 | clear_dvec(sum); |
804 | totmass = 0; |
805 | a = 0; |
806 | snew(hadAtom, natoms)(hadAtom) = save_calloc("hadAtom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 806, (natoms), sizeof(*(hadAtom))); |
807 | for (mb = 0; mb < mtop->nmolblock; mb++) |
808 | { |
809 | molb = &mtop->molblock[mb]; |
810 | nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr; |
811 | pr = &(molinfo[molb->type].plist[F_POSRES]); |
812 | prfb = &(molinfo[molb->type].plist[F_FBPOSRES]); |
813 | if (pr->nr > 0 || prfb->nr > 0) |
814 | { |
815 | atom = mtop->moltype[molb->type].atoms.atom; |
816 | for (i = 0; (i < pr->nr); i++) |
817 | { |
818 | ai = pr->param[i].AIa[0]; |
819 | if (ai >= natoms) |
820 | { |
821 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 821, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n", |
822 | ai+1, *molinfo[molb->type].name, fn, natoms); |
823 | } |
824 | hadAtom[ai] = TRUE1; |
825 | if (rc_scaling == erscCOM) |
826 | { |
827 | /* Determine the center of mass of the posres reference coordinates */ |
828 | for (j = 0; j < npbcdim; j++) |
829 | { |
830 | sum[j] += atom[ai].m*x[a+ai][j]; |
831 | } |
832 | totmass += atom[ai].m; |
833 | } |
834 | } |
835 | /* Same for flat-bottomed posres, but do not count an atom twice for COM */ |
836 | for (i = 0; (i < prfb->nr); i++) |
837 | { |
838 | ai = prfb->param[i].AIa[0]; |
839 | if (ai >= natoms) |
840 | { |
841 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 841, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n", |
842 | ai+1, *molinfo[molb->type].name, fn, natoms); |
843 | } |
844 | if (rc_scaling == erscCOM && hadAtom[ai] == FALSE0) |
845 | { |
846 | /* Determine the center of mass of the posres reference coordinates */ |
847 | for (j = 0; j < npbcdim; j++) |
848 | { |
849 | sum[j] += atom[ai].m*x[a+ai][j]; |
850 | } |
851 | totmass += atom[ai].m; |
852 | } |
853 | } |
854 | if (!bTopB) |
855 | { |
856 | molb->nposres_xA = nat_molb; |
857 | snew(molb->posres_xA, molb->nposres_xA)(molb->posres_xA) = save_calloc("molb->posres_xA", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 857, (molb->nposres_xA), sizeof(*(molb->posres_xA))); |
858 | for (i = 0; i < nat_molb; i++) |
859 | { |
860 | copy_rvec(x[a+i], molb->posres_xA[i]); |
861 | } |
862 | } |
863 | else |
864 | { |
865 | molb->nposres_xB = nat_molb; |
866 | snew(molb->posres_xB, molb->nposres_xB)(molb->posres_xB) = save_calloc("molb->posres_xB", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 866, (molb->nposres_xB), sizeof(*(molb->posres_xB))); |
867 | for (i = 0; i < nat_molb; i++) |
868 | { |
869 | copy_rvec(x[a+i], molb->posres_xB[i]); |
870 | } |
871 | } |
872 | } |
873 | a += nat_molb; |
874 | } |
875 | if (rc_scaling == erscCOM) |
876 | { |
877 | if (totmass == 0) |
878 | { |
879 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 879, "The total mass of the position restraint atoms is 0"); |
880 | } |
881 | for (j = 0; j < npbcdim; j++) |
882 | { |
883 | com[j] = sum[j]/totmass; |
884 | } |
885 | fprintf(stderrstderr, "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n", com[XX0], com[YY1], com[ZZ2]); |
886 | } |
887 | |
888 | if (rc_scaling != erscNO) |
889 | { |
890 | assert(npbcdim <= DIM)((void) (0)); |
891 | |
892 | for (mb = 0; mb < mtop->nmolblock; mb++) |
893 | { |
894 | molb = &mtop->molblock[mb]; |
895 | nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr; |
896 | if (molb->nposres_xA > 0 || molb->nposres_xB > 0) |
897 | { |
898 | xp = (!bTopB ? molb->posres_xA : molb->posres_xB); |
899 | for (i = 0; i < nat_molb; i++) |
900 | { |
901 | for (j = 0; j < npbcdim; j++) |
902 | { |
903 | if (rc_scaling == erscALL) |
904 | { |
905 | /* Convert from Cartesian to crystal coordinates */ |
906 | xp[i][j] *= invbox[j][j]; |
907 | for (k = j+1; k < npbcdim; k++) |
908 | { |
909 | xp[i][j] += invbox[k][j]*xp[i][k]; |
910 | } |
911 | } |
912 | else if (rc_scaling == erscCOM) |
913 | { |
914 | /* Subtract the center of mass */ |
915 | xp[i][j] -= com[j]; |
916 | } |
917 | } |
918 | } |
919 | } |
920 | } |
921 | |
922 | if (rc_scaling == erscCOM) |
923 | { |
924 | /* Convert the COM from Cartesian to crystal coordinates */ |
925 | for (j = 0; j < npbcdim; j++) |
926 | { |
927 | com[j] *= invbox[j][j]; |
928 | for (k = j+1; k < npbcdim; k++) |
929 | { |
930 | com[j] += invbox[k][j]*com[k]; |
931 | } |
932 | } |
933 | } |
934 | } |
935 | |
936 | free_t_atoms(&dumat, TRUE1); |
937 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 937, (x)); |
938 | sfree(v)save_free("v", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 938, (v)); |
939 | sfree(hadAtom)save_free("hadAtom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 939, (hadAtom)); |
940 | } |
941 | |
942 | static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi, |
943 | char *fnA, char *fnB, |
944 | int rc_scaling, int ePBC, |
945 | rvec com, rvec comB, |
946 | warninp_t wi) |
947 | { |
948 | int i, j; |
949 | |
950 | read_posres (mtop, mi, FALSE0, fnA, rc_scaling, ePBC, com, wi); |
951 | if (strcmp(fnA, fnB)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (fnA) && __builtin_constant_p (fnB) && (__s1_len = strlen (fnA), __s2_len = strlen (fnB), (!((size_t)(const void *)((fnA) + 1) - (size_t)(const void *)(fnA) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((fnB) + 1) - ( size_t)(const void *)(fnB) == 1) || __s2_len >= 4)) ? __builtin_strcmp (fnA, fnB) : (__builtin_constant_p (fnA) && ((size_t )(const void *)((fnA) + 1) - (size_t)(const void *)(fnA) == 1 ) && (__s1_len = strlen (fnA), __s1_len < 4) ? (__builtin_constant_p (fnB) && ((size_t)(const void *)((fnB) + 1) - (size_t )(const void *)(fnB) == 1) ? __builtin_strcmp (fnA, fnB) : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (fnB); int __result = (((const unsigned char *) (const char *) (fnA))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (fnA))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (fnA))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (fnA))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (fnB) && ((size_t)(const void *)((fnB) + 1) - (size_t )(const void *)(fnB) == 1) && (__s2_len = strlen (fnB ), __s2_len < 4) ? (__builtin_constant_p (fnA) && ( (size_t)(const void *)((fnA) + 1) - (size_t)(const void *)(fnA ) == 1) ? __builtin_strcmp (fnA, fnB) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (fnA); int __result = (((const unsigned char *) (const char * ) (fnB))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ( fnB))[1] - __s2[1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (fnB ))[2] - __s2[2]); if (__s2_len > 2 && __result == 0 ) __result = (((const unsigned char *) (const char *) (fnB))[ 3] - __s2[3]); } } __result; })))) : __builtin_strcmp (fnA, fnB )))); }) != 0) |
952 | { |
953 | read_posres(mtop, mi, TRUE1, fnB, rc_scaling, ePBC, comB, wi); |
954 | } |
955 | } |
956 | |
957 | static void set_wall_atomtype(gpp_atomtype_t at, t_gromppopts *opts, |
958 | t_inputrec *ir, warninp_t wi) |
959 | { |
960 | int i; |
961 | char warn_buf[STRLEN4096]; |
962 | |
963 | if (ir->nwall > 0) |
964 | { |
965 | fprintf(stderrstderr, "Searching the wall atom type(s)\n"); |
966 | } |
967 | for (i = 0; i < ir->nwall; i++) |
968 | { |
969 | ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i], at); |
970 | if (ir->wall_atomtype[i] == NOTSET-12345) |
971 | { |
972 | sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]); |
973 | warning_error(wi, warn_buf); |
974 | } |
975 | } |
976 | } |
977 | |
978 | static int nrdf_internal(t_atoms *atoms) |
979 | { |
980 | int i, nmass, nrdf; |
981 | |
982 | nmass = 0; |
983 | for (i = 0; i < atoms->nr; i++) |
984 | { |
985 | /* Vsite ptype might not be set here yet, so also check the mass */ |
986 | if ((atoms->atom[i].ptype == eptAtom || |
987 | atoms->atom[i].ptype == eptNucleus) |
988 | && atoms->atom[i].m > 0) |
989 | { |
990 | nmass++; |
991 | } |
992 | } |
993 | switch (nmass) |
994 | { |
995 | case 0: nrdf = 0; break; |
996 | case 1: nrdf = 0; break; |
997 | case 2: nrdf = 1; break; |
998 | default: nrdf = nmass*3 - 6; break; |
999 | } |
1000 | |
1001 | return nrdf; |
1002 | } |
1003 | |
1004 | void |
1005 | spline1d( double dx, |
1006 | double * y, |
1007 | int n, |
1008 | double * u, |
1009 | double * y2 ) |
1010 | { |
1011 | int i; |
1012 | double p, q; |
1013 | |
1014 | y2[0] = 0.0; |
1015 | u[0] = 0.0; |
1016 | |
1017 | for (i = 1; i < n-1; i++) |
1018 | { |
1019 | p = 0.5*y2[i-1]+2.0; |
1020 | y2[i] = -0.5/p; |
1021 | q = (y[i+1]-2.0*y[i]+y[i-1])/dx; |
1022 | u[i] = (3.0*q/dx-0.5*u[i-1])/p; |
1023 | } |
1024 | |
1025 | y2[n-1] = 0.0; |
1026 | |
1027 | for (i = n-2; i >= 0; i--) |
1028 | { |
1029 | y2[i] = y2[i]*y2[i+1]+u[i]; |
1030 | } |
1031 | } |
1032 | |
1033 | |
1034 | void |
1035 | interpolate1d( double xmin, |
1036 | double dx, |
1037 | double * ya, |
1038 | double * y2a, |
1039 | double x, |
1040 | double * y, |
1041 | double * y1) |
1042 | { |
1043 | int ix; |
1044 | double a, b; |
1045 | |
1046 | ix = (x-xmin)/dx; |
1047 | |
1048 | a = (xmin+(ix+1)*dx-x)/dx; |
1049 | b = (x-xmin-ix*dx)/dx; |
1050 | |
1051 | *y = a*ya[ix]+b*ya[ix+1]+((a*a*a-a)*y2a[ix]+(b*b*b-b)*y2a[ix+1])*(dx*dx)/6.0; |
1052 | *y1 = (ya[ix+1]-ya[ix])/dx-(3.0*a*a-1.0)/6.0*dx*y2a[ix]+(3.0*b*b-1.0)/6.0*dx*y2a[ix+1]; |
1053 | } |
1054 | |
1055 | |
1056 | void |
1057 | setup_cmap (int grid_spacing, |
1058 | int nc, |
1059 | real * grid, |
1060 | gmx_cmap_t * cmap_grid) |
1061 | { |
1062 | double *tmp_u, *tmp_u2, *tmp_yy, *tmp_y1, *tmp_t2, *tmp_grid; |
1063 | |
1064 | int i, j, k, ii, jj, kk, idx; |
1065 | int offset; |
1066 | double dx, xmin, v, v1, v2, v12; |
1067 | double phi, psi; |
1068 | |
1069 | snew(tmp_u, 2*grid_spacing)(tmp_u) = save_calloc("tmp_u", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1069, (2*grid_spacing), sizeof(*(tmp_u))); |
1070 | snew(tmp_u2, 2*grid_spacing)(tmp_u2) = save_calloc("tmp_u2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1070, (2*grid_spacing), sizeof(*(tmp_u2))); |
1071 | snew(tmp_yy, 2*grid_spacing)(tmp_yy) = save_calloc("tmp_yy", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1071, (2*grid_spacing), sizeof(*(tmp_yy))); |
1072 | snew(tmp_y1, 2*grid_spacing)(tmp_y1) = save_calloc("tmp_y1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1072, (2*grid_spacing), sizeof(*(tmp_y1))); |
1073 | snew(tmp_t2, 2*grid_spacing*2*grid_spacing)(tmp_t2) = save_calloc("tmp_t2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1073, (2*grid_spacing*2*grid_spacing), sizeof(*(tmp_t2))); |
1074 | snew(tmp_grid, 2*grid_spacing*2*grid_spacing)(tmp_grid) = save_calloc("tmp_grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1074, (2*grid_spacing*2*grid_spacing), sizeof(*(tmp_grid))); |
1075 | |
1076 | dx = 360.0/grid_spacing; |
1077 | xmin = -180.0-dx*grid_spacing/2; |
1078 | |
1079 | for (kk = 0; kk < nc; kk++) |
1080 | { |
1081 | /* Compute an offset depending on which cmap we are using |
1082 | * Offset will be the map number multiplied with the |
1083 | * grid_spacing * grid_spacing * 2 |
1084 | */ |
1085 | offset = kk * grid_spacing * grid_spacing * 2; |
1086 | |
1087 | for (i = 0; i < 2*grid_spacing; i++) |
1088 | { |
1089 | ii = (i+grid_spacing-grid_spacing/2)%grid_spacing; |
1090 | |
1091 | for (j = 0; j < 2*grid_spacing; j++) |
1092 | { |
1093 | jj = (j+grid_spacing-grid_spacing/2)%grid_spacing; |
1094 | tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj]; |
1095 | } |
1096 | } |
1097 | |
1098 | for (i = 0; i < 2*grid_spacing; i++) |
1099 | { |
1100 | spline1d(dx, &(tmp_grid[2*grid_spacing*i]), 2*grid_spacing, tmp_u, &(tmp_t2[2*grid_spacing*i])); |
1101 | } |
1102 | |
1103 | for (i = grid_spacing/2; i < grid_spacing+grid_spacing/2; i++) |
1104 | { |
1105 | ii = i-grid_spacing/2; |
1106 | phi = ii*dx-180.0; |
1107 | |
1108 | for (j = grid_spacing/2; j < grid_spacing+grid_spacing/2; j++) |
1109 | { |
1110 | jj = j-grid_spacing/2; |
1111 | psi = jj*dx-180.0; |
1112 | |
1113 | for (k = 0; k < 2*grid_spacing; k++) |
1114 | { |
1115 | interpolate1d(xmin, dx, &(tmp_grid[2*grid_spacing*k]), |
1116 | &(tmp_t2[2*grid_spacing*k]), psi, &tmp_yy[k], &tmp_y1[k]); |
1117 | } |
1118 | |
1119 | spline1d(dx, tmp_yy, 2*grid_spacing, tmp_u, tmp_u2); |
1120 | interpolate1d(xmin, dx, tmp_yy, tmp_u2, phi, &v, &v1); |
1121 | spline1d(dx, tmp_y1, 2*grid_spacing, tmp_u, tmp_u2); |
1122 | interpolate1d(xmin, dx, tmp_y1, tmp_u2, phi, &v2, &v12); |
1123 | |
1124 | idx = ii*grid_spacing+jj; |
1125 | cmap_grid->cmapdata[kk].cmap[idx*4] = grid[offset+ii*grid_spacing+jj]; |
1126 | cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1; |
1127 | cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2; |
1128 | cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12; |
1129 | } |
1130 | } |
1131 | } |
1132 | } |
1133 | |
1134 | void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing) |
1135 | { |
1136 | int i, k, nelem; |
1137 | |
1138 | cmap_grid->ngrid = ngrid; |
1139 | cmap_grid->grid_spacing = grid_spacing; |
1140 | nelem = cmap_grid->grid_spacing*cmap_grid->grid_spacing; |
1141 | |
1142 | snew(cmap_grid->cmapdata, ngrid)(cmap_grid->cmapdata) = save_calloc("cmap_grid->cmapdata" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1142, (ngrid), sizeof(*(cmap_grid->cmapdata))); |
1143 | |
1144 | for (i = 0; i < cmap_grid->ngrid; i++) |
1145 | { |
1146 | snew(cmap_grid->cmapdata[i].cmap, 4*nelem)(cmap_grid->cmapdata[i].cmap) = save_calloc("cmap_grid->cmapdata[i].cmap" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1146, (4*nelem), sizeof(*(cmap_grid->cmapdata[i].cmap))); |
1147 | } |
1148 | } |
1149 | |
1150 | |
1151 | static int count_constraints(gmx_mtop_t *mtop, t_molinfo *mi, warninp_t wi) |
1152 | { |
1153 | int count, count_mol, i, mb; |
1154 | gmx_molblock_t *molb; |
1155 | t_params *plist; |
1156 | char buf[STRLEN4096]; |
1157 | |
1158 | count = 0; |
1159 | for (mb = 0; mb < mtop->nmolblock; mb++) |
1160 | { |
1161 | count_mol = 0; |
1162 | molb = &mtop->molblock[mb]; |
1163 | plist = mi[molb->type].plist; |
1164 | |
1165 | for (i = 0; i < F_NRE; i++) |
1166 | { |
1167 | if (i == F_SETTLE) |
1168 | { |
1169 | count_mol += 3*plist[i].nr; |
1170 | } |
1171 | else if (interaction_function[i].flags & IF_CONSTRAINT1<<2) |
1172 | { |
1173 | count_mol += plist[i].nr; |
1174 | } |
1175 | } |
1176 | |
1177 | if (count_mol > nrdf_internal(&mi[molb->type].atoms)) |
1178 | { |
1179 | sprintf(buf, |
1180 | "Molecule type '%s' has %d constraints.\n" |
1181 | "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n", |
1182 | *mi[molb->type].name, count_mol, |
1183 | nrdf_internal(&mi[molb->type].atoms)); |
1184 | warning(wi, buf); |
1185 | } |
1186 | count += molb->nmol*count_mol; |
1187 | } |
1188 | |
1189 | return count; |
1190 | } |
1191 | |
1192 | static void check_gbsa_params_charged(gmx_mtop_t *sys, gpp_atomtype_t atype) |
1193 | { |
1194 | int i, nmiss, natoms, mt; |
1195 | real q; |
1196 | const t_atoms *atoms; |
1197 | |
1198 | nmiss = 0; |
1199 | for (mt = 0; mt < sys->nmoltype; mt++) |
1200 | { |
1201 | atoms = &sys->moltype[mt].atoms; |
1202 | natoms = atoms->nr; |
1203 | |
1204 | for (i = 0; i < natoms; i++) |
1205 | { |
1206 | q = atoms->atom[i].q; |
1207 | if ((get_atomtype_radius(atoms->atom[i].type, atype) == 0 || |
1208 | get_atomtype_vol(atoms->atom[i].type, atype) == 0 || |
1209 | get_atomtype_surftens(atoms->atom[i].type, atype) == 0 || |
1210 | get_atomtype_gb_radius(atoms->atom[i].type, atype) == 0 || |
1211 | get_atomtype_S_hct(atoms->atom[i].type, atype) == 0) && |
1212 | q != 0) |
1213 | { |
1214 | fprintf(stderrstderr, "\nGB parameter(s) zero for atom type '%s' while charge is %g\n", |
1215 | get_atomtype_name(atoms->atom[i].type, atype), q); |
1216 | nmiss++; |
1217 | } |
1218 | } |
1219 | } |
1220 | |
1221 | if (nmiss > 0) |
1222 | { |
1223 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1223, "Can't do GB electrostatics; the implicit_genborn_params section of the forcefield has parameters with value zero for %d atomtypes that occur as charged atoms.", nmiss); |
1224 | } |
1225 | } |
1226 | |
1227 | |
1228 | static void check_gbsa_params(gpp_atomtype_t atype) |
1229 | { |
1230 | int nmiss, i; |
1231 | |
1232 | /* If we are doing GBSA, check that we got the parameters we need |
1233 | * This checking is to see if there are GBSA paratmeters for all |
1234 | * atoms in the force field. To go around this for testing purposes |
1235 | * comment out the nerror++ counter temporarily |
1236 | */ |
1237 | nmiss = 0; |
1238 | for (i = 0; i < get_atomtype_ntypes(atype); i++) |
1239 | { |
1240 | if (get_atomtype_radius(i, atype) < 0 || |
1241 | get_atomtype_vol(i, atype) < 0 || |
1242 | get_atomtype_surftens(i, atype) < 0 || |
1243 | get_atomtype_gb_radius(i, atype) < 0 || |
1244 | get_atomtype_S_hct(i, atype) < 0) |
1245 | { |
1246 | fprintf(stderrstderr, "\nGB parameter(s) missing or negative for atom type '%s'\n", |
1247 | get_atomtype_name(i, atype)); |
1248 | nmiss++; |
1249 | } |
1250 | } |
1251 | |
1252 | if (nmiss > 0) |
1253 | { |
1254 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1254, "Can't do GB electrostatics; the implicit_genborn_params section of the forcefield is missing parameters for %d atomtypes or they might be negative.", nmiss); |
1255 | } |
1256 | |
1257 | } |
1258 | |
1259 | static real calc_temp(const gmx_mtop_t *mtop, |
1260 | const t_inputrec *ir, |
1261 | rvec *v) |
1262 | { |
1263 | double sum_mv2; |
1264 | gmx_mtop_atomloop_all_t aloop; |
1265 | t_atom *atom; |
1266 | int a; |
1267 | int nrdf, g; |
1268 | |
1269 | sum_mv2 = 0; |
1270 | |
1271 | aloop = gmx_mtop_atomloop_all_init(mtop); |
1272 | while (gmx_mtop_atomloop_all_next(aloop, &a, &atom)) |
1273 | { |
1274 | sum_mv2 += atom->m*norm2(v[a]); |
1275 | } |
1276 | |
1277 | nrdf = 0; |
1278 | for (g = 0; g < ir->opts.ngtc; g++) |
1279 | { |
1280 | nrdf += ir->opts.nrdf[g]; |
1281 | } |
1282 | |
1283 | return sum_mv2/(nrdf*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))); |
1284 | } |
1285 | |
1286 | static real get_max_reference_temp(const t_inputrec *ir, |
1287 | warninp_t wi) |
1288 | { |
1289 | real ref_t; |
1290 | int i; |
1291 | gmx_bool bNoCoupl; |
1292 | |
1293 | ref_t = 0; |
1294 | bNoCoupl = FALSE0; |
1295 | for (i = 0; i < ir->opts.ngtc; i++) |
1296 | { |
1297 | if (ir->opts.tau_t[i] < 0) |
1298 | { |
1299 | bNoCoupl = TRUE1; |
1300 | } |
1301 | else |
1302 | { |
1303 | ref_t = max(ref_t, ir->opts.ref_t[i])(((ref_t) > (ir->opts.ref_t[i])) ? (ref_t) : (ir->opts .ref_t[i]) ); |
1304 | } |
1305 | } |
1306 | |
1307 | if (bNoCoupl) |
1308 | { |
1309 | char buf[STRLEN4096]; |
1310 | |
1311 | sprintf(buf, "Some temperature coupling groups do not use temperature coupling. We will assume their temperature is not more than %.3f K. If their temperature is higher, the energy error and the Verlet buffer might be underestimated.", |
1312 | ref_t); |
1313 | warning(wi, buf); |
1314 | } |
1315 | |
1316 | return ref_t; |
1317 | } |
1318 | |
1319 | static void set_verlet_buffer(const gmx_mtop_t *mtop, |
1320 | t_inputrec *ir, |
1321 | real buffer_temp, |
1322 | matrix box, |
1323 | warninp_t wi) |
1324 | { |
1325 | int i; |
1326 | verletbuf_list_setup_t ls; |
1327 | real rlist_1x1; |
1328 | int n_nonlin_vsite; |
1329 | char warn_buf[STRLEN4096]; |
1330 | |
1331 | printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp); |
1332 | |
1333 | /* Calculate the buffer size for simple atom vs atoms list */ |
1334 | ls.cluster_size_i = 1; |
1335 | ls.cluster_size_j = 1; |
1336 | calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp, |
1337 | &ls, &n_nonlin_vsite, &rlist_1x1); |
1338 | |
1339 | /* Set the pair-list buffer size in ir */ |
1340 | verletbuf_get_list_setup(FALSE0, &ls); |
1341 | calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp, |
1342 | &ls, &n_nonlin_vsite, &ir->rlist); |
1343 | |
1344 | if (n_nonlin_vsite > 0) |
1345 | { |
1346 | sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy error is approximated. In most cases this does not affect the error significantly.", n_nonlin_vsite); |
1347 | warning_note(wi, warn_buf); |
1348 | } |
1349 | |
1350 | printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n", |
1351 | 1, 1, rlist_1x1, rlist_1x1-max(ir->rvdw, ir->rcoulomb)(((ir->rvdw) > (ir->rcoulomb)) ? (ir->rvdw) : (ir ->rcoulomb) )); |
1352 | |
1353 | ir->rlistlong = ir->rlist; |
1354 | printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n", |
1355 | ls.cluster_size_i, ls.cluster_size_j, |
1356 | ir->rlist, ir->rlist-max(ir->rvdw, ir->rcoulomb)(((ir->rvdw) > (ir->rcoulomb)) ? (ir->rvdw) : (ir ->rcoulomb) )); |
1357 | |
1358 | if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box)) |
1359 | { |
1360 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1360, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, sqrt(max_cutoff2(ir->ePBC, box))); |
1361 | } |
1362 | } |
1363 | |
1364 | int gmx_grompp(int argc, char *argv[]) |
1365 | { |
1366 | static const char *desc[] = { |
1367 | "[THISMODULE] (the gromacs preprocessor)", |
1368 | "reads a molecular topology file, checks the validity of the", |
1369 | "file, expands the topology from a molecular description to an atomic", |
1370 | "description. The topology file contains information about", |
1371 | "molecule types and the number of molecules, the preprocessor", |
1372 | "copies each molecule as needed. ", |
1373 | "There is no limitation on the number of molecule types. ", |
1374 | "Bonds and bond-angles can be converted into constraints, separately", |
1375 | "for hydrogens and heavy atoms.", |
1376 | "Then a coordinate file is read and velocities can be generated", |
1377 | "from a Maxwellian distribution if requested.", |
1378 | "[THISMODULE] also reads parameters for [gmx-mdrun] ", |
1379 | "(eg. number of MD steps, time step, cut-off), and others such as", |
1380 | "NEMD parameters, which are corrected so that the net acceleration", |
1381 | "is zero.", |
1382 | "Eventually a binary file is produced that can serve as the sole input", |
1383 | "file for the MD program.[PAR]", |
1384 | |
1385 | "[THISMODULE] uses the atom names from the topology file. The atom names", |
1386 | "in the coordinate file (option [TT]-c[tt]) are only read to generate", |
1387 | "warnings when they do not match the atom names in the topology.", |
1388 | "Note that the atom names are irrelevant for the simulation as", |
1389 | "only the atom types are used for generating interaction parameters.[PAR]", |
1390 | |
1391 | "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ", |
1392 | "etc. The preprocessor supports the following keywords:[PAR]", |
1393 | "#ifdef VARIABLE[BR]", |
1394 | "#ifndef VARIABLE[BR]", |
1395 | "#else[BR]", |
1396 | "#endif[BR]", |
1397 | "#define VARIABLE[BR]", |
1398 | "#undef VARIABLE[BR]" |
1399 | "#include \"filename\"[BR]", |
1400 | "#include <filename>[PAR]", |
1401 | "The functioning of these statements in your topology may be modulated by", |
1402 | "using the following two flags in your [TT].mdp[tt] file:[PAR]", |
1403 | "[TT]define = -DVARIABLE1 -DVARIABLE2[BR]", |
1404 | "include = -I/home/john/doe[tt][BR]", |
1405 | "For further information a C-programming textbook may help you out.", |
1406 | "Specifying the [TT]-pp[tt] flag will get the pre-processed", |
1407 | "topology file written out so that you can verify its contents.[PAR]", |
1408 | |
1409 | /* cpp has been unnecessary for some time, hasn't it? |
1410 | "If your system does not have a C-preprocessor, you can still", |
1411 | "use [TT]grompp[tt], but you do not have access to the features ", |
1412 | "from the cpp. Command line options to the C-preprocessor can be given", |
1413 | "in the [TT].mdp[tt] file. See your local manual (man cpp).[PAR]", |
1414 | */ |
1415 | |
1416 | "When using position restraints a file with restraint coordinates", |
1417 | "can be supplied with [TT]-r[tt], otherwise restraining will be done", |
1418 | "with respect to the conformation from the [TT]-c[tt] option.", |
1419 | "For free energy calculation the the coordinates for the B topology", |
1420 | "can be supplied with [TT]-rb[tt], otherwise they will be equal to", |
1421 | "those of the A topology.[PAR]", |
1422 | |
1423 | "Starting coordinates can be read from trajectory with [TT]-t[tt].", |
1424 | "The last frame with coordinates and velocities will be read,", |
1425 | "unless the [TT]-time[tt] option is used. Only if this information", |
1426 | "is absent will the coordinates in the [TT]-c[tt] file be used.", |
1427 | "Note that these velocities will not be used when [TT]gen_vel = yes[tt]", |
1428 | "in your [TT].mdp[tt] file. An energy file can be supplied with", |
1429 | "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling", |
1430 | "variables.[PAR]", |
1431 | |
1432 | "[THISMODULE] can be used to restart simulations (preserving", |
1433 | "continuity) by supplying just a checkpoint file with [TT]-t[tt].", |
1434 | "However, for simply changing the number of run steps to extend", |
1435 | "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].", |
1436 | "You then supply the old checkpoint file directly to [gmx-mdrun]", |
1437 | "with [TT]-cpi[tt]. If you wish to change the ensemble or things", |
1438 | "like output frequency, then supplying the checkpoint file to", |
1439 | "[THISMODULE] with [TT]-t[tt] along with a new [TT].mdp[tt] file", |
1440 | "with [TT]-f[tt] is the recommended procedure.[PAR]", |
1441 | |
1442 | "By default, all bonded interactions which have constant energy due to", |
1443 | "virtual site constructions will be removed. If this constant energy is", |
1444 | "not zero, this will result in a shift in the total energy. All bonded", |
1445 | "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,", |
1446 | "all constraints for distances which will be constant anyway because", |
1447 | "of virtual site constructions will be removed. If any constraints remain", |
1448 | "which involve virtual sites, a fatal error will result.[PAR]" |
1449 | |
1450 | "To verify your run input file, please take note of all warnings", |
1451 | "on the screen, and correct where necessary. Do also look at the contents", |
1452 | "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as", |
1453 | "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]", |
1454 | "with the [TT]-debug[tt] option which will give you more information", |
1455 | "in a file called [TT]grompp.log[tt] (along with real debug info). You", |
1456 | "can see the contents of the run input file with the [gmx-dump]", |
1457 | "program. [gmx-check] can be used to compare the contents of two", |
1458 | "run input files.[PAR]" |
1459 | |
1460 | "The [TT]-maxwarn[tt] option can be used to override warnings printed", |
1461 | "by [THISMODULE] that otherwise halt output. In some cases, warnings are", |
1462 | "harmless, but usually they are not. The user is advised to carefully", |
1463 | "interpret the output messages before attempting to bypass them with", |
1464 | "this option." |
1465 | }; |
1466 | t_gromppopts *opts; |
1467 | gmx_mtop_t *sys; |
1468 | int nmi; |
1469 | t_molinfo *mi; |
1470 | gpp_atomtype_t atype; |
1471 | t_inputrec *ir; |
1472 | int natoms, nvsite, comb, mt; |
1473 | t_params *plist; |
1474 | t_state state; |
1475 | matrix box; |
1476 | real max_spacing, fudgeQQ; |
1477 | double reppow; |
1478 | char fn[STRLEN4096], fnB[STRLEN4096]; |
1479 | const char *mdparin; |
1480 | int ntype; |
1481 | gmx_bool bNeedVel, bGenVel; |
1482 | gmx_bool have_atomnumber; |
1483 | int n12, n13, n14; |
1484 | t_params *gb_plist = NULL((void*)0); |
1485 | gmx_genborn_t *born = NULL((void*)0); |
1486 | output_env_t oenv; |
1487 | gmx_bool bVerbose = FALSE0; |
1488 | warninp_t wi; |
1489 | char warn_buf[STRLEN4096]; |
1490 | unsigned int useed; |
1491 | t_atoms IMDatoms; /* Atoms to be operated on interactively (IMD) */ |
1492 | |
1493 | t_filenm fnm[] = { |
1494 | { efMDP, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
1495 | { efMDP, "-po", "mdout", ffWRITE1<<2 }, |
1496 | { efSTX, "-c", NULL((void*)0), ffREAD1<<1 }, |
1497 | { efSTX, "-r", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1498 | { efSTX, "-rb", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1499 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1500 | { efTOP, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
1501 | { efTOP, "-pp", "processed", ffOPTWR(1<<2| 1<<3) }, |
1502 | { efTPX, "-o", NULL((void*)0), ffWRITE1<<2 }, |
1503 | { efTRN, "-t", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1504 | { efEDR, "-e", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1505 | /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */ |
1506 | { efGRO, "-imd", "imdgroup", ffOPTWR(1<<2| 1<<3) }, |
1507 | { efTRN, "-ref", "rotref", ffOPTRW((1<<1 | 1<<2) | 1<<3) } |
1508 | }; |
1509 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
1510 | |
1511 | /* Command line options */ |
1512 | static gmx_bool bRenum = TRUE1; |
1513 | static gmx_bool bRmVSBds = TRUE1, bZero = FALSE0; |
1514 | static int i, maxwarn = 0; |
1515 | static real fr_time = -1; |
1516 | t_pargs pa[] = { |
1517 | { "-v", FALSE0, etBOOL, {&bVerbose}, |
1518 | "Be loud and noisy" }, |
1519 | { "-time", FALSE0, etREAL, {&fr_time}, |
1520 | "Take frame at or first after this time." }, |
1521 | { "-rmvsbds", FALSE0, etBOOL, {&bRmVSBds}, |
1522 | "Remove constant bonded interactions with virtual sites" }, |
1523 | { "-maxwarn", FALSE0, etINT, {&maxwarn}, |
1524 | "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" }, |
1525 | { "-zero", FALSE0, etBOOL, {&bZero}, |
1526 | "Set parameters for bonded interactions without defaults to zero instead of generating an error" }, |
1527 | { "-renum", FALSE0, etBOOL, {&bRenum}, |
1528 | "Renumber atomtypes and minimize number of atomtypes" } |
1529 | }; |
1530 | |
1531 | /* Initiate some variables */ |
1532 | snew(ir, 1)(ir) = save_calloc("ir", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1532, (1), sizeof(*(ir))); |
1533 | snew(opts, 1)(opts) = save_calloc("opts", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1533, (1), sizeof(*(opts))); |
1534 | init_ir(ir, opts); |
1535 | |
1536 | /* Parse the command line */ |
1537 | if (!parse_common_args(&argc, argv, 0, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, |
1538 | asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
1539 | { |
1540 | return 0; |
1541 | } |
1542 | |
1543 | wi = init_warning(TRUE1, maxwarn); |
1544 | |
1545 | /* PARAMETER file processing */ |
1546 | mdparin = opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1547 | set_warning_line(wi, mdparin, -1); |
1548 | get_ir(mdparin, opt2fn("-po", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ir, opts, wi); |
1549 | |
1550 | if (bVerbose) |
1551 | { |
1552 | fprintf(stderrstderr, "checking input for internal consistency...\n"); |
1553 | } |
1554 | check_ir(mdparin, ir, opts, wi); |
1555 | |
1556 | if (ir->ld_seed == -1) |
1557 | { |
1558 | ir->ld_seed = (gmx_int64_t)gmx_rng_make_seed(); |
1559 | fprintf(stderrstderr, "Setting the LD random seed to %"GMX_PRId64"l" "d" "\n", ir->ld_seed); |
1560 | } |
1561 | |
1562 | if (ir->expandedvals->lmc_seed == -1) |
1563 | { |
1564 | ir->expandedvals->lmc_seed = (int)gmx_rng_make_seed(); |
1565 | fprintf(stderrstderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed); |
1566 | } |
1567 | |
1568 | bNeedVel = EI_STATE_VELOCITY(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2)); |
1569 | bGenVel = (bNeedVel && opts->bGenVel); |
1570 | if (bGenVel && ir->bContinuation) |
1571 | { |
1572 | sprintf(warn_buf, |
1573 | "Generating velocities is inconsistent with attempting " |
1574 | "to continue a previous run. Choose only one of " |
1575 | "gen-vel = yes and continuation = yes."); |
1576 | warning_error(wi, warn_buf); |
1577 | } |
1578 | |
1579 | snew(plist, F_NRE)(plist) = save_calloc("plist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1579, (F_NRE), sizeof(*(plist))); |
1580 | init_plist(plist); |
1581 | snew(sys, 1)(sys) = save_calloc("sys", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1581, (1), sizeof(*(sys))); |
1582 | atype = init_atomtype(); |
1583 | if (debug) |
1584 | { |
1585 | pr_symtab(debug, 0, "Just opened", &sys->symtab); |
1586 | } |
1587 | |
1588 | strcpy(fn, ftp2fn(efTOP, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
1589 | if (!gmx_fexist(fn)) |
1590 | { |
1591 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1591, "%s does not exist", fn); |
1592 | } |
1593 | new_status(fn, opt2fn_null("-pp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-c", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
1594 | opts, ir, bZero, bGenVel, bVerbose, &state, |
1595 | atype, sys, &nmi, &mi, plist, &comb, &reppow, &fudgeQQ, |
1596 | opts->bMorse, |
1597 | wi); |
1598 | |
1599 | if (debug) |
1600 | { |
1601 | pr_symtab(debug, 0, "After new_status", &sys->symtab); |
1602 | } |
1603 | |
1604 | nvsite = 0; |
1605 | /* set parameters for virtual site construction (not for vsiten) */ |
1606 | for (mt = 0; mt < sys->nmoltype; mt++) |
1607 | { |
1608 | nvsite += |
1609 | set_vsites(bVerbose, &sys->moltype[mt].atoms, atype, mi[mt].plist); |
1610 | } |
1611 | /* now throw away all obsolete bonds, angles and dihedrals: */ |
1612 | /* note: constraints are ALWAYS removed */ |
1613 | if (nvsite) |
1614 | { |
1615 | for (mt = 0; mt < sys->nmoltype; mt++) |
1616 | { |
1617 | clean_vsite_bondeds(mi[mt].plist, sys->moltype[mt].atoms.nr, bRmVSBds); |
1618 | } |
1619 | } |
1620 | |
1621 | if (ir->cutoff_scheme == ecutsVERLET) |
1622 | { |
1623 | fprintf(stderrstderr, "Removing all charge groups because cutoff-scheme=%s\n", |
1624 | ecutscheme_names[ir->cutoff_scheme]); |
1625 | |
1626 | /* Remove all charge groups */ |
1627 | gmx_mtop_remove_chargegroups(sys); |
1628 | } |
1629 | |
1630 | if (count_constraints(sys, mi, wi) && (ir->eConstrAlg == econtSHAKE)) |
1631 | { |
1632 | if (ir->eI == eiCG || ir->eI == eiLBFGS) |
1633 | { |
1634 | sprintf(warn_buf, "Can not do %s with %s, use %s", |
1635 | EI(ir->eI)((((ir->eI) < 0) || ((ir->eI) >= (eiNR))) ? "UNDEFINED" : (ei_names)[ir->eI]), econstr_names[econtSHAKE], econstr_names[econtLINCS]); |
1636 | warning_error(wi, warn_buf); |
1637 | } |
1638 | if (ir->bPeriodicMols) |
1639 | { |
1640 | sprintf(warn_buf, "Can not do periodic molecules with %s, use %s", |
1641 | econstr_names[econtSHAKE], econstr_names[econtLINCS]); |
1642 | warning_error(wi, warn_buf); |
1643 | } |
1644 | } |
1645 | |
1646 | if (EI_SD (ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2) && ir->etc != etcNO) |
1647 | { |
1648 | warning_note(wi, "Temperature coupling is ignored with SD integrators."); |
1649 | } |
1650 | |
1651 | /* If we are doing QM/MM, check that we got the atom numbers */ |
1652 | have_atomnumber = TRUE1; |
1653 | for (i = 0; i < get_atomtype_ntypes(atype); i++) |
1654 | { |
1655 | have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i, atype) >= 0); |
1656 | } |
1657 | if (!have_atomnumber && ir->bQMMM) |
1658 | { |
1659 | warning_error(wi, |
1660 | "\n" |
1661 | "It appears as if you are trying to run a QM/MM calculation, but the force\n" |
1662 | "field you are using does not contain atom numbers fields. This is an\n" |
1663 | "optional field (introduced in Gromacs 3.3) for general runs, but mandatory\n" |
1664 | "for QM/MM. The good news is that it is easy to add - put the atom number as\n" |
1665 | "an integer just before the mass column in ffXXXnb.itp.\n" |
1666 | "NB: United atoms have the same atom numbers as normal ones.\n\n"); |
1667 | } |
1668 | |
1669 | if (ir->bAdress) |
1670 | { |
1671 | if ((ir->adress->const_wf > 1) || (ir->adress->const_wf < 0)) |
1672 | { |
1673 | warning_error(wi, "AdResS contant weighting function should be between 0 and 1\n\n"); |
1674 | } |
1675 | /** TODO check size of ex+hy width against box size */ |
1676 | } |
1677 | |
1678 | /* Check for errors in the input now, since they might cause problems |
1679 | * during processing further down. |
1680 | */ |
1681 | check_warning_error(wi, FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1681); |
1682 | |
1683 | if (opt2bSet("-r", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
1684 | { |
1685 | sprintf(fn, "%s", opt2fn("-r", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
1686 | } |
1687 | else |
1688 | { |
1689 | sprintf(fn, "%s", opt2fn("-c", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
1690 | } |
1691 | if (opt2bSet("-rb", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
1692 | { |
1693 | sprintf(fnB, "%s", opt2fn("-rb", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
1694 | } |
1695 | else |
1696 | { |
1697 | strcpy(fnB, fn); |
1698 | } |
1699 | |
1700 | if (nint_ftype(sys, mi, F_POSRES) > 0 || nint_ftype(sys, mi, F_FBPOSRES) > 0) |
1701 | { |
1702 | if (bVerbose) |
1703 | { |
1704 | fprintf(stderrstderr, "Reading position restraint coords from %s", fn); |
1705 | if (strcmp(fn, fnB)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (fn) && __builtin_constant_p (fnB) && (__s1_len = strlen (fn), __s2_len = strlen (fnB), (!((size_t)(const void *)((fn) + 1) - (size_t)(const void *)(fn) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((fnB) + 1) - (size_t )(const void *)(fnB) == 1) || __s2_len >= 4)) ? __builtin_strcmp (fn, fnB) : (__builtin_constant_p (fn) && ((size_t)( const void *)((fn) + 1) - (size_t)(const void *)(fn) == 1) && (__s1_len = strlen (fn), __s1_len < 4) ? (__builtin_constant_p (fnB) && ((size_t)(const void *)((fnB) + 1) - (size_t )(const void *)(fnB) == 1) ? __builtin_strcmp (fn, fnB) : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (fnB); int __result = (((const unsigned char *) (const char *) (fn))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ( fn))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (fn ))[2] - __s2[2]); if (__s1_len > 2 && __result == 0 ) __result = (((const unsigned char *) (const char *) (fn))[3 ] - __s2[3]); } } __result; }))) : (__builtin_constant_p (fnB ) && ((size_t)(const void *)((fnB) + 1) - (size_t)(const void *)(fnB) == 1) && (__s2_len = strlen (fnB), __s2_len < 4) ? (__builtin_constant_p (fn) && ((size_t)(const void *)((fn) + 1) - (size_t)(const void *)(fn) == 1) ? __builtin_strcmp (fn, fnB) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (fn); int __result = ( ((const unsigned char *) (const char *) (fnB))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = ( ((const unsigned char *) (const char *) (fnB))[1] - __s2[1]); if (__s2_len > 1 && __result == 0) { __result = ( ((const unsigned char *) (const char *) (fnB))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = ((( const unsigned char *) (const char *) (fnB))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (fn, fnB)))); }) == 0) |
1706 | { |
1707 | fprintf(stderrstderr, "\n"); |
1708 | } |
1709 | else |
1710 | { |
1711 | fprintf(stderrstderr, " and %s\n", fnB); |
1712 | } |
1713 | } |
1714 | gen_posres(sys, mi, fn, fnB, |
1715 | ir->refcoord_scaling, ir->ePBC, |
1716 | ir->posres_com, ir->posres_comB, |
1717 | wi); |
1718 | } |
1719 | |
1720 | /* If we are using CMAP, setup the pre-interpolation grid */ |
1721 | if (plist->ncmap > 0) |
1722 | { |
1723 | init_cmap_grid(&sys->ffparams.cmap_grid, plist->nc, plist->grid_spacing); |
1724 | setup_cmap(plist->grid_spacing, plist->nc, plist->cmap, &sys->ffparams.cmap_grid); |
1725 | } |
1726 | |
1727 | set_wall_atomtype(atype, opts, ir, wi); |
1728 | if (bRenum) |
1729 | { |
1730 | renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose); |
1731 | ntype = get_atomtype_ntypes(atype); |
1732 | } |
1733 | |
1734 | if (ir->implicit_solvent != eisNO) |
1735 | { |
1736 | /* Now we have renumbered the atom types, we can check the GBSA params */ |
1737 | check_gbsa_params(atype); |
1738 | |
1739 | /* Check that all atoms that have charge and/or LJ-parameters also have |
1740 | * sensible GB-parameters |
1741 | */ |
1742 | check_gbsa_params_charged(sys, atype); |
1743 | } |
1744 | |
1745 | /* PELA: Copy the atomtype data to the topology atomtype list */ |
1746 | copy_atomtype_atomtypes(atype, &(sys->atomtypes)); |
1747 | |
1748 | if (debug) |
1749 | { |
1750 | pr_symtab(debug, 0, "After renum_atype", &sys->symtab); |
1751 | } |
1752 | |
1753 | if (bVerbose) |
1754 | { |
1755 | fprintf(stderrstderr, "converting bonded parameters...\n"); |
1756 | } |
1757 | |
1758 | ntype = get_atomtype_ntypes(atype); |
1759 | convert_params(ntype, plist, mi, comb, reppow, fudgeQQ, sys); |
1760 | |
1761 | if (debug) |
1762 | { |
1763 | pr_symtab(debug, 0, "After convert_params", &sys->symtab); |
1764 | } |
1765 | |
1766 | /* set ptype to VSite for virtual sites */ |
1767 | for (mt = 0; mt < sys->nmoltype; mt++) |
1768 | { |
1769 | set_vsites_ptype(FALSE0, &sys->moltype[mt]); |
1770 | } |
1771 | if (debug) |
1772 | { |
1773 | pr_symtab(debug, 0, "After virtual sites", &sys->symtab); |
1774 | } |
1775 | /* Check velocity for virtual sites and shells */ |
1776 | if (bGenVel) |
1777 | { |
1778 | check_vel(sys, state.v); |
1779 | } |
1780 | |
1781 | /* check masses */ |
1782 | check_mol(sys, wi); |
1783 | |
1784 | for (i = 0; i < sys->nmoltype; i++) |
1785 | { |
1786 | check_cg_sizes(ftp2fn(efTOP, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &sys->moltype[i].cgs, wi); |
1787 | } |
1788 | |
1789 | if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || (ir->eI) == eiBD) && ir->eI != eiBD) |
1790 | { |
1791 | check_bonds_timestep(sys, ir->delta_t, wi); |
1792 | } |
1793 | |
1794 | if (EI_ENERGY_MINIMIZATION(ir->eI)((ir->eI) == eiSteep || (ir->eI) == eiCG || (ir->eI) == eiLBFGS) && 0 == ir->nsteps) |
1795 | { |
1796 | warning_note(wi, "Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun."); |
1797 | } |
1798 | |
1799 | check_warning_error(wi, FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1799); |
1800 | |
1801 | if (bVerbose) |
1802 | { |
1803 | fprintf(stderrstderr, "initialising group options...\n"); |
1804 | } |
1805 | do_index(mdparin, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
1806 | sys, bVerbose, ir, |
1807 | bGenVel ? state.v : NULL((void*)0), |
1808 | wi); |
1809 | |
1810 | if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && |
1811 | ir->nstlist > 1) |
1812 | { |
1813 | if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || (ir->eI) == eiBD) && inputrec2nboundeddim(ir) == 3) |
1814 | { |
1815 | real buffer_temp; |
1816 | |
1817 | if (EI_MD(ir->eI)((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) && ir->etc == etcNO) |
1818 | { |
1819 | if (bGenVel) |
1820 | { |
1821 | buffer_temp = opts->tempi; |
1822 | } |
1823 | else |
1824 | { |
1825 | buffer_temp = calc_temp(sys, ir, state.v); |
1826 | } |
1827 | if (buffer_temp > 0) |
1828 | { |
1829 | sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp); |
1830 | warning_note(wi, warn_buf); |
1831 | } |
1832 | else |
1833 | { |
1834 | sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!", |
1835 | (int)(verlet_buffer_ratio_NVE_T0*100 + 0.5)); |
1836 | warning_note(wi, warn_buf); |
1837 | } |
1838 | } |
1839 | else |
1840 | { |
1841 | buffer_temp = get_max_reference_temp(ir, wi); |
1842 | } |
1843 | |
1844 | if (EI_MD(ir->eI)((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) && ir->etc == etcNO && buffer_temp == 0) |
1845 | { |
1846 | /* NVE with initial T=0: we add a fixed ratio to rlist. |
1847 | * Since we don't actually use verletbuf_tol, we set it to -1 |
1848 | * so it can't be misused later. |
1849 | */ |
1850 | ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0; |
1851 | ir->verletbuf_tol = -1; |
1852 | } |
1853 | else |
1854 | { |
1855 | /* We warn for NVE simulations with >1(.1)% drift tolerance */ |
1856 | const real drift_tol = 0.01; |
1857 | real ener_runtime; |
1858 | |
1859 | /* We use 2 DOF per atom = 2kT pot+kin energy, to be on |
1860 | * the safe side with constraints (without constraints: 3 DOF). |
1861 | */ |
1862 | ener_runtime = 2*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*buffer_temp/(ir->nsteps*ir->delta_t); |
1863 | |
1864 | if (EI_MD(ir->eI)((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) && ir->etc == etcNO && ir->nstlist > 1 && |
1865 | ir->nsteps > 0 && |
1866 | ir->verletbuf_tol > 1.1*drift_tol*ener_runtime) |
1867 | { |
1868 | sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%%, you might need to set verlet-buffer-tolerance to %.1e.", |
1869 | ir->verletbuf_tol, ir->nsteps*ir->delta_t, |
1870 | (int)(ir->verletbuf_tol/ener_runtime*100 + 0.5), |
1871 | (int)(100*drift_tol + 0.5), |
1872 | drift_tol*ener_runtime); |
1873 | warning_note(wi, warn_buf); |
1874 | } |
1875 | |
1876 | set_verlet_buffer(sys, ir, buffer_temp, state.box, wi); |
1877 | } |
1878 | } |
1879 | } |
1880 | |
1881 | /* Init the temperature coupling state */ |
1882 | init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */ |
1883 | |
1884 | if (bVerbose) |
1885 | { |
1886 | fprintf(stderrstderr, "Checking consistency between energy and charge groups...\n"); |
1887 | } |
1888 | check_eg_vs_cg(sys); |
1889 | |
1890 | if (debug) |
1891 | { |
1892 | pr_symtab(debug, 0, "After index", &sys->symtab); |
1893 | } |
1894 | |
1895 | triple_check(mdparin, ir, sys, wi); |
1896 | close_symtab(&sys->symtab); |
1897 | if (debug) |
1898 | { |
1899 | pr_symtab(debug, 0, "After close", &sys->symtab); |
1900 | } |
1901 | |
1902 | /* make exclusions between QM atoms */ |
1903 | if (ir->bQMMM) |
1904 | { |
1905 | if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE) |
1906 | { |
1907 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1907, "electrostatic embedding only works with grid neighboursearching, use ns-type=grid instead\n"); |
1908 | } |
1909 | else |
1910 | { |
1911 | generate_qmexcl(sys, ir, wi); |
1912 | } |
1913 | } |
1914 | |
1915 | if (ftp2bSet(efTRN, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
1916 | { |
1917 | if (bVerbose) |
1918 | { |
1919 | fprintf(stderrstderr, "getting data from old trajectory ...\n"); |
1920 | } |
1921 | cont_status(ftp2fn(efTRN, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ftp2fn_null(efEDR, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
1922 | bNeedVel, bGenVel, fr_time, ir, &state, sys, oenv); |
1923 | } |
1924 | |
1925 | if (ir->ePBC == epbcXY && ir->nwall != 2) |
1926 | { |
1927 | clear_rvec(state.box[ZZ2]); |
1928 | } |
1929 | |
1930 | if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0) |
1931 | { |
1932 | set_warning_line(wi, mdparin, -1); |
1933 | check_chargegroup_radii(sys, ir, state.x, wi); |
1934 | } |
1935 | |
1936 | if (EEL_FULL(ir->coulombtype)((((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir->coulombtype ) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD) || (ir->coulombtype) == eelEWALD) || (ir->coulombtype) == eelPOISSON) || EVDW_PME(ir->vdwtype)((ir->vdwtype) == evdwPME)) |
1937 | { |
1938 | /* Calculate the optimal grid dimensions */ |
1939 | copy_mat(state.box, box); |
1940 | if (ir->ePBC == epbcXY && ir->nwall == 2) |
1941 | { |
1942 | svmul(ir->wall_ewald_zfac, box[ZZ2], box[ZZ2]); |
1943 | } |
1944 | if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0) |
1945 | { |
1946 | /* Mark fourier_spacing as not used */ |
1947 | ir->fourier_spacing = 0; |
1948 | } |
1949 | else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0) |
1950 | { |
1951 | set_warning_line(wi, mdparin, -1); |
1952 | warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set."); |
1953 | } |
1954 | max_spacing = calc_grid(stdoutstdout, box, ir->fourier_spacing, |
Value stored to 'max_spacing' is never read | |
1955 | &(ir->nkx), &(ir->nky), &(ir->nkz)); |
1956 | } |
1957 | |
1958 | /* MRS: eventually figure out better logic for initializing the fep |
1959 | values that makes declaring the lambda and declaring the state not |
1960 | potentially conflict if not handled correctly. */ |
1961 | if (ir->efep != efepNO) |
1962 | { |
1963 | state.fep_state = ir->fepvals->init_fep_state; |
1964 | for (i = 0; i < efptNR; i++) |
1965 | { |
1966 | /* init_lambda trumps state definitions*/ |
1967 | if (ir->fepvals->init_lambda >= 0) |
1968 | { |
1969 | state.lambda[i] = ir->fepvals->init_lambda; |
1970 | } |
1971 | else |
1972 | { |
1973 | if (ir->fepvals->all_lambda[i] == NULL((void*)0)) |
1974 | { |
1975 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 1975, "Values of lambda not set for a free energy calculation!"); |
1976 | } |
1977 | else |
1978 | { |
1979 | state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state]; |
1980 | } |
1981 | } |
1982 | } |
1983 | } |
1984 | |
1985 | if (ir->ePull != epullNO) |
1986 | { |
1987 | set_pull_init(ir, sys, state.x, state.box, state.lambda[efptMASS], oenv, opts->pull_start); |
1988 | } |
1989 | |
1990 | if (ir->bRot) |
1991 | { |
1992 | set_reference_positions(ir->rot, state.x, state.box, |
1993 | opt2fn("-ref", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2bSet("-ref", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
1994 | wi); |
1995 | } |
1996 | |
1997 | /* reset_multinr(sys); */ |
1998 | |
1999 | if (EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir->coulombtype ) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD)) |
2000 | { |
2001 | float ratio = pme_load_estimate(sys, ir, state.box); |
2002 | fprintf(stderrstderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio); |
2003 | /* With free energy we might need to do PME both for the A and B state |
2004 | * charges. This will double the cost, but the optimal performance will |
2005 | * then probably be at a slightly larger cut-off and grid spacing. |
2006 | */ |
2007 | if ((ir->efep == efepNO && ratio > 1.0/2.0) || |
2008 | (ir->efep != efepNO && ratio > 2.0/3.0)) |
2009 | { |
2010 | warning_note(wi, |
2011 | "The optimal PME mesh load for parallel simulations is below 0.5\n" |
2012 | "and for highly parallel simulations between 0.25 and 0.33,\n" |
2013 | "for higher performance, increase the cut-off and the PME grid spacing.\n"); |
2014 | if (ir->efep != efepNO) |
2015 | { |
2016 | warning_note(wi, |
2017 | "For free energy simulations, the optimal load limit increases from 0.5 to 0.667\n"); |
2018 | } |
2019 | } |
2020 | } |
2021 | |
2022 | { |
2023 | char warn_buf[STRLEN4096]; |
2024 | double cio = compute_io(ir, sys->natoms, &sys->groups, F_NRE, 1); |
2025 | sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio); |
2026 | if (cio > 2000) |
2027 | { |
2028 | set_warning_line(wi, mdparin, -1); |
2029 | warning_note(wi, warn_buf); |
2030 | } |
2031 | else |
2032 | { |
2033 | printf("%s\n", warn_buf); |
2034 | } |
2035 | } |
2036 | |
2037 | if (bVerbose) |
2038 | { |
2039 | fprintf(stderrstderr, "writing run input file...\n"); |
2040 | } |
2041 | |
2042 | done_warning(wi, FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/grompp.c" , 2042); |
2043 | write_tpx_state(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ir, &state, sys); |
2044 | |
2045 | /* Output IMD group, if bIMD is TRUE */ |
2046 | write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2047 | |
2048 | done_atomtype(atype); |
2049 | done_mtop(sys, TRUE1); |
2050 | done_inputrec_strings(); |
2051 | |
2052 | return 0; |
2053 | } |