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.
45 #include "gromacs/domdec/domdec_struct.h"
46 #include "gromacs/gmxlib/nrnb.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/mdlib/splitter.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/mdtypes/mdatom.h"
53 #include "gromacs/topology/invblock.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/smalloc.h"
57 typedef struct gmx_shakedata
60 real
*half_of_reduced_mass
;
61 real
*distance_squared_tolerance
;
62 real
*constraint_distance_squared
;
68 int nblocks
; /* The number of SHAKE blocks */
69 int *sblock
; /* The SHAKE blocks */
70 int sblock_nalloc
; /* The allocation size of sblock */
71 /*! \brief Scaled Lagrange multiplier for each constraint.
73 * Value is -2 * eta from p. 336 of the paper, divided by the
74 * constraint distance. */
75 real
*scaled_lagrange_multiplier
;
76 int lagr_nalloc
; /* The allocation size of scaled_lagrange_multiplier */
79 gmx_shakedata_t
shake_init()
87 d
->half_of_reduced_mass
= nullptr;
88 d
->distance_squared_tolerance
= nullptr;
89 d
->constraint_distance_squared
= nullptr;
91 /* SOR initialization */
104 static int pcomp(const void *p1
, const void *p2
)
107 int min1
, min2
, max1
, max2
;
108 t_sortblock
*a1
= (t_sortblock
*)p1
;
109 t_sortblock
*a2
= (t_sortblock
*)p2
;
111 db
= a1
->blocknr
-a2
->blocknr
;
118 min1
= std::min(a1
->iatom
[1], a1
->iatom
[2]);
119 max1
= std::max(a1
->iatom
[1], a1
->iatom
[2]);
120 min2
= std::min(a2
->iatom
[1], a2
->iatom
[2]);
121 max2
= std::max(a2
->iatom
[1], a2
->iatom
[2]);
133 static void pr_sortblock(FILE *fp
, const char *title
, int nsb
, t_sortblock sb
[])
137 fprintf(fp
, "%s\n", title
);
138 for (i
= 0; (i
< nsb
); i
++)
140 fprintf(fp
, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
141 i
, sb
[i
].iatom
[0], sb
[i
].iatom
[1], sb
[i
].iatom
[2],
146 static void resizeLagrangianData(gmx_shakedata
*shaked
, int ncons
)
148 if (ncons
> shaked
->lagr_nalloc
)
150 shaked
->lagr_nalloc
= over_alloc_dd(ncons
);
151 srenew(shaked
->scaled_lagrange_multiplier
, shaked
->lagr_nalloc
);
156 make_shake_sblock_serial(gmx_shakedata
*shaked
,
157 const t_idef
*idef
, const t_mdatoms
*md
)
166 /* Since we are processing the local topology,
167 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
169 ncons
= idef
->il
[F_CONSTR
].nr
/3;
171 init_blocka(&sblocks
);
172 gen_sblocks(nullptr, 0, md
->homenr
, idef
, &sblocks
, FALSE
);
175 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
176 nblocks=blocks->multinr[idef->nodeid] - bstart;
179 shaked
->nblocks
= sblocks
.nr
;
182 fprintf(debug
, "ncons: %d, bstart: %d, nblocks: %d\n",
183 ncons
, bstart
, shaked
->nblocks
);
186 /* Calculate block number for each atom */
187 inv_sblock
= make_invblocka(&sblocks
, md
->nr
);
189 done_blocka(&sblocks
);
191 /* Store the block number in temp array and
192 * sort the constraints in order of the sblock number
193 * and the atom numbers, really sorting a segment of the array!
195 iatom
= idef
->il
[F_CONSTR
].iatoms
;
197 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
199 for (m
= 0; (m
< 3); m
++)
201 sb
[i
].iatom
[m
] = iatom
[m
];
203 sb
[i
].blocknr
= inv_sblock
[iatom
[1]];
206 /* Now sort the blocks */
209 pr_sortblock(debug
, "Before sorting", ncons
, sb
);
210 fprintf(debug
, "Going to sort constraints\n");
213 qsort(sb
, ncons
, (size_t)sizeof(*sb
), pcomp
);
217 pr_sortblock(debug
, "After sorting", ncons
, sb
);
220 iatom
= idef
->il
[F_CONSTR
].iatoms
;
221 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
223 for (m
= 0; (m
< 3); m
++)
225 iatom
[m
] = sb
[i
].iatom
[m
];
230 snew(shaked
->sblock
, shaked
->nblocks
+1);
232 for (i
= 0; (i
< ncons
); i
++)
234 if (sb
[i
].blocknr
!= bnr
)
237 shaked
->sblock
[j
++] = 3*i
;
241 shaked
->sblock
[j
++] = 3*ncons
;
243 if (j
!= (shaked
->nblocks
+1))
245 fprintf(stderr
, "bstart: %d\n", bstart
);
246 fprintf(stderr
, "j: %d, nblocks: %d, ncons: %d\n",
247 j
, shaked
->nblocks
, ncons
);
248 for (i
= 0; (i
< ncons
); i
++)
250 fprintf(stderr
, "i: %5d sb[i].blocknr: %5d\n", i
, sb
[i
].blocknr
);
252 for (j
= 0; (j
<= shaked
->nblocks
); j
++)
254 fprintf(stderr
, "sblock[%3d]=%5d\n", j
, (int)shaked
->sblock
[j
]);
256 gmx_fatal(FARGS
, "DEATH HORROR: "
257 "sblocks does not match idef->il[F_CONSTR]");
261 resizeLagrangianData(shaked
, ncons
);
265 make_shake_sblock_dd(gmx_shakedata
*shaked
,
266 const t_ilist
*ilcon
, const t_block
*cgs
,
267 const gmx_domdec_t
*dd
)
272 if (dd
->ncg_home
+1 > shaked
->sblock_nalloc
)
274 shaked
->sblock_nalloc
= over_alloc_dd(dd
->ncg_home
+1);
275 srenew(shaked
->sblock
, shaked
->sblock_nalloc
);
279 iatom
= ilcon
->iatoms
;
282 for (c
= 0; c
< ncons
; c
++)
284 if (c
== 0 || iatom
[1] >= cgs
->index
[cg
+1])
286 shaked
->sblock
[shaked
->nblocks
++] = 3*c
;
287 while (iatom
[1] >= cgs
->index
[cg
+1])
294 shaked
->sblock
[shaked
->nblocks
] = 3*ncons
;
295 resizeLagrangianData(shaked
, ncons
);
298 /*! \brief Inner kernel for SHAKE constraints
300 * Original implementation from R.C. van Schaik and W.F. van Gunsteren
301 * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
302 * Spoel November 1992.
304 * The algorithm here is based section five of Ryckaert, Ciccotti and
305 * Berendsen, J Comp Phys, 23, 327, 1977.
307 * \param[in] iatom Mini-topology of triples of constraint type (unused in this
308 * function) and indices of the two atoms involved
309 * \param[in] ncon Number of constraints
310 * \param[out] nnit Number of iterations performed
311 * \param[in] maxnit Maximum number of iterations permitted
312 * \param[in] constraint_distance_squared The objective value for each constraint
313 * \param[inout] positions The initial (and final) values of the positions of all atoms
314 * \param[in] initial_displacements The initial displacements of each constraint
315 * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
316 * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
317 * using shake-sor=yes in the .mdp, but there is no documentation anywhere)
318 * \param[in] invmass Inverse mass of each atom
319 * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
320 * square of the constrained distance (see code)
321 * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
322 * of the paper, divided by the constraint distance)
323 * \param[out] nerror Zero upon success, returns one more than the index of the
324 * problematic constraint if the input was malformed
326 * \todo Make SHAKE use better data structures, in particular for iatom. */
327 void cshake(const int iatom
[], int ncon
, int *nnit
, int maxnit
,
328 const real constraint_distance_squared
[], real positions
[],
329 const real initial_displacements
[], const real half_of_reduced_mass
[], real omega
,
330 const real invmass
[], const real distance_squared_tolerance
[],
331 real scaled_lagrange_multiplier
[], int *nerror
)
333 /* default should be increased! MRS 8/4/2009 */
334 const real mytol
= 1e-10;
336 int ll
, i
, j
, i3
, j3
, l3
;
337 int ix
, iy
, iz
, jx
, jy
, jz
;
339 real constraint_distance_squared_ll
;
340 real r_prime_squared
;
341 real scaled_lagrange_multiplier_ll
;
342 real r_prime_x
, r_prime_y
, r_prime_z
, diff
, im
, jm
;
343 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
344 int nit
, error
, nconv
;
347 // TODO nconv is used solely as a boolean, so we should write the
351 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
354 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
357 rijx
= initial_displacements
[l3
+XX
];
358 rijy
= initial_displacements
[l3
+YY
];
359 rijz
= initial_displacements
[l3
+ZZ
];
371 /* Compute r prime between atoms i and j, which is the
372 displacement *before* this update stage */
373 r_prime_x
= positions
[ix
]-positions
[jx
];
374 r_prime_y
= positions
[iy
]-positions
[jy
];
375 r_prime_z
= positions
[iz
]-positions
[jz
];
376 r_prime_squared
= (r_prime_x
* r_prime_x
+
377 r_prime_y
* r_prime_y
+
378 r_prime_z
* r_prime_z
);
379 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
380 diff
= constraint_distance_squared_ll
- r_prime_squared
;
382 /* iconvf is less than 1 when the error is smaller than a bound */
383 iconvf
= fabs(diff
) * distance_squared_tolerance
[ll
];
387 nconv
= static_cast<int>(iconvf
);
388 r_dot_r_prime
= (rijx
* r_prime_x
+
392 if (r_dot_r_prime
< constraint_distance_squared_ll
* mytol
)
398 /* The next line solves equation 5.6 (neglecting
399 the term in g^2), for g */
400 scaled_lagrange_multiplier_ll
= omega
*diff
*half_of_reduced_mass
[ll
]/r_dot_r_prime
;
401 scaled_lagrange_multiplier
[ll
] += scaled_lagrange_multiplier_ll
;
402 xh
= rijx
* scaled_lagrange_multiplier_ll
;
403 yh
= rijy
* scaled_lagrange_multiplier_ll
;
404 zh
= rijz
* scaled_lagrange_multiplier_ll
;
407 positions
[ix
] += xh
*im
;
408 positions
[iy
] += yh
*im
;
409 positions
[iz
] += zh
*im
;
410 positions
[jx
] -= xh
*jm
;
411 positions
[jy
] -= yh
*jm
;
412 positions
[jz
] -= zh
*jm
;
422 crattle(int iatom
[], int ncon
, int *nnit
, int maxnit
,
423 real constraint_distance_squared
[], real vp
[], real rij
[], real m2
[], real omega
,
424 const real invmass
[], real distance_squared_tolerance
[], real scaled_lagrange_multiplier
[],
425 int *nerror
, real invdt
)
428 * r.c. van schaik and w.f. van gunsteren
431 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
432 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
433 * second part of rattle algorithm
436 int ll
, i
, j
, i3
, j3
, l3
;
437 int ix
, iy
, iz
, jx
, jy
, jz
;
438 real constraint_distance_squared_ll
;
439 real vpijd
, vx
, vy
, vz
, acor
, fac
, im
, jm
;
440 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
441 int nit
, error
, nconv
;
444 // TODO nconv is used solely as a boolean, so we should write the
448 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
451 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
471 vpijd
= vx
*rijx
+vy
*rijy
+vz
*rijz
;
472 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
474 /* iconv is zero when the error is smaller than a bound */
475 iconvf
= fabs(vpijd
)*(distance_squared_tolerance
[ll
]/invdt
);
479 nconv
= static_cast<int>(iconvf
);
480 fac
= omega
*2.0*m2
[ll
]/constraint_distance_squared_ll
;
482 scaled_lagrange_multiplier
[ll
] += acor
;
503 static int vec_shakef(FILE *fplog
, gmx_shakedata_t shaked
,
504 const real invmass
[], int ncon
,
505 t_iparams ip
[], t_iatom
*iatom
,
506 real tol
, rvec x
[], rvec prime
[], real omega
,
507 gmx_bool bFEP
, real lambda
, real scaled_lagrange_multiplier
[],
509 gmx_bool bCalcVir
, tensor vir_r_m_dr
, int econq
)
512 real
*half_of_reduced_mass
, *distance_squared_tolerance
, *constraint_distance_squared
;
514 int nit
= 0, ll
, i
, j
, d
, d2
, type
;
519 real constraint_distance
;
521 if (ncon
> shaked
->nalloc
)
523 shaked
->nalloc
= over_alloc_dd(ncon
);
524 srenew(shaked
->rij
, shaked
->nalloc
);
525 srenew(shaked
->half_of_reduced_mass
, shaked
->nalloc
);
526 srenew(shaked
->distance_squared_tolerance
, shaked
->nalloc
);
527 srenew(shaked
->constraint_distance_squared
, shaked
->nalloc
);
530 half_of_reduced_mass
= shaked
->half_of_reduced_mass
;
531 distance_squared_tolerance
= shaked
->distance_squared_tolerance
;
532 constraint_distance_squared
= shaked
->constraint_distance_squared
;
536 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
542 mm
= 2.0*(invmass
[i
]+invmass
[j
]);
543 rij
[ll
][XX
] = x
[i
][XX
]-x
[j
][XX
];
544 rij
[ll
][YY
] = x
[i
][YY
]-x
[j
][YY
];
545 rij
[ll
][ZZ
] = x
[i
][ZZ
]-x
[j
][ZZ
];
546 half_of_reduced_mass
[ll
] = 1.0/mm
;
549 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
553 constraint_distance
= ip
[type
].constr
.dA
;
555 constraint_distance_squared
[ll
] = gmx::square(constraint_distance
);
556 distance_squared_tolerance
[ll
] = 0.5/(constraint_distance_squared
[ll
]*tol
);
562 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
);
565 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
);
573 fprintf(fplog
, "Shake did not converge in %d steps\n", maxnit
);
575 fprintf(stderr
, "Shake did not converge in %d steps\n", maxnit
);
582 fprintf(fplog
, "Inner product between old and new vector <= 0.0!\n"
583 "constraint #%d atoms %d and %d\n",
584 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
586 fprintf(stderr
, "Inner product between old and new vector <= 0.0!\n"
587 "constraint #%d atoms %d and %d\n",
588 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
592 /* Constraint virial and correct the Lagrange multipliers for the length */
596 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
602 if ((econq
== econqCoord
) && v
!= nullptr)
604 /* Correct the velocities */
605 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[i
]*invdt
;
606 for (d
= 0; d
< DIM
; d
++)
608 v
[ia
[1]][d
] += mm
*rij
[ll
][d
];
610 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[j
]*invdt
;
611 for (d
= 0; d
< DIM
; d
++)
613 v
[ia
[2]][d
] -= mm
*rij
[ll
][d
];
618 /* constraint virial */
621 mm
= scaled_lagrange_multiplier
[ll
];
622 for (d
= 0; d
< DIM
; d
++)
625 for (d2
= 0; d2
< DIM
; d2
++)
627 vir_r_m_dr
[d
][d2
] -= tmp
*rij
[ll
][d2
];
633 /* cshake and crattle produce Lagrange multipliers scaled by
634 the reciprocal of the constraint length, so fix that */
637 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
641 constraint_distance
= ip
[type
].constr
.dA
;
643 scaled_lagrange_multiplier
[ll
] *= constraint_distance
;
649 static void check_cons(FILE *log
, int nc
, rvec x
[], rvec prime
[], rvec v
[],
650 t_iparams ip
[], t_iatom
*iatom
,
651 const real invmass
[], int econq
)
660 " i mi j mj before after should be\n");
662 for (i
= 0; (i
< nc
); i
++, ia
+= 3)
666 rvec_sub(x
[ai
], x
[aj
], dx
);
672 rvec_sub(prime
[ai
], prime
[aj
], dx
);
674 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
675 ai
+1, 1.0/invmass
[ai
],
676 aj
+1, 1.0/invmass
[aj
], d
, dp
, ip
[ia
[0]].constr
.dA
);
679 rvec_sub(v
[ai
], v
[aj
], dv
);
681 rvec_sub(prime
[ai
], prime
[aj
], dv
);
683 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
684 ai
+1, 1.0/invmass
[ai
],
685 aj
+1, 1.0/invmass
[aj
], d
, dp
, 0.);
692 bshakef(FILE *log
, gmx_shakedata_t shaked
,
693 const real invmass
[],
694 const t_idef
*idef
, const t_inputrec
*ir
, rvec x_s
[], rvec prime
[],
695 t_nrnb
*nrnb
, real lambda
, real
*dvdlambda
,
696 real invdt
, rvec
*v
, gmx_bool bCalcVir
, tensor vir_r_m_dr
,
697 gmx_bool bDumpOnError
, int econq
)
700 real
*lam
, dt_2
, dvdl
;
701 int i
, n0
, ncon
, blen
, type
, ll
;
702 int tnit
= 0, trij
= 0;
704 ncon
= idef
->il
[F_CONSTR
].nr
/3;
706 for (ll
= 0; ll
< ncon
; ll
++)
708 shaked
->scaled_lagrange_multiplier
[ll
] = 0;
711 // TODO Rewrite this block so that it is obvious that i, iatoms
712 // and lam are all iteration variables. Is this easier if the
713 // sblock data structure is organized differently?
714 iatoms
= &(idef
->il
[F_CONSTR
].iatoms
[shaked
->sblock
[0]]);
715 lam
= shaked
->scaled_lagrange_multiplier
;
716 for (i
= 0; (i
< shaked
->nblocks
); )
718 blen
= (shaked
->sblock
[i
+1]-shaked
->sblock
[i
]);
720 n0
= vec_shakef(log
, shaked
, invmass
, blen
, idef
->iparams
,
721 iatoms
, ir
->shake_tol
, x_s
, prime
, shaked
->omega
,
722 ir
->efep
!= efepNO
, lambda
, lam
, invdt
, v
, bCalcVir
, vir_r_m_dr
,
727 if (bDumpOnError
&& log
)
730 check_cons(log
, blen
, x_s
, prime
, v
, idef
->iparams
, iatoms
, invmass
, econq
);
737 iatoms
+= 3*blen
; /* Increment pointer! */
741 /* only for position part? */
742 if (econq
== econqCoord
)
744 if (ir
->efep
!= efepNO
)
747 /* TODO This should probably use invdt, so that sd integrator scaling works properly */
748 dt_2
= 1/gmx::square(ir
->delta_t
);
750 for (ll
= 0; ll
< ncon
; ll
++)
752 type
= idef
->il
[F_CONSTR
].iatoms
[3*ll
];
754 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
755 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
756 estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
757 so the pre-factors are already present. */
758 bondA
= idef
->iparams
[type
].constr
.dA
;
759 bondB
= idef
->iparams
[type
].constr
.dB
;
760 dvdl
+= shaked
->scaled_lagrange_multiplier
[ll
] * dt_2
* (bondB
- bondA
);
767 if (tnit
> shaked
->gamma
)
769 shaked
->delta
*= -0.5;
771 shaked
->omega
+= shaked
->delta
;
772 shaked
->gamma
= tnit
;
774 inc_nrnb(nrnb
, eNR_SHAKE
, tnit
);
775 inc_nrnb(nrnb
, eNR_SHAKE_RIJ
, trij
);
778 inc_nrnb(nrnb
, eNR_CONSTR_V
, trij
*2);
782 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, trij
);
789 constrain_shake(FILE *log
,
790 gmx_shakedata_t shaked
,
791 const real invmass
[],
793 const t_inputrec
*ir
,
804 gmx_bool bDumpOnError
,
807 if (shaked
->nblocks
== 0)
815 bOK
= bshakef(log
, shaked
,
817 idef
, ir
, x_s
, xprime
, nrnb
,
819 invdt
, v
, bCalcVir
, vir_r_m_dr
,
820 bDumpOnError
, econq
);
823 bOK
= bshakef(log
, shaked
,
825 idef
, ir
, x_s
, vprime
, nrnb
,
827 invdt
, nullptr, bCalcVir
, vir_r_m_dr
,
828 bDumpOnError
, econq
);
831 gmx_fatal(FARGS
, "Internal error, SHAKE called for constraining something else than coordinates");