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, 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.
41 #include "gromacs/gmxlib/nrnb.h"
42 #include "gromacs/math/functions.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/mdlib/constr.h"
45 #include "gromacs/mdtypes/inputrec.h"
46 #include "gromacs/mdtypes/md_enums.h"
47 #include "gromacs/utility/smalloc.h"
49 typedef struct gmx_shakedata
52 real
*half_of_reduced_mass
;
53 real
*distance_squared_tolerance
;
54 real
*constraint_distance_squared
;
62 gmx_shakedata_t
shake_init()
70 d
->half_of_reduced_mass
= NULL
;
71 d
->distance_squared_tolerance
= NULL
;
72 d
->constraint_distance_squared
= NULL
;
74 /* SOR initialization */
82 /*! \brief Inner kernel for SHAKE constraints
84 * Original implementation from R.C. van Schaik and W.F. van Gunsteren
85 * (ETH Zuerich, June 1992), adapted for GROMACS by David van der
86 * Spoel November 1992.
88 * The algorithm here is based section five of Ryckaert, Ciccotti and
89 * Berendsen, J Comp Phys, 23, 327, 1977.
91 * \param[in] iatom Mini-topology of triples of constraint type (unused in this
92 * function) and indices of the two atoms involved
93 * \param[in] ncon Number of constraints
94 * \param[out] nnit Number of iterations performed
95 * \param[in] maxnit Maximum number of iterations permitted
96 * \param[in] constraint_distance_squared The objective value for each constraint
97 * \param[inout] positions The initial (and final) values of the positions of all atoms
98 * \param[in] initial_displacements The initial displacements of each constraint
99 * \param[in] half_of_reduced_mass Half of the reduced mass for each constraint
100 * \param[in] omega SHAKE over-relaxation factor (set non-1.0 by
101 * using shake-sor=yes in the .mdp, but there is no documentation anywhere)
102 * \param[in] invmass Inverse mass of each atom
103 * \param[in] distance_squared_tolerance Multiplicative tolerance on the difference in the
104 * square of the constrained distance (see code)
105 * \param[out] scaled_lagrange_multiplier Scaled Lagrange multiplier for each constraint (-2 * eta from p. 336
106 * of the paper, divided by the constraint distance)
107 * \param[out] nerror Zero upon success, returns one more than the index of the
108 * problematic constraint if the input was malformed
110 * \todo Make SHAKE use better data structures, in particular for iatom. */
111 void cshake(const int iatom
[], int ncon
, int *nnit
, int maxnit
,
112 const real constraint_distance_squared
[], real positions
[],
113 const real initial_displacements
[], const real half_of_reduced_mass
[], real omega
,
114 const real invmass
[], const real distance_squared_tolerance
[],
115 real scaled_lagrange_multiplier
[], int *nerror
)
117 /* default should be increased! MRS 8/4/2009 */
118 const real mytol
= 1e-10;
120 int ll
, i
, j
, i3
, j3
, l3
;
121 int ix
, iy
, iz
, jx
, jy
, jz
;
123 real constraint_distance_squared_ll
;
124 real r_prime_squared
;
125 real scaled_lagrange_multiplier_ll
;
126 real r_prime_x
, r_prime_y
, r_prime_z
, diff
, im
, jm
;
127 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
128 int nit
, error
, nconv
;
131 // TODO nconv is used solely as a boolean, so we should write the
135 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
138 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
141 rijx
= initial_displacements
[l3
+XX
];
142 rijy
= initial_displacements
[l3
+YY
];
143 rijz
= initial_displacements
[l3
+ZZ
];
155 /* Compute r prime between atoms i and j, which is the
156 displacement *before* this update stage */
157 r_prime_x
= positions
[ix
]-positions
[jx
];
158 r_prime_y
= positions
[iy
]-positions
[jy
];
159 r_prime_z
= positions
[iz
]-positions
[jz
];
160 r_prime_squared
= (r_prime_x
* r_prime_x
+
161 r_prime_y
* r_prime_y
+
162 r_prime_z
* r_prime_z
);
163 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
164 diff
= constraint_distance_squared_ll
- r_prime_squared
;
166 /* iconvf is less than 1 when the error is smaller than a bound */
167 iconvf
= fabs(diff
) * distance_squared_tolerance
[ll
];
171 nconv
= static_cast<int>(iconvf
);
172 r_dot_r_prime
= (rijx
* r_prime_x
+
176 if (r_dot_r_prime
< constraint_distance_squared_ll
* mytol
)
182 /* The next line solves equation 5.6 (neglecting
183 the term in g^2), for g */
184 scaled_lagrange_multiplier_ll
= omega
*diff
*half_of_reduced_mass
[ll
]/r_dot_r_prime
;
185 scaled_lagrange_multiplier
[ll
] += scaled_lagrange_multiplier_ll
;
186 xh
= rijx
* scaled_lagrange_multiplier_ll
;
187 yh
= rijy
* scaled_lagrange_multiplier_ll
;
188 zh
= rijz
* scaled_lagrange_multiplier_ll
;
191 positions
[ix
] += xh
*im
;
192 positions
[iy
] += yh
*im
;
193 positions
[iz
] += zh
*im
;
194 positions
[jx
] -= xh
*jm
;
195 positions
[jy
] -= yh
*jm
;
196 positions
[jz
] -= zh
*jm
;
205 int vec_shakef(FILE *fplog
, gmx_shakedata_t shaked
,
206 real invmass
[], int ncon
,
207 t_iparams ip
[], t_iatom
*iatom
,
208 real tol
, rvec x
[], rvec prime
[], real omega
,
209 gmx_bool bFEP
, real lambda
, real scaled_lagrange_multiplier
[],
211 gmx_bool bCalcVir
, tensor vir_r_m_dr
, int econq
)
214 real
*half_of_reduced_mass
, *distance_squared_tolerance
, *constraint_distance_squared
;
216 int nit
= 0, ll
, i
, j
, d
, d2
, type
;
221 real constraint_distance
;
223 if (ncon
> shaked
->nalloc
)
225 shaked
->nalloc
= over_alloc_dd(ncon
);
226 srenew(shaked
->rij
, shaked
->nalloc
);
227 srenew(shaked
->half_of_reduced_mass
, shaked
->nalloc
);
228 srenew(shaked
->distance_squared_tolerance
, shaked
->nalloc
);
229 srenew(shaked
->constraint_distance_squared
, shaked
->nalloc
);
232 half_of_reduced_mass
= shaked
->half_of_reduced_mass
;
233 distance_squared_tolerance
= shaked
->distance_squared_tolerance
;
234 constraint_distance_squared
= shaked
->constraint_distance_squared
;
238 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
244 mm
= 2.0*(invmass
[i
]+invmass
[j
]);
245 rij
[ll
][XX
] = x
[i
][XX
]-x
[j
][XX
];
246 rij
[ll
][YY
] = x
[i
][YY
]-x
[j
][YY
];
247 rij
[ll
][ZZ
] = x
[i
][ZZ
]-x
[j
][ZZ
];
248 half_of_reduced_mass
[ll
] = 1.0/mm
;
251 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
255 constraint_distance
= ip
[type
].constr
.dA
;
257 constraint_distance_squared
[ll
] = gmx::square(constraint_distance
);
258 distance_squared_tolerance
[ll
] = 0.5/(constraint_distance_squared
[ll
]*tol
);
264 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
);
267 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
);
275 fprintf(fplog
, "Shake did not converge in %d steps\n", maxnit
);
277 fprintf(stderr
, "Shake did not converge in %d steps\n", maxnit
);
284 fprintf(fplog
, "Inner product between old and new vector <= 0.0!\n"
285 "constraint #%d atoms %d and %d\n",
286 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
288 fprintf(stderr
, "Inner product between old and new vector <= 0.0!\n"
289 "constraint #%d atoms %d and %d\n",
290 error
-1, iatom
[3*(error
-1)+1]+1, iatom
[3*(error
-1)+2]+1);
294 /* Constraint virial and correct the Lagrange multipliers for the length */
298 for (ll
= 0; (ll
< ncon
); ll
++, ia
+= 3)
304 if ((econq
== econqCoord
) && v
!= NULL
)
306 /* Correct the velocities */
307 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[i
]*invdt
;
308 for (d
= 0; d
< DIM
; d
++)
310 v
[ia
[1]][d
] += mm
*rij
[ll
][d
];
312 mm
= scaled_lagrange_multiplier
[ll
]*invmass
[j
]*invdt
;
313 for (d
= 0; d
< DIM
; d
++)
315 v
[ia
[2]][d
] -= mm
*rij
[ll
][d
];
320 /* constraint virial */
323 mm
= scaled_lagrange_multiplier
[ll
];
324 for (d
= 0; d
< DIM
; d
++)
327 for (d2
= 0; d2
< DIM
; d2
++)
329 vir_r_m_dr
[d
][d2
] -= tmp
*rij
[ll
][d2
];
335 /* cshake and crattle produce Lagrange multipliers scaled by
336 the reciprocal of the constraint length, so fix that */
339 constraint_distance
= L1
*ip
[type
].constr
.dA
+ lambda
*ip
[type
].constr
.dB
;
343 constraint_distance
= ip
[type
].constr
.dA
;
345 scaled_lagrange_multiplier
[ll
] *= constraint_distance
;
351 static void check_cons(FILE *log
, int nc
, rvec x
[], rvec prime
[], rvec v
[],
352 t_iparams ip
[], t_iatom
*iatom
,
353 real invmass
[], int econq
)
362 " i mi j mj before after should be\n");
364 for (i
= 0; (i
< nc
); i
++, ia
+= 3)
368 rvec_sub(x
[ai
], x
[aj
], dx
);
374 rvec_sub(prime
[ai
], prime
[aj
], dx
);
376 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
377 ai
+1, 1.0/invmass
[ai
],
378 aj
+1, 1.0/invmass
[aj
], d
, dp
, ip
[ia
[0]].constr
.dA
);
381 rvec_sub(v
[ai
], v
[aj
], dv
);
383 rvec_sub(prime
[ai
], prime
[aj
], dv
);
385 fprintf(log
, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
386 ai
+1, 1.0/invmass
[ai
],
387 aj
+1, 1.0/invmass
[aj
], d
, dp
, 0.);
393 gmx_bool
bshakef(FILE *log
, gmx_shakedata_t shaked
,
394 real invmass
[], int nblocks
, int sblock
[],
395 t_idef
*idef
, t_inputrec
*ir
, rvec x_s
[], rvec prime
[],
396 t_nrnb
*nrnb
, real
*scaled_lagrange_multiplier
, real lambda
, real
*dvdlambda
,
397 real invdt
, rvec
*v
, gmx_bool bCalcVir
, tensor vir_r_m_dr
,
398 gmx_bool bDumpOnError
, int econq
)
402 int i
, n0
, ncon
, blen
, type
, ll
;
403 int tnit
= 0, trij
= 0;
406 fprintf(log
, "nblocks=%d, sblock[0]=%d\n", nblocks
, sblock
[0]);
409 ncon
= idef
->il
[F_CONSTR
].nr
/3;
411 for (ll
= 0; ll
< ncon
; ll
++)
413 scaled_lagrange_multiplier
[ll
] = 0;
416 iatoms
= &(idef
->il
[F_CONSTR
].iatoms
[sblock
[0]]);
417 for (i
= 0; (i
< nblocks
); )
419 blen
= (sblock
[i
+1]-sblock
[i
]);
421 n0
= vec_shakef(log
, shaked
, invmass
, blen
, idef
->iparams
,
422 iatoms
, ir
->shake_tol
, x_s
, prime
, shaked
->omega
,
423 ir
->efep
!= efepNO
, lambda
, scaled_lagrange_multiplier
, invdt
, v
, bCalcVir
, vir_r_m_dr
,
427 check_cons(log
, blen
, x_s
, prime
, v
, idef
->iparams
, iatoms
, invmass
, econq
);
432 if (bDumpOnError
&& log
)
435 check_cons(log
, blen
, x_s
, prime
, v
, idef
->iparams
, iatoms
, invmass
, econq
);
442 iatoms
+= 3*blen
; /* Increment pointer! */
443 scaled_lagrange_multiplier
+= blen
;
446 /* only for position part? */
447 if (econq
== econqCoord
)
449 if (ir
->efep
!= efepNO
)
452 /* TODO This should probably use invdt, so that sd integrator scaling works properly */
453 dt_2
= 1/gmx::square(ir
->delta_t
);
455 for (ll
= 0; ll
< ncon
; ll
++)
457 type
= idef
->il
[F_CONSTR
].iatoms
[3*ll
];
459 /* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
460 /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
461 estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
462 so the pre-factors are already present. */
463 bondA
= idef
->iparams
[type
].constr
.dA
;
464 bondB
= idef
->iparams
[type
].constr
.dB
;
465 dvdl
+= scaled_lagrange_multiplier
[ll
] * dt_2
* (bondB
- bondA
);
471 fprintf(log
, "tnit: %5d omega: %10.5f\n", tnit
, omega
);
475 if (tnit
> shaked
->gamma
)
477 shaked
->delta
*= -0.5;
479 shaked
->omega
+= shaked
->delta
;
480 shaked
->gamma
= tnit
;
482 inc_nrnb(nrnb
, eNR_SHAKE
, tnit
);
483 inc_nrnb(nrnb
, eNR_SHAKE_RIJ
, trij
);
486 inc_nrnb(nrnb
, eNR_CONSTR_V
, trij
*2);
490 inc_nrnb(nrnb
, eNR_CONSTR_VIR
, trij
);
496 void crattle(int iatom
[], int ncon
, int *nnit
, int maxnit
,
497 real constraint_distance_squared
[], real vp
[], real rij
[], real m2
[], real omega
,
498 real invmass
[], real distance_squared_tolerance
[], real scaled_lagrange_multiplier
[],
499 int *nerror
, real invdt
)
502 * r.c. van schaik and w.f. van gunsteren
505 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
506 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
507 * second part of rattle algorithm
510 int ll
, i
, j
, i3
, j3
, l3
;
511 int ix
, iy
, iz
, jx
, jy
, jz
;
512 real constraint_distance_squared_ll
;
513 real vpijd
, vx
, vy
, vz
, acor
, fac
, im
, jm
;
514 real xh
, yh
, zh
, rijx
, rijy
, rijz
;
515 int nit
, error
, nconv
;
518 // TODO nconv is used solely as a boolean, so we should write the
522 for (nit
= 0; (nit
< maxnit
) && (nconv
!= 0) && (error
== 0); nit
++)
525 for (ll
= 0; (ll
< ncon
) && (error
== 0); ll
++)
545 vpijd
= vx
*rijx
+vy
*rijy
+vz
*rijz
;
546 constraint_distance_squared_ll
= constraint_distance_squared
[ll
];
548 /* iconv is zero when the error is smaller than a bound */
549 iconvf
= fabs(vpijd
)*(distance_squared_tolerance
[ll
]/invdt
);
553 nconv
= static_cast<int>(iconvf
);
554 fac
= omega
*2.0*m2
[ll
]/constraint_distance_squared_ll
;
556 scaled_lagrange_multiplier
[ll
] += acor
;