Bug Summary

File:gromacs/mdlib/constr.c
Location:line 350, column 5
Description:Value stored to 'nrend' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <stdlib.h>
42
43#include "gromacs/fileio/confio.h"
44#include "types/commrec.h"
45#include "constr.h"
46#include "copyrite.h"
47#include "invblock.h"
48#include "mdrun.h"
49#include "nrnb.h"
50#include "gromacs/utility/smalloc.h"
51#include "gromacs/math/vec.h"
52#include "physics.h"
53#include "names.h"
54#include "txtdump.h"
55#include "domdec.h"
56#include "gromacs/fileio/pdbio.h"
57#include "splitter.h"
58#include "mtop_util.h"
59#include "gromacs/fileio/gmxfio.h"
60#include "macros.h"
61#include "gmx_omp_nthreads.h"
62#include "gromacs/essentialdynamics/edsam.h"
63#include "gromacs/pulling/pull.h"
64
65#include "gromacs/utility/fatalerror.h"
66
67typedef struct gmx_constr {
68 int ncon_tot; /* The total number of constraints */
69 int nflexcon; /* The number of flexible constraints */
70 int n_at2con_mt; /* The size of at2con = #moltypes */
71 t_blocka *at2con_mt; /* A list of atoms to constraints */
72 int n_at2settle_mt; /* The size of at2settle = #moltypes */
73 int **at2settle_mt; /* A list of atoms to settles */
74 gmx_bool bInterCGsettles;
75 gmx_lincsdata_t lincsd; /* LINCS data */
76 gmx_shakedata_t shaked; /* SHAKE data */
77 gmx_settledata_t settled; /* SETTLE data */
78 int nblocks; /* The number of SHAKE blocks */
79 int *sblock; /* The SHAKE blocks */
80 int sblock_nalloc; /* The allocation size of sblock */
81 real *lagr; /* Lagrange multipliers for SHAKE */
82 int lagr_nalloc; /* The allocation size of lagr */
83 int maxwarn; /* The maximum number of warnings */
84 int warncount_lincs;
85 int warncount_settle;
86 gmx_edsam_t ed; /* The essential dynamics data */
87
88 tensor *vir_r_m_dr_th; /* Thread local working data */
89 int *settle_error; /* Thread local working data */
90
91 gmx_mtop_t *warn_mtop; /* Only used for printing warnings */
92} t_gmx_constr;
93
94typedef struct {
95 atom_id iatom[3];
96 atom_id blocknr;
97} t_sortblock;
98
99static void *init_vetavars(t_vetavars *vars,
100 gmx_bool constr_deriv,
101 real veta, real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
102{
103 double g;
104 int i;
105
106 /* first, set the alpha integrator variable */
107 if ((ir->opts.nrdf[0] > 0) && bPscal)
108 {
109 vars->alpha = 1.0 + DIM3/((double)ir->opts.nrdf[0]);
110 }
111 else
112 {
113 vars->alpha = 1.0;
114 }
115 g = 0.5*veta*ir->delta_t;
116 vars->rscale = exp(g)*series_sinhx(g);
117 g = -0.25*vars->alpha*veta*ir->delta_t;
118 vars->vscale = exp(g)*series_sinhx(g);
119 vars->rvscale = vars->vscale*vars->rscale;
120 vars->veta = vetanew;
121
122 if (constr_deriv)
123 {
124 snew(vars->vscale_nhc, ir->opts.ngtc)(vars->vscale_nhc) = save_calloc("vars->vscale_nhc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 124, (ir->opts.ngtc), sizeof(*(vars->vscale_nhc)));
125 if ((ekind == NULL((void*)0)) || (!bPscal))
126 {
127 for (i = 0; i < ir->opts.ngtc; i++)
128 {
129 vars->vscale_nhc[i] = 1;
130 }
131 }
132 else
133 {
134 for (i = 0; i < ir->opts.ngtc; i++)
135 {
136 vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
137 }
138 }
139 }
140 else
141 {
142 vars->vscale_nhc = NULL((void*)0);
143 }
144
145 return vars;
146}
147
148static void free_vetavars(t_vetavars *vars)
149{
150 if (vars->vscale_nhc != NULL((void*)0))
151 {
152 sfree(vars->vscale_nhc)save_free("vars->vscale_nhc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 152, (vars->vscale_nhc));
153 }
154}
155
156static int pcomp(const void *p1, const void *p2)
157{
158 int db;
159 atom_id min1, min2, max1, max2;
160 t_sortblock *a1 = (t_sortblock *)p1;
161 t_sortblock *a2 = (t_sortblock *)p2;
162
163 db = a1->blocknr-a2->blocknr;
164
165 if (db != 0)
166 {
167 return db;
168 }
169
170 min1 = min(a1->iatom[1], a1->iatom[2])(((a1->iatom[1]) < (a1->iatom[2])) ? (a1->iatom[1
]) : (a1->iatom[2]) );
171 max1 = max(a1->iatom[1], a1->iatom[2])(((a1->iatom[1]) > (a1->iatom[2])) ? (a1->iatom[1
]) : (a1->iatom[2]) );
172 min2 = min(a2->iatom[1], a2->iatom[2])(((a2->iatom[1]) < (a2->iatom[2])) ? (a2->iatom[1
]) : (a2->iatom[2]) );
173 max2 = max(a2->iatom[1], a2->iatom[2])(((a2->iatom[1]) > (a2->iatom[2])) ? (a2->iatom[1
]) : (a2->iatom[2]) );
174
175 if (min1 == min2)
176 {
177 return max1-max2;
178 }
179 else
180 {
181 return min1-min2;
182 }
183}
184
185int n_flexible_constraints(struct gmx_constr *constr)
186{
187 int nflexcon;
188
189 if (constr)
190 {
191 nflexcon = constr->nflexcon;
192 }
193 else
194 {
195 nflexcon = 0;
196 }
197
198 return nflexcon;
199}
200
201void too_many_constraint_warnings(int eConstrAlg, int warncount)
202{
203 const char *abort = "- aborting to avoid logfile runaway.\n"
204 "This normally happens when your system is not sufficiently equilibrated,"
205 "or if you are changing lambda too fast in free energy simulations.\n";
206
207 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
207,
208 "Too many %s warnings (%d)\n"
209 "If you know what you are doing you can %s"
210 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
211 "but normally it is better to fix the problem",
212 (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE", warncount,
213 (eConstrAlg == econtLINCS) ?
214 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
215}
216
217static void write_constr_pdb(const char *fn, const char *title,
218 gmx_mtop_t *mtop,
219 int start, int homenr, t_commrec *cr,
220 rvec x[], matrix box)
221{
222 char fname[STRLEN4096], format[STRLEN4096];
223 FILE *out;
224 int dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
225 gmx_domdec_t *dd;
226 char *anm, *resnm;
227
228 dd = NULL((void*)0);
229 if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes >
1)))
230 {
231 dd = cr->dd;
232 dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
233 start = 0;
234 homenr = dd_ac1;
235 }
236
237 if (PAR(cr)((cr)->nnodes > 1))
238 {
239 sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
240 }
241 else
242 {
243 sprintf(fname, "%s.pdb", fn);
244 }
245 sprintf(format, "%s\n", get_pdbformat());
246
247 out = gmx_fio_fopen(fname, "w");
248
249 fprintf(out, "TITLE %s\n", title);
250 gmx_write_pdb_box(out, -1, box);
251 for (i = start; i < start+homenr; i++)
252 {
253 if (dd != NULL((void*)0))
254 {
255 if (i >= dd->nat_home && i < dd_ac0)
256 {
257 continue;
258 }
259 ii = dd->gatindex[i];
260 }
261 else
262 {
263 ii = i;
264 }
265 gmx_mtop_atominfo_global(mtop, ii, &anm, &resnr, &resnm);
266 fprintf(out, format, "ATOM", (ii+1)%100000,
267 anm, resnm, ' ', resnr%10000, ' ',
268 10*x[i][XX0], 10*x[i][YY1], 10*x[i][ZZ2]);
269 }
270 fprintf(out, "TER\n");
271
272 gmx_fio_fclose(out);
273}
274
275static void dump_confs(FILE *fplog, gmx_int64_t step, gmx_mtop_t *mtop,
276 int start, int homenr, t_commrec *cr,
277 rvec x[], rvec xprime[], matrix box)
278{
279 char buf[256], buf2[22];
280
281 char *env = getenv("GMX_SUPPRESS_DUMP");
282 if (env)
283 {
284 return;
285 }
286
287 sprintf(buf, "step%sb", gmx_step_str(step, buf2));
288 write_constr_pdb(buf, "initial coordinates",
289 mtop, start, homenr, cr, x, box);
290 sprintf(buf, "step%sc", gmx_step_str(step, buf2));
291 write_constr_pdb(buf, "coordinates after constraining",
292 mtop, start, homenr, cr, xprime, box);
293 if (fplog)
294 {
295 fprintf(fplog, "Wrote pdb files with previous and current coordinates\n");
296 }
297 fprintf(stderrstderr, "Wrote pdb files with previous and current coordinates\n");
298}
299
300static void pr_sortblock(FILE *fp, const char *title, int nsb, t_sortblock sb[])
301{
302 int i;
303
304 fprintf(fp, "%s\n", title);
305 for (i = 0; (i < nsb); i++)
306 {
307 fprintf(fp, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
308 i, sb[i].iatom[0], sb[i].iatom[1], sb[i].iatom[2],
309 sb[i].blocknr);
310 }
311}
312
313gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
314 struct gmx_constr *constr,
315 t_idef *idef, t_inputrec *ir, gmx_ekindata_t *ekind,
316 t_commrec *cr,
317 gmx_int64_t step, int delta_step,
318 t_mdatoms *md,
319 rvec *x, rvec *xprime, rvec *min_proj,
320 gmx_bool bMolPBC, matrix box,
321 real lambda, real *dvdlambda,
322 rvec *v, tensor *vir,
323 t_nrnb *nrnb, int econq, gmx_bool bPscal,
324 real veta, real vetanew)
325{
326 gmx_bool bOK, bDump;
327 int start, homenr, nrend;
328 int i, j, d;
329 int ncons, settle_error;
330 tensor vir_r_m_dr;
331 rvec *vstor;
332 real invdt, vir_fac, t;
333 t_ilist *settle;
334 int nsettle;
335 t_pbc pbc, *pbc_null;
336 char buf[22];
337 t_vetavars vetavar;
338 int nth, th;
339
340 if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI)((ir->eI) == eiSteep || (ir->eI) == eiCG || (ir->eI)
== eiLBFGS))
341 {
342 gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted")_gmx_error("incons", "constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 342
);
343 }
344
345 bOK = TRUE1;
346 bDump = FALSE0;
347
348 start = 0;
349 homenr = md->homenr;
350 nrend = start+homenr;
Value stored to 'nrend' is never read
351
352 /* set constants for pressure control integration */
353 init_vetavars(&vetavar, econq != econqCoord,
354 veta, vetanew, ir, ekind, bPscal);
355
356 if (ir->delta_t == 0)
357 {
358 invdt = 0;
359 }
360 else
361 {
362 invdt = 1/ir->delta_t;
363 }
364
365 if (ir->efep != efepNO && EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD))
366 {
367 /* Set the constraint lengths for the step at which this configuration
368 * is meant to be. The invmasses should not be changed.
369 */
370 lambda += delta_step*ir->fepvals->delta_lambda;
371 }
372
373 if (vir != NULL((void*)0))
374 {
375 clear_mat(vir_r_m_dr);
376 }
377
378 where()_where("/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 378);
379
380 settle = &idef->il[F_SETTLE];
381 nsettle = settle->nr/(1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms));
382
383 if (nsettle > 0)
384 {
385 nth = gmx_omp_nthreads_get(emntSETTLE);
386 }
387 else
388 {
389 nth = 1;
390 }
391
392 if (nth > 1 && constr->vir_r_m_dr_th == NULL((void*)0))
393 {
394 snew(constr->vir_r_m_dr_th, nth)(constr->vir_r_m_dr_th) = save_calloc("constr->vir_r_m_dr_th"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 394
, (nth), sizeof(*(constr->vir_r_m_dr_th)));
395 snew(constr->settle_error, nth)(constr->settle_error) = save_calloc("constr->settle_error"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 395
, (nth), sizeof(*(constr->settle_error)));
396 }
397
398 settle_error = -1;
399
400 /* We do not need full pbc when constraints do not cross charge groups,
401 * i.e. when dd->constraint_comm==NULL.
402 * Note that PBC for constraints is different from PBC for bondeds.
403 * For constraints there is both forward and backward communication.
404 */
405 if (ir->ePBC != epbcNONE &&
406 (cr->dd || bMolPBC) && !(cr->dd && cr->dd->constraint_comm == NULL((void*)0)))
407 {
408 /* With pbc=screw the screw has been changed to a shift
409 * by the constraint coordinate communication routine,
410 * so that here we can use normal pbc.
411 */
412 pbc_null = set_pbc_dd(&pbc, ir->ePBC, cr->dd, FALSE0, box);
413 }
414 else
415 {
416 pbc_null = NULL((void*)0);
417 }
418
419 /* Communicate the coordinates required for the non-local constraints
420 * for LINCS and/or SETTLE.
421 */
422 if (cr->dd)
423 {
424 dd_move_x_constraints(cr->dd, box, x, xprime, econq == econqCoord);
425 }
426
427 if (constr->lincsd != NULL((void*)0))
428 {
429 bOK = constrain_lincs(fplog, bLog, bEner, ir, step, constr->lincsd, md, cr,
430 x, xprime, min_proj,
431 box, pbc_null, lambda, dvdlambda,
432 invdt, v, vir != NULL((void*)0), vir_r_m_dr,
433 econq, nrnb,
434 constr->maxwarn, &constr->warncount_lincs);
435 if (!bOK && constr->maxwarn >= 0)
436 {
437 if (fplog != NULL((void*)0))
438 {
439 fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
440 econstr_names[econtLINCS], gmx_step_str(step, buf));
441 }
442 bDump = TRUE1;
443 }
444 }
445
446 if (constr->nblocks > 0)
447 {
448 switch (econq)
449 {
450 case (econqCoord):
451 bOK = bshakef(fplog, constr->shaked,
452 md->invmass, constr->nblocks, constr->sblock,
453 idef, ir, x, xprime, nrnb,
454 constr->lagr, lambda, dvdlambda,
455 invdt, v, vir != NULL((void*)0), vir_r_m_dr,
456 constr->maxwarn >= 0, econq, &vetavar);
457 break;
458 case (econqVeloc):
459 bOK = bshakef(fplog, constr->shaked,
460 md->invmass, constr->nblocks, constr->sblock,
461 idef, ir, x, min_proj, nrnb,
462 constr->lagr, lambda, dvdlambda,
463 invdt, NULL((void*)0), vir != NULL((void*)0), vir_r_m_dr,
464 constr->maxwarn >= 0, econq, &vetavar);
465 break;
466 default:
467 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
467, "Internal error, SHAKE called for constraining something else than coordinates");
468 break;
469 }
470
471 if (!bOK && constr->maxwarn >= 0)
472 {
473 if (fplog != NULL((void*)0))
474 {
475 fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
476 econstr_names[econtSHAKE], gmx_step_str(step, buf));
477 }
478 bDump = TRUE1;
479 }
480 }
481
482 if (nsettle > 0)
483 {
484 int calcvir_atom_end;
485
486 if (vir == NULL((void*)0))
487 {
488 calcvir_atom_end = 0;
489 }
490 else
491 {
492 calcvir_atom_end = md->homenr;
493 }
494
495 switch (econq)
496 {
497 case econqCoord:
498#pragma omp parallel for num_threads(nth) schedule(static)
499 for (th = 0; th < nth; th++)
500 {
501 int start_th, end_th;
502
503 if (th > 0)
504 {
505 clear_mat(constr->vir_r_m_dr_th[th]);
506 }
507
508 start_th = (nsettle* th )/nth;
509 end_th = (nsettle*(th+1))/nth;
510 if (start_th >= 0 && end_th - start_th > 0)
511 {
512 csettle(constr->settled,
513 end_th-start_th,
514 settle->iatoms+start_th*(1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms)),
515 pbc_null,
516 x[0], xprime[0],
517 invdt, v ? v[0] : NULL((void*)0), calcvir_atom_end,
518 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
519 th == 0 ? &settle_error : &constr->settle_error[th],
520 &vetavar);
521 }
522 }
523 inc_nrnb(nrnb, eNR_SETTLE, nsettle)(nrnb)->n[eNR_SETTLE] += nsettle;
524 if (v != NULL((void*)0))
525 {
526 inc_nrnb(nrnb, eNR_CONSTR_V, nsettle*3)(nrnb)->n[eNR_CONSTR_V] += nsettle*3;
527 }
528 if (vir != NULL((void*)0))
529 {
530 inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle*3)(nrnb)->n[eNR_CONSTR_VIR] += nsettle*3;
531 }
532 break;
533 case econqVeloc:
534 case econqDeriv:
535 case econqForce:
536 case econqForceDispl:
537#pragma omp parallel for num_threads(nth) schedule(static)
538 for (th = 0; th < nth; th++)
539 {
540 int start_th, end_th;
541
542 if (th > 0)
543 {
544 clear_mat(constr->vir_r_m_dr_th[th]);
545 }
546
547 start_th = (nsettle* th )/nth;
548 end_th = (nsettle*(th+1))/nth;
549
550 if (start_th >= 0 && end_th - start_th > 0)
551 {
552 settle_proj(constr->settled, econq,
553 end_th-start_th,
554 settle->iatoms+start_th*(1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms)),
555 pbc_null,
556 x,
557 xprime, min_proj, calcvir_atom_end,
558 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
559 &vetavar);
560 }
561 }
562 /* This is an overestimate */
563 inc_nrnb(nrnb, eNR_SETTLE, nsettle)(nrnb)->n[eNR_SETTLE] += nsettle;
564 break;
565 case econqDeriv_FlexCon:
566 /* Nothing to do, since the are no flexible constraints in settles */
567 break;
568 default:
569 gmx_incons("Unknown constraint quantity for settle")_gmx_error("incons", "Unknown constraint quantity for settle"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 569
);
570 }
571 }
572
573 if (settle->nr > 0)
574 {
575 /* Combine virial and error info of the other threads */
576 for (i = 1; i < nth; i++)
577 {
578 m_add(vir_r_m_dr, constr->vir_r_m_dr_th[i], vir_r_m_dr);
579 settle_error = constr->settle_error[i];
580 }
581
582 if (econq == econqCoord && settle_error >= 0)
583 {
584 bOK = FALSE0;
585 if (constr->maxwarn >= 0)
586 {
587 char buf[256];
588 sprintf(buf,
589 "\nstep " "%"GMX_PRId64"l" "d" ": Water molecule starting at atom %d can not be "
590 "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
591 step, ddglatnr(cr->dd, settle->iatoms[settle_error*(1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms))+1]));
592 if (fplog)
593 {
594 fprintf(fplog, "%s", buf);
595 }
596 fprintf(stderrstderr, "%s", buf);
597 constr->warncount_settle++;
598 if (constr->warncount_settle > constr->maxwarn)
599 {
600 too_many_constraint_warnings(-1, constr->warncount_settle);
601 }
602 bDump = TRUE1;
603 }
604 }
605 }
606
607 free_vetavars(&vetavar);
608
609 if (vir != NULL((void*)0))
610 {
611 switch (econq)
612 {
613 case econqCoord:
614 vir_fac = 0.5/(ir->delta_t*ir->delta_t);
615 break;
616 case econqVeloc:
617 vir_fac = 0.5/ir->delta_t;
618 break;
619 case econqForce:
620 case econqForceDispl:
621 vir_fac = 0.5;
622 break;
623 default:
624 vir_fac = 0;
625 gmx_incons("Unsupported constraint quantity for virial")_gmx_error("incons", "Unsupported constraint quantity for virial"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 625
);
626 }
627
628 if (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK))
629 {
630 vir_fac *= 2; /* only constraining over half the distance here */
631 }
632 for (i = 0; i < DIM3; i++)
633 {
634 for (j = 0; j < DIM3; j++)
635 {
636 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
637 }
638 }
639 }
640
641 if (bDump)
642 {
643 dump_confs(fplog, step, constr->warn_mtop, start, homenr, cr, x, xprime, box);
644 }
645
646 if (econq == econqCoord)
647 {
648 if (ir->ePull == epullCONSTRAINT)
649 {
650 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD))
651 {
652 t = ir->init_t + (step + delta_step)*ir->delta_t;
653 }
654 else
655 {
656 t = ir->init_t;
657 }
658 set_pbc(&pbc, ir->ePBC, box);
659 pull_constraint(ir->pull, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir);
660 }
661 if (constr->ed && delta_step > 0)
662 {
663 /* apply the essential dynamcs constraints here */
664 do_edsam(ir, step, cr, xprime, v, box, constr->ed);
665 }
666 }
667
668 return bOK;
669}
670
671real *constr_rmsd_data(struct gmx_constr *constr)
672{
673 if (constr->lincsd)
674 {
675 return lincs_rmsd_data(constr->lincsd);
676 }
677 else
678 {
679 return NULL((void*)0);
680 }
681}
682
683real constr_rmsd(struct gmx_constr *constr, gmx_bool bSD2)
684{
685 if (constr->lincsd)
686 {
687 return lincs_rmsd(constr->lincsd, bSD2);
688 }
689 else
690 {
691 return 0;
692 }
693}
694
695static void make_shake_sblock_serial(struct gmx_constr *constr,
696 t_idef *idef, t_mdatoms *md)
697{
698 int i, j, m, ncons;
699 int bstart, bnr;
700 t_blocka sblocks;
701 t_sortblock *sb;
702 t_iatom *iatom;
703 atom_id *inv_sblock;
704
705 /* Since we are processing the local topology,
706 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
707 */
708 ncons = idef->il[F_CONSTR].nr/3;
709
710 init_blocka(&sblocks);
711 gen_sblocks(NULL((void*)0), 0, md->homenr, idef, &sblocks, FALSE0);
712
713 /*
714 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
715 nblocks=blocks->multinr[idef->nodeid] - bstart;
716 */
717 bstart = 0;
718 constr->nblocks = sblocks.nr;
719 if (debug)
720 {
721 fprintf(debug, "ncons: %d, bstart: %d, nblocks: %d\n",
722 ncons, bstart, constr->nblocks);
723 }
724
725 /* Calculate block number for each atom */
726 inv_sblock = make_invblocka(&sblocks, md->nr);
727
728 done_blocka(&sblocks);
729
730 /* Store the block number in temp array and
731 * sort the constraints in order of the sblock number
732 * and the atom numbers, really sorting a segment of the array!
733 */
734#ifdef DEBUGIDEF
735 pr_idef(fplog, 0, "Before Sort", idef);
736#endif
737 iatom = idef->il[F_CONSTR].iatoms;
738 snew(sb, ncons)(sb) = save_calloc("sb", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 738, (ncons), sizeof(*(sb)));
739 for (i = 0; (i < ncons); i++, iatom += 3)
740 {
741 for (m = 0; (m < 3); m++)
742 {
743 sb[i].iatom[m] = iatom[m];
744 }
745 sb[i].blocknr = inv_sblock[iatom[1]];
746 }
747
748 /* Now sort the blocks */
749 if (debug)
750 {
751 pr_sortblock(debug, "Before sorting", ncons, sb);
752 fprintf(debug, "Going to sort constraints\n");
753 }
754
755 qsort(sb, ncons, (size_t)sizeof(*sb), pcomp);
756
757 if (debug)
758 {
759 pr_sortblock(debug, "After sorting", ncons, sb);
760 }
761
762 iatom = idef->il[F_CONSTR].iatoms;
763 for (i = 0; (i < ncons); i++, iatom += 3)
764 {
765 for (m = 0; (m < 3); m++)
766 {
767 iatom[m] = sb[i].iatom[m];
768 }
769 }
770#ifdef DEBUGIDEF
771 pr_idef(fplog, 0, "After Sort", idef);
772#endif
773
774 j = 0;
775 snew(constr->sblock, constr->nblocks+1)(constr->sblock) = save_calloc("constr->sblock", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 775, (constr->nblocks+1), sizeof(*(constr->sblock)));
776 bnr = -2;
777 for (i = 0; (i < ncons); i++)
778 {
779 if (sb[i].blocknr != bnr)
780 {
781 bnr = sb[i].blocknr;
782 constr->sblock[j++] = 3*i;
783 }
784 }
785 /* Last block... */
786 constr->sblock[j++] = 3*ncons;
787
788 if (j != (constr->nblocks+1))
789 {
790 fprintf(stderrstderr, "bstart: %d\n", bstart);
791 fprintf(stderrstderr, "j: %d, nblocks: %d, ncons: %d\n",
792 j, constr->nblocks, ncons);
793 for (i = 0; (i < ncons); i++)
794 {
795 fprintf(stderrstderr, "i: %5d sb[i].blocknr: %5u\n", i, sb[i].blocknr);
796 }
797 for (j = 0; (j <= constr->nblocks); j++)
798 {
799 fprintf(stderrstderr, "sblock[%3d]=%5d\n", j, (int)constr->sblock[j]);
800 }
801 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
801, "DEATH HORROR: "
802 "sblocks does not match idef->il[F_CONSTR]");
803 }
804 sfree(sb)save_free("sb", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 804, (sb));
805 sfree(inv_sblock)save_free("inv_sblock", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 805, (inv_sblock));
806}
807
808static void make_shake_sblock_dd(struct gmx_constr *constr,
809 t_ilist *ilcon, t_block *cgs,
810 gmx_domdec_t *dd)
811{
812 int ncons, c, cg;
813 t_iatom *iatom;
814
815 if (dd->ncg_home+1 > constr->sblock_nalloc)
816 {
817 constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
818 srenew(constr->sblock, constr->sblock_nalloc)(constr->sblock) = save_realloc("constr->sblock", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 818, (constr->sblock), (constr->sblock_nalloc), sizeof
(*(constr->sblock)));
819 }
820
821 ncons = ilcon->nr/3;
822 iatom = ilcon->iatoms;
823 constr->nblocks = 0;
824 cg = 0;
825 for (c = 0; c < ncons; c++)
826 {
827 if (c == 0 || iatom[1] >= cgs->index[cg+1])
828 {
829 constr->sblock[constr->nblocks++] = 3*c;
830 while (iatom[1] >= cgs->index[cg+1])
831 {
832 cg++;
833 }
834 }
835 iatom += 3;
836 }
837 constr->sblock[constr->nblocks] = 3*ncons;
838}
839
840t_blocka make_at2con(int start, int natoms,
841 t_ilist *ilist, t_iparams *iparams,
842 gmx_bool bDynamics, int *nflexiblecons)
843{
844 int *count, ncon, con, con_tot, nflexcon, ftype, i, a;
845 t_iatom *ia;
846 t_blocka at2con;
847 gmx_bool bFlexCon;
848
849 snew(count, natoms)(count) = save_calloc("count", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 849, (natoms), sizeof(*(count)));
850 nflexcon = 0;
851 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
852 {
853 ncon = ilist[ftype].nr/3;
854 ia = ilist[ftype].iatoms;
855 for (con = 0; con < ncon; con++)
856 {
857 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
858 iparams[ia[0]].constr.dB == 0);
859 if (bFlexCon)
860 {
861 nflexcon++;
862 }
863 if (bDynamics || !bFlexCon)
864 {
865 for (i = 1; i < 3; i++)
866 {
867 a = ia[i] - start;
868 count[a]++;
869 }
870 }
871 ia += 3;
872 }
873 }
874 *nflexiblecons = nflexcon;
875
876 at2con.nr = natoms;
877 at2con.nalloc_index = at2con.nr+1;
878 snew(at2con.index, at2con.nalloc_index)(at2con.index) = save_calloc("at2con.index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 878, (at2con.nalloc_index), sizeof(*(at2con.index)));
879 at2con.index[0] = 0;
880 for (a = 0; a < natoms; a++)
881 {
882 at2con.index[a+1] = at2con.index[a] + count[a];
883 count[a] = 0;
884 }
885 at2con.nra = at2con.index[natoms];
886 at2con.nalloc_a = at2con.nra;
887 snew(at2con.a, at2con.nalloc_a)(at2con.a) = save_calloc("at2con.a", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 887, (at2con.nalloc_a), sizeof(*(at2con.a)));
888
889 /* The F_CONSTRNC constraints have constraint numbers
890 * that continue after the last F_CONSTR constraint.
891 */
892 con_tot = 0;
893 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
894 {
895 ncon = ilist[ftype].nr/3;
896 ia = ilist[ftype].iatoms;
897 for (con = 0; con < ncon; con++)
898 {
899 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
900 iparams[ia[0]].constr.dB == 0);
901 if (bDynamics || !bFlexCon)
902 {
903 for (i = 1; i < 3; i++)
904 {
905 a = ia[i] - start;
906 at2con.a[at2con.index[a]+count[a]++] = con_tot;
907 }
908 }
909 con_tot++;
910 ia += 3;
911 }
912 }
913
914 sfree(count)save_free("count", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 914, (count));
915
916 return at2con;
917}
918
919static int *make_at2settle(int natoms, const t_ilist *ilist)
920{
921 int *at2s;
922 int a, stride, s;
923
924 snew(at2s, natoms)(at2s) = save_calloc("at2s", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 924, (natoms), sizeof(*(at2s)));
925 /* Set all to no settle */
926 for (a = 0; a < natoms; a++)
927 {
928 at2s[a] = -1;
929 }
930
931 stride = 1 + NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms);
932
933 for (s = 0; s < ilist->nr; s += stride)
934 {
935 at2s[ilist->iatoms[s+1]] = s/stride;
936 at2s[ilist->iatoms[s+2]] = s/stride;
937 at2s[ilist->iatoms[s+3]] = s/stride;
938 }
939
940 return at2s;
941}
942
943void set_constraints(struct gmx_constr *constr,
944 gmx_localtop_t *top, t_inputrec *ir,
945 t_mdatoms *md, t_commrec *cr)
946{
947 t_idef *idef;
948 int ncons;
949 t_ilist *settle;
950 int iO, iH;
951
952 idef = &top->idef;
953
954 if (constr->ncon_tot > 0)
955 {
956 /* We are using the local topology,
957 * so there are only F_CONSTR constraints.
958 */
959 ncons = idef->il[F_CONSTR].nr/3;
960
961 /* With DD we might also need to call LINCS with ncons=0 for
962 * communicating coordinates to other nodes that do have constraints.
963 */
964 if (ir->eConstrAlg == econtLINCS)
965 {
966 set_lincs(idef, md, EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD), cr, constr->lincsd);
967 }
968 if (ir->eConstrAlg == econtSHAKE)
969 {
970 if (cr->dd)
971 {
972 make_shake_sblock_dd(constr, &idef->il[F_CONSTR], &top->cgs, cr->dd);
973 }
974 else
975 {
976 make_shake_sblock_serial(constr, idef, md);
977 }
978 if (ncons > constr->lagr_nalloc)
979 {
980 constr->lagr_nalloc = over_alloc_dd(ncons);
981 srenew(constr->lagr, constr->lagr_nalloc)(constr->lagr) = save_realloc("constr->lagr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 981, (constr->lagr), (constr->lagr_nalloc), sizeof(*(
constr->lagr)));
982 }
983 }
984 }
985
986 if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL((void*)0))
987 {
988 settle = &idef->il[F_SETTLE];
989 iO = settle->iatoms[1];
990 iH = settle->iatoms[2];
991 constr->settled =
992 settle_init(md->massT[iO], md->massT[iH],
993 md->invmass[iO], md->invmass[iH],
994 idef->iparams[settle->iatoms[0]].settle.doh,
995 idef->iparams[settle->iatoms[0]].settle.dhh);
996 }
997
998 /* Make a selection of the local atoms for essential dynamics */
999 if (constr->ed && cr->dd)
1000 {
1001 dd_make_local_ed_indices(cr->dd, constr->ed);
1002 }
1003}
1004
1005static void constr_recur(t_blocka *at2con,
1006 t_ilist *ilist, t_iparams *iparams, gmx_bool bTopB,
1007 int at, int depth, int nc, int *path,
1008 real r0, real r1, real *r2max,
1009 int *count)
1010{
1011 int ncon1;
1012 t_iatom *ia1, *ia2;
1013 int c, con, a1;
1014 gmx_bool bUse;
1015 t_iatom *ia;
1016 real len, rn0, rn1;
1017
1018 (*count)++;
1019
1020 ncon1 = ilist[F_CONSTR].nr/3;
1021 ia1 = ilist[F_CONSTR].iatoms;
1022 ia2 = ilist[F_CONSTRNC].iatoms;
1023
1024 /* Loop over all constraints connected to this atom */
1025 for (c = at2con->index[at]; c < at2con->index[at+1]; c++)
1026 {
1027 con = at2con->a[c];
1028 /* Do not walk over already used constraints */
1029 bUse = TRUE1;
1030 for (a1 = 0; a1 < depth; a1++)
1031 {
1032 if (con == path[a1])
1033 {
1034 bUse = FALSE0;
1035 }
1036 }
1037 if (bUse)
1038 {
1039 ia = constr_iatomptr(ncon1, ia1, ia2, con)((con) < (ncon1) ? (ia1)+(con)*3 : (ia2)+(con-ncon1)*3);
1040 /* Flexible constraints currently have length 0, which is incorrect */
1041 if (!bTopB)
1042 {
1043 len = iparams[ia[0]].constr.dA;
1044 }
1045 else
1046 {
1047 len = iparams[ia[0]].constr.dB;
1048 }
1049 /* In the worst case the bond directions alternate */
1050 if (nc % 2 == 0)
1051 {
1052 rn0 = r0 + len;
1053 rn1 = r1;
1054 }
1055 else
1056 {
1057 rn0 = r0;
1058 rn1 = r1 + len;
1059 }
1060 /* Assume angles of 120 degrees between all bonds */
1061 if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max)
1062 {
1063 *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
1064 if (debug)
1065 {
1066 fprintf(debug, "Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0, rn1, sqrt(*r2max));
1067 for (a1 = 0; a1 < depth; a1++)
1068 {
1069 fprintf(debug, " %d %5.3f",
1070 path[a1],
1071 iparams[constr_iatomptr(ncon1, ia1, ia2, con)((con) < (ncon1) ? (ia1)+(con)*3 : (ia2)+(con-ncon1)*3)[0]].constr.dA);
1072 }
1073 fprintf(debug, " %d %5.3f\n", con, len);
1074 }
1075 }
1076 /* Limit the number of recursions to 1000*nc,
1077 * so a call does not take more than a second,
1078 * even for highly connected systems.
1079 */
1080 if (depth + 1 < nc && *count < 1000*nc)
1081 {
1082 if (ia[1] == at)
1083 {
1084 a1 = ia[2];
1085 }
1086 else
1087 {
1088 a1 = ia[1];
1089 }
1090 /* Recursion */
1091 path[depth] = con;
1092 constr_recur(at2con, ilist, iparams,
1093 bTopB, a1, depth+1, nc, path, rn0, rn1, r2max, count);
1094 path[depth] = -1;
1095 }
1096 }
1097 }
1098}
1099
1100static real constr_r_max_moltype(gmx_moltype_t *molt, t_iparams *iparams,
1101 t_inputrec *ir)
1102{
1103 int natoms, nflexcon, *path, at, count;
1104
1105 t_blocka at2con;
1106 real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
1107
1108 if (molt->ilist[F_CONSTR].nr == 0 &&
1109 molt->ilist[F_CONSTRNC].nr == 0)
1110 {
1111 return 0;
1112 }
1113
1114 natoms = molt->atoms.nr;
1115
1116 at2con = make_at2con(0, natoms, molt->ilist, iparams,
1117 EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD), &nflexcon);
1118 snew(path, 1+ir->nProjOrder)(path) = save_calloc("path", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1118, (1+ir->nProjOrder), sizeof(*(path)));
1119 for (at = 0; at < 1+ir->nProjOrder; at++)
1120 {
1121 path[at] = -1;
1122 }
1123
1124 r2maxA = 0;
1125 for (at = 0; at < natoms; at++)
1126 {
1127 r0 = 0;
1128 r1 = 0;
1129
1130 count = 0;
1131 constr_recur(&at2con, molt->ilist, iparams,
1132 FALSE0, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxA, &count);
1133 }
1134 if (ir->efep == efepNO)
1135 {
1136 rmax = sqrt(r2maxA);
1137 }
1138 else
1139 {
1140 r2maxB = 0;
1141 for (at = 0; at < natoms; at++)
1142 {
1143 r0 = 0;
1144 r1 = 0;
1145 count = 0;
1146 constr_recur(&at2con, molt->ilist, iparams,
1147 TRUE1, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxB, &count);
1148 }
1149 lam0 = ir->fepvals->init_lambda;
1150 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD))
1151 {
1152 lam0 += ir->init_step*ir->fepvals->delta_lambda;
1153 }
1154 rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
1155 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD))
1156 {
1157 lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
1158 rmax = max(rmax, (1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB))(((rmax) > ((1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB))) ?
(rmax) : ((1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB)) );
1159 }
1160 }
1161
1162 done_blocka(&at2con);
1163 sfree(path)save_free("path", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1163, (path));
1164
1165 return rmax;
1166}
1167
1168real constr_r_max(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir)
1169{
1170 int mt;
1171 real rmax;
1172
1173 rmax = 0;
1174 for (mt = 0; mt < mtop->nmoltype; mt++)
1175 {
1176 rmax = max(rmax,(((rmax) > (constr_r_max_moltype(&mtop->moltype[mt]
, mtop->ffparams.iparams, ir))) ? (rmax) : (constr_r_max_moltype
(&mtop->moltype[mt], mtop->ffparams.iparams, ir)) )
1177 constr_r_max_moltype(&mtop->moltype[mt],(((rmax) > (constr_r_max_moltype(&mtop->moltype[mt]
, mtop->ffparams.iparams, ir))) ? (rmax) : (constr_r_max_moltype
(&mtop->moltype[mt], mtop->ffparams.iparams, ir)) )
1178 mtop->ffparams.iparams, ir))(((rmax) > (constr_r_max_moltype(&mtop->moltype[mt]
, mtop->ffparams.iparams, ir))) ? (rmax) : (constr_r_max_moltype
(&mtop->moltype[mt], mtop->ffparams.iparams, ir)) );
1179 }
1180
1181 if (fplog)
1182 {
1183 fprintf(fplog, "Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n", 1+ir->nProjOrder, rmax);
1184 }
1185
1186 return rmax;
1187}
1188
1189gmx_constr_t init_constraints(FILE *fplog,
1190 gmx_mtop_t *mtop, t_inputrec *ir,
1191 gmx_edsam_t ed, t_state *state,
1192 t_commrec *cr)
1193{
1194 int ncon, nset, nmol, settle_type, i, natoms, mt, nflexcon;
1195 struct gmx_constr *constr;
1196 char *env;
1197 t_ilist *ilist;
1198 gmx_mtop_ilistloop_t iloop;
1199
1200 ncon =
1201 gmx_mtop_ftype_count(mtop, F_CONSTR) +
1202 gmx_mtop_ftype_count(mtop, F_CONSTRNC);
1203 nset = gmx_mtop_ftype_count(mtop, F_SETTLE);
1204
1205 if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL((void*)0))
1206 {
1207 return NULL((void*)0);
1208 }
1209
1210 snew(constr, 1)(constr) = save_calloc("constr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1210, (1), sizeof(*(constr)));
1211
1212 constr->ncon_tot = ncon;
1213 constr->nflexcon = 0;
1214 if (ncon > 0)
1215 {
1216 constr->n_at2con_mt = mtop->nmoltype;
1217 snew(constr->at2con_mt, constr->n_at2con_mt)(constr->at2con_mt) = save_calloc("constr->at2con_mt", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1217, (constr->n_at2con_mt), sizeof(*(constr->at2con_mt
)));
1218 for (mt = 0; mt < mtop->nmoltype; mt++)
1219 {
1220 constr->at2con_mt[mt] = make_at2con(0, mtop->moltype[mt].atoms.nr,
1221 mtop->moltype[mt].ilist,
1222 mtop->ffparams.iparams,
1223 EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD), &nflexcon);
1224 for (i = 0; i < mtop->nmolblock; i++)
1225 {
1226 if (mtop->molblock[i].type == mt)
1227 {
1228 constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
1229 }
1230 }
1231 }
1232
1233 if (constr->nflexcon > 0)
1234 {
1235 if (fplog)
1236 {
1237 fprintf(fplog, "There are %d flexible constraints\n",
1238 constr->nflexcon);
1239 if (ir->fc_stepsize == 0)
1240 {
1241 fprintf(fplog, "\n"
1242 "WARNING: step size for flexible constraining = 0\n"
1243 " All flexible constraints will be rigid.\n"
1244 " Will try to keep all flexible constraints at their original length,\n"
1245 " but the lengths may exhibit some drift.\n\n");
1246 constr->nflexcon = 0;
1247 }
1248 }
1249 if (constr->nflexcon > 0)
1250 {
1251 please_cite(fplog, "Hess2002");
1252 }
1253 }
1254
1255 if (ir->eConstrAlg == econtLINCS)
1256 {
1257 constr->lincsd = init_lincs(fplog, mtop,
1258 constr->nflexcon, constr->at2con_mt,
1259 DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes >
1)) && cr->dd->bInterCGcons,
1260 ir->nLincsIter, ir->nProjOrder);
1261 }
1262
1263 if (ir->eConstrAlg == econtSHAKE)
1264 {
1265 if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes >
1)) && cr->dd->bInterCGcons)
1266 {
1267 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
1267, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1268 }
1269 if (constr->nflexcon)
1270 {
1271 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
1271, "For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1272 }
1273 please_cite(fplog, "Ryckaert77a");
1274 if (ir->bShakeSOR)
1275 {
1276 please_cite(fplog, "Barth95a");
1277 }
1278
1279 constr->shaked = shake_init();
1280 }
1281 }
1282
1283 if (nset > 0)
1284 {
1285 please_cite(fplog, "Miyamoto92a");
1286
1287 constr->bInterCGsettles = inter_charge_group_settles(mtop);
1288
1289 /* Check that we have only one settle type */
1290 settle_type = -1;
1291 iloop = gmx_mtop_ilistloop_init(mtop);
1292 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
1293 {
1294 for (i = 0; i < ilist[F_SETTLE].nr; i += 4)
1295 {
1296 if (settle_type == -1)
1297 {
1298 settle_type = ilist[F_SETTLE].iatoms[i];
1299 }
1300 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1301 {
1302 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c",
1302,
1303 "The [molecules] section of your topology specifies more than one block of\n"
1304 "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1305 "are trying to partition your solvent into different *groups* (e.g. for\n"
1306 "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1307 "files specify groups. Otherwise, you may wish to change the least-used\n"
1308 "block of molecules with SETTLE constraints into 3 normal constraints.");
1309 }
1310 }
1311 }
1312
1313 constr->n_at2settle_mt = mtop->nmoltype;
1314 snew(constr->at2settle_mt, constr->n_at2settle_mt)(constr->at2settle_mt) = save_calloc("constr->at2settle_mt"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c", 1314
, (constr->n_at2settle_mt), sizeof(*(constr->at2settle_mt
)));
1315 for (mt = 0; mt < mtop->nmoltype; mt++)
1316 {
1317 constr->at2settle_mt[mt] =
1318 make_at2settle(mtop->moltype[mt].atoms.nr,
1319 &mtop->moltype[mt].ilist[F_SETTLE]);
1320 }
1321 }
1322
1323 constr->maxwarn = 999;
1324 env = getenv("GMX_MAXCONSTRWARN");
1325 if (env)
1326 {
1327 constr->maxwarn = 0;
1328 sscanf(env, "%d", &constr->maxwarn);
1329 if (fplog)
1330 {
1331 fprintf(fplog,
1332 "Setting the maximum number of constraint warnings to %d\n",
1333 constr->maxwarn);
1334 }
1335 if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
1336 {
1337 fprintf(stderrstderr,
1338 "Setting the maximum number of constraint warnings to %d\n",
1339 constr->maxwarn);
1340 }
1341 }
1342 if (constr->maxwarn < 0 && fplog)
1343 {
1344 fprintf(fplog, "maxwarn < 0, will not stop on constraint errors\n");
1345 }
1346 constr->warncount_lincs = 0;
1347 constr->warncount_settle = 0;
1348
1349 /* Initialize the essential dynamics sampling.
1350 * Put the pointer to the ED struct in constr */
1351 constr->ed = ed;
1352 if (ed != NULL((void*)0) || state->edsamstate.nED > 0)
1353 {
1354 init_edsam(mtop, ir, cr, ed, state->x, state->box, &state->edsamstate);
1355 }
1356
1357 constr->warn_mtop = mtop;
1358
1359 return constr;
1360}
1361
1362const t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1363{
1364 return constr->at2con_mt;
1365}
1366
1367const int **atom2settle_moltype(gmx_constr_t constr)
1368{
1369 return (const int **)constr->at2settle_mt;
1370}
1371
1372
1373gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
1374{
1375 const gmx_moltype_t *molt;
1376 const t_block *cgs;
1377 const t_ilist *il;
1378 int mb;
1379 int nat, *at2cg, cg, a, ftype, i;
1380 gmx_bool bInterCG;
1381
1382 bInterCG = FALSE0;
1383 for (mb = 0; mb < mtop->nmolblock && !bInterCG; mb++)
1384 {
1385 molt = &mtop->moltype[mtop->molblock[mb].type];
1386
1387 if (molt->ilist[F_CONSTR].nr > 0 ||
1388 molt->ilist[F_CONSTRNC].nr > 0 ||
1389 molt->ilist[F_SETTLE].nr > 0)
1390 {
1391 cgs = &molt->cgs;
1392 snew(at2cg, molt->atoms.nr)(at2cg) = save_calloc("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1392, (molt->atoms.nr), sizeof(*(at2cg)));
1393 for (cg = 0; cg < cgs->nr; cg++)
1394 {
1395 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1396 {
1397 at2cg[a] = cg;
1398 }
1399 }
1400
1401 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
1402 {
1403 il = &molt->ilist[ftype];
1404 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype)(interaction_function[(ftype)].nratoms))
1405 {
1406 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1407 {
1408 bInterCG = TRUE1;
1409 }
1410 }
1411 }
1412
1413 sfree(at2cg)save_free("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1413, (at2cg));
1414 }
1415 }
1416
1417 return bInterCG;
1418}
1419
1420gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop)
1421{
1422 const gmx_moltype_t *molt;
1423 const t_block *cgs;
1424 const t_ilist *il;
1425 int mb;
1426 int nat, *at2cg, cg, a, ftype, i;
1427 gmx_bool bInterCG;
1428
1429 bInterCG = FALSE0;
1430 for (mb = 0; mb < mtop->nmolblock && !bInterCG; mb++)
1431 {
1432 molt = &mtop->moltype[mtop->molblock[mb].type];
1433
1434 if (molt->ilist[F_SETTLE].nr > 0)
1435 {
1436 cgs = &molt->cgs;
1437 snew(at2cg, molt->atoms.nr)(at2cg) = save_calloc("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1437, (molt->atoms.nr), sizeof(*(at2cg)));
1438 for (cg = 0; cg < cgs->nr; cg++)
1439 {
1440 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1441 {
1442 at2cg[a] = cg;
1443 }
1444 }
1445
1446 for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
1447 {
1448 il = &molt->ilist[ftype];
1449 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms))
1450 {
1451 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1452 at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1453 {
1454 bInterCG = TRUE1;
1455 }
1456 }
1457 }
1458
1459 sfree(at2cg)save_free("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/constr.c"
, 1459, (at2cg));
1460 }
1461 }
1462
1463 return bInterCG;
1464}
1465
1466/* helper functions for andersen temperature control, because the
1467 * gmx_constr construct is only defined in constr.c. Return the list
1468 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
1469
1470extern int *get_sblock(struct gmx_constr *constr)
1471{
1472 return constr->sblock;
1473}
1474
1475extern int get_nblocks(struct gmx_constr *constr)
1476{
1477 return constr->nblocks;
1478}