2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
42 #include "gromacs/fileio/confio.h"
43 #include "types/commrec.h"
48 #include "gromacs/math/vec.h"
52 #include "gromacs/fileio/pdbio.h"
54 #include "gromacs/topology/mtop_util.h"
55 #include "gromacs/fileio/gmxfio.h"
57 #include "gmx_omp_nthreads.h"
58 #include "gromacs/essentialdynamics/edsam.h"
59 #include "gromacs/pulling/pull.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/topology/block.h"
63 #include "gromacs/topology/invblock.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/smalloc.h"
67 typedef 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 */
86 gmx_edsam_t ed
; /* The essential dynamics data */
88 tensor
*vir_r_m_dr_th
; /* Thread local working data */
89 int *settle_error
; /* Thread local working data */
91 gmx_mtop_t
*warn_mtop
; /* Only used for printing warnings */
99 /* delta_t should be used instead of ir->delta_t, to permit the time
100 step to be scaled by the calling code */
101 static void *init_vetavars(t_vetavars
*vars
,
102 gmx_bool constr_deriv
,
103 real veta
, real vetanew
,
104 t_inputrec
*ir
, real delta_t
,
105 gmx_ekindata_t
*ekind
, gmx_bool bPscal
)
110 /* first, set the alpha integrator variable */
111 if ((ir
->opts
.nrdf
[0] > 0) && bPscal
)
113 vars
->alpha
= 1.0 + DIM
/((double)ir
->opts
.nrdf
[0]);
119 g
= 0.5*veta
*delta_t
;
120 vars
->rscale
= exp(g
)*series_sinhx(g
);
121 g
= -0.25*vars
->alpha
*veta
*delta_t
;
122 vars
->vscale
= exp(g
)*series_sinhx(g
);
123 vars
->rvscale
= vars
->vscale
*vars
->rscale
;
124 vars
->veta
= vetanew
;
128 snew(vars
->vscale_nhc
, ir
->opts
.ngtc
);
129 if ((ekind
== NULL
) || (!bPscal
))
131 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
133 vars
->vscale_nhc
[i
] = 1;
138 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
140 vars
->vscale_nhc
[i
] = ekind
->tcstat
[i
].vscale_nhc
;
146 vars
->vscale_nhc
= NULL
;
152 static void free_vetavars(t_vetavars
*vars
)
154 if (vars
->vscale_nhc
!= NULL
)
156 sfree(vars
->vscale_nhc
);
160 static int pcomp(const void *p1
, const void *p2
)
163 atom_id min1
, min2
, max1
, max2
;
164 t_sortblock
*a1
= (t_sortblock
*)p1
;
165 t_sortblock
*a2
= (t_sortblock
*)p2
;
167 db
= a1
->blocknr
-a2
->blocknr
;
174 min1
= min(a1
->iatom
[1], a1
->iatom
[2]);
175 max1
= max(a1
->iatom
[1], a1
->iatom
[2]);
176 min2
= min(a2
->iatom
[1], a2
->iatom
[2]);
177 max2
= max(a2
->iatom
[1], a2
->iatom
[2]);
189 int n_flexible_constraints(struct gmx_constr
*constr
)
195 nflexcon
= constr
->nflexcon
;
205 void too_many_constraint_warnings(int eConstrAlg
, int warncount
)
207 const char *abort
= "- aborting to avoid logfile runaway.\n"
208 "This normally happens when your system is not sufficiently equilibrated,"
209 "or if you are changing lambda too fast in free energy simulations.\n";
212 "Too many %s warnings (%d)\n"
213 "If you know what you are doing you can %s"
214 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
215 "but normally it is better to fix the problem",
216 (eConstrAlg
== econtLINCS
) ? "LINCS" : "SETTLE", warncount
,
217 (eConstrAlg
== econtLINCS
) ?
218 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
221 static void write_constr_pdb(const char *fn
, const char *title
,
223 int start
, int homenr
, t_commrec
*cr
,
224 rvec x
[], matrix box
)
228 int dd_ac0
= 0, dd_ac1
= 0, i
, ii
, resnr
;
233 if (DOMAINDECOMP(cr
))
236 dd_get_constraint_range(dd
, &dd_ac0
, &dd_ac1
);
243 sprintf(fname
, "%s_n%d.pdb", fn
, cr
->sim_nodeid
);
247 sprintf(fname
, "%s.pdb", fn
);
250 out
= gmx_fio_fopen(fname
, "w");
252 fprintf(out
, "TITLE %s\n", title
);
253 gmx_write_pdb_box(out
, -1, box
);
254 for (i
= start
; i
< start
+homenr
; i
++)
258 if (i
>= dd
->nat_home
&& i
< dd_ac0
)
262 ii
= dd
->gatindex
[i
];
268 gmx_mtop_atominfo_global(mtop
, ii
, &anm
, &resnr
, &resnm
);
269 gmx_fprintf_pdb_atomline(out
, epdbATOM
, ii
+1, anm
, ' ', resnm
, ' ', resnr
, ' ',
270 10*x
[i
][XX
], 10*x
[i
][YY
], 10*x
[i
][ZZ
], 1.0, 0.0, "");
272 fprintf(out
, "TER\n");
277 static void dump_confs(FILE *fplog
, gmx_int64_t step
, gmx_mtop_t
*mtop
,
278 int start
, int homenr
, t_commrec
*cr
,
279 rvec x
[], rvec xprime
[], matrix box
)
281 char buf
[256], buf2
[22];
283 char *env
= getenv("GMX_SUPPRESS_DUMP");
289 sprintf(buf
, "step%sb", gmx_step_str(step
, buf2
));
290 write_constr_pdb(buf
, "initial coordinates",
291 mtop
, start
, homenr
, cr
, x
, box
);
292 sprintf(buf
, "step%sc", gmx_step_str(step
, buf2
));
293 write_constr_pdb(buf
, "coordinates after constraining",
294 mtop
, start
, homenr
, cr
, xprime
, box
);
297 fprintf(fplog
, "Wrote pdb files with previous and current coordinates\n");
299 fprintf(stderr
, "Wrote pdb files with previous and current coordinates\n");
302 static void pr_sortblock(FILE *fp
, const char *title
, int nsb
, t_sortblock sb
[])
306 fprintf(fp
, "%s\n", title
);
307 for (i
= 0; (i
< nsb
); i
++)
309 fprintf(fp
, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
310 i
, sb
[i
].iatom
[0], sb
[i
].iatom
[1], sb
[i
].iatom
[2],
315 gmx_bool
constrain(FILE *fplog
, gmx_bool bLog
, gmx_bool bEner
,
316 struct gmx_constr
*constr
,
317 t_idef
*idef
, t_inputrec
*ir
, gmx_ekindata_t
*ekind
,
319 gmx_int64_t step
, int delta_step
,
322 rvec
*x
, rvec
*xprime
, rvec
*min_proj
,
323 gmx_bool bMolPBC
, matrix box
,
324 real lambda
, real
*dvdlambda
,
325 rvec
*v
, tensor
*vir
,
326 t_nrnb
*nrnb
, int econq
, gmx_bool bPscal
,
327 real veta
, real vetanew
)
330 int start
, homenr
, nrend
;
332 int ncons
, settle_error
;
336 real invdt
, vir_fac
, t
;
339 t_pbc pbc
, *pbc_null
;
344 if (econq
== econqForceDispl
&& !EI_ENERGY_MINIMIZATION(ir
->eI
))
346 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");
354 nrend
= start
+homenr
;
356 scaled_delta_t
= step_scaling
* ir
->delta_t
;
358 /* set constants for pressure control integration */
359 init_vetavars(&vetavar
, econq
!= econqCoord
,
360 veta
, vetanew
, ir
, scaled_delta_t
, ekind
, bPscal
);
362 /* Prepare time step for use in constraint implementations, and
363 avoid generating inf when ir->delta_t = 0. */
364 if (ir
->delta_t
== 0)
370 invdt
= 1.0/scaled_delta_t
;
373 if (ir
->efep
!= efepNO
&& EI_DYNAMICS(ir
->eI
))
375 /* Set the constraint lengths for the step at which this configuration
376 * is meant to be. The invmasses should not be changed.
378 lambda
+= delta_step
*ir
->fepvals
->delta_lambda
;
383 clear_mat(vir_r_m_dr
);
388 settle
= &idef
->il
[F_SETTLE
];
389 nsettle
= settle
->nr
/(1+NRAL(F_SETTLE
));
393 nth
= gmx_omp_nthreads_get(emntSETTLE
);
400 if (nth
> 1 && constr
->vir_r_m_dr_th
== NULL
)
402 snew(constr
->vir_r_m_dr_th
, nth
);
403 snew(constr
->settle_error
, nth
);
408 /* We do not need full pbc when constraints do not cross charge groups,
409 * i.e. when dd->constraint_comm==NULL.
410 * Note that PBC for constraints is different from PBC for bondeds.
411 * For constraints there is both forward and backward communication.
413 if (ir
->ePBC
!= epbcNONE
&&
414 (cr
->dd
|| bMolPBC
) && !(cr
->dd
&& cr
->dd
->constraint_comm
== NULL
))
416 /* With pbc=screw the screw has been changed to a shift
417 * by the constraint coordinate communication routine,
418 * so that here we can use normal pbc.
420 pbc_null
= set_pbc_dd(&pbc
, ir
->ePBC
, cr
->dd
, FALSE
, box
);
427 /* Communicate the coordinates required for the non-local constraints
428 * for LINCS and/or SETTLE.
432 dd_move_x_constraints(cr
->dd
, box
, x
, xprime
, econq
== econqCoord
);
435 if (constr
->lincsd
!= NULL
)
437 bOK
= constrain_lincs(fplog
, bLog
, bEner
, ir
, step
, constr
->lincsd
, md
, cr
,
439 box
, pbc_null
, lambda
, dvdlambda
,
440 invdt
, v
, vir
!= NULL
, vir_r_m_dr
,
442 constr
->maxwarn
, &constr
->warncount_lincs
);
443 if (!bOK
&& constr
->maxwarn
>= 0)
447 fprintf(fplog
, "Constraint error in algorithm %s at step %s\n",
448 econstr_names
[econtLINCS
], gmx_step_str(step
, buf
));
454 if (constr
->nblocks
> 0)
459 bOK
= bshakef(fplog
, constr
->shaked
,
460 md
->invmass
, constr
->nblocks
, constr
->sblock
,
461 idef
, ir
, x
, xprime
, nrnb
,
462 constr
->lagr
, lambda
, dvdlambda
,
463 invdt
, v
, vir
!= NULL
, vir_r_m_dr
,
464 constr
->maxwarn
>= 0, econq
, &vetavar
);
467 bOK
= bshakef(fplog
, constr
->shaked
,
468 md
->invmass
, constr
->nblocks
, constr
->sblock
,
469 idef
, ir
, x
, min_proj
, nrnb
,
470 constr
->lagr
, lambda
, dvdlambda
,
471 invdt
, NULL
, vir
!= NULL
, vir_r_m_dr
,
472 constr
->maxwarn
>= 0, econq
, &vetavar
);
475 gmx_fatal(FARGS
, "Internal error, SHAKE called for constraining something else than coordinates");
479 if (!bOK
&& constr
->maxwarn
>= 0)
483 fprintf(fplog
, "Constraint error in algorithm %s at step %s\n",
484 econstr_names
[econtSHAKE
], gmx_step_str(step
, buf
));
492 int calcvir_atom_end
;
496 calcvir_atom_end
= 0;
500 calcvir_atom_end
= md
->homenr
;
506 #pragma omp parallel for num_threads(nth) schedule(static)
507 for (th
= 0; th
< nth
; th
++)
509 int start_th
, end_th
;
513 clear_mat(constr
->vir_r_m_dr_th
[th
]);
516 start_th
= (nsettle
* th
)/nth
;
517 end_th
= (nsettle
*(th
+1))/nth
;
518 if (start_th
>= 0 && end_th
- start_th
> 0)
520 csettle(constr
->settled
,
522 settle
->iatoms
+start_th
*(1+NRAL(F_SETTLE
)),
525 invdt
, v
? v
[0] : NULL
, calcvir_atom_end
,
526 th
== 0 ? vir_r_m_dr
: constr
->vir_r_m_dr_th
[th
],
527 th
== 0 ? &settle_error
: &constr
->settle_error
[th
],
531 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
534 inc_nrnb(nrnb
, eNR_CONSTR_V
, nsettle
*3);
538 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, nsettle
*3);
544 case econqForceDispl
:
545 #pragma omp parallel for num_threads(nth) schedule(static)
546 for (th
= 0; th
< nth
; th
++)
548 int start_th
, end_th
;
552 clear_mat(constr
->vir_r_m_dr_th
[th
]);
555 start_th
= (nsettle
* th
)/nth
;
556 end_th
= (nsettle
*(th
+1))/nth
;
558 if (start_th
>= 0 && end_th
- start_th
> 0)
560 settle_proj(constr
->settled
, econq
,
562 settle
->iatoms
+start_th
*(1+NRAL(F_SETTLE
)),
565 xprime
, min_proj
, calcvir_atom_end
,
566 th
== 0 ? vir_r_m_dr
: constr
->vir_r_m_dr_th
[th
],
570 /* This is an overestimate */
571 inc_nrnb(nrnb
, eNR_SETTLE
, nsettle
);
573 case econqDeriv_FlexCon
:
574 /* Nothing to do, since the are no flexible constraints in settles */
577 gmx_incons("Unknown constraint quantity for settle");
583 /* Combine virial and error info of the other threads */
584 for (i
= 1; i
< nth
; i
++)
586 m_add(vir_r_m_dr
, constr
->vir_r_m_dr_th
[i
], vir_r_m_dr
);
587 settle_error
= constr
->settle_error
[i
];
590 if (econq
== econqCoord
&& settle_error
>= 0)
593 if (constr
->maxwarn
>= 0)
597 "\nstep " "%"GMX_PRId64
": Water molecule starting at atom %d can not be "
598 "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
599 step
, ddglatnr(cr
->dd
, settle
->iatoms
[settle_error
*(1+NRAL(F_SETTLE
))+1]));
602 fprintf(fplog
, "%s", buf
);
604 fprintf(stderr
, "%s", buf
);
605 constr
->warncount_settle
++;
606 if (constr
->warncount_settle
> constr
->maxwarn
)
608 too_many_constraint_warnings(-1, constr
->warncount_settle
);
615 free_vetavars(&vetavar
);
619 /* The normal uses of constrain() pass step_scaling = 1.0.
620 * The call to constrain() for SD1 that passes step_scaling =
621 * 0.5 also passes vir = NULL, so cannot reach this
622 * assertion. This assertion should remain until someone knows
623 * that this path works for their intended purpose, and then
624 * they can use scaled_delta_t instead of ir->delta_t
626 assert(gmx_within_tol(step_scaling
, 1.0, GMX_REAL_EPS
));
630 vir_fac
= 0.5/(ir
->delta_t
*ir
->delta_t
);
633 vir_fac
= 0.5/ir
->delta_t
;
636 case econqForceDispl
:
641 gmx_incons("Unsupported constraint quantity for virial");
646 vir_fac
*= 2; /* only constraining over half the distance here */
648 for (i
= 0; i
< DIM
; i
++)
650 for (j
= 0; j
< DIM
; j
++)
652 (*vir
)[i
][j
] = vir_fac
*vir_r_m_dr
[i
][j
];
659 dump_confs(fplog
, step
, constr
->warn_mtop
, start
, homenr
, cr
, x
, xprime
, box
);
662 if (econq
== econqCoord
)
664 if (ir
->ePull
== epullCONSTRAINT
)
666 if (EI_DYNAMICS(ir
->eI
))
668 t
= ir
->init_t
+ (step
+ delta_step
)*ir
->delta_t
;
674 set_pbc(&pbc
, ir
->ePBC
, box
);
675 pull_constraint(ir
->pull
, md
, &pbc
, cr
, ir
->delta_t
, t
, x
, xprime
, v
, *vir
);
677 if (constr
->ed
&& delta_step
> 0)
679 /* apply the essential dynamcs constraints here */
680 do_edsam(ir
, step
, cr
, xprime
, v
, box
, constr
->ed
);
687 real
*constr_rmsd_data(struct gmx_constr
*constr
)
691 return lincs_rmsd_data(constr
->lincsd
);
699 real
constr_rmsd(struct gmx_constr
*constr
, gmx_bool bSD2
)
703 return lincs_rmsd(constr
->lincsd
, bSD2
);
711 static void make_shake_sblock_serial(struct gmx_constr
*constr
,
712 t_idef
*idef
, t_mdatoms
*md
)
721 /* Since we are processing the local topology,
722 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
724 ncons
= idef
->il
[F_CONSTR
].nr
/3;
726 init_blocka(&sblocks
);
727 gen_sblocks(NULL
, 0, md
->homenr
, idef
, &sblocks
, FALSE
);
730 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
731 nblocks=blocks->multinr[idef->nodeid] - bstart;
734 constr
->nblocks
= sblocks
.nr
;
737 fprintf(debug
, "ncons: %d, bstart: %d, nblocks: %d\n",
738 ncons
, bstart
, constr
->nblocks
);
741 /* Calculate block number for each atom */
742 inv_sblock
= make_invblocka(&sblocks
, md
->nr
);
744 done_blocka(&sblocks
);
746 /* Store the block number in temp array and
747 * sort the constraints in order of the sblock number
748 * and the atom numbers, really sorting a segment of the array!
751 pr_idef(fplog
, 0, "Before Sort", idef
);
753 iatom
= idef
->il
[F_CONSTR
].iatoms
;
755 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
757 for (m
= 0; (m
< 3); m
++)
759 sb
[i
].iatom
[m
] = iatom
[m
];
761 sb
[i
].blocknr
= inv_sblock
[iatom
[1]];
764 /* Now sort the blocks */
767 pr_sortblock(debug
, "Before sorting", ncons
, sb
);
768 fprintf(debug
, "Going to sort constraints\n");
771 qsort(sb
, ncons
, (size_t)sizeof(*sb
), pcomp
);
775 pr_sortblock(debug
, "After sorting", ncons
, sb
);
778 iatom
= idef
->il
[F_CONSTR
].iatoms
;
779 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
781 for (m
= 0; (m
< 3); m
++)
783 iatom
[m
] = sb
[i
].iatom
[m
];
787 pr_idef(fplog
, 0, "After Sort", idef
);
791 snew(constr
->sblock
, constr
->nblocks
+1);
793 for (i
= 0; (i
< ncons
); i
++)
795 if (sb
[i
].blocknr
!= bnr
)
798 constr
->sblock
[j
++] = 3*i
;
802 constr
->sblock
[j
++] = 3*ncons
;
804 if (j
!= (constr
->nblocks
+1))
806 fprintf(stderr
, "bstart: %d\n", bstart
);
807 fprintf(stderr
, "j: %d, nblocks: %d, ncons: %d\n",
808 j
, constr
->nblocks
, ncons
);
809 for (i
= 0; (i
< ncons
); i
++)
811 fprintf(stderr
, "i: %5d sb[i].blocknr: %5d\n", i
, sb
[i
].blocknr
);
813 for (j
= 0; (j
<= constr
->nblocks
); j
++)
815 fprintf(stderr
, "sblock[%3d]=%5d\n", j
, (int)constr
->sblock
[j
]);
817 gmx_fatal(FARGS
, "DEATH HORROR: "
818 "sblocks does not match idef->il[F_CONSTR]");
824 static void make_shake_sblock_dd(struct gmx_constr
*constr
,
825 t_ilist
*ilcon
, t_block
*cgs
,
831 if (dd
->ncg_home
+1 > constr
->sblock_nalloc
)
833 constr
->sblock_nalloc
= over_alloc_dd(dd
->ncg_home
+1);
834 srenew(constr
->sblock
, constr
->sblock_nalloc
);
838 iatom
= ilcon
->iatoms
;
841 for (c
= 0; c
< ncons
; c
++)
843 if (c
== 0 || iatom
[1] >= cgs
->index
[cg
+1])
845 constr
->sblock
[constr
->nblocks
++] = 3*c
;
846 while (iatom
[1] >= cgs
->index
[cg
+1])
853 constr
->sblock
[constr
->nblocks
] = 3*ncons
;
856 t_blocka
make_at2con(int start
, int natoms
,
857 t_ilist
*ilist
, t_iparams
*iparams
,
858 gmx_bool bDynamics
, int *nflexiblecons
)
860 int *count
, ncon
, con
, con_tot
, nflexcon
, ftype
, i
, a
;
867 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
869 ncon
= ilist
[ftype
].nr
/3;
870 ia
= ilist
[ftype
].iatoms
;
871 for (con
= 0; con
< ncon
; con
++)
873 bFlexCon
= (iparams
[ia
[0]].constr
.dA
== 0 &&
874 iparams
[ia
[0]].constr
.dB
== 0);
879 if (bDynamics
|| !bFlexCon
)
881 for (i
= 1; i
< 3; i
++)
890 *nflexiblecons
= nflexcon
;
893 at2con
.nalloc_index
= at2con
.nr
+1;
894 snew(at2con
.index
, at2con
.nalloc_index
);
896 for (a
= 0; a
< natoms
; a
++)
898 at2con
.index
[a
+1] = at2con
.index
[a
] + count
[a
];
901 at2con
.nra
= at2con
.index
[natoms
];
902 at2con
.nalloc_a
= at2con
.nra
;
903 snew(at2con
.a
, at2con
.nalloc_a
);
905 /* The F_CONSTRNC constraints have constraint numbers
906 * that continue after the last F_CONSTR constraint.
909 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
911 ncon
= ilist
[ftype
].nr
/3;
912 ia
= ilist
[ftype
].iatoms
;
913 for (con
= 0; con
< ncon
; con
++)
915 bFlexCon
= (iparams
[ia
[0]].constr
.dA
== 0 &&
916 iparams
[ia
[0]].constr
.dB
== 0);
917 if (bDynamics
|| !bFlexCon
)
919 for (i
= 1; i
< 3; i
++)
922 at2con
.a
[at2con
.index
[a
]+count
[a
]++] = con_tot
;
935 static int *make_at2settle(int natoms
, const t_ilist
*ilist
)
941 /* Set all to no settle */
942 for (a
= 0; a
< natoms
; a
++)
947 stride
= 1 + NRAL(F_SETTLE
);
949 for (s
= 0; s
< ilist
->nr
; s
+= stride
)
951 at2s
[ilist
->iatoms
[s
+1]] = s
/stride
;
952 at2s
[ilist
->iatoms
[s
+2]] = s
/stride
;
953 at2s
[ilist
->iatoms
[s
+3]] = s
/stride
;
959 void set_constraints(struct gmx_constr
*constr
,
960 gmx_localtop_t
*top
, t_inputrec
*ir
,
961 t_mdatoms
*md
, t_commrec
*cr
)
970 if (constr
->ncon_tot
> 0)
972 /* We are using the local topology,
973 * so there are only F_CONSTR constraints.
975 ncons
= idef
->il
[F_CONSTR
].nr
/3;
977 /* With DD we might also need to call LINCS with ncons=0 for
978 * communicating coordinates to other nodes that do have constraints.
980 if (ir
->eConstrAlg
== econtLINCS
)
982 set_lincs(idef
, md
, EI_DYNAMICS(ir
->eI
), cr
, constr
->lincsd
);
984 if (ir
->eConstrAlg
== econtSHAKE
)
988 make_shake_sblock_dd(constr
, &idef
->il
[F_CONSTR
], &top
->cgs
, cr
->dd
);
992 make_shake_sblock_serial(constr
, idef
, md
);
994 if (ncons
> constr
->lagr_nalloc
)
996 constr
->lagr_nalloc
= over_alloc_dd(ncons
);
997 srenew(constr
->lagr
, constr
->lagr_nalloc
);
1002 if (idef
->il
[F_SETTLE
].nr
> 0 && constr
->settled
== NULL
)
1004 settle
= &idef
->il
[F_SETTLE
];
1005 iO
= settle
->iatoms
[1];
1006 iH
= settle
->iatoms
[2];
1008 settle_init(md
->massT
[iO
], md
->massT
[iH
],
1009 md
->invmass
[iO
], md
->invmass
[iH
],
1010 idef
->iparams
[settle
->iatoms
[0]].settle
.doh
,
1011 idef
->iparams
[settle
->iatoms
[0]].settle
.dhh
);
1014 /* Make a selection of the local atoms for essential dynamics */
1015 if (constr
->ed
&& cr
->dd
)
1017 dd_make_local_ed_indices(cr
->dd
, constr
->ed
);
1021 static void constr_recur(t_blocka
*at2con
,
1022 t_ilist
*ilist
, t_iparams
*iparams
, gmx_bool bTopB
,
1023 int at
, int depth
, int nc
, int *path
,
1024 real r0
, real r1
, real
*r2max
,
1036 ncon1
= ilist
[F_CONSTR
].nr
/3;
1037 ia1
= ilist
[F_CONSTR
].iatoms
;
1038 ia2
= ilist
[F_CONSTRNC
].iatoms
;
1040 /* Loop over all constraints connected to this atom */
1041 for (c
= at2con
->index
[at
]; c
< at2con
->index
[at
+1]; c
++)
1044 /* Do not walk over already used constraints */
1046 for (a1
= 0; a1
< depth
; a1
++)
1048 if (con
== path
[a1
])
1055 ia
= constr_iatomptr(ncon1
, ia1
, ia2
, con
);
1056 /* Flexible constraints currently have length 0, which is incorrect */
1059 len
= iparams
[ia
[0]].constr
.dA
;
1063 len
= iparams
[ia
[0]].constr
.dB
;
1065 /* In the worst case the bond directions alternate */
1076 /* Assume angles of 120 degrees between all bonds */
1077 if (rn0
*rn0
+ rn1
*rn1
+ rn0
*rn1
> *r2max
)
1079 *r2max
= rn0
*rn0
+ rn1
*rn1
+ r0
*rn1
;
1082 fprintf(debug
, "Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0
, rn1
, sqrt(*r2max
));
1083 for (a1
= 0; a1
< depth
; a1
++)
1085 fprintf(debug
, " %d %5.3f",
1087 iparams
[constr_iatomptr(ncon1
, ia1
, ia2
, con
)[0]].constr
.dA
);
1089 fprintf(debug
, " %d %5.3f\n", con
, len
);
1092 /* Limit the number of recursions to 1000*nc,
1093 * so a call does not take more than a second,
1094 * even for highly connected systems.
1096 if (depth
+ 1 < nc
&& *count
< 1000*nc
)
1108 constr_recur(at2con
, ilist
, iparams
,
1109 bTopB
, a1
, depth
+1, nc
, path
, rn0
, rn1
, r2max
, count
);
1116 static real
constr_r_max_moltype(gmx_moltype_t
*molt
, t_iparams
*iparams
,
1119 int natoms
, nflexcon
, *path
, at
, count
;
1122 real r0
, r1
, r2maxA
, r2maxB
, rmax
, lam0
, lam1
;
1124 if (molt
->ilist
[F_CONSTR
].nr
== 0 &&
1125 molt
->ilist
[F_CONSTRNC
].nr
== 0)
1130 natoms
= molt
->atoms
.nr
;
1132 at2con
= make_at2con(0, natoms
, molt
->ilist
, iparams
,
1133 EI_DYNAMICS(ir
->eI
), &nflexcon
);
1134 snew(path
, 1+ir
->nProjOrder
);
1135 for (at
= 0; at
< 1+ir
->nProjOrder
; at
++)
1141 for (at
= 0; at
< natoms
; at
++)
1147 constr_recur(&at2con
, molt
->ilist
, iparams
,
1148 FALSE
, at
, 0, 1+ir
->nProjOrder
, path
, r0
, r1
, &r2maxA
, &count
);
1150 if (ir
->efep
== efepNO
)
1152 rmax
= sqrt(r2maxA
);
1157 for (at
= 0; at
< natoms
; at
++)
1162 constr_recur(&at2con
, molt
->ilist
, iparams
,
1163 TRUE
, at
, 0, 1+ir
->nProjOrder
, path
, r0
, r1
, &r2maxB
, &count
);
1165 lam0
= ir
->fepvals
->init_lambda
;
1166 if (EI_DYNAMICS(ir
->eI
))
1168 lam0
+= ir
->init_step
*ir
->fepvals
->delta_lambda
;
1170 rmax
= (1 - lam0
)*sqrt(r2maxA
) + lam0
*sqrt(r2maxB
);
1171 if (EI_DYNAMICS(ir
->eI
))
1173 lam1
= ir
->fepvals
->init_lambda
+ (ir
->init_step
+ ir
->nsteps
)*ir
->fepvals
->delta_lambda
;
1174 rmax
= max(rmax
, (1 - lam1
)*sqrt(r2maxA
) + lam1
*sqrt(r2maxB
));
1178 done_blocka(&at2con
);
1184 real
constr_r_max(FILE *fplog
, gmx_mtop_t
*mtop
, t_inputrec
*ir
)
1190 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1193 constr_r_max_moltype(&mtop
->moltype
[mt
],
1194 mtop
->ffparams
.iparams
, ir
));
1199 fprintf(fplog
, "Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n", 1+ir
->nProjOrder
, rmax
);
1205 gmx_constr_t
init_constraints(FILE *fplog
,
1206 gmx_mtop_t
*mtop
, t_inputrec
*ir
,
1207 gmx_edsam_t ed
, t_state
*state
,
1210 int ncon
, nset
, nmol
, settle_type
, i
, natoms
, mt
, nflexcon
;
1211 struct gmx_constr
*constr
;
1214 gmx_mtop_ilistloop_t iloop
;
1217 gmx_mtop_ftype_count(mtop
, F_CONSTR
) +
1218 gmx_mtop_ftype_count(mtop
, F_CONSTRNC
);
1219 nset
= gmx_mtop_ftype_count(mtop
, F_SETTLE
);
1221 if (ncon
+nset
== 0 && ir
->ePull
!= epullCONSTRAINT
&& ed
== NULL
)
1228 constr
->ncon_tot
= ncon
;
1229 constr
->nflexcon
= 0;
1232 constr
->n_at2con_mt
= mtop
->nmoltype
;
1233 snew(constr
->at2con_mt
, constr
->n_at2con_mt
);
1234 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1236 constr
->at2con_mt
[mt
] = make_at2con(0, mtop
->moltype
[mt
].atoms
.nr
,
1237 mtop
->moltype
[mt
].ilist
,
1238 mtop
->ffparams
.iparams
,
1239 EI_DYNAMICS(ir
->eI
), &nflexcon
);
1240 for (i
= 0; i
< mtop
->nmolblock
; i
++)
1242 if (mtop
->molblock
[i
].type
== mt
)
1244 constr
->nflexcon
+= mtop
->molblock
[i
].nmol
*nflexcon
;
1249 if (constr
->nflexcon
> 0)
1253 fprintf(fplog
, "There are %d flexible constraints\n",
1255 if (ir
->fc_stepsize
== 0)
1258 "WARNING: step size for flexible constraining = 0\n"
1259 " All flexible constraints will be rigid.\n"
1260 " Will try to keep all flexible constraints at their original length,\n"
1261 " but the lengths may exhibit some drift.\n\n");
1262 constr
->nflexcon
= 0;
1265 if (constr
->nflexcon
> 0)
1267 please_cite(fplog
, "Hess2002");
1271 if (ir
->eConstrAlg
== econtLINCS
)
1273 constr
->lincsd
= init_lincs(fplog
, mtop
,
1274 constr
->nflexcon
, constr
->at2con_mt
,
1275 DOMAINDECOMP(cr
) && cr
->dd
->bInterCGcons
,
1276 ir
->nLincsIter
, ir
->nProjOrder
);
1279 if (ir
->eConstrAlg
== econtSHAKE
)
1281 if (DOMAINDECOMP(cr
) && cr
->dd
->bInterCGcons
)
1283 gmx_fatal(FARGS
, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1285 if (constr
->nflexcon
)
1287 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");
1289 please_cite(fplog
, "Ryckaert77a");
1292 please_cite(fplog
, "Barth95a");
1295 constr
->shaked
= shake_init();
1301 please_cite(fplog
, "Miyamoto92a");
1303 constr
->bInterCGsettles
= inter_charge_group_settles(mtop
);
1305 /* Check that we have only one settle type */
1307 iloop
= gmx_mtop_ilistloop_init(mtop
);
1308 while (gmx_mtop_ilistloop_next(iloop
, &ilist
, &nmol
))
1310 for (i
= 0; i
< ilist
[F_SETTLE
].nr
; i
+= 4)
1312 if (settle_type
== -1)
1314 settle_type
= ilist
[F_SETTLE
].iatoms
[i
];
1316 else if (ilist
[F_SETTLE
].iatoms
[i
] != settle_type
)
1319 "The [molecules] section of your topology specifies more than one block of\n"
1320 "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1321 "are trying to partition your solvent into different *groups* (e.g. for\n"
1322 "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1323 "files specify groups. Otherwise, you may wish to change the least-used\n"
1324 "block of molecules with SETTLE constraints into 3 normal constraints.");
1329 constr
->n_at2settle_mt
= mtop
->nmoltype
;
1330 snew(constr
->at2settle_mt
, constr
->n_at2settle_mt
);
1331 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1333 constr
->at2settle_mt
[mt
] =
1334 make_at2settle(mtop
->moltype
[mt
].atoms
.nr
,
1335 &mtop
->moltype
[mt
].ilist
[F_SETTLE
]);
1339 constr
->maxwarn
= 999;
1340 env
= getenv("GMX_MAXCONSTRWARN");
1343 constr
->maxwarn
= 0;
1344 sscanf(env
, "%d", &constr
->maxwarn
);
1348 "Setting the maximum number of constraint warnings to %d\n",
1354 "Setting the maximum number of constraint warnings to %d\n",
1358 if (constr
->maxwarn
< 0 && fplog
)
1360 fprintf(fplog
, "maxwarn < 0, will not stop on constraint errors\n");
1362 constr
->warncount_lincs
= 0;
1363 constr
->warncount_settle
= 0;
1365 /* Initialize the essential dynamics sampling.
1366 * Put the pointer to the ED struct in constr */
1368 if (ed
!= NULL
|| state
->edsamstate
.nED
> 0)
1370 init_edsam(mtop
, ir
, cr
, ed
, state
->x
, state
->box
, &state
->edsamstate
);
1373 constr
->warn_mtop
= mtop
;
1378 const t_blocka
*atom2constraints_moltype(gmx_constr_t constr
)
1380 return constr
->at2con_mt
;
1383 const int **atom2settle_moltype(gmx_constr_t constr
)
1385 return (const int **)constr
->at2settle_mt
;
1389 gmx_bool
inter_charge_group_constraints(const gmx_mtop_t
*mtop
)
1391 const gmx_moltype_t
*molt
;
1395 int nat
, *at2cg
, cg
, a
, ftype
, i
;
1399 for (mb
= 0; mb
< mtop
->nmolblock
&& !bInterCG
; mb
++)
1401 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
1403 if (molt
->ilist
[F_CONSTR
].nr
> 0 ||
1404 molt
->ilist
[F_CONSTRNC
].nr
> 0 ||
1405 molt
->ilist
[F_SETTLE
].nr
> 0)
1408 snew(at2cg
, molt
->atoms
.nr
);
1409 for (cg
= 0; cg
< cgs
->nr
; cg
++)
1411 for (a
= cgs
->index
[cg
]; a
< cgs
->index
[cg
+1]; a
++)
1417 for (ftype
= F_CONSTR
; ftype
<= F_CONSTRNC
; ftype
++)
1419 il
= &molt
->ilist
[ftype
];
1420 for (i
= 0; i
< il
->nr
&& !bInterCG
; i
+= 1+NRAL(ftype
))
1422 if (at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+2]])
1436 gmx_bool
inter_charge_group_settles(const gmx_mtop_t
*mtop
)
1438 const gmx_moltype_t
*molt
;
1442 int nat
, *at2cg
, cg
, a
, ftype
, i
;
1446 for (mb
= 0; mb
< mtop
->nmolblock
&& !bInterCG
; mb
++)
1448 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
1450 if (molt
->ilist
[F_SETTLE
].nr
> 0)
1453 snew(at2cg
, molt
->atoms
.nr
);
1454 for (cg
= 0; cg
< cgs
->nr
; cg
++)
1456 for (a
= cgs
->index
[cg
]; a
< cgs
->index
[cg
+1]; a
++)
1462 for (ftype
= F_SETTLE
; ftype
<= F_SETTLE
; ftype
++)
1464 il
= &molt
->ilist
[ftype
];
1465 for (i
= 0; i
< il
->nr
&& !bInterCG
; i
+= 1+NRAL(F_SETTLE
))
1467 if (at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+2]] ||
1468 at2cg
[il
->iatoms
[i
+1]] != at2cg
[il
->iatoms
[i
+3]])
1482 /* helper functions for andersen temperature control, because the
1483 * gmx_constr construct is only defined in constr.c. Return the list
1484 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
1486 extern int *get_sblock(struct gmx_constr
*constr
)
1488 return constr
->sblock
;
1491 extern int get_nblocks(struct gmx_constr
*constr
)
1493 return constr
->nblocks
;