2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief This file defines low-level functions necessary for
38 * computing energies and forces for position restraints.
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_listed-forces
47 #include "position-restraints.h"
53 #include "gromacs/gmxlib/nrnb.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdtypes/forcerec.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/mdtypes/md_enums.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/timing/wallcycle.h"
61 #include "gromacs/topology/idef.h"
62 #include "gromacs/utility/basedefinitions.h"
63 #include "gromacs/utility/fatalerror.h"
70 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
72 void posres_dx(const rvec x
, const rvec pos0A
, const rvec pos0B
,
73 const rvec comA_sc
, const rvec comB_sc
,
75 const t_pbc
*pbc
, int refcoord_scaling
, int npbcdim
,
76 rvec dx
, rvec rdist
, rvec dpdl
)
79 real posA
, posB
, L1
, ref
= 0.;
84 for (m
= 0; m
< DIM
; m
++)
90 switch (refcoord_scaling
)
94 rdist
[m
] = L1
*posA
+ lambda
*posB
;
95 dpdl
[m
] = posB
- posA
;
98 /* Box relative coordinates are stored for dimensions with pbc */
99 posA
*= pbc
->box
[m
][m
];
100 posB
*= pbc
->box
[m
][m
];
101 assert(npbcdim
<= DIM
);
102 for (d
= m
+1; d
< npbcdim
; d
++)
104 posA
+= pos0A
[d
]*pbc
->box
[d
][m
];
105 posB
+= pos0B
[d
]*pbc
->box
[d
][m
];
107 ref
= L1
*posA
+ lambda
*posB
;
109 dpdl
[m
] = posB
- posA
;
112 ref
= L1
*comA_sc
[m
] + lambda
*comB_sc
[m
];
113 rdist
[m
] = L1
*posA
+ lambda
*posB
;
114 dpdl
[m
] = comB_sc
[m
] - comA_sc
[m
] + posB
- posA
;
117 gmx_fatal(FARGS
, "No such scaling method implemented");
122 ref
= L1
*posA
+ lambda
*posB
;
124 dpdl
[m
] = posB
- posA
;
127 /* We do pbc_dx with ref+rdist,
128 * since with only ref we can be up to half a box vector wrong.
130 pos
[m
] = ref
+ rdist
[m
];
135 pbc_dx(pbc
, x
, pos
, dx
);
139 rvec_sub(x
, pos
, dx
);
143 /*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
144 * Returns the flat-bottom potential. */
145 real
do_fbposres_cylinder(int fbdim
, rvec fm
, rvec dx
, real rfb
, real kk
, gmx_bool bInvert
)
148 real dr
, dr2
, invdr
, v
, rfb2
;
151 rfb2
= gmx::square(rfb
);
154 for (d
= 0; d
< DIM
; d
++)
158 dr2
+= gmx::square(dx
[d
]);
163 ( (dr2
> rfb2
&& bInvert
== FALSE
) || (dr2
< rfb2
&& bInvert
== TRUE
) )
168 v
= 0.5*kk
*gmx::square(dr
- rfb
);
169 for (d
= 0; d
< DIM
; d
++)
173 fm
[d
] = -kk
*(dr
-rfb
)*dx
[d
]*invdr
; /* Force pointing to the center */
181 /*! \brief Adds forces of flat-bottomed positions restraints to f[]
182 * and fixes vir_diag.
184 * Returns the flat-bottomed potential. Same PBC treatment as in
185 * normal position restraints */
186 real
fbposres(int nbonds
,
187 const t_iatom forceatoms
[], const t_iparams forceparams
[],
188 const rvec x
[], rvec f
[], rvec vir_diag
,
190 int refcoord_scaling
, int ePBC
, rvec com
)
191 /* compute flat-bottomed positions restraints */
193 int i
, ai
, m
, d
, type
, npbcdim
= 0, fbdim
;
196 real dr
, dr2
, rfb
, rfb2
, fact
;
197 rvec com_sc
, rdist
, dx
, dpdl
, fm
;
200 npbcdim
= ePBC2npbcdim(ePBC
);
202 if (refcoord_scaling
== erscCOM
)
205 for (m
= 0; m
< npbcdim
; m
++)
207 assert(npbcdim
<= DIM
);
208 for (d
= m
; d
< npbcdim
; d
++)
210 com_sc
[m
] += com
[d
]*pbc
->box
[d
][m
];
216 for (i
= 0; (i
< nbonds
); )
218 type
= forceatoms
[i
++];
219 ai
= forceatoms
[i
++];
220 pr
= &forceparams
[type
];
222 /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
223 posres_dx(x
[ai
], forceparams
[type
].fbposres
.pos0
, forceparams
[type
].fbposres
.pos0
,
225 pbc
, refcoord_scaling
, npbcdim
,
232 rfb
= pr
->fbposres
.r
;
233 rfb2
= gmx::square(rfb
);
235 /* with rfb<0, push particle out of the sphere/cylinder/layer */
243 switch (pr
->fbposres
.geom
)
245 case efbposresSPHERE
:
246 /* spherical flat-bottom posres */
249 ( (dr2
> rfb2
&& bInvert
== FALSE
) || (dr2
< rfb2
&& bInvert
== TRUE
) )
253 v
= 0.5*kk
*gmx::square(dr
- rfb
);
254 fact
= -kk
*(dr
-rfb
)/dr
; /* Force pointing to the center pos0 */
258 case efbposresCYLINDERX
:
259 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
261 v
= do_fbposres_cylinder(fbdim
, fm
, dx
, rfb
, kk
, bInvert
);
263 case efbposresCYLINDERY
:
264 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
266 v
= do_fbposres_cylinder(fbdim
, fm
, dx
, rfb
, kk
, bInvert
);
268 case efbposresCYLINDER
:
269 /* equivalent to efbposresCYLINDERZ for backwards compatibility */
270 case efbposresCYLINDERZ
:
271 /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
273 v
= do_fbposres_cylinder(fbdim
, fm
, dx
, rfb
, kk
, bInvert
);
275 case efbposresX
: /* fbdim=XX */
276 case efbposresY
: /* fbdim=YY */
277 case efbposresZ
: /* fbdim=ZZ */
278 /* 1D flat-bottom potential */
279 fbdim
= pr
->fbposres
.geom
- efbposresX
;
281 if ( ( dr
> rfb
&& bInvert
== FALSE
) || ( 0 < dr
&& dr
< rfb
&& bInvert
== TRUE
) )
283 v
= 0.5*kk
*gmx::square(dr
- rfb
);
284 fm
[fbdim
] = -kk
*(dr
- rfb
);
286 else if ( (dr
< (-rfb
) && bInvert
== FALSE
) || ( (-rfb
) < dr
&& dr
< 0 && bInvert
== TRUE
))
288 v
= 0.5*kk
*gmx::square(dr
+ rfb
);
289 fm
[fbdim
] = -kk
*(dr
+ rfb
);
296 for (m
= 0; (m
< DIM
); m
++)
299 /* Here we correct for the pbc_dx which included rdist */
300 vir_diag
[m
] -= 0.5*(dx
[m
] + rdist
[m
])*fm
[m
];
308 /*! \brief Compute energies and forces for position restraints
310 * Note that position restraints require a different pbc treatment
311 * from other bondeds */
312 real
posres(int nbonds
,
313 const t_iatom forceatoms
[], const t_iparams forceparams
[],
314 const rvec x
[], rvec f
[], rvec vir_diag
,
315 const struct t_pbc
*pbc
,
316 real lambda
, real
*dvdlambda
,
317 int refcoord_scaling
, int ePBC
, rvec comA
, rvec comB
)
319 int i
, ai
, m
, d
, type
, npbcdim
= 0;
323 rvec comA_sc
, comB_sc
, rdist
, dpdl
, dx
;
324 gmx_bool bForceValid
= TRUE
;
326 if ((f
== nullptr) || (vir_diag
== nullptr)) /* should both be null together! */
331 npbcdim
= ePBC2npbcdim(ePBC
);
333 if (refcoord_scaling
== erscCOM
)
337 for (m
= 0; m
< npbcdim
; m
++)
339 assert(npbcdim
<= DIM
);
340 for (d
= m
; d
< npbcdim
; d
++)
342 comA_sc
[m
] += comA
[d
]*pbc
->box
[d
][m
];
343 comB_sc
[m
] += comB
[d
]*pbc
->box
[d
][m
];
351 for (i
= 0; (i
< nbonds
); )
353 type
= forceatoms
[i
++];
354 ai
= forceatoms
[i
++];
355 pr
= &forceparams
[type
];
357 /* return dx, rdist, and dpdl */
358 posres_dx(x
[ai
], forceparams
[type
].posres
.pos0A
, forceparams
[type
].posres
.pos0B
,
359 comA_sc
, comB_sc
, lambda
,
360 pbc
, refcoord_scaling
, npbcdim
,
363 for (m
= 0; (m
< DIM
); m
++)
365 kk
= L1
*pr
->posres
.fcA
[m
] + lambda
*pr
->posres
.fcB
[m
];
367 vtot
+= 0.5*kk
*dx
[m
]*dx
[m
];
369 0.5*(pr
->posres
.fcB
[m
] - pr
->posres
.fcA
[m
])*dx
[m
]*dx
[m
]
372 /* Here we correct for the pbc_dx which included rdist */
376 vir_diag
[m
] -= 0.5*(dx
[m
] + rdist
[m
])*fm
;
387 posres_wrapper(t_nrnb
*nrnb
,
389 const struct t_pbc
*pbc
,
391 gmx_enerdata_t
*enerd
,
398 v
= posres(idef
->il
[F_POSRES
].nr
, idef
->il
[F_POSRES
].iatoms
,
399 idef
->iparams_posres
,
400 x
, as_rvec_array(fr
->f_novirsum
->data()), fr
->vir_diag_posres
,
401 fr
->ePBC
== epbcNONE
? nullptr : pbc
,
402 lambda
[efptRESTRAINT
], &dvdl
,
403 fr
->rc_scaling
, fr
->ePBC
, fr
->posres_com
, fr
->posres_comB
);
404 enerd
->term
[F_POSRES
] += v
;
405 /* If just the force constant changes, the FEP term is linear,
406 * but if k changes, it is not.
408 enerd
->dvdl_nonlin
[efptRESTRAINT
] += dvdl
;
409 inc_nrnb(nrnb
, eNR_POSRES
, idef
->il
[F_POSRES
].nr
/2);
413 posres_wrapper_lambda(struct gmx_wallcycle
*wcycle
,
414 const t_lambda
*fepvals
,
416 const struct t_pbc
*pbc
,
418 gmx_enerdata_t
*enerd
,
425 if (0 == idef
->il
[F_POSRES
].nr
)
430 wallcycle_sub_start_nocount(wcycle
, ewcsRESTRAINTS
);
431 for (i
= 0; i
< enerd
->n_lambda
; i
++)
433 real dvdl_dum
= 0, lambda_dum
;
435 lambda_dum
= (i
== 0 ? lambda
[efptRESTRAINT
] : fepvals
->all_lambda
[efptRESTRAINT
][i
-1]);
436 v
= posres(idef
->il
[F_POSRES
].nr
, idef
->il
[F_POSRES
].iatoms
,
437 idef
->iparams_posres
,
439 fr
->ePBC
== epbcNONE
? nullptr : pbc
, lambda_dum
, &dvdl_dum
,
440 fr
->rc_scaling
, fr
->ePBC
, fr
->posres_com
, fr
->posres_comB
);
441 enerd
->enerpart_lambda
[i
] += v
;
443 wallcycle_sub_stop(wcycle
, ewcsRESTRAINTS
);
446 /*! \brief Helper function that wraps calls to fbposres for
447 free-energy perturbation */
448 void fbposres_wrapper(t_nrnb
*nrnb
,
450 const struct t_pbc
*pbc
,
452 gmx_enerdata_t
*enerd
,
457 v
= fbposres(idef
->il
[F_FBPOSRES
].nr
, idef
->il
[F_FBPOSRES
].iatoms
,
458 idef
->iparams_fbposres
,
459 x
, as_rvec_array(fr
->f_novirsum
->data()), fr
->vir_diag_posres
,
460 fr
->ePBC
== epbcNONE
? nullptr : pbc
,
461 fr
->rc_scaling
, fr
->ePBC
, fr
->posres_com
);
462 enerd
->term
[F_FBPOSRES
] += v
;
463 inc_nrnb(nrnb
, eNR_FBPOSRES
, idef
->il
[F_FBPOSRES
].nr
/2);