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, 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.
48 #include "nonbonded.h"
62 #include "gmx_omp_nthreads.h"
64 #include "gromacs/legacyheaders/types/commrec.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/timing/wallcycle.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/smalloc.h"
79 gmx_bool bDoLongRangeNS
)
85 if (!fr
->ns
.nblist_initialized
)
87 init_neighbor_list(fp
, fr
, md
->homenr
);
95 nsearch
= search_neighbours(fp
, fr
, box
, top
, groups
, cr
, nrnb
, md
,
96 bFillGrid
, bDoLongRangeNS
);
99 fprintf(debug
, "nsearch = %d\n", nsearch
);
102 /* Check whether we have to do dynamic load balancing */
103 /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
104 count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
105 &(top->idef),opts->ngener);
107 if (fr
->ns
.dump_nl
> 0)
109 dump_nblist(fp
, cr
, fr
, fr
->ns
.dump_nl
);
113 static void reduce_thread_forces(int n
, rvec
*f
,
114 tensor vir_q
, tensor vir_lj
,
115 real
*Vcorr_q
, real
*Vcorr_lj
,
116 real
*dvdl_q
, real
*dvdl_lj
,
117 int nthreads
, f_thread_t
*f_t
)
121 /* This reduction can run over any number of threads */
122 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
123 for (i
= 0; i
< n
; i
++)
125 for (t
= 1; t
< nthreads
; t
++)
127 rvec_inc(f
[i
], f_t
[t
].f
[i
]);
130 for (t
= 1; t
< nthreads
; t
++)
132 *Vcorr_q
+= f_t
[t
].Vcorr_q
;
133 *Vcorr_lj
+= f_t
[t
].Vcorr_lj
;
134 *dvdl_q
+= f_t
[t
].dvdl
[efptCOUL
];
135 *dvdl_lj
+= f_t
[t
].dvdl
[efptVDW
];
136 m_add(vir_q
, f_t
[t
].vir_q
, vir_q
);
137 m_add(vir_lj
, f_t
[t
].vir_lj
, vir_lj
);
141 void gmx_print_sepdvdl(FILE *fplog
, const char *s
, real v
, real dvdlambda
)
143 fprintf(fplog
, " %-30s V %12.5e dVdl %12.5e\n", s
, v
, dvdlambda
);
146 void do_force_lowlevel(FILE *fplog
, gmx_int64_t step
,
147 t_forcerec
*fr
, t_inputrec
*ir
,
148 t_idef
*idef
, t_commrec
*cr
,
149 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
151 rvec x
[], history_t
*hist
,
154 gmx_enerdata_t
*enerd
,
171 gmx_bool bDoEpot
, bSepDVDL
, bSB
;
177 double clam_i
, vlam_i
;
178 real dvdl_dum
[efptNR
], dvdl_nb
[efptNR
], lam_i
[efptNR
];
179 real dvdl_q
, dvdl_lj
;
182 double t0
= 0.0, t1
, t2
, t3
; /* time measurement for coarse load balancing */
185 #define PRINT_SEPDVDL(s, v, dvdlambda) if (bSepDVDL) { gmx_print_sepdvdl(fplog, s, v, dvdlambda); }
187 set_pbc(&pbc
, fr
->ePBC
, box
);
189 /* reset free energy components */
190 for (i
= 0; i
< efptNR
; i
++)
197 for (i
= 0; (i
< DIM
); i
++)
199 box_size
[i
] = box
[i
][i
];
202 bSepDVDL
= (fr
->bSepDVDL
&& do_per_step(step
, ir
->nstlog
));
205 /* do QMMM first if requested */
208 enerd
->term
[F_EQM
] = calculate_QMMM(cr
, x
, f
, fr
);
213 fprintf(fplog
, "Step %s: non-bonded V and dVdl for node %d:\n",
214 gmx_step_str(step
, buf
), cr
->nodeid
);
217 /* Call the short range functions all in one go. */
220 /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
221 #define TAKETIME FALSE
224 MPI_Barrier(cr
->mpi_comm_mygroup
);
231 /* foreign lambda component for walls */
232 real dvdl_walls
= do_walls(ir
, fr
, box
, md
, x
, f
, lambda
[efptVDW
],
233 enerd
->grpp
.ener
[egLJSR
], nrnb
);
234 PRINT_SEPDVDL("Walls", 0.0, dvdl_walls
);
235 enerd
->dvdl_lin
[efptVDW
] += dvdl_walls
;
238 /* If doing GB, reset dvda and calculate the Born radii */
239 if (ir
->implicit_solvent
)
241 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
243 for (i
= 0; i
< born
->nr
; i
++)
250 calc_gb_rad(cr
, fr
, ir
, top
, x
, &(fr
->gblist
), born
, md
, nrnb
);
253 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
257 /* We only do non-bonded calculation with group scheme here, the verlet
258 * calls are done from do_force_cutsVERLET(). */
259 if (fr
->cutoff_scheme
== ecutsGROUP
&& (flags
& GMX_FORCE_NONBONDED
))
262 /* Add short-range interactions */
263 donb_flags
|= GMX_NONBONDED_DO_SR
;
265 /* Currently all group scheme kernels always calculate (shift-)forces */
266 if (flags
& GMX_FORCE_FORCES
)
268 donb_flags
|= GMX_NONBONDED_DO_FORCE
;
270 if (flags
& GMX_FORCE_VIRIAL
)
272 donb_flags
|= GMX_NONBONDED_DO_SHIFTFORCE
;
274 if (flags
& GMX_FORCE_ENERGY
)
276 donb_flags
|= GMX_NONBONDED_DO_POTENTIAL
;
278 if (flags
& GMX_FORCE_DO_LR
)
280 donb_flags
|= GMX_NONBONDED_DO_LR
;
283 wallcycle_sub_start(wcycle
, ewcsNONBONDED
);
284 do_nonbonded(fr
, x
, f
, f_longrange
, md
, excl
,
286 lambda
, dvdl_nb
, -1, -1, donb_flags
);
288 /* If we do foreign lambda and we have soft-core interactions
289 * we have to recalculate the (non-linear) energies contributions.
291 if (fepvals
->n_lambda
> 0 && (flags
& GMX_FORCE_DHDL
) && fepvals
->sc_alpha
!= 0)
293 for (i
= 0; i
< enerd
->n_lambda
; i
++)
295 for (j
= 0; j
< efptNR
; j
++)
297 lam_i
[j
] = (i
== 0 ? lambda
[j
] : fepvals
->all_lambda
[j
][i
-1]);
299 reset_foreign_enerdata(enerd
);
300 do_nonbonded(fr
, x
, f
, f_longrange
, md
, excl
,
301 &(enerd
->foreign_grpp
), nrnb
,
302 lam_i
, dvdl_dum
, -1, -1,
303 (donb_flags
& ~GMX_NONBONDED_DO_FORCE
) | GMX_NONBONDED_DO_FOREIGNLAMBDA
);
304 sum_epot(&(enerd
->foreign_grpp
), enerd
->foreign_term
);
305 enerd
->enerpart_lambda
[i
] += enerd
->foreign_term
[F_EPOT
];
308 wallcycle_sub_stop(wcycle
, ewcsNONBONDED
);
312 /* If we are doing GB, calculate bonded forces and apply corrections
313 * to the solvation forces */
314 /* MRS: Eventually, many need to include free energy contribution here! */
315 if (ir
->implicit_solvent
)
317 wallcycle_sub_start(wcycle
, ewcsBONDED
);
318 calc_gb_forces(cr
, md
, born
, top
, x
, f
, fr
, idef
,
319 ir
->gb_algorithm
, ir
->sa_algorithm
, nrnb
, &pbc
, graph
, enerd
);
320 wallcycle_sub_stop(wcycle
, ewcsBONDED
);
331 if (fepvals
->sc_alpha
!= 0)
333 enerd
->dvdl_nonlin
[efptVDW
] += dvdl_nb
[efptVDW
];
337 enerd
->dvdl_lin
[efptVDW
] += dvdl_nb
[efptVDW
];
340 if (fepvals
->sc_alpha
!= 0)
342 /* even though coulomb part is linear, we already added it, beacuse we
343 need to go through the vdw calculation anyway */
345 enerd
->dvdl_nonlin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
349 enerd
->dvdl_lin
[efptCOUL
] += dvdl_nb
[efptCOUL
];
354 real V_short_range
= 0;
355 real dvdl_short_range
= 0;
357 for (i
= 0; i
< enerd
->grpp
.nener
; i
++)
361 enerd
->grpp
.ener
[egBHAMSR
][i
] :
362 enerd
->grpp
.ener
[egLJSR
][i
])
363 + enerd
->grpp
.ener
[egCOULSR
][i
] + enerd
->grpp
.ener
[egGB
][i
];
365 dvdl_short_range
= dvdl_nb
[efptVDW
] + dvdl_nb
[efptCOUL
];
366 PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",
375 pr_rvecs(debug
, 0, "fshift after SR", fr
->fshift
, SHIFTS
);
378 /* Shift the coordinates. Must be done before bonded forces and PPPM,
379 * but is also necessary for SHAKE and update, therefore it can NOT
380 * go when no bonded forces have to be evaluated.
383 /* Here sometimes we would not need to shift with NBFonly,
384 * but we do so anyhow for consistency of the returned coordinates.
388 shift_self(graph
, box
, x
);
391 inc_nrnb(nrnb
, eNR_SHIFTX
, 2*graph
->nnodes
);
395 inc_nrnb(nrnb
, eNR_SHIFTX
, graph
->nnodes
);
398 /* Check whether we need to do bondeds or correct for exclusions */
400 ((flags
& GMX_FORCE_BONDED
)
401 || EEL_RF(fr
->eeltype
) || EEL_FULL(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
)))
403 /* Since all atoms are in the rectangular or triclinic unit-cell,
404 * only single box vector shifts (2 in x) are required.
406 set_pbc_dd(&pbc
, fr
->ePBC
, cr
->dd
, TRUE
, box
);
410 if (flags
& GMX_FORCE_BONDED
)
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(&(enerd
->foreign_grpp
), enerd
->foreign_term
);
439 enerd
->enerpart_lambda
[i
] += enerd
->foreign_term
[F_EPOT
];
444 wallcycle_sub_stop(wcycle
, ewcsBONDED
);
450 clear_mat(fr
->vir_el_recip
);
451 clear_mat(fr
->vir_lj_recip
);
453 /* Do long-range electrostatics and/or LJ-PME, including related short-range
456 if (EEL_FULL(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
))
458 real Vlr
= 0, Vcorr
= 0;
459 real dvdl_long_range
= 0;
461 real Vlr_q
= 0, Vlr_lj
= 0, Vcorr_q
= 0, Vcorr_lj
= 0;
462 real dvdl_long_range_q
= 0, dvdl_long_range_lj
= 0;
464 bSB
= (ir
->nwall
== 2);
468 svmul(ir
->wall_ewald_zfac
, boxs
[ZZ
], boxs
[ZZ
]);
469 box_size
[ZZ
] *= ir
->wall_ewald_zfac
;
472 if (EEL_PME_EWALD(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
))
474 real dvdl_long_range_correction_q
= 0;
475 real dvdl_long_range_correction_lj
= 0;
476 /* With the Verlet scheme exclusion forces are calculated
477 * in the non-bonded kernel.
479 /* The TPI molecule does not have exclusions with the rest
480 * of the system and no intra-molecular PME grid
481 * contributions will be calculated in
482 * gmx_pme_calc_energy.
484 if ((ir
->cutoff_scheme
== ecutsGROUP
&& fr
->n_tpi
== 0) ||
485 ir
->ewald_geometry
!= eewg3D
||
486 ir
->epsilon_surface
!= 0)
490 wallcycle_sub_start(wcycle
, ewcsEWALD_CORRECTION
);
494 gmx_fatal(FARGS
, "TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
497 nthreads
= gmx_omp_nthreads_get(emntBonded
);
498 #pragma omp parallel for num_threads(nthreads) schedule(static)
499 for (t
= 0; t
< nthreads
; t
++)
503 tensor
*vir_q
, *vir_lj
;
504 real
*Vcorrt_q
, *Vcorrt_lj
, *dvdlt_q
, *dvdlt_lj
;
507 fnv
= fr
->f_novirsum
;
508 vir_q
= &fr
->vir_el_recip
;
509 vir_lj
= &fr
->vir_lj_recip
;
511 Vcorrt_lj
= &Vcorr_lj
;
512 dvdlt_q
= &dvdl_long_range_correction_q
;
513 dvdlt_lj
= &dvdl_long_range_correction_lj
;
518 vir_q
= &fr
->f_t
[t
].vir_q
;
519 vir_lj
= &fr
->f_t
[t
].vir_lj
;
520 Vcorrt_q
= &fr
->f_t
[t
].Vcorr_q
;
521 Vcorrt_lj
= &fr
->f_t
[t
].Vcorr_lj
;
522 dvdlt_q
= &fr
->f_t
[t
].dvdl
[efptCOUL
];
523 dvdlt_lj
= &fr
->f_t
[t
].dvdl
[efptVDW
];
524 for (i
= 0; i
< fr
->natoms_force
; i
++)
534 ewald_LRcorrection(fr
->excl_load
[t
], fr
->excl_load
[t
+1],
537 md
->nChargePerturbed
? md
->chargeB
: NULL
,
539 md
->nTypePerturbed
? md
->sqrt_c6B
: NULL
,
541 md
->nTypePerturbed
? md
->sigmaB
: NULL
,
543 md
->nTypePerturbed
? md
->sigma3B
: NULL
,
544 ir
->cutoff_scheme
!= ecutsVERLET
,
545 excl
, x
, bSB
? boxs
: box
, mu_tot
,
548 fnv
, *vir_q
, *vir_lj
,
550 lambda
[efptCOUL
], lambda
[efptVDW
],
555 reduce_thread_forces(fr
->natoms_force
, fr
->f_novirsum
,
556 fr
->vir_el_recip
, fr
->vir_lj_recip
,
558 &dvdl_long_range_correction_q
,
559 &dvdl_long_range_correction_lj
,
562 wallcycle_sub_stop(wcycle
, ewcsEWALD_CORRECTION
);
565 if (EEL_PME_EWALD(fr
->eeltype
) && fr
->n_tpi
== 0)
567 Vcorr_q
+= ewald_charge_correction(cr
, fr
, lambda
[efptCOUL
], box
,
568 &dvdl_long_range_correction_q
,
572 PRINT_SEPDVDL("Ewald excl./charge/dip. corr.", Vcorr_q
, dvdl_long_range_correction_q
);
573 PRINT_SEPDVDL("Ewald excl. corr. LJ", Vcorr_lj
, dvdl_long_range_correction_lj
);
574 enerd
->dvdl_lin
[efptCOUL
] += dvdl_long_range_correction_q
;
575 enerd
->dvdl_lin
[efptVDW
] += dvdl_long_range_correction_lj
;
577 if ((EEL_PME(fr
->eeltype
) || EVDW_PME(fr
->vdwtype
)) && (cr
->duty
& DUTY_PME
))
579 /* Do reciprocal PME for Coulomb and/or LJ. */
580 assert(fr
->n_tpi
>= 0);
581 if (fr
->n_tpi
== 0 || (flags
& GMX_FORCE_STATECHANGED
))
583 pme_flags
= GMX_PME_SPREAD
| GMX_PME_SOLVE
;
584 if (EEL_PME(fr
->eeltype
))
586 pme_flags
|= GMX_PME_DO_COULOMB
;
588 if (EVDW_PME(fr
->vdwtype
))
590 pme_flags
|= GMX_PME_DO_LJ
;
592 if (flags
& GMX_FORCE_FORCES
)
594 pme_flags
|= GMX_PME_CALC_F
;
596 if (flags
& GMX_FORCE_VIRIAL
)
598 pme_flags
|= GMX_PME_CALC_ENER_VIR
;
602 /* We don't calculate f, but we do want the potential */
603 pme_flags
|= GMX_PME_CALC_POT
;
605 wallcycle_start(wcycle
, ewcPMEMESH
);
606 status
= gmx_pme_do(fr
->pmedata
,
607 0, md
->homenr
- fr
->n_tpi
,
609 md
->chargeA
, md
->chargeB
,
610 md
->sqrt_c6A
, md
->sqrt_c6B
,
611 md
->sigmaA
, md
->sigmaB
,
612 bSB
? boxs
: box
, cr
,
613 DOMAINDECOMP(cr
) ? dd_pme_maxshift_x(cr
->dd
) : 0,
614 DOMAINDECOMP(cr
) ? dd_pme_maxshift_y(cr
->dd
) : 0,
616 fr
->vir_el_recip
, fr
->ewaldcoeff_q
,
617 fr
->vir_lj_recip
, fr
->ewaldcoeff_lj
,
619 lambda
[efptCOUL
], lambda
[efptVDW
],
620 &dvdl_long_range_q
, &dvdl_long_range_lj
, pme_flags
);
621 *cycles_pme
= wallcycle_stop(wcycle
, ewcPMEMESH
);
624 gmx_fatal(FARGS
, "Error %d in reciprocal PME routine", status
);
626 /* We should try to do as little computation after
627 * this as possible, because parallel PME synchronizes
628 * the nodes, so we want all load imbalance of the
629 * rest of the force calculation to be before the PME
630 * call. DD load balancing is done on the whole time
631 * of the force call (without PME).
636 if (EVDW_PME(ir
->vdwtype
))
639 gmx_fatal(FARGS
, "Test particle insertion not implemented with LJ-PME");
641 /* Determine the PME grid energy of the test molecule
642 * with the PME grid potential of the other charges.
644 gmx_pme_calc_energy(fr
->pmedata
, fr
->n_tpi
,
645 x
+ md
->homenr
- fr
->n_tpi
,
646 md
->chargeA
+ md
->homenr
- fr
->n_tpi
,
649 PRINT_SEPDVDL("PME mesh", Vlr_q
+ Vlr_lj
, dvdl_long_range_q
+dvdl_long_range_lj
);
653 if (!EEL_PME(fr
->eeltype
) && EEL_PME_EWALD(fr
->eeltype
))
655 Vlr_q
= do_ewald(ir
, x
, fr
->f_novirsum
,
656 md
->chargeA
, md
->chargeB
,
657 box_size
, cr
, md
->homenr
,
658 fr
->vir_el_recip
, fr
->ewaldcoeff_q
,
659 lambda
[efptCOUL
], &dvdl_long_range_q
, fr
->ewald_table
);
660 PRINT_SEPDVDL("Ewald long-range", Vlr_q
, dvdl_long_range_q
);
663 /* Note that with separate PME nodes we get the real energies later */
664 enerd
->dvdl_lin
[efptCOUL
] += dvdl_long_range_q
;
665 enerd
->dvdl_lin
[efptVDW
] += dvdl_long_range_lj
;
666 enerd
->term
[F_COUL_RECIP
] = Vlr_q
+ Vcorr_q
;
667 enerd
->term
[F_LJ_RECIP
] = Vlr_lj
+ Vcorr_lj
;
670 fprintf(debug
, "Vlr_q = %g, Vcorr_q = %g, Vlr_corr_q = %g\n",
671 Vlr_q
, Vcorr_q
, enerd
->term
[F_COUL_RECIP
]);
672 pr_rvecs(debug
, 0, "vir_el_recip after corr", fr
->vir_el_recip
, DIM
);
673 pr_rvecs(debug
, 0, "fshift after LR Corrections", fr
->fshift
, SHIFTS
);
674 fprintf(debug
, "Vlr_lj: %g, Vcorr_lj = %g, Vlr_corr_lj = %g\n",
675 Vlr_lj
, Vcorr_lj
, enerd
->term
[F_LJ_RECIP
]);
676 pr_rvecs(debug
, 0, "vir_lj_recip after corr", fr
->vir_lj_recip
, DIM
);
681 /* Is there a reaction-field exclusion correction needed? */
682 if (EEL_RF(fr
->eeltype
) && eelRF_NEC
!= fr
->eeltype
)
684 /* With the Verlet scheme, exclusion forces are calculated
685 * in the non-bonded kernel.
687 if (ir
->cutoff_scheme
!= ecutsVERLET
)
689 real dvdl_rf_excl
= 0;
690 enerd
->term
[F_RF_EXCL
] =
691 RF_excl_correction(fr
, graph
, md
, excl
, x
, f
,
692 fr
->fshift
, &pbc
, lambda
[efptCOUL
], &dvdl_rf_excl
);
694 enerd
->dvdl_lin
[efptCOUL
] += dvdl_rf_excl
;
695 PRINT_SEPDVDL("RF exclusion correction",
696 enerd
->term
[F_RF_EXCL
], dvdl_rf_excl
);
705 print_nrnb(debug
, nrnb
);
713 MPI_Barrier(cr
->mpi_comm_mygroup
);
716 if (fr
->timesteps
== 11)
718 fprintf(stderr
, "* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
719 cr
->nodeid
, gmx_step_str(fr
->timesteps
, buf
),
720 100*fr
->t_wait
/(fr
->t_wait
+fr
->t_fnbf
),
721 (fr
->t_fnbf
+fr
->t_wait
)/fr
->t_fnbf
);
729 pr_rvecs(debug
, 0, "fshift after bondeds", fr
->fshift
, SHIFTS
);
734 void init_enerdata(int ngener
, int n_lambda
, gmx_enerdata_t
*enerd
)
738 for (i
= 0; i
< F_NRE
; i
++)
741 enerd
->foreign_term
[i
] = 0;
745 for (i
= 0; i
< efptNR
; i
++)
747 enerd
->dvdl_lin
[i
] = 0;
748 enerd
->dvdl_nonlin
[i
] = 0;
754 fprintf(debug
, "Creating %d sized group matrix for energies\n", n2
);
756 enerd
->grpp
.nener
= n2
;
757 enerd
->foreign_grpp
.nener
= n2
;
758 for (i
= 0; (i
< egNR
); i
++)
760 snew(enerd
->grpp
.ener
[i
], n2
);
761 snew(enerd
->foreign_grpp
.ener
[i
], n2
);
766 enerd
->n_lambda
= 1 + n_lambda
;
767 snew(enerd
->enerpart_lambda
, enerd
->n_lambda
);
775 void destroy_enerdata(gmx_enerdata_t
*enerd
)
779 for (i
= 0; (i
< egNR
); i
++)
781 sfree(enerd
->grpp
.ener
[i
]);
784 for (i
= 0; (i
< egNR
); i
++)
786 sfree(enerd
->foreign_grpp
.ener
[i
]);
791 sfree(enerd
->enerpart_lambda
);
795 static real
sum_v(int n
, real v
[])
801 for (i
= 0; (i
< n
); i
++)
809 void sum_epot(gmx_grppairener_t
*grpp
, real
*epot
)
813 /* Accumulate energies */
814 epot
[F_COUL_SR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULSR
]);
815 epot
[F_LJ
] = sum_v(grpp
->nener
, grpp
->ener
[egLJSR
]);
816 epot
[F_LJ14
] = sum_v(grpp
->nener
, grpp
->ener
[egLJ14
]);
817 epot
[F_COUL14
] = sum_v(grpp
->nener
, grpp
->ener
[egCOUL14
]);
818 epot
[F_COUL_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egCOULLR
]);
819 epot
[F_LJ_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egLJLR
]);
820 /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
821 epot
[F_GBPOL
] += sum_v(grpp
->nener
, grpp
->ener
[egGB
]);
823 /* lattice part of LR doesnt belong to any group
824 * and has been added earlier
826 epot
[F_BHAM
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMSR
]);
827 epot
[F_BHAM_LR
] = sum_v(grpp
->nener
, grpp
->ener
[egBHAMLR
]);
830 for (i
= 0; (i
< F_EPOT
); i
++)
832 if (i
!= F_DISRESVIOL
&& i
!= F_ORIRESDEV
)
834 epot
[F_EPOT
] += epot
[i
];
839 void sum_dhdl(gmx_enerdata_t
*enerd
, real
*lambda
, t_lambda
*fepvals
)
844 enerd
->dvdl_lin
[efptVDW
] += enerd
->term
[F_DVDL_VDW
]; /* include dispersion correction */
845 enerd
->term
[F_DVDL
] = 0.0;
846 for (i
= 0; i
< efptNR
; i
++)
848 if (fepvals
->separate_dvdl
[i
])
850 /* could this be done more readably/compactly? */
863 index
= F_DVDL_BONDED
;
865 case (efptRESTRAINT
):
866 index
= F_DVDL_RESTRAINT
;
872 enerd
->term
[index
] = enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
875 fprintf(debug
, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
876 efpt_names
[i
], i
, enerd
->term
[index
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
881 enerd
->term
[F_DVDL
] += enerd
->dvdl_lin
[i
] + enerd
->dvdl_nonlin
[i
];
884 fprintf(debug
, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
885 efpt_names
[0], i
, enerd
->term
[F_DVDL
], enerd
->dvdl_nonlin
[i
], enerd
->dvdl_lin
[i
]);
890 /* Notes on the foreign lambda free energy difference evaluation:
891 * Adding the potential and ekin terms that depend linearly on lambda
892 * as delta lam * dvdl to the energy differences is exact.
893 * For the constraints this is not exact, but we have no other option
894 * without literally changing the lengths and reevaluating the energies at each step.
895 * (try to remedy this post 4.6 - MRS)
896 * For the non-bonded LR term we assume that the soft-core (if present)
897 * no longer affects the energy beyond the short-range cut-off,
898 * which is a very good approximation (except for exotic settings).
899 * (investigate how to overcome this post 4.6 - MRS)
901 if (fepvals
->separate_dvdl
[efptBONDED
])
903 enerd
->term
[F_DVDL_BONDED
] += enerd
->term
[F_DVDL_CONSTR
];
907 enerd
->term
[F_DVDL
] += enerd
->term
[F_DVDL_CONSTR
];
909 enerd
->term
[F_DVDL_CONSTR
] = 0;
911 for (i
= 0; i
< fepvals
->n_lambda
; i
++)
913 /* note we are iterating over fepvals here!
914 For the current lam, dlam = 0 automatically,
915 so we don't need to add anything to the
916 enerd->enerpart_lambda[0] */
918 /* we don't need to worry about dvdl_lin contributions to dE at
919 current lambda, because the contributions to the current
920 lambda are automatically zeroed */
922 for (j
= 0; j
< efptNR
; j
++)
924 /* Note that this loop is over all dhdl components, not just the separated ones */
925 dlam
= (fepvals
->all_lambda
[j
][i
]-lambda
[j
]);
926 enerd
->enerpart_lambda
[i
+1] += dlam
*enerd
->dvdl_lin
[j
];
929 fprintf(debug
, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
930 fepvals
->all_lambda
[j
][i
], efpt_names
[j
],
931 (enerd
->enerpart_lambda
[i
+1] - enerd
->enerpart_lambda
[0]),
932 dlam
, enerd
->dvdl_lin
[j
]);
939 void reset_foreign_enerdata(gmx_enerdata_t
*enerd
)
943 /* First reset all foreign energy components. Foreign energies always called on
944 neighbor search steps */
945 for (i
= 0; (i
< egNR
); i
++)
947 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
949 enerd
->foreign_grpp
.ener
[i
][j
] = 0.0;
953 /* potential energy components */
954 for (i
= 0; (i
<= F_EPOT
); i
++)
956 enerd
->foreign_term
[i
] = 0.0;
960 void reset_enerdata(t_forcerec
*fr
, gmx_bool bNS
,
961 gmx_enerdata_t
*enerd
,
967 /* First reset all energy components, except for the long range terms
968 * on the master at non neighbor search steps, since the long range
969 * terms have already been summed at the last neighbor search step.
971 bKeepLR
= (fr
->bTwinRange
&& !bNS
);
972 for (i
= 0; (i
< egNR
); i
++)
974 if (!(bKeepLR
&& bMaster
&& (i
== egCOULLR
|| i
== egLJLR
)))
976 for (j
= 0; (j
< enerd
->grpp
.nener
); j
++)
978 enerd
->grpp
.ener
[i
][j
] = 0.0;
982 for (i
= 0; i
< efptNR
; i
++)
984 enerd
->dvdl_lin
[i
] = 0.0;
985 enerd
->dvdl_nonlin
[i
] = 0.0;
988 /* Normal potential energy components */
989 for (i
= 0; (i
<= F_EPOT
); i
++)
991 enerd
->term
[i
] = 0.0;
993 /* Initialize the dVdlambda term with the long range contribution */
994 /* Initialize the dvdl term with the long range contribution */
995 enerd
->term
[F_DVDL
] = 0.0;
996 enerd
->term
[F_DVDL_COUL
] = 0.0;
997 enerd
->term
[F_DVDL_VDW
] = 0.0;
998 enerd
->term
[F_DVDL_BONDED
] = 0.0;
999 enerd
->term
[F_DVDL_RESTRAINT
] = 0.0;
1000 enerd
->term
[F_DKDL
] = 0.0;
1001 if (enerd
->n_lambda
> 0)
1003 for (i
= 0; i
< enerd
->n_lambda
; i
++)
1005 enerd
->enerpart_lambda
[i
] = 0.0;
1008 /* reset foreign energy data - separate function since we also call it elsewhere */
1009 reset_foreign_enerdata(enerd
);