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.
7 * Copyright (c) 2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
39 * \brief Defines SHAKE code.
41 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
42 * \author Berk Hess <hess@kth.se>
43 * \author Mark Abraham <mark.j.abraham@gmail.com>
44 * \ingroup module_mdlib
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/pbcutil/pbc.h"
62 #include "gromacs/topology/invblock.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.h"
75 //! Compares sort blocks.
76 static int pcomp(const void* p1
, const void* p2
)
79 int min1
, min2
, max1
, max2
;
80 const t_sortblock
* a1
= reinterpret_cast<const t_sortblock
*>(p1
);
81 const t_sortblock
* a2
= reinterpret_cast<const t_sortblock
*>(p2
);
83 db
= a1
->blocknr
- a2
->blocknr
;
90 min1
= std::min(a1
->iatom
[1], a1
->iatom
[2]);
91 max1
= std::max(a1
->iatom
[1], a1
->iatom
[2]);
92 min2
= std::min(a2
->iatom
[1], a2
->iatom
[2]);
93 max2
= std::max(a2
->iatom
[1], a2
->iatom
[2]);
105 //! Prints sortblocks
106 static void pr_sortblock(FILE* fp
, const char* title
, int nsb
, t_sortblock sb
[])
110 fprintf(fp
, "%s\n", title
);
111 for (i
= 0; (i
< nsb
); i
++)
113 fprintf(fp
, "i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n", i
, sb
[i
].iatom
[0],
114 sb
[i
].iatom
[1], sb
[i
].iatom
[2], sb
[i
].blocknr
);
118 //! Reallocates a vector.
119 static void resizeLagrangianData(shakedata
* shaked
, int ncons
)
121 shaked
->scaled_lagrange_multiplier
.resize(ncons
);
124 void make_shake_sblock_serial(shakedata
* shaked
, InteractionDefinitions
* idef
, const int numAtoms
)
132 /* Since we are processing the local topology,
133 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
135 ncons
= idef
->il
[F_CONSTR
].size() / 3;
137 init_blocka(&sblocks
);
138 sfree(sblocks
.index
); // To solve memory leak
139 gen_sblocks(nullptr, numAtoms
, *idef
, &sblocks
, FALSE
);
142 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
143 nblocks=blocks->multinr[idef->nodeid] - bstart;
148 fprintf(debug
, "ncons: %d, bstart: %d, nblocks: %d\n", ncons
, bstart
, sblocks
.nr
);
151 /* Calculate block number for each atom */
152 inv_sblock
= make_invblocka(&sblocks
, numAtoms
);
154 done_blocka(&sblocks
);
156 /* Store the block number in temp array and
157 * sort the constraints in order of the sblock number
158 * and the atom numbers, really sorting a segment of the array!
160 int* iatom
= idef
->il
[F_CONSTR
].iatoms
.data();
162 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
164 for (m
= 0; (m
< 3); m
++)
166 sb
[i
].iatom
[m
] = iatom
[m
];
168 sb
[i
].blocknr
= inv_sblock
[iatom
[1]];
171 /* Now sort the blocks */
174 pr_sortblock(debug
, "Before sorting", ncons
, sb
);
175 fprintf(debug
, "Going to sort constraints\n");
178 std::qsort(sb
, ncons
, sizeof(*sb
), pcomp
);
182 pr_sortblock(debug
, "After sorting", ncons
, sb
);
185 iatom
= idef
->il
[F_CONSTR
].iatoms
.data();
186 for (i
= 0; (i
< ncons
); i
++, iatom
+= 3)
188 for (m
= 0; (m
< 3); m
++)
190 iatom
[m
] = sb
[i
].iatom
[m
];
194 shaked
->sblock
.clear();
196 for (i
= 0; (i
< ncons
); i
++)
198 if (sb
[i
].blocknr
!= bnr
)
201 shaked
->sblock
.push_back(3 * i
);
205 shaked
->sblock
.push_back(3 * ncons
);
209 resizeLagrangianData(shaked
, ncons
);
212 void make_shake_sblock_dd(shakedata
* shaked
, const InteractionList
& ilcon
)
216 ncons
= ilcon
.size() / 3;
217 const int* iatom
= ilcon
.iatoms
.data();
218 shaked
->sblock
.clear();
220 for (c
= 0; c
< ncons
; c
++)
222 if (c
== 0 || iatom
[1] >= cg
+ 1)
224 shaked
->sblock
.push_back(3 * c
);
225 while (iatom
[1] >= cg
+ 1)
232 shaked
->sblock
.push_back(3 * ncons
);
233 resizeLagrangianData(shaked
, ncons
);
236 /*! \brief Inner kernel for SHAKE constraints
238 * Original implementation from R.C. van Schaik and W.F. van Gunsteren
239 * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
240 * Spoel November 1992.
242 * The algorithm here is based section five of Ryckaert, Ciccotti and
243 * Berendsen, J Comp Phys, 23, 327, 1977.
245 * \param[in] iatom Mini-topology of triplets of constraint type (unused
246 * in this function) and indices of two atoms involved
247 * \param[in] ncon Number of constraints
248 * \param[out] nnit Number of iterations performed
249 * \param[in] maxnit Maximum number of iterations permitted
250 * \param[in] constraint_distance_squared The objective value for each constraint
251 * \param[inout] positions The initial (and final) values of the positions
253 * \param[in] pbc PBC information
254 * \param[in] initial_displacements The initial displacements of each constraint
255 * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
256 * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
257 * using shake-sor=yes in the .mdp,
258 * but there is no documentation anywhere)
259 * \param[in] invmass Inverse mass of each atom
260 * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
261 * square of the constrained distance (see code)
262 * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint
263 * (-2 * eta from p. 336 of the paper, divided by
264 * the constraint distance)
265 * \param[out] nerror Zero upon success, returns one more than the index of
266 * the problematic constraint if the input was malformed
268 * \todo Make SHAKE use better data structures, in particular for iatom. */
269 void cshake(const int iatom
[],
273 ArrayRef
<const real
> constraint_distance_squared
,
274 ArrayRef
<RVec
> positions
,
276 ArrayRef
<const RVec
> initial_displacements
,
277 ArrayRef
<const real
> half_of_reduced_mass
,
279 const real invmass
[],
280 ArrayRef
<const real
> distance_squared_tolerance
,
281 ArrayRef
<real
> scaled_lagrange_multiplier
,
284 /* default should be increased! MRS 8/4/2009 */
285 const real mytol
= 1e-10;
289 real constraint_distance_squared_ll
;
290 real scaled_lagrange_multiplier_ll
;
292 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
293 int nit
, error
, nconv
;
296 // TODO nconv is used solely as a boolean, so we should write the
300 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
303 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
306 rijx
= initial_displacements
[ll
][XX
];
307 rijy
= initial_displacements
[ll
][YY
];
308 rijz
= initial_displacements
[ll
][ZZ
];
312 /* Compute r prime between atoms i and j, which is the
313 displacement *before* this update stage */
317 pbc_dx(pbc
, positions
[i
], positions
[j
], r_prime
);
321 rvec_sub(positions
[i
], positions
[j
], r_prime
);
323 const real r_prime_squared
= norm2(r_prime
);
324 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
325 diff
= constraint_distance_squared_ll
- r_prime_squared
;
327 /* iconvf is less than 1 when the error is smaller than a bound */
328 iconvf
= fabs(diff
) * distance_squared_tolerance
[ll
];
332 nconv
= static_cast<int>(iconvf
);
333 r_dot_r_prime
= (rijx
* r_prime
[XX
] + rijy
* r_prime
[YY
] + rijz
* r_prime
[ZZ
]);
335 if (r_dot_r_prime
< constraint_distance_squared_ll
* mytol
)
341 /* The next line solves equation 5.6 (neglecting
342 the term in g^2), for g */
343 scaled_lagrange_multiplier_ll
= omega
* diff
* half_of_reduced_mass
[ll
] / r_dot_r_prime
;
344 scaled_lagrange_multiplier
[ll
] += scaled_lagrange_multiplier_ll
;
345 xh
= rijx
* scaled_lagrange_multiplier_ll
;
346 yh
= rijy
* scaled_lagrange_multiplier_ll
;
347 zh
= rijz
* scaled_lagrange_multiplier_ll
;
350 positions
[i
][XX
] += xh
* im
;
351 positions
[i
][YY
] += yh
* im
;
352 positions
[i
][ZZ
] += zh
* im
;
353 positions
[j
][XX
] -= xh
* jm
;
354 positions
[j
][YY
] -= yh
* jm
;
355 positions
[j
][ZZ
] -= zh
* jm
;
364 //! Implements RATTLE (ie. SHAKE for velocity verlet integrators)
365 static void crattle(const int iatom
[],
369 ArrayRef
<const real
> constraint_distance_squared
,
371 ArrayRef
<const RVec
> rij
,
372 ArrayRef
<const real
> m2
,
374 const real invmass
[],
375 ArrayRef
<const real
> distance_squared_tolerance
,
376 ArrayRef
<real
> scaled_lagrange_multiplier
,
381 * r.c. van schaik and w.f. van gunsteren
384 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
385 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
386 * second part of rattle algorithm
390 real constraint_distance_squared_ll
;
391 real vpijd
, acor
, fac
, im
, jm
;
392 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
393 int nit
, error
, nconv
;
396 // TODO nconv is used solely as a boolean, so we should write the
400 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
403 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
412 rvec_sub(vp
[i
], vp
[j
], v
);
414 vpijd
= v
[XX
] * rijx
+ v
[YY
] * rijy
+ v
[ZZ
] * rijz
;
415 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
417 /* iconv is zero when the error is smaller than a bound */
418 iconvf
= fabs(vpijd
) * (distance_squared_tolerance
[ll
] / invdt
);
422 nconv
= static_cast<int>(iconvf
);
423 fac
= omega
* 2.0 * m2
[ll
] / constraint_distance_squared_ll
;
425 scaled_lagrange_multiplier
[ll
] += acor
;
433 vp
[i
][XX
] += xh
* im
;
434 vp
[i
][YY
] += yh
* im
;
435 vp
[i
][ZZ
] += zh
* im
;
436 vp
[j
][XX
] -= xh
* jm
;
437 vp
[j
][YY
] -= yh
* jm
;
438 vp
[j
][ZZ
] -= zh
* jm
;
447 static int vec_shakef(FILE* fplog
,
449 const real invmass
[],
451 ArrayRef
<const t_iparams
> ip
,
454 ArrayRef
<const RVec
> x
,
455 ArrayRef
<RVec
> prime
,
460 ArrayRef
<real
> scaled_lagrange_multiplier
,
465 ConstraintVariable econq
)
468 int nit
= 0, ll
, i
, j
, d
, d2
, type
;
472 real constraint_distance
;
474 shaked
->rij
.resize(ncon
);
475 shaked
->half_of_reduced_mass
.resize(ncon
);
476 shaked
->distance_squared_tolerance
.resize(ncon
);
477 shaked
->constraint_distance_squared
.resize(ncon
);
479 ArrayRef
<RVec
> rij
= shaked
->rij
;
480 ArrayRef
<real
> half_of_reduced_mass
= shaked
->half_of_reduced_mass
;
481 ArrayRef
<real
> distance_squared_tolerance
= shaked
->distance_squared_tolerance
;
482 ArrayRef
<real
> constraint_distance_squared
= shaked
->constraint_distance_squared
;
485 const int* ia
= iatom
;
486 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
492 mm
= 2.0 * (invmass
[i
] + invmass
[j
]);
493 rij
[ll
][XX
] = x
[i
][XX
] - x
[j
][XX
];
494 rij
[ll
][YY
] = x
[i
][YY
] - x
[j
][YY
];
495 rij
[ll
][ZZ
] = x
[i
][ZZ
] - x
[j
][ZZ
];
496 half_of_reduced_mass
[ll
] = 1.0 / mm
;
499 constraint_distance
= L1
* ip
[type
].constr
.dA
+ lambda
* ip
[type
].constr
.dB
;
503 constraint_distance
= ip
[type
].constr
.dA
;
505 constraint_distance_squared
[ll
] = gmx::square(constraint_distance
);
506 distance_squared_tolerance
[ll
] = 0.5 / (constraint_distance_squared
[ll
] * tol
);
511 case ConstraintVariable::Positions
:
512 cshake(iatom
, ncon
, &nit
, maxnit
, constraint_distance_squared
, prime
, pbc
, rij
,
513 half_of_reduced_mass
, omega
, invmass
, distance_squared_tolerance
,
514 scaled_lagrange_multiplier
, &error
);
516 case ConstraintVariable::Velocities
:
517 crattle(iatom
, ncon
, &nit
, maxnit
, constraint_distance_squared
, prime
, rij
,
518 half_of_reduced_mass
, omega
, invmass
, distance_squared_tolerance
,
519 scaled_lagrange_multiplier
, &error
, invdt
);
521 default: gmx_incons("Unknown constraint quantity for SHAKE");
528 fprintf(fplog
, "Shake did not converge in %d steps\n", maxnit
);
530 fprintf(stderr
, "Shake did not converge in %d steps\n", maxnit
);
538 "Inner product between old and new vector <= 0.0!\n"
539 "constraint #%d atoms %d and %d\n",
540 error
- 1, iatom
[3 * (error
- 1) + 1] + 1, iatom
[3 * (error
- 1) + 2] + 1);
543 "Inner product between old and new vector <= 0.0!\n"
544 "constraint #%d atoms %d and %d\n",
545 error
- 1, iatom
[3 * (error
- 1) + 1] + 1, iatom
[3 * (error
- 1) + 2] + 1);
549 /* Constraint virial and correct the Lagrange multipliers for the length */
553 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
559 if ((econq
== ConstraintVariable::Positions
) && !v
.empty())
561 /* Correct the velocities */
562 mm
= scaled_lagrange_multiplier
[ll
] * invmass
[i
] * invdt
;
563 for (d
= 0; d
< DIM
; d
++)
565 v
[ia
[1]][d
] += mm
* rij
[ll
][d
];
567 mm
= scaled_lagrange_multiplier
[ll
] * invmass
[j
] * invdt
;
568 for (d
= 0; d
< DIM
; d
++)
570 v
[ia
[2]][d
] -= mm
* rij
[ll
][d
];
575 /* constraint virial */
578 mm
= scaled_lagrange_multiplier
[ll
];
579 for (d
= 0; d
< DIM
; d
++)
581 tmp
= mm
* rij
[ll
][d
];
582 for (d2
= 0; d2
< DIM
; d2
++)
584 vir_r_m_dr
[d
][d2
] -= tmp
* rij
[ll
][d2
];
590 /* cshake and crattle produce Lagrange multipliers scaled by
591 the reciprocal of the constraint length, so fix that */
594 constraint_distance
= L1
* ip
[type
].constr
.dA
+ lambda
* ip
[type
].constr
.dB
;
598 constraint_distance
= ip
[type
].constr
.dA
;
600 scaled_lagrange_multiplier
[ll
] *= constraint_distance
;
606 //! Check that constraints are satisfied.
607 static void check_cons(FILE* log
,
609 ArrayRef
<const RVec
> x
,
610 ArrayRef
<const RVec
> prime
,
611 ArrayRef
<const RVec
> v
,
613 ArrayRef
<const t_iparams
> ip
,
615 const real invmass
[],
616 ConstraintVariable econq
)
623 GMX_ASSERT(!v
.empty(), "Input has to be non-null");
624 fprintf(log
, " i mi j mj before after should be\n");
625 const int* ia
= iatom
;
626 for (i
= 0; (i
< nc
); i
++, ia
+= 3)
630 rvec_sub(x
[ai
], x
[aj
], dx
);
635 case ConstraintVariable::Positions
:
638 pbc_dx(pbc
, prime
[ai
], prime
[aj
], dx
);
642 rvec_sub(prime
[ai
], prime
[aj
], dx
);
645 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n", ai
+ 1,
646 1.0 / invmass
[ai
], aj
+ 1, 1.0 / invmass
[aj
], d
, dp
, ip
[ia
[0]].constr
.dA
);
648 case ConstraintVariable::Velocities
:
649 rvec_sub(v
[ai
], v
[aj
], dv
);
651 rvec_sub(prime
[ai
], prime
[aj
], dv
);
653 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n", ai
+ 1,
654 1.0 / invmass
[ai
], aj
+ 1, 1.0 / invmass
[aj
], d
, dp
, 0.);
656 default: gmx_incons("Unknown constraint quantity for SHAKE");
662 static bool bshakef(FILE* log
,
664 const real invmass
[],
665 const InteractionDefinitions
& idef
,
666 const t_inputrec
& ir
,
667 ArrayRef
<const RVec
> x_s
,
668 ArrayRef
<RVec
> prime
,
678 ConstraintVariable econq
)
681 int i
, n0
, ncon
, blen
, type
, ll
;
682 int tnit
= 0, trij
= 0;
684 ncon
= idef
.il
[F_CONSTR
].size() / 3;
686 for (ll
= 0; ll
< ncon
; ll
++)
688 shaked
->scaled_lagrange_multiplier
[ll
] = 0;
691 // TODO Rewrite this block so that it is obvious that i, iatoms
692 // and lam are all iteration variables. Is this easier if the
693 // sblock data structure is organized differently?
694 const int* iatoms
= &(idef
.il
[F_CONSTR
].iatoms
[shaked
->sblock
[0]]);
695 ArrayRef
<real
> lam
= shaked
->scaled_lagrange_multiplier
;
696 for (i
= 0; (i
< shaked
->numShakeBlocks());)
698 blen
= (shaked
->sblock
[i
+ 1] - shaked
->sblock
[i
]);
700 n0
= vec_shakef(log
, shaked
, invmass
, blen
, idef
.iparams
, iatoms
, ir
.shake_tol
, x_s
, prime
,
701 pbc
, shaked
->omega
, ir
.efep
!= efepNO
, lambda
, lam
, invdt
, v
, bCalcVir
,
706 if (bDumpOnError
&& log
)
709 check_cons(log
, blen
, x_s
, prime
, v
, pbc
, idef
.iparams
, iatoms
, invmass
, econq
);
716 iatoms
+= 3 * blen
; /* Increment pointer! */
717 lam
= lam
.subArray(blen
, lam
.ssize() - blen
);
720 /* only for position part? */
721 if (econq
== ConstraintVariable::Positions
)
723 if (ir
.efep
!= efepNO
)
725 ArrayRef
<const t_iparams
> iparams
= idef
.iparams
;
727 /* TODO This should probably use invdt, so that sd integrator scaling works properly */
728 dt_2
= 1 / gmx::square(ir
.delta_t
);
730 for (ll
= 0; ll
< ncon
; ll
++)
732 type
= idef
.il
[F_CONSTR
].iatoms
[3 * ll
];
734 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
735 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll
736 (eta_ll is the estimate of the Langrangian, definition on page 336 of Ryckaert et
737 al 1977), so the pre-factors are already present. */
738 const real bondA
= iparams
[type
].constr
.dA
;
739 const real bondB
= iparams
[type
].constr
.dB
;
740 dvdl
+= shaked
->scaled_lagrange_multiplier
[ll
] * dt_2
* (bondB
- bondA
);
747 if (tnit
> shaked
->gamma
)
749 shaked
->delta
*= -0.5;
751 shaked
->omega
+= shaked
->delta
;
752 shaked
->gamma
= tnit
;
754 inc_nrnb(nrnb
, eNR_SHAKE
, tnit
);
755 inc_nrnb(nrnb
, eNR_SHAKE_RIJ
, trij
);
758 inc_nrnb(nrnb
, eNR_CONSTR_V
, trij
* 2);
762 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, trij
);
768 bool constrain_shake(FILE* log
,
770 const real invmass
[],
771 const InteractionDefinitions
& idef
,
772 const t_inputrec
& ir
,
773 ArrayRef
<const RVec
> x_s
,
774 ArrayRef
<RVec
> xprime
,
775 ArrayRef
<RVec
> vprime
,
785 ConstraintVariable econq
)
787 if (shaked
->numShakeBlocks() == 0)
794 case (ConstraintVariable::Positions
):
795 bOK
= bshakef(log
, shaked
, invmass
, idef
, ir
, x_s
, xprime
, pbc
, nrnb
, lambda
, dvdlambda
,
796 invdt
, v
, bCalcVir
, vir_r_m_dr
, bDumpOnError
, econq
);
798 case (ConstraintVariable::Velocities
):
799 bOK
= bshakef(log
, shaked
, invmass
, idef
, ir
, x_s
, vprime
, pbc
, nrnb
, lambda
, dvdlambda
,
800 invdt
, {}, bCalcVir
, vir_r_m_dr
, bDumpOnError
, econq
);
804 "Internal error, SHAKE called for constraining something else than "