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 anneal
[STRLEN
],anneal_npoints
[STRLEN
],
79 anneal_time
[STRLEN
],anneal_temp
[STRLEN
];
80 static char QMmethod
[STRLEN
],QMbasis
[STRLEN
],QMcharge
[STRLEN
],QMmult
[STRLEN
],
81 bSH
[STRLEN
],CASorbitals
[STRLEN
], CASelectrons
[STRLEN
],SAon
[STRLEN
],
82 SAoff
[STRLEN
],SAsteps
[STRLEN
],bTS
[STRLEN
],bOPT
[STRLEN
];
83 static char efield_x
[STRLEN
],efield_xt
[STRLEN
],efield_y
[STRLEN
],
84 efield_yt
[STRLEN
],efield_z
[STRLEN
],efield_zt
[STRLEN
];
86 enum { egrptpALL
, egrptpALL_GENREST
, egrptpPART
, egrptpONE
};
89 void init_ir(t_inputrec
*ir
, t_gromppopts
*opts
)
91 snew(opts
->include
,STRLEN
);
92 snew(opts
->define
,STRLEN
);
95 static void _low_check(bool b
,char *s
,int *n
)
98 fprintf(stderr
,"\nERROR: %s\n\n",s
);
103 static void check_nst(const char *desc_nst
,int nst
,
104 const char *desc_p
,int *p
)
108 if (*p
> 0 && *p
% nst
!= 0)
110 /* Round up to the next multiple of nst */
111 *p
= ((*p
)/nst
+ 1)*nst
;
112 sprintf(buf
,"%s should be a multiple of %s, changing %s to %d\n",
113 desc_p
,desc_nst
,desc_p
,*p
);
118 void check_ir(const char *mdparin
,t_inputrec
*ir
, t_gromppopts
*opts
,
120 /* Check internal consistency */
122 /* Strange macro: first one fills the err_buf, and then one can check
123 * the condition, which will print the message and increase the error
126 #define CHECK(b) _low_check(b,err_buf,nerror)
131 set_warning_line(mdparin
,-1);
133 /* BASIC CUT-OFF STUFF */
134 if (ir
->rlist
== 0 ||
135 !(( EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) && ir
->rcoulomb
> ir
->rlist
) ||
136 (EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
) && ir
->rvdw
> ir
->rlist
))) {
137 /* No switched potential and/or no twin-range:
138 * we can set the long-range cut-off to the maximum of the other cut-offs.
140 ir
->rlistlong
= max_cutoff(ir
->rlist
,max_cutoff(ir
->rvdw
,ir
->rcoulomb
));
141 } else if (ir
->rlistlong
< 0) {
142 ir
->rlistlong
= max_cutoff(ir
->rlist
,max_cutoff(ir
->rvdw
,ir
->rcoulomb
));
143 sprintf(warn_buf
,"rlistlong was not set, setting it to %g (no buffer)",
147 if (ir
->rlistlong
== 0 && ir
->ePBC
!= epbcNONE
) {
148 warning_error("Can not have an infinite cut-off with PBC");
150 if (ir
->rlistlong
> 0 && (ir
->rlist
== 0 || ir
->rlistlong
< ir
->rlist
)) {
151 warning_error("rlistlong can not be shorter than rlist");
153 if (IR_TWINRANGE(*ir
) && ir
->nstlist
<= 0) {
154 warning_error("Can not have nstlist<=0 with twin-range interactions");
157 /* GENERAL INTEGRATOR STUFF */
158 if (!(ir
->eI
== eiMD
|| EI_VV(ir
->eI
))) {
161 if (!EI_DYNAMICS(ir
->eI
)) {
164 if (EI_DYNAMICS(ir
->eI
)) {
165 if (ir
->nstcalcenergy
< 0) {
166 gmx_fatal(FARGS
,"Can not have nstcalcenergy < 0");
168 if ((ir
->etc
!= etcNO
|| ir
->epc
!= epcNO
) && ir
->nstcalcenergy
== 0) {
169 gmx_fatal(FARGS
,"Can not have nstcalcenergy=0 with global T/P-coupling");
171 if (IR_TWINRANGE(*ir
)) {
172 check_nst("nstlist",ir
->nstlist
,"nstcalcenergy",&ir
->nstcalcenergy
);
174 dt_coupl
= ir
->nstcalcenergy
*ir
->delta_t
;
176 if (ir
->nstcalcenergy
> 1) {
177 /* Energy and log file writing trigger energy calculation,
178 * so we need some checks.
180 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstenergy",&ir
->nstenergy
);
181 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstlog",&ir
->nstlog
);
182 if (ir
->efep
!= efepNO
) {
183 check_nst("nstcalcenergy",ir
->nstcalcenergy
,"nstdhdl",&ir
->nstdhdl
);
189 if ((EI_SD(ir
->eI
) || ir
->eI
== eiBD
) &&
190 ir
->bContinuation
&& ir
->ld_seed
!= -1) {
191 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)");
195 if (EI_TPI(ir
->eI
)) {
196 sprintf(err_buf
,"TPI only works with pbc = %s",epbc_names
[epbcXYZ
]);
197 CHECK(ir
->ePBC
!= epbcXYZ
);
198 sprintf(err_buf
,"TPI only works with ns = %s",ens_names
[ensGRID
]);
199 CHECK(ir
->ns_type
!= ensGRID
);
200 sprintf(err_buf
,"with TPI nstlist should be larger than zero");
201 CHECK(ir
->nstlist
<= 0);
202 sprintf(err_buf
,"TPI does not work with full electrostatics other than PME");
203 CHECK(EEL_FULL(ir
->coulombtype
) && !EEL_PME(ir
->coulombtype
));
207 if ( (opts
->nshake
> 0) && (opts
->bMorse
) ) {
209 "Using morse bond-potentials while constraining bonds is useless");
213 sprintf(err_buf
,"shake_tol must be > 0 instead of %g while using shake",
215 CHECK(((ir
->shake_tol
<= 0.0) && (opts
->nshake
>0) &&
216 (ir
->eConstrAlg
== econtSHAKE
)));
219 sprintf(err_buf
,"walls only work with pbc=%s",epbc_names
[epbcXY
]);
220 CHECK(ir
->nwall
&& ir
->ePBC
!=epbcXY
);
223 if (ir
->ePBC
!= epbcXYZ
&& ir
->nwall
!= 2) {
224 if (ir
->ePBC
== epbcNONE
) {
225 if (ir
->epc
!= epcNO
) {
226 warning("Turning off pressure coupling for vacuum system");
230 sprintf(err_buf
,"Can not have pressure coupling with pbc=%s",
231 epbc_names
[ir
->ePBC
]);
232 CHECK(ir
->epc
!= epcNO
);
234 sprintf(err_buf
,"Can not have Ewald with pbc=%s",epbc_names
[ir
->ePBC
]);
235 CHECK(EEL_FULL(ir
->coulombtype
));
237 sprintf(err_buf
,"Can not have dispersion correction with pbc=%s",
238 epbc_names
[ir
->ePBC
]);
239 CHECK(ir
->eDispCorr
!= edispcNO
);
242 if (ir
->rlist
== 0.0) {
243 sprintf(err_buf
,"can only have neighborlist cut-off zero (=infinite)\n"
244 "with coulombtype = %s or coulombtype = %s\n"
245 "without periodic boundary conditions (pbc = %s) and\n"
246 "rcoulomb and rvdw set to zero",
247 eel_names
[eelCUT
],eel_names
[eelUSER
],epbc_names
[epbcNONE
]);
248 CHECK(((ir
->coulombtype
!= eelCUT
) && (ir
->coulombtype
!= eelUSER
)) ||
249 (ir
->ePBC
!= epbcNONE
) ||
250 (ir
->rcoulomb
!= 0.0) || (ir
->rvdw
!= 0.0));
252 if (ir
->nstlist
< 0) {
253 warning_error("Can not have heuristic neighborlist updates without cut-off");
255 if (ir
->nstlist
> 0) {
256 warning_note("Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
261 if (ir
->nstcomm
== 0) {
262 ir
->comm_mode
= ecmNO
;
264 if (ir
->comm_mode
!= ecmNO
) {
265 if (ir
->nstcomm
< 0) {
266 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");
267 ir
->nstcomm
= abs(ir
->nstcomm
);
270 if (ir
->nstcalcenergy
> 0 && ir
->nstcomm
< ir
->nstcalcenergy
) {
271 warning_note("nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
272 ir
->nstcomm
= ir
->nstcalcenergy
;
275 if (ir
->comm_mode
== ecmANGULAR
) {
276 sprintf(err_buf
,"Can not remove the rotation around the center of mass with periodic molecules");
277 CHECK(ir
->bPeriodicMols
);
278 if (ir
->ePBC
!= epbcNONE
)
279 warning("Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
283 if (EI_STATE_VELOCITY(ir
->eI
) && ir
->ePBC
== epbcNONE
&& ir
->comm_mode
!= ecmANGULAR
) {
284 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.");
287 sprintf(err_buf
,"Free-energy not implemented for Ewald and PPPM");
288 CHECK((ir
->coulombtype
==eelEWALD
|| ir
->coulombtype
==eelPPPM
)
289 && (ir
->efep
!=efepNO
));
291 sprintf(err_buf
,"Twin-range neighbour searching (NS) with simple NS"
292 " algorithm not implemented");
293 CHECK(((ir
->rcoulomb
> ir
->rlist
) || (ir
->rvdw
> ir
->rlist
))
294 && (ir
->ns_type
== ensSIMPLE
));
296 /* TEMPERATURE COUPLING */
297 if(ir
->etc
== etcYES
) {
298 ir
->etc
= etcBERENDSEN
;
299 warning_note("Old option for temperature coupling given: "
300 "changing \"yes\" to \"Berendsen\"\n");
303 if (ir
->etc
== etcBERENDSEN
) {
304 sprintf(warn_buf
,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
305 ETCOUPLTYPE(ir
->etc
),ETCOUPLTYPE(etcVRESCALE
));
309 if((ir
->etc
==etcNOSEHOOVER
|| ir
->etc
==etcANDERSEN
|| ir
->etc
==etcANDERSENINTERVAL
)
310 && ir
->epc
==epcBERENDSEN
) {
311 sprintf(warn_buf
,"Using Berendsen pressure coupling invalidates the "
312 "true ensemble for the thermostat");
316 /* PRESSURE COUPLING */
317 if (ir
->epc
== epcISOTROPIC
) {
318 ir
->epc
= epcBERENDSEN
;
319 warning_note("Old option for pressure coupling given: "
320 "changing \"Isotropic\" to \"Berendsen\"\n");
323 if (ir
->epc
!= epcNO
) {
324 sprintf(err_buf
,"tau_p must be > 0 instead of %g\n",ir
->tau_p
);
325 CHECK(ir
->tau_p
<= 0);
327 if (ir
->tau_p
< 100*dt_coupl
) {
328 sprintf(warn_buf
,"For proper barostat integration tau_p (%g) should be more than two orders of magnitude larger than nstcalcenergy*dt (%g)",
333 sprintf(err_buf
,"compressibility must be > 0 when using pressure"
334 " coupling %s\n",EPCOUPLTYPE(ir
->epc
));
335 CHECK(ir
->compress
[XX
][XX
] < 0 || ir
->compress
[YY
][YY
] < 0 ||
336 ir
->compress
[ZZ
][ZZ
] < 0 ||
337 (trace(ir
->compress
) == 0 && ir
->compress
[YY
][XX
] <= 0 &&
338 ir
->compress
[ZZ
][XX
] <= 0 && ir
->compress
[ZZ
][YY
] <= 0));
340 sprintf(err_buf
,"pressure coupling with PPPM not implemented, use PME");
341 CHECK(ir
->coulombtype
== eelPPPM
);
343 } else if (ir
->coulombtype
== eelPPPM
) {
344 sprintf(warn_buf
,"The pressure with PPPM is incorrect, if you need the pressure use PME");
349 if (ir
->epc
> epcNO
) {
350 if (ir
->epc
!=epcMTTK
) {
351 warning_error("NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
357 /* More checks are in triple check (grompp.c) */
358 if (ir
->coulombtype
== eelSWITCH
) {
359 sprintf(warn_buf
,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
360 eel_names
[ir
->coulombtype
],
361 eel_names
[eelRF_ZERO
]);
365 if (EEL_RF(ir
->coulombtype
) && ir
->epsilon_rf
==1 && ir
->epsilon_r
!=1) {
366 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
);
368 ir
->epsilon_rf
= ir
->epsilon_r
;
372 if (getenv("GALACTIC_DYNAMICS") == NULL
) {
373 sprintf(err_buf
,"epsilon_r must be >= 0 instead of %g\n",ir
->epsilon_r
);
374 CHECK(ir
->epsilon_r
< 0);
377 if (EEL_RF(ir
->coulombtype
)) {
378 /* reaction field (at the cut-off) */
380 if (ir
->coulombtype
== eelRF_ZERO
) {
381 sprintf(err_buf
,"With coulombtype = %s, epsilon_rf must be 0",
382 eel_names
[ir
->coulombtype
]);
383 CHECK(ir
->epsilon_rf
!= 0);
386 sprintf(err_buf
,"epsilon_rf must be >= epsilon_r");
387 CHECK((ir
->epsilon_rf
< ir
->epsilon_r
&& ir
->epsilon_rf
!= 0) ||
388 (ir
->epsilon_r
== 0));
389 if (ir
->epsilon_rf
== ir
->epsilon_r
) {
390 sprintf(warn_buf
,"Using epsilon_rf = epsilon_r with %s does not make sense",
391 eel_names
[ir
->coulombtype
]);
395 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
396 * means the interaction is zero outside rcoulomb, but it helps to
397 * provide accurate energy conservation.
399 if (EEL_ZERO_AT_CUTOFF(ir
->coulombtype
)) {
400 if (EEL_SWITCHED(ir
->coulombtype
)) {
402 "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
403 eel_names
[ir
->coulombtype
]);
404 CHECK(ir
->rcoulomb_switch
>= ir
->rcoulomb
);
406 } else if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
)) {
407 sprintf(err_buf
,"With coulombtype = %s, rcoulomb must be >= rlist",
408 eel_names
[ir
->coulombtype
]);
409 CHECK(ir
->rlist
> ir
->rcoulomb
);
412 if (EEL_FULL(ir
->coulombtype
)) {
413 if (ir
->coulombtype
==eelPMESWITCH
|| ir
->coulombtype
==eelPMEUSER
) {
414 sprintf(err_buf
,"With coulombtype = %s, rcoulomb must be <= rlist",
415 eel_names
[ir
->coulombtype
]);
416 CHECK(ir
->rcoulomb
> ir
->rlist
);
418 if (ir
->coulombtype
== eelPME
) {
420 "With coulombtype = %s, rcoulomb must be equal to rlist\n"
421 "If you want optimal energy conservation or exact integration use %s",
422 eel_names
[ir
->coulombtype
],eel_names
[eelPMESWITCH
]);
425 "With coulombtype = %s, rcoulomb must be equal to rlist",
426 eel_names
[ir
->coulombtype
]);
428 CHECK(ir
->rcoulomb
!= ir
->rlist
);
432 if (EEL_PME(ir
->coulombtype
)) {
433 if (ir
->pme_order
< 3) {
434 warning_error("pme_order can not be smaller than 3");
438 if (ir
->nwall
==2 && EEL_FULL(ir
->coulombtype
)) {
439 if (ir
->ewald_geometry
== eewg3D
) {
440 sprintf(warn_buf
,"With pbc=%s you should use ewald_geometry=%s",
441 epbc_names
[ir
->ePBC
],eewg_names
[eewg3DC
]);
444 /* This check avoids extra pbc coding for exclusion corrections */
445 sprintf(err_buf
,"wall_ewald_zfac should be >= 2");
446 CHECK(ir
->wall_ewald_zfac
< 2);
449 if (EVDW_SWITCHED(ir
->vdwtype
)) {
450 sprintf(err_buf
,"With vdwtype = %s rvdw_switch must be < rvdw",
451 evdw_names
[ir
->vdwtype
]);
452 CHECK(ir
->rvdw_switch
>= ir
->rvdw
);
453 } else if (ir
->vdwtype
== evdwCUT
) {
454 sprintf(err_buf
,"With vdwtype = %s, rvdw must be >= rlist",evdw_names
[ir
->vdwtype
]);
455 CHECK(ir
->rlist
> ir
->rvdw
);
457 if ((EEL_SWITCHED(ir
->coulombtype
) || ir
->coulombtype
== eelRF_ZERO
)
458 && (ir
->rlistlong
<= ir
->rcoulomb
)) {
459 sprintf(warn_buf
,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
460 IR_TWINRANGE(*ir
) ? "rlistlong" : "rlist");
463 if (EVDW_SWITCHED(ir
->vdwtype
) && (ir
->rlistlong
<= ir
->rvdw
)) {
464 sprintf(warn_buf
,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
465 IR_TWINRANGE(*ir
) ? "rlistlong" : "rlist");
469 if (ir
->vdwtype
== evdwUSER
&& ir
->eDispCorr
!= edispcNO
) {
470 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.");
473 if (ir
->nstlist
== -1) {
475 "nstlist=-1 only works with switched or shifted potentials,\n"
476 "suggestion: use vdw-type=%s and coulomb-type=%s",
477 evdw_names
[evdwSHIFT
],eel_names
[eelPMESWITCH
]);
478 CHECK(!(EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) &&
479 EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
)));
481 sprintf(err_buf
,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
482 CHECK(ir
->rvdw
>= ir
->rlist
|| ir
->rcoulomb
>= ir
->rlist
);
484 sprintf(err_buf
,"nstlist can not be smaller than -1");
485 CHECK(ir
->nstlist
< -1);
487 if (ir
->eI
== eiLBFGS
&& (ir
->coulombtype
==eelCUT
|| ir
->vdwtype
==evdwCUT
)
489 warning("For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
492 if (ir
->eI
== eiLBFGS
&& ir
->nbfgscorr
<= 0) {
493 warning("Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
497 if (ir
->efep
!= efepNO
) {
498 sprintf(err_buf
,"The soft-core power is %d and can only be 1 or 2",
500 CHECK(ir
->sc_alpha
!=0 && ir
->sc_power
!=1 && ir
->sc_power
!=2);
503 /* ENERGY CONSERVATION */
504 if (ir
->eI
== eiMD
&& ir
->etc
== etcNO
) {
505 if (!EVDW_ZERO_AT_CUTOFF(ir
->vdwtype
) && ir
->rvdw
> 0) {
506 sprintf(warn_buf
,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
507 evdw_names
[evdwSHIFT
]);
510 if (!EEL_ZERO_AT_CUTOFF(ir
->coulombtype
) && ir
->rcoulomb
> 0) {
511 sprintf(warn_buf
,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
512 eel_names
[eelPMESWITCH
],eel_names
[eelRF_ZERO
]);
517 if(ir
->coulombtype
==eelGB_NOTUSED
)
519 ir
->coulombtype
=eelCUT
;
520 ir
->implicit_solvent
=eisGBSA
;
521 fprintf(stderr
,"Note: Old option for generalized born electrostatics given:\n"
522 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
523 "setting implicit_solvent value to \"GBSA\" in input section.\n");
526 if(ir
->implicit_solvent
==eisGBSA
)
528 sprintf(err_buf
,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
529 CHECK(ir
->rgbradii
!= ir
->rlist
);
533 static int str_nelem(const char *str
,int maxptr
,char *ptr
[])
541 while (*copy
!= '\0') {
543 gmx_fatal(FARGS
,"Too many groups on line: '%s' (max is %d)",
548 while ((*copy
!= '\0') && !isspace(*copy
))
562 static void parse_n_double(char *str
,int *n
,double **r
)
567 *n
= str_nelem(str
,MAXPTR
,ptr
);
570 for(i
=0; i
<*n
; i
++) {
571 (*r
)[i
] = strtod(ptr
[i
],NULL
);
575 static void do_wall_params(t_inputrec
*ir
,
576 char *wall_atomtype
, char *wall_density
,
583 opts
->wall_atomtype
[0] = NULL
;
584 opts
->wall_atomtype
[1] = NULL
;
586 ir
->wall_atomtype
[0] = -1;
587 ir
->wall_atomtype
[1] = -1;
588 ir
->wall_density
[0] = 0;
589 ir
->wall_density
[1] = 0;
592 nstr
= str_nelem(wall_atomtype
,MAXPTR
,names
);
593 if (nstr
!= ir
->nwall
)
594 gmx_fatal(FARGS
,"Expected %d elements for wall_atomtype, found %d",
596 for(i
=0; i
<ir
->nwall
; i
++)
597 opts
->wall_atomtype
[i
] = strdup(names
[i
]);
599 if (ir
->wall_type
!= ewtTABLE
) {
600 nstr
= str_nelem(wall_density
,MAXPTR
,names
);
601 if (nstr
!= ir
->nwall
)
602 gmx_fatal(FARGS
,"Expected %d elements for wall_density, found %d",
604 for(i
=0; i
<ir
->nwall
; i
++) {
605 sscanf(names
[i
],"%lf",&dbl
);
607 gmx_fatal(FARGS
,"wall_density[%d] = %f\n",i
,dbl
);
608 ir
->wall_density
[i
] = dbl
;
614 static void add_wall_energrps(gmx_groups_t
*groups
,int nwall
,t_symtab
*symtab
)
621 srenew(groups
->grpname
,groups
->ngrpname
+nwall
);
622 grps
= &(groups
->grps
[egcENER
]);
623 srenew(grps
->nm_ind
,grps
->nr
+nwall
);
624 for(i
=0; i
<nwall
; i
++) {
625 sprintf(str
,"wall%d",i
);
626 groups
->grpname
[groups
->ngrpname
] = put_symtab(symtab
,str
);
627 grps
->nm_ind
[grps
->nr
++] = groups
->ngrpname
++;
632 void get_ir(const char *mdparin
,const char *mdparout
,
633 t_inputrec
*ir
,t_gromppopts
*opts
, int *nerror
)
641 inp
=read_inpfile(mdparin
,&ninp
, NULL
);
643 snew(dumstr
[0],STRLEN
);
644 snew(dumstr
[1],STRLEN
);
648 REM_TYPE("domain-decomposition");
649 REPL_TYPE("unconstrained-start","continuation");
650 REM_TYPE("dihre-tau");
651 REM_TYPE("nstdihreout");
652 REM_TYPE("nstcheckpoint");
654 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
655 CTYPE ("Preprocessor information: use cpp syntax.");
656 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
657 STYPE ("include", opts
->include
, NULL
);
658 CTYPE ("e.g.: -DI_Want_Cookies -DMe_Too");
659 STYPE ("define", opts
->define
, NULL
);
661 CCTYPE ("RUN CONTROL PARAMETERS");
662 EETYPE("integrator", ir
->eI
, ei_names
, nerror
, TRUE
);
663 CTYPE ("Start time and timestep in ps");
664 RTYPE ("tinit", ir
->init_t
, 0.0);
665 RTYPE ("dt", ir
->delta_t
, 0.001);
666 STEPTYPE ("nsteps", ir
->nsteps
, 0);
667 CTYPE ("For exact run continuation or redoing part of a run");
668 STEPTYPE ("init_step",ir
->init_step
, 0);
669 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
670 ITYPE ("simulation_part", ir
->simulation_part
, 1);
671 CTYPE ("mode for center of mass motion removal");
672 CTYPE ("energy calculation and T/P-coupling frequency");
673 ITYPE ("nstcalcenergy",ir
->nstcalcenergy
, 1);
674 EETYPE("comm-mode", ir
->comm_mode
, ecm_names
, nerror
, TRUE
);
675 CTYPE ("number of steps for center of mass motion removal");
676 ITYPE ("nstcomm", ir
->nstcomm
, 1);
677 CTYPE ("group(s) for center of mass motion removal");
678 STYPE ("comm-grps", vcm
, NULL
);
680 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
681 CTYPE ("Friction coefficient (amu/ps) and random seed");
682 RTYPE ("bd-fric", ir
->bd_fric
, 0.0);
683 ITYPE ("ld-seed", ir
->ld_seed
, 1993);
686 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
687 CTYPE ("Force tolerance and initial step-size");
688 RTYPE ("emtol", ir
->em_tol
, 10.0);
689 RTYPE ("emstep", ir
->em_stepsize
,0.01);
690 CTYPE ("Max number of iterations in relax_shells");
691 ITYPE ("niter", ir
->niter
, 20);
692 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
693 RTYPE ("fcstep", ir
->fc_stepsize
, 0);
694 CTYPE ("Frequency of steepest descents steps when doing CG");
695 ITYPE ("nstcgsteep", ir
->nstcgsteep
, 1000);
696 ITYPE ("nbfgscorr", ir
->nbfgscorr
, 10);
698 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
699 RTYPE ("rtpi", ir
->rtpi
, 0.05);
702 CCTYPE ("OUTPUT CONTROL OPTIONS");
703 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
704 ITYPE ("nstxout", ir
->nstxout
, 100);
705 ITYPE ("nstvout", ir
->nstvout
, 100);
706 ITYPE ("nstfout", ir
->nstfout
, 0);
707 ir
->nstcheckpoint
= 1000;
708 CTYPE ("Output frequency for energies to log file and energy file");
709 ITYPE ("nstlog", ir
->nstlog
, 100);
710 ITYPE ("nstenergy", ir
->nstenergy
, 100);
711 CTYPE ("Output frequency and precision for xtc file");
712 ITYPE ("nstxtcout", ir
->nstxtcout
, 0);
713 RTYPE ("xtc-precision",ir
->xtcprec
, 1000.0);
714 CTYPE ("This selects the subset of atoms for the xtc file. You can");
715 CTYPE ("select multiple groups. By default all atoms will be written.");
716 STYPE ("xtc-grps", xtc_grps
, NULL
);
717 CTYPE ("Selection of energy groups");
718 STYPE ("energygrps", energy
, NULL
);
720 /* Neighbor searching */
721 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
722 CTYPE ("nblist update frequency");
723 ITYPE ("nstlist", ir
->nstlist
, 10);
724 CTYPE ("ns algorithm (simple or grid)");
725 EETYPE("ns-type", ir
->ns_type
, ens_names
, nerror
, TRUE
);
726 /* set ndelta to the optimal value of 2 */
728 CTYPE ("Periodic boundary conditions: xyz, no, xy");
729 EETYPE("pbc", ir
->ePBC
, epbc_names
, nerror
, TRUE
);
730 EETYPE("periodic_molecules", ir
->bPeriodicMols
, yesno_names
, nerror
, TRUE
);
731 CTYPE ("nblist cut-off");
732 RTYPE ("rlist", ir
->rlist
, 1.0);
733 CTYPE ("long-range cut-off for switched potentials");
734 RTYPE ("rlistlong", ir
->rlistlong
, -1);
737 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
738 CTYPE ("Method for doing electrostatics");
739 EETYPE("coulombtype", ir
->coulombtype
, eel_names
, nerror
, TRUE
);
740 CTYPE ("cut-off lengths");
741 RTYPE ("rcoulomb-switch", ir
->rcoulomb_switch
, 0.0);
742 RTYPE ("rcoulomb", ir
->rcoulomb
, 1.0);
743 CTYPE ("Relative dielectric constant for the medium and the reaction field");
744 RTYPE ("epsilon_r", ir
->epsilon_r
, 1.0);
745 RTYPE ("epsilon_rf", ir
->epsilon_rf
, 1.0);
746 CTYPE ("Method for doing Van der Waals");
747 EETYPE("vdw-type", ir
->vdwtype
, evdw_names
, nerror
, TRUE
);
748 CTYPE ("cut-off lengths");
749 RTYPE ("rvdw-switch", ir
->rvdw_switch
, 0.0);
750 RTYPE ("rvdw", ir
->rvdw
, 1.0);
751 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
752 EETYPE("DispCorr", ir
->eDispCorr
, edispc_names
, nerror
, TRUE
);
753 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
754 RTYPE ("table-extension", ir
->tabext
, 1.0);
755 CTYPE ("Seperate tables between energy group pairs");
756 STYPE ("energygrp_table", egptable
, NULL
);
757 CTYPE ("Spacing for the PME/PPPM FFT grid");
758 RTYPE ("fourierspacing", opts
->fourierspacing
,0.12);
759 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
760 ITYPE ("fourier_nx", ir
->nkx
, 0);
761 ITYPE ("fourier_ny", ir
->nky
, 0);
762 ITYPE ("fourier_nz", ir
->nkz
, 0);
763 CTYPE ("EWALD/PME/PPPM parameters");
764 ITYPE ("pme_order", ir
->pme_order
, 4);
765 RTYPE ("ewald_rtol", ir
->ewald_rtol
, 0.00001);
766 EETYPE("ewald_geometry", ir
->ewald_geometry
, eewg_names
, nerror
, TRUE
);
767 RTYPE ("epsilon_surface", ir
->epsilon_surface
, 0.0);
768 EETYPE("optimize_fft",ir
->bOptFFT
, yesno_names
, nerror
, TRUE
);
770 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
771 EETYPE("implicit_solvent", ir
->implicit_solvent
, eis_names
, nerror
, TRUE
);
773 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
774 CTYPE ("Algorithm for calculating Born radii");
775 EETYPE("gb_algorithm", ir
->gb_algorithm
, egb_names
, nerror
, TRUE
);
776 CTYPE ("Frequency of calculating the Born radii inside rlist");
777 ITYPE ("nstgbradii", ir
->nstgbradii
, 1);
778 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
779 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
780 RTYPE ("rgbradii", ir
->rgbradii
, 1.0);
781 CTYPE ("Dielectric coefficient of the implicit solvent");
782 RTYPE ("gb_epsilon_solvent",ir
->gb_epsilon_solvent
, 80.0);
783 CTYPE ("Salt concentration in M for Generalized Born models");
784 RTYPE ("gb_saltconc", ir
->gb_saltconc
, 0.0);
785 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
786 RTYPE ("gb_obc_alpha", ir
->gb_obc_alpha
, 1.0);
787 RTYPE ("gb_obc_beta", ir
->gb_obc_beta
, 0.8);
788 RTYPE ("gb_obc_gamma", ir
->gb_obc_gamma
, 4.85);
789 RTYPE ("gb_dielectric_offset", ir
->gb_dielectric_offset
, 0.009);
790 EETYPE("sa_algorithm", ir
->sa_algorithm
, esa_names
, nerror
, TRUE
);
791 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
792 CTYPE ("The default value (2.092) corresponds to 0.005 kcal/mol/Angstrom^2.");
793 RTYPE ("sa_surface_tension", ir
->sa_surface_tension
, 2.092);
796 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
797 CTYPE ("Temperature coupling");
798 EETYPE("tcoupl", ir
->etc
, etcoupl_names
, nerror
, TRUE
);
799 CTYPE ("Groups to couple separately");
800 STYPE ("tc-grps", tcgrps
, NULL
);
801 CTYPE ("Time constant (ps) and reference temperature (K)");
802 STYPE ("tau-t", tau_t
, NULL
);
803 STYPE ("ref-t", ref_t
, NULL
);
804 CTYPE ("Pressure coupling");
805 EETYPE("Pcoupl", ir
->epc
, epcoupl_names
, nerror
, TRUE
);
806 EETYPE("Pcoupltype", ir
->epct
, epcoupltype_names
, nerror
, TRUE
);
807 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
808 RTYPE ("tau-p", ir
->tau_p
, 1.0);
809 STYPE ("compressibility", dumstr
[0], NULL
);
810 STYPE ("ref-p", dumstr
[1], NULL
);
811 CTYPE ("Scaling of reference coordinates, No, All or COM");
812 EETYPE ("refcoord_scaling",ir
->refcoord_scaling
,erefscaling_names
, nerror
, TRUE
);
814 CTYPE ("Random seed for Andersen thermostat");
815 ITYPE ("andersen_seed", ir
->andersen_seed
, 815131);
818 CCTYPE ("OPTIONS FOR QMMM calculations");
819 EETYPE("QMMM", ir
->bQMMM
, yesno_names
, nerror
, TRUE
);
820 CTYPE ("Groups treated Quantum Mechanically");
821 STYPE ("QMMM-grps", QMMM
, NULL
);
823 STYPE("QMmethod", QMmethod
, NULL
);
824 CTYPE ("QMMM scheme");
825 EETYPE("QMMMscheme", ir
->QMMMscheme
, eQMMMscheme_names
, nerror
, TRUE
);
826 CTYPE ("QM basisset");
827 STYPE("QMbasis", QMbasis
, NULL
);
829 STYPE ("QMcharge", QMcharge
,NULL
);
830 CTYPE ("QM multiplicity");
831 STYPE ("QMmult", QMmult
,NULL
);
832 CTYPE ("Surface Hopping");
833 STYPE ("SH", bSH
, NULL
);
834 CTYPE ("CAS space options");
835 STYPE ("CASorbitals", CASorbitals
, NULL
);
836 STYPE ("CASelectrons", CASelectrons
, NULL
);
837 STYPE ("SAon", SAon
, NULL
);
838 STYPE ("SAoff",SAoff
,NULL
);
839 STYPE ("SAsteps", SAsteps
, NULL
);
840 CTYPE ("Scale factor for MM charges");
841 RTYPE ("MMChargeScaleFactor", ir
->scalefactor
, 1.0);
842 CTYPE ("Optimization of QM subsystem");
843 STYPE ("bOPT", bOPT
, NULL
);
844 STYPE ("bTS", bTS
, NULL
);
846 /* Simulated annealing */
847 CCTYPE("SIMULATED ANNEALING");
848 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
849 STYPE ("annealing", anneal
, NULL
);
850 CTYPE ("Number of time points to use for specifying annealing in each group");
851 STYPE ("annealing_npoints", anneal_npoints
, NULL
);
852 CTYPE ("List of times at the annealing points for each group");
853 STYPE ("annealing_time", anneal_time
, NULL
);
854 CTYPE ("Temp. at each annealing point, for each group.");
855 STYPE ("annealing_temp", anneal_temp
, NULL
);
858 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
859 EETYPE("gen-vel", opts
->bGenVel
, yesno_names
, nerror
, TRUE
);
860 RTYPE ("gen-temp", opts
->tempi
, 300.0);
861 ITYPE ("gen-seed", opts
->seed
, 173529);
864 CCTYPE ("OPTIONS FOR BONDS");
865 EETYPE("constraints", opts
->nshake
, constraints
, nerror
, TRUE
);
866 CTYPE ("Type of constraint algorithm");
867 EETYPE("constraint-algorithm", ir
->eConstrAlg
, econstr_names
, nerror
, TRUE
);
868 CTYPE ("Do not constrain the start configuration");
869 EETYPE("continuation", ir
->bContinuation
, yesno_names
, nerror
, TRUE
);
870 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
871 EETYPE("Shake-SOR", ir
->bShakeSOR
, yesno_names
, nerror
, TRUE
);
872 CTYPE ("Relative tolerance of shake");
873 RTYPE ("shake-tol", ir
->shake_tol
, 0.0001);
874 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
875 ITYPE ("lincs-order", ir
->nProjOrder
, 4);
876 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
877 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
878 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
879 ITYPE ("lincs-iter", ir
->nLincsIter
, 1);
880 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
881 CTYPE ("rotates over more degrees than");
882 RTYPE ("lincs-warnangle", ir
->LincsWarnAngle
, 30.0);
883 CTYPE ("Convert harmonic bonds to morse potentials");
884 EETYPE("morse", opts
->bMorse
,yesno_names
, nerror
, TRUE
);
886 /* Energy group exclusions */
887 CCTYPE ("ENERGY GROUP EXCLUSIONS");
888 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
889 STYPE ("energygrp_excl", egpexcl
, NULL
);
893 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
894 ITYPE ("nwall", ir
->nwall
, 0);
895 EETYPE("wall_type", ir
->wall_type
, ewt_names
, nerror
, TRUE
);
896 RTYPE ("wall_r_linpot", ir
->wall_r_linpot
, -1);
897 STYPE ("wall_atomtype", wall_atomtype
, NULL
);
898 STYPE ("wall_density", wall_density
, NULL
);
899 RTYPE ("wall_ewald_zfac", ir
->wall_ewald_zfac
, 3);
902 CCTYPE("COM PULLING");
903 CTYPE("Pull type: no, umbrella, constraint or constant_force");
904 EETYPE("pull", ir
->ePull
, epull_names
, nerror
, TRUE
);
905 if (ir
->ePull
!= epullNO
) {
907 pull_grp
= read_pullparams(&ninp
,&inp
,ir
->pull
,&opts
->pull_start
,nerror
);
911 CCTYPE("NMR refinement stuff");
912 CTYPE ("Distance restraints type: No, Simple or Ensemble");
913 EETYPE("disre", ir
->eDisre
, edisre_names
, nerror
, TRUE
);
914 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
915 EETYPE("disre-weighting", ir
->eDisreWeighting
, edisreweighting_names
, nerror
, TRUE
);
916 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
917 EETYPE("disre-mixed", ir
->bDisreMixed
, yesno_names
, nerror
, TRUE
);
918 RTYPE ("disre-fc", ir
->dr_fc
, 1000.0);
919 RTYPE ("disre-tau", ir
->dr_tau
, 0.0);
920 CTYPE ("Output frequency for pair distances to energy file");
921 ITYPE ("nstdisreout", ir
->nstdisreout
, 100);
922 CTYPE ("Orientation restraints: No or Yes");
923 EETYPE("orire", opts
->bOrire
, yesno_names
, nerror
, TRUE
);
924 CTYPE ("Orientation restraints force constant and tau for time averaging");
925 RTYPE ("orire-fc", ir
->orires_fc
, 0.0);
926 RTYPE ("orire-tau", ir
->orires_tau
, 0.0);
927 STYPE ("orire-fitgrp",orirefitgrp
, NULL
);
928 CTYPE ("Output frequency for trace(SD) and S to energy file");
929 ITYPE ("nstorireout", ir
->nstorireout
, 100);
930 CTYPE ("Dihedral angle restraints: No or Yes");
931 EETYPE("dihre", opts
->bDihre
, yesno_names
, nerror
, TRUE
);
932 RTYPE ("dihre-fc", ir
->dihre_fc
, 1000.0);
934 /* Free energy stuff */
935 CCTYPE ("Free energy control stuff");
936 EETYPE("free-energy", ir
->efep
, efep_names
, nerror
, TRUE
);
937 RTYPE ("init-lambda", ir
->init_lambda
,0.0);
938 RTYPE ("delta-lambda",ir
->delta_lambda
,0.0);
939 STYPE ("foreign_lambda", foreign_lambda
, NULL
);
940 RTYPE ("sc-alpha",ir
->sc_alpha
,0.0);
941 ITYPE ("sc-power",ir
->sc_power
,0);
942 RTYPE ("sc-sigma",ir
->sc_sigma
,0.3);
943 ITYPE ("nstdhdl", ir
->nstdhdl
, 10);
944 STYPE ("couple-moltype", couple_moltype
, NULL
);
945 EETYPE("couple-lambda0", opts
->couple_lam0
, couple_lam
, nerror
, TRUE
);
946 EETYPE("couple-lambda1", opts
->couple_lam1
, couple_lam
, nerror
, TRUE
);
947 EETYPE("couple-intramol", opts
->bCoupleIntra
, yesno_names
, nerror
, TRUE
);
949 /* Non-equilibrium MD stuff */
950 CCTYPE("Non-equilibrium MD stuff");
951 STYPE ("acc-grps", accgrps
, NULL
);
952 STYPE ("accelerate", acc
, NULL
);
953 STYPE ("freezegrps", freeze
, NULL
);
954 STYPE ("freezedim", frdim
, NULL
);
955 RTYPE ("cos-acceleration", ir
->cos_accel
, 0);
956 STYPE ("deform", deform
, NULL
);
958 /* Electric fields */
959 CCTYPE("Electric fields");
960 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
961 CTYPE ("and a phase angle (real)");
962 STYPE ("E-x", efield_x
, NULL
);
963 STYPE ("E-xt", efield_xt
, NULL
);
964 STYPE ("E-y", efield_y
, NULL
);
965 STYPE ("E-yt", efield_yt
, NULL
);
966 STYPE ("E-z", efield_z
, NULL
);
967 STYPE ("E-zt", efield_zt
, NULL
);
969 /* User defined thingies */
970 CCTYPE ("User defined thingies");
971 STYPE ("user1-grps", user1
, NULL
);
972 STYPE ("user2-grps", user2
, NULL
);
973 ITYPE ("userint1", ir
->userint1
, 0);
974 ITYPE ("userint2", ir
->userint2
, 0);
975 ITYPE ("userint3", ir
->userint3
, 0);
976 ITYPE ("userint4", ir
->userint4
, 0);
977 RTYPE ("userreal1", ir
->userreal1
, 0);
978 RTYPE ("userreal2", ir
->userreal2
, 0);
979 RTYPE ("userreal3", ir
->userreal3
, 0);
980 RTYPE ("userreal4", ir
->userreal4
, 0);
983 write_inpfile(mdparout
,ninp
,inp
,FALSE
);
984 for (i
=0; (i
<ninp
); i
++) {
990 /* Process options if necessary */
992 for(i
=0; i
<2*DIM
; i
++)
997 if (sscanf(dumstr
[m
],"%lf",&(dumdub
[m
][XX
]))!=1) {
999 "ERROR: pressure coupling not enough values (I need 1)\n");
1002 dumdub
[m
][YY
]=dumdub
[m
][ZZ
]=dumdub
[m
][XX
];
1004 case epctSEMIISOTROPIC
:
1005 case epctSURFACETENSION
:
1006 if (sscanf(dumstr
[m
],"%lf%lf",
1007 &(dumdub
[m
][XX
]),&(dumdub
[m
][ZZ
]))!=2) {
1009 "ERROR: pressure coupling not enough values (I need 2)\n");
1012 dumdub
[m
][YY
]=dumdub
[m
][XX
];
1014 case epctANISOTROPIC
:
1015 if (sscanf(dumstr
[m
],"%lf%lf%lf%lf%lf%lf",
1016 &(dumdub
[m
][XX
]),&(dumdub
[m
][YY
]),&(dumdub
[m
][ZZ
]),
1017 &(dumdub
[m
][3]),&(dumdub
[m
][4]),&(dumdub
[m
][5]))!=6) {
1019 "ERROR: pressure coupling not enough values (I need 6)\n");
1024 gmx_fatal(FARGS
,"Pressure coupling type %s not implemented yet",
1025 epcoupltype_names
[ir
->epct
]);
1029 clear_mat(ir
->ref_p
);
1030 clear_mat(ir
->compress
);
1031 for(i
=0; i
<DIM
; i
++) {
1032 ir
->ref_p
[i
][i
] = dumdub
[1][i
];
1033 ir
->compress
[i
][i
] = dumdub
[0][i
];
1035 if (ir
->epct
== epctANISOTROPIC
) {
1036 ir
->ref_p
[XX
][YY
] = dumdub
[1][3];
1037 ir
->ref_p
[XX
][ZZ
] = dumdub
[1][4];
1038 ir
->ref_p
[YY
][ZZ
] = dumdub
[1][5];
1039 if (ir
->ref_p
[XX
][YY
]!=0 && ir
->ref_p
[XX
][ZZ
]!=0 && ir
->ref_p
[YY
][ZZ
]!=0) {
1040 warning("All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1042 ir
->compress
[XX
][YY
] = dumdub
[0][3];
1043 ir
->compress
[XX
][ZZ
] = dumdub
[0][4];
1044 ir
->compress
[YY
][ZZ
] = dumdub
[0][5];
1045 for(i
=0; i
<DIM
; i
++) {
1046 for(m
=0; m
<i
; m
++) {
1047 ir
->ref_p
[i
][m
] = ir
->ref_p
[m
][i
];
1048 ir
->compress
[i
][m
] = ir
->compress
[m
][i
];
1053 if (ir
->comm_mode
== ecmNO
)
1056 opts
->couple_moltype
= NULL
;
1057 if (strlen(couple_moltype
) > 0) {
1058 if (ir
->efep
!= efepNO
) {
1059 opts
->couple_moltype
= strdup(couple_moltype
);
1060 if (opts
->couple_lam0
== opts
->couple_lam1
)
1061 warning("The lambda=0 and lambda=1 states for coupling are identical");
1062 if (ir
->eI
== eiMD
&& (opts
->couple_lam0
== ecouplamNONE
||
1063 opts
->couple_lam1
== ecouplamNONE
)) {
1064 warning("For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
1067 warning("Can not couple a molecule with free_energy = no");
1071 do_wall_params(ir
,wall_atomtype
,wall_density
,opts
);
1073 if (opts
->bOrire
&& str_nelem(orirefitgrp
,MAXPTR
,NULL
)!=1) {
1074 fprintf(stderr
,"ERROR: Need one orientation restraint fit group\n");
1078 clear_mat(ir
->deform
);
1081 m
= sscanf(deform
,"%lf %lf %lf %lf %lf %lf",
1082 &(dumdub
[0][0]),&(dumdub
[0][1]),&(dumdub
[0][2]),
1083 &(dumdub
[0][3]),&(dumdub
[0][4]),&(dumdub
[0][5]));
1085 ir
->deform
[i
][i
] = dumdub
[0][i
];
1086 ir
->deform
[YY
][XX
] = dumdub
[0][3];
1087 ir
->deform
[ZZ
][XX
] = dumdub
[0][4];
1088 ir
->deform
[ZZ
][YY
] = dumdub
[0][5];
1089 if (ir
->epc
!= epcNO
) {
1092 if (ir
->deform
[i
][j
]!=0 && ir
->compress
[i
][j
]!=0) {
1093 fprintf(stderr
,"ERROR: A box element has deform set and compressibility > 0\n");
1098 if (ir
->deform
[i
][j
]!=0) {
1099 for(m
=j
; m
<DIM
; m
++)
1100 if (ir
->compress
[m
][j
]!=0) {
1101 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.");
1107 if (ir
->efep
!= efepNO
) {
1108 parse_n_double(foreign_lambda
,&ir
->n_flambda
,&ir
->flambda
);
1109 if (ir
->n_flambda
> 0 && ir
->rlist
< max(ir
->rvdw
,ir
->rcoulomb
)) {
1110 warning_note("For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
1120 static int search_QMstring(char *s
,int ng
,const char *gn
[])
1122 /* same as normal search_string, but this one searches QM strings */
1125 for(i
=0; (i
<ng
); i
++)
1126 if (strcasecmp(s
,gn
[i
]) == 0)
1129 gmx_fatal(FARGS
,"this QM method or basisset (%s) is not implemented\n!",s
);
1133 } /* search_QMstring */
1136 int search_string(char *s
,int ng
,char *gn
[])
1140 for(i
=0; (i
<ng
); i
++)
1141 if (strcasecmp(s
,gn
[i
]) == 0)
1144 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
);
1149 static bool do_numbering(int natoms
,gmx_groups_t
*groups
,int ng
,char *ptrs
[],
1150 t_blocka
*block
,char *gnames
[],
1151 int gtype
,int restnm
,
1152 int grptp
,bool bVerbose
)
1154 unsigned short *cbuf
;
1155 t_grps
*grps
=&(groups
->grps
[gtype
]);
1156 int i
,j
,gid
,aj
,ognr
,ntot
=0;
1161 fprintf(debug
,"Starting numbering %d groups of type %d\n",ng
,gtype
);
1163 title
= gtypes
[gtype
];
1166 for(i
=0; (i
<natoms
); i
++)
1169 snew(grps
->nm_ind
,ng
+1); /* +1 for possible rest group */
1170 for(i
=0; (i
<ng
); i
++) {
1171 /* Lookup the group name in the block structure */
1172 gid
= search_string(ptrs
[i
],block
->nr
,gnames
);
1173 if ((grptp
!= egrptpONE
) || (i
== 0))
1174 grps
->nm_ind
[grps
->nr
++]=gid
;
1176 fprintf(debug
,"Found gid %d for group %s\n",gid
,ptrs
[i
]);
1178 /* Now go over the atoms in the group */
1179 for(j
=block
->index
[gid
]; (j
<block
->index
[gid
+1]); j
++) {
1182 /* Range checking */
1183 if ((aj
< 0) || (aj
>= natoms
))
1184 gmx_fatal(FARGS
,"Invalid atom number %d in indexfile",aj
);
1186 /* Lookup up the old group number */
1189 gmx_fatal(FARGS
,"Atom %d in multiple %s groups (%d and %d)",
1190 aj
+1,title
,ognr
+1,i
+1);
1192 /* Store the group number in buffer */
1193 if (grptp
== egrptpONE
)
1202 /* Now check whether we have done all atoms */
1204 if (ntot
!= natoms
) {
1205 if (grptp
== egrptpALL
) {
1206 gmx_fatal(FARGS
,"%d atoms are not part of any of the %s groups",
1208 } else if (grptp
== egrptpPART
) {
1209 sprintf(warn_buf
,"%d atoms are not part of any of the %s groups",
1213 /* Assign all atoms currently unassigned to a rest group */
1214 for(j
=0; (j
<natoms
); j
++) {
1215 if (cbuf
[j
] == NOGID
) {
1220 if (grptp
!= egrptpPART
) {
1223 "Making dummy/rest group for %s containing %d elements\n",
1225 /* Add group name "rest" */
1226 grps
->nm_ind
[grps
->nr
] = restnm
;
1228 /* Assign the rest name to all atoms not currently assigned to a group */
1229 for(j
=0; (j
<natoms
); j
++) {
1230 if (cbuf
[j
] == NOGID
)
1237 if (grps
->nr
== 1) {
1238 groups
->ngrpnr
[gtype
] = 0;
1239 groups
->grpnr
[gtype
] = NULL
;
1241 groups
->ngrpnr
[gtype
] = natoms
;
1242 snew(groups
->grpnr
[gtype
],natoms
);
1243 for(j
=0; (j
<natoms
); j
++) {
1244 groups
->grpnr
[gtype
][j
] = cbuf
[j
];
1250 return (bRest
&& grptp
== egrptpPART
);
1253 static void calc_nrdf(gmx_mtop_t
*mtop
,t_inputrec
*ir
,char **gnames
)
1256 gmx_groups_t
*groups
;
1258 int natoms
,ai
,aj
,i
,j
,d
,g
,imin
,jmin
,nc
;
1260 int *nrdf
,*na_vcm
,na_tot
;
1261 double *nrdf_tc
,*nrdf_vcm
,nrdf_uc
,n_sub
=0;
1262 gmx_mtop_atomloop_all_t aloop
;
1264 int mb
,mol
,ftype
,as
;
1265 gmx_molblock_t
*molb
;
1266 gmx_moltype_t
*molt
;
1269 * First calc 3xnr-atoms for each group
1270 * then subtract half a degree of freedom for each constraint
1272 * Only atoms and nuclei contribute to the degrees of freedom...
1277 groups
= &mtop
->groups
;
1278 natoms
= mtop
->natoms
;
1280 /* Allocate one more for a possible rest group */
1281 /* We need to sum degrees of freedom into doubles,
1282 * since floats give too low nrdf's above 3 million atoms.
1284 snew(nrdf_tc
,groups
->grps
[egcTC
].nr
+1);
1285 snew(nrdf_vcm
,groups
->grps
[egcVCM
].nr
+1);
1286 snew(na_vcm
,groups
->grps
[egcVCM
].nr
+1);
1288 for(i
=0; i
<groups
->grps
[egcTC
].nr
; i
++)
1290 for(i
=0; i
<groups
->grps
[egcVCM
].nr
+1; i
++)
1294 aloop
= gmx_mtop_atomloop_all_init(mtop
);
1295 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
1297 if (atom
->ptype
== eptAtom
|| atom
->ptype
== eptNucleus
) {
1298 g
= ggrpnr(groups
,egcFREEZE
,i
);
1299 /* Double count nrdf for particle i */
1300 for(d
=0; d
<DIM
; d
++) {
1301 if (opts
->nFreeze
[g
][d
] == 0) {
1305 nrdf_tc
[ggrpnr(groups
,egcTC
,i
)] += 0.5*nrdf
[i
];
1306 nrdf_vcm
[ggrpnr(groups
,egcVCM
,i
)] += 0.5*nrdf
[i
];
1311 for(mb
=0; mb
<mtop
->nmolblock
; mb
++) {
1312 molb
= &mtop
->molblock
[mb
];
1313 molt
= &mtop
->moltype
[molb
->type
];
1314 atom
= molt
->atoms
.atom
;
1315 for(mol
=0; mol
<molb
->nmol
; mol
++) {
1316 for (ftype
=F_CONSTR
; ftype
<=F_CONSTRNC
; ftype
++) {
1317 ia
= molt
->ilist
[ftype
].iatoms
;
1318 for(i
=0; i
<molt
->ilist
[ftype
].nr
; ) {
1319 /* Subtract degrees of freedom for the constraints,
1320 * if the particles still have degrees of freedom left.
1321 * If one of the particles is a vsite or a shell, then all
1322 * constraint motion will go there, but since they do not
1323 * contribute to the constraints the degrees of freedom do not
1328 if (((atom
[ia
[1]].ptype
== eptNucleus
) ||
1329 (atom
[ia
[1]].ptype
== eptAtom
)) &&
1330 ((atom
[ia
[2]].ptype
== eptNucleus
) ||
1331 (atom
[ia
[2]].ptype
== eptAtom
))) {
1340 imin
= min(imin
,nrdf
[ai
]);
1341 jmin
= min(jmin
,nrdf
[aj
]);
1344 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5*imin
;
1345 nrdf_tc
[ggrpnr(groups
,egcTC
,aj
)] -= 0.5*jmin
;
1346 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5*imin
;
1347 nrdf_vcm
[ggrpnr(groups
,egcVCM
,aj
)] -= 0.5*jmin
;
1349 ia
+= interaction_function
[ftype
].nratoms
+1;
1350 i
+= interaction_function
[ftype
].nratoms
+1;
1353 ia
= molt
->ilist
[F_SETTLE
].iatoms
;
1354 for(i
=0; i
<molt
->ilist
[F_SETTLE
].nr
; ) {
1355 /* Subtract 1 dof from every atom in the SETTLE */
1356 for(ai
=as
+ia
[1]; ai
<as
+ia
[1]+3; ai
++) {
1357 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 1;
1358 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 1;
1363 as
+= molt
->atoms
.nr
;
1367 if (ir
->ePull
== epullCONSTRAINT
) {
1368 /* Correct nrdf for the COM constraints.
1369 * We correct using the TC and VCM group of the first atom
1370 * in the reference and pull group. If atoms in one pull group
1371 * belong to different TC or VCM groups it is anyhow difficult
1372 * to determine the optimal nrdf assignment.
1375 if (pull
->eGeom
== epullgPOS
) {
1377 for(i
=0; i
<DIM
; i
++) {
1384 for(i
=0; i
<pull
->ngrp
; i
++) {
1386 if (pull
->grp
[0].nat
> 0) {
1387 /* Subtract 1/2 dof from the reference group */
1388 ai
= pull
->grp
[0].ind
[0];
1389 if (nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] > 1) {
1390 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5;
1391 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5;
1395 /* Subtract 1/2 dof from the pulled group */
1396 ai
= pull
->grp
[1+i
].ind
[0];
1397 nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] -= 0.5*imin
;
1398 nrdf_vcm
[ggrpnr(groups
,egcVCM
,ai
)] -= 0.5*imin
;
1399 if (nrdf_tc
[ggrpnr(groups
,egcTC
,ai
)] < 0)
1400 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
)]]);
1404 if (ir
->nstcomm
!= 0) {
1405 /* Subtract 3 from the number of degrees of freedom in each vcm group
1406 * when com translation is removed and 6 when rotation is removed
1409 switch (ir
->comm_mode
) {
1411 n_sub
= ndof_com(ir
);
1418 gmx_incons("Checking comm_mode");
1421 for(i
=0; i
<groups
->grps
[egcTC
].nr
; i
++) {
1422 /* Count the number of atoms of TC group i for every VCM group */
1423 for(j
=0; j
<groups
->grps
[egcVCM
].nr
+1; j
++)
1426 for(ai
=0; ai
<natoms
; ai
++)
1427 if (ggrpnr(groups
,egcTC
,ai
) == i
) {
1428 na_vcm
[ggrpnr(groups
,egcVCM
,ai
)]++;
1431 /* Correct for VCM removal according to the fraction of each VCM
1432 * group present in this TC group.
1434 nrdf_uc
= nrdf_tc
[i
];
1436 fprintf(debug
,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
1440 for(j
=0; j
<groups
->grps
[egcVCM
].nr
+1; j
++) {
1441 if (nrdf_vcm
[j
] > n_sub
) {
1442 nrdf_tc
[i
] += nrdf_uc
*((double)na_vcm
[j
]/(double)na_tot
)*
1443 (nrdf_vcm
[j
] - n_sub
)/nrdf_vcm
[j
];
1446 fprintf(debug
," nrdf_vcm[%d] = %g, nrdf = %g\n",
1447 j
,nrdf_vcm
[j
],nrdf_tc
[i
]);
1452 for(i
=0; (i
<groups
->grps
[egcTC
].nr
); i
++) {
1453 opts
->nrdf
[i
] = nrdf_tc
[i
];
1454 if (opts
->nrdf
[i
] < 0)
1457 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
1458 gnames
[groups
->grps
[egcTC
].nm_ind
[i
]],opts
->nrdf
[i
]);
1467 static void decode_cos(char *s
,t_cosines
*cosine
,bool bTime
)
1470 char format
[STRLEN
],f1
[STRLEN
];
1481 sscanf(t
,"%d",&(cosine
->n
));
1482 if (cosine
->n
<= 0) {
1485 snew(cosine
->a
,cosine
->n
);
1486 snew(cosine
->phi
,cosine
->n
);
1488 sprintf(format
,"%%*d");
1489 for(i
=0; (i
<cosine
->n
); i
++) {
1491 strcat(f1
,"%lf%lf");
1492 if (sscanf(t
,f1
,&a
,&phi
) < 2)
1493 gmx_fatal(FARGS
,"Invalid input for electric field shift: '%s'",t
);
1496 strcat(format
,"%*lf%*lf");
1503 static bool do_egp_flag(t_inputrec
*ir
,gmx_groups_t
*groups
,
1504 const char *option
,const char *val
,int flag
)
1506 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
1507 * But since this is much larger than STRLEN, such a line can not be parsed.
1508 * The real maximum is the number of names that fit in a string: STRLEN/2.
1510 #define EGP_MAX (STRLEN/2)
1512 char *names
[EGP_MAX
];
1516 gnames
= groups
->grpname
;
1518 nelem
= str_nelem(val
,EGP_MAX
,names
);
1520 gmx_fatal(FARGS
,"The number of groups for %s is odd",option
);
1521 nr
= groups
->grps
[egcENER
].nr
;
1523 for(i
=0; i
<nelem
/2; i
++) {
1526 strcasecmp(names
[2*i
],*(gnames
[groups
->grps
[egcENER
].nm_ind
[j
]])))
1529 gmx_fatal(FARGS
,"%s in %s is not an energy group\n",
1533 strcasecmp(names
[2*i
+1],*(gnames
[groups
->grps
[egcENER
].nm_ind
[k
]])))
1536 gmx_fatal(FARGS
,"%s in %s is not an energy group\n",
1537 names
[2*i
+1],option
);
1538 if ((j
< nr
) && (k
< nr
)) {
1539 ir
->opts
.egp_flags
[nr
*j
+k
] |= flag
;
1540 ir
->opts
.egp_flags
[nr
*k
+j
] |= flag
;
1548 void do_index(const char* mdparin
, const char *ndx
,
1551 t_inputrec
*ir
,rvec
*v
)
1554 gmx_groups_t
*groups
;
1558 char warnbuf
[STRLEN
],**gnames
;
1559 int nr
,ntcg
,ntau_t
,nref_t
,nacc
,nofg
,nSA
,nSA_points
,nSA_time
,nSA_temp
;
1560 int nacg
,nfreeze
,nfrdim
,nenergy
,nvcm
,nuser
;
1561 char *ptr1
[MAXPTR
],*ptr2
[MAXPTR
],*ptr3
[MAXPTR
];
1564 bool bExcl
,bTable
,bSetTCpar
,bAnneal
,bRest
;
1565 int nQMmethod
,nQMbasis
,nQMcharge
,nQMmult
,nbSH
,nCASorb
,nCASelec
,
1566 nSAon
,nSAoff
,nSAsteps
,nQMg
,nbOPT
,nbTS
;
1569 fprintf(stderr
,"processing index file...\n");
1573 snew(grps
->index
,1);
1575 atoms_all
= gmx_mtop_global_atoms(mtop
);
1576 analyse(&atoms_all
,grps
,&gnames
,FALSE
,TRUE
);
1577 free_t_atoms(&atoms_all
,FALSE
);
1579 grps
= init_index(ndx
,&gnames
);
1582 groups
= &mtop
->groups
;
1583 natoms
= mtop
->natoms
;
1584 symtab
= &mtop
->symtab
;
1586 snew(groups
->grpname
,grps
->nr
+1);
1588 for(i
=0; (i
<grps
->nr
); i
++) {
1589 groups
->grpname
[i
] = put_symtab(symtab
,gnames
[i
]);
1591 groups
->grpname
[i
] = put_symtab(symtab
,"rest");
1593 srenew(gnames
,grps
->nr
+1);
1594 gnames
[restnm
] = *(groups
->grpname
[i
]);
1595 groups
->ngrpname
= grps
->nr
+1;
1597 set_warning_line(mdparin
,-1);
1599 ntau_t
= str_nelem(tau_t
,MAXPTR
,ptr1
);
1600 nref_t
= str_nelem(ref_t
,MAXPTR
,ptr2
);
1601 ntcg
= str_nelem(tcgrps
,MAXPTR
,ptr3
);
1602 if ((ntau_t
!= ntcg
) || (nref_t
!= ntcg
)) {
1603 gmx_fatal(FARGS
,"Invalid T coupling input: %d groups, %d ref_t values and "
1604 "%d tau_t values",ntcg
,nref_t
,ntau_t
);
1607 bSetTCpar
= (ir
->etc
|| EI_SD(ir
->eI
) || ir
->eI
==eiBD
|| EI_TPI(ir
->eI
));
1608 do_numbering(natoms
,groups
,ntcg
,ptr3
,grps
,gnames
,egcTC
,
1609 restnm
,bSetTCpar
? egrptpALL
: egrptpALL_GENREST
,bVerbose
);
1610 nr
= groups
->grps
[egcTC
].nr
;
1612 snew(ir
->opts
.nrdf
,nr
);
1613 snew(ir
->opts
.tau_t
,nr
);
1614 snew(ir
->opts
.ref_t
,nr
);
1615 if (ir
->eI
==eiBD
&& ir
->bd_fric
==0) {
1616 fprintf(stderr
,"bd_fric=0, so tau_t will be used as the inverse friction constant(s)\n");
1620 gmx_fatal(FARGS
,"Not enough ref_t and tau_t values!");
1621 for(i
=0; (i
<nr
); i
++) {
1622 ir
->opts
.tau_t
[i
]=strtod(ptr1
[i
],NULL
);
1623 if (ir
->opts
.tau_t
[i
] < 0) {
1624 gmx_fatal(FARGS
,"tau_t for group %d negative",i
);
1626 /* We check the relative magnitude of the coupling time tau_t.
1627 * V-rescale works correctly, even for tau_t=0.
1629 if ((ir
->etc
== etcBERENDSEN
|| ir
->etc
== etcNOSEHOOVER
) &&
1630 ir
->opts
.tau_t
[i
] != 0 &&
1631 ir
->opts
.tau_t
[i
] < 10*ir
->nstcalcenergy
*ir
->delta_t
) {
1632 sprintf(warn_buf
,"For proper thermostat integration tau_t (%g) should be more than an order of magnitude larger than nstcalcenergy*dt (%g)",
1633 ir
->opts
.tau_t
[i
],ir
->nstcalcenergy
*ir
->delta_t
);
1637 for(i
=0; (i
<nr
); i
++) {
1638 ir
->opts
.ref_t
[i
]=strtod(ptr2
[i
],NULL
);
1639 if (ir
->opts
.ref_t
[i
] < 0)
1640 gmx_fatal(FARGS
,"ref_t for group %d negative",i
);
1644 /* Simulated annealing for each group. There are nr groups */
1645 nSA
= str_nelem(anneal
,MAXPTR
,ptr1
);
1646 if (nSA
== 1 && (ptr1
[0][0]=='n' || ptr1
[0][0]=='N'))
1648 if(nSA
>0 && nSA
!= nr
)
1649 gmx_fatal(FARGS
,"Not enough annealing values: %d (for %d groups)\n",nSA
,nr
);
1651 snew(ir
->opts
.annealing
,nr
);
1652 snew(ir
->opts
.anneal_npoints
,nr
);
1653 snew(ir
->opts
.anneal_time
,nr
);
1654 snew(ir
->opts
.anneal_temp
,nr
);
1656 ir
->opts
.annealing
[i
]=eannNO
;
1657 ir
->opts
.anneal_npoints
[i
]=0;
1658 ir
->opts
.anneal_time
[i
]=NULL
;
1659 ir
->opts
.anneal_temp
[i
]=NULL
;
1664 if(ptr1
[i
][0]=='n' || ptr1
[i
][0]=='N') {
1665 ir
->opts
.annealing
[i
]=eannNO
;
1666 } else if(ptr1
[i
][0]=='s'|| ptr1
[i
][0]=='S') {
1667 ir
->opts
.annealing
[i
]=eannSINGLE
;
1669 } else if(ptr1
[i
][0]=='p'|| ptr1
[i
][0]=='P') {
1670 ir
->opts
.annealing
[i
]=eannPERIODIC
;
1675 /* Read the other fields too */
1676 nSA_points
= str_nelem(anneal_npoints
,MAXPTR
,ptr1
);
1678 gmx_fatal(FARGS
,"Found %d annealing_npoints values for %d groups\n",nSA_points
,nSA
);
1679 for(k
=0,i
=0;i
<nr
;i
++) {
1680 ir
->opts
.anneal_npoints
[i
]=strtol(ptr1
[i
],NULL
,10);
1681 if(ir
->opts
.anneal_npoints
[i
]==1)
1682 gmx_fatal(FARGS
,"Please specify at least a start and an end point for annealing\n");
1683 snew(ir
->opts
.anneal_time
[i
],ir
->opts
.anneal_npoints
[i
]);
1684 snew(ir
->opts
.anneal_temp
[i
],ir
->opts
.anneal_npoints
[i
]);
1685 k
+= ir
->opts
.anneal_npoints
[i
];
1688 nSA_time
= str_nelem(anneal_time
,MAXPTR
,ptr1
);
1690 gmx_fatal(FARGS
,"Found %d annealing_time values, wanter %d\n",nSA_time
,k
);
1691 nSA_temp
= str_nelem(anneal_temp
,MAXPTR
,ptr2
);
1693 gmx_fatal(FARGS
,"Found %d annealing_temp values, wanted %d\n",nSA_temp
,k
);
1695 for(i
=0,k
=0;i
<nr
;i
++) {
1697 for(j
=0;j
<ir
->opts
.anneal_npoints
[i
];j
++) {
1698 ir
->opts
.anneal_time
[i
][j
]=strtod(ptr1
[k
],NULL
);
1699 ir
->opts
.anneal_temp
[i
][j
]=strtod(ptr2
[k
],NULL
);
1701 if(ir
->opts
.anneal_time
[i
][0] > (ir
->init_t
+GMX_REAL_EPS
))
1702 gmx_fatal(FARGS
,"First time point for annealing > init_t.\n");
1705 if(ir
->opts
.anneal_time
[i
][j
]<ir
->opts
.anneal_time
[i
][j
-1])
1706 gmx_fatal(FARGS
,"Annealing timepoints out of order: t=%f comes after t=%f\n",
1707 ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_time
[i
][j
-1]);
1709 if(ir
->opts
.anneal_temp
[i
][j
]<0)
1710 gmx_fatal(FARGS
,"Found negative temperature in annealing: %f\n",ir
->opts
.anneal_temp
[i
][j
]);
1714 /* Print out some summary information, to make sure we got it right */
1715 for(i
=0,k
=0;i
<nr
;i
++) {
1716 if(ir
->opts
.annealing
[i
]!=eannNO
) {
1717 j
= groups
->grps
[egcTC
].nm_ind
[i
];
1718 fprintf(stderr
,"Simulated annealing for group %s: %s, %d timepoints\n",
1719 *(groups
->grpname
[j
]),eann_names
[ir
->opts
.annealing
[i
]],
1720 ir
->opts
.anneal_npoints
[i
]);
1721 fprintf(stderr
,"Time (ps) Temperature (K)\n");
1722 /* All terms except the last one */
1723 for(j
=0;j
<(ir
->opts
.anneal_npoints
[i
]-1);j
++)
1724 fprintf(stderr
,"%9.1f %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1726 /* Finally the last one */
1727 j
= ir
->opts
.anneal_npoints
[i
]-1;
1728 if(ir
->opts
.annealing
[i
]==eannSINGLE
)
1729 fprintf(stderr
,"%9.1f- %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1731 fprintf(stderr
,"%9.1f %5.1f\n",ir
->opts
.anneal_time
[i
][j
],ir
->opts
.anneal_temp
[i
][j
]);
1732 if(fabs(ir
->opts
.anneal_temp
[i
][j
]-ir
->opts
.anneal_temp
[i
][0])>GMX_REAL_EPS
)
1733 warning_note("There is a temperature jump when your annealing loops back.\n");
1741 if (ir
->ePull
!= epullNO
) {
1742 make_pull_groups(ir
->pull
,pull_grp
,grps
,gnames
);
1745 nacc
= str_nelem(acc
,MAXPTR
,ptr1
);
1746 nacg
= str_nelem(accgrps
,MAXPTR
,ptr2
);
1747 if (nacg
*DIM
!= nacc
)
1748 gmx_fatal(FARGS
,"Invalid Acceleration input: %d groups and %d acc. values",
1750 do_numbering(natoms
,groups
,nacg
,ptr2
,grps
,gnames
,egcACC
,
1751 restnm
,egrptpALL_GENREST
,bVerbose
);
1752 nr
= groups
->grps
[egcACC
].nr
;
1753 snew(ir
->opts
.acc
,nr
);
1756 for(i
=k
=0; (i
<nacg
); i
++)
1757 for(j
=0; (j
<DIM
); j
++,k
++)
1758 ir
->opts
.acc
[i
][j
]=strtod(ptr1
[k
],NULL
);
1760 for(j
=0; (j
<DIM
); j
++)
1761 ir
->opts
.acc
[i
][j
]=0;
1763 nfrdim
= str_nelem(frdim
,MAXPTR
,ptr1
);
1764 nfreeze
= str_nelem(freeze
,MAXPTR
,ptr2
);
1765 if (nfrdim
!= DIM
*nfreeze
)
1766 gmx_fatal(FARGS
,"Invalid Freezing input: %d groups and %d freeze values",
1768 do_numbering(natoms
,groups
,nfreeze
,ptr2
,grps
,gnames
,egcFREEZE
,
1769 restnm
,egrptpALL_GENREST
,bVerbose
);
1770 nr
= groups
->grps
[egcFREEZE
].nr
;
1772 snew(ir
->opts
.nFreeze
,nr
);
1773 for(i
=k
=0; (i
<nfreeze
); i
++)
1774 for(j
=0; (j
<DIM
); j
++,k
++) {
1775 ir
->opts
.nFreeze
[i
][j
]=(strncasecmp(ptr1
[k
],"Y",1)==0);
1776 if (!ir
->opts
.nFreeze
[i
][j
]) {
1777 if (strncasecmp(ptr1
[k
],"N",1) != 0) {
1778 sprintf(warnbuf
,"Please use Y(ES) or N(O) for freezedim only "
1779 "(not %s)", ptr1
[k
]);
1785 for(j
=0; (j
<DIM
); j
++)
1786 ir
->opts
.nFreeze
[i
][j
]=0;
1788 nenergy
=str_nelem(energy
,MAXPTR
,ptr1
);
1789 do_numbering(natoms
,groups
,nenergy
,ptr1
,grps
,gnames
,egcENER
,
1790 restnm
,egrptpALL_GENREST
,bVerbose
);
1791 add_wall_energrps(groups
,ir
->nwall
,symtab
);
1792 ir
->opts
.ngener
= groups
->grps
[egcENER
].nr
;
1793 nvcm
=str_nelem(vcm
,MAXPTR
,ptr1
);
1795 do_numbering(natoms
,groups
,nvcm
,ptr1
,grps
,gnames
,egcVCM
,
1796 restnm
,nvcm
==0 ? egrptpALL_GENREST
: egrptpPART
,bVerbose
);
1798 warning("Some atoms are not part of any center of mass motion removal group.\n"
1799 "This may lead to artifacts.\n"
1800 "In most cases one should use one group for the whole system.");
1803 /* Now we have filled the freeze struct, so we can calculate NRDF */
1804 calc_nrdf(mtop
,ir
,gnames
);
1809 /* Must check per group! */
1810 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1811 ntot
+= ir
->opts
.nrdf
[i
];
1812 if (ntot
!= (DIM
*natoms
)) {
1813 fac
= sqrt(ntot
/(DIM
*natoms
));
1815 fprintf(stderr
,"Scaling velocities by a factor of %.3f to account for constraints\n"
1816 "and removal of center of mass motion\n",fac
);
1817 for(i
=0; (i
<natoms
); i
++)
1818 svmul(fac
,v
[i
],v
[i
]);
1822 nuser
=str_nelem(user1
,MAXPTR
,ptr1
);
1823 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcUser1
,
1824 restnm
,egrptpALL_GENREST
,bVerbose
);
1825 nuser
=str_nelem(user2
,MAXPTR
,ptr1
);
1826 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcUser2
,
1827 restnm
,egrptpALL_GENREST
,bVerbose
);
1828 nuser
=str_nelem(xtc_grps
,MAXPTR
,ptr1
);
1829 do_numbering(natoms
,groups
,nuser
,ptr1
,grps
,gnames
,egcXTC
,
1830 restnm
,egrptpONE
,bVerbose
);
1831 nofg
= str_nelem(orirefitgrp
,MAXPTR
,ptr1
);
1832 do_numbering(natoms
,groups
,nofg
,ptr1
,grps
,gnames
,egcORFIT
,
1833 restnm
,egrptpALL_GENREST
,bVerbose
);
1835 /* QMMM input processing */
1836 nQMg
= str_nelem(QMMM
,MAXPTR
,ptr1
);
1837 nQMmethod
= str_nelem(QMmethod
,MAXPTR
,ptr2
);
1838 nQMbasis
= str_nelem(QMbasis
,MAXPTR
,ptr3
);
1839 if((nQMmethod
!= nQMg
)||(nQMbasis
!= nQMg
)){
1840 gmx_fatal(FARGS
,"Invalid QMMM input: %d groups %d basissets"
1841 " and %d methods\n",nQMg
,nQMbasis
,nQMmethod
);
1843 /* group rest, if any, is always MM! */
1844 do_numbering(natoms
,groups
,nQMg
,ptr1
,grps
,gnames
,egcQMMM
,
1845 restnm
,egrptpALL_GENREST
,bVerbose
);
1846 nr
= nQMg
; /*atoms->grps[egcQMMM].nr;*/
1847 ir
->opts
.ngQM
= nQMg
;
1848 snew(ir
->opts
.QMmethod
,nr
);
1849 snew(ir
->opts
.QMbasis
,nr
);
1851 /* input consists of strings: RHF CASSCF PM3 .. These need to be
1852 * converted to the corresponding enum in names.c
1854 ir
->opts
.QMmethod
[i
] = search_QMstring(ptr2
[i
],eQMmethodNR
,
1856 ir
->opts
.QMbasis
[i
] = search_QMstring(ptr3
[i
],eQMbasisNR
,
1860 nQMmult
= str_nelem(QMmult
,MAXPTR
,ptr1
);
1861 nQMcharge
= str_nelem(QMcharge
,MAXPTR
,ptr2
);
1862 nbSH
= str_nelem(bSH
,MAXPTR
,ptr3
);
1863 snew(ir
->opts
.QMmult
,nr
);
1864 snew(ir
->opts
.QMcharge
,nr
);
1865 snew(ir
->opts
.bSH
,nr
);
1868 ir
->opts
.QMmult
[i
] = strtol(ptr1
[i
],NULL
,10);
1869 ir
->opts
.QMcharge
[i
] = strtol(ptr2
[i
],NULL
,10);
1870 ir
->opts
.bSH
[i
] = (strncasecmp(ptr3
[i
],"Y",1)==0);
1873 nCASelec
= str_nelem(CASelectrons
,MAXPTR
,ptr1
);
1874 nCASorb
= str_nelem(CASorbitals
,MAXPTR
,ptr2
);
1875 snew(ir
->opts
.CASelectrons
,nr
);
1876 snew(ir
->opts
.CASorbitals
,nr
);
1878 ir
->opts
.CASelectrons
[i
]= strtol(ptr1
[i
],NULL
,10);
1879 ir
->opts
.CASorbitals
[i
] = strtol(ptr2
[i
],NULL
,10);
1881 /* special optimization options */
1883 nbOPT
= str_nelem(bOPT
,MAXPTR
,ptr1
);
1884 nbTS
= str_nelem(bTS
,MAXPTR
,ptr2
);
1885 snew(ir
->opts
.bOPT
,nr
);
1886 snew(ir
->opts
.bTS
,nr
);
1888 ir
->opts
.bOPT
[i
] = (strncasecmp(ptr1
[i
],"Y",1)==0);
1889 ir
->opts
.bTS
[i
] = (strncasecmp(ptr2
[i
],"Y",1)==0);
1891 nSAon
= str_nelem(SAon
,MAXPTR
,ptr1
);
1892 nSAoff
= str_nelem(SAoff
,MAXPTR
,ptr2
);
1893 nSAsteps
= str_nelem(SAsteps
,MAXPTR
,ptr3
);
1894 snew(ir
->opts
.SAon
,nr
);
1895 snew(ir
->opts
.SAoff
,nr
);
1896 snew(ir
->opts
.SAsteps
,nr
);
1899 ir
->opts
.SAon
[i
] = strtod(ptr1
[i
],NULL
);
1900 ir
->opts
.SAoff
[i
] = strtod(ptr2
[i
],NULL
);
1901 ir
->opts
.SAsteps
[i
] = strtol(ptr3
[i
],NULL
,10);
1903 /* end of QMMM input */
1906 for(i
=0; (i
<egcNR
); i
++) {
1907 fprintf(stderr
,"%-16s has %d element(s):",gtypes
[i
],groups
->grps
[i
].nr
);
1908 for(j
=0; (j
<groups
->grps
[i
].nr
); j
++)
1909 fprintf(stderr
," %s",*(groups
->grpname
[groups
->grps
[i
].nm_ind
[j
]]));
1910 fprintf(stderr
,"\n");
1913 nr
= groups
->grps
[egcENER
].nr
;
1914 snew(ir
->opts
.egp_flags
,nr
*nr
);
1916 bExcl
= do_egp_flag(ir
,groups
,"energygrp_excl",egpexcl
,EGP_EXCL
);
1917 if (bExcl
&& EEL_FULL(ir
->coulombtype
))
1918 warning("Can not exclude the lattice Coulomb energy between energy groups");
1920 bTable
= do_egp_flag(ir
,groups
,"energygrp_table",egptable
,EGP_TABLE
);
1921 if (bTable
&& !(ir
->vdwtype
== evdwUSER
) &&
1922 !(ir
->coulombtype
== eelUSER
) &&!(ir
->coulombtype
== eelPMEUSER
))
1923 gmx_fatal(FARGS
,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
1925 decode_cos(efield_x
,&(ir
->ex
[XX
]),FALSE
);
1926 decode_cos(efield_xt
,&(ir
->et
[XX
]),TRUE
);
1927 decode_cos(efield_y
,&(ir
->ex
[YY
]),FALSE
);
1928 decode_cos(efield_yt
,&(ir
->et
[YY
]),TRUE
);
1929 decode_cos(efield_z
,&(ir
->ex
[ZZ
]),FALSE
);
1930 decode_cos(efield_zt
,&(ir
->et
[ZZ
]),TRUE
);
1932 for(i
=0; (i
<grps
->nr
); i
++)
1942 static void check_disre(gmx_mtop_t
*mtop
)
1944 gmx_ffparams_t
*ffparams
;
1945 t_functype
*functype
;
1947 int i
,ndouble
,ftype
;
1948 int label
,old_label
;
1950 if (gmx_mtop_ftype_count(mtop
,F_DISRES
) > 0) {
1951 ffparams
= &mtop
->ffparams
;
1952 functype
= ffparams
->functype
;
1953 ip
= ffparams
->iparams
;
1956 for(i
=0; i
<ffparams
->ntypes
; i
++) {
1957 ftype
= functype
[i
];
1958 if (ftype
== F_DISRES
) {
1959 label
= ip
[i
].disres
.label
;
1960 if (label
== old_label
) {
1961 fprintf(stderr
,"Distance restraint index %d occurs twice\n",label
);
1968 gmx_fatal(FARGS
,"Found %d double distance restraint indices,\n"
1969 "probably the parameters for multiple pairs in one restraint "
1970 "are not identical\n",ndouble
);
1974 static bool absolute_reference(t_inputrec
*ir
,gmx_mtop_t
*sys
,ivec AbsRef
)
1977 gmx_mtop_ilistloop_t iloop
;
1983 for(d
=0; d
<DIM
; d
++) {
1984 AbsRef
[d
] = (d
< ndof_com(ir
) ? 0 : 1);
1986 /* Check for freeze groups */
1987 for(g
=0; g
<ir
->opts
.ngfrz
; g
++) {
1988 for(d
=0; d
<DIM
; d
++) {
1989 if (ir
->opts
.nFreeze
[g
][d
] != 0) {
1994 /* Check for position restraints */
1995 iloop
= gmx_mtop_ilistloop_init(sys
);
1996 while (gmx_mtop_ilistloop_next(iloop
,&ilist
,&nmol
)) {
1998 for(i
=0; i
<ilist
[F_POSRES
].nr
; i
+=2) {
1999 pr
= &sys
->ffparams
.iparams
[ilist
[F_POSRES
].iatoms
[i
]];
2000 for(d
=0; d
<DIM
; d
++) {
2001 if (pr
->posres
.fcA
[d
] != 0) {
2009 return (AbsRef
[XX
] != 0 && AbsRef
[YY
] != 0 && AbsRef
[ZZ
] != 0);
2012 void triple_check(const char *mdparin
,t_inputrec
*ir
,gmx_mtop_t
*sys
,int *nerror
)
2015 int i
,m
,g
,nmol
,npct
;
2017 real gdt_max
,*mgrp
,mt
;
2019 gmx_mtop_atomloop_block_t aloopb
;
2020 gmx_mtop_atomloop_all_t aloop
;
2024 if (EI_DYNAMICS(ir
->eI
) && !EI_SD(ir
->eI
) && ir
->eI
!= eiBD
&&
2025 ir
->comm_mode
== ecmNO
&&
2026 !(absolute_reference(ir
,sys
,AbsRef
) || ir
->nsteps
<= 10)) {
2027 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");
2030 if (ir
->coulombtype
== eelCUT
&& ir
->rcoulomb
> 0) {
2032 aloopb
= gmx_mtop_atomloop_block_init(sys
);
2033 while (gmx_mtop_atomloop_block_next(aloopb
,&atom
,&nmol
)) {
2034 if (atom
->q
!= 0 || atom
->qB
!= 0) {
2039 set_warning_line(mdparin
,-1);
2041 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
2042 "You might want to consider using %s electrostatics.\n",
2044 warning_note(err_buf
);
2048 /* Generalized reaction field */
2049 if (ir
->opts
.ngtc
== 0) {
2050 sprintf(err_buf
,"No temperature coupling while using coulombtype %s",
2052 CHECK(ir
->coulombtype
== eelGRF
);
2055 sprintf(err_buf
,"When using coulombtype = %s"
2056 " ref_t for temperature coupling should be > 0",
2058 CHECK((ir
->coulombtype
== eelGRF
) && (ir
->opts
.ref_t
[0] <= 0));
2061 if (ir
->eI
== eiSD1
) {
2063 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
2064 gdt_max
= max(gdt_max
,ir
->delta_t
/ir
->opts
.tau_t
[i
]);
2065 if (0.5*gdt_max
> 0.0015) {
2066 set_warning_line(mdparin
,-1);
2067 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",
2068 ei_names
[ir
->eI
],0.5*gdt_max
,ei_names
[eiSD2
]);
2074 for(i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++) {
2075 for(m
=0; (m
<DIM
); m
++) {
2076 if (fabs(ir
->opts
.acc
[i
][m
]) > 1e-6) {
2083 snew(mgrp
,sys
->groups
.grps
[egcACC
].nr
);
2084 aloop
= gmx_mtop_atomloop_all_init(sys
);
2085 while (gmx_mtop_atomloop_all_next(aloop
,&i
,&atom
)) {
2086 mgrp
[ggrpnr(&sys
->groups
,egcACC
,i
)] += atom
->m
;
2089 for(i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++) {
2090 for(m
=0; (m
<DIM
); m
++)
2091 acc
[m
] += ir
->opts
.acc
[i
][m
]*mgrp
[i
];
2094 for(m
=0; (m
<DIM
); m
++) {
2095 if (fabs(acc
[m
]) > 1e-6) {
2096 const char *dim
[DIM
] = { "X", "Y", "Z" };
2098 "Net Acceleration in %s direction, will %s be corrected\n",
2099 dim
[m
],ir
->nstcomm
!= 0 ? "" : "not");
2100 if (ir
->nstcomm
!= 0 && m
< ndof_com(ir
)) {
2102 for (i
=0; (i
<sys
->groups
.grps
[egcACC
].nr
); i
++)
2103 ir
->opts
.acc
[i
][m
] -= acc
[m
];
2110 if (ir
->efep
!= efepNO
&& ir
->sc_alpha
!= 0 &&
2111 !gmx_within_tol(sys
->ffparams
.reppow
,12.0,10*GMX_DOUBLE_EPS
)) {
2112 gmx_fatal(FARGS
,"Soft-core interactions are only supported with VdW repulsion power 12");
2115 if (ir
->ePull
!= epullNO
) {
2116 if (ir
->pull
->grp
[0].nat
== 0) {
2117 absolute_reference(ir
,sys
,AbsRef
);
2118 for(m
=0; m
<DIM
; m
++) {
2119 if (ir
->pull
->dim
[m
] && !AbsRef
[m
]) {
2120 set_warning_line(mdparin
,-1);
2121 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.");
2127 if (ir
->pull
->eGeom
== epullgDIRPBC
) {
2128 for(i
=0; i
<3; i
++) {
2129 for(m
=0; m
<=i
; m
++) {
2130 if ((ir
->epc
!= epcNO
&& ir
->compress
[i
][m
] != 0) ||
2131 ir
->deform
[i
][m
] != 0) {
2132 for(g
=1; g
<ir
->pull
->ngrp
; g
++) {
2133 if (ir
->pull
->grp
[g
].vec
[m
] != 0) {
2134 gmx_fatal(FARGS
,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir
->pull
->eGeom
),'x'+m
);
2146 void double_check(t_inputrec
*ir
,matrix box
,bool bConstr
,int *nerror
)
2152 ptr
= check_box(ir
->ePBC
,box
);
2159 if (bConstr
&& ir
->eConstrAlg
== econtSHAKE
) {
2160 if (ir
->shake_tol
<= 0.0) {
2161 fprintf(stderr
,"ERROR: shake_tol must be > 0 instead of %g\n",
2166 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1) {
2167 sprintf(warn_buf
,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
2168 if (ir
->epc
== epcNO
) {
2171 fprintf(stderr
,"ERROR: %s\n",warn_buf
);
2177 if( (ir
->eConstrAlg
== econtLINCS
) && bConstr
) {
2178 /* If we have Lincs constraints: */
2179 if(ir
->eI
==eiMD
&& ir
->etc
==etcNO
&&
2180 ir
->eConstrAlg
==econtLINCS
&& ir
->nLincsIter
==1) {
2181 sprintf(warn_buf
,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
2185 if ((ir
->eI
== eiCG
|| ir
->eI
== eiLBFGS
) && (ir
->nProjOrder
<8)) {
2186 sprintf(warn_buf
,"For accurate %s with LINCS constraints, lincs_order should be 8 or more.",ei_names
[ir
->eI
]);
2191 if (ir
->LincsWarnAngle
> 90.0) {
2192 sprintf(warn_buf
,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
2194 ir
->LincsWarnAngle
= 90.0;
2197 if (ir
->ePBC
!= epbcNONE
) {
2198 if (ir
->nstlist
== 0) {
2199 warning("With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
2201 bTWIN
= (ir
->rlistlong
> ir
->rlist
);
2202 if (ir
->ns_type
== ensGRID
) {
2203 if (sqr(ir
->rlistlong
) >= max_cutoff2(ir
->ePBC
,box
)) {
2204 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",
2205 bTWIN
? (ir
->rcoulomb
==ir
->rlistlong
? "rcoulomb" : "rvdw"):"rlist");
2209 min_size
= min(box
[XX
][XX
],min(box
[YY
][YY
],box
[ZZ
][ZZ
]));
2210 if (2*ir
->rlistlong
>= min_size
) {
2211 fprintf(stderr
,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
2214 fprintf(stderr
,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");