Bug Summary

File:gromacs/gmxpreprocess/readir.c
Location:line 3584, column 5
Description:Value stored to 'nQMmult' is never read

Annotated Source Code

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#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <ctype.h>
42#include <stdlib.h>
43#include <limits.h>
44#include "gromacs/utility/smalloc.h"
45#include "typedefs.h"
46#include "physics.h"
47#include "names.h"
48#include "gromacs/utility/fatalerror.h"
49#include "macros.h"
50#include "index.h"
51#include "symtab.h"
52#include "gromacs/utility/cstringutil.h"
53#include "readinp.h"
54#include "warninp.h"
55#include "readir.h"
56#include "toputil.h"
57#include "index.h"
58#include "network.h"
59#include "gromacs/math/vec.h"
60#include "pbc.h"
61#include "mtop_util.h"
62#include "chargegroup.h"
63#include "inputrec.h"
64#include "calc_verletbuf.h"
65
66#define MAXPTR254 254
67#define NOGID255 255
68
69/* Resource parameters
70 * Do not change any of these until you read the instruction
71 * in readinp.h. Some cpp's do not take spaces after the backslash
72 * (like the c-shell), which will give you a very weird compiler
73 * message.
74 */
75
76typedef struct t_inputrec_strings
77{
78 char tcgrps[STRLEN4096], tau_t[STRLEN4096], ref_t[STRLEN4096],
79 acc[STRLEN4096], accgrps[STRLEN4096], freeze[STRLEN4096], frdim[STRLEN4096],
80 energy[STRLEN4096], user1[STRLEN4096], user2[STRLEN4096], vcm[STRLEN4096], x_compressed_groups[STRLEN4096],
81 couple_moltype[STRLEN4096], orirefitgrp[STRLEN4096], egptable[STRLEN4096], egpexcl[STRLEN4096],
82 wall_atomtype[STRLEN4096], wall_density[STRLEN4096], deform[STRLEN4096], QMMM[STRLEN4096],
83 imd_grp[STRLEN4096];
84 char fep_lambda[efptNR][STRLEN4096];
85 char lambda_weights[STRLEN4096];
86 char **pull_grp;
87 char **rot_grp;
88 char anneal[STRLEN4096], anneal_npoints[STRLEN4096],
89 anneal_time[STRLEN4096], anneal_temp[STRLEN4096];
90 char QMmethod[STRLEN4096], QMbasis[STRLEN4096], QMcharge[STRLEN4096], QMmult[STRLEN4096],
91 bSH[STRLEN4096], CASorbitals[STRLEN4096], CASelectrons[STRLEN4096], SAon[STRLEN4096],
92 SAoff[STRLEN4096], SAsteps[STRLEN4096], bTS[STRLEN4096], bOPT[STRLEN4096];
93 char efield_x[STRLEN4096], efield_xt[STRLEN4096], efield_y[STRLEN4096],
94 efield_yt[STRLEN4096], efield_z[STRLEN4096], efield_zt[STRLEN4096];
95
96} gmx_inputrec_strings;
97
98static gmx_inputrec_strings *is = NULL((void*)0);
99
100void init_inputrec_strings()
101{
102 if (is)
103 {
104 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.")_gmx_error("incons", "Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time."
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 104)
;
105 }
106 snew(is, 1)(is) = save_calloc("is", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 106, (1), sizeof(*(is)))
;
107}
108
109void done_inputrec_strings()
110{
111 sfree(is)save_free("is", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 111, (is))
;
112 is = NULL((void*)0);
113}
114
115static char swapgrp[STRLEN4096], splitgrp0[STRLEN4096], splitgrp1[STRLEN4096], solgrp[STRLEN4096];
116
117enum {
118 egrptpALL, /* All particles have to be a member of a group. */
119 egrptpALL_GENREST, /* A rest group with name is generated for particles *
120 * that are not part of any group. */
121 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
122 * for the rest group. */
123 egrptpONE /* Merge all selected groups into one group, *
124 * make a rest group for the remaining particles. */
125};
126
127static const char *constraints[eshNR+1] = {
128 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL((void*)0)
129};
130
131static const char *couple_lam[ecouplamNR+1] = {
132 "vdw-q", "vdw", "q", "none", NULL((void*)0)
133};
134
135void init_ir(t_inputrec *ir, t_gromppopts *opts)
136{
137 snew(opts->include, STRLEN)(opts->include) = save_calloc("opts->include", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 137, (4096), sizeof(*(opts->include)))
;
138 snew(opts->define, STRLEN)(opts->define) = save_calloc("opts->define", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 138, (4096), sizeof(*(opts->define)))
;
139 snew(ir->fepvals, 1)(ir->fepvals) = save_calloc("ir->fepvals", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 139, (1), sizeof(*(ir->fepvals)))
;
140 snew(ir->expandedvals, 1)(ir->expandedvals) = save_calloc("ir->expandedvals", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 140, (1), sizeof(*(ir->expandedvals)))
;
141 snew(ir->simtempvals, 1)(ir->simtempvals) = save_calloc("ir->simtempvals", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 141, (1), sizeof(*(ir->simtempvals)))
;
142}
143
144static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
145{
146
147 int i;
148
149 for (i = 0; i < ntemps; i++)
150 {
151 /* simple linear scaling -- allows more control */
152 if (simtemp->eSimTempScale == esimtempLINEAR)
153 {
154 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
155 }
156 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
157 {
158 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
159 }
160 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
161 {
162 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
163 }
164 else
165 {
166 char errorstr[128];
167 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
168 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 168
, errorstr);
169 }
170 }
171}
172
173
174
175static void _low_check(gmx_bool b, char *s, warninp_t wi)
176{
177 if (b)
178 {
179 warning_error(wi, s);
180 }
181}
182
183static void check_nst(const char *desc_nst, int nst,
184 const char *desc_p, int *p,
185 warninp_t wi)
186{
187 char buf[STRLEN4096];
188
189 if (*p > 0 && *p % nst != 0)
190 {
191 /* Round up to the next multiple of nst */
192 *p = ((*p)/nst + 1)*nst;
193 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
194 desc_p, desc_nst, desc_p, *p);
195 warning(wi, buf);
196 }
197}
198
199static gmx_bool ir_NVE(const t_inputrec *ir)
200{
201 return ((ir->eI == eiMD || EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK)) && ir->etc == etcNO);
202}
203
204static int lcd(int n1, int n2)
205{
206 int d, i;
207
208 d = 1;
209 for (i = 2; (i <= n1 && i <= n2); i++)
210 {
211 if (n1 % i == 0 && n2 % i == 0)
212 {
213 d = i;
214 }
215 }
216
217 return d;
218}
219
220static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
221{
222 if (*eintmod == eintmodPOTSHIFT_VERLET)
223 {
224 if (ir->cutoff_scheme == ecutsVERLET)
225 {
226 *eintmod = eintmodPOTSHIFT;
227 }
228 else
229 {
230 *eintmod = eintmodNONE;
231 }
232 }
233}
234
235void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
236 warninp_t wi)
237/* Check internal consistency */
238{
239 /* Strange macro: first one fills the err_buf, and then one can check
240 * the condition, which will print the message and increase the error
241 * counter.
242 */
243#define CHECK(b)_low_check(b, err_buf, wi) _low_check(b, err_buf, wi)
244 char err_buf[256], warn_buf[STRLEN4096];
245 int i, j;
246 int ns_type = 0;
247 real dt_coupl = 0;
248 real dt_pcoupl;
249 int nstcmin;
250 t_lambda *fep = ir->fepvals;
251 t_expanded *expand = ir->expandedvals;
252
253 set_warning_line(wi, mdparin, -1);
254
255 /* BASIC CUT-OFF STUFF */
256 if (ir->rcoulomb < 0)
257 {
258 warning_error(wi, "rcoulomb should be >= 0");
259 }
260 if (ir->rvdw < 0)
261 {
262 warning_error(wi, "rvdw should be >= 0");
263 }
264 if (ir->rlist < 0 &&
265 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
266 {
267 warning_error(wi, "rlist should be >= 0");
268 }
269
270 process_interaction_modifier(ir, &ir->coulomb_modifier);
271 process_interaction_modifier(ir, &ir->vdw_modifier);
272
273 if (ir->cutoff_scheme == ecutsGROUP)
274 {
275 warning_note(wi,
276 "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
277 "release when all interaction forms are supported for the verlet scheme. The verlet "
278 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
279
280 /* BASIC CUT-OFF STUFF */
281 if (ir->rlist == 0 ||
282 !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
283 (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
284 {
285 /* No switched potential and/or no twin-range:
286 * we can set the long-range cut-off to the maximum of the other cut-offs.
287 */
288 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
289 }
290 else if (ir->rlistlong < 0)
291 {
292 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
293 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
294 ir->rlistlong);
295 warning(wi, warn_buf);
296 }
297 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
298 {
299 warning_error(wi, "Can not have an infinite cut-off with PBC");
300 }
301 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
302 {
303 warning_error(wi, "rlistlong can not be shorter than rlist");
304 }
305 if (IR_TWINRANGE(*ir)((*ir).rlist > 0 && ((*ir).rlistlong == 0 || (*ir)
.rlistlong > (*ir).rlist))
&& ir->nstlist <= 0)
306 {
307 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
308 }
309 }
310
311 if (ir->rlistlong == ir->rlist)
312 {
313 ir->nstcalclr = 0;
314 }
315 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
316 {
317 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
318 }
319
320 if (ir->cutoff_scheme == ecutsVERLET)
321 {
322 real rc_max;
323
324 /* Normal Verlet type neighbor-list, currently only limited feature support */
325 if (inputrec2nboundeddim(ir) < 3)
326 {
327 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
328 }
329 if (ir->rcoulomb != ir->rvdw)
330 {
331 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
332 }
333 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
334 {
335 if (ir->vdw_modifier == eintmodNONE ||
336 ir->vdw_modifier == eintmodPOTSHIFT)
337 {
338 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
339
340 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
341 warning_note(wi, warn_buf);
342
343 ir->vdwtype = evdwCUT;
344 }
345 else
346 {
347 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
348 warning_error(wi, warn_buf);
349 }
350 }
351
352 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
353 {
354 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
355 }
356 if (!(ir->coulombtype == eelCUT ||
357 (EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF
|| (ir->coulombtype) == eelRF_NEC || (ir->coulombtype)
== eelRF_ZERO )
&& ir->coulombtype != eelRF_NEC) ||
358 EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH
|| (ir->coulombtype) == eelPMEUSER || (ir->coulombtype
) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD)
|| ir->coulombtype == eelEWALD))
359 {
360 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
361 }
362 if (!(ir->coulomb_modifier == eintmodNONE ||
363 ir->coulomb_modifier == eintmodPOTSHIFT))
364 {
365 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
366 warning_error(wi, warn_buf);
367 }
368
369 if (ir->nstlist <= 0)
370 {
371 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
372 }
373
374 if (ir->nstlist < 10)
375 {
376 warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
377 }
378
379 rc_max = max(ir->rvdw, ir->rcoulomb)(((ir->rvdw) > (ir->rcoulomb)) ? (ir->rvdw) : (ir
->rcoulomb) )
;
380
381 if (ir->verletbuf_tol <= 0)
382 {
383 if (ir->verletbuf_tol == 0)
384 {
385 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
386 }
387
388 if (ir->rlist < rc_max)
389 {
390 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
391 }
392
393 if (ir->rlist == rc_max && ir->nstlist > 1)
394 {
395 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
396 }
397 }
398 else
399 {
400 if (ir->rlist > rc_max)
401 {
402 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
403 }
404
405 if (ir->nstlist == 1)
406 {
407 /* No buffer required */
408 ir->rlist = rc_max;
409 }
410 else
411 {
412 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD)
)
413 {
414 if (inputrec2nboundeddim(ir) < 3)
415 {
416 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
417 }
418 /* Set rlist temporarily so we can continue processing */
419 ir->rlist = rc_max;
420 }
421 else
422 {
423 /* Set the buffer to 5% of the cut-off */
424 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
425 }
426 }
427 }
428
429 /* No twin-range calculations with Verlet lists */
430 ir->rlistlong = ir->rlist;
431 }
432
433 if (ir->nstcalclr == -1)
434 {
435 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
436 ir->nstcalclr = ir->nstlist;
437 }
438 else if (ir->nstcalclr > 0)
439 {
440 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
441 {
442 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
443 }
444 }
445 else if (ir->nstcalclr < -1)
446 {
447 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
448 }
449
450 if (EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH
|| (ir->coulombtype) == eelPMEUSER || (ir->coulombtype
) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD)
&& ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
451 {
452 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
453 }
454
455 /* GENERAL INTEGRATOR STUFF */
456 if (!(ir->eI == eiMD || EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK)))
457 {
458 ir->etc = etcNO;
459 }
460 if (ir->eI == eiVVAK)
461 {
462 sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
463 warning_note(wi, warn_buf);
464 }
465 if (!EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD)
)
466 {
467 ir->epc = epcNO;
468 }
469 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD)
)
470 {
471 if (ir->nstcalcenergy < 0)
472 {
473 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
474 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
475 {
476 /* nstcalcenergy larger than nstener does not make sense.
477 * We ideally want nstcalcenergy=nstener.
478 */
479 if (ir->nstlist > 0)
480 {
481 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
482 }
483 else
484 {
485 ir->nstcalcenergy = ir->nstenergy;
486 }
487 }
488 }
489 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
490 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
491 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
492
493 {
494 const char *nsten = "nstenergy";
495 const char *nstdh = "nstdhdl";
496 const char *min_name = nsten;
497 int min_nst = ir->nstenergy;
498
499 /* find the smallest of ( nstenergy, nstdhdl ) */
500 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
501 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
502 {
503 min_nst = ir->fepvals->nstdhdl;
504 min_name = nstdh;
505 }
506 /* If the user sets nstenergy small, we should respect that */
507 sprintf(warn_buf,
508 "Setting nstcalcenergy (%d) equal to %s (%d)",
509 ir->nstcalcenergy, min_name, min_nst);
510 warning_note(wi, warn_buf);
511 ir->nstcalcenergy = min_nst;
512 }
513
514 if (ir->epc != epcNO)
515 {
516 if (ir->nstpcouple < 0)
517 {
518 ir->nstpcouple = ir_optimal_nstpcouple(ir);
519 }
520 }
521 if (IR_TWINRANGE(*ir)((*ir).rlist > 0 && ((*ir).rlistlong == 0 || (*ir)
.rlistlong > (*ir).rlist))
)
522 {
523 check_nst("nstlist", ir->nstlist,
524 "nstcalcenergy", &ir->nstcalcenergy, wi);
525 if (ir->epc != epcNO)
526 {
527 check_nst("nstlist", ir->nstlist,
528 "nstpcouple", &ir->nstpcouple, wi);
529 }
530 }
531
532 if (ir->nstcalcenergy > 0)
533 {
534 if (ir->efep != efepNO)
535 {
536 /* nstdhdl should be a multiple of nstcalcenergy */
537 check_nst("nstcalcenergy", ir->nstcalcenergy,
538 "nstdhdl", &ir->fepvals->nstdhdl, wi);
539 /* nstexpanded should be a multiple of nstcalcenergy */
540 check_nst("nstcalcenergy", ir->nstcalcenergy,
541 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
542 }
543 /* for storing exact averages nstenergy should be
544 * a multiple of nstcalcenergy
545 */
546 check_nst("nstcalcenergy", ir->nstcalcenergy,
547 "nstenergy", &ir->nstenergy, wi);
548 }
549 }
550
551 /* LD STUFF */
552 if ((EI_SD(ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || ir->eI == eiBD) &&
553 ir->bContinuation && ir->ld_seed != -1)
554 {
555 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
556 }
557
558 /* TPI STUFF */
559 if (EI_TPI(ir->eI)((ir->eI) == eiTPI || (ir->eI) == eiTPIC))
560 {
561 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
562 CHECK(ir->ePBC != epbcXYZ)_low_check(ir->ePBC != epbcXYZ, err_buf, wi);
563 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
564 CHECK(ir->ns_type != ensGRID)_low_check(ir->ns_type != ensGRID, err_buf, wi);
565 sprintf(err_buf, "with TPI nstlist should be larger than zero");
566 CHECK(ir->nstlist <= 0)_low_check(ir->nstlist <= 0, err_buf, wi);
567 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
568 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype))_low_check(((((ir->coulombtype) == eelPME || (ir->coulombtype
) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir
->coulombtype) == eelPMEUSERSWITCH || (ir->coulombtype)
== eelP3M_AD) || (ir->coulombtype) == eelEWALD) || (ir->
coulombtype) == eelPOISSON) && !((ir->coulombtype)
== eelPME || (ir->coulombtype) == eelPMESWITCH || (ir->
coulombtype) == eelPMEUSER || (ir->coulombtype) == eelPMEUSERSWITCH
|| (ir->coulombtype) == eelP3M_AD), err_buf, wi)
;
569 }
570
571 /* SHAKE / LINCS */
572 if ( (opts->nshake > 0) && (opts->bMorse) )
573 {
574 sprintf(warn_buf,
575 "Using morse bond-potentials while constraining bonds is useless");
576 warning(wi, warn_buf);
577 }
578
579 if ((EI_SD(ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || ir->eI == eiBD) &&
580 ir->bContinuation && ir->ld_seed != -1)
581 {
582 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
583 }
584 /* verify simulated tempering options */
585
586 if (ir->bSimTemp)
587 {
588 gmx_bool bAllTempZero = TRUE1;
589 for (i = 0; i < fep->n_lambda; i++)
590 {
591 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
592 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1))_low_check((fep->all_lambda[efptTEMPERATURE][i] < 0) ||
(fep->all_lambda[efptTEMPERATURE][i] > 1), err_buf, wi
)
;
593 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
594 {
595 bAllTempZero = FALSE0;
596 }
597 }
598 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
599 CHECK(bAllTempZero == TRUE)_low_check(bAllTempZero == 1, err_buf, wi);
600
601 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
602 CHECK(ir->eI != eiVV)_low_check(ir->eI != eiVV, err_buf, wi);
603
604 /* check compatability of the temperature coupling with simulated tempering */
605
606 if (ir->etc == etcNOSEHOOVER)
607 {
608 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
609 warning_note(wi, warn_buf);
610 }
611
612 /* check that the temperatures make sense */
613
614 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
615 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low)_low_check(ir->simtempvals->simtemp_high <= ir->simtempvals
->simtemp_low, err_buf, wi)
;
616
617 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
618 CHECK(ir->simtempvals->simtemp_high <= 0)_low_check(ir->simtempvals->simtemp_high <= 0, err_buf
, wi)
;
619
620 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
621 CHECK(ir->simtempvals->simtemp_low <= 0)_low_check(ir->simtempvals->simtemp_low <= 0, err_buf
, wi)
;
622 }
623
624 /* verify free energy options */
625
626 if (ir->efep != efepNO)
627 {
628 fep = ir->fepvals;
629 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
630 fep->sc_power);
631 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2)_low_check(fep->sc_alpha != 0 && fep->sc_power !=
1 && fep->sc_power != 2, err_buf, wi)
;
632
633 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
634 (int)fep->sc_r_power);
635 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0)_low_check(fep->sc_alpha != 0 && fep->sc_r_power
!= 6.0 && fep->sc_r_power != 48.0, err_buf, wi)
;
636
637 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
638 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)))_low_check(fep->delta_lambda > 0 && ((fep->init_fep_state
> 0) || (fep->init_lambda > 0)), err_buf, wi)
;
639
640 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
641 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED))_low_check(fep->delta_lambda > 0 && (ir->efep
== efepEXPANDED), err_buf, wi)
;
642
643 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
644 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED))_low_check(!(((ir->eI) == eiVV || (ir->eI) == eiVVAK)) &&
(ir->efep == efepEXPANDED), err_buf, wi)
;
645
646 sprintf(err_buf, "Free-energy not implemented for Ewald");
647 CHECK(ir->coulombtype == eelEWALD)_low_check(ir->coulombtype == eelEWALD, err_buf, wi);
648
649 /* check validty of lambda inputs */
650 if (fep->n_lambda == 0)
651 {
652 /* Clear output in case of no states:*/
653 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
654 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0))_low_check((fep->init_fep_state >= 0) && (fep->
n_lambda == 0), err_buf, wi)
;
655 }
656 else
657 {
658 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
659 CHECK((fep->init_fep_state >= fep->n_lambda))_low_check((fep->init_fep_state >= fep->n_lambda), err_buf
, wi)
;
660 }
661
662 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
663 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0))_low_check((fep->init_fep_state < 0) && (fep->
init_lambda < 0), err_buf, wi)
;
664
665 sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
666 fep->init_lambda, fep->init_fep_state);
667 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0))_low_check((fep->init_fep_state >= 0) && (fep->
init_lambda >= 0), err_buf, wi)
;
668
669
670
671 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
672 {
673 int n_lambda_terms;
674 n_lambda_terms = 0;
675 for (i = 0; i < efptNR; i++)
676 {
677 if (fep->separate_dvdl[i])
678 {
679 n_lambda_terms++;
680 }
681 }
682 if (n_lambda_terms > 1)
683 {
684 sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
685 warning(wi, warn_buf);
686 }
687
688 if (n_lambda_terms < 2 && fep->n_lambda > 0)
689 {
690 warning_note(wi,
691 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
692 }
693 }
694
695 for (j = 0; j < efptNR; j++)
696 {
697 for (i = 0; i < fep->n_lambda; i++)
698 {
699 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
700 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1))_low_check((fep->all_lambda[j][i] < 0) || (fep->all_lambda
[j][i] > 1), err_buf, wi)
;
701 }
702 }
703
704 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
705 {
706 for (i = 0; i < fep->n_lambda; i++)
707 {
708 sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
709 fep->all_lambda[efptCOUL][i]);
710 CHECK((fep->sc_alpha > 0) &&_low_check((fep->sc_alpha > 0) && (((fep->all_lambda
[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL
][i] < 1.0)) && ((fep->all_lambda[efptVDW][i] >
0.0) && (fep->all_lambda[efptVDW][i] < 1.0))),
err_buf, wi)
711 (((fep->all_lambda[efptCOUL][i] > 0.0) &&_low_check((fep->sc_alpha > 0) && (((fep->all_lambda
[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL
][i] < 1.0)) && ((fep->all_lambda[efptVDW][i] >
0.0) && (fep->all_lambda[efptVDW][i] < 1.0))),
err_buf, wi)
712 (fep->all_lambda[efptCOUL][i] < 1.0)) &&_low_check((fep->sc_alpha > 0) && (((fep->all_lambda
[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL
][i] < 1.0)) && ((fep->all_lambda[efptVDW][i] >
0.0) && (fep->all_lambda[efptVDW][i] < 1.0))),
err_buf, wi)
713 ((fep->all_lambda[efptVDW][i] > 0.0) &&_low_check((fep->sc_alpha > 0) && (((fep->all_lambda
[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL
][i] < 1.0)) && ((fep->all_lambda[efptVDW][i] >
0.0) && (fep->all_lambda[efptVDW][i] < 1.0))),
err_buf, wi)
714 (fep->all_lambda[efptVDW][i] < 1.0))))_low_check((fep->sc_alpha > 0) && (((fep->all_lambda
[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL
][i] < 1.0)) && ((fep->all_lambda[efptVDW][i] >
0.0) && (fep->all_lambda[efptVDW][i] < 1.0))),
err_buf, wi)
;
715 }
716 }
717
718 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH
|| (ir->coulombtype) == eelPMEUSER || (ir->coulombtype
) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD)
))
719 {
720 real sigma, lambda, r_sc;
721
722 sigma = 0.34;
723 /* Maximum estimate for A and B charges equal with lambda power 1 */
724 lambda = 0.5;
725 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
726 sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
727 fep->sc_r_power,
728 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
729 warning_note(wi, warn_buf);
730 }
731
732 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
733 be treated differently, but that's the next step */
734
735 for (i = 0; i < efptNR; i++)
736 {
737 for (j = 0; j < fep->n_lambda; j++)
738 {
739 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
740 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1))_low_check((fep->all_lambda[i][j] < 0) || (fep->all_lambda
[i][j] > 1), err_buf, wi)
;
741 }
742 }
743 }
744
745 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
746 {
747 fep = ir->fepvals;
748 expand = ir->expandedvals;
749
750 /* checking equilibration of weights inputs for validity */
751
752 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
753 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
754 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM))_low_check((expand->equil_n_at_lam > 0) && (expand
->elmceq != elmceqNUMATLAM), err_buf, wi)
;
755
756 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
757 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
758 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES))_low_check((expand->equil_samples > 0) && (expand
->elmceq != elmceqSAMPLES), err_buf, wi)
;
759
760 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
761 expand->equil_steps, elmceq_names[elmceqSTEPS]);
762 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS))_low_check((expand->equil_steps > 0) && (expand
->elmceq != elmceqSTEPS), err_buf, wi)
;
763
764 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
765 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
766 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA))_low_check((expand->equil_wl_delta > 0) && (expand
->elmceq != elmceqWLDELTA), err_buf, wi)
;
767
768 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
769 expand->equil_ratio, elmceq_names[elmceqRATIO]);
770 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO))_low_check((expand->equil_ratio > 0) && (expand
->elmceq != elmceqRATIO), err_buf, wi)
;
771
772 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
773 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
774 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM))_low_check((expand->equil_n_at_lam <= 0) && (expand
->elmceq == elmceqNUMATLAM), err_buf, wi)
;
775
776 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
777 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
778 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES))_low_check((expand->equil_samples <= 0) && (expand
->elmceq == elmceqSAMPLES), err_buf, wi)
;
779
780 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
781 expand->equil_steps, elmceq_names[elmceqSTEPS]);
782 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS))_low_check((expand->equil_steps <= 0) && (expand
->elmceq == elmceqSTEPS), err_buf, wi)
;
783
784 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
785 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
786 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA))_low_check((expand->equil_wl_delta <= 0) && (expand
->elmceq == elmceqWLDELTA), err_buf, wi)
;
787
788 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
789 expand->equil_ratio, elmceq_names[elmceqRATIO]);
790 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO))_low_check((expand->equil_ratio <= 0) && (expand
->elmceq == elmceqRATIO), err_buf, wi)
;
791
792 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
793 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
794 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)))_low_check((expand->elmceq == elmceqWLDELTA) && (!
((expand->elamstats) == elamstatsWL || (expand->elamstats
) == elamstatsWWL)), err_buf, wi)
;
795
796 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
797 CHECK((expand->lmc_repeats <= 0))_low_check((expand->lmc_repeats <= 0), err_buf, wi);
798 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
799 CHECK((expand->minvarmin <= 0))_low_check((expand->minvarmin <= 0), err_buf, wi);
800 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
801 CHECK((expand->c_range < 0))_low_check((expand->c_range < 0), err_buf, wi);
802 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
803 fep->init_fep_state, expand->lmc_forced_nstart);
804 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO))_low_check((fep->init_fep_state != 0) && (expand->
lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO
), err_buf, wi)
;
805 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
806 CHECK((expand->lmc_forced_nstart < 0))_low_check((expand->lmc_forced_nstart < 0), err_buf, wi
)
;
807 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
808 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda))_low_check((fep->init_fep_state < 0) || (fep->init_fep_state
>= fep->n_lambda), err_buf, wi)
;
809
810 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
811 CHECK((expand->init_wl_delta < 0))_low_check((expand->init_wl_delta < 0), err_buf, wi);
812 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
813 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1))_low_check((expand->wl_ratio <= 0) || (expand->wl_ratio
>= 1), err_buf, wi)
;
814 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
815 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1))_low_check((expand->wl_scale <= 0) || (expand->wl_scale
>= 1), err_buf, wi)
;
816
817 /* if there is no temperature control, we need to specify an MC temperature */
818 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
819 if (expand->nstTij > 0)
820 {
821 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
822 expand->nstTij, ir->nstlog);
823 CHECK((mod(expand->nstTij, ir->nstlog) != 0))_low_check((_mod((expand->nstTij), (ir->nstlog), "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 823) != 0), err_buf, wi)
;
824 }
825 }
826
827 /* PBC/WALLS */
828 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
829 CHECK(ir->nwall && ir->ePBC != epbcXY)_low_check(ir->nwall && ir->ePBC != epbcXY, err_buf
, wi)
;
830
831 /* VACUUM STUFF */
832 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
833 {
834 if (ir->ePBC == epbcNONE)
835 {
836 if (ir->epc != epcNO)
837 {
838 warning(wi, "Turning off pressure coupling for vacuum system");
839 ir->epc = epcNO;
840 }
841 }
842 else
843 {
844 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
845 epbc_names[ir->ePBC]);
846 CHECK(ir->epc != epcNO)_low_check(ir->epc != epcNO, err_buf, wi);
847 }
848 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
849 CHECK(EEL_FULL(ir->coulombtype))_low_check(((((ir->coulombtype) == eelPME || (ir->coulombtype
) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir
->coulombtype) == eelPMEUSERSWITCH || (ir->coulombtype)
== eelP3M_AD) || (ir->coulombtype) == eelEWALD) || (ir->
coulombtype) == eelPOISSON), err_buf, wi)
;
850
851 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
852 epbc_names[ir->ePBC]);
853 CHECK(ir->eDispCorr != edispcNO)_low_check(ir->eDispCorr != edispcNO, err_buf, wi);
854 }
855
856 if (ir->rlist == 0.0)
857 {
858 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
859 "with coulombtype = %s or coulombtype = %s\n"
860 "without periodic boundary conditions (pbc = %s) and\n"
861 "rcoulomb and rvdw set to zero",
862 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
863 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||_low_check(((ir->coulombtype != eelCUT) && (ir->
coulombtype != eelUSER)) || (ir->ePBC != epbcNONE) || (ir->
rcoulomb != 0.0) || (ir->rvdw != 0.0), err_buf, wi)
864 (ir->ePBC != epbcNONE) ||_low_check(((ir->coulombtype != eelCUT) && (ir->
coulombtype != eelUSER)) || (ir->ePBC != epbcNONE) || (ir->
rcoulomb != 0.0) || (ir->rvdw != 0.0), err_buf, wi)
865 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0))_low_check(((ir->coulombtype != eelCUT) && (ir->
coulombtype != eelUSER)) || (ir->ePBC != epbcNONE) || (ir->
rcoulomb != 0.0) || (ir->rvdw != 0.0), err_buf, wi)
;
866
867 if (ir->nstlist < 0)
868 {
869 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
870 }
871 if (ir->nstlist > 0)
872 {
873 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
874 }
875 }
876
877 /* COMM STUFF */
878 if (ir->nstcomm == 0)
879 {
880 ir->comm_mode = ecmNO;
881 }
882 if (ir->comm_mode != ecmNO)
883 {
884 if (ir->nstcomm < 0)
885 {
886 warning(wi, "If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
887 ir->nstcomm = abs(ir->nstcomm);
888 }
889
890 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
891 {
892 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
893 ir->nstcomm = ir->nstcalcenergy;
894 }
895
896 if (ir->comm_mode == ecmANGULAR)
897 {
898 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
899 CHECK(ir->bPeriodicMols)_low_check(ir->bPeriodicMols, err_buf, wi);
900 if (ir->ePBC != epbcNONE)
901 {
902 warning(wi, "Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
903 }
904 }
905 }
906
907 if (EI_STATE_VELOCITY(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2))
&& ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
908 {
909 warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
910 }
911
912 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
913 " algorithm not implemented");
914 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))_low_check(((ir->rcoulomb > ir->rlist) || (ir->rvdw
> ir->rlist)) && (ir->ns_type == ensSIMPLE)
, err_buf, wi)
915 && (ir->ns_type == ensSIMPLE))_low_check(((ir->rcoulomb > ir->rlist) || (ir->rvdw
> ir->rlist)) && (ir->ns_type == ensSIMPLE)
, err_buf, wi)
;
916
917 /* TEMPERATURE COUPLING */
918 if (ir->etc == etcYES)
919 {
920 ir->etc = etcBERENDSEN;
921 warning_note(wi, "Old option for temperature coupling given: "
922 "changing \"yes\" to \"Berendsen\"\n");
923 }
924
925 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
926 {
927 if (ir->opts.nhchainlength < 1)
928 {
929 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
930 ir->opts.nhchainlength = 1;
931 warning(wi, warn_buf);
932 }
933
934 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) && ir->opts.nhchainlength > 1)
935 {
936 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
937 ir->opts.nhchainlength = 1;
938 }
939 }
940 else
941 {
942 ir->opts.nhchainlength = 0;
943 }
944
945 if (ir->eI == eiVVAK)
946 {
947 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
948 ei_names[eiVVAK]);
949 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1))_low_check((ir->nsttcouple != 1) || (ir->nstpcouple != 1
), err_buf, wi)
;
950 }
951
952 if (ETC_ANDERSEN(ir->etc)(((ir->etc) == etcANDERSENMASSIVE) || ((ir->etc) == etcANDERSEN
))
)
953 {
954 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
955 CHECK(!(EI_VV(ir->eI)))_low_check(!(((ir->eI) == eiVV || (ir->eI) == eiVVAK)),
err_buf, wi)
;
956
957 for (i = 0; i < ir->opts.ngtc; i++)
958 {
959 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
960 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i])_low_check(ir->opts.tau_t[0] != ir->opts.tau_t[i], err_buf
, wi)
;
961 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
962 i, ir->opts.tau_t[i]);
963 CHECK(ir->opts.tau_t[i] < 0)_low_check(ir->opts.tau_t[i] < 0, err_buf, wi);
964 }
965 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
966 {
967 sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
968 warning_note(wi, warn_buf);
969 }
970
971 sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
972 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN))_low_check(ir->nstcomm > 1 && (ir->etc == etcANDERSEN
), err_buf, wi)
;
973
974 for (i = 0; i < ir->opts.ngtc; i++)
975 {
976 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
977 sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
978 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE))_low_check((nsteps % ir->nstcomm) && (ir->etc ==
etcANDERSENMASSIVE), err_buf, wi)
;
979 }
980 }
981 if (ir->etc == etcBERENDSEN)
982 {
983 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
984 ETCOUPLTYPE(ir->etc)((((ir->etc) < 0) || ((ir->etc) >= (etcNR))) ? "UNDEFINED"
: (etcoupl_names)[ir->etc])
, ETCOUPLTYPE(etcVRESCALE)((((etcVRESCALE) < 0) || ((etcVRESCALE) >= (etcNR))) ? "UNDEFINED"
: (etcoupl_names)[etcVRESCALE])
);
985 warning_note(wi, warn_buf);
986 }
987
988 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)(((ir->etc) == etcANDERSENMASSIVE) || ((ir->etc) == etcANDERSEN
))
)
989 && ir->epc == epcBERENDSEN)
990 {
991 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
992 "true ensemble for the thermostat");
993 warning(wi, warn_buf);
994 }
995
996 /* PRESSURE COUPLING */
997 if (ir->epc == epcISOTROPIC)
998 {
999 ir->epc = epcBERENDSEN;
1000 warning_note(wi, "Old option for pressure coupling given: "
1001 "changing \"Isotropic\" to \"Berendsen\"\n");
1002 }
1003
1004 if (ir->epc != epcNO)
1005 {
1006 dt_pcoupl = ir->nstpcouple*ir->delta_t;
1007
1008 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1009 CHECK(ir->tau_p <= 0)_low_check(ir->tau_p <= 0, err_buf, wi);
1010
1011 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
1012 {
1013 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
1014 EPCOUPLTYPE(ir->epc)((((ir->epc) < 0) || ((ir->epc) >= (epcNR))) ? "UNDEFINED"
: (epcoupl_names)[ir->epc])
, ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1015 warning(wi, warn_buf);
1016 }
1017
1018 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1019 " coupling %s\n", EPCOUPLTYPE(ir->epc)((((ir->epc) < 0) || ((ir->epc) >= (epcNR))) ? "UNDEFINED"
: (epcoupl_names)[ir->epc])
);
1020 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||_low_check(ir->compress[0][0] < 0 || ir->compress[1]
[1] < 0 || ir->compress[2][2] < 0 || (trace(ir->compress
) == 0 && ir->compress[1][0] <= 0 && ir
->compress[2][0] <= 0 && ir->compress[2][1] <=
0), err_buf, wi)
1021 ir->compress[ZZ][ZZ] < 0 ||_low_check(ir->compress[0][0] < 0 || ir->compress[1]
[1] < 0 || ir->compress[2][2] < 0 || (trace(ir->compress
) == 0 && ir->compress[1][0] <= 0 && ir
->compress[2][0] <= 0 && ir->compress[2][1] <=
0), err_buf, wi)
1022 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&_low_check(ir->compress[0][0] < 0 || ir->compress[1]
[1] < 0 || ir->compress[2][2] < 0 || (trace(ir->compress
) == 0 && ir->compress[1][0] <= 0 && ir
->compress[2][0] <= 0 && ir->compress[2][1] <=
0), err_buf, wi)
1023 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0))_low_check(ir->compress[0][0] < 0 || ir->compress[1]
[1] < 0 || ir->compress[2][2] < 0 || (trace(ir->compress
) == 0 && ir->compress[1][0] <= 0 && ir
->compress[2][0] <= 0 && ir->compress[2][1] <=
0), err_buf, wi)
;
1024
1025 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1026 {
1027 sprintf(warn_buf,
1028 "You are generating velocities so I am assuming you "
1029 "are equilibrating a system. You are using "
1030 "%s pressure coupling, but this can be "
1031 "unstable for equilibration. If your system crashes, try "
1032 "equilibrating first with Berendsen pressure coupling. If "
1033 "you are not equilibrating the system, you can probably "
1034 "ignore this warning.",
1035 epcoupl_names[ir->epc]);
1036 warning(wi, warn_buf);
1037 }
1038 }
1039
1040 if (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK))
1041 {
1042 if (ir->epc > epcNO)
1043 {
1044 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1045 {
1046 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
1047 }
1048 }
1049 }
1050 else
1051 {
1052 if (ir->epc == epcMTTK)
1053 {
1054 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1055 }
1056 }
1057
1058 /* ELECTROSTATICS */
1059 /* More checks are in triple check (grompp.c) */
1060
1061 if (ir->coulombtype == eelSWITCH)
1062 {
1063 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1064 "artifacts, advice: use coulombtype = %s",
1065 eel_names[ir->coulombtype],
1066 eel_names[eelRF_ZERO]);
1067 warning(wi, warn_buf);
1068 }
1069
1070 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1071 {
1072 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1073 warning_note(wi, warn_buf);
1074 }
1075
1076 if (EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF
|| (ir->coulombtype) == eelRF_NEC || (ir->coulombtype)
== eelRF_ZERO )
&& ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1077 {
1078 sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
1079 warning(wi, warn_buf);
1080 ir->epsilon_rf = ir->epsilon_r;
1081 ir->epsilon_r = 1.0;
1082 }
1083
1084 if (getenv("GALACTIC_DYNAMICS") == NULL((void*)0))
1085 {
1086 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1087 CHECK(ir->epsilon_r < 0)_low_check(ir->epsilon_r < 0, err_buf, wi);
1088 }
1089
1090 if (EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF
|| (ir->coulombtype) == eelRF_NEC || (ir->coulombtype)
== eelRF_ZERO )
)
1091 {
1092 /* reaction field (at the cut-off) */
1093
1094 if (ir->coulombtype == eelRF_ZERO)
1095 {
1096 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1097 eel_names[ir->coulombtype]);
1098 CHECK(ir->epsilon_rf != 0)_low_check(ir->epsilon_rf != 0, err_buf, wi);
1099 ir->epsilon_rf = 0.0;
1100 }
1101
1102 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1103 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||_low_check((ir->epsilon_rf < ir->epsilon_r &&
ir->epsilon_rf != 0) || (ir->epsilon_r == 0), err_buf,
wi)
1104 (ir->epsilon_r == 0))_low_check((ir->epsilon_rf < ir->epsilon_r &&
ir->epsilon_rf != 0) || (ir->epsilon_r == 0), err_buf,
wi)
;
1105 if (ir->epsilon_rf == ir->epsilon_r)
1106 {
1107 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1108 eel_names[ir->coulombtype]);
1109 warning(wi, warn_buf);
1110 }
1111 }
1112 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1113 * means the interaction is zero outside rcoulomb, but it helps to
1114 * provide accurate energy conservation.
1115 */
1116 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1117 {
1118 if (ir_coulomb_switched(ir))
1119 {
1120 sprintf(err_buf,
1121 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1122 eel_names[ir->coulombtype]);
1123 CHECK(ir->rcoulomb_switch >= ir->rcoulomb)_low_check(ir->rcoulomb_switch >= ir->rcoulomb, err_buf
, wi)
;
1124 }
1125 }
1126 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF
|| (ir->coulombtype) == eelRF_NEC || (ir->coulombtype)
== eelRF_ZERO )
)
1127 {
1128 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1129 {
1130 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1131 eel_names[ir->coulombtype]);
1132 CHECK(ir->rlist > ir->rcoulomb)_low_check(ir->rlist > ir->rcoulomb, err_buf, wi);
1133 }
1134 }
1135
1136 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1137 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1138 {
1139 sprintf(warn_buf,
1140 "The switch/shift interaction settings are just for compatibility; you will get better "
1141 "performance from applying potential modifiers to your interactions!\n");
1142 warning_note(wi, warn_buf);
1143 }
1144
1145 if (ir->coulombtype == eelPMESWITCH)
1146 {
1147 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1148 {
1149 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1150 eel_names[ir->coulombtype],
1151 ir->ewald_rtol);
1152 warning(wi, warn_buf);
1153 }
1154 }
1155
1156 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)
)
1157 {
1158 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1159 ir->coulombtype == eelPMEUSERSWITCH)
1160 {
1161 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1162 eel_names[ir->coulombtype]);
1163 CHECK(ir->rcoulomb > ir->rlist)_low_check(ir->rcoulomb > ir->rlist, err_buf, wi);
1164 }
1165 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1166 {
1167 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1168 {
1169 sprintf(err_buf,
1170 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1171 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1172 "a potential modifier.", eel_names[ir->coulombtype]);
1173 if (ir->nstcalclr == 1)
1174 {
1175 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong)_low_check(ir->rcoulomb != ir->rlist && ir->
rcoulomb != ir->rlistlong, err_buf, wi)
;
1176 }
1177 else
1178 {
1179 CHECK(ir->rcoulomb != ir->rlist)_low_check(ir->rcoulomb != ir->rlist, err_buf, wi);
1180 }
1181 }
1182 }
1183 }
1184
1185 if (EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH
|| (ir->coulombtype) == eelPMEUSER || (ir->coulombtype
) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD)
|| EVDW_PME(ir->vdwtype)((ir->vdwtype) == evdwPME))
1186 {
1187 if (ir->pme_order < 3)
1188 {
1189 warning_error(wi, "pme-order can not be smaller than 3");
1190 }
1191 }
1192
1193 if (ir->nwall == 2 && 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)
)
1194 {
1195 if (ir->ewald_geometry == eewg3D)
1196 {
1197 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1198 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1199 warning(wi, warn_buf);
1200 }
1201 /* This check avoids extra pbc coding for exclusion corrections */
1202 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1203 CHECK(ir->wall_ewald_zfac < 2)_low_check(ir->wall_ewald_zfac < 2, err_buf, wi);
1204 }
1205
1206 if (ir_vdw_switched(ir))
1207 {
1208 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1209 CHECK(ir->rvdw_switch >= ir->rvdw)_low_check(ir->rvdw_switch >= ir->rvdw, err_buf, wi);
1210
1211 if (ir->rvdw_switch < 0.5*ir->rvdw)
1212 {
1213 sprintf(warn_buf, "You are applying a switch function to vdw forces or potentials from %g to %g nm, which is more than half the interaction range, whereas switch functions are intended to act only close to the cut-off.",
1214 ir->rvdw_switch, ir->rvdw);
1215 warning_note(wi, warn_buf);
1216 }
1217 }
1218 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1219 {
1220 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1221 {
1222 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1223 CHECK(ir->rlist > ir->rvdw)_low_check(ir->rlist > ir->rvdw, err_buf, wi);
1224 }
1225 }
1226
1227 if (ir->vdwtype == evdwPME)
1228 {
1229 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1230 {
1231 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
1232nd %s",
1233 evdw_names[ir->vdwtype],
1234 eintmod_names[eintmodPOTSHIFT],
1235 eintmod_names[eintmodNONE]);
1236 }
1237 }
1238
1239 if (ir->cutoff_scheme == ecutsGROUP)
1240 {
1241 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1242 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1243 ir->nstlist != 1)
1244 {
1245 warning_note(wi, "With exact cut-offs, rlist should be "
1246 "larger than rcoulomb and rvdw, so that there "
1247 "is a buffer region for particle motion "
1248 "between neighborsearch steps");
1249 }
1250
1251 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
1252 {
1253 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1254 IR_TWINRANGE(*ir)((*ir).rlist > 0 && ((*ir).rlistlong == 0 || (*ir)
.rlistlong > (*ir).rlist))
? "rlistlong" : "rlist");
1255 warning_note(wi, warn_buf);
1256 }
1257 if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
1258 {
1259 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1260 IR_TWINRANGE(*ir)((*ir).rlist > 0 && ((*ir).rlistlong == 0 || (*ir)
.rlistlong > (*ir).rlist))
? "rlistlong" : "rlist");
1261 warning_note(wi, warn_buf);
1262 }
1263 }
1264
1265 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1266 {
1267 warning_note(wi, "You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
1268 }
1269
1270 if (ir->nstlist == -1)
1271 {
1272 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1273 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist)_low_check(ir->rvdw >= ir->rlist || ir->rcoulomb >=
ir->rlist, err_buf, wi)
;
1274 }
1275 sprintf(err_buf, "nstlist can not be smaller than -1");
1276 CHECK(ir->nstlist < -1)_low_check(ir->nstlist < -1, err_buf, wi);
1277
1278 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1279 && ir->rvdw != 0)
1280 {
1281 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1282 }
1283
1284 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1285 {
1286 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1287 }
1288
1289 /* ENERGY CONSERVATION */
1290 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1291 {
1292 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1293 {
1294 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1295 evdw_names[evdwSHIFT]);
1296 warning_note(wi, warn_buf);
1297 }
1298 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1299 {
1300 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1301 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1302 warning_note(wi, warn_buf);
1303 }
1304 }
1305
1306 /* IMPLICIT SOLVENT */
1307 if (ir->coulombtype == eelGB_NOTUSED)
1308 {
1309 ir->coulombtype = eelCUT;
1310 ir->implicit_solvent = eisGBSA;
1311 fprintf(stderrstderr, "Note: Old option for generalized born electrostatics given:\n"
1312 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1313 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1314 }
1315
1316 if (ir->sa_algorithm == esaSTILL)
1317 {
1318 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1319 CHECK(ir->sa_algorithm == esaSTILL)_low_check(ir->sa_algorithm == esaSTILL, err_buf, wi);
1320 }
1321
1322 if (ir->implicit_solvent == eisGBSA)
1323 {
1324 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1325 CHECK(ir->rgbradii != ir->rlist)_low_check(ir->rgbradii != ir->rlist, err_buf, wi);
1326
1327 if (ir->coulombtype != eelCUT)
1328 {
1329 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1330 CHECK(ir->coulombtype != eelCUT)_low_check(ir->coulombtype != eelCUT, err_buf, wi);
1331 }
1332 if (ir->vdwtype != evdwCUT)
1333 {
1334 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1335 CHECK(ir->vdwtype != evdwCUT)_low_check(ir->vdwtype != evdwCUT, err_buf, wi);
1336 }
1337 if (ir->nstgbradii < 1)
1338 {
1339 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1340 warning_note(wi, warn_buf);
1341 ir->nstgbradii = 1;
1342 }
1343 if (ir->sa_algorithm == esaNO)
1344 {
1345 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1346 warning_note(wi, warn_buf);
1347 }
1348 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1349 {
1350 sprintf(warn_buf, "Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
1351 warning_note(wi, warn_buf);
1352
1353 if (ir->gb_algorithm == egbSTILL)
1354 {
1355 ir->sa_surface_tension = 0.0049 * CAL2JOULE(4.184) * 100;
1356 }
1357 else
1358 {
1359 ir->sa_surface_tension = 0.0054 * CAL2JOULE(4.184) * 100;
1360 }
1361 }
1362 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1363 {
1364 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1365 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)_low_check(ir->sa_surface_tension == 0 && ir->sa_algorithm
!= esaNO, err_buf, wi)
;
1366 }
1367
1368 }
1369
1370 if (ir->bAdress)
1371 {
1372 if (ir->cutoff_scheme != ecutsGROUP)
1373 {
1374 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1375 }
1376 if (!EI_SD(ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2))
1377 {
1378 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1379 }
1380 if (ir->epc != epcNO)
1381 {
1382 warning_error(wi, "AdresS simulation does not support pressure coupling");
1383 }
1384 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)
)
1385 {
1386 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1387 }
1388 }
1389}
1390
1391/* count the number of text elemets separated by whitespace in a string.
1392 str = the input string
1393 maxptr = the maximum number of allowed elements
1394 ptr = the output array of pointers to the first character of each element
1395 returns: the number of elements. */
1396int str_nelem(const char *str, int maxptr, char *ptr[])
1397{
1398 int np = 0;
1399 char *copy0, *copy;
1400
1401 copy0 = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t
)(const void *)((str) + 1) - (size_t)(const void *)(str) == 1
) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, str, __len); __retval
; })) : __strdup (str)))
;
1402 copy = copy0;
1403 ltrim(copy);
1404 while (*copy != '\0')
1405 {
1406 if (np >= maxptr)
1407 {
1408 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1408
, "Too many groups on line: '%s' (max is %d)",
1409 str, maxptr);
1410 }
1411 if (ptr)
1412 {
1413 ptr[np] = copy;
1414 }
1415 np++;
1416 while ((*copy != '\0') && !isspace(*copy)((*__ctype_b_loc ())[(int) ((*copy))] & (unsigned short int
) _ISspace)
)
1417 {
1418 copy++;
1419 }
1420 if (*copy != '\0')
1421 {
1422 *copy = '\0';
1423 copy++;
1424 }
1425 ltrim(copy);
1426 }
1427 if (ptr == NULL((void*)0))
1428 {
1429 sfree(copy0)save_free("copy0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1429, (copy0))
;
1430 }
1431
1432 return np;
1433}
1434
1435/* interpret a number of doubles from a string and put them in an array,
1436 after allocating space for them.
1437 str = the input string
1438 n = the (pre-allocated) number of doubles read
1439 r = the output array of doubles. */
1440static void parse_n_real(char *str, int *n, real **r)
1441{
1442 char *ptr[MAXPTR254];
1443 int i;
1444
1445 *n = str_nelem(str, MAXPTR254, ptr);
1446
1447 snew(*r, *n)(*r) = save_calloc("*r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1447, (*n), sizeof(*(*r)))
;
1448 for (i = 0; i < *n; i++)
1449 {
1450 (*r)[i] = strtod(ptr[i], NULL((void*)0));
1451 }
1452}
1453
1454static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN4096], char weights[STRLEN4096])
1455{
1456
1457 int i, j, max_n_lambda, nweights, nfep[efptNR];
1458 t_lambda *fep = ir->fepvals;
1459 t_expanded *expand = ir->expandedvals;
1460 real **count_fep_lambdas;
1461 gmx_bool bOneLambda = TRUE1;
1462
1463 snew(count_fep_lambdas, efptNR)(count_fep_lambdas) = save_calloc("count_fep_lambdas", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1463, (efptNR), sizeof(*(count_fep_lambdas)))
;
1464
1465 /* FEP input processing */
1466 /* first, identify the number of lambda values for each type.
1467 All that are nonzero must have the same number */
1468
1469 for (i = 0; i < efptNR; i++)
1470 {
1471 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1472 }
1473
1474 /* now, determine the number of components. All must be either zero, or equal. */
1475
1476 max_n_lambda = 0;
1477 for (i = 0; i < efptNR; i++)
1478 {
1479 if (nfep[i] > max_n_lambda)
1480 {
1481 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1482 must have the same number if its not zero.*/
1483 break;
1484 }
1485 }
1486
1487 for (i = 0; i < efptNR; i++)
1488 {
1489 if (nfep[i] == 0)
1490 {
1491 ir->fepvals->separate_dvdl[i] = FALSE0;
1492 }
1493 else if (nfep[i] == max_n_lambda)
1494 {
1495 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1496 respect to the temperature currently */
1497 {
1498 ir->fepvals->separate_dvdl[i] = TRUE1;
1499 }
1500 }
1501 else
1502 {
1503 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1503
, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1504 nfep[i], efpt_names[i], max_n_lambda);
1505 }
1506 }
1507 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1508 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE0;
1509
1510 /* the number of lambdas is the number we've read in, which is either zero
1511 or the same for all */
1512 fep->n_lambda = max_n_lambda;
1513
1514 /* allocate space for the array of lambda values */
1515 snew(fep->all_lambda, efptNR)(fep->all_lambda) = save_calloc("fep->all_lambda", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1515, (efptNR), sizeof(*(fep->all_lambda)))
;
1516 /* if init_lambda is defined, we need to set lambda */
1517 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1518 {
1519 ir->fepvals->separate_dvdl[efptFEP] = TRUE1;
1520 }
1521 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1522 for (i = 0; i < efptNR; i++)
1523 {
1524 snew(fep->all_lambda[i], fep->n_lambda)(fep->all_lambda[i]) = save_calloc("fep->all_lambda[i]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1524, (fep->n_lambda), sizeof(*(fep->all_lambda[i])))
;
1525 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1526 are zero */
1527 {
1528 for (j = 0; j < fep->n_lambda; j++)
1529 {
1530 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1531 }
1532 sfree(count_fep_lambdas[i])save_free("count_fep_lambdas[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1532, (count_fep_lambdas[i]))
;
1533 }
1534 }
1535 sfree(count_fep_lambdas)save_free("count_fep_lambdas", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1535, (count_fep_lambdas))
;
1536
1537 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1538 bookkeeping -- for now, init_lambda */
1539
1540 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1541 {
1542 for (i = 0; i < fep->n_lambda; i++)
1543 {
1544 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1545 }
1546 }
1547
1548 /* check to see if only a single component lambda is defined, and soft core is defined.
1549 In this case, turn on coulomb soft core */
1550
1551 if (max_n_lambda == 0)
1552 {
1553 bOneLambda = TRUE1;
1554 }
1555 else
1556 {
1557 for (i = 0; i < efptNR; i++)
1558 {
1559 if ((nfep[i] != 0) && (i != efptFEP))
1560 {
1561 bOneLambda = FALSE0;
1562 }
1563 }
1564 }
1565 if ((bOneLambda) && (fep->sc_alpha > 0))
1566 {
1567 fep->bScCoul = TRUE1;
1568 }
1569
1570 /* Fill in the others with the efptFEP if they are not explicitly
1571 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1572 they are all zero. */
1573
1574 for (i = 0; i < efptNR; i++)
1575 {
1576 if ((nfep[i] == 0) && (i != efptFEP))
1577 {
1578 for (j = 0; j < fep->n_lambda; j++)
1579 {
1580 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1581 }
1582 }
1583 }
1584
1585
1586 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1587 if (fep->sc_r_power == 48)
1588 {
1589 if (fep->sc_alpha > 0.1)
1590 {
1591 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1591
, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1592 }
1593 }
1594
1595 expand = ir->expandedvals;
1596 /* now read in the weights */
1597 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1598 if (nweights == 0)
1599 {
1600 snew(expand->init_lambda_weights, fep->n_lambda)(expand->init_lambda_weights) = save_calloc("expand->init_lambda_weights"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1600, (fep->n_lambda), sizeof(*(expand->init_lambda_weights
)))
; /* initialize to zero */
1601 }
1602 else if (nweights != fep->n_lambda)
1603 {
1604 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1604
, "Number of weights (%d) is not equal to number of lambda values (%d)",
1605 nweights, fep->n_lambda);
1606 }
1607 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1608 {
1609 expand->nstexpanded = fep->nstdhdl;
1610 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1611 }
1612 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1613 {
1614 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1615 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1616 2*tau_t just to be careful so it's not to frequent */
1617 }
1618}
1619
1620
1621static void do_simtemp_params(t_inputrec *ir)
1622{
1623
1624 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda)(ir->simtempvals->temperatures) = save_calloc("ir->simtempvals->temperatures"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1624, (ir->fepvals->n_lambda), sizeof(*(ir->simtempvals
->temperatures)))
;
1625 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1626
1627 return;
1628}
1629
1630static void do_wall_params(t_inputrec *ir,
1631 char *wall_atomtype, char *wall_density,
1632 t_gromppopts *opts)
1633{
1634 int nstr, i;
1635 char *names[MAXPTR254];
1636 double dbl;
1637
1638 opts->wall_atomtype[0] = NULL((void*)0);
1639 opts->wall_atomtype[1] = NULL((void*)0);
1640
1641 ir->wall_atomtype[0] = -1;
1642 ir->wall_atomtype[1] = -1;
1643 ir->wall_density[0] = 0;
1644 ir->wall_density[1] = 0;
1645
1646 if (ir->nwall > 0)
1647 {
1648 nstr = str_nelem(wall_atomtype, MAXPTR254, names);
1649 if (nstr != ir->nwall)
1650 {
1651 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1651
, "Expected %d elements for wall_atomtype, found %d",
1652 ir->nwall, nstr);
1653 }
1654 for (i = 0; i < ir->nwall; i++)
1655 {
1656 opts->wall_atomtype[i] = strdup(names[i])(__extension__ (__builtin_constant_p (names[i]) && ((
size_t)(const void *)((names[i]) + 1) - (size_t)(const void *
)(names[i]) == 1) ? (((const char *) (names[i]))[0] == '\0' ?
(char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len =
strlen (names[i]) + 1; char *__retval = (char *) malloc (__len
); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval
, names[i], __len); __retval; })) : __strdup (names[i])))
;
1657 }
1658
1659 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1660 {
1661 nstr = str_nelem(wall_density, MAXPTR254, names);
1662 if (nstr != ir->nwall)
1663 {
1664 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1664
, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1665 }
1666 for (i = 0; i < ir->nwall; i++)
1667 {
1668 sscanf(names[i], "%lf", &dbl);
1669 if (dbl <= 0)
1670 {
1671 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1671
, "wall-density[%d] = %f\n", i, dbl);
1672 }
1673 ir->wall_density[i] = dbl;
1674 }
1675 }
1676 }
1677}
1678
1679static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1680{
1681 int i;
1682 t_grps *grps;
1683 char str[STRLEN4096];
1684
1685 if (nwall > 0)
1686 {
1687 srenew(groups->grpname, groups->ngrpname+nwall)(groups->grpname) = save_realloc("groups->grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1687, (groups->grpname), (groups->ngrpname+nwall), sizeof
(*(groups->grpname)))
;
1688 grps = &(groups->grps[egcENER]);
1689 srenew(grps->nm_ind, grps->nr+nwall)(grps->nm_ind) = save_realloc("grps->nm_ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1689, (grps->nm_ind), (grps->nr+nwall), sizeof(*(grps
->nm_ind)))
;
1690 for (i = 0; i < nwall; i++)
1691 {
1692 sprintf(str, "wall%d", i);
1693 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1694 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1695 }
1696 }
1697}
1698
1699void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1700 t_expanded *expand, warninp_t wi)
1701{
1702 int ninp, nerror = 0;
1703 t_inpfile *inp;
1704
1705 ninp = *ninp_p;
1706 inp = *inp_p;
1707
1708 /* read expanded ensemble parameters */
1709 CCTYPE ("expanded ensemble variables")get_estr(&ninp, &inp, "\n; " "expanded ensemble variables"
, ((void*)0))
;
1710 ITYPE ("nstexpanded", expand->nstexpanded, -1)expand->nstexpanded = get_eint(&ninp, &inp, "nstexpanded"
, -1, wi)
;
1711 EETYPE("lmc-stats", expand->elamstats, elamstats_names)expand->elamstats = get_eeenum(&ninp, &inp, "lmc-stats"
, elamstats_names, wi)
;
1712 EETYPE("lmc-move", expand->elmcmove, elmcmove_names)expand->elmcmove = get_eeenum(&ninp, &inp, "lmc-move"
, elmcmove_names, wi)
;
1713 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names)expand->elmceq = get_eeenum(&ninp, &inp, "lmc-weights-equil"
, elmceq_names, wi)
;
1714 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1)expand->equil_n_at_lam = get_eint(&ninp, &inp, "weight-equil-number-all-lambda"
, -1, wi)
;
1715 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1)expand->equil_samples = get_eint(&ninp, &inp, "weight-equil-number-samples"
, -1, wi)
;
1716 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1)expand->equil_steps = get_eint(&ninp, &inp, "weight-equil-number-steps"
, -1, wi)
;
1717 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1)expand->equil_wl_delta = get_ereal(&ninp, &inp, "weight-equil-wl-delta"
, -1, wi)
;
1718 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1)expand->equil_ratio = get_ereal(&ninp, &inp, "weight-equil-count-ratio"
, -1, wi)
;
1719 CCTYPE("Seed for Monte Carlo in lambda space")get_estr(&ninp, &inp, "\n; " "Seed for Monte Carlo in lambda space"
, ((void*)0))
;
1720 ITYPE ("lmc-seed", expand->lmc_seed, -1)expand->lmc_seed = get_eint(&ninp, &inp, "lmc-seed"
, -1, wi)
;
1721 RTYPE ("mc-temperature", expand->mc_temp, -1)expand->mc_temp = get_ereal(&ninp, &inp, "mc-temperature"
, -1, wi)
;
1722 ITYPE ("lmc-repeats", expand->lmc_repeats, 1)expand->lmc_repeats = get_eint(&ninp, &inp, "lmc-repeats"
, 1, wi)
;
1723 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1)expand->gibbsdeltalam = get_eint(&ninp, &inp, "lmc-gibbsdelta"
, -1, wi)
;
1724 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0)expand->lmc_forced_nstart = get_eint(&ninp, &inp, "lmc-forced-nstart"
, 0, wi)
;
1725 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names)expand->bSymmetrizedTMatrix = get_eeenum(&ninp, &inp
, "symmetrized-transition-matrix", yesno_names, wi)
;
1726 ITYPE("nst-transition-matrix", expand->nstTij, -1)expand->nstTij = get_eint(&ninp, &inp, "nst-transition-matrix"
, -1, wi)
;
1727 ITYPE ("mininum-var-min", expand->minvarmin, 100)expand->minvarmin = get_eint(&ninp, &inp, "mininum-var-min"
, 100, wi)
; /*default is reasonable */
1728 ITYPE ("weight-c-range", expand->c_range, 0)expand->c_range = get_eint(&ninp, &inp, "weight-c-range"
, 0, wi)
; /* default is just C=0 */
1729 RTYPE ("wl-scale", expand->wl_scale, 0.8)expand->wl_scale = get_ereal(&ninp, &inp, "wl-scale"
, 0.8, wi)
;
1730 RTYPE ("wl-ratio", expand->wl_ratio, 0.8)expand->wl_ratio = get_ereal(&ninp, &inp, "wl-ratio"
, 0.8, wi)
;
1731 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0)expand->init_wl_delta = get_ereal(&ninp, &inp, "init-wl-delta"
, 1.0, wi)
;
1732 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names)expand->bWLoneovert = get_eeenum(&ninp, &inp, "wl-oneovert"
, yesno_names, wi)
;
1733
1734 *ninp_p = ninp;
1735 *inp_p = inp;
1736
1737 return;
1738}
1739
1740void get_ir(const char *mdparin, const char *mdparout,
1741 t_inputrec *ir, t_gromppopts *opts,
1742 warninp_t wi)
1743{
1744 char *dumstr[2];
1745 double dumdub[2][6];
1746 t_inpfile *inp;
1747 const char *tmp;
1748 int i, j, m, ninp;
1749 char warn_buf[STRLEN4096];
1750 t_lambda *fep = ir->fepvals;
1751 t_expanded *expand = ir->expandedvals;
1752
1753 init_inputrec_strings();
1754 inp = read_inpfile(mdparin, &ninp, wi);
1755
1756 snew(dumstr[0], STRLEN)(dumstr[0]) = save_calloc("dumstr[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1756, (4096), sizeof(*(dumstr[0])))
;
1757 snew(dumstr[1], STRLEN)(dumstr[1]) = save_calloc("dumstr[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 1757, (4096), sizeof(*(dumstr[1])))
;
1758
1759 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1760 {
1761 sprintf(warn_buf,
1762 "%s did not specify a value for the .mdp option "
1763 "\"cutoff-scheme\". Probably it was first intended for use "
1764 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1765 "introduced, but the group scheme was still the default. "
1766 "The default is now the Verlet scheme, so you will observe "
1767 "different behaviour.", mdparin);
1768 warning_note(wi, warn_buf);
1769 }
1770
1771 /* remove the following deprecated commands */
1772 REM_TYPE("title")replace_inp_entry(ninp, inp, "title", ((void*)0));
1773 REM_TYPE("cpp")replace_inp_entry(ninp, inp, "cpp", ((void*)0));
1774 REM_TYPE("domain-decomposition")replace_inp_entry(ninp, inp, "domain-decomposition", ((void*)
0))
;
1775 REM_TYPE("andersen-seed")replace_inp_entry(ninp, inp, "andersen-seed", ((void*)0));
1776 REM_TYPE("dihre")replace_inp_entry(ninp, inp, "dihre", ((void*)0));
1777 REM_TYPE("dihre-fc")replace_inp_entry(ninp, inp, "dihre-fc", ((void*)0));
1778 REM_TYPE("dihre-tau")replace_inp_entry(ninp, inp, "dihre-tau", ((void*)0));
1779 REM_TYPE("nstdihreout")replace_inp_entry(ninp, inp, "nstdihreout", ((void*)0));
1780 REM_TYPE("nstcheckpoint")replace_inp_entry(ninp, inp, "nstcheckpoint", ((void*)0));
1781
1782 /* replace the following commands with the clearer new versions*/
1783 REPL_TYPE("unconstrained-start", "continuation")replace_inp_entry(ninp, inp, "unconstrained-start", "continuation"
)
;
1784 REPL_TYPE("foreign-lambda", "fep-lambdas")replace_inp_entry(ninp, inp, "foreign-lambda", "fep-lambdas");
1785 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance")replace_inp_entry(ninp, inp, "verlet-buffer-drift", "verlet-buffer-tolerance"
)
;
1786 REPL_TYPE("nstxtcout", "nstxout-compressed")replace_inp_entry(ninp, inp, "nstxtcout", "nstxout-compressed"
)
;
1787 REPL_TYPE("xtc-grps", "compressed-x-grps")replace_inp_entry(ninp, inp, "xtc-grps", "compressed-x-grps");
1788 REPL_TYPE("xtc-precision", "compressed-x-precision")replace_inp_entry(ninp, inp, "xtc-precision", "compressed-x-precision"
)
;
1789
1790 CCTYPE ("VARIOUS PREPROCESSING OPTIONS")get_estr(&ninp, &inp, "\n; " "VARIOUS PREPROCESSING OPTIONS"
, ((void*)0))
;
1791 CTYPE ("Preprocessor information: use cpp syntax.");
1792 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1793 STYPE ("include", opts->include, NULL)if ((tmp = get_estr(&ninp, &inp, "include", ((void*)0
))) != ((void*)0)) strcpy(opts->include, tmp)
;
1794 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1795 STYPE ("define", opts->define, NULL)if ((tmp = get_estr(&ninp, &inp, "define", ((void*)0)
)) != ((void*)0)) strcpy(opts->define, tmp)
;
1796
1797 CCTYPE ("RUN CONTROL PARAMETERS")get_estr(&ninp, &inp, "\n; " "RUN CONTROL PARAMETERS"
, ((void*)0))
;
1798 EETYPE("integrator", ir->eI, ei_names)ir->eI = get_eeenum(&ninp, &inp, "integrator", ei_names
, wi)
;
1799 CTYPE ("Start time and timestep in ps");
1800 RTYPE ("tinit", ir->init_t, 0.0)ir->init_t = get_ereal(&ninp, &inp, "tinit", 0.0, wi
)
;
1801 RTYPE ("dt", ir->delta_t, 0.001)ir->delta_t = get_ereal(&ninp, &inp, "dt", 0.001, wi
)
;
1802 STEPTYPE ("nsteps", ir->nsteps, 0)ir->nsteps = get_eint64(&ninp, &inp, "nsteps", 0, wi
)
;
1803 CTYPE ("For exact run continuation or redoing part of a run");
1804 STEPTYPE ("init-step", ir->init_step, 0)ir->init_step = get_eint64(&ninp, &inp, "init-step"
, 0, wi)
;
1805 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1806 ITYPE ("simulation-part", ir->simulation_part, 1)ir->simulation_part = get_eint(&ninp, &inp, "simulation-part"
, 1, wi)
;
1807 CTYPE ("mode for center of mass motion removal");
1808 EETYPE("comm-mode", ir->comm_mode, ecm_names)ir->comm_mode = get_eeenum(&ninp, &inp, "comm-mode"
, ecm_names, wi)
;
1809 CTYPE ("number of steps for center of mass motion removal");
1810 ITYPE ("nstcomm", ir->nstcomm, 100)ir->nstcomm = get_eint(&ninp, &inp, "nstcomm", 100
, wi)
;
1811 CTYPE ("group(s) for center of mass motion removal");
1812 STYPE ("comm-grps", is->vcm, NULL)if ((tmp = get_estr(&ninp, &inp, "comm-grps", ((void*
)0))) != ((void*)0)) strcpy(is->vcm, tmp)
;
1813
1814 CCTYPE ("LANGEVIN DYNAMICS OPTIONS")get_estr(&ninp, &inp, "\n; " "LANGEVIN DYNAMICS OPTIONS"
, ((void*)0))
;
1815 CTYPE ("Friction coefficient (amu/ps) and random seed");
1816 RTYPE ("bd-fric", ir->bd_fric, 0.0)ir->bd_fric = get_ereal(&ninp, &inp, "bd-fric", 0.0
, wi)
;
1817 STEPTYPE ("ld-seed", ir->ld_seed, -1)ir->ld_seed = get_eint64(&ninp, &inp, "ld-seed", -
1, wi)
;
1818
1819 /* Em stuff */
1820 CCTYPE ("ENERGY MINIMIZATION OPTIONS")get_estr(&ninp, &inp, "\n; " "ENERGY MINIMIZATION OPTIONS"
, ((void*)0))
;
1821 CTYPE ("Force tolerance and initial step-size");
1822 RTYPE ("emtol", ir->em_tol, 10.0)ir->em_tol = get_ereal(&ninp, &inp, "emtol", 10.0,
wi)
;
1823 RTYPE ("emstep", ir->em_stepsize, 0.01)ir->em_stepsize = get_ereal(&ninp, &inp, "emstep",
0.01, wi)
;
1824 CTYPE ("Max number of iterations in relax-shells");
1825 ITYPE ("niter", ir->niter, 20)ir->niter = get_eint(&ninp, &inp, "niter", 20, wi);
1826 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1827 RTYPE ("fcstep", ir->fc_stepsize, 0)ir->fc_stepsize = get_ereal(&ninp, &inp, "fcstep",
0, wi)
;
1828 CTYPE ("Frequency of steepest descents steps when doing CG");
1829 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000)ir->nstcgsteep = get_eint(&ninp, &inp, "nstcgsteep"
, 1000, wi)
;
1830 ITYPE ("nbfgscorr", ir->nbfgscorr, 10)ir->nbfgscorr = get_eint(&ninp, &inp, "nbfgscorr",
10, wi)
;
1831
1832 CCTYPE ("TEST PARTICLE INSERTION OPTIONS")get_estr(&ninp, &inp, "\n; " "TEST PARTICLE INSERTION OPTIONS"
, ((void*)0))
;
1833 RTYPE ("rtpi", ir->rtpi, 0.05)ir->rtpi = get_ereal(&ninp, &inp, "rtpi", 0.05, wi
)
;
1834
1835 /* Output options */
1836 CCTYPE ("OUTPUT CONTROL OPTIONS")get_estr(&ninp, &inp, "\n; " "OUTPUT CONTROL OPTIONS"
, ((void*)0))
;
1837 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1838 ITYPE ("nstxout", ir->nstxout, 0)ir->nstxout = get_eint(&ninp, &inp, "nstxout", 0, wi
)
;
1839 ITYPE ("nstvout", ir->nstvout, 0)ir->nstvout = get_eint(&ninp, &inp, "nstvout", 0, wi
)
;
1840 ITYPE ("nstfout", ir->nstfout, 0)ir->nstfout = get_eint(&ninp, &inp, "nstfout", 0, wi
)
;
1841 ir->nstcheckpoint = 1000;
1842 CTYPE ("Output frequency for energies to log file and energy file");
1843 ITYPE ("nstlog", ir->nstlog, 1000)ir->nstlog = get_eint(&ninp, &inp, "nstlog", 1000,
wi)
;
1844 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100)ir->nstcalcenergy = get_eint(&ninp, &inp, "nstcalcenergy"
, 100, wi)
;
1845 ITYPE ("nstenergy", ir->nstenergy, 1000)ir->nstenergy = get_eint(&ninp, &inp, "nstenergy",
1000, wi)
;
1846 CTYPE ("Output frequency and precision for .xtc file");
1847 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0)ir->nstxout_compressed = get_eint(&ninp, &inp, "nstxout-compressed"
, 0, wi)
;
1848 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0)ir->x_compression_precision = get_ereal(&ninp, &inp
, "compressed-x-precision", 1000.0, wi)
;
1849 CTYPE ("This selects the subset of atoms for the compressed");
1850 CTYPE ("trajectory file. You can select multiple groups. By");
1851 CTYPE ("default, all atoms will be written.");
1852 STYPE ("compressed-x-grps", is->x_compressed_groups, NULL)if ((tmp = get_estr(&ninp, &inp, "compressed-x-grps",
((void*)0))) != ((void*)0)) strcpy(is->x_compressed_groups
, tmp)
;
1853 CTYPE ("Selection of energy groups");
1854 STYPE ("energygrps", is->energy, NULL)if ((tmp = get_estr(&ninp, &inp, "energygrps", ((void
*)0))) != ((void*)0)) strcpy(is->energy, tmp)
;
1855
1856 /* Neighbor searching */
1857 CCTYPE ("NEIGHBORSEARCHING PARAMETERS")get_estr(&ninp, &inp, "\n; " "NEIGHBORSEARCHING PARAMETERS"
, ((void*)0))
;
1858 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1859 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names)ir->cutoff_scheme = get_eeenum(&ninp, &inp, "cutoff-scheme"
, ecutscheme_names, wi)
;
1860 CTYPE ("nblist update frequency");
1861 ITYPE ("nstlist", ir->nstlist, 10)ir->nstlist = get_eint(&ninp, &inp, "nstlist", 10,
wi)
;
1862 CTYPE ("ns algorithm (simple or grid)");
1863 EETYPE("ns-type", ir->ns_type, ens_names)ir->ns_type = get_eeenum(&ninp, &inp, "ns-type", ens_names
, wi)
;
1864 /* set ndelta to the optimal value of 2 */
1865 ir->ndelta = 2;
1866 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1867 EETYPE("pbc", ir->ePBC, epbc_names)ir->ePBC = get_eeenum(&ninp, &inp, "pbc", epbc_names
, wi)
;
1868 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names)ir->bPeriodicMols = get_eeenum(&ninp, &inp, "periodic-molecules"
, yesno_names, wi)
;
1869 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1870 CTYPE ("a value of -1 means: use rlist");
1871 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005)ir->verletbuf_tol = get_ereal(&ninp, &inp, "verlet-buffer-tolerance"
, 0.005, wi)
;
1872 CTYPE ("nblist cut-off");
1873 RTYPE ("rlist", ir->rlist, 1.0)ir->rlist = get_ereal(&ninp, &inp, "rlist", 1.0, wi
)
;
1874 CTYPE ("long-range cut-off for switched potentials");
1875 RTYPE ("rlistlong", ir->rlistlong, -1)ir->rlistlong = get_ereal(&ninp, &inp, "rlistlong"
, -1, wi)
;
1876 ITYPE ("nstcalclr", ir->nstcalclr, -1)ir->nstcalclr = get_eint(&ninp, &inp, "nstcalclr",
-1, wi)
;
1877
1878 /* Electrostatics */
1879 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW")get_estr(&ninp, &inp, "\n; " "OPTIONS FOR ELECTROSTATICS AND VDW"
, ((void*)0))
;
1880 CTYPE ("Method for doing electrostatics");
1881 EETYPE("coulombtype", ir->coulombtype, eel_names)ir->coulombtype = get_eeenum(&ninp, &inp, "coulombtype"
, eel_names, wi)
;
1882 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names)ir->coulomb_modifier = get_eeenum(&ninp, &inp, "coulomb-modifier"
, eintmod_names, wi)
;
1883 CTYPE ("cut-off lengths");
1884 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0)ir->rcoulomb_switch = get_ereal(&ninp, &inp, "rcoulomb-switch"
, 0.0, wi)
;
1885 RTYPE ("rcoulomb", ir->rcoulomb, 1.0)ir->rcoulomb = get_ereal(&ninp, &inp, "rcoulomb", 1.0
, wi)
;
1886 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1887 RTYPE ("epsilon-r", ir->epsilon_r, 1.0)ir->epsilon_r = get_ereal(&ninp, &inp, "epsilon-r"
, 1.0, wi)
;
1888 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0)ir->epsilon_rf = get_ereal(&ninp, &inp, "epsilon-rf"
, 0.0, wi)
;
1889 CTYPE ("Method for doing Van der Waals");
1890 EETYPE("vdw-type", ir->vdwtype, evdw_names)ir->vdwtype = get_eeenum(&ninp, &inp, "vdw-type", evdw_names
, wi)
;
1891 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names)ir->vdw_modifier = get_eeenum(&ninp, &inp, "vdw-modifier"
, eintmod_names, wi)
;
1892 CTYPE ("cut-off lengths");
1893 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0)ir->rvdw_switch = get_ereal(&ninp, &inp, "rvdw-switch"
, 0.0, wi)
;
1894 RTYPE ("rvdw", ir->rvdw, 1.0)ir->rvdw = get_ereal(&ninp, &inp, "rvdw", 1.0, wi);
1895 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1896 EETYPE("DispCorr", ir->eDispCorr, edispc_names)ir->eDispCorr = get_eeenum(&ninp, &inp, "DispCorr"
, edispc_names, wi)
;
1897 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1898 RTYPE ("table-extension", ir->tabext, 1.0)ir->tabext = get_ereal(&ninp, &inp, "table-extension"
, 1.0, wi)
;
1899 CTYPE ("Separate tables between energy group pairs");
1900 STYPE ("energygrp-table", is->egptable, NULL)if ((tmp = get_estr(&ninp, &inp, "energygrp-table", (
(void*)0))) != ((void*)0)) strcpy(is->egptable, tmp)
;
1901 CTYPE ("Spacing for the PME/PPPM FFT grid");
1902 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12)ir->fourier_spacing = get_ereal(&ninp, &inp, "fourierspacing"
, 0.12, wi)
;
1903 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1904 ITYPE ("fourier-nx", ir->nkx, 0)ir->nkx = get_eint(&ninp, &inp, "fourier-nx", 0, wi
)
;
1905 ITYPE ("fourier-ny", ir->nky, 0)ir->nky = get_eint(&ninp, &inp, "fourier-ny", 0, wi
)
;
1906 ITYPE ("fourier-nz", ir->nkz, 0)ir->nkz = get_eint(&ninp, &inp, "fourier-nz", 0, wi
)
;
1907 CTYPE ("EWALD/PME/PPPM parameters");
1908 ITYPE ("pme-order", ir->pme_order, 4)ir->pme_order = get_eint(&ninp, &inp, "pme-order",
4, wi)
;
1909 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001)ir->ewald_rtol = get_ereal(&ninp, &inp, "ewald-rtol"
, 0.00001, wi)
;
1910 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001)ir->ewald_rtol_lj = get_ereal(&ninp, &inp, "ewald-rtol-lj"
, 0.001, wi)
;
1911 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names)ir->ljpme_combination_rule = get_eeenum(&ninp, &inp
, "lj-pme-comb-rule", eljpme_names, wi)
;
1912 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names)ir->ewald_geometry = get_eeenum(&ninp, &inp, "ewald-geometry"
, eewg_names, wi)
;
1913 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0)ir->epsilon_surface = get_ereal(&ninp, &inp, "epsilon-surface"
, 0.0, wi)
;
1914 EETYPE("optimize-fft", ir->bOptFFT, yesno_names)ir->bOptFFT = get_eeenum(&ninp, &inp, "optimize-fft"
, yesno_names, wi)
;
1915
1916 CCTYPE("IMPLICIT SOLVENT ALGORITHM")get_estr(&ninp, &inp, "\n; " "IMPLICIT SOLVENT ALGORITHM"
, ((void*)0))
;
1917 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names)ir->implicit_solvent = get_eeenum(&ninp, &inp, "implicit-solvent"
, eis_names, wi)
;
1918
1919 CCTYPE ("GENERALIZED BORN ELECTROSTATICS")get_estr(&ninp, &inp, "\n; " "GENERALIZED BORN ELECTROSTATICS"
, ((void*)0))
;
1920 CTYPE ("Algorithm for calculating Born radii");
1921 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names)ir->gb_algorithm = get_eeenum(&ninp, &inp, "gb-algorithm"
, egb_names, wi)
;
1922 CTYPE ("Frequency of calculating the Born radii inside rlist");
1923 ITYPE ("nstgbradii", ir->nstgbradii, 1)ir->nstgbradii = get_eint(&ninp, &inp, "nstgbradii"
, 1, wi)
;
1924 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1925 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1926 RTYPE ("rgbradii", ir->rgbradii, 1.0)ir->rgbradii = get_ereal(&ninp, &inp, "rgbradii", 1.0
, wi)
;
1927 CTYPE ("Dielectric coefficient of the implicit solvent");
1928 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0)ir->gb_epsilon_solvent = get_ereal(&ninp, &inp, "gb-epsilon-solvent"
, 80.0, wi)
;
1929 CTYPE ("Salt concentration in M for Generalized Born models");
1930 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0)ir->gb_saltconc = get_ereal(&ninp, &inp, "gb-saltconc"
, 0.0, wi)
;
1931 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1932 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0)ir->gb_obc_alpha = get_ereal(&ninp, &inp, "gb-obc-alpha"
, 1.0, wi)
;
1933 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8)ir->gb_obc_beta = get_ereal(&ninp, &inp, "gb-obc-beta"
, 0.8, wi)
;
1934 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85)ir->gb_obc_gamma = get_ereal(&ninp, &inp, "gb-obc-gamma"
, 4.85, wi)
;
1935 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009)ir->gb_dielectric_offset = get_ereal(&ninp, &inp, "gb-dielectric-offset"
, 0.009, wi)
;
1936 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names)ir->sa_algorithm = get_eeenum(&ninp, &inp, "sa-algorithm"
, esa_names, wi)
;
1937 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1938 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1939 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1)ir->sa_surface_tension = get_ereal(&ninp, &inp, "sa-surface-tension"
, -1, wi)
;
1940
1941 /* Coupling stuff */
1942 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS")get_estr(&ninp, &inp, "\n; " "OPTIONS FOR WEAK COUPLING ALGORITHMS"
, ((void*)0))
;
1943 CTYPE ("Temperature coupling");
1944 EETYPE("tcoupl", ir->etc, etcoupl_names)ir->etc = get_eeenum(&ninp, &inp, "tcoupl", etcoupl_names
, wi)
;
1945 ITYPE ("nsttcouple", ir->nsttcouple, -1)ir->nsttcouple = get_eint(&ninp, &inp, "nsttcouple"
, -1, wi)
;
1946 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10)ir->opts.nhchainlength = get_eint(&ninp, &inp, "nh-chain-length"
, 10, wi)
;
1947 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names)ir->bPrintNHChains = get_eeenum(&ninp, &inp, "print-nose-hoover-chain-variables"
, yesno_names, wi)
;
1948 CTYPE ("Groups to couple separately");
1949 STYPE ("tc-grps", is->tcgrps, NULL)if ((tmp = get_estr(&ninp, &inp, "tc-grps", ((void*)0
))) != ((void*)0)) strcpy(is->tcgrps, tmp)
;
1950 CTYPE ("Time constant (ps) and reference temperature (K)");
1951 STYPE ("tau-t", is->tau_t, NULL)if ((tmp = get_estr(&ninp, &inp, "tau-t", ((void*)0))
) != ((void*)0)) strcpy(is->tau_t, tmp)
;
1952 STYPE ("ref-t", is->ref_t, NULL)if ((tmp = get_estr(&ninp, &inp, "ref-t", ((void*)0))
) != ((void*)0)) strcpy(is->ref_t, tmp)
;
1953 CTYPE ("pressure coupling");
1954 EETYPE("pcoupl", ir->epc, epcoupl_names)ir->epc = get_eeenum(&ninp, &inp, "pcoupl", epcoupl_names
, wi)
;
1955 EETYPE("pcoupltype", ir->epct, epcoupltype_names)ir->epct = get_eeenum(&ninp, &inp, "pcoupltype", epcoupltype_names
, wi)
;
1956 ITYPE ("nstpcouple", ir->nstpcouple, -1)ir->nstpcouple = get_eint(&ninp, &inp, "nstpcouple"
, -1, wi)
;
1957 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1958 RTYPE ("tau-p", ir->tau_p, 1.0)ir->tau_p = get_ereal(&ninp, &inp, "tau-p", 1.0, wi
)
;
1959 STYPE ("compressibility", dumstr[0], NULL)if ((tmp = get_estr(&ninp, &inp, "compressibility", (
(void*)0))) != ((void*)0)) strcpy(dumstr[0], tmp)
;
1960 STYPE ("ref-p", dumstr[1], NULL)if ((tmp = get_estr(&ninp, &inp, "ref-p", ((void*)0))
) != ((void*)0)) strcpy(dumstr[1], tmp)
;
1961 CTYPE ("Scaling of reference coordinates, No, All or COM");
1962 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names)ir->refcoord_scaling = get_eeenum(&ninp, &inp, "refcoord-scaling"
, erefscaling_names, wi)
;
1963
1964 /* QMMM */
1965 CCTYPE ("OPTIONS FOR QMMM calculations")get_estr(&ninp, &inp, "\n; " "OPTIONS FOR QMMM calculations"
, ((void*)0))
;
1966 EETYPE("QMMM", ir->bQMMM, yesno_names)ir->bQMMM = get_eeenum(&ninp, &inp, "QMMM", yesno_names
, wi)
;
1967 CTYPE ("Groups treated Quantum Mechanically");
1968 STYPE ("QMMM-grps", is->QMMM, NULL)if ((tmp = get_estr(&ninp, &inp, "QMMM-grps", ((void*
)0))) != ((void*)0)) strcpy(is->QMMM, tmp)
;
1969 CTYPE ("QM method");
1970 STYPE("QMmethod", is->QMmethod, NULL)if ((tmp = get_estr(&ninp, &inp, "QMmethod", ((void*)
0))) != ((void*)0)) strcpy(is->QMmethod, tmp)
;
1971 CTYPE ("QMMM scheme");
1972 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names)ir->QMMMscheme = get_eeenum(&ninp, &inp, "QMMMscheme"
, eQMMMscheme_names, wi)
;
1973 CTYPE ("QM basisset");
1974 STYPE("QMbasis", is->QMbasis, NULL)if ((tmp = get_estr(&ninp, &inp, "QMbasis", ((void*)0
))) != ((void*)0)) strcpy(is->QMbasis, tmp)
;
1975 CTYPE ("QM charge");
1976 STYPE ("QMcharge", is->QMcharge, NULL)if ((tmp = get_estr(&ninp, &inp, "QMcharge", ((void*)
0))) != ((void*)0)) strcpy(is->QMcharge, tmp)
;
1977 CTYPE ("QM multiplicity");
1978 STYPE ("QMmult", is->QMmult, NULL)if ((tmp = get_estr(&ninp, &inp, "QMmult", ((void*)0)
)) != ((void*)0)) strcpy(is->QMmult, tmp)
;
1979 CTYPE ("Surface Hopping");
1980 STYPE ("SH", is->bSH, NULL)if ((tmp = get_estr(&ninp, &inp, "SH", ((void*)0))) !=
((void*)0)) strcpy(is->bSH, tmp)
;
1981 CTYPE ("CAS space options");
1982 STYPE ("CASorbitals", is->CASorbitals, NULL)if ((tmp = get_estr(&ninp, &inp, "CASorbitals", ((void
*)0))) != ((void*)0)) strcpy(is->CASorbitals, tmp)
;
1983 STYPE ("CASelectrons", is->CASelectrons, NULL)if ((tmp = get_estr(&ninp, &inp, "CASelectrons", ((void
*)0))) != ((void*)0)) strcpy(is->CASelectrons, tmp)
;
1984 STYPE ("SAon", is->SAon, NULL)if ((tmp = get_estr(&ninp, &inp, "SAon", ((void*)0)))
!= ((void*)0)) strcpy(is->SAon, tmp)
;
1985 STYPE ("SAoff", is->SAoff, NULL)if ((tmp = get_estr(&ninp, &inp, "SAoff", ((void*)0))
) != ((void*)0)) strcpy(is->SAoff, tmp)
;
1986 STYPE ("SAsteps", is->SAsteps, NULL)if ((tmp = get_estr(&ninp, &inp, "SAsteps", ((void*)0
))) != ((void*)0)) strcpy(is->SAsteps, tmp)
;
1987 CTYPE ("Scale factor for MM charges");
1988 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0)ir->scalefactor = get_ereal(&ninp, &inp, "MMChargeScaleFactor"
, 1.0, wi)
;
1989 CTYPE ("Optimization of QM subsystem");
1990 STYPE ("bOPT", is->bOPT, NULL)if ((tmp = get_estr(&ninp, &inp, "bOPT", ((void*)0)))
!= ((void*)0)) strcpy(is->bOPT, tmp)
;
1991 STYPE ("bTS", is->bTS, NULL)if ((tmp = get_estr(&ninp, &inp, "bTS", ((void*)0))) !=
((void*)0)) strcpy(is->bTS, tmp)
;
1992
1993 /* Simulated annealing */
1994 CCTYPE("SIMULATED ANNEALING")get_estr(&ninp, &inp, "\n; " "SIMULATED ANNEALING", (
(void*)0))
;
1995 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1996 STYPE ("annealing", is->anneal, NULL)if ((tmp = get_estr(&ninp, &inp, "annealing", ((void*
)0))) != ((void*)0)) strcpy(is->anneal, tmp)
;
1997 CTYPE ("Number of time points to use for specifying annealing in each group");
1998 STYPE ("annealing-npoints", is->anneal_npoints, NULL)if ((tmp = get_estr(&ninp, &inp, "annealing-npoints",
((void*)0))) != ((void*)0)) strcpy(is->anneal_npoints, tmp
)
;
1999 CTYPE ("List of times at the annealing points for each group");
2000 STYPE ("annealing-time", is->anneal_time, NULL)if ((tmp = get_estr(&ninp, &inp, "annealing-time", ((
void*)0))) != ((void*)0)) strcpy(is->anneal_time, tmp)
;
2001 CTYPE ("Temp. at each annealing point, for each group.");
2002 STYPE ("annealing-temp", is->anneal_temp, NULL)if ((tmp = get_estr(&ninp, &inp, "annealing-temp", ((
void*)0))) != ((void*)0)) strcpy(is->anneal_temp, tmp)
;
2003
2004 /* Startup run */
2005 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN")get_estr(&ninp, &inp, "\n; " "GENERATE VELOCITIES FOR STARTUP RUN"
, ((void*)0))
;
2006 EETYPE("gen-vel", opts->bGenVel, yesno_names)opts->bGenVel = get_eeenum(&ninp, &inp, "gen-vel",
yesno_names, wi)
;
2007 RTYPE ("gen-temp", opts->tempi, 300.0)opts->tempi = get_ereal(&ninp, &inp, "gen-temp", 300.0
, wi)
;
2008 ITYPE ("gen-seed", opts->seed, -1)opts->seed = get_eint(&ninp, &inp, "gen-seed", -1,
wi)
;
2009
2010 /* Shake stuff */
2011 CCTYPE ("OPTIONS FOR BONDS")get_estr(&ninp, &inp, "\n; " "OPTIONS FOR BONDS", ((void
*)0))
;
2012 EETYPE("constraints", opts->nshake, constraints)opts->nshake = get_eeenum(&ninp, &inp, "constraints"
, constraints, wi)
;
2013 CTYPE ("Type of constraint algorithm");
2014 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names)ir->eConstrAlg = get_eeenum(&ninp, &inp, "constraint-algorithm"
, econstr_names, wi)
;
2015 CTYPE ("Do not constrain the start configuration");
2016 EETYPE("continuation", ir->bContinuation, yesno_names)ir->bContinuation = get_eeenum(&ninp, &inp, "continuation"
, yesno_names, wi)
;
2017 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2018 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names)ir->bShakeSOR = get_eeenum(&ninp, &inp, "Shake-SOR"
, yesno_names, wi)
;
2019 CTYPE ("Relative tolerance of shake");
2020 RTYPE ("shake-tol", ir->shake_tol, 0.0001)ir->shake_tol = get_ereal(&ninp, &inp, "shake-tol"
, 0.0001, wi)
;
2021 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2022 ITYPE ("lincs-order", ir->nProjOrder, 4)ir->nProjOrder = get_eint(&ninp, &inp, "lincs-order"
, 4, wi)
;
2023 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2024 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2025 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2026 ITYPE ("lincs-iter", ir->nLincsIter, 1)ir->nLincsIter = get_eint(&ninp, &inp, "lincs-iter"
, 1, wi)
;
2027 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2028 CTYPE ("rotates over more degrees than");
2029 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0)ir->LincsWarnAngle = get_ereal(&ninp, &inp, "lincs-warnangle"
, 30.0, wi)
;
2030 CTYPE ("Convert harmonic bonds to morse potentials");
2031 EETYPE("morse", opts->bMorse, yesno_names)opts->bMorse = get_eeenum(&ninp, &inp, "morse", yesno_names
, wi)
;
2032
2033 /* Energy group exclusions */
2034 CCTYPE ("ENERGY GROUP EXCLUSIONS")get_estr(&ninp, &inp, "\n; " "ENERGY GROUP EXCLUSIONS"
, ((void*)0))
;
2035 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2036 STYPE ("energygrp-excl", is->egpexcl, NULL)if ((tmp = get_estr(&ninp, &inp, "energygrp-excl", ((
void*)0))) != ((void*)0)) strcpy(is->egpexcl, tmp)
;
2037
2038 /* Walls */
2039 CCTYPE ("WALLS")get_estr(&ninp, &inp, "\n; " "WALLS", ((void*)0));
2040 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2041 ITYPE ("nwall", ir->nwall, 0)ir->nwall = get_eint(&ninp, &inp, "nwall", 0, wi);
2042 EETYPE("wall-type", ir->wall_type, ewt_names)ir->wall_type = get_eeenum(&ninp, &inp, "wall-type"
, ewt_names, wi)
;
2043 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1)ir->wall_r_linpot = get_ereal(&ninp, &inp, "wall-r-linpot"
, -1, wi)
;
2044 STYPE ("wall-atomtype", is->wall_atomtype, NULL)if ((tmp = get_estr(&ninp, &inp, "wall-atomtype", ((void
*)0))) != ((void*)0)) strcpy(is->wall_atomtype, tmp)
;
2045 STYPE ("wall-density", is->wall_density, NULL)if ((tmp = get_estr(&ninp, &inp, "wall-density", ((void
*)0))) != ((void*)0)) strcpy(is->wall_density, tmp)
;
2046 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3)ir->wall_ewald_zfac = get_ereal(&ninp, &inp, "wall-ewald-zfac"
, 3, wi)
;
2047
2048 /* COM pulling */
2049 CCTYPE("COM PULLING")get_estr(&ninp, &inp, "\n; " "COM PULLING", ((void*)0
))
;
2050 CTYPE("Pull type: no, umbrella, constraint or constant-force");
2051 EETYPE("pull", ir->ePull, epull_names)ir->ePull = get_eeenum(&ninp, &inp, "pull", epull_names
, wi)
;
2052 if (ir->ePull != epullNO)
2053 {
2054 snew(ir->pull, 1)(ir->pull) = save_calloc("ir->pull", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2054, (1), sizeof(*(ir->pull)))
;
2055 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2056 }
2057
2058 /* Enforced rotation */
2059 CCTYPE("ENFORCED ROTATION")get_estr(&ninp, &inp, "\n; " "ENFORCED ROTATION", ((void
*)0))
;
2060 CTYPE("Enforced rotation: No or Yes");
2061 EETYPE("rotation", ir->bRot, yesno_names)ir->bRot = get_eeenum(&ninp, &inp, "rotation", yesno_names
, wi)
;
2062 if (ir->bRot)
2063 {
2064 snew(ir->rot, 1)(ir->rot) = save_calloc("ir->rot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2064, (1), sizeof(*(ir->rot)))
;
2065 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2066 }
2067
2068 /* Interactive MD */
2069 ir->bIMD = FALSE0;
2070 CCTYPE("Group to display and/or manipulate in interactive MD session")get_estr(&ninp, &inp, "\n; " "Group to display and/or manipulate in interactive MD session"
, ((void*)0))
;
2071 STYPE ("IMD-group", is->imd_grp, NULL)if ((tmp = get_estr(&ninp, &inp, "IMD-group", ((void*
)0))) != ((void*)0)) strcpy(is->imd_grp, tmp)
;
2072 if (is->imd_grp[0] != '\0')
2073 {
2074 snew(ir->imd, 1)(ir->imd) = save_calloc("ir->imd", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2074, (1), sizeof(*(ir->imd)))
;
2075 ir->bIMD = TRUE1;
2076 }
2077
2078 /* Refinement */
2079 CCTYPE("NMR refinement stuff")get_estr(&ninp, &inp, "\n; " "NMR refinement stuff", (
(void*)0))
;
2080 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2081 EETYPE("disre", ir->eDisre, edisre_names)ir->eDisre = get_eeenum(&ninp, &inp, "disre", edisre_names
, wi)
;
2082 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2083 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names)ir->eDisreWeighting = get_eeenum(&ninp, &inp, "disre-weighting"
, edisreweighting_names, wi)
;
2084 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2085 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names)ir->bDisreMixed = get_eeenum(&ninp, &inp, "disre-mixed"
, yesno_names, wi)
;
2086 RTYPE ("disre-fc", ir->dr_fc, 1000.0)ir->dr_fc = get_ereal(&ninp, &inp, "disre-fc", 1000.0
, wi)
;
2087 RTYPE ("disre-tau", ir->dr_tau, 0.0)ir->dr_tau = get_ereal(&ninp, &inp, "disre-tau", 0.0
, wi)
;
2088 CTYPE ("Output frequency for pair distances to energy file");
2089 ITYPE ("nstdisreout", ir->nstdisreout, 100)ir->nstdisreout = get_eint(&ninp, &inp, "nstdisreout"
, 100, wi)
;
2090 CTYPE ("Orientation restraints: No or Yes");
2091 EETYPE("orire", opts->bOrire, yesno_names)opts->bOrire = get_eeenum(&ninp, &inp, "orire", yesno_names
, wi)
;
2092 CTYPE ("Orientation restraints force constant and tau for time averaging");
2093 RTYPE ("orire-fc", ir->orires_fc, 0.0)ir->orires_fc = get_ereal(&ninp, &inp, "orire-fc",
0.0, wi)
;
2094 RTYPE ("orire-tau", ir->orires_tau, 0.0)ir->orires_tau = get_ereal(&ninp, &inp, "orire-tau"
, 0.0, wi)
;
2095 STYPE ("orire-fitgrp", is->orirefitgrp, NULL)if ((tmp = get_estr(&ninp, &inp, "orire-fitgrp", ((void
*)0))) != ((void*)0)) strcpy(is->orirefitgrp, tmp)
;
2096 CTYPE ("Output frequency for trace(SD) and S to energy file");
2097 ITYPE ("nstorireout", ir->nstorireout, 100)ir->nstorireout = get_eint(&ninp, &inp, "nstorireout"
, 100, wi)
;
2098
2099 /* free energy variables */
2100 CCTYPE ("Free energy variables")get_estr(&ninp, &inp, "\n; " "Free energy variables",
((void*)0))
;
2101 EETYPE("free-energy", ir->efep, efep_names)ir->efep = get_eeenum(&ninp, &inp, "free-energy", efep_names
, wi)
;
2102 STYPE ("couple-moltype", is->couple_moltype, NULL)if ((tmp = get_estr(&ninp, &inp, "couple-moltype", ((
void*)0))) != ((void*)0)) strcpy(is->couple_moltype, tmp)
;
2103 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam)opts->couple_lam0 = get_eeenum(&ninp, &inp, "couple-lambda0"
, couple_lam, wi)
;
2104 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam)opts->couple_lam1 = get_eeenum(&ninp, &inp, "couple-lambda1"
, couple_lam, wi)
;
2105 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names)opts->bCoupleIntra = get_eeenum(&ninp, &inp, "couple-intramol"
, yesno_names, wi)
;
2106
2107 RTYPE ("init-lambda", fep->init_lambda, -1)fep->init_lambda = get_ereal(&ninp, &inp, "init-lambda"
, -1, wi)
; /* start with -1 so
2108 we can recognize if
2109 it was not entered */
2110 ITYPE ("init-lambda-state", fep->init_fep_state, -1)fep->init_fep_state = get_eint(&ninp, &inp, "init-lambda-state"
, -1, wi)
;
2111 RTYPE ("delta-lambda", fep->delta_lambda, 0.0)fep->delta_lambda = get_ereal(&ninp, &inp, "delta-lambda"
, 0.0, wi)
;
2112 ITYPE ("nstdhdl", fep->nstdhdl, 50)fep->nstdhdl = get_eint(&ninp, &inp, "nstdhdl", 50
, wi)
;
2113 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL)if ((tmp = get_estr(&ninp, &inp, "fep-lambdas", ((void
*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptFEP], tmp)
;
2114 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL)if ((tmp = get_estr(&ninp, &inp, "mass-lambdas", ((void
*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptMASS], tmp
)
;
2115 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL)if ((tmp = get_estr(&ninp, &inp, "coul-lambdas", ((void
*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptCOUL], tmp
)
;
2116 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL)if ((tmp = get_estr(&ninp, &inp, "vdw-lambdas", ((void
*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptVDW], tmp)
;
2117 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL)if ((tmp = get_estr(&ninp, &inp, "bonded-lambdas", ((
void*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptBONDED
], tmp)
;
2118 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL)if ((tmp = get_estr(&ninp, &inp, "restraint-lambdas",
((void*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptRESTRAINT
], tmp)
;
2119 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL)if ((tmp = get_estr(&ninp, &inp, "temperature-lambdas"
, ((void*)0))) != ((void*)0)) strcpy(is->fep_lambda[efptTEMPERATURE
], tmp)
;
2120 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1)fep->lambda_neighbors = get_eint(&ninp, &inp, "calc-lambda-neighbors"
, 1, wi)
;
2121 STYPE ("init-lambda-weights", is->lambda_weights, NULL)if ((tmp = get_estr(&ninp, &inp, "init-lambda-weights"
, ((void*)0))) != ((void*)0)) strcpy(is->lambda_weights, tmp
)
;
2122 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names)fep->bPrintEnergy = get_eeenum(&ninp, &inp, "dhdl-print-energy"
, yesno_names, wi)
;
2123 RTYPE ("sc-alpha", fep->sc_alpha, 0.0)fep->sc_alpha = get_ereal(&ninp, &inp, "sc-alpha",
0.0, wi)
;
2124 ITYPE ("sc-power", fep->sc_power, 1)fep->sc_power = get_eint(&ninp, &inp, "sc-power", 1
, wi)
;
2125 RTYPE ("sc-r-power", fep->sc_r_power, 6.0)fep->sc_r_power = get_ereal(&ninp, &inp, "sc-r-power"
, 6.0, wi)
;
2126 RTYPE ("sc-sigma", fep->sc_sigma, 0.3)fep->sc_sigma = get_ereal(&ninp, &inp, "sc-sigma",
0.3, wi)
;
2127 EETYPE("sc-coul", fep->bScCoul, yesno_names)fep->bScCoul = get_eeenum(&ninp, &inp, "sc-coul", yesno_names
, wi)
;
2128 ITYPE ("dh_hist_size", fep->dh_hist_size, 0)fep->dh_hist_size = get_eint(&ninp, &inp, "dh_hist_size"
, 0, wi)
;
2129 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1)fep->dh_hist_spacing = get_ereal(&ninp, &inp, "dh_hist_spacing"
, 0.1, wi)
;
2130 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,fep->separate_dhdl_file = get_eeenum(&ninp, &inp, "separate-dhdl-file"
, separate_dhdl_file_names, wi)
2131 separate_dhdl_file_names)fep->separate_dhdl_file = get_eeenum(&ninp, &inp, "separate-dhdl-file"
, separate_dhdl_file_names, wi)
;
2132 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names)fep->dhdl_derivatives = get_eeenum(&ninp, &inp, "dhdl-derivatives"
, dhdl_derivatives_names, wi)
;
2133 ITYPE ("dh_hist_size", fep->dh_hist_size, 0)fep->dh_hist_size = get_eint(&ninp, &inp, "dh_hist_size"
, 0, wi)
;
2134 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1)fep->dh_hist_spacing = get_ereal(&ninp, &inp, "dh_hist_spacing"
, 0.1, wi)
;
2135
2136 /* Non-equilibrium MD stuff */
2137 CCTYPE("Non-equilibrium MD stuff")get_estr(&ninp, &inp, "\n; " "Non-equilibrium MD stuff"
, ((void*)0))
;
2138 STYPE ("acc-grps", is->accgrps, NULL)if ((tmp = get_estr(&ninp, &inp, "acc-grps", ((void*)
0))) != ((void*)0)) strcpy(is->accgrps, tmp)
;
2139 STYPE ("accelerate", is->acc, NULL)if ((tmp = get_estr(&ninp, &inp, "accelerate", ((void
*)0))) != ((void*)0)) strcpy(is->acc, tmp)
;
2140 STYPE ("freezegrps", is->freeze, NULL)if ((tmp = get_estr(&ninp, &inp, "freezegrps", ((void
*)0))) != ((void*)0)) strcpy(is->freeze, tmp)
;
2141 STYPE ("freezedim", is->frdim, NULL)if ((tmp = get_estr(&ninp, &inp, "freezedim", ((void*
)0))) != ((void*)0)) strcpy(is->frdim, tmp)
;
2142 RTYPE ("cos-acceleration", ir->cos_accel, 0)ir->cos_accel = get_ereal(&ninp, &inp, "cos-acceleration"
, 0, wi)
;
2143 STYPE ("deform", is->deform, NULL)if ((tmp = get_estr(&ninp, &inp, "deform", ((void*)0)
)) != ((void*)0)) strcpy(is->deform, tmp)
;
2144
2145 /* simulated tempering variables */
2146 CCTYPE("simulated tempering variables")get_estr(&ninp, &inp, "\n; " "simulated tempering variables"
, ((void*)0))
;
2147 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names)ir->bSimTemp = get_eeenum(&ninp, &inp, "simulated-tempering"
, yesno_names, wi)
;
2148 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names)ir->simtempvals->eSimTempScale = get_eeenum(&ninp, &
inp, "simulated-tempering-scaling", esimtemp_names, wi)
;
2149 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0)ir->simtempvals->simtemp_low = get_ereal(&ninp, &
inp, "sim-temp-low", 300.0, wi)
;
2150 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0)ir->simtempvals->simtemp_high = get_ereal(&ninp, &
inp, "sim-temp-high", 300.0, wi)
;
2151
2152 /* expanded ensemble variables */
2153 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2154 {
2155 read_expandedparams(&ninp, &inp, expand, wi);
2156 }
2157
2158 /* Electric fields */
2159 CCTYPE("Electric fields")get_estr(&ninp, &inp, "\n; " "Electric fields", ((void
*)0))
;
2160 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2161 CTYPE ("and a phase angle (real)");
2162 STYPE ("E-x", is->efield_x, NULL)if ((tmp = get_estr(&ninp, &inp, "E-x", ((void*)0))) !=
((void*)0)) strcpy(is->efield_x, tmp)
;
2163 STYPE ("E-xt", is->efield_xt, NULL)if ((tmp = get_estr(&ninp, &inp, "E-xt", ((void*)0)))
!= ((void*)0)) strcpy(is->efield_xt, tmp)
;
2164 STYPE ("E-y", is->efield_y, NULL)if ((tmp = get_estr(&ninp, &inp, "E-y", ((void*)0))) !=
((void*)0)) strcpy(is->efield_y, tmp)
;
2165 STYPE ("E-yt", is->efield_yt, NULL)if ((tmp = get_estr(&ninp, &inp, "E-yt", ((void*)0)))
!= ((void*)0)) strcpy(is->efield_yt, tmp)
;
2166 STYPE ("E-z", is->efield_z, NULL)if ((tmp = get_estr(&ninp, &inp, "E-z", ((void*)0))) !=
((void*)0)) strcpy(is->efield_z, tmp)
;
2167 STYPE ("E-zt", is->efield_zt, NULL)if ((tmp = get_estr(&ninp, &inp, "E-zt", ((void*)0)))
!= ((void*)0)) strcpy(is->efield_zt, tmp)
;
2168
2169 CCTYPE("Ion/water position swapping for computational electrophysiology setups")get_estr(&ninp, &inp, "\n; " "Ion/water position swapping for computational electrophysiology setups"
, ((void*)0))
;
2170 CTYPE("Swap positions along direction: no, X, Y, Z");
2171 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names)ir->eSwapCoords = get_eeenum(&ninp, &inp, "swapcoords"
, eSwapTypes_names, wi)
;
2172 if (ir->eSwapCoords != eswapNO)
2173 {
2174 snew(ir->swap, 1)(ir->swap) = save_calloc("ir->swap", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2174, (1), sizeof(*(ir->swap)))
;
2175 CTYPE("Swap attempt frequency");
2176 ITYPE("swap-frequency", ir->swap->nstswap, 1)ir->swap->nstswap = get_eint(&ninp, &inp, "swap-frequency"
, 1, wi)
;
2177 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2178 STYPE("split-group0", splitgrp0, NULL)if ((tmp = get_estr(&ninp, &inp, "split-group0", ((void
*)0))) != ((void*)0)) strcpy(splitgrp0, tmp)
;
2179 STYPE("split-group1", splitgrp1, NULL)if ((tmp = get_estr(&ninp, &inp, "split-group1", ((void
*)0))) != ((void*)0)) strcpy(splitgrp1, tmp)
;
2180 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2181 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names)ir->swap->massw_split[0] = get_eeenum(&ninp, &inp
, "massw-split0", yesno_names, wi)
;
2182 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names)ir->swap->massw_split[1] = get_eeenum(&ninp, &inp
, "massw-split1", yesno_names, wi)
;
2183
2184 CTYPE("Group name of ions that can be exchanged with solvent molecules");
2185 STYPE("swap-group", swapgrp, NULL)if ((tmp = get_estr(&ninp, &inp, "swap-group", ((void
*)0))) != ((void*)0)) strcpy(swapgrp, tmp)
;
2186 CTYPE("Group name of solvent molecules");
2187 STYPE("solvent-group", solgrp, NULL)if ((tmp = get_estr(&ninp, &inp, "solvent-group", ((void
*)0))) != ((void*)0)) strcpy(solgrp, tmp)
;
2188
2189 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2190 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2191 CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2192 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0)ir->swap->cyl0r = get_ereal(&ninp, &inp, "cyl0-r"
, 2.0, wi)
;
2193 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0)ir->swap->cyl0u = get_ereal(&ninp, &inp, "cyl0-up"
, 1.0, wi)
;
2194 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0)ir->swap->cyl0l = get_ereal(&ninp, &inp, "cyl0-down"
, 1.0, wi)
;
2195 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0)ir->swap->cyl1r = get_ereal(&ninp, &inp, "cyl1-r"
, 2.0, wi)
;
2196 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0)ir->swap->cyl1u = get_ereal(&ninp, &inp, "cyl1-up"
, 1.0, wi)
;
2197 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0)ir->swap->cyl1l = get_ereal(&ninp, &inp, "cyl1-down"
, 1.0, wi)
;
2198
2199 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2200 ITYPE("coupl-steps", ir->swap->nAverage, 10)ir->swap->nAverage = get_eint(&ninp, &inp, "coupl-steps"
, 10, wi)
;
2201 CTYPE("Requested number of anions and cations for each of the two compartments");
2202 CTYPE("-1 means fix the numbers as found in time step 0");
2203 ITYPE("anionsA", ir->swap->nanions[0], -1)ir->swap->nanions[0] = get_eint(&ninp, &inp, "anionsA"
, -1, wi)
;
2204 ITYPE("cationsA", ir->swap->ncations[0], -1)ir->swap->ncations[0] = get_eint(&ninp, &inp, "cationsA"
, -1, wi)
;
2205 ITYPE("anionsB", ir->swap->nanions[1], -1)ir->swap->nanions[1] = get_eint(&ninp, &inp, "anionsB"
, -1, wi)
;
2206 ITYPE("cationsB", ir->swap->ncations[1], -1)ir->swap->ncations[1] = get_eint(&ninp, &inp, "cationsB"
, -1, wi)
;
2207 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2208 RTYPE("threshold", ir->swap->threshold, 1.0)ir->swap->threshold = get_ereal(&ninp, &inp, "threshold"
, 1.0, wi)
;
2209 }
2210
2211 /* AdResS defined thingies */
2212 CCTYPE ("AdResS parameters")get_estr(&ninp, &inp, "\n; " "AdResS parameters", ((void
*)0))
;
2213 EETYPE("adress", ir->bAdress, yesno_names)ir->bAdress = get_eeenum(&ninp, &inp, "adress", yesno_names
, wi)
;
2214 if (ir->bAdress)
2215 {
2216 snew(ir->adress, 1)(ir->adress) = save_calloc("ir->adress", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2216, (1), sizeof(*(ir->adress)))
;
2217 read_adressparams(&ninp, &inp, ir->adress, wi);
2218 }
2219
2220 /* User defined thingies */
2221 CCTYPE ("User defined thingies")get_estr(&ninp, &inp, "\n; " "User defined thingies",
((void*)0))
;
2222 STYPE ("user1-grps", is->user1, NULL)if ((tmp = get_estr(&ninp, &inp, "user1-grps", ((void
*)0))) != ((void*)0)) strcpy(is->user1, tmp)
;
2223 STYPE ("user2-grps", is->user2, NULL)if ((tmp = get_estr(&ninp, &inp, "user2-grps", ((void
*)0))) != ((void*)0)) strcpy(is->user2, tmp)
;
2224 ITYPE ("userint1", ir->userint1, 0)ir->userint1 = get_eint(&ninp, &inp, "userint1", 0
, wi)
;
2225 ITYPE ("userint2", ir->userint2, 0)ir->userint2 = get_eint(&ninp, &inp, "userint2", 0
, wi)
;
2226 ITYPE ("userint3", ir->userint3, 0)ir->userint3 = get_eint(&ninp, &inp, "userint3", 0
, wi)
;
2227 ITYPE ("userint4", ir->userint4, 0)ir->userint4 = get_eint(&ninp, &inp, "userint4", 0
, wi)
;
2228 RTYPE ("userreal1", ir->userreal1, 0)ir->userreal1 = get_ereal(&ninp, &inp, "userreal1"
, 0, wi)
;
2229 RTYPE ("userreal2", ir->userreal2, 0)ir->userreal2 = get_ereal(&ninp, &inp, "userreal2"
, 0, wi)
;
2230 RTYPE ("userreal3", ir->userreal3, 0)ir->userreal3 = get_ereal(&ninp, &inp, "userreal3"
, 0, wi)
;
2231 RTYPE ("userreal4", ir->userreal4, 0)ir->userreal4 = get_ereal(&ninp, &inp, "userreal4"
, 0, wi)
;
2232#undef CTYPE
2233
2234 write_inpfile(mdparout, ninp, inp, FALSE0, wi);
2235 for (i = 0; (i < ninp); i++)
2236 {
2237 sfree(inp[i].name)save_free("inp[i].name", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2237, (inp[i].name))
;
2238 sfree(inp[i].value)save_free("inp[i].value", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2238, (inp[i].value))
;
2239 }
2240 sfree(inp)save_free("inp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2240, (inp))
;
2241
2242 /* Process options if necessary */
2243 for (m = 0; m < 2; m++)
2244 {
2245 for (i = 0; i < 2*DIM3; i++)
2246 {
2247 dumdub[m][i] = 0.0;
2248 }
2249 if (ir->epc)
2250 {
2251 switch (ir->epct)
2252 {
2253 case epctISOTROPIC:
2254 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX0])) != 1)
2255 {
2256 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2257 }
2258 dumdub[m][YY1] = dumdub[m][ZZ2] = dumdub[m][XX0];
2259 break;
2260 case epctSEMIISOTROPIC:
2261 case epctSURFACETENSION:
2262 if (sscanf(dumstr[m], "%lf%lf",
2263 &(dumdub[m][XX0]), &(dumdub[m][ZZ2])) != 2)
2264 {
2265 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2266 }
2267 dumdub[m][YY1] = dumdub[m][XX0];
2268 break;
2269 case epctANISOTROPIC:
2270 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2271 &(dumdub[m][XX0]), &(dumdub[m][YY1]), &(dumdub[m][ZZ2]),
2272 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2273 {
2274 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2275 }
2276 break;
2277 default:
2278 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2278
, "Pressure coupling type %s not implemented yet",
2279 epcoupltype_names[ir->epct]);
2280 }
2281 }
2282 }
2283 clear_mat(ir->ref_p);
2284 clear_mat(ir->compress);
2285 for (i = 0; i < DIM3; i++)
2286 {
2287 ir->ref_p[i][i] = dumdub[1][i];
2288 ir->compress[i][i] = dumdub[0][i];
2289 }
2290 if (ir->epct == epctANISOTROPIC)
2291 {
2292 ir->ref_p[XX0][YY1] = dumdub[1][3];
2293 ir->ref_p[XX0][ZZ2] = dumdub[1][4];
2294 ir->ref_p[YY1][ZZ2] = dumdub[1][5];
2295 if (ir->ref_p[XX0][YY1] != 0 && ir->ref_p[XX0][ZZ2] != 0 && ir->ref_p[YY1][ZZ2] != 0)
2296 {
2297 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2298 }
2299 ir->compress[XX0][YY1] = dumdub[0][3];
2300 ir->compress[XX0][ZZ2] = dumdub[0][4];
2301 ir->compress[YY1][ZZ2] = dumdub[0][5];
2302 for (i = 0; i < DIM3; i++)
2303 {
2304 for (m = 0; m < i; m++)
2305 {
2306 ir->ref_p[i][m] = ir->ref_p[m][i];
2307 ir->compress[i][m] = ir->compress[m][i];
2308 }
2309 }
2310 }
2311
2312 if (ir->comm_mode == ecmNO)
2313 {
2314 ir->nstcomm = 0;
2315 }
2316
2317 opts->couple_moltype = NULL((void*)0);
2318 if (strlen(is->couple_moltype) > 0)
2319 {
2320 if (ir->efep != efepNO)
2321 {
2322 opts->couple_moltype = strdup(is->couple_moltype)(__extension__ (__builtin_constant_p (is->couple_moltype) &&
((size_t)(const void *)((is->couple_moltype) + 1) - (size_t
)(const void *)(is->couple_moltype) == 1) ? (((const char *
) (is->couple_moltype))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (is->couple_moltype
) + 1; char *__retval = (char *) malloc (__len); if (__retval
!= ((void*)0)) __retval = (char *) memcpy (__retval, is->
couple_moltype, __len); __retval; })) : __strdup (is->couple_moltype
)))
;
2323 if (opts->couple_lam0 == opts->couple_lam1)
2324 {
2325 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2326 }
2327 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2328 opts->couple_lam1 == ecouplamNONE))
2329 {
2330 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2331 }
2332 }
2333 else
2334 {
2335 warning(wi, "Can not couple a molecule with free_energy = no");
2336 }
2337 }
2338 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2339 if (ir->efep != efepNO)
2340 {
2341 if (fep->delta_lambda > 0)
2342 {
2343 ir->efep = efepSLOWGROWTH;
2344 }
2345 }
2346
2347 if (ir->bSimTemp)
2348 {
2349 fep->bPrintEnergy = TRUE1;
2350 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2351 if the temperature is changing. */
2352 }
2353
2354 if ((ir->efep != efepNO) || ir->bSimTemp)
2355 {
2356 ir->bExpanded = FALSE0;
2357 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2358 {
2359 ir->bExpanded = TRUE1;
2360 }
2361 do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2362 if (ir->bSimTemp) /* done after fep params */
2363 {
2364 do_simtemp_params(ir);
2365 }
2366 }
2367 else
2368 {
2369 ir->fepvals->n_lambda = 0;
2370 }
2371
2372 /* WALL PARAMETERS */
2373
2374 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2375
2376 /* ORIENTATION RESTRAINT PARAMETERS */
2377
2378 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR254, NULL((void*)0)) != 1)
2379 {
2380 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2381 }
2382
2383 /* DEFORMATION PARAMETERS */
2384
2385 clear_mat(ir->deform);
2386 for (i = 0; i < 6; i++)
2387 {
2388 dumdub[0][i] = 0;
2389 }
2390 m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2391 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2392 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2393 for (i = 0; i < 3; i++)
2394 {
2395 ir->deform[i][i] = dumdub[0][i];
2396 }
2397 ir->deform[YY1][XX0] = dumdub[0][3];
2398 ir->deform[ZZ2][XX0] = dumdub[0][4];
2399 ir->deform[ZZ2][YY1] = dumdub[0][5];
2400 if (ir->epc != epcNO)
2401 {
2402 for (i = 0; i < 3; i++)
2403 {
2404 for (j = 0; j <= i; j++)
2405 {
2406 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2407 {
2408 warning_error(wi, "A box element has deform set and compressibility > 0");
2409 }
2410 }
2411 }
2412 for (i = 0; i < 3; i++)
2413 {
2414 for (j = 0; j < i; j++)
2415 {
2416 if (ir->deform[i][j] != 0)
2417 {
2418 for (m = j; m < DIM3; m++)
2419 {
2420 if (ir->compress[m][j] != 0)
2421 {
2422 sprintf(warn_buf, "An off-diagonal box element has deform set while compressibility > 0 for the same component of another box vector, this might lead to spurious periodicity effects.");
2423 warning(wi, warn_buf);
2424 }
2425 }
2426 }
2427 }
2428 }
2429 }
2430
2431 /* Ion/water position swapping checks */
2432 if (ir->eSwapCoords != eswapNO)
2433 {
2434 if (ir->swap->nstswap < 1)
2435 {
2436 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2437 }
2438 if (ir->swap->nAverage < 1)
2439 {
2440 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2441 }
2442 if (ir->swap->threshold < 1.0)
2443 {
2444 warning_error(wi, "Ion count threshold must be at least 1.\n");
2445 }
2446 }
2447
2448 sfree(dumstr[0])save_free("dumstr[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2448, (dumstr[0]))
;
2449 sfree(dumstr[1])save_free("dumstr[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2449, (dumstr[1]))
;
2450}
2451
2452static int search_QMstring(const char *s, int ng, const char *gn[])
2453{
2454 /* same as normal search_string, but this one searches QM strings */
2455 int i;
2456
2457 for (i = 0; (i < ng); i++)
2458 {
2459 if (gmx_strcasecmp(s, gn[i]) == 0)
2460 {
2461 return i;
2462 }
2463 }
2464
2465 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2465
, "this QM method or basisset (%s) is not implemented\n!", s);
2466
2467 return -1;
2468
2469} /* search_QMstring */
2470
2471/* We would like gn to be const as well, but C doesn't allow this */
2472int search_string(const char *s, int ng, char *gn[])
2473{
2474 int i;
2475
2476 for (i = 0; (i < ng); i++)
2477 {
2478 if (gmx_strcasecmp(s, gn[i]) == 0)
2479 {
2480 return i;
2481 }
2482 }
2483
2484 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2484
,
2485 "Group %s referenced in the .mdp file was not found in the index file.\n"
2486 "Group names must match either [moleculetype] names or custom index group\n"
2487 "names, in which case you must supply an index file to the '-n' option\n"
2488 "of grompp.",
2489 s);
2490
2491 return -1;
2492}
2493
2494static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2495 t_blocka *block, char *gnames[],
2496 int gtype, int restnm,
2497 int grptp, gmx_bool bVerbose,
2498 warninp_t wi)
2499{
2500 unsigned short *cbuf;
2501 t_grps *grps = &(groups->grps[gtype]);
2502 int i, j, gid, aj, ognr, ntot = 0;
2503 const char *title;
2504 gmx_bool bRest;
2505 char warn_buf[STRLEN4096];
2506
2507 if (debug)
2508 {
2509 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2510 }
2511
2512 title = gtypes[gtype];
2513
2514 snew(cbuf, natoms)(cbuf) = save_calloc("cbuf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2514, (natoms), sizeof(*(cbuf)))
;
2515 /* Mark all id's as not set */
2516 for (i = 0; (i < natoms); i++)
2517 {
2518 cbuf[i] = NOGID255;
2519 }
2520
2521 snew(grps->nm_ind, ng+1)(grps->nm_ind) = save_calloc("grps->nm_ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2521, (ng+1), sizeof(*(grps->nm_ind)))
; /* +1 for possible rest group */
2522 for (i = 0; (i < ng); i++)
2523 {
2524 /* Lookup the group name in the block structure */
2525 gid = search_string(ptrs[i], block->nr, gnames);
2526 if ((grptp != egrptpONE) || (i == 0))
2527 {
2528 grps->nm_ind[grps->nr++] = gid;
2529 }
2530 if (debug)
2531 {
2532 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2533 }
2534
2535 /* Now go over the atoms in the group */
2536 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2537 {
2538
2539 aj = block->a[j];
2540
2541 /* Range checking */
2542 if ((aj < 0) || (aj >= natoms))
2543 {
2544 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2544
, "Invalid atom number %d in indexfile", aj);
2545 }
2546 /* Lookup up the old group number */
2547 ognr = cbuf[aj];
2548 if (ognr != NOGID255)
2549 {
2550 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2550
, "Atom %d in multiple %s groups (%d and %d)",
2551 aj+1, title, ognr+1, i+1);
2552 }
2553 else
2554 {
2555 /* Store the group number in buffer */
2556 if (grptp == egrptpONE)
2557 {
2558 cbuf[aj] = 0;
2559 }
2560 else
2561 {
2562 cbuf[aj] = i;
2563 }
2564 ntot++;
2565 }
2566 }
2567 }
2568
2569 /* Now check whether we have done all atoms */
2570 bRest = FALSE0;
2571 if (ntot != natoms)
2572 {
2573 if (grptp == egrptpALL)
2574 {
2575 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2575
, "%d atoms are not part of any of the %s groups",
2576 natoms-ntot, title);
2577 }
2578 else if (grptp == egrptpPART)
2579 {
2580 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2581 natoms-ntot, title);
2582 warning_note(wi, warn_buf);
2583 }
2584 /* Assign all atoms currently unassigned to a rest group */
2585 for (j = 0; (j < natoms); j++)
2586 {
2587 if (cbuf[j] == NOGID255)
2588 {
2589 cbuf[j] = grps->nr;
2590 bRest = TRUE1;
2591 }
2592 }
2593 if (grptp != egrptpPART)
2594 {
2595 if (bVerbose)
2596 {
2597 fprintf(stderrstderr,
2598 "Making dummy/rest group for %s containing %d elements\n",
2599 title, natoms-ntot);
2600 }
2601 /* Add group name "rest" */
2602 grps->nm_ind[grps->nr] = restnm;
2603
2604 /* Assign the rest name to all atoms not currently assigned to a group */
2605 for (j = 0; (j < natoms); j++)
2606 {
2607 if (cbuf[j] == NOGID255)
2608 {
2609 cbuf[j] = grps->nr;
2610 }
2611 }
2612 grps->nr++;
2613 }
2614 }
2615
2616 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2617 {
2618 /* All atoms are part of one (or no) group, no index required */
2619 groups->ngrpnr[gtype] = 0;
2620 groups->grpnr[gtype] = NULL((void*)0);
2621 }
2622 else
2623 {
2624 groups->ngrpnr[gtype] = natoms;
2625 snew(groups->grpnr[gtype], natoms)(groups->grpnr[gtype]) = save_calloc("groups->grpnr[gtype]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2625, (natoms), sizeof(*(groups->grpnr[gtype])))
;
2626 for (j = 0; (j < natoms); j++)
2627 {
2628 groups->grpnr[gtype][j] = cbuf[j];
2629 }
2630 }
2631
2632 sfree(cbuf)save_free("cbuf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2632, (cbuf))
;
2633
2634 return (bRest && grptp == egrptpPART);
2635}
2636
2637static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2638{
2639 t_grpopts *opts;
2640 gmx_groups_t *groups;
2641 t_pull *pull;
2642 int natoms, ai, aj, i, j, d, g, imin, jmin;
2643 t_iatom *ia;
2644 int *nrdf2, *na_vcm, na_tot;
2645 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2646 gmx_mtop_atomloop_all_t aloop;
2647 t_atom *atom;
2648 int mb, mol, ftype, as;
2649 gmx_molblock_t *molb;
2650 gmx_moltype_t *molt;
2651
2652 /* Calculate nrdf.
2653 * First calc 3xnr-atoms for each group
2654 * then subtract half a degree of freedom for each constraint
2655 *
2656 * Only atoms and nuclei contribute to the degrees of freedom...
2657 */
2658
2659 opts = &ir->opts;
2660
2661 groups = &mtop->groups;
2662 natoms = mtop->natoms;
2663
2664 /* Allocate one more for a possible rest group */
2665 /* We need to sum degrees of freedom into doubles,
2666 * since floats give too low nrdf's above 3 million atoms.
2667 */
2668 snew(nrdf_tc, groups->grps[egcTC].nr+1)(nrdf_tc) = save_calloc("nrdf_tc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2668, (groups->grps[egcTC].nr+1), sizeof(*(nrdf_tc)))
;
2669 snew(nrdf_vcm, groups->grps[egcVCM].nr+1)(nrdf_vcm) = save_calloc("nrdf_vcm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2669, (groups->grps[egcVCM].nr+1), sizeof(*(nrdf_vcm)))
;
2670 snew(na_vcm, groups->grps[egcVCM].nr+1)(na_vcm) = save_calloc("na_vcm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2670, (groups->grps[egcVCM].nr+1), sizeof(*(na_vcm)))
;
2671
2672 for (i = 0; i < groups->grps[egcTC].nr; i++)
2673 {
2674 nrdf_tc[i] = 0;
2675 }
2676 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2677 {
2678 nrdf_vcm[i] = 0;
2679 }
2680
2681 snew(nrdf2, natoms)(nrdf2) = save_calloc("nrdf2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2681, (natoms), sizeof(*(nrdf2)))
;
2682 aloop = gmx_mtop_atomloop_all_init(mtop);
2683 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2684 {
2685 nrdf2[i] = 0;
2686 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2687 {
2688 g = ggrpnr(groups, egcFREEZE, i)((groups)->grpnr[egcFREEZE] ? (groups)->grpnr[egcFREEZE
][i] : 0)
;
2689 /* Double count nrdf for particle i */
2690 for (d = 0; d < DIM3; d++)
2691 {
2692 if (opts->nFreeze[g][d] == 0)
2693 {
2694 nrdf2[i] += 2;
2695 }
2696 }
2697 nrdf_tc [ggrpnr(groups, egcTC, i)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][i] : 0
)
] += 0.5*nrdf2[i];
2698 nrdf_vcm[ggrpnr(groups, egcVCM, i)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][i] :
0)
] += 0.5*nrdf2[i];
2699 }
2700 }
2701
2702 as = 0;
2703 for (mb = 0; mb < mtop->nmolblock; mb++)
2704 {
2705 molb = &mtop->molblock[mb];
2706 molt = &mtop->moltype[molb->type];
2707 atom = molt->atoms.atom;
2708 for (mol = 0; mol < molb->nmol; mol++)
2709 {
2710 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2711 {
2712 ia = molt->ilist[ftype].iatoms;
2713 for (i = 0; i < molt->ilist[ftype].nr; )
2714 {
2715 /* Subtract degrees of freedom for the constraints,
2716 * if the particles still have degrees of freedom left.
2717 * If one of the particles is a vsite or a shell, then all
2718 * constraint motion will go there, but since they do not
2719 * contribute to the constraints the degrees of freedom do not
2720 * change.
2721 */
2722 ai = as + ia[1];
2723 aj = as + ia[2];
2724 if (((atom[ia[1]].ptype == eptNucleus) ||
2725 (atom[ia[1]].ptype == eptAtom)) &&
2726 ((atom[ia[2]].ptype == eptNucleus) ||
2727 (atom[ia[2]].ptype == eptAtom)))
2728 {
2729 if (nrdf2[ai] > 0)
2730 {
2731 jmin = 1;
2732 }
2733 else
2734 {
2735 jmin = 2;
2736 }
2737 if (nrdf2[aj] > 0)
2738 {
2739 imin = 1;
2740 }
2741 else
2742 {
2743 imin = 2;
2744 }
2745 imin = min(imin, nrdf2[ai])(((imin) < (nrdf2[ai])) ? (imin) : (nrdf2[ai]) );
2746 jmin = min(jmin, nrdf2[aj])(((jmin) < (nrdf2[aj])) ? (jmin) : (nrdf2[aj]) );
2747 nrdf2[ai] -= imin;
2748 nrdf2[aj] -= jmin;
2749 nrdf_tc [ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
] -= 0.5*imin;
2750 nrdf_tc [ggrpnr(groups, egcTC, aj)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][aj] : 0
)
] -= 0.5*jmin;
2751 nrdf_vcm[ggrpnr(groups, egcVCM, ai)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][ai] :
0)
] -= 0.5*imin;
2752 nrdf_vcm[ggrpnr(groups, egcVCM, aj)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][aj] :
0)
] -= 0.5*jmin;
2753 }
2754 ia += interaction_function[ftype].nratoms+1;
2755 i += interaction_function[ftype].nratoms+1;
2756 }
2757 }
2758 ia = molt->ilist[F_SETTLE].iatoms;
2759 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2760 {
2761 /* Subtract 1 dof from every atom in the SETTLE */
2762 for (j = 0; j < 3; j++)
2763 {
2764 ai = as + ia[1+j];
2765 imin = min(2, nrdf2[ai])(((2) < (nrdf2[ai])) ? (2) : (nrdf2[ai]) );
2766 nrdf2[ai] -= imin;
2767 nrdf_tc [ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
] -= 0.5*imin;
2768 nrdf_vcm[ggrpnr(groups, egcVCM, ai)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][ai] :
0)
] -= 0.5*imin;
2769 }
2770 ia += 4;
2771 i += 4;
2772 }
2773 as += molt->atoms.nr;
2774 }
2775 }
2776
2777 if (ir->ePull == epullCONSTRAINT)
2778 {
2779 /* Correct nrdf for the COM constraints.
2780 * We correct using the TC and VCM group of the first atom
2781 * in the reference and pull group. If atoms in one pull group
2782 * belong to different TC or VCM groups it is anyhow difficult
2783 * to determine the optimal nrdf assignment.
2784 */
2785 pull = ir->pull;
2786
2787 for (i = 0; i < pull->ncoord; i++)
2788 {
2789 imin = 1;
2790
2791 for (j = 0; j < 2; j++)
2792 {
2793 const t_pull_group *pgrp;
2794
2795 pgrp = &pull->group[pull->coord[i].group[j]];
2796
2797 if (pgrp->nat > 0)
2798 {
2799 /* Subtract 1/2 dof from each group */
2800 ai = pgrp->ind[0];
2801 nrdf_tc [ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
] -= 0.5*imin;
2802 nrdf_vcm[ggrpnr(groups, egcVCM, ai)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][ai] :
0)
] -= 0.5*imin;
2803 if (nrdf_tc[ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
] < 0)
2804 {
2805 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2805
, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
]]);
2806 }
2807 }
2808 else
2809 {
2810 /* We need to subtract the whole DOF from group j=1 */
2811 imin += 1;
2812 }
2813 }
2814 }
2815 }
2816
2817 if (ir->nstcomm != 0)
2818 {
2819 /* Subtract 3 from the number of degrees of freedom in each vcm group
2820 * when com translation is removed and 6 when rotation is removed
2821 * as well.
2822 */
2823 switch (ir->comm_mode)
2824 {
2825 case ecmLINEAR:
2826 n_sub = ndof_com(ir);
2827 break;
2828 case ecmANGULAR:
2829 n_sub = 6;
2830 break;
2831 default:
2832 n_sub = 0;
2833 gmx_incons("Checking comm_mode")_gmx_error("incons", "Checking comm_mode", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2833)
;
2834 }
2835
2836 for (i = 0; i < groups->grps[egcTC].nr; i++)
2837 {
2838 /* Count the number of atoms of TC group i for every VCM group */
2839 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2840 {
2841 na_vcm[j] = 0;
2842 }
2843 na_tot = 0;
2844 for (ai = 0; ai < natoms; ai++)
2845 {
2846 if (ggrpnr(groups, egcTC, ai)((groups)->grpnr[egcTC] ? (groups)->grpnr[egcTC][ai] : 0
)
== i)
2847 {
2848 na_vcm[ggrpnr(groups, egcVCM, ai)((groups)->grpnr[egcVCM] ? (groups)->grpnr[egcVCM][ai] :
0)
]++;
2849 na_tot++;
2850 }
2851 }
2852 /* Correct for VCM removal according to the fraction of each VCM
2853 * group present in this TC group.
2854 */
2855 nrdf_uc = nrdf_tc[i];
2856 if (debug)
2857 {
2858 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2859 i, nrdf_uc, n_sub);
2860 }
2861 nrdf_tc[i] = 0;
2862 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2863 {
2864 if (nrdf_vcm[j] > n_sub)
2865 {
2866 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2867 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2868 }
2869 if (debug)
2870 {
2871 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2872 j, nrdf_vcm[j], nrdf_tc[i]);
2873 }
2874 }
2875 }
2876 }
2877 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2878 {
2879 opts->nrdf[i] = nrdf_tc[i];
2880 if (opts->nrdf[i] < 0)
2881 {
2882 opts->nrdf[i] = 0;
2883 }
2884 fprintf(stderrstderr,
2885 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2886 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2887 }
2888
2889 sfree(nrdf2)save_free("nrdf2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2889, (nrdf2))
;
2890 sfree(nrdf_tc)save_free("nrdf_tc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2890, (nrdf_tc))
;
2891 sfree(nrdf_vcm)save_free("nrdf_vcm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2891, (nrdf_vcm))
;
2892 sfree(na_vcm)save_free("na_vcm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2892, (na_vcm))
;
2893}
2894
2895static void decode_cos(char *s, t_cosines *cosine)
2896{
2897 char *t;
2898 char format[STRLEN4096], f1[STRLEN4096];
2899 double a, phi;
2900 int i;
2901
2902 t = strdup(s)(__extension__ (__builtin_constant_p (s) && ((size_t)
(const void *)((s) + 1) - (size_t)(const void *)(s) == 1) ? (
((const char *) (s))[0] == '\0' ? (char *) calloc ((size_t) 1
, (size_t) 1) : ({ size_t __len = strlen (s) + 1; char *__retval
= (char *) malloc (__len); if (__retval != ((void*)0)) __retval
= (char *) memcpy (__retval, s, __len); __retval; })) : __strdup
(s)))
;
2903 trim(t);
2904
2905 cosine->n = 0;
2906 cosine->a = NULL((void*)0);
2907 cosine->phi = NULL((void*)0);
2908 if (strlen(t))
2909 {
2910 sscanf(t, "%d", &(cosine->n));
2911 if (cosine->n <= 0)
2912 {
2913 cosine->n = 0;
2914 }
2915 else
2916 {
2917 snew(cosine->a, cosine->n)(cosine->a) = save_calloc("cosine->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2917, (cosine->n), sizeof(*(cosine->a)))
;
2918 snew(cosine->phi, cosine->n)(cosine->phi) = save_calloc("cosine->phi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2918, (cosine->n), sizeof(*(cosine->phi)))
;
2919
2920 sprintf(format, "%%*d");
2921 for (i = 0; (i < cosine->n); i++)
2922 {
2923 strcpy(f1, format);
2924 strcat(f1, "%lf%lf");
2925 if (sscanf(t, f1, &a, &phi) < 2)
2926 {
2927 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2927
, "Invalid input for electric field shift: '%s'", t);
2928 }
2929 cosine->a[i] = a;
2930 cosine->phi[i] = phi;
2931 strcat(format, "%*lf%*lf");
2932 }
2933 }
2934 }
2935 sfree(t)save_free("t", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2935, (t))
;
2936}
2937
2938static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2939 const char *option, const char *val, int flag)
2940{
2941 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2942 * But since this is much larger than STRLEN, such a line can not be parsed.
2943 * The real maximum is the number of names that fit in a string: STRLEN/2.
2944 */
2945#define EGP_MAX(4096/2) (STRLEN4096/2)
2946 int nelem, i, j, k, nr;
2947 char *names[EGP_MAX(4096/2)];
2948 char ***gnames;
2949 gmx_bool bSet;
2950
2951 gnames = groups->grpname;
2952
2953 nelem = str_nelem(val, EGP_MAX(4096/2), names);
2954 if (nelem % 2 != 0)
2955 {
2956 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2956
, "The number of groups for %s is odd", option);
2957 }
2958 nr = groups->grps[egcENER].nr;
2959 bSet = FALSE0;
2960 for (i = 0; i < nelem/2; i++)
2961 {
2962 j = 0;
2963 while ((j < nr) &&
2964 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2965 {
2966 j++;
2967 }
2968 if (j == nr)
2969 {
2970 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2970
, "%s in %s is not an energy group\n",
2971 names[2*i], option);
2972 }
2973 k = 0;
2974 while ((k < nr) &&
2975 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2976 {
2977 k++;
2978 }
2979 if (k == nr)
2980 {
2981 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 2981
, "%s in %s is not an energy group\n",
2982 names[2*i+1], option);
2983 }
2984 if ((j < nr) && (k < nr))
2985 {
2986 ir->opts.egp_flags[nr*j+k] |= flag;
2987 ir->opts.egp_flags[nr*k+j] |= flag;
2988 bSet = TRUE1;
2989 }
2990 }
2991
2992 return bSet;
2993}
2994
2995
2996static void make_swap_groups(
2997 t_swapcoords *swap,
2998 char *swapgname,
2999 char *splitg0name,
3000 char *splitg1name,
3001 char *solgname,
3002 t_blocka *grps,
3003 char **gnames)
3004{
3005 int ig = -1, i = 0, j;
3006 char *splitg;
3007
3008
3009 /* Just a quick check here, more thorough checks are in mdrun */
3010 if (strcmp(splitg0name, splitg1name)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(splitg0name) && __builtin_constant_p (splitg1name) &&
(__s1_len = strlen (splitg0name), __s2_len = strlen (splitg1name
), (!((size_t)(const void *)((splitg0name) + 1) - (size_t)(const
void *)(splitg0name) == 1) || __s1_len >= 4) && (
!((size_t)(const void *)((splitg1name) + 1) - (size_t)(const void
*)(splitg1name) == 1) || __s2_len >= 4)) ? __builtin_strcmp
(splitg0name, splitg1name) : (__builtin_constant_p (splitg0name
) && ((size_t)(const void *)((splitg0name) + 1) - (size_t
)(const void *)(splitg0name) == 1) && (__s1_len = strlen
(splitg0name), __s1_len < 4) ? (__builtin_constant_p (splitg1name
) && ((size_t)(const void *)((splitg1name) + 1) - (size_t
)(const void *)(splitg1name) == 1) ? __builtin_strcmp (splitg0name
, splitg1name) : (__extension__ ({ const unsigned char *__s2 =
(const unsigned char *) (const char *) (splitg1name); int __result
= (((const unsigned char *) (const char *) (splitg0name))[0]
- __s2[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (splitg0name))[1]
- __s2[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (splitg0name))[2]
- __s2[2]); if (__s1_len > 2 && __result == 0) __result
= (((const unsigned char *) (const char *) (splitg0name))[3]
- __s2[3]); } } __result; }))) : (__builtin_constant_p (splitg1name
) && ((size_t)(const void *)((splitg1name) + 1) - (size_t
)(const void *)(splitg1name) == 1) && (__s2_len = strlen
(splitg1name), __s2_len < 4) ? (__builtin_constant_p (splitg0name
) && ((size_t)(const void *)((splitg0name) + 1) - (size_t
)(const void *)(splitg0name) == 1) ? __builtin_strcmp (splitg0name
, splitg1name) : (- (__extension__ ({ const unsigned char *__s2
= (const unsigned char *) (const char *) (splitg0name); int __result
= (((const unsigned char *) (const char *) (splitg1name))[0]
- __s2[0]); if (__s2_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (splitg1name))[1]
- __s2[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (splitg1name))[2]
- __s2[2]); if (__s2_len > 2 && __result == 0) __result
= (((const unsigned char *) (const char *) (splitg1name))[3]
- __s2[3]); } } __result; })))) : __builtin_strcmp (splitg0name
, splitg1name)))); })
== 0)
3011 {
3012 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3012
, "The split groups can not both be '%s'.", splitg0name);
3013 }
3014
3015 /* First get the swap group index atoms */
3016 ig = search_string(swapgname, grps->nr, gnames);
3017 swap->nat = grps->index[ig+1] - grps->index[ig];
3018 if (swap->nat > 0)
3019 {
3020 fprintf(stderrstderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
3021 snew(swap->ind, swap->nat)(swap->ind) = save_calloc("swap->ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3021, (swap->nat), sizeof(*(swap->ind)))
;
3022 for (i = 0; i < swap->nat; i++)
3023 {
3024 swap->ind[i] = grps->a[grps->index[ig]+i];
3025 }
3026 }
3027 else
3028 {
3029 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3029
, "You defined an empty group of atoms for swapping.");
3030 }
3031
3032 /* Now do so for the split groups */
3033 for (j = 0; j < 2; j++)
3034 {
3035 if (j == 0)
3036 {
3037 splitg = splitg0name;
3038 }
3039 else
3040 {
3041 splitg = splitg1name;
3042 }
3043
3044 ig = search_string(splitg, grps->nr, gnames);
3045 swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3046 if (swap->nat_split[j] > 0)
3047 {
3048 fprintf(stderrstderr, "Split group %d '%s' contains %d atom%s.\n",
3049 j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3050 snew(swap->ind_split[j], swap->nat_split[j])(swap->ind_split[j]) = save_calloc("swap->ind_split[j]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3050, (swap->nat_split[j]), sizeof(*(swap->ind_split[
j])))
;
3051 for (i = 0; i < swap->nat_split[j]; i++)
3052 {
3053 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3054 }
3055 }
3056 else
3057 {
3058 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3058
, "Split group %d has to contain at least 1 atom!", j);
3059 }
3060 }
3061
3062 /* Now get the solvent group index atoms */
3063 ig = search_string(solgname, grps->nr, gnames);
3064 swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3065 if (swap->nat_sol > 0)
3066 {
3067 fprintf(stderrstderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3068 snew(swap->ind_sol, swap->nat_sol)(swap->ind_sol) = save_calloc("swap->ind_sol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3068, (swap->nat_sol), sizeof(*(swap->ind_sol)))
;
3069 for (i = 0; i < swap->nat_sol; i++)
3070 {
3071 swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3072 }
3073 }
3074 else
3075 {
3076 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3076
, "You defined an empty group of solvent. Cannot exchange ions.");
3077 }
3078}
3079
3080
3081void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3082{
3083 int ig = -1, i;
3084
3085
3086 ig = search_string(IMDgname, grps->nr, gnames);
3087 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3088
3089 if (IMDgroup->nat > 0)
3090 {
3091 fprintf(stderrstderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3092 IMDgname, IMDgroup->nat);
3093 snew(IMDgroup->ind, IMDgroup->nat)(IMDgroup->ind) = save_calloc("IMDgroup->ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3093, (IMDgroup->nat), sizeof(*(IMDgroup->ind)))
;
3094 for (i = 0; i < IMDgroup->nat; i++)
3095 {
3096 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3097 }
3098 }
3099}
3100
3101
3102void do_index(const char* mdparin, const char *ndx,
3103 gmx_mtop_t *mtop,
3104 gmx_bool bVerbose,
3105 t_inputrec *ir, rvec *v,
3106 warninp_t wi)
3107{
3108 t_blocka *grps;
3109 gmx_groups_t *groups;
3110 int natoms;
3111 t_symtab *symtab;
3112 t_atoms atoms_all;
3113 char warnbuf[STRLEN4096], **gnames;
3114 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3115 real tau_min;
3116 int nstcmin;
3117 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3118 char *ptr1[MAXPTR254], *ptr2[MAXPTR254], *ptr3[MAXPTR254];
3119 int i, j, k, restnm;
3120 real SAtime;
3121 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3122 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3123 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3124 char warn_buf[STRLEN4096];
3125
3126 if (bVerbose)
3127 {
3128 fprintf(stderrstderr, "processing index file...\n");
3129 }
3130 debug_gmx();
3131 if (ndx == NULL((void*)0))
3132 {
3133 snew(grps, 1)(grps) = save_calloc("grps", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3133, (1), sizeof(*(grps)))
;
3134 snew(grps->index, 1)(grps->index) = save_calloc("grps->index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3134, (1), sizeof(*(grps->index)))
;
3135 snew(gnames, 1)(gnames) = save_calloc("gnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3135, (1), sizeof(*(gnames)))
;
3136 atoms_all = gmx_mtop_global_atoms(mtop);
3137 analyse(&atoms_all, grps, &gnames, FALSE0, TRUE1);
3138 free_t_atoms(&atoms_all, FALSE0);
3139 }
3140 else
3141 {
3142 grps = init_index(ndx, &gnames);
3143 }
3144
3145 groups = &mtop->groups;
3146 natoms = mtop->natoms;
3147 symtab = &mtop->symtab;
3148
3149 snew(groups->grpname, grps->nr+1)(groups->grpname) = save_calloc("groups->grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3149, (grps->nr+1), sizeof(*(groups->grpname)))
;
3150
3151 for (i = 0; (i < grps->nr); i++)
3152 {
3153 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3154 }
3155 groups->grpname[i] = put_symtab(symtab, "rest");
3156 restnm = i;
3157 srenew(gnames, grps->nr+1)(gnames) = save_realloc("gnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3157, (gnames), (grps->nr+1), sizeof(*(gnames)))
;
3158 gnames[restnm] = *(groups->grpname[i]);
3159 groups->ngrpname = grps->nr+1;
3160
3161 set_warning_line(wi, mdparin, -1);
3162
3163 ntau_t = str_nelem(is->tau_t, MAXPTR254, ptr1);
3164 nref_t = str_nelem(is->ref_t, MAXPTR254, ptr2);
3165 ntcg = str_nelem(is->tcgrps, MAXPTR254, ptr3);
3166 if ((ntau_t != ntcg) || (nref_t != ntcg))
3167 {
3168 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3168
, "Invalid T coupling input: %d groups, %d ref-t values and "
3169 "%d tau-t values", ntcg, nref_t, ntau_t);
3170 }
3171
3172 bSetTCpar = (ir->etc || EI_SD(ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || ir->eI == eiBD || EI_TPI(ir->eI)((ir->eI) == eiTPI || (ir->eI) == eiTPIC));
3173 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3174 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3175 nr = groups->grps[egcTC].nr;
3176 ir->opts.ngtc = nr;
3177 snew(ir->opts.nrdf, nr)(ir->opts.nrdf) = save_calloc("ir->opts.nrdf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3177, (nr), sizeof(*(ir->opts.nrdf)))
;
3178 snew(ir->opts.tau_t, nr)(ir->opts.tau_t) = save_calloc("ir->opts.tau_t", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3178, (nr), sizeof(*(ir->opts.tau_t)))
;
3179 snew(ir->opts.ref_t, nr)(ir->opts.ref_t) = save_calloc("ir->opts.ref_t", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3179, (nr), sizeof(*(ir->opts.ref_t)))
;
3180 if (ir->eI == eiBD && ir->bd_fric == 0)
3181 {
3182 fprintf(stderrstderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3183 }
3184
3185 if (bSetTCpar)
3186 {
3187 if (nr != nref_t)
3188 {
3189 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3189
, "Not enough ref-t and tau-t values!");
3190 }
3191
3192 tau_min = 1e20;
3193 for (i = 0; (i < nr); i++)
3194 {
3195 ir->opts.tau_t[i] = strtod(ptr1[i], NULL((void*)0));
3196 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3197 {
3198 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3199 warning_error(wi, warn_buf);
3200 }
3201
3202 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3203 {
3204 warning_note(wi, "tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3205 }
3206
3207 if (ir->opts.tau_t[i] >= 0)
3208 {
3209 tau_min = min(tau_min, ir->opts.tau_t[i])(((tau_min) < (ir->opts.tau_t[i])) ? (tau_min) : (ir->
opts.tau_t[i]) )
;
3210 }
3211 }
3212 if (ir->etc != etcNO && ir->nsttcouple == -1)
3213 {
3214 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3215 }
3216
3217 if (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK))
3218 {
3219 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3220 {
3221 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3221
, "Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
3222 }
3223 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3224 {
3225 if (ir->nstpcouple != ir->nsttcouple)
3226 {
3227 int mincouple = min(ir->nstpcouple, ir->nsttcouple)(((ir->nstpcouple) < (ir->nsttcouple)) ? (ir->nstpcouple
) : (ir->nsttcouple) )
;
3228 ir->nstpcouple = ir->nsttcouple = mincouple;
3229 sprintf(warn_buf, "for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d", mincouple);
3230 warning_note(wi, warn_buf);
3231 }
3232 }
3233 }
3234 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3235 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3236
3237 if (ETC_ANDERSEN(ir->etc)(((ir->etc) == etcANDERSENMASSIVE) || ((ir->etc) == etcANDERSEN
))
)
3238 {
3239 if (ir->nsttcouple != 1)
3240 {
3241 ir->nsttcouple = 1;
3242 sprintf(warn_buf, "Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
3243 warning_note(wi, warn_buf);
3244 }
3245 }
3246 nstcmin = tcouple_min_integration_steps(ir->etc);
3247 if (nstcmin > 1)
3248 {
3249 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3250 {
3251 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3252 ETCOUPLTYPE(ir->etc)((((ir->etc) < 0) || ((ir->etc) >= (etcNR))) ? "UNDEFINED"
: (etcoupl_names)[ir->etc])
,
3253 tau_min, nstcmin,
3254 ir->nsttcouple*ir->delta_t);
3255 warning(wi, warn_buf);
3256 }
3257 }
3258 for (i = 0; (i < nr); i++)
3259 {
3260 ir->opts.ref_t[i] = strtod(ptr2[i], NULL((void*)0));
3261 if (ir->opts.ref_t[i] < 0)
3262 {
3263 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3263
, "ref-t for group %d negative", i);
3264 }
3265 }
3266 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3267 if we are in this conditional) if mc_temp is negative */
3268 if (ir->expandedvals->mc_temp < 0)
3269 {
3270 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3271 }
3272 }
3273
3274 /* Simulated annealing for each group. There are nr groups */
3275 nSA = str_nelem(is->anneal, MAXPTR254, ptr1);
3276 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3277 {
3278 nSA = 0;
3279 }
3280 if (nSA > 0 && nSA != nr)
3281 {
3282 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3282
, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3283 }
3284 else
3285 {
3286 snew(ir->opts.annealing, nr)(ir->opts.annealing) = save_calloc("ir->opts.annealing"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3286, (nr), sizeof(*(ir->opts.annealing)))
;
3287 snew(ir->opts.anneal_npoints, nr)(ir->opts.anneal_npoints) = save_calloc("ir->opts.anneal_npoints"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3287, (nr), sizeof(*(ir->opts.anneal_npoints)))
;
3288 snew(ir->opts.anneal_time, nr)(ir->opts.anneal_time) = save_calloc("ir->opts.anneal_time"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3288, (nr), sizeof(*(ir->opts.anneal_time)))
;
3289 snew(ir->opts.anneal_temp, nr)(ir->opts.anneal_temp) = save_calloc("ir->opts.anneal_temp"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3289, (nr), sizeof(*(ir->opts.anneal_temp)))
;
3290 for (i = 0; i < nr; i++)
3291 {
3292 ir->opts.annealing[i] = eannNO;
3293 ir->opts.anneal_npoints[i] = 0;
3294 ir->opts.anneal_time[i] = NULL((void*)0);
3295 ir->opts.anneal_temp[i] = NULL((void*)0);
3296 }
3297 if (nSA > 0)
3298 {
3299 bAnneal = FALSE0;
3300 for (i = 0; i < nr; i++)
3301 {
3302 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3303 {
3304 ir->opts.annealing[i] = eannNO;
3305 }
3306 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3307 {
3308 ir->opts.annealing[i] = eannSINGLE;
3309 bAnneal = TRUE1;
3310 }
3311 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3312 {
3313 ir->opts.annealing[i] = eannPERIODIC;
3314 bAnneal = TRUE1;
3315 }
3316 }
3317 if (bAnneal)
3318 {
3319 /* Read the other fields too */
3320 nSA_points = str_nelem(is->anneal_npoints, MAXPTR254, ptr1);
3321 if (nSA_points != nSA)
3322 {
3323 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3323
, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3324 }
3325 for (k = 0, i = 0; i < nr; i++)
3326 {
3327 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL((void*)0), 10);
3328 if (ir->opts.anneal_npoints[i] == 1)
3329 {
3330 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3330
, "Please specify at least a start and an end point for annealing\n");
3331 }
3332 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i])(ir->opts.anneal_time[i]) = save_calloc("ir->opts.anneal_time[i]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3332, (ir->opts.anneal_npoints[i]), sizeof(*(ir->opts
.anneal_time[i])))
;
3333 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i])(ir->opts.anneal_temp[i]) = save_calloc("ir->opts.anneal_temp[i]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3333, (ir->opts.anneal_npoints[i]), sizeof(*(ir->opts
.anneal_temp[i])))
;
3334 k += ir->opts.anneal_npoints[i];
3335 }
3336
3337 nSA_time = str_nelem(is->anneal_time, MAXPTR254, ptr1);
3338 if (nSA_time != k)
3339 {
3340 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3340
, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3341 }
3342 nSA_temp = str_nelem(is->anneal_temp, MAXPTR254, ptr2);
3343 if (nSA_temp != k)
3344 {
3345 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3345
, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3346 }
3347
3348 for (i = 0, k = 0; i < nr; i++)
3349 {
3350
3351 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3352 {
3353 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL((void*)0));
3354 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL((void*)0));
3355 if (j == 0)
3356 {
3357 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS5.96046448E-08))
3358 {
3359 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3359
, "First time point for annealing > init_t.\n");
3360 }
3361 }
3362 else
3363 {
3364 /* j>0 */
3365 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3366 {
3367 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3367
, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3368 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3369 }
3370 }
3371 if (ir->opts.anneal_temp[i][j] < 0)
3372 {
3373 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3373
, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3374 }
3375 k++;
3376 }
3377 }
3378 /* Print out some summary information, to make sure we got it right */
3379 for (i = 0, k = 0; i < nr; i++)
3380 {
3381 if (ir->opts.annealing[i] != eannNO)
3382 {
3383 j = groups->grps[egcTC].nm_ind[i];
3384 fprintf(stderrstderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3385 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3386 ir->opts.anneal_npoints[i]);
3387 fprintf(stderrstderr, "Time (ps) Temperature (K)\n");
3388 /* All terms except the last one */
3389 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3390 {
3391 fprintf(stderrstderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3392 }
3393
3394 /* Finally the last one */
3395 j = ir->opts.anneal_npoints[i]-1;
3396 if (ir->opts.annealing[i] == eannSINGLE)
3397 {
3398 fprintf(stderrstderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3399 }
3400 else
3401 {
3402 fprintf(stderrstderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3403 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS5.96046448E-08)
3404 {
3405 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3406 }
3407 }
3408 }
3409 }
3410 }
3411 }
3412 }
3413
3414 if (ir->ePull != epullNO)
3415 {
3416 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3417
3418 make_pull_coords(ir->pull);
3419 }
3420
3421 if (ir->bRot)
3422 {
3423 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3424 }
3425
3426 if (ir->eSwapCoords != eswapNO)
3427 {
3428 make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3429 }
3430
3431 /* Make indices for IMD session */
3432 if (ir->bIMD)
3433 {
3434 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3435 }
3436
3437 nacc = str_nelem(is->acc, MAXPTR254, ptr1);
3438 nacg = str_nelem(is->accgrps, MAXPTR254, ptr2);
3439 if (nacg*DIM3 != nacc)
3440 {
3441 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3441
, "Invalid Acceleration input: %d groups and %d acc. values",
3442 nacg, nacc);
3443 }
3444 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3445 restnm, egrptpALL_GENREST, bVerbose, wi);
3446 nr = groups->grps[egcACC].nr;
3447 snew(ir->opts.acc, nr)(ir->opts.acc) = save_calloc("ir->opts.acc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3447, (nr), sizeof(*(ir->opts.acc)))
;
3448 ir->opts.ngacc = nr;
3449
3450 for (i = k = 0; (i < nacg); i++)
3451 {
3452 for (j = 0; (j < DIM3); j++, k++)
3453 {
3454 ir->opts.acc[i][j] = strtod(ptr1[k], NULL((void*)0));
3455 }
3456 }
3457 for (; (i < nr); i++)
3458 {
3459 for (j = 0; (j < DIM3); j++)
3460 {
3461 ir->opts.acc[i][j] = 0;
3462 }
3463 }
3464
3465 nfrdim = str_nelem(is->frdim, MAXPTR254, ptr1);
3466 nfreeze = str_nelem(is->freeze, MAXPTR254, ptr2);
3467 if (nfrdim != DIM3*nfreeze)
3468 {
3469 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3469
, "Invalid Freezing input: %d groups and %d freeze values",
3470 nfreeze, nfrdim);
3471 }
3472 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3473 restnm, egrptpALL_GENREST, bVerbose, wi);
3474 nr = groups->grps[egcFREEZE].nr;
3475 ir->opts.ngfrz = nr;
3476 snew(ir->opts.nFreeze, nr)(ir->opts.nFreeze) = save_calloc("ir->opts.nFreeze", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3476, (nr), sizeof(*(ir->opts.nFreeze)))
;
3477 for (i = k = 0; (i < nfreeze); i++)
3478 {
3479 for (j = 0; (j < DIM3); j++, k++)
3480 {
3481 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3482 if (!ir->opts.nFreeze[i][j])
3483 {
3484 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3485 {
3486 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3487 "(not %s)", ptr1[k]);
3488 warning(wi, warn_buf);
3489 }
3490 }
3491 }
3492 }
3493 for (; (i < nr); i++)
3494 {
3495 for (j = 0; (j < DIM3); j++)
3496 {
3497 ir->opts.nFreeze[i][j] = 0;
3498 }
3499 }
3500
3501 nenergy = str_nelem(is->energy, MAXPTR254, ptr1);
3502 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3503 restnm, egrptpALL_GENREST, bVerbose, wi);
3504 add_wall_energrps(groups, ir->nwall, symtab);
3505 ir->opts.ngener = groups->grps[egcENER].nr;
3506 nvcm = str_nelem(is->vcm, MAXPTR254, ptr1);
3507 bRest =
3508 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3509 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3510 if (bRest)
3511 {
3512 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3513 "This may lead to artifacts.\n"
3514 "In most cases one should use one group for the whole system.");
3515 }
3516
3517 /* Now we have filled the freeze struct, so we can calculate NRDF */
3518 calc_nrdf(mtop, ir, gnames);
3519
3520 if (v && NULL((void*)0))
3521 {
3522 real fac, ntot = 0;
3523
3524 /* Must check per group! */
3525 for (i = 0; (i < ir->opts.ngtc); i++)
3526 {
3527 ntot += ir->opts.nrdf[i];
3528 }
3529 if (ntot != (DIM3*natoms))
3530 {
3531 fac = sqrt(ntot/(DIM3*natoms));
3532 if (bVerbose)
3533 {
3534 fprintf(stderrstderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3535 "and removal of center of mass motion\n", fac);
3536 }
3537 for (i = 0; (i < natoms); i++)
3538 {
3539 svmul(fac, v[i], v[i]);
3540 }
3541 }
3542 }
3543
3544 nuser = str_nelem(is->user1, MAXPTR254, ptr1);
3545 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3546 restnm, egrptpALL_GENREST, bVerbose, wi);
3547 nuser = str_nelem(is->user2, MAXPTR254, ptr1);
3548 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3549 restnm, egrptpALL_GENREST, bVerbose, wi);
3550 nuser = str_nelem(is->x_compressed_groups, MAXPTR254, ptr1);
3551 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3552 restnm, egrptpONE, bVerbose, wi);
3553 nofg = str_nelem(is->orirefitgrp, MAXPTR254, ptr1);
3554 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3555 restnm, egrptpALL_GENREST, bVerbose, wi);
3556
3557 /* QMMM input processing */
3558 nQMg = str_nelem(is->QMMM, MAXPTR254, ptr1);
3559 nQMmethod = str_nelem(is->QMmethod, MAXPTR254, ptr2);
3560 nQMbasis = str_nelem(is->QMbasis, MAXPTR254, ptr3);
3561 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3562 {
3563 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3563
, "Invalid QMMM input: %d groups %d basissets"
3564 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3565 }
3566 /* group rest, if any, is always MM! */
3567 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3568 restnm, egrptpALL_GENREST, bVerbose, wi);
3569 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3570 ir->opts.ngQM = nQMg;
3571 snew(ir->opts.QMmethod, nr)(ir->opts.QMmethod) = save_calloc("ir->opts.QMmethod", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3571, (nr), sizeof(*(ir->opts.QMmethod)))
;
3572 snew(ir->opts.QMbasis, nr)(ir->opts.QMbasis) = save_calloc("ir->opts.QMbasis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3572, (nr), sizeof(*(ir->opts.QMbasis)))
;
3573 for (i = 0; i < nr; i++)
3574 {
3575 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3576 * converted to the corresponding enum in names.c
3577 */
3578 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3579 eQMmethod_names);
3580 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3581 eQMbasis_names);
3582
3583 }
3584 nQMmult = str_nelem(is->QMmult, MAXPTR254, ptr1);
Value stored to 'nQMmult' is never read
3585 nQMcharge = str_nelem(is->QMcharge, MAXPTR254, ptr2);
3586 nbSH = str_nelem(is->bSH, MAXPTR254, ptr3);
3587 snew(ir->opts.QMmult, nr)(ir->opts.QMmult) = save_calloc("ir->opts.QMmult", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3587, (nr), sizeof(*(ir->opts.QMmult)))
;
3588 snew(ir->opts.QMcharge, nr)(ir->opts.QMcharge) = save_calloc("ir->opts.QMcharge", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3588, (nr), sizeof(*(ir->opts.QMcharge)))
;
3589 snew(ir->opts.bSH, nr)(ir->opts.bSH) = save_calloc("ir->opts.bSH", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3589, (nr), sizeof(*(ir->opts.bSH)))
;
3590
3591 for (i = 0; i < nr; i++)
3592 {
3593 ir->opts.QMmult[i] = strtol(ptr1[i], NULL((void*)0), 10);
3594 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL((void*)0), 10);
3595 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3596 }
3597
3598 nCASelec = str_nelem(is->CASelectrons, MAXPTR254, ptr1);
3599 nCASorb = str_nelem(is->CASorbitals, MAXPTR254, ptr2);
3600 snew(ir->opts.CASelectrons, nr)(ir->opts.CASelectrons) = save_calloc("ir->opts.CASelectrons"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3600, (nr), sizeof(*(ir->opts.CASelectrons)))
;
3601 snew(ir->opts.CASorbitals, nr)(ir->opts.CASorbitals) = save_calloc("ir->opts.CASorbitals"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3601, (nr), sizeof(*(ir->opts.CASorbitals)))
;
3602 for (i = 0; i < nr; i++)
3603 {
3604 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL((void*)0), 10);
3605 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL((void*)0), 10);
3606 }
3607 /* special optimization options */
3608
3609 nbOPT = str_nelem(is->bOPT, MAXPTR254, ptr1);
3610 nbTS = str_nelem(is->bTS, MAXPTR254, ptr2);
3611 snew(ir->opts.bOPT, nr)(ir->opts.bOPT) = save_calloc("ir->opts.bOPT", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3611, (nr), sizeof(*(ir->opts.bOPT)))
;
3612 snew(ir->opts.bTS, nr)(ir->opts.bTS) = save_calloc("ir->opts.bTS", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3612, (nr), sizeof(*(ir->opts.bTS)))
;
3613 for (i = 0; i < nr; i++)
3614 {
3615 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3616 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3617 }
3618 nSAon = str_nelem(is->SAon, MAXPTR254, ptr1);
3619 nSAoff = str_nelem(is->SAoff, MAXPTR254, ptr2);
3620 nSAsteps = str_nelem(is->SAsteps, MAXPTR254, ptr3);
3621 snew(ir->opts.SAon, nr)(ir->opts.SAon) = save_calloc("ir->opts.SAon", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3621, (nr), sizeof(*(ir->opts.SAon)))
;
3622 snew(ir->opts.SAoff, nr)(ir->opts.SAoff) = save_calloc("ir->opts.SAoff", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3622, (nr), sizeof(*(ir->opts.SAoff)))
;
3623 snew(ir->opts.SAsteps, nr)(ir->opts.SAsteps) = save_calloc("ir->opts.SAsteps", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3623, (nr), sizeof(*(ir->opts.SAsteps)))
;
3624
3625 for (i = 0; i < nr; i++)
3626 {
3627 ir->opts.SAon[i] = strtod(ptr1[i], NULL((void*)0));
3628 ir->opts.SAoff[i] = strtod(ptr2[i], NULL((void*)0));
3629 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL((void*)0), 10);
3630 }
3631 /* end of QMMM input */
3632
3633 if (bVerbose)
3634 {
3635 for (i = 0; (i < egcNR); i++)
3636 {
3637 fprintf(stderrstderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3638 for (j = 0; (j < groups->grps[i].nr); j++)
3639 {
3640 fprintf(stderrstderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3641 }
3642 fprintf(stderrstderr, "\n");
3643 }
3644 }
3645
3646 nr = groups->grps[egcENER].nr;
3647 snew(ir->opts.egp_flags, nr*nr)(ir->opts.egp_flags) = save_calloc("ir->opts.egp_flags"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3647, (nr*nr), sizeof(*(ir->opts.egp_flags)))
;
3648
3649 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL(1<<0));
3650 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3651 {
3652 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3653 }
3654 if (bExcl && 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)
)
3655 {
3656 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3657 }
3658
3659 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE(1<<1));
3660 if (bTable && !(ir->vdwtype == evdwUSER) &&
3661 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3662 !(ir->coulombtype == eelPMEUSERSWITCH))
3663 {
3664 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3664
, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3665 }
3666
3667 decode_cos(is->efield_x, &(ir->ex[XX0]));
3668 decode_cos(is->efield_xt, &(ir->et[XX0]));
3669 decode_cos(is->efield_y, &(ir->ex[YY1]));
3670 decode_cos(is->efield_yt, &(ir->et[YY1]));
3671 decode_cos(is->efield_z, &(ir->ex[ZZ2]));
3672 decode_cos(is->efield_zt, &(ir->et[ZZ2]));
3673
3674 if (ir->bAdress)
3675 {
3676 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3677 }
3678
3679 for (i = 0; (i < grps->nr); i++)
3680 {
3681 sfree(gnames[i])save_free("gnames[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3681, (gnames[i]))
;
3682 }
3683 sfree(gnames)save_free("gnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3683, (gnames))
;
3684 done_blocka(grps);
3685 sfree(grps)save_free("grps", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3685, (grps))
;
3686
3687}
3688
3689
3690
3691static void check_disre(gmx_mtop_t *mtop)
3692{
3693 gmx_ffparams_t *ffparams;
3694 t_functype *functype;
3695 t_iparams *ip;
3696 int i, ndouble, ftype;
3697 int label, old_label;
3698
3699 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3700 {
3701 ffparams = &mtop->ffparams;
3702 functype = ffparams->functype;
3703 ip = ffparams->iparams;
3704 ndouble = 0;
3705 old_label = -1;
3706 for (i = 0; i < ffparams->ntypes; i++)
3707 {
3708 ftype = functype[i];
3709 if (ftype == F_DISRES)
3710 {
3711 label = ip[i].disres.label;
3712 if (label == old_label)
3713 {
3714 fprintf(stderrstderr, "Distance restraint index %d occurs twice\n", label);
3715 ndouble++;
3716 }
3717 old_label = label;
3718 }
3719 }
3720 if (ndouble > 0)
3721 {
3722 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3722
, "Found %d double distance restraint indices,\n"
3723 "probably the parameters for multiple pairs in one restraint "
3724 "are not identical\n", ndouble);
3725 }
3726 }
3727}
3728
3729static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3730 gmx_bool posres_only,
3731 ivec AbsRef)
3732{
3733 int d, g, i;
3734 gmx_mtop_ilistloop_t iloop;
3735 t_ilist *ilist;
3736 int nmol;
3737 t_iparams *pr;
3738
3739 clear_ivec(AbsRef);
3740
3741 if (!posres_only)
3742 {
3743 /* Check the COM */
3744 for (d = 0; d < DIM3; d++)
3745 {
3746 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3747 }
3748 /* Check for freeze groups */
3749 for (g = 0; g < ir->opts.ngfrz; g++)
3750 {
3751 for (d = 0; d < DIM3; d++)
3752 {
3753 if (ir->opts.nFreeze[g][d] != 0)
3754 {
3755 AbsRef[d] = 1;
3756 }
3757 }
3758 }
3759 }
3760
3761 /* Check for position restraints */
3762 iloop = gmx_mtop_ilistloop_init(sys);
3763 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3764 {
3765 if (nmol > 0 &&
3766 (AbsRef[XX0] == 0 || AbsRef[YY1] == 0 || AbsRef[ZZ2] == 0))
3767 {
3768 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3769 {
3770 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3771 for (d = 0; d < DIM3; d++)
3772 {
3773 if (pr->posres.fcA[d] != 0)
3774 {
3775 AbsRef[d] = 1;
3776 }
3777 }
3778 }
3779 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3780 {
3781 /* Check for flat-bottom posres */
3782 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3783 if (pr->fbposres.k != 0)
3784 {
3785 switch (pr->fbposres.geom)
3786 {
3787 case efbposresSPHERE:
3788 AbsRef[XX0] = AbsRef[YY1] = AbsRef[ZZ2] = 1;
3789 break;
3790 case efbposresCYLINDER:
3791 AbsRef[XX0] = AbsRef[YY1] = 1;
3792 break;
3793 case efbposresX: /* d=XX */
3794 case efbposresY: /* d=YY */
3795 case efbposresZ: /* d=ZZ */
3796 d = pr->fbposres.geom - efbposresX;
3797 AbsRef[d] = 1;
3798 break;
3799 default:
3800 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3800
, " Invalid geometry for flat-bottom position restraint.\n"
3801 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3802 pr->fbposres.geom);
3803 }
3804 }
3805 }
3806 }
3807 }
3808
3809 return (AbsRef[XX0] != 0 && AbsRef[YY1] != 0 && AbsRef[ZZ2] != 0);
3810}
3811
3812static void
3813check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3814 gmx_bool *bC6ParametersWorkWithGeometricRules,
3815 gmx_bool *bC6ParametersWorkWithLBRules,
3816 gmx_bool *bLBRulesPossible)
3817{
3818 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3819 int *typecount;
3820 real tol;
3821 double geometricdiff, LBdiff;
3822 double c6i, c6j, c12i, c12j;
3823 double c6, c6_geometric, c6_LB;
3824 double sigmai, sigmaj, epsi, epsj;
3825 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3826 const char *ptr;
3827
3828 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3829 * force-field floating point parameters.
3830 */
3831 tol = 1e-5;
3832 ptr = getenv("GMX_LJCOMB_TOL");
3833 if (ptr != NULL((void*)0))
3834 {
3835 double dbl;
3836
3837 sscanf(ptr, "%lf", &dbl);
3838 tol = dbl;
3839 }
3840
3841 *bC6ParametersWorkWithLBRules = TRUE1;
3842 *bC6ParametersWorkWithGeometricRules = TRUE1;
3843 bCanDoLBRules = TRUE1;
3844 bCanDoGeometricRules = TRUE1;
3845 ntypes = mtop->ffparams.atnr;
3846 snew(typecount, ntypes)(typecount) = save_calloc("typecount", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3846, (ntypes), sizeof(*(typecount)))
;
3847 gmx_mtop_count_atomtypes(mtop, state, typecount);
3848 geometricdiff = LBdiff = 0.0;
3849 *bLBRulesPossible = TRUE1;
3850 for (tpi = 0; tpi < ntypes; ++tpi)
3851 {
3852 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3853 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3854 for (tpj = tpi; tpj < ntypes; ++tpj)
3855 {
3856 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3857 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3858 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3859 c6_geometric = sqrt(c6i * c6j);
3860 if (!gmx_numzero(c6_geometric))
3861 {
3862 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3863 {
3864 sigmai = pow(c12i / c6i, 1.0/6.0);
3865 sigmaj = pow(c12j / c6j, 1.0/6.0);
3866 epsi = c6i * c6i /(4.0 * c12i);
3867 epsj = c6j * c6j /(4.0 * c12j);
3868 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3869 }
3870 else
3871 {
3872 *bLBRulesPossible = FALSE0;
3873 c6_LB = c6_geometric;
3874 }
3875 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3876 }
3877
3878 if (FALSE0 == bCanDoLBRules)
3879 {
3880 *bC6ParametersWorkWithLBRules = FALSE0;
3881 }
3882
3883 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3884
3885 if (FALSE0 == bCanDoGeometricRules)
3886 {
3887 *bC6ParametersWorkWithGeometricRules = FALSE0;
3888 }
3889 }
3890 }
3891 sfree(typecount)save_free("typecount", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 3891, (typecount))
;
3892}
3893
3894static void
3895check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3896 warninp_t wi)
3897{
3898 char err_buf[256];
3899 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3900
3901 check_combination_rule_differences(mtop, 0,
3902 &bC6ParametersWorkWithGeometricRules,
3903 &bC6ParametersWorkWithLBRules,
3904 &bLBRulesPossible);
3905 if (ir->ljpme_combination_rule == eljpmeLB)
3906 {
3907 if (FALSE0 == bC6ParametersWorkWithLBRules || FALSE0 == bLBRulesPossible)
3908 {
3909 warning(wi, "You are using arithmetic-geometric combination rules "
3910 "in LJ-PME, but your non-bonded C6 parameters do not "
3911 "follow these rules.");
3912 }
3913 }
3914 else
3915 {
3916 if (FALSE0 == bC6ParametersWorkWithGeometricRules)
3917 {
3918 if (ir->eDispCorr != edispcNO)
3919 {
3920 warning_note(wi, "You are using geometric combination rules in "
3921 "LJ-PME, but your non-bonded C6 parameters do "
3922 "not follow these rules. "
3923 "This will introduce very small errors in the forces and energies in "
3924 "your simulations. Dispersion correction will correct total energy "
3925 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3926 }
3927 else
3928 {
3929 warning_note(wi, "You are using geometric combination rules in "
3930 "LJ-PME, but your non-bonded C6 parameters do "
3931 "not follow these rules. "
3932 "This will introduce very small errors in the forces and energies in "
3933 "your simulations. If your system is homogeneous, consider using dispersion correction "
3934 "for the total energy and pressure.");
3935 }
3936 }
3937 }
3938}
3939
3940void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3941 warninp_t wi)
3942{
3943 char err_buf[256];
3944 int i, m, c, nmol, npct;
3945 gmx_bool bCharge, bAcc;
3946 real gdt_max, *mgrp, mt;
3947 rvec acc;
3948 gmx_mtop_atomloop_block_t aloopb;
3949 gmx_mtop_atomloop_all_t aloop;
3950 t_atom *atom;
3951 ivec AbsRef;
3952 char warn_buf[STRLEN4096];
3953
3954 set_warning_line(wi, mdparin, -1);
3955
3956 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD)
&& !EI_SD(ir->eI)((ir->eI) == eiSD1 || (ir->eI) == eiSD2) && ir->eI != eiBD &&
3957 ir->comm_mode == ecmNO &&
3958 !(absolute_reference(ir, sys, FALSE0, AbsRef) || ir->nsteps <= 10) &&
3959 !ETC_ANDERSEN(ir->etc)(((ir->etc) == etcANDERSENMASSIVE) || ((ir->etc) == etcANDERSEN
))
)
3960 {
3961 warning(wi, "You are not using center of mass motion removal (mdp option comm-mode), numerical rounding errors can lead to build up of kinetic energy of the center of mass");
3962 }
3963
3964 /* Check for pressure coupling with absolute position restraints */
3965 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3966 {
3967 absolute_reference(ir, sys, TRUE1, AbsRef);
3968 {
3969 for (m = 0; m < DIM3; m++)
3970 {
3971 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3972 {
3973 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3974 break;
3975 }
3976 }
3977 }
3978 }
3979
3980 bCharge = FALSE0;
3981 aloopb = gmx_mtop_atomloop_block_init(sys);
3982 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3983 {
3984 if (atom->q != 0 || atom->qB != 0)
3985 {
3986 bCharge = TRUE1;
3987 }
3988 }
3989
3990 if (!bCharge)
3991 {
3992 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)
)
3993 {
3994 sprintf(err_buf,
3995 "You are using full electrostatics treatment %s for a system without charges.\n"
3996 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3997 EELTYPE(ir->coulombtype)((((ir->coulombtype) < 0) || ((ir->coulombtype) >=
(eelNR))) ? "UNDEFINED" : (eel_names)[ir->coulombtype])
, EELTYPE(eelCUT)((((eelCUT) < 0) || ((eelCUT) >= (eelNR))) ? "UNDEFINED"
: (eel_names)[eelCUT])
);
3998 warning(wi, err_buf);
3999 }
4000 }
4001 else
4002 {
4003 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4004 {
4005 sprintf(err_buf,
4006 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4007 "You might want to consider using %s electrostatics.\n",
4008 EELTYPE(eelPME)((((eelPME) < 0) || ((eelPME) >= (eelNR))) ? "UNDEFINED"
: (eel_names)[eelPME])
);
4009 warning_note(wi, err_buf);
4010 }
4011 }
4012
4013 /* Check if combination rules used in LJ-PME are the same as in the force field */
4014 if (EVDW_PME(ir->vdwtype)((ir->vdwtype) == evdwPME))
4015 {
4016 check_combination_rules(ir, sys, wi);
4017 }
4018
4019 /* Generalized reaction field */
4020 if (ir->opts.ngtc == 0)
4021 {
4022 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4023 eel_names[eelGRF]);
4024 CHECK(ir->coulombtype == eelGRF)_low_check(ir->coulombtype == eelGRF, err_buf, wi);
4025 }
4026 else
4027 {
4028 sprintf(err_buf, "When using coulombtype = %s"
4029 " ref-t for temperature coupling should be > 0",
4030 eel_names[eelGRF]);
4031 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0))_low_check((ir->coulombtype == eelGRF) && (ir->
opts.ref_t[0] <= 0), err_buf, wi)
;
4032 }
4033
4034 if (ir->eI == eiSD1 &&
4035 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
4036 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
4037 {
4038 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
4039 warning_note(wi, warn_buf);
4040 }
4041
4042 bAcc = FALSE0;
4043 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4044 {
4045 for (m = 0; (m < DIM3); m++)
4046 {
4047 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4048 {
4049 bAcc = TRUE1;
4050 }
4051 }
4052 }
4053 if (bAcc)
4054 {
4055 clear_rvec(acc);
4056 snew(mgrp, sys->groups.grps[egcACC].nr)(mgrp) = save_calloc("mgrp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 4056, (sys->groups.grps[egcACC].nr), sizeof(*(mgrp)))
;
4057 aloop = gmx_mtop_atomloop_all_init(sys);
4058 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4059 {
4060 mgrp[ggrpnr(&sys->groups, egcACC, i)((&sys->groups)->grpnr[egcACC] ? (&sys->groups
)->grpnr[egcACC][i] : 0)
] += atom->m;
4061 }
4062 mt = 0.0;
4063 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4064 {
4065 for (m = 0; (m < DIM3); m++)
4066 {
4067 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4068 }
4069 mt += mgrp[i];
4070 }
4071 for (m = 0; (m < DIM3); m++)
4072 {
4073 if (fabs(acc[m]) > 1e-6)
4074 {
4075 const char *dim[DIM3] = { "X", "Y", "Z" };
4076 fprintf(stderrstderr,
4077 "Net Acceleration in %s direction, will %s be corrected\n",
4078 dim[m], ir->nstcomm != 0 ? "" : "not");
4079 if (ir->nstcomm != 0 && m < ndof_com(ir))
4080 {
4081 acc[m] /= mt;
4082 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4083 {
4084 ir->opts.acc[i][m] -= acc[m];
4085 }
4086 }
4087 }
4088 }
4089 sfree(mgrp)save_free("mgrp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 4089, (mgrp))
;
4090 }
4091
4092 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4093 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS1.11022302E-16))
4094 {
4095 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 4095
, "Soft-core interactions are only supported with VdW repulsion power 12");
4096 }
4097
4098 if (ir->ePull != epullNO)
4099 {
4100 gmx_bool bPullAbsoluteRef;
4101
4102 bPullAbsoluteRef = FALSE0;
4103 for (i = 0; i < ir->pull->ncoord; i++)
4104 {
4105 bPullAbsoluteRef = bPullAbsoluteRef ||
4106 ir->pull->coord[i].group[0] == 0 ||
4107 ir->pull->coord[i].group[1] == 0;
4108 }
4109 if (bPullAbsoluteRef)
4110 {
4111 absolute_reference(ir, sys, FALSE0, AbsRef);
4112 for (m = 0; m < DIM3; m++)
4113 {
4114 if (ir->pull->dim[m] && !AbsRef[m])
4115 {
4116 warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
4117 break;
4118 }
4119 }
4120 }
4121
4122 if (ir->pull->eGeom == epullgDIRPBC)
4123 {
4124 for (i = 0; i < 3; i++)
4125 {
4126 for (m = 0; m <= i; m++)
4127 {
4128 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4129 ir->deform[i][m] != 0)
4130 {
4131 for (c = 0; c < ir->pull->ncoord; c++)
4132 {
4133 if (ir->pull->coord[c].vec[m] != 0)
4134 {
4135 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/readir.c"
, 4135
, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom)((((ir->pull->eGeom) < 0) || ((ir->pull->eGeom
) >= (epullgNR))) ? "UNDEFINED" : (epullg_names)[ir->pull
->eGeom])
, 'x'+m);
4136 }
4137 }
4138 }
4139 }
4140 }
4141 }
4142 }
4143
4144 check_disre(sys);
4145}
4146
4147void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4148{
4149 real min_size;
4150 gmx_bool bTWIN;
4151 char warn_buf[STRLEN4096];
4152 const char *ptr;
4153
4154 ptr = check_box(ir->ePBC, box);
4155 if (ptr)
4156 {
4157 warning_error(wi, ptr);
4158 }
4159
4160 if (bConstr && ir->eConstrAlg == econtSHAKE)
4161 {
4162 if (ir->shake_tol <= 0.0)
4163 {
4164 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4165 ir->shake_tol);
4166 warning_error(wi, warn_buf);
4167 }
4168
4169 if (IR_TWINRANGE(*ir)((*ir).rlist > 0 && ((*ir).rlistlong == 0 || (*ir)
.rlistlong > (*ir).rlist))
&& ir->nstlist > 1)
4170 {
4171 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4172 if (ir->epc == epcNO)
4173 {
4174 warning(wi, warn_buf);
4175 }
4176 else
4177 {
4178 warning_error(wi, warn_buf);
4179 }
4180 }
4181 }
4182
4183 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4184 {
4185 /* If we have Lincs constraints: */
4186 if (ir->eI == eiMD && ir->etc == etcNO &&
4187 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4188 {
4189 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4190 warning_note(wi, warn_buf);
4191 }
4192
4193 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4194 {
4195 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4196 warning_note(wi, warn_buf);
4197 }
4198 if (ir->epc == epcMTTK)
4199 {
4200 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4201 }
4202 }
4203
4204 if (bConstr && ir->epc == epcMTTK)
4205 {
4206 warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
4207 }
4208
4209 if (ir->LincsWarnAngle > 90.0)
4210 {
4211 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4212 warning(wi, warn_buf);
4213 ir->LincsWarnAngle = 90.0;
4214 }
4215
4216 if (ir->ePBC != epbcNONE)
4217 {
4218 if (ir->nstlist == 0)
4219 {
4220 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4221 }
4222 bTWIN = (ir->rlistlong > ir->rlist);
4223 if (ir->ns_type == ensGRID)
4224 {
4225 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4226 {
4227 sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease %s.\n",
4228 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4229 warning_error(wi, warn_buf);
4230 }
4231 }
4232 else
4233 {
4234 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]))(((box[0][0]) < ((((box[1][1]) < (box[2][2])) ? (box[1]
[1]) : (box[2][2]) ))) ? (box[0][0]) : ((((box[1][1]) < (box
[2][2])) ? (box[1][1]) : (box[2][2]) )) )
;
4235 if (2*ir->rlistlong >= min_size)
4236 {
4237 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4238 warning_error(wi, warn_buf);
4239 if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0))
4240 {
4241 fprintf(stderrstderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4242 }
4243 }
4244 }
4245 }
4246}
4247
4248void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4249 rvec *x,
4250 warninp_t wi)
4251{
4252 real rvdw1, rvdw2, rcoul1, rcoul2;
4253 char warn_buf[STRLEN4096];
4254
4255 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4256
4257 if (rvdw1 > 0)
4258 {
4259 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4260 rvdw1, rvdw2);
4261 }
4262 if (rcoul1 > 0)
4263 {
4264 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4265 rcoul1, rcoul2);
4266 }
4267
4268 if (ir->rlist > 0)
4269 {
4270 if (rvdw1 + rvdw2 > ir->rlist ||
4271 rcoul1 + rcoul2 > ir->rlist)
4272 {
4273 sprintf(warn_buf,
4274 "The sum of the two largest charge group radii (%f) "
4275 "is larger than rlist (%f)\n",
4276 max(rvdw1+rvdw2, rcoul1+rcoul2)(((rvdw1+rvdw2) > (rcoul1+rcoul2)) ? (rvdw1+rvdw2) : (rcoul1
+rcoul2) )
, ir->rlist);
4277 warning(wi, warn_buf);
4278 }
4279 else
4280 {
4281 /* Here we do not use the zero at cut-off macro,
4282 * since user defined interactions might purposely
4283 * not be zero at the cut-off.
4284 */
4285 if (ir_vdw_is_zero_at_cutoff(ir) &&
4286 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4287 {
4288 sprintf(warn_buf, "The sum of the two largest charge group "
4289 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4290 "With exact cut-offs, better performance can be "
4291 "obtained with cutoff-scheme = %s, because it "
4292 "does not use charge groups at all.",
4293 rvdw1+rvdw2,
4294 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4295 ir->rlistlong, ir->rvdw,
4296 ecutscheme_names[ecutsVERLET]);
4297 if (ir_NVE(ir))
4298 {
4299 warning(wi, warn_buf);
4300 }
4301 else
4302 {
4303 warning_note(wi, warn_buf);
4304 }
4305 }
4306 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4307 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4308 {
4309 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4310 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4311 rcoul1+rcoul2,
4312 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4313 ir->rlistlong, ir->rcoulomb,
4314 ecutscheme_names[ecutsVERLET]);
4315 if (ir_NVE(ir))
4316 {
4317 warning(wi, warn_buf);
4318 }
4319 else
4320 {
4321 warning_note(wi, warn_buf);
4322 }
4323 }
4324 }
4325 }
4326}