1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
40 #include "gromacs/fileio/confio.h"
53 #include "gromacs/fileio/pdbio.h"
56 #include "mtop_util.h"
57 #include "gromacs/fileio/gmxfio.h"
59 #include "gmx_omp_nthreads.h"
61 typedef struct gmx_constr {
62 int ncon_tot; /* The total number of constraints */
63 int nflexcon; /* The number of flexible constraints */
64 int n_at2con_mt; /* The size of at2con = #moltypes */
65 t_blocka *at2con_mt; /* A list of atoms to constraints */
66 int n_at2settle_mt; /* The size of at2settle = #moltypes */
67 int **at2settle_mt; /* A list of atoms to settles */
68 gmx_bool bInterCGsettles;
69 gmx_lincsdata_t lincsd; /* LINCS data */
70 gmx_shakedata_t shaked; /* SHAKE data */
71 gmx_settledata_t settled; /* SETTLE data */
72 int nblocks; /* The number of SHAKE blocks */
73 int *sblock; /* The SHAKE blocks */
74 int sblock_nalloc; /* The allocation size of sblock */
75 real *lagr; /* Lagrange multipliers for SHAKE */
76 int lagr_nalloc; /* The allocation size of lagr */
77 int maxwarn; /* The maximum number of warnings */
80 gmx_edsam_t ed; /* The essential dynamics data */
82 tensor *vir_r_m_dr_th; /* Thread local working data */
83 int *settle_error; /* Thread local working data */
85 gmx_mtop_t *warn_mtop; /* Only used for printing warnings */
93 static void *init_vetavars(t_vetavars *vars,
94 gmx_bool constr_deriv,
95 real veta, real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
100 /* first, set the alpha integrator variable */
101 if ((ir->opts.nrdf[0] > 0) && bPscal)
103 vars->alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
109 g = 0.5*veta*ir->delta_t;
110 vars->rscale = exp(g)*series_sinhx(g);
111 g = -0.25*vars->alpha*veta*ir->delta_t;
112 vars->vscale = exp(g)*series_sinhx(g);
113 vars->rvscale = vars->vscale*vars->rscale;
114 vars->veta = vetanew;
118 snew(vars->vscale_nhc, ir->opts.ngtc);
119 if ((ekind == NULL) || (!bPscal))
121 for (i = 0; i < ir->opts.ngtc; i++)
123 vars->vscale_nhc[i] = 1;
128 for (i = 0; i < ir->opts.ngtc; i++)
130 vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
136 vars->vscale_nhc = NULL;
142 static void free_vetavars(t_vetavars *vars)
144 if (vars->vscale_nhc != NULL)
146 sfree(vars->vscale_nhc);
150 static int pcomp(const void *p1, const void *p2)
153 atom_id min1, min2, max1, max2;
154 t_sortblock *a1 = (t_sortblock *)p1;
155 t_sortblock *a2 = (t_sortblock *)p2;
157 db = a1->blocknr-a2->blocknr;
164 min1 = min(a1->iatom[1], a1->iatom[2]);
165 max1 = max(a1->iatom[1], a1->iatom[2]);
166 min2 = min(a2->iatom[1], a2->iatom[2]);
167 max2 = max(a2->iatom[1], a2->iatom[2]);
179 int n_flexible_constraints(struct gmx_constr *constr)
185 nflexcon = constr->nflexcon;
195 void too_many_constraint_warnings(int eConstrAlg, int warncount)
197 const char *abort = "- aborting to avoid logfile runaway.\n"
198 "This normally happens when your system is not sufficiently equilibrated,"
199 "or if you are changing lambda too fast in free energy simulations.\n";
202 "Too many %s warnings (%d)\n"
203 "If you know what you are doing you can %s"
204 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
205 "but normally it is better to fix the problem",
206 (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE", warncount,
207 (eConstrAlg == econtLINCS) ?
208 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
211 static void write_constr_pdb(const char *fn, const char *title,
213 int start, int homenr, t_commrec *cr,
214 rvec x[], matrix box)
216 char fname[STRLEN], format[STRLEN];
218 int dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
225 sprintf(fname, "%s_n%d.pdb", fn, cr->sim_nodeid);
226 if (DOMAINDECOMP(cr))
229 dd_get_constraint_range(dd, &dd_ac0, &dd_ac1);
236 sprintf(fname, "%s.pdb", fn);
238 sprintf(format, "%s\n", get_pdbformat());
240 out = gmx_fio_fopen(fname, "w");
242 fprintf(out, "TITLE %s\n", title);
243 gmx_write_pdb_box(out, -1, box);
244 for (i = start; i < start+homenr; i++)
248 if (i >= dd->nat_home && i < dd_ac0)
252 ii = dd->gatindex[i];
258 gmx_mtop_atominfo_global(mtop, ii, &anm, &resnr, &resnm);
259 fprintf(out, format, "ATOM", (ii+1)%100000,
260 anm, resnm, ' ', resnr%10000, ' ',
261 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ]);
263 fprintf(out, "TER\n");
268 static void dump_confs(FILE *fplog, gmx_int64_t step, gmx_mtop_t *mtop,
269 int start, int homenr, t_commrec *cr,
270 rvec x[], rvec xprime[], matrix box)
272 char buf[256], buf2[22];
274 char *env = getenv("GMX_SUPPRESS_DUMP");
280 sprintf(buf, "step%sb", gmx_step_str(step, buf2));
281 write_constr_pdb(buf, "initial coordinates",
282 mtop, start, homenr, cr, x, box);
283 sprintf(buf, "step%sc", gmx_step_str(step, buf2));
284 write_constr_pdb(buf, "coordinates after constraining",
285 mtop, start, homenr, cr, xprime, box);
288 fprintf(fplog, "Wrote pdb files with previous and current coordinates\n");
290 fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
293 static void pr_sortblock(FILE *fp, const char *title, int nsb, t_sortblock sb[])
297 fprintf(fp, "%s\n", title);
298 for (i = 0; (i < nsb); i++)
300 fprintf(fp, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
301 i, sb[i].iatom[0], sb[i].iatom[1], sb[i].iatom[2],
306 gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
307 struct gmx_constr *constr,
308 t_idef *idef, t_inputrec *ir, gmx_ekindata_t *ekind,
310 gmx_int64_t step, int delta_step,
312 rvec *x, rvec *xprime, rvec *min_proj,
313 gmx_bool bMolPBC, matrix box,
314 real lambda, real *dvdlambda,
315 rvec *v, tensor *vir,
316 t_nrnb *nrnb, int econq, gmx_bool bPscal,
317 real veta, real vetanew)
320 int start, homenr, nrend;
322 int ncons, settle_error;
325 real invdt, vir_fac, t;
328 t_pbc pbc, *pbc_null;
333 if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
335 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");
343 nrend = start+homenr;
345 /* set constants for pressure control integration */
346 init_vetavars(&vetavar, econq != econqCoord,
347 veta, vetanew, ir, ekind, bPscal);
349 if (ir->delta_t == 0)
355 invdt = 1/ir->delta_t;
358 if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
360 /* Set the constraint lengths for the step at which this configuration
361 * is meant to be. The invmasses should not be changed.
363 lambda += delta_step*ir->fepvals->delta_lambda;
368 clear_mat(vir_r_m_dr);
373 settle = &idef->il[F_SETTLE];
374 nsettle = settle->nr/(1+NRAL(F_SETTLE));
378 nth = gmx_omp_nthreads_get(emntSETTLE);
385 if (nth > 1 && constr->vir_r_m_dr_th == NULL)
387 snew(constr->vir_r_m_dr_th, nth);
388 snew(constr->settle_error, nth);
393 /* We do not need full pbc when constraints do not cross charge groups,
394 * i.e. when dd->constraint_comm==NULL.
395 * Note that PBC for constraints is different from PBC for bondeds.
396 * For constraints there is both forward and backward communication.
398 if (ir->ePBC != epbcNONE &&
399 (cr->dd || bMolPBC) && !(cr->dd && cr->dd->constraint_comm == NULL))
401 /* With pbc=screw the screw has been changed to a shift
402 * by the constraint coordinate communication routine,
403 * so that here we can use normal pbc.
405 pbc_null = set_pbc_dd(&pbc, ir->ePBC, cr->dd, FALSE, box);
412 /* Communicate the coordinates required for the non-local constraints
413 * for LINCS and/or SETTLE.
417 dd_move_x_constraints(cr->dd, box, x, xprime);
419 else if (PARTDECOMP(cr))
421 pd_move_x_constraints(cr, x, xprime);
424 if (constr->lincsd != NULL)
426 bOK = constrain_lincs(fplog, bLog, bEner, ir, step, constr->lincsd, md, cr,
428 box, pbc_null, lambda, dvdlambda,
429 invdt, v, vir != NULL, vir_r_m_dr,
431 constr->maxwarn, &constr->warncount_lincs);
432 if (!bOK && constr->maxwarn >= 0)
436 fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
437 econstr_names[econtLINCS], gmx_step_str(step, buf));
443 if (constr->nblocks > 0)
448 bOK = bshakef(fplog, constr->shaked,
449 md->invmass, constr->nblocks, constr->sblock,
450 idef, ir, x, xprime, nrnb,
451 constr->lagr, lambda, dvdlambda,
452 invdt, v, vir != NULL, vir_r_m_dr,
453 constr->maxwarn >= 0, econq, &vetavar);
456 bOK = bshakef(fplog, constr->shaked,
457 md->invmass, constr->nblocks, constr->sblock,
458 idef, ir, x, min_proj, nrnb,
459 constr->lagr, lambda, dvdlambda,
460 invdt, NULL, vir != NULL, vir_r_m_dr,
461 constr->maxwarn >= 0, econq, &vetavar);
464 gmx_fatal(FARGS, "Internal error, SHAKE called for constraining something else than coordinates");
468 if (!bOK && constr->maxwarn >= 0)
472 fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
473 econstr_names[econtSHAKE], gmx_step_str(step, buf));
481 int calcvir_atom_end;
485 calcvir_atom_end = 0;
489 calcvir_atom_end = md->start + md->homenr;
495 #pragma omp parallel for num_threads(nth) schedule(static)
496 for (th = 0; th < nth; th++)
498 int start_th, end_th;
502 clear_mat(constr->vir_r_m_dr_th[th]);
505 start_th = (nsettle* th )/nth;
506 end_th = (nsettle*(th+1))/nth;
507 if (start_th >= 0 && end_th - start_th > 0)
509 csettle(constr->settled,
511 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
514 invdt, v ? v[0] : NULL, calcvir_atom_end,
515 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
516 th == 0 ? &settle_error : &constr->settle_error[th],
520 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
523 inc_nrnb(nrnb, eNR_CONSTR_V, nsettle*3);
527 inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle*3);
533 case econqForceDispl:
534 #pragma omp parallel for num_threads(nth) schedule(static)
535 for (th = 0; th < nth; th++)
537 int start_th, end_th;
541 clear_mat(constr->vir_r_m_dr_th[th]);
544 start_th = (nsettle* th )/nth;
545 end_th = (nsettle*(th+1))/nth;
547 if (start_th >= 0 && end_th - start_th > 0)
549 settle_proj(constr->settled, econq,
551 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
554 xprime, min_proj, calcvir_atom_end,
555 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
559 /* This is an overestimate */
560 inc_nrnb(nrnb, eNR_SETTLE, nsettle);
562 case econqDeriv_FlexCon:
563 /* Nothing to do, since the are no flexible constraints in settles */
566 gmx_incons("Unknown constraint quantity for settle");
572 /* Combine virial and error info of the other threads */
573 for (i = 1; i < nth; i++)
575 m_add(vir_r_m_dr, constr->vir_r_m_dr_th[i], vir_r_m_dr);
576 settle_error = constr->settle_error[i];
579 if (econq == econqCoord && settle_error >= 0)
582 if (constr->maxwarn >= 0)
586 "\nstep " "%"GMX_PRId64 ": Water molecule starting at atom %d can not be "
587 "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
588 step, ddglatnr(cr->dd, settle->iatoms[settle_error*(1+NRAL(F_SETTLE))+1]));
591 fprintf(fplog, "%s", buf);
593 fprintf(stderr, "%s", buf);
594 constr->warncount_settle++;
595 if (constr->warncount_settle > constr->maxwarn)
597 too_many_constraint_warnings(-1, constr->warncount_settle);
604 free_vetavars(&vetavar);
611 vir_fac = 0.5/(ir->delta_t*ir->delta_t);
614 vir_fac = 0.5/ir->delta_t;
617 case econqForceDispl:
622 gmx_incons("Unsupported constraint quantity for virial");
627 vir_fac *= 2; /* only constraining over half the distance here */
629 for (i = 0; i < DIM; i++)
631 for (j = 0; j < DIM; j++)
633 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
640 dump_confs(fplog, step, constr->warn_mtop, start, homenr, cr, x, xprime, box);
643 if (econq == econqCoord)
645 if (ir->ePull == epullCONSTRAINT)
647 if (EI_DYNAMICS(ir->eI))
649 t = ir->init_t + (step + delta_step)*ir->delta_t;
655 set_pbc(&pbc, ir->ePBC, box);
656 pull_constraint(ir->pull, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir);
658 if (constr->ed && delta_step > 0)
660 /* apply the essential dynamcs constraints here */
661 do_edsam(ir, step, cr, xprime, v, box, constr->ed);
668 real *constr_rmsd_data(struct gmx_constr *constr)
672 return lincs_rmsd_data(constr->lincsd);
680 real constr_rmsd(struct gmx_constr *constr, gmx_bool bSD2)
684 return lincs_rmsd(constr->lincsd, bSD2);
692 static void make_shake_sblock_pd(struct gmx_constr *constr,
693 t_idef *idef, t_mdatoms *md)
702 /* Since we are processing the local topology,
703 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
705 ncons = idef->il[F_CONSTR].nr/3;
707 init_blocka(&sblocks);
708 gen_sblocks(NULL, md->start, md->start+md->homenr, idef, &sblocks, FALSE);
711 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
712 nblocks=blocks->multinr[idef->nodeid] - bstart;
715 constr->nblocks = sblocks.nr;
718 fprintf(debug, "ncons: %d, bstart: %d, nblocks: %d\n",
719 ncons, bstart, constr->nblocks);
722 /* Calculate block number for each atom */
723 inv_sblock = make_invblocka(&sblocks, md->nr);
725 done_blocka(&sblocks);
727 /* Store the block number in temp array and
728 * sort the constraints in order of the sblock number
729 * and the atom numbers, really sorting a segment of the array!
732 pr_idef(fplog, 0, "Before Sort", idef);
734 iatom = idef->il[F_CONSTR].iatoms;
736 for (i = 0; (i < ncons); i++, iatom += 3)
738 for (m = 0; (m < 3); m++)
740 sb[i].iatom[m] = iatom[m];
742 sb[i].blocknr = inv_sblock[iatom[1]];
745 /* Now sort the blocks */
748 pr_sortblock(debug, "Before sorting", ncons, sb);
749 fprintf(debug, "Going to sort constraints\n");
752 qsort(sb, ncons, (size_t)sizeof(*sb), pcomp);
756 pr_sortblock(debug, "After sorting", ncons, sb);
759 iatom = idef->il[F_CONSTR].iatoms;
760 for (i = 0; (i < ncons); i++, iatom += 3)
762 for (m = 0; (m < 3); m++)
764 iatom[m] = sb[i].iatom[m];
768 pr_idef(fplog, 0, "After Sort", idef);
772 snew(constr->sblock, constr->nblocks+1);
774 for (i = 0; (i < ncons); i++)
776 if (sb[i].blocknr != bnr)
779 constr->sblock[j++] = 3*i;
783 constr->sblock[j++] = 3*ncons;
785 if (j != (constr->nblocks+1))
787 fprintf(stderr, "bstart: %d\n", bstart);
788 fprintf(stderr, "j: %d, nblocks: %d, ncons: %d\n",
789 j, constr->nblocks, ncons);
790 for (i = 0; (i < ncons); i++)
792 fprintf(stderr, "i: %5d sb[i].blocknr: %5u\n", i, sb[i].blocknr);
794 for (j = 0; (j <= constr->nblocks); j++)
796 fprintf(stderr, "sblock[%3d]=%5d\n", j, (int)constr->sblock[j]);
798 gmx_fatal(FARGS, "DEATH HORROR: "
799 "sblocks does not match idef->il[F_CONSTR]");
805 static void make_shake_sblock_dd(struct gmx_constr *constr,
806 t_ilist *ilcon, t_block *cgs,
812 if (dd->ncg_home+1 > constr->sblock_nalloc)
814 constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
815 srenew(constr->sblock, constr->sblock_nalloc);
819 iatom = ilcon->iatoms;
822 for (c = 0; c < ncons; c++)
824 if (c == 0 || iatom[1] >= cgs->index[cg+1])
826 constr->sblock[constr->nblocks++] = 3*c;
827 while (iatom[1] >= cgs->index[cg+1])
834 constr->sblock[constr->nblocks] = 3*ncons;
837 t_blocka make_at2con(int start, int natoms,
838 t_ilist *ilist, t_iparams *iparams,
839 gmx_bool bDynamics, int *nflexiblecons)
841 int *count, ncon, con, con_tot, nflexcon, ftype, i, a;
848 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
850 ncon = ilist[ftype].nr/3;
851 ia = ilist[ftype].iatoms;
852 for (con = 0; con < ncon; con++)
854 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
855 iparams[ia[0]].constr.dB == 0);
860 if (bDynamics || !bFlexCon)
862 for (i = 1; i < 3; i++)
871 *nflexiblecons = nflexcon;
874 at2con.nalloc_index = at2con.nr+1;
875 snew(at2con.index, at2con.nalloc_index);
877 for (a = 0; a < natoms; a++)
879 at2con.index[a+1] = at2con.index[a] + count[a];
882 at2con.nra = at2con.index[natoms];
883 at2con.nalloc_a = at2con.nra;
884 snew(at2con.a, at2con.nalloc_a);
886 /* The F_CONSTRNC constraints have constraint numbers
887 * that continue after the last F_CONSTR constraint.
890 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
892 ncon = ilist[ftype].nr/3;
893 ia = ilist[ftype].iatoms;
894 for (con = 0; con < ncon; con++)
896 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
897 iparams[ia[0]].constr.dB == 0);
898 if (bDynamics || !bFlexCon)
900 for (i = 1; i < 3; i++)
903 at2con.a[at2con.index[a]+count[a]++] = con_tot;
916 static int *make_at2settle(int natoms, const t_ilist *ilist)
922 /* Set all to no settle */
923 for (a = 0; a < natoms; a++)
928 stride = 1 + NRAL(F_SETTLE);
930 for (s = 0; s < ilist->nr; s += stride)
932 at2s[ilist->iatoms[s+1]] = s/stride;
933 at2s[ilist->iatoms[s+2]] = s/stride;
934 at2s[ilist->iatoms[s+3]] = s/stride;
940 void set_constraints(struct gmx_constr *constr,
941 gmx_localtop_t *top, t_inputrec *ir,
942 t_mdatoms *md, t_commrec *cr)
951 if (constr->ncon_tot > 0)
953 /* We are using the local topology,
954 * so there are only F_CONSTR constraints.
956 ncons = idef->il[F_CONSTR].nr/3;
958 /* With DD we might also need to call LINCS with ncons=0 for
959 * communicating coordinates to other nodes that do have constraints.
961 if (ir->eConstrAlg == econtLINCS)
963 set_lincs(idef, md, EI_DYNAMICS(ir->eI), cr, constr->lincsd);
965 if (ir->eConstrAlg == econtSHAKE)
969 make_shake_sblock_dd(constr, &idef->il[F_CONSTR], &top->cgs, cr->dd);
973 make_shake_sblock_pd(constr, idef, md);
975 if (ncons > constr->lagr_nalloc)
977 constr->lagr_nalloc = over_alloc_dd(ncons);
978 srenew(constr->lagr, constr->lagr_nalloc);
983 if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
985 settle = &idef->il[F_SETTLE];
986 iO = settle->iatoms[1];
987 iH = settle->iatoms[2];
989 settle_init(md->massT[iO], md->massT[iH],
990 md->invmass[iO], md->invmass[iH],
991 idef->iparams[settle->iatoms[0]].settle.doh,
992 idef->iparams[settle->iatoms[0]].settle.dhh);
995 /* Make a selection of the local atoms for essential dynamics */
996 if (constr->ed && cr->dd)
998 dd_make_local_ed_indices(cr->dd, constr->ed);
1002 static void constr_recur(t_blocka *at2con,
1003 t_ilist *ilist, t_iparams *iparams, gmx_bool bTopB,
1004 int at, int depth, int nc, int *path,
1005 real r0, real r1, real *r2max,
1017 ncon1 = ilist[F_CONSTR].nr/3;
1018 ia1 = ilist[F_CONSTR].iatoms;
1019 ia2 = ilist[F_CONSTRNC].iatoms;
1021 /* Loop over all constraints connected to this atom */
1022 for (c = at2con->index[at]; c < at2con->index[at+1]; c++)
1025 /* Do not walk over already used constraints */
1027 for (a1 = 0; a1 < depth; a1++)
1029 if (con == path[a1])
1036 ia = constr_iatomptr(ncon1, ia1, ia2, con);
1037 /* Flexible constraints currently have length 0, which is incorrect */
1040 len = iparams[ia[0]].constr.dA;
1044 len = iparams[ia[0]].constr.dB;
1046 /* In the worst case the bond directions alternate */
1057 /* Assume angles of 120 degrees between all bonds */
1058 if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max)
1060 *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
1063 fprintf(debug, "Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0, rn1, sqrt(*r2max));
1064 for (a1 = 0; a1 < depth; a1++)
1066 fprintf(debug, " %d %5.3f",
1068 iparams[constr_iatomptr(ncon1, ia1, ia2, con)[0]].constr.dA);
1070 fprintf(debug, " %d %5.3f\n", con, len);
1073 /* Limit the number of recursions to 1000*nc,
1074 * so a call does not take more than a second,
1075 * even for highly connected systems.
1077 if (depth + 1 < nc && *count < 1000*nc)
1089 constr_recur(at2con, ilist, iparams,
1090 bTopB, a1, depth+1, nc, path, rn0, rn1, r2max, count);
1097 static real constr_r_max_moltype(gmx_moltype_t *molt, t_iparams *iparams,
1100 int natoms, nflexcon, *path, at, count;
1103 real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
1105 if (molt->ilist[F_CONSTR].nr == 0 &&
1106 molt->ilist[F_CONSTRNC].nr == 0)
1111 natoms = molt->atoms.nr;
1113 at2con = make_at2con(0, natoms, molt->ilist, iparams,
1114 EI_DYNAMICS(ir->eI), &nflexcon);
1115 snew(path, 1+ir->nProjOrder);
1116 for (at = 0; at < 1+ir->nProjOrder; at++)
1122 for (at = 0; at < natoms; at++)
1128 constr_recur(&at2con, molt->ilist, iparams,
1129 FALSE, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxA, &count);
1131 if (ir->efep == efepNO)
1133 rmax = sqrt(r2maxA);
1138 for (at = 0; at < natoms; at++)
1143 constr_recur(&at2con, molt->ilist, iparams,
1144 TRUE, at, 0, 1+ir->nProjOrder, path, r0, r1, &r2maxB, &count);
1146 lam0 = ir->fepvals->init_lambda;
1147 if (EI_DYNAMICS(ir->eI))
1149 lam0 += ir->init_step*ir->fepvals->delta_lambda;
1151 rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
1152 if (EI_DYNAMICS(ir->eI))
1154 lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
1155 rmax = max(rmax, (1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
1159 done_blocka(&at2con);
1165 real constr_r_max(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir)
1171 for (mt = 0; mt < mtop->nmoltype; mt++)
1174 constr_r_max_moltype(&mtop->moltype[mt],
1175 mtop->ffparams.iparams, ir));
1180 fprintf(fplog, "Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n", 1+ir->nProjOrder, rmax);
1186 gmx_constr_t init_constraints(FILE *fplog,
1187 gmx_mtop_t *mtop, t_inputrec *ir,
1188 gmx_edsam_t ed, t_state *state,
1191 int ncon, nset, nmol, settle_type, i, natoms, mt, nflexcon;
1192 struct gmx_constr *constr;
1195 gmx_mtop_ilistloop_t iloop;
1198 gmx_mtop_ftype_count(mtop, F_CONSTR) +
1199 gmx_mtop_ftype_count(mtop, F_CONSTRNC);
1200 nset = gmx_mtop_ftype_count(mtop, F_SETTLE);
1202 if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
1209 constr->ncon_tot = ncon;
1210 constr->nflexcon = 0;
1213 constr->n_at2con_mt = mtop->nmoltype;
1214 snew(constr->at2con_mt, constr->n_at2con_mt);
1215 for (mt = 0; mt < mtop->nmoltype; mt++)
1217 constr->at2con_mt[mt] = make_at2con(0, mtop->moltype[mt].atoms.nr,
1218 mtop->moltype[mt].ilist,
1219 mtop->ffparams.iparams,
1220 EI_DYNAMICS(ir->eI), &nflexcon);
1221 for (i = 0; i < mtop->nmolblock; i++)
1223 if (mtop->molblock[i].type == mt)
1225 constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
1230 if (constr->nflexcon > 0)
1234 fprintf(fplog, "There are %d flexible constraints\n",
1236 if (ir->fc_stepsize == 0)
1239 "WARNING: step size for flexible constraining = 0\n"
1240 " All flexible constraints will be rigid.\n"
1241 " Will try to keep all flexible constraints at their original length,\n"
1242 " but the lengths may exhibit some drift.\n\n");
1243 constr->nflexcon = 0;
1246 if (constr->nflexcon > 0)
1248 please_cite(fplog, "Hess2002");
1252 if (ir->eConstrAlg == econtLINCS)
1254 constr->lincsd = init_lincs(fplog, mtop,
1255 constr->nflexcon, constr->at2con_mt,
1256 DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1257 ir->nLincsIter, ir->nProjOrder);
1260 if (ir->eConstrAlg == econtSHAKE)
1262 if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1264 gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1266 if (constr->nflexcon)
1268 gmx_fatal(FARGS, "For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1270 please_cite(fplog, "Ryckaert77a");
1273 please_cite(fplog, "Barth95a");
1276 constr->shaked = shake_init();
1282 please_cite(fplog, "Miyamoto92a");
1284 constr->bInterCGsettles = inter_charge_group_settles(mtop);
1286 /* Check that we have only one settle type */
1288 iloop = gmx_mtop_ilistloop_init(mtop);
1289 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
1291 for (i = 0; i < ilist[F_SETTLE].nr; i += 4)
1293 if (settle_type == -1)
1295 settle_type = ilist[F_SETTLE].iatoms[i];
1297 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1300 "The [molecules] section of your topology specifies more than one block of\n"
1301 "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1302 "are trying to partition your solvent into different *groups* (e.g. for\n"
1303 "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1304 "files specify groups. Otherwise, you may wish to change the least-used\n"
1305 "block of molecules with SETTLE constraints into 3 normal constraints.");
1310 constr->n_at2settle_mt = mtop->nmoltype;
1311 snew(constr->at2settle_mt, constr->n_at2settle_mt);
1312 for (mt = 0; mt < mtop->nmoltype; mt++)
1314 constr->at2settle_mt[mt] =
1315 make_at2settle(mtop->moltype[mt].atoms.nr,
1316 &mtop->moltype[mt].ilist[F_SETTLE]);
1320 constr->maxwarn = 999;
1321 env = getenv("GMX_MAXCONSTRWARN");
1324 constr->maxwarn = 0;
1325 sscanf(env, "%d", &constr->maxwarn);
1329 "Setting the maximum number of constraint warnings to %d\n",
1335 "Setting the maximum number of constraint warnings to %d\n",
1339 if (constr->maxwarn < 0 && fplog)
1341 fprintf(fplog, "maxwarn < 0, will not stop on constraint errors\n");
1343 constr->warncount_lincs = 0;
1344 constr->warncount_settle = 0;
1346 /* Initialize the essential dynamics sampling.
1347 * Put the pointer to the ED struct in constr */
1349 if (ed != NULL || state->edsamstate.nED > 0)
1351 init_edsam(mtop, ir, cr, ed, state->x, state->box, &state->edsamstate);
1354 constr->warn_mtop = mtop;
1359 const t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1361 return constr->at2con_mt;
1364 const int **atom2settle_moltype(gmx_constr_t constr)
1366 return (const int **)constr->at2settle_mt;
1370 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
1372 const gmx_moltype_t *molt;
1376 int nat, *at2cg, cg, a, ftype, i;
1380 for (mb = 0; mb < mtop->nmolblock && !bInterCG; mb++)
1382 molt = &mtop->moltype[mtop->molblock[mb].type];
1384 if (molt->ilist[F_CONSTR].nr > 0 ||
1385 molt->ilist[F_CONSTRNC].nr > 0 ||
1386 molt->ilist[F_SETTLE].nr > 0)
1389 snew(at2cg, molt->atoms.nr);
1390 for (cg = 0; cg < cgs->nr; cg++)
1392 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1398 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
1400 il = &molt->ilist[ftype];
1401 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype))
1403 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1417 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop)
1419 const gmx_moltype_t *molt;
1423 int nat, *at2cg, cg, a, ftype, i;
1427 for (mb = 0; mb < mtop->nmolblock && !bInterCG; mb++)
1429 molt = &mtop->moltype[mtop->molblock[mb].type];
1431 if (molt->ilist[F_SETTLE].nr > 0)
1434 snew(at2cg, molt->atoms.nr);
1435 for (cg = 0; cg < cgs->nr; cg++)
1437 for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
1443 for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
1445 il = &molt->ilist[ftype];
1446 for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE))
1448 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1449 at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1463 /* helper functions for andersen temperature control, because the
1464 * gmx_constr construct is only defined in constr.c. Return the list
1465 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
1467 extern int *get_sblock(struct gmx_constr *constr)
1469 return constr->sblock;
1472 extern int get_nblocks(struct gmx_constr *constr)
1474 return constr->nblocks;