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,2016,2017, 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.
47 #include "gromacs/domdec/domdec.h"
48 #include "gromacs/domdec/domdec_struct.h"
49 #include "gromacs/ewald/ewald.h"
50 #include "gromacs/ewald/long-range-correction.h"
51 #include "gromacs/ewald/pme.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/gmxlib/nrnb.h"
54 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
55 #include "gromacs/listed-forces/listed-forces.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vecdump.h"
58 #include "gromacs/mdlib/forcerec-threading.h"
59 #include "gromacs/mdlib/genborn.h"
60 #include "gromacs/mdlib/mdrun.h"
61 #include "gromacs/mdlib/ns.h"
62 #include "gromacs/mdlib/qmmm.h"
63 #include "gromacs/mdtypes/commrec.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/pbcutil/mshift.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/timing/wallcycle.h"
70 #include "gromacs/utility/cstringutil.h"
71 #include "gromacs/utility/exceptions.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/smalloc.h"
88 if (!fr
->ns
->nblist_initialized
)
90 init_neighbor_list(fp
, fr
, md
->homenr
);
93 nsearch
= search_neighbours(fp
, fr
, box
, top
, groups
, cr
, nrnb
, md
,
97 fprintf(debug
, "nsearch = %d\n", nsearch
);
100 /* Check whether we have to do dynamic load balancing */
101 /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
102 count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
103 &(top->idef),opts->ngener);
105 if (fr
->ns
->dump_nl
> 0)
107 dump_nblist(fp
, cr
, fr
, fr
->ns
->dump_nl
);
111 static void reduce_thread_energies(tensor vir_q
, tensor vir_lj
,
112 real
*Vcorr_q
, real
*Vcorr_lj
,
113 real
*dvdl_q
, real
*dvdl_lj
,
115 ewald_corr_thread_t
*ewc_t
)
119 for (t
= 1; t
< nthreads
; t
++)
121 *Vcorr_q
+= ewc_t
[t
].Vcorr_q
;
122 *Vcorr_lj
+= ewc_t
[t
].Vcorr_lj
;
123 *dvdl_q
+= ewc_t
[t
].dvdl
[efptCOUL
];
124 *dvdl_lj
+= ewc_t
[t
].dvdl
[efptVDW
];
125 m_add(vir_q
, ewc_t
[t
].vir_q
, vir_q
);
126 m_add(vir_lj
, ewc_t
[t
].vir_lj
, vir_lj
);
130 void do_force_lowlevel(t_forcerec
*fr
, t_inputrec
*ir
,
131 t_idef
*idef
, t_commrec
*cr
,
132 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
134 rvec x
[], history_t
*hist
,
136 gmx_enerdata_t
*enerd
,
157 real dvdl_dum
[efptNR
], dvdl_nb
[efptNR
];
160 double t0
= 0.0, t1
, t2
, t3
; /* time measurement for coarse load balancing */
163 set_pbc(&pbc
, fr
->ePBC
, box
);
165 /* reset free energy components */
166 for (i
= 0; i
< efptNR
; i
++)
173 for (i
= 0; (i
< DIM
); i
++)
175 box_size
[i
] = box
[i
][i
];
178 /* do QMMM first if requested */
181 enerd
->term
[F_EQM
] = calculate_QMMM(cr
, x
, f
, fr
);
184 /* Call the short range functions all in one go. */
187 /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
188 #define TAKETIME FALSE
191 MPI_Barrier(cr
->mpi_comm_mygroup
);
198 /* foreign lambda component for walls */
199 real dvdl_walls
= do_walls(ir
, fr
, box
, md
, x
, f
, lambda
[efptVDW
],
200 enerd
->grpp
.ener
[egLJSR
], nrnb
);
201 enerd
->dvdl_lin
[efptVDW
] += dvdl_walls
;
204 /* If doing GB, reset dvda and calculate the Born radii */
205 if (ir
->implicit_solvent
)
207 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
209 for (i
= 0; i
< born
->nr
; i
++)
216 calc_gb_rad(cr
, fr
, ir
, top
, x
, fr
->gblist
, born
, md
, nrnb
);
219 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
223 /* We only do non-bonded calculation with group scheme here, the verlet
224 * calls are done from do_force_cutsVERLET(). */
225 if (fr
->cutoff_scheme
== ecutsGROUP
&& (flags
& GMX_FORCE_NONBONDED
))
228 /* Add short-range interactions */
229 donb_flags
|= GMX_NONBONDED_DO_SR
;
231 /* Currently all group scheme kernels always calculate (shift-)forces */
232 if (flags
& GMX_FORCE_FORCES
)
234 donb_flags
|= GMX_NONBONDED_DO_FORCE
;
236 if (flags
& GMX_FORCE_VIRIAL
)
238 donb_flags
|= GMX_NONBONDED_DO_SHIFTFORCE
;
240 if (flags
& GMX_FORCE_ENERGY
)
242 donb_flags
|= GMX_NONBONDED_DO_POTENTIAL
;
245 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
246 do_nonbonded(fr
, x
, f
, md
, excl
,
248 lambda
, dvdl_nb
, -1, -1, donb_flags
);
250 /* If we do foreign lambda and we have soft-core interactions
251 * we have to recalculate the (non-linear) energies contributions.
253 if (fepvals
->n_lambda
> 0 && (flags
& GMX_FORCE_DHDL
) && fepvals
->sc_alpha
!= 0)
255 for (i
= 0; i
< enerd
->n_lambda
; i
++)
259 for (j
= 0; j
< efptNR
; j
++)
261 lam_i
[j
] = (i
== 0 ? lambda
[j
] : fepvals
->all_lambda
[j
][i
-1]);
263 reset_foreign_enerdata(enerd
);
264 do_nonbonded(fr
, x
, f
, md
, excl
,
265 &(enerd
->foreign_grpp
), nrnb
,
266 lam_i
, dvdl_dum
, -1, -1,
267 (donb_flags
& ~GMX_NONBONDED_DO_FORCE
) | GMX_NONBONDED_DO_FOREIGNLAMBDA
);
268 sum_epot(&(enerd
->foreign_grpp
), enerd
->foreign_term
);
269 enerd
->enerpart_lambda
[i
] += enerd
->foreign_term
[F_EPOT
];
272 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
276 /* If we are doing GB, calculate bonded forces and apply corrections
277 * to the solvation forces */
278 /* MRS: Eventually, many need to include free energy contribution here! */
279 if (ir
->implicit_solvent
)
281 wallcycle_sub_start(wcycle
, ewcsLISTED
);
282 calc_gb_forces(cr
, md
, born
, top
, x
, f
, fr
, idef
,
283 ir
->gb_algorithm
, ir
->sa_algorithm
, nrnb
, &pbc
, graph
, enerd
);
284 wallcycle_sub_stop(wcycle
, ewcsLISTED
);
295 if (fepvals
->sc_alpha
!= 0)
297 enerd
->dvdl_nonlin
[efptVDW
] += dvdl_nb
[efptVDW
];
301 enerd
->dvdl_lin
[efptVDW
] += dvdl_nb
[efptVDW
];
304 if (fepvals
->sc_alpha
!= 0)
306 /* even though coulomb part is linear, we already added it, beacuse we
307 need to go through the vdw calculation anyway */
309 enerd
->dvdl_nonlin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
313 enerd
->dvdl_lin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
318 pr_rvecs(debug
, 0, "fshift after SR", fr
->fshift
, SHIFTS
);
321 /* Shift the coordinates. Must be done before listed forces and PPPM,
322 * but is also necessary for SHAKE and update, therefore it can NOT
323 * go when no listed forces have to be evaluated.
325 * The shifting and PBC code is deliberately not timed, since with
326 * the Verlet scheme it only takes non-zero time with triclinic
327 * boxes, and even then the time is around a factor of 100 less
328 * than the next smallest counter.
332 /* Here sometimes we would not need to shift with NBFonly,
333 * but we do so anyhow for consistency of the returned coordinates.
337 shift_self(graph
, box
, x
);
340 inc_nrnb(nrnb
, eNR_SHIFTX
, 2*graph
->nnodes
);
344 inc_nrnb(nrnb
, eNR_SHIFTX
, graph
->nnodes
);
347 /* Check whether we need to do listed interactions or correct for exclusions */
349 ((flags
& GMX_FORCE_LISTED
)
350 || EEL_RF(fr
->eeltype
) || EEL_FULL(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
)))
352 /* TODO There are no electrostatics methods that require this
353 transformation, when using the Verlet scheme, so update the
354 above conditional. */
355 /* Since all atoms are in the rectangular or triclinic unit-cell,
356 * only single box vector shifts (2 in x) are required.
358 set_pbc_dd(&pbc
, fr
->ePBC
, DOMAINDECOMP(cr
) ? cr
->dd
->nc
: nullptr,
362 do_force_listed(wcycle
, box
, ir
->fepvals
, cr
,
363 idef
, (const rvec
*) x
, hist
, f
, fr
,
364 &pbc
, graph
, enerd
, nrnb
, lambda
, md
, fcd
,
365 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: nullptr,
371 clear_mat(fr
->vir_el_recip
);
372 clear_mat(fr
->vir_lj_recip
);
374 /* Do long-range electrostatics and/or LJ-PME, including related short-range
377 if (EEL_FULL(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
))
380 real Vlr_q
= 0, Vlr_lj
= 0, Vcorr_q
= 0, Vcorr_lj
= 0;
381 real dvdl_long_range_q
= 0, dvdl_long_range_lj
= 0;
383 bSB
= (ir
->nwall
== 2);
387 svmul(ir
->wall_ewald_zfac
, boxs
[ZZ
], boxs
[ZZ
]);
388 box_size
[ZZ
] *= ir
->wall_ewald_zfac
;
391 if (EEL_PME_EWALD(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
))
393 real dvdl_long_range_correction_q
= 0;
394 real dvdl_long_range_correction_lj
= 0;
395 /* With the Verlet scheme exclusion forces are calculated
396 * in the non-bonded kernel.
398 /* The TPI molecule does not have exclusions with the rest
399 * of the system and no intra-molecular PME grid
400 * contributions will be calculated in
401 * gmx_pme_calc_energy.
403 if ((ir
->cutoff_scheme
== ecutsGROUP
&& fr
->n_tpi
== 0) ||
404 ir
->ewald_geometry
!= eewg3D
||
405 ir
->epsilon_surface
!= 0)
409 wallcycle_sub_start(wcycle
, ewcsEWALD_CORRECTION
);
413 gmx_fatal(FARGS
, "TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
416 nthreads
= fr
->nthread_ewc
;
417 #pragma omp parallel for num_threads(nthreads) schedule(static)
418 for (t
= 0; t
< nthreads
; t
++)
422 tensor
*vir_q
, *vir_lj
;
423 real
*Vcorrt_q
, *Vcorrt_lj
, *dvdlt_q
, *dvdlt_lj
;
426 vir_q
= &fr
->vir_el_recip
;
427 vir_lj
= &fr
->vir_lj_recip
;
429 Vcorrt_lj
= &Vcorr_lj
;
430 dvdlt_q
= &dvdl_long_range_correction_q
;
431 dvdlt_lj
= &dvdl_long_range_correction_lj
;
435 vir_q
= &fr
->ewc_t
[t
].vir_q
;
436 vir_lj
= &fr
->ewc_t
[t
].vir_lj
;
437 Vcorrt_q
= &fr
->ewc_t
[t
].Vcorr_q
;
438 Vcorrt_lj
= &fr
->ewc_t
[t
].Vcorr_lj
;
439 dvdlt_q
= &fr
->ewc_t
[t
].dvdl
[efptCOUL
];
440 dvdlt_lj
= &fr
->ewc_t
[t
].dvdl
[efptVDW
];
447 /* Threading is only supported with the Verlet cut-off
448 * scheme and then only single particle forces (no
449 * exclusion forces) are calculated, so we can store
450 * the forces in the normal, single fr->f_novirsum array.
452 ewald_LRcorrection(md
->homenr
, cr
, nthreads
, t
, fr
,
453 md
->chargeA
, md
->chargeB
,
454 md
->sqrt_c6A
, md
->sqrt_c6B
,
455 md
->sigmaA
, md
->sigmaB
,
456 md
->sigma3A
, md
->sigma3B
,
457 md
->nChargePerturbed
|| md
->nTypePerturbed
,
458 ir
->cutoff_scheme
!= ecutsVERLET
,
459 excl
, x
, bSB
? boxs
: box
, mu_tot
,
462 as_rvec_array(fr
->f_novirsum
->data()),
465 lambda
[efptCOUL
], lambda
[efptVDW
],
468 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
472 reduce_thread_energies(fr
->vir_el_recip
, fr
->vir_lj_recip
,
474 &dvdl_long_range_correction_q
,
475 &dvdl_long_range_correction_lj
,
476 nthreads
, fr
->ewc_t
);
478 wallcycle_sub_stop(wcycle
, ewcsEWALD_CORRECTION
);
481 if (EEL_PME_EWALD(fr
->eeltype
) && fr
->n_tpi
== 0)
483 /* This is not in a subcounter because it takes a
484 negligible and constant-sized amount of time */
485 Vcorr_q
+= ewald_charge_correction(cr
, fr
, lambda
[efptCOUL
], box
,
486 &dvdl_long_range_correction_q
,
490 enerd
->dvdl_lin
[efptCOUL
] += dvdl_long_range_correction_q
;
491 enerd
->dvdl_lin
[efptVDW
] += dvdl_long_range_correction_lj
;
493 if ((EEL_PME(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
)) && (cr
->duty
& DUTY_PME
))
495 /* Do reciprocal PME for Coulomb and/or LJ. */
496 assert(fr
->n_tpi
>= 0);
497 if (fr
->n_tpi
== 0 || (flags
& GMX_FORCE_STATECHANGED
))
499 pme_flags
= GMX_PME_SPREAD
| GMX_PME_SOLVE
;
501 if (flags
& GMX_FORCE_FORCES
)
503 pme_flags
|= GMX_PME_CALC_F
;
505 if (flags
& GMX_FORCE_VIRIAL
)
507 pme_flags
|= GMX_PME_CALC_ENER_VIR
;
511 /* We don't calculate f, but we do want the potential */
512 pme_flags
|= GMX_PME_CALC_POT
;
514 wallcycle_start(wcycle
, ewcPMEMESH
);
515 status
= gmx_pme_do(fr
->pmedata
,
516 0, md
->homenr
- fr
->n_tpi
,
518 as_rvec_array(fr
->f_novirsum
->data()),
519 md
->chargeA
, md
->chargeB
,
520 md
->sqrt_c6A
, md
->sqrt_c6B
,
521 md
->sigmaA
, md
->sigmaB
,
522 bSB
? boxs
: box
, cr
,
523 DOMAINDECOMP(cr
) ? dd_pme_maxshift_x(cr
->dd
) : 0,
524 DOMAINDECOMP(cr
) ? dd_pme_maxshift_y(cr
->dd
) : 0,
526 fr
->vir_el_recip
, fr
->vir_lj_recip
,
528 lambda
[efptCOUL
], lambda
[efptVDW
],
529 &dvdl_long_range_q
, &dvdl_long_range_lj
, pme_flags
);
530 *cycles_pme
= wallcycle_stop(wcycle
, ewcPMEMESH
);
533 gmx_fatal(FARGS
, "Error %d in reciprocal PME routine", status
);
535 /* We should try to do as little computation after
536 * this as possible, because parallel PME synchronizes
537 * the nodes, so we want all load imbalance of the
538 * rest of the force calculation to be before the PME
539 * call. DD load balancing is done on the whole time
540 * of the force call (without PME).
545 if (EVDW_PME(ir
->vdwtype
))
548 gmx_fatal(FARGS
, "Test particle insertion not implemented with LJ-PME");
550 /* Determine the PME grid energy of the test molecule
551 * with the PME grid potential of the other charges.
553 gmx_pme_calc_energy(fr
->pmedata
, fr
->n_tpi
,
554 x
+ md
->homenr
- fr
->n_tpi
,
555 md
->chargeA
+ md
->homenr
- fr
->n_tpi
,
561 if (!EEL_PME(fr
->eeltype
) && EEL_PME_EWALD(fr
->eeltype
))
563 Vlr_q
= do_ewald(ir
, x
, as_rvec_array(fr
->f_novirsum
->data()),
564 md
->chargeA
, md
->chargeB
,
565 box_size
, cr
, md
->homenr
,
566 fr
->vir_el_recip
, fr
->ewaldcoeff_q
,
567 lambda
[efptCOUL
], &dvdl_long_range_q
, fr
->ewald_table
);
570 /* Note that with separate PME nodes we get the real energies later */
571 enerd
->dvdl_lin
[efptCOUL
] += dvdl_long_range_q
;
572 enerd
->dvdl_lin
[efptVDW
] += dvdl_long_range_lj
;
573 enerd
->term
[F_COUL_RECIP
] = Vlr_q
+ Vcorr_q
;
574 enerd
->term
[F_LJ_RECIP
] = Vlr_lj
+ Vcorr_lj
;
577 fprintf(debug
, "Vlr_q = %g, Vcorr_q = %g, Vlr_corr_q = %g\n",
578 Vlr_q
, Vcorr_q
, enerd
->term
[F_COUL_RECIP
]);
579 pr_rvecs(debug
, 0, "vir_el_recip after corr", fr
->vir_el_recip
, DIM
);
580 pr_rvecs(debug
, 0, "fshift after LR Corrections", fr
->fshift
, SHIFTS
);
581 fprintf(debug
, "Vlr_lj: %g, Vcorr_lj = %g, Vlr_corr_lj = %g\n",
582 Vlr_lj
, Vcorr_lj
, enerd
->term
[F_LJ_RECIP
]);
583 pr_rvecs(debug
, 0, "vir_lj_recip after corr", fr
->vir_lj_recip
, DIM
);
588 /* Is there a reaction-field exclusion correction needed?
589 * With the Verlet scheme, exclusion forces are calculated
590 * in the non-bonded kernel.
592 if (ir
->cutoff_scheme
!= ecutsVERLET
&& EEL_RF(fr
->eeltype
))
594 real dvdl_rf_excl
= 0;
595 enerd
->term
[F_RF_EXCL
] =
596 RF_excl_correction(fr
, graph
, md
, excl
, x
, f
,
597 fr
->fshift
, &pbc
, lambda
[efptCOUL
], &dvdl_rf_excl
);
599 enerd
->dvdl_lin
[efptCOUL
] += dvdl_rf_excl
;
606 print_nrnb(debug
, nrnb
);
613 MPI_Barrier(cr
->mpi_comm_mygroup
);
616 if (fr
->timesteps
== 11)
619 fprintf(stderr
, "* PP load balancing info: rank %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
620 cr
->nodeid
, gmx_step_str(fr
->timesteps
, buf
),
621 100*fr
->t_wait
/(fr
->t_wait
+fr
->t_fnbf
),
622 (fr
->t_fnbf
+fr
->t_wait
)/fr
->t_fnbf
);
630 pr_rvecs(debug
, 0, "fshift after bondeds", fr
->fshift
, SHIFTS
);
635 void init_enerdata(int ngener
, int n_lambda
, gmx_enerdata_t
*enerd
)
639 for (i
= 0; i
< F_NRE
; i
++)
642 enerd
->foreign_term
[i
] = 0;
646 for (i
= 0; i
< efptNR
; i
++)
648 enerd
->dvdl_lin
[i
] = 0;
649 enerd
->dvdl_nonlin
[i
] = 0;
655 fprintf(debug
, "Creating %d sized group matrix for energies\n", n2
);
657 enerd
->grpp
.nener
= n2
;
658 enerd
->foreign_grpp
.nener
= n2
;
659 for (i
= 0; (i
< egNR
); i
++)
661 snew(enerd
->grpp
.ener
[i
], n2
);
662 snew(enerd
->foreign_grpp
.ener
[i
], n2
);
667 enerd
->n_lambda
= 1 + n_lambda
;
668 snew(enerd
->enerpart_lambda
, enerd
->n_lambda
);
676 void destroy_enerdata(gmx_enerdata_t
*enerd
)
680 for (i
= 0; (i
< egNR
); i
++)
682 sfree(enerd
->grpp
.ener
[i
]);
685 for (i
= 0; (i
< egNR
); i
++)
687 sfree(enerd
->foreign_grpp
.ener
[i
]);
692 sfree(enerd
->enerpart_lambda
);
696 static real
sum_v(int n
, real v
[])
702 for (i
= 0; (i
< n
); i
++)
710 void sum_epot(gmx_grppairener_t
*grpp
, real
*epot
)
714 /* Accumulate energies */
715 epot
[F_COUL_SR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULSR
]);
716 epot
[F_LJ
] = sum_v(grpp
->nener
, grpp
->ener
[egLJSR
]);
717 epot
[F_LJ14
] = sum_v(grpp
->nener
, grpp
->ener
[egLJ14
]);
718 epot
[F_COUL14
] = sum_v(grpp
->nener
, grpp
->ener
[egCOUL14
]);
719 /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
720 epot
[F_GBPOL
] += sum_v(grpp
->nener
, grpp
->ener
[egGB
]);
722 /* lattice part of LR doesnt belong to any group
723 * and has been added earlier
725 epot
[F_BHAM
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMSR
]);
728 for (i
= 0; (i
< F_EPOT
); i
++)
730 if (i
!= F_DISRESVIOL
&& i
!= F_ORIRESDEV
)
732 epot
[F_EPOT
] += epot
[i
];
737 void sum_dhdl(gmx_enerdata_t
*enerd
, const std::vector
<real
> *lambda
, t_lambda
*fepvals
)
742 enerd
->dvdl_lin
[efptVDW
] += enerd
->term
[F_DVDL_VDW
]; /* include dispersion correction */
743 enerd
->term
[F_DVDL
] = 0.0;
744 for (i
= 0; i
< efptNR
; i
++)
746 if (fepvals
->separate_dvdl
[i
])
748 /* could this be done more readably/compactly? */
761 index
= F_DVDL_BONDED
;
763 case (efptRESTRAINT
):
764 index
= F_DVDL_RESTRAINT
;
770 enerd
->term
[index
] = enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
773 fprintf(debug
, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
774 efpt_names
[i
], i
, enerd
->term
[index
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
779 enerd
->term
[F_DVDL
] += enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
782 fprintf(debug
, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
783 efpt_names
[0], i
, enerd
->term
[F_DVDL
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
788 /* Notes on the foreign lambda free energy difference evaluation:
789 * Adding the potential and ekin terms that depend linearly on lambda
790 * as delta lam * dvdl to the energy differences is exact.
791 * For the constraints this is not exact, but we have no other option
792 * without literally changing the lengths and reevaluating the energies at each step.
793 * (try to remedy this post 4.6 - MRS)
795 if (fepvals
->separate_dvdl
[efptBONDED
])
797 enerd
->term
[F_DVDL_BONDED
] += enerd
->term
[F_DVDL_CONSTR
];
801 enerd
->term
[F_DVDL
] += enerd
->term
[F_DVDL_CONSTR
];
803 enerd
->term
[F_DVDL_CONSTR
] = 0;
805 for (i
= 0; i
< fepvals
->n_lambda
; i
++)
807 /* note we are iterating over fepvals here!
808 For the current lam, dlam = 0 automatically,
809 so we don't need to add anything to the
810 enerd->enerpart_lambda[0] */
812 /* we don't need to worry about dvdl_lin contributions to dE at
813 current lambda, because the contributions to the current
814 lambda are automatically zeroed */
816 for (j
= 0; j
< efptNR
; j
++)
818 /* Note that this loop is over all dhdl components, not just the separated ones */
819 dlam
= (fepvals
->all_lambda
[j
][i
] - (*lambda
)[j
]);
820 enerd
->enerpart_lambda
[i
+1] += dlam
*enerd
->dvdl_lin
[j
];
823 fprintf(debug
, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
824 fepvals
->all_lambda
[j
][i
], efpt_names
[j
],
825 (enerd
->enerpart_lambda
[i
+1] - enerd
->enerpart_lambda
[0]),
826 dlam
, enerd
->dvdl_lin
[j
]);
833 void reset_foreign_enerdata(gmx_enerdata_t
*enerd
)
837 /* First reset all foreign energy components. Foreign energies always called on
838 neighbor search steps */
839 for (i
= 0; (i
< egNR
); i
++)
841 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
843 enerd
->foreign_grpp
.ener
[i
][j
] = 0.0;
847 /* potential energy components */
848 for (i
= 0; (i
<= F_EPOT
); i
++)
850 enerd
->foreign_term
[i
] = 0.0;
854 void reset_enerdata(gmx_enerdata_t
*enerd
)
858 /* First reset all energy components. */
859 for (i
= 0; (i
< egNR
); i
++)
861 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
863 enerd
->grpp
.ener
[i
][j
] = 0.0;
866 for (i
= 0; i
< efptNR
; i
++)
868 enerd
->dvdl_lin
[i
] = 0.0;
869 enerd
->dvdl_nonlin
[i
] = 0.0;
872 /* Normal potential energy components */
873 for (i
= 0; (i
<= F_EPOT
); i
++)
875 enerd
->term
[i
] = 0.0;
877 enerd
->term
[F_DVDL
] = 0.0;
878 enerd
->term
[F_DVDL_COUL
] = 0.0;
879 enerd
->term
[F_DVDL_VDW
] = 0.0;
880 enerd
->term
[F_DVDL_BONDED
] = 0.0;
881 enerd
->term
[F_DVDL_RESTRAINT
] = 0.0;
882 enerd
->term
[F_DKDL
] = 0.0;
883 if (enerd
->n_lambda
> 0)
885 for (i
= 0; i
< enerd
->n_lambda
; i
++)
887 enerd
->enerpart_lambda
[i
] = 0.0;
890 /* reset foreign energy data - separate function since we also call it elsewhere */
891 reset_foreign_enerdata(enerd
);