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, 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 */
116 //! Compares sort blocks.
117 static int pcomp(const void *p1
, const void *p2
)
120 int min1
, min2
, max1
, max2
;
121 t_sortblock
*a1
= (t_sortblock
*)p1
;
122 t_sortblock
*a2
= (t_sortblock
*)p2
;
124 db
= a1
->blocknr
-a2
->blocknr
;
131 min1
= std::min(a1
->iatom
[1], a1
->iatom
[2]);
132 max1
= std::max(a1
->iatom
[1], a1
->iatom
[2]);
133 min2
= std::min(a2
->iatom
[1], a2
->iatom
[2]);
134 max2
= std::max(a2
->iatom
[1], a2
->iatom
[2]);
146 //! Prints sortblocks
147 static void pr_sortblock(FILE *fp
, const char *title
, int nsb
, t_sortblock sb
[])
151 fprintf(fp
, "%s\n", title
);
152 for (i
= 0; (i
< nsb
); i
++)
154 fprintf(fp
, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
155 i
, sb
[i
].iatom
[0], sb
[i
].iatom
[1], sb
[i
].iatom
[2],
160 //! Reallocates a vector.
161 static void resizeLagrangianData(shakedata
*shaked
, int ncons
)
163 if (ncons
> shaked
->lagr_nalloc
)
165 shaked
->lagr_nalloc
= over_alloc_dd(ncons
);
166 srenew(shaked
->scaled_lagrange_multiplier
, shaked
->lagr_nalloc
);
171 make_shake_sblock_serial(shakedata
*shaked
,
172 const t_idef
*idef
, const t_mdatoms
&md
)
181 /* Since we are processing the local topology,
182 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
184 ncons
= idef
->il
[F_CONSTR
].nr
/3;
186 init_blocka(&sblocks
);
187 gen_sblocks(nullptr, 0, md
.homenr
, idef
, &sblocks
, FALSE
);
190 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
191 nblocks=blocks->multinr[idef->nodeid] - bstart;
194 shaked
->nblocks
= sblocks
.nr
;
197 fprintf(debug
, "ncons: %d, bstart: %d, nblocks: %d\n",
198 ncons
, bstart
, shaked
->nblocks
);
201 /* Calculate block number for each atom */
202 inv_sblock
= make_invblocka(&sblocks
, md
.nr
);
204 done_blocka(&sblocks
);
206 /* Store the block number in temp array and
207 * sort the constraints in order of the sblock number
208 * and the atom numbers, really sorting a segment of the array!
210 iatom
= idef
->il
[F_CONSTR
].iatoms
;
212 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
214 for (m
= 0; (m
< 3); m
++)
216 sb
[i
].iatom
[m
] = iatom
[m
];
218 sb
[i
].blocknr
= inv_sblock
[iatom
[1]];
221 /* Now sort the blocks */
224 pr_sortblock(debug
, "Before sorting", ncons
, sb
);
225 fprintf(debug
, "Going to sort constraints\n");
228 qsort(sb
, ncons
, (size_t)sizeof(*sb
), pcomp
);
232 pr_sortblock(debug
, "After sorting", ncons
, sb
);
235 iatom
= idef
->il
[F_CONSTR
].iatoms
;
236 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
238 for (m
= 0; (m
< 3); m
++)
240 iatom
[m
] = sb
[i
].iatom
[m
];
245 snew(shaked
->sblock
, shaked
->nblocks
+1);
247 for (i
= 0; (i
< ncons
); i
++)
249 if (sb
[i
].blocknr
!= bnr
)
252 shaked
->sblock
[j
++] = 3*i
;
256 shaked
->sblock
[j
++] = 3*ncons
;
258 if (j
!= (shaked
->nblocks
+1))
260 fprintf(stderr
, "bstart: %d\n", bstart
);
261 fprintf(stderr
, "j: %d, nblocks: %d, ncons: %d\n",
262 j
, shaked
->nblocks
, ncons
);
263 for (i
= 0; (i
< ncons
); i
++)
265 fprintf(stderr
, "i: %5d sb[i].blocknr: %5d\n", i
, sb
[i
].blocknr
);
267 for (j
= 0; (j
<= shaked
->nblocks
); j
++)
269 fprintf(stderr
, "sblock[%3d]=%5d\n", j
, (int)shaked
->sblock
[j
]);
271 gmx_fatal(FARGS
, "DEATH HORROR: "
272 "sblocks does not match idef->il[F_CONSTR]");
276 resizeLagrangianData(shaked
, ncons
);
280 make_shake_sblock_dd(shakedata
*shaked
,
281 const t_ilist
*ilcon
, const t_block
*cgs
,
282 const gmx_domdec_t
*dd
)
287 if (dd
->ncg_home
+1 > shaked
->sblock_nalloc
)
289 shaked
->sblock_nalloc
= over_alloc_dd(dd
->ncg_home
+1);
290 srenew(shaked
->sblock
, shaked
->sblock_nalloc
);
294 iatom
= ilcon
->iatoms
;
297 for (c
= 0; c
< ncons
; c
++)
299 if (c
== 0 || iatom
[1] >= cgs
->index
[cg
+1])
301 shaked
->sblock
[shaked
->nblocks
++] = 3*c
;
302 while (iatom
[1] >= cgs
->index
[cg
+1])
309 shaked
->sblock
[shaked
->nblocks
] = 3*ncons
;
310 resizeLagrangianData(shaked
, ncons
);
313 /*! \brief Inner kernel for SHAKE constraints
315 * Original implementation from R.C. van Schaik and W.F. van Gunsteren
316 * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
317 * Spoel November 1992.
319 * The algorithm here is based section five of Ryckaert, Ciccotti and
320 * Berendsen, J Comp Phys, 23, 327, 1977.
322 * \param[in] iatom Mini-topology of triples of constraint type (unused in this
323 * function) and indices of the two atoms involved
324 * \param[in] ncon Number of constraints
325 * \param[out] nnit Number of iterations performed
326 * \param[in] maxnit Maximum number of iterations permitted
327 * \param[in] constraint_distance_squared The objective value for each constraint
328 * \param[inout] positions The initial (and final) values of the positions of all atoms
329 * \param[in] initial_displacements The initial displacements of each constraint
330 * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
331 * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
332 * using shake-sor=yes in the .mdp, but there is no documentation anywhere)
333 * \param[in] invmass Inverse mass of each atom
334 * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
335 * square of the constrained distance (see code)
336 * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
337 * of the paper, divided by the constraint distance)
338 * \param[out] nerror Zero upon success, returns one more than the index of the
339 * problematic constraint if the input was malformed
341 * \todo Make SHAKE use better data structures, in particular for iatom. */
342 void cshake(const int iatom
[], int ncon
, int *nnit
, int maxnit
,
343 const real constraint_distance_squared
[], real positions
[],
344 const real initial_displacements
[], const real half_of_reduced_mass
[], real omega
,
345 const real invmass
[], const real distance_squared_tolerance
[],
346 real scaled_lagrange_multiplier
[], int *nerror
)
348 /* default should be increased! MRS 8/4/2009 */
349 const real mytol
= 1e-10;
351 int ll
, i
, j
, i3
, j3
, l3
;
352 int ix
, iy
, iz
, jx
, jy
, jz
;
354 real constraint_distance_squared_ll
;
355 real r_prime_squared
;
356 real scaled_lagrange_multiplier_ll
;
357 real r_prime_x
, r_prime_y
, r_prime_z
, diff
, im
, jm
;
358 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
359 int nit
, error
, nconv
;
362 // TODO nconv is used solely as a boolean, so we should write the
366 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
369 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
372 rijx
= initial_displacements
[l3
+XX
];
373 rijy
= initial_displacements
[l3
+YY
];
374 rijz
= initial_displacements
[l3
+ZZ
];
386 /* Compute r prime between atoms i and j, which is the
387 displacement *before* this update stage */
388 r_prime_x
= positions
[ix
]-positions
[jx
];
389 r_prime_y
= positions
[iy
]-positions
[jy
];
390 r_prime_z
= positions
[iz
]-positions
[jz
];
391 r_prime_squared
= (r_prime_x
* r_prime_x
+
392 r_prime_y
* r_prime_y
+
393 r_prime_z
* r_prime_z
);
394 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
395 diff
= constraint_distance_squared_ll
- r_prime_squared
;
397 /* iconvf is less than 1 when the error is smaller than a bound */
398 iconvf
= fabs(diff
) * distance_squared_tolerance
[ll
];
402 nconv
= static_cast<int>(iconvf
);
403 r_dot_r_prime
= (rijx
* r_prime_x
+
407 if (r_dot_r_prime
< constraint_distance_squared_ll
* mytol
)
413 /* The next line solves equation 5.6 (neglecting
414 the term in g^2), for g */
415 scaled_lagrange_multiplier_ll
= omega
*diff
*half_of_reduced_mass
[ll
]/r_dot_r_prime
;
416 scaled_lagrange_multiplier
[ll
] += scaled_lagrange_multiplier_ll
;
417 xh
= rijx
* scaled_lagrange_multiplier_ll
;
418 yh
= rijy
* scaled_lagrange_multiplier_ll
;
419 zh
= rijz
* scaled_lagrange_multiplier_ll
;
422 positions
[ix
] += xh
*im
;
423 positions
[iy
] += yh
*im
;
424 positions
[iz
] += zh
*im
;
425 positions
[jx
] -= xh
*jm
;
426 positions
[jy
] -= yh
*jm
;
427 positions
[jz
] -= zh
*jm
;
436 //! Implements RATTLE (ie. SHAKE for velocity verlet integrators)
438 crattle(int iatom
[], int ncon
, int *nnit
, int maxnit
,
439 real constraint_distance_squared
[], real vp
[], real rij
[], real m2
[], real omega
,
440 const real invmass
[], real distance_squared_tolerance
[], real scaled_lagrange_multiplier
[],
441 int *nerror
, real invdt
)
444 * r.c. van schaik and w.f. van gunsteren
447 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
448 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
449 * second part of rattle algorithm
452 int ll
, i
, j
, i3
, j3
, l3
;
453 int ix
, iy
, iz
, jx
, jy
, jz
;
454 real constraint_distance_squared_ll
;
455 real vpijd
, vx
, vy
, vz
, acor
, fac
, im
, jm
;
456 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
457 int nit
, error
, nconv
;
460 // TODO nconv is used solely as a boolean, so we should write the
464 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
467 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
487 vpijd
= vx
*rijx
+vy
*rijy
+vz
*rijz
;
488 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
490 /* iconv is zero when the error is smaller than a bound */
491 iconvf
= fabs(vpijd
)*(distance_squared_tolerance
[ll
]/invdt
);
495 nconv
= static_cast<int>(iconvf
);
496 fac
= omega
*2.0*m2
[ll
]/constraint_distance_squared_ll
;
498 scaled_lagrange_multiplier
[ll
] += acor
;
520 static int vec_shakef(FILE *fplog
, shakedata
*shaked
,
521 const real invmass
[], int ncon
,
522 t_iparams ip
[], t_iatom
*iatom
,
523 real tol
, const rvec x
[], rvec prime
[], real omega
,
524 bool bFEP
, real lambda
, real scaled_lagrange_multiplier
[],
526 bool bCalcVir
, tensor vir_r_m_dr
, ConstraintVariable econq
)
529 real
*half_of_reduced_mass
, *distance_squared_tolerance
, *constraint_distance_squared
;
531 int nit
= 0, ll
, i
, j
, d
, d2
, type
;
536 real constraint_distance
;
538 if (ncon
> shaked
->nalloc
)
540 shaked
->nalloc
= over_alloc_dd(ncon
);
541 srenew(shaked
->rij
, shaked
->nalloc
);
542 srenew(shaked
->half_of_reduced_mass
, shaked
->nalloc
);
543 srenew(shaked
->distance_squared_tolerance
, shaked
->nalloc
);
544 srenew(shaked
->constraint_distance_squared
, shaked
->nalloc
);
547 half_of_reduced_mass
= shaked
->half_of_reduced_mass
;
548 distance_squared_tolerance
= shaked
->distance_squared_tolerance
;
549 constraint_distance_squared
= shaked
->constraint_distance_squared
;
553 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
559 mm
= 2.0*(invmass
[i
]+invmass
[j
]);
560 rij
[ll
][XX
] = x
[i
][XX
]-x
[j
][XX
];
561 rij
[ll
][YY
] = x
[i
][YY
]-x
[j
][YY
];
562 rij
[ll
][ZZ
] = x
[i
][ZZ
]-x
[j
][ZZ
];
563 half_of_reduced_mass
[ll
] = 1.0/mm
;
566 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
570 constraint_distance
= ip
[type
].constr
.dA
;
572 constraint_distance_squared
[ll
] = gmx::square(constraint_distance
);
573 distance_squared_tolerance
[ll
] = 0.5/(constraint_distance_squared
[ll
]*tol
);
578 case ConstraintVariable::Positions
:
579 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
);
581 case ConstraintVariable::Velocities
:
582 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
);
585 gmx_incons("Unknown constraint quantity for SHAKE");
592 fprintf(fplog
, "Shake did not converge in %d steps\n", maxnit
);
594 fprintf(stderr
, "Shake did not converge in %d steps\n", maxnit
);
601 fprintf(fplog
, "Inner product between old and new vector <= 0.0!\n"
602 "constraint #%d atoms %d and %d\n",
603 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
605 fprintf(stderr
, "Inner product between old and new vector <= 0.0!\n"
606 "constraint #%d atoms %d and %d\n",
607 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
611 /* Constraint virial and correct the Lagrange multipliers for the length */
615 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
621 if ((econq
== ConstraintVariable::Positions
) && v
!= nullptr)
623 /* Correct the velocities */
624 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[i
]*invdt
;
625 for (d
= 0; d
< DIM
; d
++)
627 v
[ia
[1]][d
] += mm
*rij
[ll
][d
];
629 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[j
]*invdt
;
630 for (d
= 0; d
< DIM
; d
++)
632 v
[ia
[2]][d
] -= mm
*rij
[ll
][d
];
637 /* constraint virial */
640 mm
= scaled_lagrange_multiplier
[ll
];
641 for (d
= 0; d
< DIM
; d
++)
644 for (d2
= 0; d2
< DIM
; d2
++)
646 vir_r_m_dr
[d
][d2
] -= tmp
*rij
[ll
][d2
];
652 /* cshake and crattle produce Lagrange multipliers scaled by
653 the reciprocal of the constraint length, so fix that */
656 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
660 constraint_distance
= ip
[type
].constr
.dA
;
662 scaled_lagrange_multiplier
[ll
] *= constraint_distance
;
668 //! Check that constraints are satisfied.
669 static void check_cons(FILE *log
, int nc
, const rvec x
[], rvec prime
[], rvec v
[],
670 t_iparams ip
[], t_iatom
*iatom
,
671 const real invmass
[], ConstraintVariable econq
)
680 " i mi j mj before after should be\n");
682 for (i
= 0; (i
< nc
); i
++, ia
+= 3)
686 rvec_sub(x
[ai
], x
[aj
], dx
);
691 case ConstraintVariable::Positions
:
692 rvec_sub(prime
[ai
], prime
[aj
], dx
);
694 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
695 ai
+1, 1.0/invmass
[ai
],
696 aj
+1, 1.0/invmass
[aj
], d
, dp
, ip
[ia
[0]].constr
.dA
);
698 case ConstraintVariable::Velocities
:
699 rvec_sub(v
[ai
], v
[aj
], dv
);
701 rvec_sub(prime
[ai
], prime
[aj
], dv
);
703 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
704 ai
+1, 1.0/invmass
[ai
],
705 aj
+1, 1.0/invmass
[aj
], d
, dp
, 0.);
708 gmx_incons("Unknown constraint quantity for SHAKE");
715 bshakef(FILE *log
, shakedata
*shaked
,
716 const real invmass
[],
717 const t_idef
&idef
, const t_inputrec
&ir
, const rvec x_s
[], rvec prime
[],
718 t_nrnb
*nrnb
, real lambda
, real
*dvdlambda
,
719 real invdt
, rvec
*v
, bool bCalcVir
, tensor vir_r_m_dr
,
720 bool bDumpOnError
, ConstraintVariable econq
)
723 real
*lam
, dt_2
, dvdl
;
724 int i
, n0
, ncon
, blen
, type
, ll
;
725 int tnit
= 0, trij
= 0;
727 ncon
= idef
.il
[F_CONSTR
].nr
/3;
729 for (ll
= 0; ll
< ncon
; ll
++)
731 shaked
->scaled_lagrange_multiplier
[ll
] = 0;
734 // TODO Rewrite this block so that it is obvious that i, iatoms
735 // and lam are all iteration variables. Is this easier if the
736 // sblock data structure is organized differently?
737 iatoms
= &(idef
.il
[F_CONSTR
].iatoms
[shaked
->sblock
[0]]);
738 lam
= shaked
->scaled_lagrange_multiplier
;
739 for (i
= 0; (i
< shaked
->nblocks
); )
741 blen
= (shaked
->sblock
[i
+1]-shaked
->sblock
[i
]);
743 n0
= vec_shakef(log
, shaked
, invmass
, blen
, idef
.iparams
,
744 iatoms
, ir
.shake_tol
, x_s
, prime
, shaked
->omega
,
745 ir
.efep
!= efepNO
, lambda
, lam
, invdt
, v
, bCalcVir
, vir_r_m_dr
,
750 if (bDumpOnError
&& log
)
753 check_cons(log
, blen
, x_s
, prime
, v
, idef
.iparams
, iatoms
, invmass
, econq
);
760 iatoms
+= 3*blen
; /* Increment pointer! */
764 /* only for position part? */
765 if (econq
== ConstraintVariable::Positions
)
767 if (ir
.efep
!= efepNO
)
770 /* TODO This should probably use invdt, so that sd integrator scaling works properly */
771 dt_2
= 1/gmx::square(ir
.delta_t
);
773 for (ll
= 0; ll
< ncon
; ll
++)
775 type
= idef
.il
[F_CONSTR
].iatoms
[3*ll
];
777 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
778 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
779 estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
780 so the pre-factors are already present. */
781 bondA
= idef
.iparams
[type
].constr
.dA
;
782 bondB
= idef
.iparams
[type
].constr
.dB
;
783 dvdl
+= shaked
->scaled_lagrange_multiplier
[ll
] * dt_2
* (bondB
- bondA
);
790 if (tnit
> shaked
->gamma
)
792 shaked
->delta
*= -0.5;
794 shaked
->omega
+= shaked
->delta
;
795 shaked
->gamma
= tnit
;
797 inc_nrnb(nrnb
, eNR_SHAKE
, tnit
);
798 inc_nrnb(nrnb
, eNR_SHAKE_RIJ
, trij
);
801 inc_nrnb(nrnb
, eNR_CONSTR_V
, trij
*2);
805 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, trij
);
812 constrain_shake(FILE *log
,
814 const real invmass
[],
816 const t_inputrec
&ir
,
828 ConstraintVariable econq
)
830 if (shaked
->nblocks
== 0)
837 case (ConstraintVariable::Positions
):
838 bOK
= bshakef(log
, shaked
,
840 idef
, ir
, x_s
, xprime
, nrnb
,
842 invdt
, v
, bCalcVir
, vir_r_m_dr
,
843 bDumpOnError
, econq
);
845 case (ConstraintVariable::Velocities
):
846 bOK
= bshakef(log
, shaked
,
848 idef
, ir
, x_s
, vprime
, nrnb
,
850 invdt
, nullptr, bCalcVir
, vir_r_m_dr
,
851 bDumpOnError
, econq
);
854 gmx_fatal(FARGS
, "Internal error, SHAKE called for constraining something else than coordinates");