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,2015,2017,2018,2019, 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.
38 * \brief Defines SHAKE code.
40 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
41 * \author Berk Hess <hess@kth.se>
42 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 * \ingroup module_mdlib
53 #include "gromacs/domdec/domdec_struct.h"
54 #include "gromacs/gmxlib/nrnb.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/constr.h"
58 #include "gromacs/mdlib/splitter.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/mdtypes/mdatom.h"
62 #include "gromacs/topology/invblock.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.h"
72 real
*half_of_reduced_mass
;
73 real
*distance_squared_tolerance
;
74 real
*constraint_distance_squared
;
80 int nblocks
; /* The number of SHAKE blocks */
81 int *sblock
; /* The SHAKE blocks */
82 int sblock_nalloc
; /* The allocation size of sblock */
83 /*! \brief Scaled Lagrange multiplier for each constraint.
85 * Value is -2 * eta from p. 336 of the paper, divided by the
86 * constraint distance. */
87 real
*scaled_lagrange_multiplier
;
88 int lagr_nalloc
; /* The allocation size of scaled_lagrange_multiplier */
91 shakedata
*shake_init()
99 d
->half_of_reduced_mass
= nullptr;
100 d
->distance_squared_tolerance
= nullptr;
101 d
->constraint_distance_squared
= nullptr;
103 /* SOR initialization */
111 void done_shake(shakedata
*d
)
114 sfree(d
->half_of_reduced_mass
);
115 sfree(d
->distance_squared_tolerance
);
116 sfree(d
->constraint_distance_squared
);
118 sfree(d
->scaled_lagrange_multiplier
);
127 //! Compares sort blocks.
128 static int pcomp(const void *p1
, const void *p2
)
131 int min1
, min2
, max1
, max2
;
132 const t_sortblock
*a1
= reinterpret_cast<const t_sortblock
*>(p1
);
133 const t_sortblock
*a2
= reinterpret_cast<const t_sortblock
*>(p2
);
135 db
= a1
->blocknr
-a2
->blocknr
;
142 min1
= std::min(a1
->iatom
[1], a1
->iatom
[2]);
143 max1
= std::max(a1
->iatom
[1], a1
->iatom
[2]);
144 min2
= std::min(a2
->iatom
[1], a2
->iatom
[2]);
145 max2
= std::max(a2
->iatom
[1], a2
->iatom
[2]);
157 //! Prints sortblocks
158 static void pr_sortblock(FILE *fp
, const char *title
, int nsb
, t_sortblock sb
[])
162 fprintf(fp
, "%s\n", title
);
163 for (i
= 0; (i
< nsb
); i
++)
165 fprintf(fp
, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
166 i
, sb
[i
].iatom
[0], sb
[i
].iatom
[1], sb
[i
].iatom
[2],
171 //! Reallocates a vector.
172 static void resizeLagrangianData(shakedata
*shaked
, int ncons
)
174 if (ncons
> shaked
->lagr_nalloc
)
176 shaked
->lagr_nalloc
= over_alloc_dd(ncons
);
177 srenew(shaked
->scaled_lagrange_multiplier
, shaked
->lagr_nalloc
);
182 make_shake_sblock_serial(shakedata
*shaked
,
183 const t_idef
*idef
, const t_mdatoms
&md
)
192 /* Since we are processing the local topology,
193 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
195 ncons
= idef
->il
[F_CONSTR
].nr
/3;
197 init_blocka(&sblocks
);
198 sfree(sblocks
.index
); // To solve memory leak
199 gen_sblocks(nullptr, 0, md
.homenr
, idef
, &sblocks
, FALSE
);
202 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
203 nblocks=blocks->multinr[idef->nodeid] - bstart;
206 shaked
->nblocks
= sblocks
.nr
;
209 fprintf(debug
, "ncons: %d, bstart: %d, nblocks: %d\n",
210 ncons
, bstart
, shaked
->nblocks
);
213 /* Calculate block number for each atom */
214 inv_sblock
= make_invblocka(&sblocks
, md
.nr
);
216 done_blocka(&sblocks
);
218 /* Store the block number in temp array and
219 * sort the constraints in order of the sblock number
220 * and the atom numbers, really sorting a segment of the array!
222 iatom
= idef
->il
[F_CONSTR
].iatoms
;
224 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
226 for (m
= 0; (m
< 3); m
++)
228 sb
[i
].iatom
[m
] = iatom
[m
];
230 sb
[i
].blocknr
= inv_sblock
[iatom
[1]];
233 /* Now sort the blocks */
236 pr_sortblock(debug
, "Before sorting", ncons
, sb
);
237 fprintf(debug
, "Going to sort constraints\n");
240 std::qsort(sb
, ncons
, sizeof(*sb
), pcomp
);
244 pr_sortblock(debug
, "After sorting", ncons
, sb
);
247 iatom
= idef
->il
[F_CONSTR
].iatoms
;
248 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
250 for (m
= 0; (m
< 3); m
++)
252 iatom
[m
] = sb
[i
].iatom
[m
];
257 snew(shaked
->sblock
, shaked
->nblocks
+1);
259 for (i
= 0; (i
< ncons
); i
++)
261 if (sb
[i
].blocknr
!= bnr
)
264 shaked
->sblock
[j
++] = 3*i
;
268 shaked
->sblock
[j
++] = 3*ncons
;
270 if (j
!= (shaked
->nblocks
+1))
272 fprintf(stderr
, "bstart: %d\n", bstart
);
273 fprintf(stderr
, "j: %d, nblocks: %d, ncons: %d\n",
274 j
, shaked
->nblocks
, ncons
);
275 for (i
= 0; (i
< ncons
); i
++)
277 fprintf(stderr
, "i: %5d sb[i].blocknr: %5d\n", i
, sb
[i
].blocknr
);
279 for (j
= 0; (j
<= shaked
->nblocks
); j
++)
281 fprintf(stderr
, "sblock[%3d]=%5d\n", j
, shaked
->sblock
[j
]);
283 gmx_fatal(FARGS
, "DEATH HORROR: "
284 "sblocks does not match idef->il[F_CONSTR]");
288 resizeLagrangianData(shaked
, ncons
);
291 // TODO: Check if this code is useful. It might never be called.
293 make_shake_sblock_dd(shakedata
*shaked
,
294 const t_ilist
*ilcon
,
295 const gmx_domdec_t
*dd
)
300 if (dd
->ncg_home
+1 > shaked
->sblock_nalloc
)
302 shaked
->sblock_nalloc
= over_alloc_dd(dd
->ncg_home
+1);
303 srenew(shaked
->sblock
, shaked
->sblock_nalloc
);
307 iatom
= ilcon
->iatoms
;
310 for (c
= 0; c
< ncons
; c
++)
312 if (c
== 0 || iatom
[1] >= cg
+ 1)
314 shaked
->sblock
[shaked
->nblocks
++] = 3*c
;
315 while (iatom
[1] >= cg
+ 1)
322 shaked
->sblock
[shaked
->nblocks
] = 3*ncons
;
323 resizeLagrangianData(shaked
, ncons
);
326 /*! \brief Inner kernel for SHAKE constraints
328 * Original implementation from R.C. van Schaik and W.F. van Gunsteren
329 * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
330 * Spoel November 1992.
332 * The algorithm here is based section five of Ryckaert, Ciccotti and
333 * Berendsen, J Comp Phys, 23, 327, 1977.
335 * \param[in] iatom Mini-topology of triples of constraint type (unused in this
336 * function) and indices of the two atoms involved
337 * \param[in] ncon Number of constraints
338 * \param[out] nnit Number of iterations performed
339 * \param[in] maxnit Maximum number of iterations permitted
340 * \param[in] constraint_distance_squared The objective value for each constraint
341 * \param[inout] positions The initial (and final) values of the positions of all atoms
342 * \param[in] initial_displacements The initial displacements of each constraint
343 * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
344 * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
345 * using shake-sor=yes in the .mdp, but there is no documentation anywhere)
346 * \param[in] invmass Inverse mass of each atom
347 * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
348 * square of the constrained distance (see code)
349 * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
350 * of the paper, divided by the constraint distance)
351 * \param[out] nerror Zero upon success, returns one more than the index of the
352 * problematic constraint if the input was malformed
354 * \todo Make SHAKE use better data structures, in particular for iatom. */
355 void cshake(const int iatom
[], int ncon
, int *nnit
, int maxnit
,
356 const real constraint_distance_squared
[], real positions
[],
357 const real initial_displacements
[], const real half_of_reduced_mass
[], real omega
,
358 const real invmass
[], const real distance_squared_tolerance
[],
359 real scaled_lagrange_multiplier
[], int *nerror
)
361 /* default should be increased! MRS 8/4/2009 */
362 const real mytol
= 1e-10;
364 int ll
, i
, j
, i3
, j3
, l3
;
365 int ix
, iy
, iz
, jx
, jy
, jz
;
367 real constraint_distance_squared_ll
;
368 real r_prime_squared
;
369 real scaled_lagrange_multiplier_ll
;
370 real r_prime_x
, r_prime_y
, r_prime_z
, diff
, im
, jm
;
371 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
372 int nit
, error
, nconv
;
375 // TODO nconv is used solely as a boolean, so we should write the
379 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
382 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
385 rijx
= initial_displacements
[l3
+XX
];
386 rijy
= initial_displacements
[l3
+YY
];
387 rijz
= initial_displacements
[l3
+ZZ
];
399 /* Compute r prime between atoms i and j, which is the
400 displacement *before* this update stage */
401 r_prime_x
= positions
[ix
]-positions
[jx
];
402 r_prime_y
= positions
[iy
]-positions
[jy
];
403 r_prime_z
= positions
[iz
]-positions
[jz
];
404 r_prime_squared
= (r_prime_x
* r_prime_x
+
405 r_prime_y
* r_prime_y
+
406 r_prime_z
* r_prime_z
);
407 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
408 diff
= constraint_distance_squared_ll
- r_prime_squared
;
410 /* iconvf is less than 1 when the error is smaller than a bound */
411 iconvf
= fabs(diff
) * distance_squared_tolerance
[ll
];
415 nconv
= static_cast<int>(iconvf
);
416 r_dot_r_prime
= (rijx
* r_prime_x
+
420 if (r_dot_r_prime
< constraint_distance_squared_ll
* mytol
)
426 /* The next line solves equation 5.6 (neglecting
427 the term in g^2), for g */
428 scaled_lagrange_multiplier_ll
= omega
*diff
*half_of_reduced_mass
[ll
]/r_dot_r_prime
;
429 scaled_lagrange_multiplier
[ll
] += scaled_lagrange_multiplier_ll
;
430 xh
= rijx
* scaled_lagrange_multiplier_ll
;
431 yh
= rijy
* scaled_lagrange_multiplier_ll
;
432 zh
= rijz
* scaled_lagrange_multiplier_ll
;
435 positions
[ix
] += xh
*im
;
436 positions
[iy
] += yh
*im
;
437 positions
[iz
] += zh
*im
;
438 positions
[jx
] -= xh
*jm
;
439 positions
[jy
] -= yh
*jm
;
440 positions
[jz
] -= zh
*jm
;
449 //! Implements RATTLE (ie. SHAKE for velocity verlet integrators)
451 crattle(const int iatom
[], int ncon
, int *nnit
, int maxnit
,
452 const real constraint_distance_squared
[], real vp
[], const real rij
[], const real m2
[], real omega
,
453 const real invmass
[], const real distance_squared_tolerance
[], real scaled_lagrange_multiplier
[],
454 int *nerror
, real invdt
)
457 * r.c. van schaik and w.f. van gunsteren
460 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
461 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
462 * second part of rattle algorithm
465 int ll
, i
, j
, i3
, j3
, l3
;
466 int ix
, iy
, iz
, jx
, jy
, jz
;
467 real constraint_distance_squared_ll
;
468 real vpijd
, vx
, vy
, vz
, acor
, fac
, im
, jm
;
469 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
470 int nit
, error
, nconv
;
473 // TODO nconv is used solely as a boolean, so we should write the
477 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
480 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
500 vpijd
= vx
*rijx
+vy
*rijy
+vz
*rijz
;
501 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
503 /* iconv is zero when the error is smaller than a bound */
504 iconvf
= fabs(vpijd
)*(distance_squared_tolerance
[ll
]/invdt
);
508 nconv
= static_cast<int>(iconvf
);
509 fac
= omega
*2.0*m2
[ll
]/constraint_distance_squared_ll
;
511 scaled_lagrange_multiplier
[ll
] += acor
;
533 static int vec_shakef(FILE *fplog
, shakedata
*shaked
,
534 const real invmass
[], int ncon
,
535 t_iparams ip
[], t_iatom
*iatom
,
536 real tol
, const rvec x
[], rvec prime
[], real omega
,
537 bool bFEP
, real lambda
, real scaled_lagrange_multiplier
[],
539 bool bCalcVir
, tensor vir_r_m_dr
, ConstraintVariable econq
)
542 real
*half_of_reduced_mass
, *distance_squared_tolerance
, *constraint_distance_squared
;
544 int nit
= 0, ll
, i
, j
, d
, d2
, type
;
549 real constraint_distance
;
551 if (ncon
> shaked
->nalloc
)
553 shaked
->nalloc
= over_alloc_dd(ncon
);
554 srenew(shaked
->rij
, shaked
->nalloc
);
555 srenew(shaked
->half_of_reduced_mass
, shaked
->nalloc
);
556 srenew(shaked
->distance_squared_tolerance
, shaked
->nalloc
);
557 srenew(shaked
->constraint_distance_squared
, shaked
->nalloc
);
560 half_of_reduced_mass
= shaked
->half_of_reduced_mass
;
561 distance_squared_tolerance
= shaked
->distance_squared_tolerance
;
562 constraint_distance_squared
= shaked
->constraint_distance_squared
;
566 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
572 mm
= 2.0*(invmass
[i
]+invmass
[j
]);
573 rij
[ll
][XX
] = x
[i
][XX
]-x
[j
][XX
];
574 rij
[ll
][YY
] = x
[i
][YY
]-x
[j
][YY
];
575 rij
[ll
][ZZ
] = x
[i
][ZZ
]-x
[j
][ZZ
];
576 half_of_reduced_mass
[ll
] = 1.0/mm
;
579 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
583 constraint_distance
= ip
[type
].constr
.dA
;
585 constraint_distance_squared
[ll
] = gmx::square(constraint_distance
);
586 distance_squared_tolerance
[ll
] = 0.5/(constraint_distance_squared
[ll
]*tol
);
591 case ConstraintVariable::Positions
:
592 cshake(iatom
, ncon
, &nit
, maxnit
, constraint_distance_squared
, prime
[0], rij
[0], half_of_reduced_mass
, omega
, invmass
, distance_squared_tolerance
, scaled_lagrange_multiplier
, &error
);
594 case ConstraintVariable::Velocities
:
595 crattle(iatom
, ncon
, &nit
, maxnit
, constraint_distance_squared
, prime
[0], rij
[0], half_of_reduced_mass
, omega
, invmass
, distance_squared_tolerance
, scaled_lagrange_multiplier
, &error
, invdt
);
598 gmx_incons("Unknown constraint quantity for SHAKE");
605 fprintf(fplog
, "Shake did not converge in %d steps\n", maxnit
);
607 fprintf(stderr
, "Shake did not converge in %d steps\n", maxnit
);
614 fprintf(fplog
, "Inner product between old and new vector <= 0.0!\n"
615 "constraint #%d atoms %d and %d\n",
616 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
618 fprintf(stderr
, "Inner product between old and new vector <= 0.0!\n"
619 "constraint #%d atoms %d and %d\n",
620 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
624 /* Constraint virial and correct the Lagrange multipliers for the length */
628 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
634 if ((econq
== ConstraintVariable::Positions
) && v
!= nullptr)
636 /* Correct the velocities */
637 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[i
]*invdt
;
638 for (d
= 0; d
< DIM
; d
++)
640 v
[ia
[1]][d
] += mm
*rij
[ll
][d
];
642 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[j
]*invdt
;
643 for (d
= 0; d
< DIM
; d
++)
645 v
[ia
[2]][d
] -= mm
*rij
[ll
][d
];
650 /* constraint virial */
653 mm
= scaled_lagrange_multiplier
[ll
];
654 for (d
= 0; d
< DIM
; d
++)
657 for (d2
= 0; d2
< DIM
; d2
++)
659 vir_r_m_dr
[d
][d2
] -= tmp
*rij
[ll
][d2
];
665 /* cshake and crattle produce Lagrange multipliers scaled by
666 the reciprocal of the constraint length, so fix that */
669 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
673 constraint_distance
= ip
[type
].constr
.dA
;
675 scaled_lagrange_multiplier
[ll
] *= constraint_distance
;
681 //! Check that constraints are satisfied.
682 static void check_cons(FILE *log
, int nc
, const rvec x
[], rvec prime
[], rvec v
[],
683 t_iparams ip
[], t_iatom
*iatom
,
684 const real invmass
[], ConstraintVariable econq
)
692 GMX_ASSERT(v
, "Input has to be non-null");
694 " i mi j mj before after should be\n");
696 for (i
= 0; (i
< nc
); i
++, ia
+= 3)
700 rvec_sub(x
[ai
], x
[aj
], dx
);
705 case ConstraintVariable::Positions
:
706 rvec_sub(prime
[ai
], prime
[aj
], dx
);
708 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
709 ai
+1, 1.0/invmass
[ai
],
710 aj
+1, 1.0/invmass
[aj
], d
, dp
, ip
[ia
[0]].constr
.dA
);
712 case ConstraintVariable::Velocities
:
713 rvec_sub(v
[ai
], v
[aj
], dv
);
715 rvec_sub(prime
[ai
], prime
[aj
], dv
);
717 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
718 ai
+1, 1.0/invmass
[ai
],
719 aj
+1, 1.0/invmass
[aj
], d
, dp
, 0.);
722 gmx_incons("Unknown constraint quantity for SHAKE");
729 bshakef(FILE *log
, shakedata
*shaked
,
730 const real invmass
[],
731 const t_idef
&idef
, const t_inputrec
&ir
, const rvec x_s
[], rvec prime
[],
732 t_nrnb
*nrnb
, real lambda
, real
*dvdlambda
,
733 real invdt
, rvec
*v
, bool bCalcVir
, tensor vir_r_m_dr
,
734 bool bDumpOnError
, ConstraintVariable econq
)
737 real
*lam
, dt_2
, dvdl
;
738 int i
, n0
, ncon
, blen
, type
, ll
;
739 int tnit
= 0, trij
= 0;
741 ncon
= idef
.il
[F_CONSTR
].nr
/3;
743 for (ll
= 0; ll
< ncon
; ll
++)
745 shaked
->scaled_lagrange_multiplier
[ll
] = 0;
748 // TODO Rewrite this block so that it is obvious that i, iatoms
749 // and lam are all iteration variables. Is this easier if the
750 // sblock data structure is organized differently?
751 iatoms
= &(idef
.il
[F_CONSTR
].iatoms
[shaked
->sblock
[0]]);
752 lam
= shaked
->scaled_lagrange_multiplier
;
753 for (i
= 0; (i
< shaked
->nblocks
); )
755 blen
= (shaked
->sblock
[i
+1]-shaked
->sblock
[i
]);
757 n0
= vec_shakef(log
, shaked
, invmass
, blen
, idef
.iparams
,
758 iatoms
, ir
.shake_tol
, x_s
, prime
, shaked
->omega
,
759 ir
.efep
!= efepNO
, lambda
, lam
, invdt
, v
, bCalcVir
, vir_r_m_dr
,
764 if (bDumpOnError
&& log
)
767 check_cons(log
, blen
, x_s
, prime
, v
, idef
.iparams
, iatoms
, invmass
, econq
);
774 iatoms
+= 3*blen
; /* Increment pointer! */
778 /* only for position part? */
779 if (econq
== ConstraintVariable::Positions
)
781 if (ir
.efep
!= efepNO
)
784 /* TODO This should probably use invdt, so that sd integrator scaling works properly */
785 dt_2
= 1/gmx::square(ir
.delta_t
);
787 for (ll
= 0; ll
< ncon
; ll
++)
789 type
= idef
.il
[F_CONSTR
].iatoms
[3*ll
];
791 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
792 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
793 estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
794 so the pre-factors are already present. */
795 bondA
= idef
.iparams
[type
].constr
.dA
;
796 bondB
= idef
.iparams
[type
].constr
.dB
;
797 dvdl
+= shaked
->scaled_lagrange_multiplier
[ll
] * dt_2
* (bondB
- bondA
);
804 if (tnit
> shaked
->gamma
)
806 shaked
->delta
*= -0.5;
808 shaked
->omega
+= shaked
->delta
;
809 shaked
->gamma
= tnit
;
811 inc_nrnb(nrnb
, eNR_SHAKE
, tnit
);
812 inc_nrnb(nrnb
, eNR_SHAKE_RIJ
, trij
);
815 inc_nrnb(nrnb
, eNR_CONSTR_V
, trij
*2);
819 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, trij
);
826 constrain_shake(FILE *log
,
828 const real invmass
[],
830 const t_inputrec
&ir
,
842 ConstraintVariable econq
)
844 if (shaked
->nblocks
== 0)
851 case (ConstraintVariable::Positions
):
852 bOK
= bshakef(log
, shaked
,
854 idef
, ir
, x_s
, xprime
, nrnb
,
856 invdt
, v
, bCalcVir
, vir_r_m_dr
,
857 bDumpOnError
, econq
);
859 case (ConstraintVariable::Velocities
):
860 bOK
= bshakef(log
, shaked
,
862 idef
, ir
, x_s
, vprime
, nrnb
,
864 invdt
, nullptr, bCalcVir
, vir_r_m_dr
,
865 bDumpOnError
, econq
);
868 gmx_fatal(FARGS
, "Internal error, SHAKE called for constraining something else than coordinates");