3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
47 #include "gmx_fatal.h"
59 #include "mtop_util.h"
64 /* Resource parameters
65 * Do not change any of these until you read the instruction
66 * in readinp.h. Some cpp's do not take spaces after the backslash
67 * (like the c-shell), which will give you a very weird compiler
71 static char tcgrps
[STRLEN
],tau_t
[STRLEN
],ref_t
[STRLEN
],
72 acc
[STRLEN
],accgrps
[STRLEN
],freeze
[STRLEN
],frdim
[STRLEN
],
73 energy
[STRLEN
],user1
[STRLEN
],user2
[STRLEN
],vcm
[STRLEN
],xtc_grps
[STRLEN
],
74 couple_moltype
[STRLEN
],orirefitgrp
[STRLEN
],egptable
[STRLEN
],egpexcl
[STRLEN
],
75 wall_atomtype
[STRLEN
],wall_density
[STRLEN
],deform
[STRLEN
],QMMM
[STRLEN
];
76 static char foreign_lambda
[STRLEN
];
77 static char **pull_grp
;
78 static char **rot_grp
;
79 static char anneal
[STRLEN
],anneal_npoints
[STRLEN
],
80 anneal_time
[STRLEN
],anneal_temp
[STRLEN
];
81 static char QMmethod
[STRLEN
],QMbasis
[STRLEN
],QMcharge
[STRLEN
],QMmult
[STRLEN
],
82 bSH
[STRLEN
],CASorbitals
[STRLEN
], CASelectrons
[STRLEN
],SAon
[STRLEN
],
83 SAoff
[STRLEN
],SAsteps
[STRLEN
],bTS
[STRLEN
],bOPT
[STRLEN
];
84 static char efield_x
[STRLEN
],efield_xt
[STRLEN
],efield_y
[STRLEN
],
85 efield_yt
[STRLEN
],efield_z
[STRLEN
],efield_zt
[STRLEN
];
87 enum { egrptpALL
, egrptpALL_GENREST
, egrptpPART
, egrptpONE
};
90 void init_ir(t_inputrec
*ir
, t_gromppopts
*opts
)
92 snew(opts
->include
,STRLEN
);
93 snew(opts
->define
,STRLEN
);
96 static void _low_check(bool b
,char *s
,int *n
)
99 fprintf(stderr
,"\nERROR: %s\n\n",s
);
104 static void check_nst(const char *desc_nst
,int nst
,
105 const char *desc_p
,int *p
)
109 if (*p
> 0 && *p
% nst
!= 0)
111 /* Round up to the next multiple of nst */
112 *p
= ((*p
)/nst
+ 1)*nst
;
113 sprintf(buf
,"%s should be a multiple of %s, changing %s to %d\n",
114 desc_p
,desc_nst
,desc_p
,*p
);
119 void check_ir(const char *mdparin
,t_inputrec
*ir
, t_gromppopts
*opts
,
121 /* Check internal consistency */
123 /* Strange macro: first one fills the err_buf, and then one can check
124 * the condition, which will print the message and increase the error
127 #define CHECK(b) _low_check(b,err_buf,nerror)
132 set_warning_line(mdparin
,-1);
134 /* BASIC CUT-OFF STUFF */
135 if (ir
->rlist
== 0 ||
136 !(( EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) && ir
->rcoulomb
> ir
->rlist
) ||
137 (EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
) && ir
->rvdw
> ir
->rlist
))) {
138 /* No switched potential and/or no twin-range:
139 * we can set the long-range cut-off to the maximum of the other cut-offs.
141 ir
->rlistlong
= max_cutoff(ir
->rlist
,max_cutoff(ir
->rvdw
,ir
->rcoulomb
));
142 } else if (ir
->rlistlong
< 0) {
143 ir
->rlistlong
= max_cutoff(ir
->rlist
,max_cutoff(ir
->rvdw
,ir
->rcoulomb
));
144 sprintf(warn_buf
,"rlistlong was not set, setting it to %g (no buffer)",
148 if (ir
->rlistlong
== 0 && ir
->ePBC
!= epbcNONE
) {
149 warning_error("Can not have an infinite cut-off with PBC");
151 if (ir
->rlistlong
> 0 && (ir
->rlist
== 0 || ir
->rlistlong
< ir
->rlist
)) {
152 warning_error("rlistlong can not be shorter than rlist");
154 if (IR_TWINRANGE(*ir
) && ir
->nstlist
<= 0) {
155 warning_error("Can not have nstlist<=0 with twin-range interactions");
158 /* GENERAL INTEGRATOR STUFF */
159 if (!(ir
->eI
== eiMD
|| EI_VV(ir
->eI
))) {
162 if (!EI_DYNAMICS(ir
->eI
)) {
165 if (EI_DYNAMICS(ir
->eI
)) {
166 if (ir
->nstcalcenergy
< 0) {
167 gmx_fatal(FARGS
,"Can not have nstcalcenergy < 0");
169 if ((ir
->etc
!= etcNO
|| ir
->epc
!= epcNO
) && ir
->nstcalcenergy
== 0) {
170 gmx_fatal(FARGS
,"Can not have nstcalcenergy=0 with global T/P-coupling");
172 if (IR_TWINRANGE(*ir
)) {
173 check_nst("nstlist",ir
->nstlist
,"nstcalcenergy",&ir
->nstcalcenergy
);
175 dt_coupl
= ir
->nstcalcenergy
*ir
->delta_t
;
177 if (ir
->nstcalcenergy
> 1) {
178 /* Energy and log file writing trigger energy calculation,
179 * so we need some checks.
181 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstenergy",&ir
->nstenergy
);
182 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstlog",&ir
->nstlog
);
183 if (ir
->efep
!= efepNO
) {
184 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstdhdl",&ir
->nstdhdl
);
190 if ((EI_SD(ir
->eI
) || ir
->eI
== eiBD
) &&
191 ir
->bContinuation
&& ir
->ld_seed
!= -1) {
192 warning_note("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)");
196 if (EI_TPI(ir
->eI
)) {
197 sprintf(err_buf
,"TPI only works with pbc = %s",epbc_names
[epbcXYZ
]);
198 CHECK(ir
->ePBC
!= epbcXYZ
);
199 sprintf(err_buf
,"TPI only works with ns = %s",ens_names
[ensGRID
]);
200 CHECK(ir
->ns_type
!= ensGRID
);
201 sprintf(err_buf
,"with TPI nstlist should be larger than zero");
202 CHECK(ir
->nstlist
<= 0);
203 sprintf(err_buf
,"TPI does not work with full electrostatics other than PME");
204 CHECK(EEL_FULL(ir
->coulombtype
) && !EEL_PME(ir
->coulombtype
));
208 if ( (opts
->nshake
> 0) && (opts
->bMorse
) ) {
210 "Using morse bond-potentials while constraining bonds is useless");
214 sprintf(err_buf
,"shake_tol must be > 0 instead of %g while using shake",
216 CHECK(((ir
->shake_tol
<= 0.0) && (opts
->nshake
>0) &&
217 (ir
->eConstrAlg
== econtSHAKE
)));
220 sprintf(err_buf
,"walls only work with pbc=%s",epbc_names
[epbcXY
]);
221 CHECK(ir
->nwall
&& ir
->ePBC
!=epbcXY
);
224 if (ir
->ePBC
!= epbcXYZ
&& ir
->nwall
!= 2) {
225 if (ir
->ePBC
== epbcNONE
) {
226 if (ir
->epc
!= epcNO
) {
227 warning("Turning off pressure coupling for vacuum system");
231 sprintf(err_buf
,"Can not have pressure coupling with pbc=%s",
232 epbc_names
[ir
->ePBC
]);
233 CHECK(ir
->epc
!= epcNO
);
235 sprintf(err_buf
,"Can not have Ewald with pbc=%s",epbc_names
[ir
->ePBC
]);
236 CHECK(EEL_FULL(ir
->coulombtype
));
238 sprintf(err_buf
,"Can not have dispersion correction with pbc=%s",
239 epbc_names
[ir
->ePBC
]);
240 CHECK(ir
->eDispCorr
!= edispcNO
);
243 if (ir
->rlist
== 0.0) {
244 sprintf(err_buf
,"can only have neighborlist cut-off zero (=infinite)\n"
245 "with coulombtype = %s or coulombtype = %s\n"
246 "without periodic boundary conditions (pbc = %s) and\n"
247 "rcoulomb and rvdw set to zero",
248 eel_names
[eelCUT
],eel_names
[eelUSER
],epbc_names
[epbcNONE
]);
249 CHECK(((ir
->coulombtype
!= eelCUT
) && (ir
->coulombtype
!= eelUSER
)) ||
250 (ir
->ePBC
!= epbcNONE
) ||
251 (ir
->rcoulomb
!= 0.0) || (ir
->rvdw
!= 0.0));
253 if (ir
->nstlist
< 0) {
254 warning_error("Can not have heuristic neighborlist updates without cut-off");
256 if (ir
->nstlist
> 0) {
257 warning_note("Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
262 if (ir
->nstcomm
== 0) {
263 ir
->comm_mode
= ecmNO
;
265 if (ir
->comm_mode
!= ecmNO
) {
266 if (ir
->nstcomm
< 0) {
267 warning("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");
268 ir
->nstcomm
= abs(ir
->nstcomm
);
271 if (ir
->nstcalcenergy
> 0 && ir
->nstcomm
< ir
->nstcalcenergy
) {
272 warning_note("nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
273 ir
->nstcomm
= ir
->nstcalcenergy
;
276 if (ir
->comm_mode
== ecmANGULAR
) {
277 sprintf(err_buf
,"Can not remove the rotation around the center of mass with periodic molecules");
278 CHECK(ir
->bPeriodicMols
);
279 if (ir
->ePBC
!= epbcNONE
)
280 warning("Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
284 if (EI_STATE_VELOCITY(ir
->eI
) && ir
->ePBC
== epbcNONE
&& ir
->comm_mode
!= ecmANGULAR
) {
285 warning_note("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.");
288 sprintf(err_buf
,"Free-energy not implemented for Ewald and PPPM");
289 CHECK((ir
->coulombtype
==eelEWALD
|| ir
->coulombtype
==eelPPPM
)
290 && (ir
->efep
!=efepNO
));
292 sprintf(err_buf
,"Twin-range neighbour searching (NS) with simple NS"
293 " algorithm not implemented");
294 CHECK(((ir
->rcoulomb
> ir
->rlist
) || (ir
->rvdw
> ir
->rlist
))
295 && (ir
->ns_type
== ensSIMPLE
));
297 /* TEMPERATURE COUPLING */
298 if(ir
->etc
== etcYES
) {
299 ir
->etc
= etcBERENDSEN
;
300 warning_note("Old option for temperature coupling given: "
301 "changing \"yes\" to \"Berendsen\"\n");
303 if (ir
->etc
== etcNOSEHOOVER
) {
304 if (ir
->opts
.nhchainlength
< 1)
306 sprintf(warn_buf
,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir
->opts
.nhchainlength
);
307 ir
->opts
.nhchainlength
=1;
311 if (ir
->etc
==etcNOSEHOOVER
&& !EI_VV(ir
->eI
) && ir
->opts
.nhchainlength
> 1) {
312 warning_note("leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
313 ir
->opts
.nhchainlength
= 1;
316 ir
->opts
.nhchainlength
= 0;
319 if (ir
->etc
== etcBERENDSEN
) {
320 sprintf(warn_buf
,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
321 ETCOUPLTYPE(ir
->etc
),ETCOUPLTYPE(etcVRESCALE
));
325 if((ir
->etc
==etcNOSEHOOVER
|| ir
->etc
==etcANDERSEN
|| ir
->etc
==etcANDERSENINTERVAL
)
326 && ir
->epc
==epcBERENDSEN
) {
327 sprintf(warn_buf
,"Using Berendsen pressure coupling invalidates the "
328 "true ensemble for the thermostat");
332 /* PRESSURE COUPLING */
333 if (ir
->epc
== epcISOTROPIC
) {
334 ir
->epc
= epcBERENDSEN
;
335 warning_note("Old option for pressure coupling given: "
336 "changing \"Isotropic\" to \"Berendsen\"\n");
339 if (ir
->epc
!= epcNO
) {
340 sprintf(err_buf
,"tau_p must be > 0 instead of %g\n",ir
->tau_p
);
341 CHECK(ir
->tau_p
<= 0);
343 if (ir
->tau_p
< 100*dt_coupl
) {
344 sprintf(warn_buf
,"For proper barostat integration tau_p (%g) should be more than two orders of magnitude larger than nstcalcenergy*dt (%g)",
349 sprintf(err_buf
,"compressibility must be > 0 when using pressure"
350 " coupling %s\n",EPCOUPLTYPE(ir
->epc
));
351 CHECK(ir
->compress
[XX
][XX
] < 0 || ir
->compress
[YY
][YY
] < 0 ||
352 ir
->compress
[ZZ
][ZZ
] < 0 ||
353 (trace(ir
->compress
) == 0 && ir
->compress
[YY
][XX
] <= 0 &&
354 ir
->compress
[ZZ
][XX
] <= 0 && ir
->compress
[ZZ
][YY
] <= 0));
356 sprintf(err_buf
,"pressure coupling with PPPM not implemented, use PME");
357 CHECK(ir
->coulombtype
== eelPPPM
);
359 } else if (ir
->coulombtype
== eelPPPM
) {
360 sprintf(warn_buf
,"The pressure with PPPM is incorrect, if you need the pressure use PME");
365 if (ir
->epc
> epcNO
) {
366 if (ir
->epc
!=epcMTTK
) {
367 warning_error("NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
373 /* More checks are in triple check (grompp.c) */
374 if (ir
->coulombtype
== eelSWITCH
) {
375 sprintf(warn_buf
,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
376 eel_names
[ir
->coulombtype
],
377 eel_names
[eelRF_ZERO
]);
381 if (EEL_RF(ir
->coulombtype
) && ir
->epsilon_rf
==1 && ir
->epsilon_r
!=1) {
382 sprintf(warn_buf
,"epsilon_r = %g and epsilon_rf = 1 with reaction field, assuming old format and exchanging epsilon_r and epsilon_rf",ir
->epsilon_r
);
384 ir
->epsilon_rf
= ir
->epsilon_r
;
388 if (getenv("GALACTIC_DYNAMICS") == NULL
) {
389 sprintf(err_buf
,"epsilon_r must be >= 0 instead of %g\n",ir
->epsilon_r
);
390 CHECK(ir
->epsilon_r
< 0);
393 if (EEL_RF(ir
->coulombtype
)) {
394 /* reaction field (at the cut-off) */
396 if (ir
->coulombtype
== eelRF_ZERO
) {
397 sprintf(err_buf
,"With coulombtype = %s, epsilon_rf must be 0",
398 eel_names
[ir
->coulombtype
]);
399 CHECK(ir
->epsilon_rf
!= 0);
402 sprintf(err_buf
,"epsilon_rf must be >= epsilon_r");
403 CHECK((ir
->epsilon_rf
< ir
->epsilon_r
&& ir
->epsilon_rf
!= 0) ||
404 (ir
->epsilon_r
== 0));
405 if (ir
->epsilon_rf
== ir
->epsilon_r
) {
406 sprintf(warn_buf
,"Using epsilon_rf = epsilon_r with %s does not make sense",
407 eel_names
[ir
->coulombtype
]);
411 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
412 * means the interaction is zero outside rcoulomb, but it helps to
413 * provide accurate energy conservation.
415 if (EEL_ZERO_AT_CUTOFF(ir
->coulombtype
)) {
416 if (EEL_SWITCHED(ir
->coulombtype
)) {
418 "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
419 eel_names
[ir
->coulombtype
]);
420 CHECK(ir
->rcoulomb_switch
>= ir
->rcoulomb
);
422 } else if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
)) {
423 sprintf(err_buf
,"With coulombtype = %s, rcoulomb must be >= rlist",
424 eel_names
[ir
->coulombtype
]);
425 CHECK(ir
->rlist
> ir
->rcoulomb
);
428 if (EEL_FULL(ir
->coulombtype
)) {
429 if (ir
->coulombtype
==eelPMESWITCH
|| ir
->coulombtype
==eelPMEUSER
) {
430 sprintf(err_buf
,"With coulombtype = %s, rcoulomb must be <= rlist",
431 eel_names
[ir
->coulombtype
]);
432 CHECK(ir
->rcoulomb
> ir
->rlist
);
434 if (ir
->coulombtype
== eelPME
) {
436 "With coulombtype = %s, rcoulomb must be equal to rlist\n"
437 "If you want optimal energy conservation or exact integration use %s",
438 eel_names
[ir
->coulombtype
],eel_names
[eelPMESWITCH
]);
441 "With coulombtype = %s, rcoulomb must be equal to rlist",
442 eel_names
[ir
->coulombtype
]);
444 CHECK(ir
->rcoulomb
!= ir
->rlist
);
448 if (EEL_PME(ir
->coulombtype
)) {
449 if (ir
->pme_order
< 3) {
450 warning_error("pme_order can not be smaller than 3");
454 if (ir
->nwall
==2 && EEL_FULL(ir
->coulombtype
)) {
455 if (ir
->ewald_geometry
== eewg3D
) {
456 sprintf(warn_buf
,"With pbc=%s you should use ewald_geometry=%s",
457 epbc_names
[ir
->ePBC
],eewg_names
[eewg3DC
]);
460 /* This check avoids extra pbc coding for exclusion corrections */
461 sprintf(err_buf
,"wall_ewald_zfac should be >= 2");
462 CHECK(ir
->wall_ewald_zfac
< 2);
465 if (EVDW_SWITCHED(ir
->vdwtype
)) {
466 sprintf(err_buf
,"With vdwtype = %s rvdw_switch must be < rvdw",
467 evdw_names
[ir
->vdwtype
]);
468 CHECK(ir
->rvdw_switch
>= ir
->rvdw
);
469 } else if (ir
->vdwtype
== evdwCUT
) {
470 sprintf(err_buf
,"With vdwtype = %s, rvdw must be >= rlist",evdw_names
[ir
->vdwtype
]);
471 CHECK(ir
->rlist
> ir
->rvdw
);
473 if ((EEL_SWITCHED(ir
->coulombtype
) || ir
->coulombtype
== eelRF_ZERO
)
474 && (ir
->rlistlong
<= ir
->rcoulomb
)) {
475 sprintf(warn_buf
,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
476 IR_TWINRANGE(*ir
) ? "rlistlong" : "rlist");
479 if (EVDW_SWITCHED(ir
->vdwtype
) && (ir
->rlistlong
<= ir
->rvdw
)) {
480 sprintf(warn_buf
,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
481 IR_TWINRANGE(*ir
) ? "rlistlong" : "rlist");
485 if (ir
->vdwtype
== evdwUSER
&& ir
->eDispCorr
!= edispcNO
) {
486 warning_note("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.");
489 if (ir
->nstlist
== -1) {
491 "nstlist=-1 only works with switched or shifted potentials,\n"
492 "suggestion: use vdw-type=%s and coulomb-type=%s",
493 evdw_names
[evdwSHIFT
],eel_names
[eelPMESWITCH
]);
494 CHECK(!(EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) &&
495 EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
)));
497 sprintf(err_buf
,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
498 CHECK(ir
->rvdw
>= ir
->rlist
|| ir
->rcoulomb
>= ir
->rlist
);
500 sprintf(err_buf
,"nstlist can not be smaller than -1");
501 CHECK(ir
->nstlist
< -1);
503 if (ir
->eI
== eiLBFGS
&& (ir
->coulombtype
==eelCUT
|| ir
->vdwtype
==evdwCUT
)
505 warning("For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
508 if (ir
->eI
== eiLBFGS
&& ir
->nbfgscorr
<= 0) {
509 warning("Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
513 if (ir
->efep
!= efepNO
) {
514 sprintf(err_buf
,"The soft-core power is %d and can only be 1 or 2",
516 CHECK(ir
->sc_alpha
!=0 && ir
->sc_power
!=1 && ir
->sc_power
!=2);
519 /* ENERGY CONSERVATION */
520 if (ir
->eI
== eiMD
&& ir
->etc
== etcNO
) {
521 if (!EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
) && ir
->rvdw
> 0) {
522 sprintf(warn_buf
,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
523 evdw_names
[evdwSHIFT
]);
526 if (!EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) && ir
->rcoulomb
> 0) {
527 sprintf(warn_buf
,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
528 eel_names
[eelPMESWITCH
],eel_names
[eelRF_ZERO
]);
533 if(ir
->coulombtype
==eelGB_NOTUSED
)
535 ir
->coulombtype
=eelCUT
;
536 ir
->implicit_solvent
=eisGBSA
;
537 fprintf(stderr
,"Note: Old option for generalized born electrostatics given:\n"
538 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
539 "setting implicit_solvent value to \"GBSA\" in input section.\n");
542 if(ir
->implicit_solvent
==eisGBSA
)
544 sprintf(err_buf
,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
545 CHECK(ir
->rgbradii
!= ir
->rlist
);
549 static int str_nelem(const char *str
,int maxptr
,char *ptr
[])
557 while (*copy
!= '\0') {
559 gmx_fatal(FARGS
,"Too many groups on line: '%s' (max is %d)",
564 while ((*copy
!= '\0') && !isspace(*copy
))
578 static void parse_n_double(char *str
,int *n
,double **r
)
583 *n
= str_nelem(str
,MAXPTR
,ptr
);
586 for(i
=0; i
<*n
; i
++) {
587 (*r
)[i
] = strtod(ptr
[i
],NULL
);
591 static void do_wall_params(t_inputrec
*ir
,
592 char *wall_atomtype
, char *wall_density
,
599 opts
->wall_atomtype
[0] = NULL
;
600 opts
->wall_atomtype
[1] = NULL
;
602 ir
->wall_atomtype
[0] = -1;
603 ir
->wall_atomtype
[1] = -1;
604 ir
->wall_density
[0] = 0;
605 ir
->wall_density
[1] = 0;
608 nstr
= str_nelem(wall_atomtype
,MAXPTR
,names
);
609 if (nstr
!= ir
->nwall
)
610 gmx_fatal(FARGS
,"Expected %d elements for wall_atomtype, found %d",
612 for(i
=0; i
<ir
->nwall
; i
++)
613 opts
->wall_atomtype
[i
] = strdup(names
[i
]);
615 if (ir
->wall_type
!= ewtTABLE
) {
616 nstr
= str_nelem(wall_density
,MAXPTR
,names
);
617 if (nstr
!= ir
->nwall
)
618 gmx_fatal(FARGS
,"Expected %d elements for wall_density, found %d",
620 for(i
=0; i
<ir
->nwall
; i
++) {
621 sscanf(names
[i
],"%lf",&dbl
);
623 gmx_fatal(FARGS
,"wall_density[%d] = %f\n",i
,dbl
);
624 ir
->wall_density
[i
] = dbl
;
630 static void add_wall_energrps(gmx_groups_t
*groups
,int nwall
,t_symtab
*symtab
)
637 srenew(groups
->grpname
,groups
->ngrpname
+nwall
);
638 grps
= &(groups
->grps
[egcENER
]);
639 srenew(grps
->nm_ind
,grps
->nr
+nwall
);
640 for(i
=0; i
<nwall
; i
++) {
641 sprintf(str
,"wall%d",i
);
642 groups
->grpname
[groups
->ngrpname
] = put_symtab(symtab
,str
);
643 grps
->nm_ind
[grps
->nr
++] = groups
->ngrpname
++;
648 void get_ir(const char *mdparin
,const char *mdparout
,
649 t_inputrec
*ir
,t_gromppopts
*opts
, int *nerror
)
657 inp
=read_inpfile(mdparin
,&ninp
, NULL
);
659 snew(dumstr
[0],STRLEN
);
660 snew(dumstr
[1],STRLEN
);
664 REM_TYPE("domain-decomposition");
665 REPL_TYPE("unconstrained-start","continuation");
666 REM_TYPE("dihre-tau");
667 REM_TYPE("nstdihreout");
668 REM_TYPE("nstcheckpoint");
670 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
671 CTYPE ("Preprocessor information: use cpp syntax.");
672 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
673 STYPE ("include", opts
->include
, NULL
);
674 CTYPE ("e.g.: -DI_Want_Cookies -DMe_Too");
675 STYPE ("define", opts
->define
, NULL
);
677 CCTYPE ("RUN CONTROL PARAMETERS");
678 EETYPE("integrator", ir
->eI
, ei_names
, nerror
, TRUE
);
679 CTYPE ("Start time and timestep in ps");
680 RTYPE ("tinit", ir
->init_t
, 0.0);
681 RTYPE ("dt", ir
->delta_t
, 0.001);
682 STEPTYPE ("nsteps", ir
->nsteps
, 0);
683 CTYPE ("For exact run continuation or redoing part of a run");
684 STEPTYPE ("init_step",ir
->init_step
, 0);
685 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
686 ITYPE ("simulation_part", ir
->simulation_part
, 1);
687 CTYPE ("mode for center of mass motion removal");
688 CTYPE ("energy calculation and T/P-coupling frequency");
689 ITYPE ("nstcalcenergy",ir
->nstcalcenergy
, 1);
690 EETYPE("comm-mode", ir
->comm_mode
, ecm_names
, nerror
, TRUE
);
691 CTYPE ("number of steps for center of mass motion removal");
692 ITYPE ("nstcomm", ir
->nstcomm
, 1);
693 CTYPE ("group(s) for center of mass motion removal");
694 STYPE ("comm-grps", vcm
, NULL
);
696 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
697 CTYPE ("Friction coefficient (amu/ps) and random seed");
698 RTYPE ("bd-fric", ir
->bd_fric
, 0.0);
699 ITYPE ("ld-seed", ir
->ld_seed
, 1993);
702 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
703 CTYPE ("Force tolerance and initial step-size");
704 RTYPE ("emtol", ir
->em_tol
, 10.0);
705 RTYPE ("emstep", ir
->em_stepsize
,0.01);
706 CTYPE ("Max number of iterations in relax_shells");
707 ITYPE ("niter", ir
->niter
, 20);
708 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
709 RTYPE ("fcstep", ir
->fc_stepsize
, 0);
710 CTYPE ("Frequency of steepest descents steps when doing CG");
711 ITYPE ("nstcgsteep", ir
->nstcgsteep
, 1000);
712 ITYPE ("nbfgscorr", ir
->nbfgscorr
, 10);
714 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
715 RTYPE ("rtpi", ir
->rtpi
, 0.05);
718 CCTYPE ("OUTPUT CONTROL OPTIONS");
719 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
720 ITYPE ("nstxout", ir
->nstxout
, 100);
721 ITYPE ("nstvout", ir
->nstvout
, 100);
722 ITYPE ("nstfout", ir
->nstfout
, 0);
723 ir
->nstcheckpoint
= 1000;
724 CTYPE ("Output frequency for energies to log file and energy file");
725 ITYPE ("nstlog", ir
->nstlog
, 100);
726 ITYPE ("nstenergy", ir
->nstenergy
, 100);
727 CTYPE ("Output frequency and precision for xtc file");
728 ITYPE ("nstxtcout", ir
->nstxtcout
, 0);
729 RTYPE ("xtc-precision",ir
->xtcprec
, 1000.0);
730 CTYPE ("This selects the subset of atoms for the xtc file. You can");
731 CTYPE ("select multiple groups. By default all atoms will be written.");
732 STYPE ("xtc-grps", xtc_grps
, NULL
);
733 CTYPE ("Selection of energy groups");
734 STYPE ("energygrps", energy
, NULL
);
736 /* Neighbor searching */
737 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
738 CTYPE ("nblist update frequency");
739 ITYPE ("nstlist", ir
->nstlist
, 10);
740 CTYPE ("ns algorithm (simple or grid)");
741 EETYPE("ns-type", ir
->ns_type
, ens_names
, nerror
, TRUE
);
742 /* set ndelta to the optimal value of 2 */
744 CTYPE ("Periodic boundary conditions: xyz, no, xy");
745 EETYPE("pbc", ir
->ePBC
, epbc_names
, nerror
, TRUE
);
746 EETYPE("periodic_molecules", ir
->bPeriodicMols
, yesno_names
, nerror
, TRUE
);
747 CTYPE ("nblist cut-off");
748 RTYPE ("rlist", ir
->rlist
, 1.0);
749 CTYPE ("long-range cut-off for switched potentials");
750 RTYPE ("rlistlong", ir
->rlistlong
, -1);
753 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
754 CTYPE ("Method for doing electrostatics");
755 EETYPE("coulombtype", ir
->coulombtype
, eel_names
, nerror
, TRUE
);
756 CTYPE ("cut-off lengths");
757 RTYPE ("rcoulomb-switch", ir
->rcoulomb_switch
, 0.0);
758 RTYPE ("rcoulomb", ir
->rcoulomb
, 1.0);
759 CTYPE ("Relative dielectric constant for the medium and the reaction field");
760 RTYPE ("epsilon_r", ir
->epsilon_r
, 1.0);
761 RTYPE ("epsilon_rf", ir
->epsilon_rf
, 1.0);
762 CTYPE ("Method for doing Van der Waals");
763 EETYPE("vdw-type", ir
->vdwtype
, evdw_names
, nerror
, TRUE
);
764 CTYPE ("cut-off lengths");
765 RTYPE ("rvdw-switch", ir
->rvdw_switch
, 0.0);
766 RTYPE ("rvdw", ir
->rvdw
, 1.0);
767 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
768 EETYPE("DispCorr", ir
->eDispCorr
, edispc_names
, nerror
, TRUE
);
769 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
770 RTYPE ("table-extension", ir
->tabext
, 1.0);
771 CTYPE ("Seperate tables between energy group pairs");
772 STYPE ("energygrp_table", egptable
, NULL
);
773 CTYPE ("Spacing for the PME/PPPM FFT grid");
774 RTYPE ("fourierspacing", opts
->fourierspacing
,0.12);
775 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
776 ITYPE ("fourier_nx", ir
->nkx
, 0);
777 ITYPE ("fourier_ny", ir
->nky
, 0);
778 ITYPE ("fourier_nz", ir
->nkz
, 0);
779 CTYPE ("EWALD/PME/PPPM parameters");
780 ITYPE ("pme_order", ir
->pme_order
, 4);
781 RTYPE ("ewald_rtol", ir
->ewald_rtol
, 0.00001);
782 EETYPE("ewald_geometry", ir
->ewald_geometry
, eewg_names
, nerror
, TRUE
);
783 RTYPE ("epsilon_surface", ir
->epsilon_surface
, 0.0);
784 EETYPE("optimize_fft",ir
->bOptFFT
, yesno_names
, nerror
, TRUE
);
786 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
787 EETYPE("implicit_solvent", ir
->implicit_solvent
, eis_names
, nerror
, TRUE
);
789 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
790 CTYPE ("Algorithm for calculating Born radii");
791 EETYPE("gb_algorithm", ir
->gb_algorithm
, egb_names
, nerror
, TRUE
);
792 CTYPE ("Frequency of calculating the Born radii inside rlist");
793 ITYPE ("nstgbradii", ir
->nstgbradii
, 1);
794 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
795 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
796 RTYPE ("rgbradii", ir
->rgbradii
, 1.0);
797 CTYPE ("Dielectric coefficient of the implicit solvent");
798 RTYPE ("gb_epsilon_solvent",ir
->gb_epsilon_solvent
, 80.0);
799 CTYPE ("Salt concentration in M for Generalized Born models");
800 RTYPE ("gb_saltconc", ir
->gb_saltconc
, 0.0);
801 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
802 RTYPE ("gb_obc_alpha", ir
->gb_obc_alpha
, 1.0);
803 RTYPE ("gb_obc_beta", ir
->gb_obc_beta
, 0.8);
804 RTYPE ("gb_obc_gamma", ir
->gb_obc_gamma
, 4.85);
805 RTYPE ("gb_dielectric_offset", ir
->gb_dielectric_offset
, 0.009);
806 EETYPE("sa_algorithm", ir
->sa_algorithm
, esa_names
, nerror
, TRUE
);
807 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
808 CTYPE ("The default value (2.092) corresponds to 0.005 kcal/mol/Angstrom^2.");
809 RTYPE ("sa_surface_tension", ir
->sa_surface_tension
, 2.092);
812 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
813 CTYPE ("Temperature coupling");
814 EETYPE("tcoupl", ir
->etc
, etcoupl_names
, nerror
, TRUE
);
815 ITYPE("nh-chain-length", ir
->opts
.nhchainlength
, NHCHAINLENGTH
);
816 CTYPE ("Groups to couple separately");
817 STYPE ("tc-grps", tcgrps
, NULL
);
818 CTYPE ("Time constant (ps) and reference temperature (K)");
819 STYPE ("tau-t", tau_t
, NULL
);
820 STYPE ("ref-t", ref_t
, NULL
);
821 CTYPE ("Pressure coupling");
822 EETYPE("Pcoupl", ir
->epc
, epcoupl_names
, nerror
, TRUE
);
823 EETYPE("Pcoupltype", ir
->epct
, epcoupltype_names
, nerror
, TRUE
);
824 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
825 RTYPE ("tau-p", ir
->tau_p
, 1.0);
826 STYPE ("compressibility", dumstr
[0], NULL
);
827 STYPE ("ref-p", dumstr
[1], NULL
);
828 CTYPE ("Scaling of reference coordinates, No, All or COM");
829 EETYPE ("refcoord_scaling",ir
->refcoord_scaling
,erefscaling_names
, nerror
, TRUE
);
831 CTYPE ("Random seed for Andersen thermostat");
832 ITYPE ("andersen_seed", ir
->andersen_seed
, 815131);
835 CCTYPE ("OPTIONS FOR QMMM calculations");
836 EETYPE("QMMM", ir
->bQMMM
, yesno_names
, nerror
, TRUE
);
837 CTYPE ("Groups treated Quantum Mechanically");
838 STYPE ("QMMM-grps", QMMM
, NULL
);
840 STYPE("QMmethod", QMmethod
, NULL
);
841 CTYPE ("QMMM scheme");
842 EETYPE("QMMMscheme", ir
->QMMMscheme
, eQMMMscheme_names
, nerror
, TRUE
);
843 CTYPE ("QM basisset");
844 STYPE("QMbasis", QMbasis
, NULL
);
846 STYPE ("QMcharge", QMcharge
,NULL
);
847 CTYPE ("QM multiplicity");
848 STYPE ("QMmult", QMmult
,NULL
);
849 CTYPE ("Surface Hopping");
850 STYPE ("SH", bSH
, NULL
);
851 CTYPE ("CAS space options");
852 STYPE ("CASorbitals", CASorbitals
, NULL
);
853 STYPE ("CASelectrons", CASelectrons
, NULL
);
854 STYPE ("SAon", SAon
, NULL
);
855 STYPE ("SAoff",SAoff
,NULL
);
856 STYPE ("SAsteps", SAsteps
, NULL
);
857 CTYPE ("Scale factor for MM charges");
858 RTYPE ("MMChargeScaleFactor", ir
->scalefactor
, 1.0);
859 CTYPE ("Optimization of QM subsystem");
860 STYPE ("bOPT", bOPT
, NULL
);
861 STYPE ("bTS", bTS
, NULL
);
863 /* Simulated annealing */
864 CCTYPE("SIMULATED ANNEALING");
865 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
866 STYPE ("annealing", anneal
, NULL
);
867 CTYPE ("Number of time points to use for specifying annealing in each group");
868 STYPE ("annealing_npoints", anneal_npoints
, NULL
);
869 CTYPE ("List of times at the annealing points for each group");
870 STYPE ("annealing_time", anneal_time
, NULL
);
871 CTYPE ("Temp. at each annealing point, for each group.");
872 STYPE ("annealing_temp", anneal_temp
, NULL
);
875 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
876 EETYPE("gen-vel", opts
->bGenVel
, yesno_names
, nerror
, TRUE
);
877 RTYPE ("gen-temp", opts
->tempi
, 300.0);
878 ITYPE ("gen-seed", opts
->seed
, 173529);
881 CCTYPE ("OPTIONS FOR BONDS");
882 EETYPE("constraints", opts
->nshake
, constraints
, nerror
, TRUE
);
883 CTYPE ("Type of constraint algorithm");
884 EETYPE("constraint-algorithm", ir
->eConstrAlg
, econstr_names
, nerror
, TRUE
);
885 CTYPE ("Do not constrain the start configuration");
886 EETYPE("continuation", ir
->bContinuation
, yesno_names
, nerror
, TRUE
);
887 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
888 EETYPE("Shake-SOR", ir
->bShakeSOR
, yesno_names
, nerror
, TRUE
);
889 CTYPE ("Relative tolerance of shake");
890 RTYPE ("shake-tol", ir
->shake_tol
, 0.0001);
891 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
892 ITYPE ("lincs-order", ir
->nProjOrder
, 4);
893 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
894 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
895 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
896 ITYPE ("lincs-iter", ir
->nLincsIter
, 1);
897 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
898 CTYPE ("rotates over more degrees than");
899 RTYPE ("lincs-warnangle", ir
->LincsWarnAngle
, 30.0);
900 CTYPE ("Convert harmonic bonds to morse potentials");
901 EETYPE("morse", opts
->bMorse
,yesno_names
, nerror
, TRUE
);
903 /* Energy group exclusions */
904 CCTYPE ("ENERGY GROUP EXCLUSIONS");
905 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
906 STYPE ("energygrp_excl", egpexcl
, NULL
);
910 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
911 ITYPE ("nwall", ir
->nwall
, 0);
912 EETYPE("wall_type", ir
->wall_type
, ewt_names
, nerror
, TRUE
);
913 RTYPE ("wall_r_linpot", ir
->wall_r_linpot
, -1);
914 STYPE ("wall_atomtype", wall_atomtype
, NULL
);
915 STYPE ("wall_density", wall_density
, NULL
);
916 RTYPE ("wall_ewald_zfac", ir
->wall_ewald_zfac
, 3);
919 CCTYPE("COM PULLING");
920 CTYPE("Pull type: no, umbrella, constraint or constant_force");
921 EETYPE("pull", ir
->ePull
, epull_names
, nerror
, TRUE
);
922 if (ir
->ePull
!= epullNO
) {
924 pull_grp
= read_pullparams(&ninp
,&inp
,ir
->pull
,&opts
->pull_start
,nerror
);
927 /* Enforced rotation */
928 CCTYPE("ENFORCED ROTATION");
929 CTYPE("Enforced rotation: No or Yes");
930 EETYPE("rotation", ir
->bRot
, yesno_names
, nerror
, TRUE
);
933 rot_grp
= read_rotparams(&ninp
,&inp
,ir
->rot
);
937 CCTYPE("NMR refinement stuff");
938 CTYPE ("Distance restraints type: No, Simple or Ensemble");
939 EETYPE("disre", ir
->eDisre
, edisre_names
, nerror
, TRUE
);
940 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
941 EETYPE("disre-weighting", ir
->eDisreWeighting
, edisreweighting_names
, nerror
, TRUE
);
942 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
943 EETYPE("disre-mixed", ir
->bDisreMixed
, yesno_names
, nerror
, TRUE
);
944 RTYPE ("disre-fc", ir
->dr_fc
, 1000.0);
945 RTYPE ("disre-tau", ir
->dr_tau
, 0.0);
946 CTYPE ("Output frequency for pair distances to energy file");
947 ITYPE ("nstdisreout", ir
->nstdisreout
, 100);
948 CTYPE ("Orientation restraints: No or Yes");
949 EETYPE("orire", opts
->bOrire
, yesno_names
, nerror
, TRUE
);
950 CTYPE ("Orientation restraints force constant and tau for time averaging");
951 RTYPE ("orire-fc", ir
->orires_fc
, 0.0);
952 RTYPE ("orire-tau", ir
->orires_tau
, 0.0);
953 STYPE ("orire-fitgrp",orirefitgrp
, NULL
);
954 CTYPE ("Output frequency for trace(SD) and S to energy file");
955 ITYPE ("nstorireout", ir
->nstorireout
, 100);
956 CTYPE ("Dihedral angle restraints: No or Yes");
957 EETYPE("dihre", opts
->bDihre
, yesno_names
, nerror
, TRUE
);
958 RTYPE ("dihre-fc", ir
->dihre_fc
, 1000.0);
960 /* Free energy stuff */
961 CCTYPE ("Free energy control stuff");
962 EETYPE("free-energy", ir
->efep
, efep_names
, nerror
, TRUE
);
963 RTYPE ("init-lambda", ir
->init_lambda
,0.0);
964 RTYPE ("delta-lambda",ir
->delta_lambda
,0.0);
965 STYPE ("foreign_lambda", foreign_lambda
, NULL
);
966 RTYPE ("sc-alpha",ir
->sc_alpha
,0.0);
967 ITYPE ("sc-power",ir
->sc_power
,0);
968 RTYPE ("sc-sigma",ir
->sc_sigma
,0.3);
969 ITYPE ("nstdhdl", ir
->nstdhdl
, 10);
970 STYPE ("couple-moltype", couple_moltype
, NULL
);
971 EETYPE("couple-lambda0", opts
->couple_lam0
, couple_lam
, nerror
, TRUE
);
972 EETYPE("couple-lambda1", opts
->couple_lam1
, couple_lam
, nerror
, TRUE
);
973 EETYPE("couple-intramol", opts
->bCoupleIntra
, yesno_names
, nerror
, TRUE
);
975 /* Non-equilibrium MD stuff */
976 CCTYPE("Non-equilibrium MD stuff");
977 STYPE ("acc-grps", accgrps
, NULL
);
978 STYPE ("accelerate", acc
, NULL
);
979 STYPE ("freezegrps", freeze
, NULL
);
980 STYPE ("freezedim", frdim
, NULL
);
981 RTYPE ("cos-acceleration", ir
->cos_accel
, 0);
982 STYPE ("deform", deform
, NULL
);
984 /* Electric fields */
985 CCTYPE("Electric fields");
986 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
987 CTYPE ("and a phase angle (real)");
988 STYPE ("E-x", efield_x
, NULL
);
989 STYPE ("E-xt", efield_xt
, NULL
);
990 STYPE ("E-y", efield_y
, NULL
);
991 STYPE ("E-yt", efield_yt
, NULL
);
992 STYPE ("E-z", efield_z
, NULL
);
993 STYPE ("E-zt", efield_zt
, NULL
);
995 /* User defined thingies */
996 CCTYPE ("User defined thingies");
997 STYPE ("user1-grps", user1
, NULL
);
998 STYPE ("user2-grps", user2
, NULL
);
999 ITYPE ("userint1", ir
->userint1
, 0);
1000 ITYPE ("userint2", ir
->userint2
, 0);
1001 ITYPE ("userint3", ir
->userint3
, 0);
1002 ITYPE ("userint4", ir
->userint4
, 0);
1003 RTYPE ("userreal1", ir
->userreal1
, 0);
1004 RTYPE ("userreal2", ir
->userreal2
, 0);
1005 RTYPE ("userreal3", ir
->userreal3
, 0);
1006 RTYPE ("userreal4", ir
->userreal4
, 0);
1009 write_inpfile(mdparout
,ninp
,inp
,FALSE
);
1010 for (i
=0; (i
<ninp
); i
++) {
1012 sfree(inp
[i
].value
);
1016 /* Process options if necessary */
1017 for(m
=0; m
<2; m
++) {
1018 for(i
=0; i
<2*DIM
; i
++)
1023 if (sscanf(dumstr
[m
],"%lf",&(dumdub
[m
][XX
]))!=1) {
1025 "ERROR: pressure coupling not enough values (I need 1)\n");
1028 dumdub
[m
][YY
]=dumdub
[m
][ZZ
]=dumdub
[m
][XX
];
1030 case epctSEMIISOTROPIC
:
1031 case epctSURFACETENSION
:
1032 if (sscanf(dumstr
[m
],"%lf%lf",
1033 &(dumdub
[m
][XX
]),&(dumdub
[m
][ZZ
]))!=2) {
1035 "ERROR: pressure coupling not enough values (I need 2)\n");
1038 dumdub
[m
][YY
]=dumdub
[m
][XX
];
1040 case epctANISOTROPIC
:
1041 if (sscanf(dumstr
[m
],"%lf%lf%lf%lf%lf%lf",
1042 &(dumdub
[m
][XX
]),&(dumdub
[m
][YY
]),&(dumdub
[m
][ZZ
]),
1043 &(dumdub
[m
][3]),&(dumdub
[m
][4]),&(dumdub
[m
][5]))!=6) {
1045 "ERROR: pressure coupling not enough values (I need 6)\n");
1050 gmx_fatal(FARGS
,"Pressure coupling type %s not implemented yet",
1051 epcoupltype_names
[ir
->epct
]);
1055 clear_mat(ir
->ref_p
);
1056 clear_mat(ir
->compress
);
1057 for(i
=0; i
<DIM
; i
++) {
1058 ir
->ref_p
[i
][i
] = dumdub
[1][i
];
1059 ir
->compress
[i
][i
] = dumdub
[0][i
];
1061 if (ir
->epct
== epctANISOTROPIC
) {
1062 ir
->ref_p
[XX
][YY
] = dumdub
[1][3];
1063 ir
->ref_p
[XX
][ZZ
] = dumdub
[1][4];
1064 ir
->ref_p
[YY
][ZZ
] = dumdub
[1][5];
1065 if (ir
->ref_p
[XX
][YY
]!=0 && ir
->ref_p
[XX
][ZZ
]!=0 && ir
->ref_p
[YY
][ZZ
]!=0) {
1066 warning("All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1068 ir
->compress
[XX
][YY
] = dumdub
[0][3];
1069 ir
->compress
[XX
][ZZ
] = dumdub
[0][4];
1070 ir
->compress
[YY
][ZZ
] = dumdub
[0][5];
1071 for(i
=0; i
<DIM
; i
++) {
1072 for(m
=0; m
<i
; m
++) {
1073 ir
->ref_p
[i
][m
] = ir
->ref_p
[m
][i
];
1074 ir
->compress
[i
][m
] = ir
->compress
[m
][i
];
1079 if (ir
->comm_mode
== ecmNO
)
1082 opts
->couple_moltype
= NULL
;
1083 if (strlen(couple_moltype
) > 0) {
1084 if (ir
->efep
!= efepNO
) {
1085 opts
->couple_moltype
= strdup(couple_moltype
);
1086 if (opts
->couple_lam0
== opts
->couple_lam1
)
1087 warning("The lambda=0 and lambda=1 states for coupling are identical");
1088 if (ir
->eI
== eiMD
&& (opts
->couple_lam0
== ecouplamNONE
||
1089 opts
->couple_lam1
== ecouplamNONE
)) {
1090 warning("For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
1093 warning("Can not couple a molecule with free_energy = no");
1097 do_wall_params(ir
,wall_atomtype
,wall_density
,opts
);
1099 if (opts
->bOrire
&& str_nelem(orirefitgrp
,MAXPTR
,NULL
)!=1) {
1100 fprintf(stderr
,"ERROR: Need one orientation restraint fit group\n");
1104 clear_mat(ir
->deform
);
1107 m
= sscanf(deform
,"%lf %lf %lf %lf %lf %lf",
1108 &(dumdub
[0][0]),&(dumdub
[0][1]),&(dumdub
[0][2]),
1109 &(dumdub
[0][3]),&(dumdub
[0][4]),&(dumdub
[0][5]));
1111 ir
->deform
[i
][i
] = dumdub
[0][i
];
1112 ir
->deform
[YY
][XX
] = dumdub
[0][3];
1113 ir
->deform
[ZZ
][XX
] = dumdub
[0][4];
1114 ir
->deform
[ZZ
][YY
] = dumdub
[0][5];
1115 if (ir
->epc
!= epcNO
) {
1118 if (ir
->deform
[i
][j
]!=0 && ir
->compress
[i
][j
]!=0) {
1119 fprintf(stderr
,"ERROR: A box element has deform set and compressibility > 0\n");
1124 if (ir
->deform
[i
][j
]!=0) {
1125 for(m
=j
; m
<DIM
; m
++)
1126 if (ir
->compress
[m
][j
]!=0) {
1127 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.");
1133 if (ir
->efep
!= efepNO
) {
1134 parse_n_double(foreign_lambda
,&ir
->n_flambda
,&ir
->flambda
);
1135 if (ir
->n_flambda
> 0 && ir
->rlist
< max(ir
->rvdw
,ir
->rcoulomb
)) {
1136 warning_note("For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
1146 static int search_QMstring(char *s
,int ng
,const char *gn
[])
1148 /* same as normal search_string, but this one searches QM strings */
1151 for(i
=0; (i
<ng
); i
++)
1152 if (strcasecmp(s
,gn
[i
]) == 0)
1155 gmx_fatal(FARGS
,"this QM method or basisset (%s) is not implemented\n!",s
);
1159 } /* search_QMstring */
1162 int search_string(char *s
,int ng
,char *gn
[])
1166 for(i
=0; (i
<ng
); i
++)
1167 if (strcasecmp(s
,gn
[i
]) == 0)
1170 gmx_fatal(FARGS
,"Group %s not found in indexfile.\nMaybe you have non-default goups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s
);
1175 static bool do_numbering(int natoms
,gmx_groups_t
*groups
,int ng
,char *ptrs
[],
1176 t_blocka
*block
,char *gnames
[],
1177 int gtype
,int restnm
,
1178 int grptp
,bool bVerbose
)
1180 unsigned short *cbuf
;
1181 t_grps
*grps
=&(groups
->grps
[gtype
]);
1182 int i
,j
,gid
,aj
,ognr
,ntot
=0;
1187 fprintf(debug
,"Starting numbering %d groups of type %d\n",ng
,gtype
);
1189 title
= gtypes
[gtype
];
1192 for(i
=0; (i
<natoms
); i
++)
1195 snew(grps
->nm_ind
,ng
+1); /* +1 for possible rest group */
1196 for(i
=0; (i
<ng
); i
++) {
1197 /* Lookup the group name in the block structure */
1198 gid
= search_string(ptrs
[i
],block
->nr
,gnames
);
1199 if ((grptp
!= egrptpONE
) || (i
== 0))
1200 grps
->nm_ind
[grps
->nr
++]=gid
;
1202 fprintf(debug
,"Found gid %d for group %s\n",gid
,ptrs
[i
]);
1204 /* Now go over the atoms in the group */
1205 for(j
=block
->index
[gid
]; (j
<block
->index
[gid
+1]); j
++) {
1208 /* Range checking */
1209 if ((aj
< 0) || (aj
>= natoms
))
1210 gmx_fatal(FARGS
,"Invalid atom number %d in indexfile",aj
);
1212 /* Lookup up the old group number */
1215 gmx_fatal(FARGS
,"Atom %d in multiple %s groups (%d and %d)",
1216 aj
+1,title
,ognr
+1,i
+1);
1218 /* Store the group number in buffer */
1219 if (grptp
== egrptpONE
)
1228 /* Now check whether we have done all atoms */
1230 if (ntot
!= natoms
) {
1231 if (grptp
== egrptpALL
) {
1232 gmx_fatal(FARGS
,"%d atoms are not part of any of the %s groups",
1234 } else if (grptp
== egrptpPART
) {
1235 sprintf(warn_buf
,"%d atoms are not part of any of the %s groups",
1239 /* Assign all atoms currently unassigned to a rest group */
1240 for(j
=0; (j
<natoms
); j
++) {
1241 if (cbuf
[j
] == NOGID
) {
1246 if (grptp
!= egrptpPART
) {
1249 "Making dummy/rest group for %s containing %d elements\n",
1251 /* Add group name "rest" */
1252 grps
->nm_ind
[grps
->nr
] = restnm
;
1254 /* Assign the rest name to all atoms not currently assigned to a group */
1255 for(j
=0; (j
<natoms
); j
++) {
1256 if (cbuf
[j
] == NOGID
)
1263 if (grps
->nr
== 1) {
1264 groups
->ngrpnr
[gtype
] = 0;
1265 groups
->grpnr
[gtype
] = NULL
;
1267 groups
->ngrpnr
[gtype
] = natoms
;
1268 snew(groups
->grpnr
[gtype
],natoms
);
1269 for(j
=0; (j
<natoms
); j
++) {
1270 groups
->grpnr
[gtype
][j
] = cbuf
[j
];
1276 return (bRest
&& grptp
== egrptpPART
);
1279 static void calc_nrdf(gmx_mtop_t
*mtop
,t_inputrec
*ir
,char **gnames
)
1282 gmx_groups_t
*groups
;
1284 int natoms
,ai
,aj
,i
,j
,d
,g
,imin
,jmin
,nc
;
1286 int *nrdf
,*na_vcm
,na_tot
;
1287 double *nrdf_tc
,*nrdf_vcm
,nrdf_uc
,n_sub
=0;
1288 gmx_mtop_atomloop_all_t aloop
;
1290 int mb
,mol
,ftype
,as
;
1291 gmx_molblock_t
*molb
;
1292 gmx_moltype_t
*molt
;
1295 * First calc 3xnr-atoms for each group
1296 * then subtract half a degree of freedom for each constraint
1298 * Only atoms and nuclei contribute to the degrees of freedom...
1303 groups
= &mtop
->groups
;
1304 natoms
= mtop
->natoms
;
1306 /* Allocate one more for a possible rest group */
1307 /* We need to sum degrees of freedom into doubles,
1308 * since floats give too low nrdf's above 3 million atoms.
1310 snew(nrdf_tc
,groups
->grps
[egcTC
].nr
+1);
1311 snew(nrdf_vcm
,groups
->grps
[egcVCM
].nr
+1);
1312 snew(na_vcm
,groups
->grps
[egcVCM
].nr
+1);
1314 for(i
=0; i
<groups
->grps
[egcTC
].nr
; i
++)
1316 for(i
=0; i
<groups
->grps
[egcVCM
].nr
+1; i
++)
1320 aloop
= gmx_mtop_atomloop_all_init(mtop
);
1321 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
1323 if (atom
->ptype
== eptAtom
|| atom
->ptype
== eptNucleus
) {
1324 g
= ggrpnr(groups
,egcFREEZE
,i
);
1325 /* Double count nrdf for particle i */
1326 for(d
=0; d
<DIM
; d
++) {
1327 if (opts
->nFreeze
[g
][d
] == 0) {
1331 nrdf_tc
[ggrpnr(groups
,egcTC
,i
)] += 0.5*nrdf
[i
];
1332 nrdf_vcm
[ggrpnr(groups
,egcVCM
,i
)] += 0.5*nrdf
[i
];
1337 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
1338 molb
= &mtop
->molblock
[mb
];
1339 molt
= &mtop
->moltype
[molb
->type
];
1340 atom
= molt
->atoms
.atom
;
1341 for(mol
=0; mol
<molb
->nmol
; mol
++) {
1342 for (ftype
=F_CONSTR
; ftype
<=F_CONSTRNC
; ftype
++) {
1343 ia
= molt
->ilist
[ftype
].iatoms
;
1344 for(i
=0; i
<molt
->ilist
[ftype
].nr
; ) {
1345 /* Subtract degrees of freedom for the constraints,
1346 * if the particles still have degrees of freedom left.
1347 * If one of the particles is a vsite or a shell, then all
1348 * constraint motion will go there, but since they do not
1349 * contribute to the constraints the degrees of freedom do not
1354 if (((atom
[ia
[1]].ptype
== eptNucleus
) ||
1355 (atom
[ia
[1]].ptype
== eptAtom
)) &&
1356 ((atom
[ia
[2]].ptype
== eptNucleus
) ||
1357 (atom
[ia
[2]].ptype
== eptAtom
))) {
1366 imin
= min(imin
,nrdf
[ai
]);
1367 jmin
= min(jmin
,nrdf
[aj
]);
1370 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5*imin
;
1371 nrdf_tc
[ggrpnr(groups
,egcTC
,aj
)] -= 0.5*jmin
;
1372 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5*imin
;
1373 nrdf_vcm
[ggrpnr(groups
,egcVCM
,aj
)] -= 0.5*jmin
;
1375 ia
+= interaction_function
[ftype
].nratoms
+1;
1376 i
+= interaction_function
[ftype
].nratoms
+1;
1379 ia
= molt
->ilist
[F_SETTLE
].iatoms
;
1380 for(i
=0; i
<molt
->ilist
[F_SETTLE
].nr
; ) {
1381 /* Subtract 1 dof from every atom in the SETTLE */
1382 for(ai
=as
+ia
[1]; ai
<as
+ia
[1]+3; ai
++) {
1383 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 1;
1384 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 1;
1389 as
+= molt
->atoms
.nr
;
1393 if (ir
->ePull
== epullCONSTRAINT
) {
1394 /* Correct nrdf for the COM constraints.
1395 * We correct using the TC and VCM group of the first atom
1396 * in the reference and pull group. If atoms in one pull group
1397 * belong to different TC or VCM groups it is anyhow difficult
1398 * to determine the optimal nrdf assignment.
1401 if (pull
->eGeom
== epullgPOS
) {
1403 for(i
=0; i
<DIM
; i
++) {
1410 for(i
=0; i
<pull
->ngrp
; i
++) {
1412 if (pull
->grp
[0].nat
> 0) {
1413 /* Subtract 1/2 dof from the reference group */
1414 ai
= pull
->grp
[0].ind
[0];
1415 if (nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] > 1) {
1416 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5;
1417 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5;
1421 /* Subtract 1/2 dof from the pulled group */
1422 ai
= pull
->grp
[1+i
].ind
[0];
1423 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5*imin
;
1424 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5*imin
;
1425 if (nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] < 0)
1426 gmx_fatal(FARGS
,"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
)]]);
1430 if (ir
->nstcomm
!= 0) {
1431 /* Subtract 3 from the number of degrees of freedom in each vcm group
1432 * when com translation is removed and 6 when rotation is removed
1435 switch (ir
->comm_mode
) {
1437 n_sub
= ndof_com(ir
);
1444 gmx_incons("Checking comm_mode");
1447 for(i
=0; i
<groups
->grps
[egcTC
].nr
; i
++) {
1448 /* Count the number of atoms of TC group i for every VCM group */
1449 for(j
=0; j
<groups
->grps
[egcVCM
].nr
+1; j
++)
1452 for(ai
=0; ai
<natoms
; ai
++)
1453 if (ggrpnr(groups
,egcTC
,ai
) == i
) {
1454 na_vcm
[ggrpnr(groups
,egcVCM
,ai
)]++;
1457 /* Correct for VCM removal according to the fraction of each VCM
1458 * group present in this TC group.
1460 nrdf_uc
= nrdf_tc
[i
];
1462 fprintf(debug
,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
1466 for(j
=0; j
<groups
->grps
[egcVCM
].nr
+1; j
++) {
1467 if (nrdf_vcm
[j
] > n_sub
) {
1468 nrdf_tc
[i
] += nrdf_uc
*((double)na_vcm
[j
]/(double)na_tot
)*
1469 (nrdf_vcm
[j
] - n_sub
)/nrdf_vcm
[j
];
1472 fprintf(debug
," nrdf_vcm[%d] = %g, nrdf = %g\n",
1473 j
,nrdf_vcm
[j
],nrdf_tc
[i
]);
1478 for(i
=0; (i
<groups
->grps
[egcTC
].nr
); i
++) {
1479 opts
->nrdf
[i
] = nrdf_tc
[i
];
1480 if (opts
->nrdf
[i
] < 0)
1483 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
1484 gnames
[groups
->grps
[egcTC
].nm_ind
[i
]],opts
->nrdf
[i
]);
1493 static void decode_cos(char *s
,t_cosines
*cosine
,bool bTime
)
1496 char format
[STRLEN
],f1
[STRLEN
];
1507 sscanf(t
,"%d",&(cosine
->n
));
1508 if (cosine
->n
<= 0) {
1511 snew(cosine
->a
,cosine
->n
);
1512 snew(cosine
->phi
,cosine
->n
);
1514 sprintf(format
,"%%*d");
1515 for(i
=0; (i
<cosine
->n
); i
++) {
1517 strcat(f1
,"%lf%lf");
1518 if (sscanf(t
,f1
,&a
,&phi
) < 2)
1519 gmx_fatal(FARGS
,"Invalid input for electric field shift: '%s'",t
);
1522 strcat(format
,"%*lf%*lf");
1529 static bool do_egp_flag(t_inputrec
*ir
,gmx_groups_t
*groups
,
1530 const char *option
,const char *val
,int flag
)
1532 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
1533 * But since this is much larger than STRLEN, such a line can not be parsed.
1534 * The real maximum is the number of names that fit in a string: STRLEN/2.
1536 #define EGP_MAX (STRLEN/2)
1538 char *names
[EGP_MAX
];
1542 gnames
= groups
->grpname
;
1544 nelem
= str_nelem(val
,EGP_MAX
,names
);
1546 gmx_fatal(FARGS
,"The number of groups for %s is odd",option
);
1547 nr
= groups
->grps
[egcENER
].nr
;
1549 for(i
=0; i
<nelem
/2; i
++) {
1552 strcasecmp(names
[2*i
],*(gnames
[groups
->grps
[egcENER
].nm_ind
[j
]])))
1555 gmx_fatal(FARGS
,"%s in %s is not an energy group\n",
1559 strcasecmp(names
[2*i
+1],*(gnames
[groups
->grps
[egcENER
].nm_ind
[k
]])))
1562 gmx_fatal(FARGS
,"%s in %s is not an energy group\n",
1563 names
[2*i
+1],option
);
1564 if ((j
< nr
) && (k
< nr
)) {
1565 ir
->opts
.egp_flags
[nr
*j
+k
] |= flag
;
1566 ir
->opts
.egp_flags
[nr
*k
+j
] |= flag
;
1574 void do_index(const char* mdparin
, const char *ndx
,
1577 t_inputrec
*ir
,rvec
*v
)
1580 gmx_groups_t
*groups
;
1584 char warnbuf
[STRLEN
],**gnames
;
1585 int nr
,ntcg
,ntau_t
,nref_t
,nacc
,nofg
,nSA
,nSA_points
,nSA_time
,nSA_temp
;
1586 int nacg
,nfreeze
,nfrdim
,nenergy
,nvcm
,nuser
;
1587 char *ptr1
[MAXPTR
],*ptr2
[MAXPTR
],*ptr3
[MAXPTR
];
1590 bool bExcl
,bTable
,bSetTCpar
,bAnneal
,bRest
;
1591 int nQMmethod
,nQMbasis
,nQMcharge
,nQMmult
,nbSH
,nCASorb
,nCASelec
,
1592 nSAon
,nSAoff
,nSAsteps
,nQMg
,nbOPT
,nbTS
;
1595 fprintf(stderr
,"processing index file...\n");
1599 snew(grps
->index
,1);
1601 atoms_all
= gmx_mtop_global_atoms(mtop
);
1602 analyse(&atoms_all
,grps
,&gnames
,FALSE
,TRUE
);
1603 free_t_atoms(&atoms_all
,FALSE
);
1605 grps
= init_index(ndx
,&gnames
);
1608 groups
= &mtop
->groups
;
1609 natoms
= mtop
->natoms
;
1610 symtab
= &mtop
->symtab
;
1612 snew(groups
->grpname
,grps
->nr
+1);
1614 for(i
=0; (i
<grps
->nr
); i
++) {
1615 groups
->grpname
[i
] = put_symtab(symtab
,gnames
[i
]);
1617 groups
->grpname
[i
] = put_symtab(symtab
,"rest");
1619 srenew(gnames
,grps
->nr
+1);
1620 gnames
[restnm
] = *(groups
->grpname
[i
]);
1621 groups
->ngrpname
= grps
->nr
+1;
1623 set_warning_line(mdparin
,-1);
1625 ntau_t
= str_nelem(tau_t
,MAXPTR
,ptr1
);
1626 nref_t
= str_nelem(ref_t
,MAXPTR
,ptr2
);
1627 ntcg
= str_nelem(tcgrps
,MAXPTR
,ptr3
);
1628 if ((ntau_t
!= ntcg
) || (nref_t
!= ntcg
)) {
1629 gmx_fatal(FARGS
,"Invalid T coupling input: %d groups, %d ref_t values and "
1630 "%d tau_t values",ntcg
,nref_t
,ntau_t
);
1633 bSetTCpar
= (ir
->etc
|| EI_SD(ir
->eI
) || ir
->eI
==eiBD
|| EI_TPI(ir
->eI
));
1634 do_numbering(natoms
,groups
,ntcg
,ptr3
,grps
,gnames
,egcTC
,
1635 restnm
,bSetTCpar
? egrptpALL
: egrptpALL_GENREST
,bVerbose
);
1636 nr
= groups
->grps
[egcTC
].nr
;
1638 snew(ir
->opts
.nrdf
,nr
);
1639 snew(ir
->opts
.tau_t
,nr
);
1640 snew(ir
->opts
.ref_t
,nr
);
1641 if (ir
->eI
==eiBD
&& ir
->bd_fric
==0) {
1642 fprintf(stderr
,"bd_fric=0, so tau_t will be used as the inverse friction constant(s)\n");
1646 gmx_fatal(FARGS
,"Not enough ref_t and tau_t values!");
1647 for(i
=0; (i
<nr
); i
++) {
1648 ir
->opts
.tau_t
[i
]=strtod(ptr1
[i
],NULL
);
1649 if (ir
->opts
.tau_t
[i
] < 0) {
1650 gmx_fatal(FARGS
,"tau_t for group %d negative",i
);
1652 /* We check the relative magnitude of the coupling time tau_t.
1653 * V-rescale works correctly, even for tau_t=0.
1655 if ((ir
->etc
== etcBERENDSEN
|| ir
->etc
== etcNOSEHOOVER
) &&
1656 ir
->opts
.tau_t
[i
] != 0 &&
1657 ir
->opts
.tau_t
[i
] < 10*ir
->nstcalcenergy
*ir
->delta_t
) {
1658 sprintf(warn_buf
,"For proper thermostat integration tau_t (%g) should be more than an order of magnitude larger than nstcalcenergy*dt (%g)",
1659 ir
->opts
.tau_t
[i
],ir
->nstcalcenergy
*ir
->delta_t
);
1663 for(i
=0; (i
<nr
); i
++) {
1664 ir
->opts
.ref_t
[i
]=strtod(ptr2
[i
],NULL
);
1665 if (ir
->opts
.ref_t
[i
] < 0)
1666 gmx_fatal(FARGS
,"ref_t for group %d negative",i
);
1670 /* Simulated annealing for each group. There are nr groups */
1671 nSA
= str_nelem(anneal
,MAXPTR
,ptr1
);
1672 if (nSA
== 1 && (ptr1
[0][0]=='n' || ptr1
[0][0]=='N'))
1674 if(nSA
>0 && nSA
!= nr
)
1675 gmx_fatal(FARGS
,"Not enough annealing values: %d (for %d groups)\n",nSA
,nr
);
1677 snew(ir
->opts
.annealing
,nr
);
1678 snew(ir
->opts
.anneal_npoints
,nr
);
1679 snew(ir
->opts
.anneal_time
,nr
);
1680 snew(ir
->opts
.anneal_temp
,nr
);
1682 ir
->opts
.annealing
[i
]=eannNO
;
1683 ir
->opts
.anneal_npoints
[i
]=0;
1684 ir
->opts
.anneal_time
[i
]=NULL
;
1685 ir
->opts
.anneal_temp
[i
]=NULL
;
1690 if(ptr1
[i
][0]=='n' || ptr1
[i
][0]=='N') {
1691 ir
->opts
.annealing
[i
]=eannNO
;
1692 } else if(ptr1
[i
][0]=='s'|| ptr1
[i
][0]=='S') {
1693 ir
->opts
.annealing
[i
]=eannSINGLE
;
1695 } else if(ptr1
[i
][0]=='p'|| ptr1
[i
][0]=='P') {
1696 ir
->opts
.annealing
[i
]=eannPERIODIC
;
1701 /* Read the other fields too */
1702 nSA_points
= str_nelem(anneal_npoints
,MAXPTR
,ptr1
);
1704 gmx_fatal(FARGS
,"Found %d annealing_npoints values for %d groups\n",nSA_points
,nSA
);
1705 for(k
=0,i
=0;i
<nr
;i
++) {
1706 ir
->opts
.anneal_npoints
[i
]=strtol(ptr1
[i
],NULL
,10);
1707 if(ir
->opts
.anneal_npoints
[i
]==1)
1708 gmx_fatal(FARGS
,"Please specify at least a start and an end point for annealing\n");
1709 snew(ir
->opts
.anneal_time
[i
],ir
->opts
.anneal_npoints
[i
]);
1710 snew(ir
->opts
.anneal_temp
[i
],ir
->opts
.anneal_npoints
[i
]);
1711 k
+= ir
->opts
.anneal_npoints
[i
];
1714 nSA_time
= str_nelem(anneal_time
,MAXPTR
,ptr1
);
1716 gmx_fatal(FARGS
,"Found %d annealing_time values, wanter %d\n",nSA_time
,k
);
1717 nSA_temp
= str_nelem(anneal_temp
,MAXPTR
,ptr2
);
1719 gmx_fatal(FARGS
,"Found %d annealing_temp values, wanted %d\n",nSA_temp
,k
);
1721 for(i
=0,k
=0;i
<nr
;i
++) {
1723 for(j
=0;j
<ir
->opts
.anneal_npoints
[i
];j
++) {
1724 ir
->opts
.anneal_time
[i
][j
]=strtod(ptr1
[k
],NULL
);
1725 ir
->opts
.anneal_temp
[i
][j
]=strtod(ptr2
[k
],NULL
);
1727 if(ir
->opts
.anneal_time
[i
][0] > (ir
->init_t
+GMX_REAL_EPS
))
1728 gmx_fatal(FARGS
,"First time point for annealing > init_t.\n");
1731 if(ir
->opts
.anneal_time
[i
][j
]<ir
->opts
.anneal_time
[i
][j
-1])
1732 gmx_fatal(FARGS
,"Annealing timepoints out of order: t=%f comes after t=%f\n",
1733 ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_time
[i
][j
-1]);
1735 if(ir
->opts
.anneal_temp
[i
][j
]<0)
1736 gmx_fatal(FARGS
,"Found negative temperature in annealing: %f\n",ir
->opts
.anneal_temp
[i
][j
]);
1740 /* Print out some summary information, to make sure we got it right */
1741 for(i
=0,k
=0;i
<nr
;i
++) {
1742 if(ir
->opts
.annealing
[i
]!=eannNO
) {
1743 j
= groups
->grps
[egcTC
].nm_ind
[i
];
1744 fprintf(stderr
,"Simulated annealing for group %s: %s, %d timepoints\n",
1745 *(groups
->grpname
[j
]),eann_names
[ir
->opts
.annealing
[i
]],
1746 ir
->opts
.anneal_npoints
[i
]);
1747 fprintf(stderr
,"Time (ps) Temperature (K)\n");
1748 /* All terms except the last one */
1749 for(j
=0;j
<(ir
->opts
.anneal_npoints
[i
]-1);j
++)
1750 fprintf(stderr
,"%9.1f %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1752 /* Finally the last one */
1753 j
= ir
->opts
.anneal_npoints
[i
]-1;
1754 if(ir
->opts
.annealing
[i
]==eannSINGLE
)
1755 fprintf(stderr
,"%9.1f- %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1757 fprintf(stderr
,"%9.1f %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1758 if(fabs(ir
->opts
.anneal_temp
[i
][j
]-ir
->opts
.anneal_temp
[i
][0])>GMX_REAL_EPS
)
1759 warning_note("There is a temperature jump when your annealing loops back.\n");
1767 if (ir
->ePull
!= epullNO
) {
1768 make_pull_groups(ir
->pull
,pull_grp
,grps
,gnames
);
1772 make_rotation_groups(ir
->rot
,rot_grp
,grps
,gnames
);
1775 nacc
= str_nelem(acc
,MAXPTR
,ptr1
);
1776 nacg
= str_nelem(accgrps
,MAXPTR
,ptr2
);
1777 if (nacg
*DIM
!= nacc
)
1778 gmx_fatal(FARGS
,"Invalid Acceleration input: %d groups and %d acc. values",
1780 do_numbering(natoms
,groups
,nacg
,ptr2
,grps
,gnames
,egcACC
,
1781 restnm
,egrptpALL_GENREST
,bVerbose
);
1782 nr
= groups
->grps
[egcACC
].nr
;
1783 snew(ir
->opts
.acc
,nr
);
1786 for(i
=k
=0; (i
<nacg
); i
++)
1787 for(j
=0; (j
<DIM
); j
++,k
++)
1788 ir
->opts
.acc
[i
][j
]=strtod(ptr1
[k
],NULL
);
1790 for(j
=0; (j
<DIM
); j
++)
1791 ir
->opts
.acc
[i
][j
]=0;
1793 nfrdim
= str_nelem(frdim
,MAXPTR
,ptr1
);
1794 nfreeze
= str_nelem(freeze
,MAXPTR
,ptr2
);
1795 if (nfrdim
!= DIM
*nfreeze
)
1796 gmx_fatal(FARGS
,"Invalid Freezing input: %d groups and %d freeze values",
1798 do_numbering(natoms
,groups
,nfreeze
,ptr2
,grps
,gnames
,egcFREEZE
,
1799 restnm
,egrptpALL_GENREST
,bVerbose
);
1800 nr
= groups
->grps
[egcFREEZE
].nr
;
1802 snew(ir
->opts
.nFreeze
,nr
);
1803 for(i
=k
=0; (i
<nfreeze
); i
++)
1804 for(j
=0; (j
<DIM
); j
++,k
++) {
1805 ir
->opts
.nFreeze
[i
][j
]=(strncasecmp(ptr1
[k
],"Y",1)==0);
1806 if (!ir
->opts
.nFreeze
[i
][j
]) {
1807 if (strncasecmp(ptr1
[k
],"N",1) != 0) {
1808 sprintf(warnbuf
,"Please use Y(ES) or N(O) for freezedim only "
1809 "(not %s)", ptr1
[k
]);
1815 for(j
=0; (j
<DIM
); j
++)
1816 ir
->opts
.nFreeze
[i
][j
]=0;
1818 nenergy
=str_nelem(energy
,MAXPTR
,ptr1
);
1819 do_numbering(natoms
,groups
,nenergy
,ptr1
,grps
,gnames
,egcENER
,
1820 restnm
,egrptpALL_GENREST
,bVerbose
);
1821 add_wall_energrps(groups
,ir
->nwall
,symtab
);
1822 ir
->opts
.ngener
= groups
->grps
[egcENER
].nr
;
1823 nvcm
=str_nelem(vcm
,MAXPTR
,ptr1
);
1825 do_numbering(natoms
,groups
,nvcm
,ptr1
,grps
,gnames
,egcVCM
,
1826 restnm
,nvcm
==0 ? egrptpALL_GENREST
: egrptpPART
,bVerbose
);
1828 warning("Some atoms are not part of any center of mass motion removal group.\n"
1829 "This may lead to artifacts.\n"
1830 "In most cases one should use one group for the whole system.");
1833 /* Now we have filled the freeze struct, so we can calculate NRDF */
1834 calc_nrdf(mtop
,ir
,gnames
);
1839 /* Must check per group! */
1840 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1841 ntot
+= ir
->opts
.nrdf
[i
];
1842 if (ntot
!= (DIM
*natoms
)) {
1843 fac
= sqrt(ntot
/(DIM
*natoms
));
1845 fprintf(stderr
,"Scaling velocities by a factor of %.3f to account for constraints\n"
1846 "and removal of center of mass motion\n",fac
);
1847 for(i
=0; (i
<natoms
); i
++)
1848 svmul(fac
,v
[i
],v
[i
]);
1852 nuser
=str_nelem(user1
,MAXPTR
,ptr1
);
1853 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcUser1
,
1854 restnm
,egrptpALL_GENREST
,bVerbose
);
1855 nuser
=str_nelem(user2
,MAXPTR
,ptr1
);
1856 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcUser2
,
1857 restnm
,egrptpALL_GENREST
,bVerbose
);
1858 nuser
=str_nelem(xtc_grps
,MAXPTR
,ptr1
);
1859 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcXTC
,
1860 restnm
,egrptpONE
,bVerbose
);
1861 nofg
= str_nelem(orirefitgrp
,MAXPTR
,ptr1
);
1862 do_numbering(natoms
,groups
,nofg
,ptr1
,grps
,gnames
,egcORFIT
,
1863 restnm
,egrptpALL_GENREST
,bVerbose
);
1865 /* QMMM input processing */
1866 nQMg
= str_nelem(QMMM
,MAXPTR
,ptr1
);
1867 nQMmethod
= str_nelem(QMmethod
,MAXPTR
,ptr2
);
1868 nQMbasis
= str_nelem(QMbasis
,MAXPTR
,ptr3
);
1869 if((nQMmethod
!= nQMg
)||(nQMbasis
!= nQMg
)){
1870 gmx_fatal(FARGS
,"Invalid QMMM input: %d groups %d basissets"
1871 " and %d methods\n",nQMg
,nQMbasis
,nQMmethod
);
1873 /* group rest, if any, is always MM! */
1874 do_numbering(natoms
,groups
,nQMg
,ptr1
,grps
,gnames
,egcQMMM
,
1875 restnm
,egrptpALL_GENREST
,bVerbose
);
1876 nr
= nQMg
; /*atoms->grps[egcQMMM].nr;*/
1877 ir
->opts
.ngQM
= nQMg
;
1878 snew(ir
->opts
.QMmethod
,nr
);
1879 snew(ir
->opts
.QMbasis
,nr
);
1881 /* input consists of strings: RHF CASSCF PM3 .. These need to be
1882 * converted to the corresponding enum in names.c
1884 ir
->opts
.QMmethod
[i
] = search_QMstring(ptr2
[i
],eQMmethodNR
,
1886 ir
->opts
.QMbasis
[i
] = search_QMstring(ptr3
[i
],eQMbasisNR
,
1890 nQMmult
= str_nelem(QMmult
,MAXPTR
,ptr1
);
1891 nQMcharge
= str_nelem(QMcharge
,MAXPTR
,ptr2
);
1892 nbSH
= str_nelem(bSH
,MAXPTR
,ptr3
);
1893 snew(ir
->opts
.QMmult
,nr
);
1894 snew(ir
->opts
.QMcharge
,nr
);
1895 snew(ir
->opts
.bSH
,nr
);
1898 ir
->opts
.QMmult
[i
] = strtol(ptr1
[i
],NULL
,10);
1899 ir
->opts
.QMcharge
[i
] = strtol(ptr2
[i
],NULL
,10);
1900 ir
->opts
.bSH
[i
] = (strncasecmp(ptr3
[i
],"Y",1)==0);
1903 nCASelec
= str_nelem(CASelectrons
,MAXPTR
,ptr1
);
1904 nCASorb
= str_nelem(CASorbitals
,MAXPTR
,ptr2
);
1905 snew(ir
->opts
.CASelectrons
,nr
);
1906 snew(ir
->opts
.CASorbitals
,nr
);
1908 ir
->opts
.CASelectrons
[i
]= strtol(ptr1
[i
],NULL
,10);
1909 ir
->opts
.CASorbitals
[i
] = strtol(ptr2
[i
],NULL
,10);
1911 /* special optimization options */
1913 nbOPT
= str_nelem(bOPT
,MAXPTR
,ptr1
);
1914 nbTS
= str_nelem(bTS
,MAXPTR
,ptr2
);
1915 snew(ir
->opts
.bOPT
,nr
);
1916 snew(ir
->opts
.bTS
,nr
);
1918 ir
->opts
.bOPT
[i
] = (strncasecmp(ptr1
[i
],"Y",1)==0);
1919 ir
->opts
.bTS
[i
] = (strncasecmp(ptr2
[i
],"Y",1)==0);
1921 nSAon
= str_nelem(SAon
,MAXPTR
,ptr1
);
1922 nSAoff
= str_nelem(SAoff
,MAXPTR
,ptr2
);
1923 nSAsteps
= str_nelem(SAsteps
,MAXPTR
,ptr3
);
1924 snew(ir
->opts
.SAon
,nr
);
1925 snew(ir
->opts
.SAoff
,nr
);
1926 snew(ir
->opts
.SAsteps
,nr
);
1929 ir
->opts
.SAon
[i
] = strtod(ptr1
[i
],NULL
);
1930 ir
->opts
.SAoff
[i
] = strtod(ptr2
[i
],NULL
);
1931 ir
->opts
.SAsteps
[i
] = strtol(ptr3
[i
],NULL
,10);
1933 /* end of QMMM input */
1936 for(i
=0; (i
<egcNR
); i
++) {
1937 fprintf(stderr
,"%-16s has %d element(s):",gtypes
[i
],groups
->grps
[i
].nr
);
1938 for(j
=0; (j
<groups
->grps
[i
].nr
); j
++)
1939 fprintf(stderr
," %s",*(groups
->grpname
[groups
->grps
[i
].nm_ind
[j
]]));
1940 fprintf(stderr
,"\n");
1943 nr
= groups
->grps
[egcENER
].nr
;
1944 snew(ir
->opts
.egp_flags
,nr
*nr
);
1946 bExcl
= do_egp_flag(ir
,groups
,"energygrp_excl",egpexcl
,EGP_EXCL
);
1947 if (bExcl
&& EEL_FULL(ir
->coulombtype
))
1948 warning("Can not exclude the lattice Coulomb energy between energy groups");
1950 bTable
= do_egp_flag(ir
,groups
,"energygrp_table",egptable
,EGP_TABLE
);
1951 if (bTable
&& !(ir
->vdwtype
== evdwUSER
) &&
1952 !(ir
->coulombtype
== eelUSER
) &&!(ir
->coulombtype
== eelPMEUSER
))
1953 gmx_fatal(FARGS
,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
1955 decode_cos(efield_x
,&(ir
->ex
[XX
]),FALSE
);
1956 decode_cos(efield_xt
,&(ir
->et
[XX
]),TRUE
);
1957 decode_cos(efield_y
,&(ir
->ex
[YY
]),FALSE
);
1958 decode_cos(efield_yt
,&(ir
->et
[YY
]),TRUE
);
1959 decode_cos(efield_z
,&(ir
->ex
[ZZ
]),FALSE
);
1960 decode_cos(efield_zt
,&(ir
->et
[ZZ
]),TRUE
);
1962 for(i
=0; (i
<grps
->nr
); i
++)
1972 static void check_disre(gmx_mtop_t
*mtop
)
1974 gmx_ffparams_t
*ffparams
;
1975 t_functype
*functype
;
1977 int i
,ndouble
,ftype
;
1978 int label
,old_label
;
1980 if (gmx_mtop_ftype_count(mtop
,F_DISRES
) > 0) {
1981 ffparams
= &mtop
->ffparams
;
1982 functype
= ffparams
->functype
;
1983 ip
= ffparams
->iparams
;
1986 for(i
=0; i
<ffparams
->ntypes
; i
++) {
1987 ftype
= functype
[i
];
1988 if (ftype
== F_DISRES
) {
1989 label
= ip
[i
].disres
.label
;
1990 if (label
== old_label
) {
1991 fprintf(stderr
,"Distance restraint index %d occurs twice\n",label
);
1998 gmx_fatal(FARGS
,"Found %d double distance restraint indices,\n"
1999 "probably the parameters for multiple pairs in one restraint "
2000 "are not identical\n",ndouble
);
2004 static bool absolute_reference(t_inputrec
*ir
,gmx_mtop_t
*sys
,ivec AbsRef
)
2007 gmx_mtop_ilistloop_t iloop
;
2013 for(d
=0; d
<DIM
; d
++) {
2014 AbsRef
[d
] = (d
< ndof_com(ir
) ? 0 : 1);
2016 /* Check for freeze groups */
2017 for(g
=0; g
<ir
->opts
.ngfrz
; g
++) {
2018 for(d
=0; d
<DIM
; d
++) {
2019 if (ir
->opts
.nFreeze
[g
][d
] != 0) {
2024 /* Check for position restraints */
2025 iloop
= gmx_mtop_ilistloop_init(sys
);
2026 while (gmx_mtop_ilistloop_next(iloop
,&ilist
,&nmol
)) {
2028 for(i
=0; i
<ilist
[F_POSRES
].nr
; i
+=2) {
2029 pr
= &sys
->ffparams
.iparams
[ilist
[F_POSRES
].iatoms
[i
]];
2030 for(d
=0; d
<DIM
; d
++) {
2031 if (pr
->posres
.fcA
[d
] != 0) {
2039 return (AbsRef
[XX
] != 0 && AbsRef
[YY
] != 0 && AbsRef
[ZZ
] != 0);
2042 void triple_check(const char *mdparin
,t_inputrec
*ir
,gmx_mtop_t
*sys
,int *nerror
)
2045 int i
,m
,g
,nmol
,npct
;
2047 real gdt_max
,*mgrp
,mt
;
2049 gmx_mtop_atomloop_block_t aloopb
;
2050 gmx_mtop_atomloop_all_t aloop
;
2054 if (EI_DYNAMICS(ir
->eI
) && !EI_SD(ir
->eI
) && ir
->eI
!= eiBD
&&
2055 ir
->comm_mode
== ecmNO
&&
2056 !(absolute_reference(ir
,sys
,AbsRef
) || ir
->nsteps
<= 10)) {
2057 warning("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");
2060 if (ir
->coulombtype
== eelCUT
&& ir
->rcoulomb
> 0) {
2062 aloopb
= gmx_mtop_atomloop_block_init(sys
);
2063 while (gmx_mtop_atomloop_block_next(aloopb
,&atom
,&nmol
)) {
2064 if (atom
->q
!= 0 || atom
->qB
!= 0) {
2069 set_warning_line(mdparin
,-1);
2071 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
2072 "You might want to consider using %s electrostatics.\n",
2074 warning_note(err_buf
);
2078 /* Generalized reaction field */
2079 if (ir
->opts
.ngtc
== 0) {
2080 sprintf(err_buf
,"No temperature coupling while using coulombtype %s",
2082 CHECK(ir
->coulombtype
== eelGRF
);
2085 sprintf(err_buf
,"When using coulombtype = %s"
2086 " ref_t for temperature coupling should be > 0",
2088 CHECK((ir
->coulombtype
== eelGRF
) && (ir
->opts
.ref_t
[0] <= 0));
2091 if (ir
->eI
== eiSD1
) {
2093 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
2094 gdt_max
= max(gdt_max
,ir
->delta_t
/ir
->opts
.tau_t
[i
]);
2095 if (0.5*gdt_max
> 0.0015) {
2096 set_warning_line(mdparin
,-1);
2097 sprintf(warn_buf
,"The relative error with integrator %s is 0.5*delta_t/tau_t = %g, you might want to switch to integrator %s\n",
2098 ei_names
[ir
->eI
],0.5*gdt_max
,ei_names
[eiSD2
]);
2104 for(i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++) {
2105 for(m
=0; (m
<DIM
); m
++) {
2106 if (fabs(ir
->opts
.acc
[i
][m
]) > 1e-6) {
2113 snew(mgrp
,sys
->groups
.grps
[egcACC
].nr
);
2114 aloop
= gmx_mtop_atomloop_all_init(sys
);
2115 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
2116 mgrp
[ggrpnr(&sys
->groups
,egcACC
,i
)] += atom
->m
;
2119 for(i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++) {
2120 for(m
=0; (m
<DIM
); m
++)
2121 acc
[m
] += ir
->opts
.acc
[i
][m
]*mgrp
[i
];
2124 for(m
=0; (m
<DIM
); m
++) {
2125 if (fabs(acc
[m
]) > 1e-6) {
2126 const char *dim
[DIM
] = { "X", "Y", "Z" };
2128 "Net Acceleration in %s direction, will %s be corrected\n",
2129 dim
[m
],ir
->nstcomm
!= 0 ? "" : "not");
2130 if (ir
->nstcomm
!= 0 && m
< ndof_com(ir
)) {
2132 for (i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++)
2133 ir
->opts
.acc
[i
][m
] -= acc
[m
];
2140 if (ir
->efep
!= efepNO
&& ir
->sc_alpha
!= 0 &&
2141 !gmx_within_tol(sys
->ffparams
.reppow
,12.0,10*GMX_DOUBLE_EPS
)) {
2142 gmx_fatal(FARGS
,"Soft-core interactions are only supported with VdW repulsion power 12");
2145 if (ir
->ePull
!= epullNO
) {
2146 if (ir
->pull
->grp
[0].nat
== 0) {
2147 absolute_reference(ir
,sys
,AbsRef
);
2148 for(m
=0; m
<DIM
; m
++) {
2149 if (ir
->pull
->dim
[m
] && !AbsRef
[m
]) {
2150 set_warning_line(mdparin
,-1);
2151 warning("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.");
2157 if (ir
->pull
->eGeom
== epullgDIRPBC
) {
2158 for(i
=0; i
<3; i
++) {
2159 for(m
=0; m
<=i
; m
++) {
2160 if ((ir
->epc
!= epcNO
&& ir
->compress
[i
][m
] != 0) ||
2161 ir
->deform
[i
][m
] != 0) {
2162 for(g
=1; g
<ir
->pull
->ngrp
; g
++) {
2163 if (ir
->pull
->grp
[g
].vec
[m
] != 0) {
2164 gmx_fatal(FARGS
,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir
->pull
->eGeom
),'x'+m
);
2176 void double_check(t_inputrec
*ir
,matrix box
,bool bConstr
,int *nerror
)
2182 ptr
= check_box(ir
->ePBC
,box
);
2189 if (bConstr
&& ir
->eConstrAlg
== econtSHAKE
) {
2190 if (ir
->shake_tol
<= 0.0) {
2191 fprintf(stderr
,"ERROR: shake_tol must be > 0 instead of %g\n",
2196 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1) {
2197 sprintf(warn_buf
,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
2198 if (ir
->epc
== epcNO
) {
2201 fprintf(stderr
,"ERROR: %s\n",warn_buf
);
2207 if( (ir
->eConstrAlg
== econtLINCS
) && bConstr
) {
2208 /* If we have Lincs constraints: */
2209 if(ir
->eI
==eiMD
&& ir
->etc
==etcNO
&&
2210 ir
->eConstrAlg
==econtLINCS
&& ir
->nLincsIter
==1) {
2211 sprintf(warn_buf
,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
2215 if ((ir
->eI
== eiCG
|| ir
->eI
== eiLBFGS
) && (ir
->nProjOrder
<8)) {
2216 sprintf(warn_buf
,"For accurate %s with LINCS constraints, lincs_order should be 8 or more.",ei_names
[ir
->eI
]);
2221 if (ir
->LincsWarnAngle
> 90.0) {
2222 sprintf(warn_buf
,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
2224 ir
->LincsWarnAngle
= 90.0;
2227 if (ir
->ePBC
!= epbcNONE
) {
2228 if (ir
->nstlist
== 0) {
2229 warning("With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
2231 bTWIN
= (ir
->rlistlong
> ir
->rlist
);
2232 if (ir
->ns_type
== ensGRID
) {
2233 if (sqr(ir
->rlistlong
) >= max_cutoff2(ir
->ePBC
,box
)) {
2234 fprintf(stderr
,"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",
2235 bTWIN
? (ir
->rcoulomb
==ir
->rlistlong
? "rcoulomb" : "rvdw"):"rlist");
2239 min_size
= min(box
[XX
][XX
],min(box
[YY
][YY
],box
[ZZ
][ZZ
]));
2240 if (2*ir
->rlistlong
>= min_size
) {
2241 fprintf(stderr
,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
2244 fprintf(stderr
,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");