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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * 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.
52 #include "nonbonded.h"
67 #include "mpelogging.h"
68 #include "gmx_omp_nthreads.h"
82 gmx_grppairener_t
*grppener
,
84 gmx_bool bDoLongRangeNS
)
89 GMX_MPE_LOG(ev_ns_start
);
90 if (!fr
->ns
.nblist_initialized
)
92 init_neighbor_list(fp
, fr
, md
->homenr
);
100 nsearch
= search_neighbours(fp
, fr
, x
, box
, top
, groups
, cr
, nrnb
, md
,
101 lambda
, dvdlambda
, grppener
,
102 bFillGrid
, bDoLongRangeNS
, TRUE
);
105 fprintf(debug
, "nsearch = %d\n", nsearch
);
108 /* Check whether we have to do dynamic load balancing */
109 /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
110 count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
111 &(top->idef),opts->ngener);
113 if (fr
->ns
.dump_nl
> 0)
115 dump_nblist(fp
, cr
, fr
, fr
->ns
.dump_nl
);
118 GMX_MPE_LOG(ev_ns_finish
);
121 static void reduce_thread_forces(int n
, rvec
*f
,
124 int efpt_ind
, real
*dvdl
,
125 int nthreads
, f_thread_t
*f_t
)
129 /* This reduction can run over any number of threads */
130 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
131 for (i
= 0; i
< n
; i
++)
133 for (t
= 1; t
< nthreads
; t
++)
135 rvec_inc(f
[i
], f_t
[t
].f
[i
]);
138 for (t
= 1; t
< nthreads
; t
++)
140 *Vcorr
+= f_t
[t
].Vcorr
;
141 *dvdl
+= f_t
[t
].dvdl
[efpt_ind
];
142 m_add(vir
, f_t
[t
].vir
, vir
);
146 void do_force_lowlevel(FILE *fplog
, gmx_large_int_t step
,
147 t_forcerec
*fr
, t_inputrec
*ir
,
148 t_idef
*idef
, t_commrec
*cr
,
149 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
152 rvec x
[], history_t
*hist
,
155 gmx_enerdata_t
*enerd
,
173 gmx_bool bDoEpot
, bSepDVDL
, bSB
;
177 real Vsr
, Vlr
, Vcorr
= 0;
181 double clam_i
, vlam_i
;
182 real dvdl_dum
[efptNR
], dvdl
, dvdl_nb
[efptNR
], lam_i
[efptNR
];
186 double t0
= 0.0, t1
, t2
, t3
; /* time measurement for coarse load balancing */
189 #define PRINT_SEPDVDL(s, v, dvdlambda) if (bSepDVDL) {fprintf(fplog, sepdvdlformat, s, v, dvdlambda); }
191 GMX_MPE_LOG(ev_force_start
);
192 set_pbc(&pbc
, fr
->ePBC
, box
);
194 /* reset free energy components */
195 for (i
= 0; i
< efptNR
; i
++)
202 for (i
= 0; (i
< DIM
); i
++)
204 box_size
[i
] = box
[i
][i
];
207 bSepDVDL
= (fr
->bSepDVDL
&& do_per_step(step
, ir
->nstlog
));
210 /* do QMMM first if requested */
213 enerd
->term
[F_EQM
] = calculate_QMMM(cr
, x
, f
, fr
, md
);
218 fprintf(fplog
, "Step %s: non-bonded V and dVdl for node %d:\n",
219 gmx_step_str(step
, buf
), cr
->nodeid
);
222 /* Call the short range functions all in one go. */
223 GMX_MPE_LOG(ev_do_fnbf_start
);
226 /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
227 #define TAKETIME FALSE
230 MPI_Barrier(cr
->mpi_comm_mygroup
);
237 /* foreign lambda component for walls */
238 dvdl
= do_walls(ir
, fr
, box
, md
, x
, f
, lambda
[efptVDW
],
239 enerd
->grpp
.ener
[egLJSR
], nrnb
);
240 PRINT_SEPDVDL("Walls", 0.0, dvdl
);
241 enerd
->dvdl_lin
[efptVDW
] += dvdl
;
244 /* If doing GB, reset dvda and calculate the Born radii */
245 if (ir
->implicit_solvent
)
247 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
249 for (i
= 0; i
< born
->nr
; i
++)
256 calc_gb_rad(cr
, fr
, ir
, top
, atype
, x
, &(fr
->gblist
), born
, md
, nrnb
);
259 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
263 /* We only do non-bonded calculation with group scheme here, the verlet
264 * calls are done from do_force_cutsVERLET(). */
265 if (fr
->cutoff_scheme
== ecutsGROUP
&& (flags
& GMX_FORCE_NONBONDED
))
268 /* Add short-range interactions */
269 donb_flags
|= GMX_NONBONDED_DO_SR
;
271 if (flags
& GMX_FORCE_FORCES
)
273 donb_flags
|= GMX_NONBONDED_DO_FORCE
;
275 if (flags
& GMX_FORCE_ENERGY
)
277 donb_flags
|= GMX_NONBONDED_DO_POTENTIAL
;
279 if (flags
& GMX_FORCE_DO_LR
)
281 donb_flags
|= GMX_NONBONDED_DO_LR
;
284 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
285 do_nonbonded(cr
, fr
, x
, f
, f_longrange
, md
, excl
,
286 &enerd
->grpp
, box_size
, nrnb
,
287 lambda
, dvdl_nb
, -1, -1, donb_flags
);
289 /* If we do foreign lambda and we have soft-core interactions
290 * we have to recalculate the (non-linear) energies contributions.
292 if (fepvals
->n_lambda
> 0 && (flags
& GMX_FORCE_DHDL
) && fepvals
->sc_alpha
!= 0)
294 for (i
= 0; i
< enerd
->n_lambda
; i
++)
296 for (j
= 0; j
< efptNR
; j
++)
298 lam_i
[j
] = (i
== 0 ? lambda
[j
] : fepvals
->all_lambda
[j
][i
-1]);
300 reset_foreign_enerdata(enerd
);
301 do_nonbonded(cr
, fr
, x
, f
, f_longrange
, md
, excl
,
302 &(enerd
->foreign_grpp
), box_size
, nrnb
,
303 lam_i
, dvdl_dum
, -1, -1,
304 (donb_flags
& ~GMX_NONBONDED_DO_FORCE
) | GMX_NONBONDED_DO_FOREIGNLAMBDA
);
305 sum_epot(&ir
->opts
, &(enerd
->foreign_grpp
), enerd
->foreign_term
);
306 enerd
->enerpart_lambda
[i
] += enerd
->foreign_term
[F_EPOT
];
309 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
313 /* If we are doing GB, calculate bonded forces and apply corrections
314 * to the solvation forces */
315 /* MRS: Eventually, many need to include free energy contribution here! */
316 if (ir
->implicit_solvent
)
318 wallcycle_sub_start(wcycle
, ewcsBONDED
);
319 calc_gb_forces(cr
, md
, born
, top
, atype
, x
, f
, fr
, idef
,
320 ir
->gb_algorithm
, ir
->sa_algorithm
, nrnb
, bBornRadii
, &pbc
, graph
, enerd
);
321 wallcycle_sub_stop(wcycle
, ewcsBONDED
);
332 if (fepvals
->sc_alpha
!= 0)
334 enerd
->dvdl_nonlin
[efptVDW
] += dvdl_nb
[efptVDW
];
338 enerd
->dvdl_lin
[efptVDW
] += dvdl_nb
[efptVDW
];
341 if (fepvals
->sc_alpha
!= 0)
343 /* even though coulomb part is linear, we already added it, beacuse we
344 need to go through the vdw calculation anyway */
346 enerd
->dvdl_nonlin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
350 enerd
->dvdl_lin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
356 for (i
= 0; i
< enerd
->grpp
.nener
; i
++)
360 enerd
->grpp
.ener
[egBHAMSR
][i
] :
361 enerd
->grpp
.ener
[egLJSR
][i
])
362 + enerd
->grpp
.ener
[egCOULSR
][i
] + enerd
->grpp
.ener
[egGB
][i
];
364 dvdlsum
= dvdl_nb
[efptVDW
] + dvdl_nb
[efptCOUL
];
365 PRINT_SEPDVDL("VdW and Coulomb SR particle-p.", Vsr
, dvdlsum
);
369 GMX_MPE_LOG(ev_do_fnbf_finish
);
373 pr_rvecs(debug
, 0, "fshift after SR", fr
->fshift
, SHIFTS
);
376 /* Shift the coordinates. Must be done before bonded forces and PPPM,
377 * but is also necessary for SHAKE and update, therefore it can NOT
378 * go when no bonded forces have to be evaluated.
381 /* Here sometimes we would not need to shift with NBFonly,
382 * but we do so anyhow for consistency of the returned coordinates.
386 shift_self(graph
, box
, x
);
389 inc_nrnb(nrnb
, eNR_SHIFTX
, 2*graph
->nnodes
);
393 inc_nrnb(nrnb
, eNR_SHIFTX
, graph
->nnodes
);
396 /* Check whether we need to do bondeds or correct for exclusions */
398 ((flags
& GMX_FORCE_BONDED
)
399 || EEL_RF(fr
->eeltype
) || EEL_FULL(fr
->eeltype
)))
401 /* Since all atoms are in the rectangular or triclinic unit-cell,
402 * only single box vector shifts (2 in x) are required.
404 set_pbc_dd(&pbc
, fr
->ePBC
, cr
->dd
, TRUE
, box
);
408 if (flags
& GMX_FORCE_BONDED
)
410 GMX_MPE_LOG(ev_calc_bonds_start
);
412 wallcycle_sub_start(wcycle
, ewcsBONDED
);
413 calc_bonds(fplog
, cr
->ms
,
414 idef
, x
, hist
, f
, fr
, &pbc
, graph
, enerd
, nrnb
, lambda
, md
, fcd
,
415 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: NULL
, atype
, born
,
417 fr
->bSepDVDL
&& do_per_step(step
, ir
->nstlog
), step
);
419 /* Check if we have to determine energy differences
420 * at foreign lambda's.
422 if (fepvals
->n_lambda
> 0 && (flags
& GMX_FORCE_DHDL
) &&
423 idef
->ilsort
!= ilsortNO_FE
)
425 if (idef
->ilsort
!= ilsortFE_SORTED
)
427 gmx_incons("The bonded interactions are not sorted for free energy");
429 for (i
= 0; i
< enerd
->n_lambda
; i
++)
431 reset_foreign_enerdata(enerd
);
432 for (j
= 0; j
< efptNR
; j
++)
434 lam_i
[j
] = (i
== 0 ? lambda
[j
] : fepvals
->all_lambda
[j
][i
-1]);
436 calc_bonds_lambda(fplog
, idef
, x
, fr
, &pbc
, graph
, &(enerd
->foreign_grpp
), enerd
->foreign_term
, nrnb
, lam_i
, md
,
437 fcd
, DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: NULL
);
438 sum_epot(&ir
->opts
, &(enerd
->foreign_grpp
), enerd
->foreign_term
);
439 enerd
->enerpart_lambda
[i
] += enerd
->foreign_term
[F_EPOT
];
443 GMX_MPE_LOG(ev_calc_bonds_finish
);
444 wallcycle_sub_stop(wcycle
, ewcsBONDED
);
450 if (EEL_FULL(fr
->eeltype
))
452 bSB
= (ir
->nwall
== 2);
456 svmul(ir
->wall_ewald_zfac
, boxs
[ZZ
], boxs
[ZZ
]);
457 box_size
[ZZ
] *= ir
->wall_ewald_zfac
;
460 clear_mat(fr
->vir_el_recip
);
467 /* With the Verlet scheme exclusion forces are calculated
468 * in the non-bonded kernel.
470 /* The TPI molecule does not have exclusions with the rest
471 * of the system and no intra-molecular PME grid contributions
472 * will be calculated in gmx_pme_calc_energy.
474 if ((ir
->cutoff_scheme
== ecutsGROUP
&& fr
->n_tpi
== 0) ||
475 ir
->ewald_geometry
!= eewg3D
||
476 ir
->epsilon_surface
!= 0)
480 wallcycle_sub_start(wcycle
, ewcsEWALD_CORRECTION
);
484 gmx_fatal(FARGS
, "TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
487 nthreads
= gmx_omp_nthreads_get(emntBonded
);
488 #pragma omp parallel for num_threads(nthreads) schedule(static)
489 for (t
= 0; t
< nthreads
; t
++)
494 real
*Vcorrt
, *dvdlt
;
497 fnv
= fr
->f_novirsum
;
498 vir
= &fr
->vir_el_recip
;
505 vir
= &fr
->f_t
[t
].vir
;
506 Vcorrt
= &fr
->f_t
[t
].Vcorr
;
507 dvdlt
= &fr
->f_t
[t
].dvdl
[efptCOUL
];
508 for (i
= 0; i
< fr
->natoms_force
; i
++)
516 ewald_LRcorrection(fplog
,
517 fr
->excl_load
[t
], fr
->excl_load
[t
+1],
520 md
->nChargePerturbed
? md
->chargeB
: NULL
,
521 ir
->cutoff_scheme
!= ecutsVERLET
,
522 excl
, x
, bSB
? boxs
: box
, mu_tot
,
526 lambda
[efptCOUL
], dvdlt
);
530 reduce_thread_forces(fr
->natoms_force
, fr
->f_novirsum
,
532 &Vcorr
, efptCOUL
, &dvdl
,
536 wallcycle_sub_stop(wcycle
, ewcsEWALD_CORRECTION
);
541 Vcorr
+= ewald_charge_correction(cr
, fr
, lambda
[efptCOUL
], box
,
542 &dvdl
, fr
->vir_el_recip
);
545 PRINT_SEPDVDL("Ewald excl./charge/dip. corr.", Vcorr
, dvdl
);
546 enerd
->dvdl_lin
[efptCOUL
] += dvdl
;
557 case eelPMEUSERSWITCH
:
559 if (cr
->duty
& DUTY_PME
)
561 assert(fr
->n_tpi
>= 0);
562 if (fr
->n_tpi
== 0 || (flags
& GMX_FORCE_STATECHANGED
))
564 pme_flags
= GMX_PME_SPREAD_Q
| GMX_PME_SOLVE
;
565 if (flags
& GMX_FORCE_FORCES
)
567 pme_flags
|= GMX_PME_CALC_F
;
569 if (flags
& (GMX_FORCE_VIRIAL
| GMX_FORCE_ENERGY
))
571 pme_flags
|= GMX_PME_CALC_ENER_VIR
;
575 /* We don't calculate f, but we do want the potential */
576 pme_flags
|= GMX_PME_CALC_POT
;
578 wallcycle_start(wcycle
, ewcPMEMESH
);
579 status
= gmx_pme_do(fr
->pmedata
,
580 md
->start
, md
->homenr
- fr
->n_tpi
,
582 md
->chargeA
, md
->chargeB
,
583 bSB
? boxs
: box
, cr
,
584 DOMAINDECOMP(cr
) ? dd_pme_maxshift_x(cr
->dd
) : 0,
585 DOMAINDECOMP(cr
) ? dd_pme_maxshift_y(cr
->dd
) : 0,
587 fr
->vir_el_recip
, fr
->ewaldcoeff
,
588 &Vlr
, lambda
[efptCOUL
], &dvdl
,
590 *cycles_pme
= wallcycle_stop(wcycle
, ewcPMEMESH
);
592 /* We should try to do as little computation after
593 * this as possible, because parallel PME synchronizes
594 * the nodes, so we want all load imbalance of the rest
595 * of the force calculation to be before the PME call.
596 * DD load balancing is done on the whole time of
597 * the force call (without PME).
602 /* Determine the PME grid energy of the test molecule
603 * with the PME grid potential of the other charges.
605 gmx_pme_calc_energy(fr
->pmedata
, fr
->n_tpi
,
606 x
+ md
->homenr
- fr
->n_tpi
,
607 md
->chargeA
+ md
->homenr
- fr
->n_tpi
,
610 PRINT_SEPDVDL("PME mesh", Vlr
, dvdl
);
614 Vlr
= do_ewald(fplog
, FALSE
, ir
, x
, fr
->f_novirsum
,
615 md
->chargeA
, md
->chargeB
,
616 box_size
, cr
, md
->homenr
,
617 fr
->vir_el_recip
, fr
->ewaldcoeff
,
618 lambda
[efptCOUL
], &dvdl
, fr
->ewald_table
);
619 PRINT_SEPDVDL("Ewald long-range", Vlr
, dvdl
);
622 gmx_fatal(FARGS
, "No such electrostatics method implemented %s",
623 eel_names
[fr
->eeltype
]);
627 gmx_fatal(FARGS
, "Error %d in long range electrostatics routine %s",
628 status
, EELTYPE(fr
->eeltype
));
630 /* Note that with separate PME nodes we get the real energies later */
631 enerd
->dvdl_lin
[efptCOUL
] += dvdl
;
632 enerd
->term
[F_COUL_RECIP
] = Vlr
+ Vcorr
;
635 fprintf(debug
, "Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
636 Vlr
, Vcorr
, enerd
->term
[F_COUL_RECIP
]);
637 pr_rvecs(debug
, 0, "vir_el_recip after corr", fr
->vir_el_recip
, DIM
);
638 pr_rvecs(debug
, 0, "fshift after LR Corrections", fr
->fshift
, SHIFTS
);
643 if (EEL_RF(fr
->eeltype
))
645 /* With the Verlet scheme exclusion forces are calculated
646 * in the non-bonded kernel.
648 if (ir
->cutoff_scheme
!= ecutsVERLET
&& fr
->eeltype
!= eelRF_NEC
)
651 enerd
->term
[F_RF_EXCL
] =
652 RF_excl_correction(fplog
, fr
, graph
, md
, excl
, x
, f
,
653 fr
->fshift
, &pbc
, lambda
[efptCOUL
], &dvdl
);
656 enerd
->dvdl_lin
[efptCOUL
] += dvdl
;
657 PRINT_SEPDVDL("RF exclusion correction",
658 enerd
->term
[F_RF_EXCL
], dvdl
);
666 print_nrnb(debug
, nrnb
);
674 MPI_Barrier(cr
->mpi_comm_mygroup
);
677 if (fr
->timesteps
== 11)
679 fprintf(stderr
, "* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
680 cr
->nodeid
, gmx_step_str(fr
->timesteps
, buf
),
681 100*fr
->t_wait
/(fr
->t_wait
+fr
->t_fnbf
),
682 (fr
->t_fnbf
+fr
->t_wait
)/fr
->t_fnbf
);
690 pr_rvecs(debug
, 0, "fshift after bondeds", fr
->fshift
, SHIFTS
);
693 GMX_MPE_LOG(ev_force_finish
);
697 void init_enerdata(int ngener
, int n_lambda
, gmx_enerdata_t
*enerd
)
701 for (i
= 0; i
< F_NRE
; i
++)
704 enerd
->foreign_term
[i
] = 0;
708 for (i
= 0; i
< efptNR
; i
++)
710 enerd
->dvdl_lin
[i
] = 0;
711 enerd
->dvdl_nonlin
[i
] = 0;
717 fprintf(debug
, "Creating %d sized group matrix for energies\n", n2
);
719 enerd
->grpp
.nener
= n2
;
720 enerd
->foreign_grpp
.nener
= n2
;
721 for (i
= 0; (i
< egNR
); i
++)
723 snew(enerd
->grpp
.ener
[i
], n2
);
724 snew(enerd
->foreign_grpp
.ener
[i
], n2
);
729 enerd
->n_lambda
= 1 + n_lambda
;
730 snew(enerd
->enerpart_lambda
, enerd
->n_lambda
);
738 void destroy_enerdata(gmx_enerdata_t
*enerd
)
742 for (i
= 0; (i
< egNR
); i
++)
744 sfree(enerd
->grpp
.ener
[i
]);
747 for (i
= 0; (i
< egNR
); i
++)
749 sfree(enerd
->foreign_grpp
.ener
[i
]);
754 sfree(enerd
->enerpart_lambda
);
758 static real
sum_v(int n
, real v
[])
764 for (i
= 0; (i
< n
); i
++)
772 void sum_epot(t_grpopts
*opts
, gmx_grppairener_t
*grpp
, real
*epot
)
776 /* Accumulate energies */
777 epot
[F_COUL_SR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULSR
]);
778 epot
[F_LJ
] = sum_v(grpp
->nener
, grpp
->ener
[egLJSR
]);
779 epot
[F_LJ14
] = sum_v(grpp
->nener
, grpp
->ener
[egLJ14
]);
780 epot
[F_COUL14
] = sum_v(grpp
->nener
, grpp
->ener
[egCOUL14
]);
781 epot
[F_COUL_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULLR
]);
782 epot
[F_LJ_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egLJLR
]);
783 /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
784 epot
[F_GBPOL
] += sum_v(grpp
->nener
, grpp
->ener
[egGB
]);
786 /* lattice part of LR doesnt belong to any group
787 * and has been added earlier
789 epot
[F_BHAM
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMSR
]);
790 epot
[F_BHAM_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMLR
]);
793 for (i
= 0; (i
< F_EPOT
); i
++)
795 if (i
!= F_DISRESVIOL
&& i
!= F_ORIRESDEV
)
797 epot
[F_EPOT
] += epot
[i
];
802 void sum_dhdl(gmx_enerdata_t
*enerd
, real
*lambda
, t_lambda
*fepvals
)
807 enerd
->dvdl_lin
[efptVDW
] += enerd
->term
[F_DVDL_VDW
]; /* include dispersion correction */
808 enerd
->term
[F_DVDL
] = 0.0;
809 for (i
= 0; i
< efptNR
; i
++)
811 if (fepvals
->separate_dvdl
[i
])
813 /* could this be done more readably/compactly? */
826 index
= F_DVDL_BONDED
;
828 case (efptRESTRAINT
):
829 index
= F_DVDL_RESTRAINT
;
835 enerd
->term
[index
] = enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
838 fprintf(debug
, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
839 efpt_names
[i
], i
, enerd
->term
[index
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
844 enerd
->term
[F_DVDL
] += enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
847 fprintf(debug
, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
848 efpt_names
[0], i
, enerd
->term
[F_DVDL
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
853 /* Notes on the foreign lambda free energy difference evaluation:
854 * Adding the potential and ekin terms that depend linearly on lambda
855 * as delta lam * dvdl to the energy differences is exact.
856 * For the constraints this is not exact, but we have no other option
857 * without literally changing the lengths and reevaluating the energies at each step.
858 * (try to remedy this post 4.6 - MRS)
859 * For the non-bonded LR term we assume that the soft-core (if present)
860 * no longer affects the energy beyond the short-range cut-off,
861 * which is a very good approximation (except for exotic settings).
862 * (investigate how to overcome this post 4.6 - MRS)
865 for (i
= 0; i
< fepvals
->n_lambda
; i
++)
866 { /* note we are iterating over fepvals here!
867 For the current lam, dlam = 0 automatically,
868 so we don't need to add anything to the
869 enerd->enerpart_lambda[0] */
871 /* we don't need to worry about dvdl_lin contributions to dE at
872 current lambda, because the contributions to the current
873 lambda are automatically zeroed */
875 for (j
= 0; j
< efptNR
; j
++)
877 /* Note that this loop is over all dhdl components, not just the separated ones */
878 dlam
= (fepvals
->all_lambda
[j
][i
]-lambda
[j
]);
879 enerd
->enerpart_lambda
[i
+1] += dlam
*enerd
->dvdl_lin
[j
];
882 fprintf(debug
, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
883 fepvals
->all_lambda
[j
][i
], efpt_names
[j
],
884 (enerd
->enerpart_lambda
[i
+1] - enerd
->enerpart_lambda
[0]),
885 dlam
, enerd
->dvdl_lin
[j
]);
892 void reset_foreign_enerdata(gmx_enerdata_t
*enerd
)
896 /* First reset all foreign energy components. Foreign energies always called on
897 neighbor search steps */
898 for (i
= 0; (i
< egNR
); i
++)
900 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
902 enerd
->foreign_grpp
.ener
[i
][j
] = 0.0;
906 /* potential energy components */
907 for (i
= 0; (i
<= F_EPOT
); i
++)
909 enerd
->foreign_term
[i
] = 0.0;
913 void reset_enerdata(t_grpopts
*opts
,
914 t_forcerec
*fr
, gmx_bool bNS
,
915 gmx_enerdata_t
*enerd
,
921 /* First reset all energy components, except for the long range terms
922 * on the master at non neighbor search steps, since the long range
923 * terms have already been summed at the last neighbor search step.
925 bKeepLR
= (fr
->bTwinRange
&& !bNS
);
926 for (i
= 0; (i
< egNR
); i
++)
928 if (!(bKeepLR
&& bMaster
&& (i
== egCOULLR
|| i
== egLJLR
)))
930 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
932 enerd
->grpp
.ener
[i
][j
] = 0.0;
936 for (i
= 0; i
< efptNR
; i
++)
938 enerd
->dvdl_lin
[i
] = 0.0;
939 enerd
->dvdl_nonlin
[i
] = 0.0;
942 /* Normal potential energy components */
943 for (i
= 0; (i
<= F_EPOT
); i
++)
945 enerd
->term
[i
] = 0.0;
947 /* Initialize the dVdlambda term with the long range contribution */
948 /* Initialize the dvdl term with the long range contribution */
949 enerd
->term
[F_DVDL
] = 0.0;
950 enerd
->term
[F_DVDL_COUL
] = 0.0;
951 enerd
->term
[F_DVDL_VDW
] = 0.0;
952 enerd
->term
[F_DVDL_BONDED
] = 0.0;
953 enerd
->term
[F_DVDL_RESTRAINT
] = 0.0;
954 enerd
->term
[F_DKDL
] = 0.0;
955 if (enerd
->n_lambda
> 0)
957 for (i
= 0; i
< enerd
->n_lambda
; i
++)
959 enerd
->enerpart_lambda
[i
] = 0.0;
962 /* reset foreign energy data - separate function since we also call it elsewhere */
963 reset_foreign_enerdata(enerd
);