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.
51 #include "gromacs/commandline/filenm.h"
52 #include "gromacs/domdec/domdec.h"
53 #include "gromacs/domdec/domdec_struct.h"
54 #include "gromacs/ewald/ewald.h"
55 #include "gromacs/ewald/ewald-utils.h"
56 #include "gromacs/fileio/filetypes.h"
57 #include "gromacs/gmxlib/network.h"
58 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
59 #include "gromacs/gpu_utils/gpu_utils.h"
60 #include "gromacs/hardware/detecthardware.h"
61 #include "gromacs/listed-forces/manage-threading.h"
62 #include "gromacs/listed-forces/pairs.h"
63 #include "gromacs/math/functions.h"
64 #include "gromacs/math/units.h"
65 #include "gromacs/math/utilities.h"
66 #include "gromacs/math/vec.h"
67 #include "gromacs/mdlib/force.h"
68 #include "gromacs/mdlib/forcerec-threading.h"
69 #include "gromacs/mdlib/gmx_omp_nthreads.h"
70 #include "gromacs/mdlib/md_support.h"
71 #include "gromacs/mdlib/nb_verlet.h"
72 #include "gromacs/mdlib/nbnxn_atomdata.h"
73 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
74 #include "gromacs/mdlib/nbnxn_search.h"
75 #include "gromacs/mdlib/nbnxn_simd.h"
76 #include "gromacs/mdlib/nbnxn_tuning.h"
77 #include "gromacs/mdlib/nbnxn_util.h"
78 #include "gromacs/mdlib/ns.h"
79 #include "gromacs/mdlib/qmmm.h"
80 #include "gromacs/mdlib/sim_util.h"
81 #include "gromacs/mdtypes/commrec.h"
82 #include "gromacs/mdtypes/fcdata.h"
83 #include "gromacs/mdtypes/group.h"
84 #include "gromacs/mdtypes/iforceprovider.h"
85 #include "gromacs/mdtypes/inputrec.h"
86 #include "gromacs/mdtypes/md_enums.h"
87 #include "gromacs/pbcutil/ishift.h"
88 #include "gromacs/pbcutil/pbc.h"
89 #include "gromacs/simd/simd.h"
90 #include "gromacs/tables/forcetable.h"
91 #include "gromacs/topology/mtop_util.h"
92 #include "gromacs/trajectory/trajectoryframe.h"
93 #include "gromacs/utility/cstringutil.h"
94 #include "gromacs/utility/exceptions.h"
95 #include "gromacs/utility/fatalerror.h"
96 #include "gromacs/utility/gmxassert.h"
97 #include "gromacs/utility/logger.h"
98 #include "gromacs/utility/pleasecite.h"
99 #include "gromacs/utility/smalloc.h"
100 #include "gromacs/utility/strconvert.h"
102 #include "nbnxn_gpu_jit_support.h"
104 const char *egrp_nm
[egNR
+1] = {
105 "Coul-SR", "LJ-SR", "Buck-SR",
106 "Coul-14", "LJ-14", nullptr
109 t_forcerec
*mk_forcerec(void)
119 static void pr_nbfp(FILE *fp
, real
*nbfp
, gmx_bool bBHAM
, int atnr
)
123 for (i
= 0; (i
< atnr
); i
++)
125 for (j
= 0; (j
< atnr
); j
++)
127 fprintf(fp
, "%2d - %2d", i
, j
);
130 fprintf(fp
, " a=%10g, b=%10g, c=%10g\n", BHAMA(nbfp
, atnr
, i
, j
),
131 BHAMB(nbfp
, atnr
, i
, j
), BHAMC(nbfp
, atnr
, i
, j
)/6.0);
135 fprintf(fp
, " c6=%10g, c12=%10g\n", C6(nbfp
, atnr
, i
, j
)/6.0,
136 C12(nbfp
, atnr
, i
, j
)/12.0);
143 static real
*mk_nbfp(const gmx_ffparams_t
*idef
, gmx_bool bBHAM
)
151 snew(nbfp
, 3*atnr
*atnr
);
152 for (i
= k
= 0; (i
< atnr
); i
++)
154 for (j
= 0; (j
< atnr
); j
++, k
++)
156 BHAMA(nbfp
, atnr
, i
, j
) = idef
->iparams
[k
].bham
.a
;
157 BHAMB(nbfp
, atnr
, i
, j
) = idef
->iparams
[k
].bham
.b
;
158 /* nbfp now includes the 6.0 derivative prefactor */
159 BHAMC(nbfp
, atnr
, i
, j
) = idef
->iparams
[k
].bham
.c
*6.0;
165 snew(nbfp
, 2*atnr
*atnr
);
166 for (i
= k
= 0; (i
< atnr
); i
++)
168 for (j
= 0; (j
< atnr
); j
++, k
++)
170 /* nbfp now includes the 6.0/12.0 derivative prefactors */
171 C6(nbfp
, atnr
, i
, j
) = idef
->iparams
[k
].lj
.c6
*6.0;
172 C12(nbfp
, atnr
, i
, j
) = idef
->iparams
[k
].lj
.c12
*12.0;
180 static real
*make_ljpme_c6grid(const gmx_ffparams_t
*idef
, t_forcerec
*fr
)
183 real c6
, c6i
, c6j
, c12i
, c12j
, epsi
, epsj
, sigmai
, sigmaj
;
186 /* For LJ-PME simulations, we correct the energies with the reciprocal space
187 * inside of the cut-off. To do this the non-bonded kernels needs to have
188 * access to the C6-values used on the reciprocal grid in pme.c
192 snew(grid
, 2*atnr
*atnr
);
193 for (i
= k
= 0; (i
< atnr
); i
++)
195 for (j
= 0; (j
< atnr
); j
++, k
++)
197 c6i
= idef
->iparams
[i
*(atnr
+1)].lj
.c6
;
198 c12i
= idef
->iparams
[i
*(atnr
+1)].lj
.c12
;
199 c6j
= idef
->iparams
[j
*(atnr
+1)].lj
.c6
;
200 c12j
= idef
->iparams
[j
*(atnr
+1)].lj
.c12
;
201 c6
= std::sqrt(c6i
* c6j
);
202 if (fr
->ljpme_combination_rule
== eljpmeLB
203 && !gmx_numzero(c6
) && !gmx_numzero(c12i
) && !gmx_numzero(c12j
))
205 sigmai
= gmx::sixthroot(c12i
/ c6i
);
206 sigmaj
= gmx::sixthroot(c12j
/ c6j
);
207 epsi
= c6i
* c6i
/ c12i
;
208 epsj
= c6j
* c6j
/ c12j
;
209 c6
= std::sqrt(epsi
* epsj
) * gmx::power6(0.5*(sigmai
+sigmaj
));
211 /* Store the elements at the same relative positions as C6 in nbfp in order
212 * to simplify access in the kernels
214 grid
[2*(atnr
*i
+j
)] = c6
*6.0;
220 static real
*mk_nbfp_combination_rule(const gmx_ffparams_t
*idef
, int comb_rule
)
224 real c6i
, c6j
, c12i
, c12j
, epsi
, epsj
, sigmai
, sigmaj
;
228 snew(nbfp
, 2*atnr
*atnr
);
229 for (i
= 0; i
< atnr
; ++i
)
231 for (j
= 0; j
< atnr
; ++j
)
233 c6i
= idef
->iparams
[i
*(atnr
+1)].lj
.c6
;
234 c12i
= idef
->iparams
[i
*(atnr
+1)].lj
.c12
;
235 c6j
= idef
->iparams
[j
*(atnr
+1)].lj
.c6
;
236 c12j
= idef
->iparams
[j
*(atnr
+1)].lj
.c12
;
237 c6
= std::sqrt(c6i
* c6j
);
238 c12
= std::sqrt(c12i
* c12j
);
239 if (comb_rule
== eCOMB_ARITHMETIC
240 && !gmx_numzero(c6
) && !gmx_numzero(c12
))
242 sigmai
= gmx::sixthroot(c12i
/ c6i
);
243 sigmaj
= gmx::sixthroot(c12j
/ c6j
);
244 epsi
= c6i
* c6i
/ c12i
;
245 epsj
= c6j
* c6j
/ c12j
;
246 c6
= std::sqrt(epsi
* epsj
) * gmx::power6(0.5*(sigmai
+sigmaj
));
247 c12
= std::sqrt(epsi
* epsj
) * gmx::power12(0.5*(sigmai
+sigmaj
));
249 C6(nbfp
, atnr
, i
, j
) = c6
*6.0;
250 C12(nbfp
, atnr
, i
, j
) = c12
*12.0;
256 /* This routine sets fr->solvent_opt to the most common solvent in the
257 * system, e.g. esolSPC or esolTIP4P. It will also mark each charge group in
258 * the fr->solvent_type array with the correct type (or esolNO).
260 * Charge groups that fulfill the conditions but are not identical to the
261 * most common one will be marked as esolNO in the solvent_type array.
263 * TIP3p is identical to SPC for these purposes, so we call it
264 * SPC in the arrays (Apologies to Bill Jorgensen ;-)
266 * NOTE: QM particle should not
267 * become an optimized solvent. Not even if there is only one charge
277 } solvent_parameters_t
;
280 check_solvent_cg(const gmx_moltype_t
*molt
,
283 const unsigned char *qm_grpnr
,
284 const t_grps
*qm_grps
,
286 int *n_solvent_parameters
,
287 solvent_parameters_t
**solvent_parameters_p
,
297 real tmp_charge
[4] = { 0.0 }; /* init to zero to make gcc4.8 happy */
298 int tmp_vdwtype
[4] = { 0 }; /* init to zero to make gcc4.8 happy */
301 solvent_parameters_t
*solvent_parameters
;
303 /* We use a list with parameters for each solvent type.
304 * Every time we discover a new molecule that fulfills the basic
305 * conditions for a solvent we compare with the previous entries
306 * in these lists. If the parameters are the same we just increment
307 * the counter for that type, and otherwise we create a new type
308 * based on the current molecule.
310 * Once we've finished going through all molecules we check which
311 * solvent is most common, and mark all those molecules while we
312 * clear the flag on all others.
315 solvent_parameters
= *solvent_parameters_p
;
317 /* Mark the cg first as non optimized */
320 /* Check if this cg has no exclusions with atoms in other charge groups
321 * and all atoms inside the charge group excluded.
322 * We only have 3 or 4 atom solvent loops.
324 if (GET_CGINFO_EXCL_INTER(cginfo
) ||
325 !GET_CGINFO_EXCL_INTRA(cginfo
))
330 /* Get the indices of the first atom in this charge group */
331 j0
= molt
->cgs
.index
[cg0
];
332 j1
= molt
->cgs
.index
[cg0
+1];
334 /* Number of atoms in our molecule */
340 "Moltype '%s': there are %d atoms in this charge group\n",
344 /* Check if it could be an SPC (3 atoms) or TIP4p (4) water,
347 if (nj
< 3 || nj
> 4)
352 /* Check if we are doing QM on this group */
354 if (qm_grpnr
!= nullptr)
356 for (j
= j0
; j
< j1
&& !qm
; j
++)
358 qm
= (qm_grpnr
[j
] < qm_grps
->nr
- 1);
361 /* Cannot use solvent optimization with QM */
367 atom
= molt
->atoms
.atom
;
369 /* Still looks like a solvent, time to check parameters */
371 /* If it is perturbed (free energy) we can't use the solvent loops,
372 * so then we just skip to the next molecule.
376 for (j
= j0
; j
< j1
&& !perturbed
; j
++)
378 perturbed
= PERTURBED(atom
[j
]);
386 /* Now it's only a question if the VdW and charge parameters
387 * are OK. Before doing the check we compare and see if they are
388 * identical to a possible previous solvent type.
389 * First we assign the current types and charges.
391 for (j
= 0; j
< nj
; j
++)
393 tmp_vdwtype
[j
] = atom
[j0
+j
].type
;
394 tmp_charge
[j
] = atom
[j0
+j
].q
;
397 /* Does it match any previous solvent type? */
398 for (k
= 0; k
< *n_solvent_parameters
; k
++)
403 /* We can only match SPC with 3 atoms and TIP4p with 4 atoms */
404 if ( (solvent_parameters
[k
].model
== esolSPC
&& nj
!= 3) ||
405 (solvent_parameters
[k
].model
== esolTIP4P
&& nj
!= 4) )
410 /* Check that types & charges match for all atoms in molecule */
411 for (j
= 0; j
< nj
&& match
== TRUE
; j
++)
413 if (tmp_vdwtype
[j
] != solvent_parameters
[k
].vdwtype
[j
])
417 if (tmp_charge
[j
] != solvent_parameters
[k
].charge
[j
])
424 /* Congratulations! We have a matched solvent.
425 * Flag it with this type for later processing.
428 solvent_parameters
[k
].count
+= nmol
;
430 /* We are done with this charge group */
435 /* If we get here, we have a tentative new solvent type.
436 * Before we add it we must check that it fulfills the requirements
437 * of the solvent optimized loops. First determine which atoms have
440 for (j
= 0; j
< nj
; j
++)
443 tjA
= tmp_vdwtype
[j
];
445 /* Go through all other tpes and see if any have non-zero
446 * VdW parameters when combined with this one.
448 for (k
= 0; k
< fr
->ntype
&& (has_vdw
[j
] == FALSE
); k
++)
450 /* We already checked that the atoms weren't perturbed,
451 * so we only need to check state A now.
455 has_vdw
[j
] = (has_vdw
[j
] ||
456 (BHAMA(fr
->nbfp
, fr
->ntype
, tjA
, k
) != 0.0) ||
457 (BHAMB(fr
->nbfp
, fr
->ntype
, tjA
, k
) != 0.0) ||
458 (BHAMC(fr
->nbfp
, fr
->ntype
, tjA
, k
) != 0.0));
463 has_vdw
[j
] = (has_vdw
[j
] ||
464 (C6(fr
->nbfp
, fr
->ntype
, tjA
, k
) != 0.0) ||
465 (C12(fr
->nbfp
, fr
->ntype
, tjA
, k
) != 0.0));
470 /* Now we know all we need to make the final check and assignment. */
474 * For this we require thatn all atoms have charge,
475 * the charges on atom 2 & 3 should be the same, and only
476 * atom 1 might have VdW.
478 if (has_vdw
[1] == FALSE
&&
479 has_vdw
[2] == FALSE
&&
480 tmp_charge
[0] != 0 &&
481 tmp_charge
[1] != 0 &&
482 tmp_charge
[2] == tmp_charge
[1])
484 srenew(solvent_parameters
, *n_solvent_parameters
+1);
485 solvent_parameters
[*n_solvent_parameters
].model
= esolSPC
;
486 solvent_parameters
[*n_solvent_parameters
].count
= nmol
;
487 for (k
= 0; k
< 3; k
++)
489 solvent_parameters
[*n_solvent_parameters
].vdwtype
[k
] = tmp_vdwtype
[k
];
490 solvent_parameters
[*n_solvent_parameters
].charge
[k
] = tmp_charge
[k
];
493 *cg_sp
= *n_solvent_parameters
;
494 (*n_solvent_parameters
)++;
499 /* Or could it be a TIP4P?
500 * For this we require thatn atoms 2,3,4 have charge, but not atom 1.
501 * Only atom 1 mght have VdW.
503 if (has_vdw
[1] == FALSE
&&
504 has_vdw
[2] == FALSE
&&
505 has_vdw
[3] == FALSE
&&
506 tmp_charge
[0] == 0 &&
507 tmp_charge
[1] != 0 &&
508 tmp_charge
[2] == tmp_charge
[1] &&
511 srenew(solvent_parameters
, *n_solvent_parameters
+1);
512 solvent_parameters
[*n_solvent_parameters
].model
= esolTIP4P
;
513 solvent_parameters
[*n_solvent_parameters
].count
= nmol
;
514 for (k
= 0; k
< 4; k
++)
516 solvent_parameters
[*n_solvent_parameters
].vdwtype
[k
] = tmp_vdwtype
[k
];
517 solvent_parameters
[*n_solvent_parameters
].charge
[k
] = tmp_charge
[k
];
520 *cg_sp
= *n_solvent_parameters
;
521 (*n_solvent_parameters
)++;
525 *solvent_parameters_p
= solvent_parameters
;
529 check_solvent(FILE * fp
,
530 const gmx_mtop_t
* mtop
,
532 cginfo_mb_t
*cginfo_mb
)
535 const gmx_moltype_t
*molt
;
536 int mb
, mol
, cg_mol
, at_offset
, am
, cgm
, i
, nmol_ch
, nmol
;
537 int n_solvent_parameters
;
538 solvent_parameters_t
*solvent_parameters
;
544 fprintf(debug
, "Going to determine what solvent types we have.\n");
547 n_solvent_parameters
= 0;
548 solvent_parameters
= nullptr;
549 /* Allocate temporary array for solvent type */
550 snew(cg_sp
, mtop
->nmolblock
);
553 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
555 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
557 /* Here we have to loop over all individual molecules
558 * because we need to check for QMMM particles.
560 snew(cg_sp
[mb
], cginfo_mb
[mb
].cg_mod
);
561 nmol_ch
= cginfo_mb
[mb
].cg_mod
/cgs
->nr
;
562 nmol
= mtop
->molblock
[mb
].nmol
/nmol_ch
;
563 for (mol
= 0; mol
< nmol_ch
; mol
++)
566 am
= mol
*cgs
->index
[cgs
->nr
];
567 for (cg_mol
= 0; cg_mol
< cgs
->nr
; cg_mol
++)
569 check_solvent_cg(molt
, cg_mol
, nmol
,
570 mtop
->groups
.grpnr
[egcQMMM
] ?
571 mtop
->groups
.grpnr
[egcQMMM
]+at_offset
+am
: nullptr,
572 &mtop
->groups
.grps
[egcQMMM
],
574 &n_solvent_parameters
, &solvent_parameters
,
575 cginfo_mb
[mb
].cginfo
[cgm
+cg_mol
],
576 &cg_sp
[mb
][cgm
+cg_mol
]);
579 at_offset
+= cgs
->index
[cgs
->nr
];
582 /* Puh! We finished going through all charge groups.
583 * Now find the most common solvent model.
586 /* Most common solvent this far */
588 for (i
= 0; i
< n_solvent_parameters
; i
++)
591 solvent_parameters
[i
].count
> solvent_parameters
[bestsp
].count
)
599 bestsol
= solvent_parameters
[bestsp
].model
;
607 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
609 cgs
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].cgs
;
610 nmol
= (mtop
->molblock
[mb
].nmol
*cgs
->nr
)/cginfo_mb
[mb
].cg_mod
;
611 for (i
= 0; i
< cginfo_mb
[mb
].cg_mod
; i
++)
613 if (cg_sp
[mb
][i
] == bestsp
)
615 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[i
], bestsol
);
620 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[i
], esolNO
);
627 if (bestsol
!= esolNO
&& fp
!= nullptr)
629 fprintf(fp
, "\nEnabling %s-like water optimization for %d molecules.\n\n",
631 solvent_parameters
[bestsp
].count
);
634 sfree(solvent_parameters
);
635 fr
->solvent_opt
= bestsol
;
639 acNONE
= 0, acCONSTRAINT
, acSETTLE
642 static cginfo_mb_t
*init_cginfo_mb(FILE *fplog
, const gmx_mtop_t
*mtop
,
643 t_forcerec
*fr
, gmx_bool bNoSolvOpt
,
644 gmx_bool
*bFEP_NonBonded
,
645 gmx_bool
*bExcl_IntraCGAll_InterCGNone
)
648 const t_blocka
*excl
;
649 const gmx_moltype_t
*molt
;
650 const gmx_molblock_t
*molb
;
651 cginfo_mb_t
*cginfo_mb
;
654 int cg_offset
, a_offset
;
655 int mb
, m
, cg
, a0
, a1
, gid
, ai
, j
, aj
, excl_nalloc
;
659 gmx_bool bId
, *bExcl
, bExclIntraAll
, bExclInter
, bHaveVDW
, bHaveQ
, bHavePerturbedAtoms
;
661 snew(cginfo_mb
, mtop
->nmolblock
);
663 snew(type_VDW
, fr
->ntype
);
664 for (ai
= 0; ai
< fr
->ntype
; ai
++)
666 type_VDW
[ai
] = FALSE
;
667 for (j
= 0; j
< fr
->ntype
; j
++)
669 type_VDW
[ai
] = type_VDW
[ai
] ||
671 C6(fr
->nbfp
, fr
->ntype
, ai
, j
) != 0 ||
672 C12(fr
->nbfp
, fr
->ntype
, ai
, j
) != 0;
676 *bFEP_NonBonded
= FALSE
;
677 *bExcl_IntraCGAll_InterCGNone
= TRUE
;
680 snew(bExcl
, excl_nalloc
);
683 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
685 molb
= &mtop
->molblock
[mb
];
686 molt
= &mtop
->moltype
[molb
->type
];
690 /* Check if the cginfo is identical for all molecules in this block.
691 * If so, we only need an array of the size of one molecule.
692 * Otherwise we make an array of #mol times #cgs per molecule.
695 for (m
= 0; m
< molb
->nmol
; m
++)
697 int am
= m
*cgs
->index
[cgs
->nr
];
698 for (cg
= 0; cg
< cgs
->nr
; cg
++)
701 a1
= cgs
->index
[cg
+1];
702 if (ggrpnr(&mtop
->groups
, egcENER
, a_offset
+am
+a0
) !=
703 ggrpnr(&mtop
->groups
, egcENER
, a_offset
+a0
))
707 if (mtop
->groups
.grpnr
[egcQMMM
] != nullptr)
709 for (ai
= a0
; ai
< a1
; ai
++)
711 if (mtop
->groups
.grpnr
[egcQMMM
][a_offset
+am
+ai
] !=
712 mtop
->groups
.grpnr
[egcQMMM
][a_offset
+ai
])
721 cginfo_mb
[mb
].cg_start
= cg_offset
;
722 cginfo_mb
[mb
].cg_end
= cg_offset
+ molb
->nmol
*cgs
->nr
;
723 cginfo_mb
[mb
].cg_mod
= (bId
? 1 : molb
->nmol
)*cgs
->nr
;
724 snew(cginfo_mb
[mb
].cginfo
, cginfo_mb
[mb
].cg_mod
);
725 cginfo
= cginfo_mb
[mb
].cginfo
;
727 /* Set constraints flags for constrained atoms */
728 snew(a_con
, molt
->atoms
.nr
);
729 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
731 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
736 for (ia
= 0; ia
< molt
->ilist
[ftype
].nr
; ia
+= 1+nral
)
740 for (a
= 0; a
< nral
; a
++)
742 a_con
[molt
->ilist
[ftype
].iatoms
[ia
+1+a
]] =
743 (ftype
== F_SETTLE
? acSETTLE
: acCONSTRAINT
);
749 for (m
= 0; m
< (bId
? 1 : molb
->nmol
); m
++)
752 int am
= m
*cgs
->index
[cgs
->nr
];
753 for (cg
= 0; cg
< cgs
->nr
; cg
++)
756 a1
= cgs
->index
[cg
+1];
758 /* Store the energy group in cginfo */
759 gid
= ggrpnr(&mtop
->groups
, egcENER
, a_offset
+am
+a0
);
760 SET_CGINFO_GID(cginfo
[cgm
+cg
], gid
);
762 /* Check the intra/inter charge group exclusions */
763 if (a1
-a0
> excl_nalloc
)
765 excl_nalloc
= a1
- a0
;
766 srenew(bExcl
, excl_nalloc
);
768 /* bExclIntraAll: all intra cg interactions excluded
769 * bExclInter: any inter cg interactions excluded
771 bExclIntraAll
= TRUE
;
775 bHavePerturbedAtoms
= FALSE
;
776 for (ai
= a0
; ai
< a1
; ai
++)
778 /* Check VDW and electrostatic interactions */
779 bHaveVDW
= bHaveVDW
|| (type_VDW
[molt
->atoms
.atom
[ai
].type
] ||
780 type_VDW
[molt
->atoms
.atom
[ai
].typeB
]);
781 bHaveQ
= bHaveQ
|| (molt
->atoms
.atom
[ai
].q
!= 0 ||
782 molt
->atoms
.atom
[ai
].qB
!= 0);
784 bHavePerturbedAtoms
= bHavePerturbedAtoms
|| (PERTURBED(molt
->atoms
.atom
[ai
]) != 0);
786 /* Clear the exclusion list for atom ai */
787 for (aj
= a0
; aj
< a1
; aj
++)
789 bExcl
[aj
-a0
] = FALSE
;
791 /* Loop over all the exclusions of atom ai */
792 for (j
= excl
->index
[ai
]; j
< excl
->index
[ai
+1]; j
++)
795 if (aj
< a0
|| aj
>= a1
)
804 /* Check if ai excludes a0 to a1 */
805 for (aj
= a0
; aj
< a1
; aj
++)
809 bExclIntraAll
= FALSE
;
816 SET_CGINFO_CONSTR(cginfo
[cgm
+cg
]);
819 SET_CGINFO_SETTLE(cginfo
[cgm
+cg
]);
827 SET_CGINFO_EXCL_INTRA(cginfo
[cgm
+cg
]);
831 SET_CGINFO_EXCL_INTER(cginfo
[cgm
+cg
]);
833 if (a1
- a0
> MAX_CHARGEGROUP_SIZE
)
835 /* The size in cginfo is currently only read with DD */
836 gmx_fatal(FARGS
, "A charge group has size %d which is larger than the limit of %d atoms", a1
-a0
, MAX_CHARGEGROUP_SIZE
);
840 SET_CGINFO_HAS_VDW(cginfo
[cgm
+cg
]);
844 SET_CGINFO_HAS_Q(cginfo
[cgm
+cg
]);
846 if (bHavePerturbedAtoms
&& fr
->efep
!= efepNO
)
848 SET_CGINFO_FEP(cginfo
[cgm
+cg
]);
849 *bFEP_NonBonded
= TRUE
;
851 /* Store the charge group size */
852 SET_CGINFO_NATOMS(cginfo
[cgm
+cg
], a1
-a0
);
854 if (!bExclIntraAll
|| bExclInter
)
856 *bExcl_IntraCGAll_InterCGNone
= FALSE
;
863 cg_offset
+= molb
->nmol
*cgs
->nr
;
864 a_offset
+= molb
->nmol
*cgs
->index
[cgs
->nr
];
868 /* the solvent optimizer is called after the QM is initialized,
869 * because we don't want to have the QM subsystemto become an
873 check_solvent(fplog
, mtop
, fr
, cginfo_mb
);
875 if (getenv("GMX_NO_SOLV_OPT"))
879 fprintf(fplog
, "Found environment variable GMX_NO_SOLV_OPT.\n"
880 "Disabling all solvent optimization\n");
882 fr
->solvent_opt
= esolNO
;
886 fr
->solvent_opt
= esolNO
;
888 if (!fr
->solvent_opt
)
890 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
892 for (cg
= 0; cg
< cginfo_mb
[mb
].cg_mod
; cg
++)
894 SET_CGINFO_SOLOPT(cginfo_mb
[mb
].cginfo
[cg
], esolNO
);
902 static int *cginfo_expand(int nmb
, cginfo_mb_t
*cgi_mb
)
907 ncg
= cgi_mb
[nmb
-1].cg_end
;
910 for (cg
= 0; cg
< ncg
; cg
++)
912 while (cg
>= cgi_mb
[mb
].cg_end
)
917 cgi_mb
[mb
].cginfo
[(cg
- cgi_mb
[mb
].cg_start
) % cgi_mb
[mb
].cg_mod
];
923 /* Sets the sum of charges (squared) and C6 in the system in fr.
924 * Returns whether the system has a net charge.
926 static bool set_chargesum(FILE *log
, t_forcerec
*fr
, const gmx_mtop_t
*mtop
)
928 /*This now calculates sum for q and c6*/
929 double qsum
, q2sum
, q
, c6sum
, c6
;
931 const t_atoms
*atoms
;
936 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
938 nmol
= mtop
->molblock
[mb
].nmol
;
939 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
940 for (i
= 0; i
< atoms
->nr
; i
++)
942 q
= atoms
->atom
[i
].q
;
945 c6
= mtop
->ffparams
.iparams
[atoms
->atom
[i
].type
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
950 fr
->q2sum
[0] = q2sum
;
951 fr
->c6sum
[0] = c6sum
;
953 if (fr
->efep
!= efepNO
)
958 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
960 nmol
= mtop
->molblock
[mb
].nmol
;
961 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
962 for (i
= 0; i
< atoms
->nr
; i
++)
964 q
= atoms
->atom
[i
].qB
;
967 c6
= mtop
->ffparams
.iparams
[atoms
->atom
[i
].typeB
*(mtop
->ffparams
.atnr
+1)].lj
.c6
;
971 fr
->q2sum
[1] = q2sum
;
972 fr
->c6sum
[1] = c6sum
;
977 fr
->qsum
[1] = fr
->qsum
[0];
978 fr
->q2sum
[1] = fr
->q2sum
[0];
979 fr
->c6sum
[1] = fr
->c6sum
[0];
983 if (fr
->efep
== efepNO
)
985 fprintf(log
, "System total charge: %.3f\n", fr
->qsum
[0]);
989 fprintf(log
, "System total charge, top. A: %.3f top. B: %.3f\n",
990 fr
->qsum
[0], fr
->qsum
[1]);
994 /* A cut-off of 1e-4 is used to catch rounding errors due to ascii input */
995 return (std::abs(fr
->qsum
[0]) > 1e-4 ||
996 std::abs(fr
->qsum
[1]) > 1e-4);
999 void update_forcerec(t_forcerec
*fr
, matrix box
)
1001 if (fr
->ic
->eeltype
== eelGRF
)
1003 calc_rffac(nullptr, fr
->ic
->eeltype
, fr
->ic
->epsilon_r
, fr
->ic
->epsilon_rf
,
1004 fr
->ic
->rcoulomb
, fr
->temp
, fr
->zsquare
, box
,
1005 &fr
->ic
->k_rf
, &fr
->ic
->c_rf
);
1009 void set_avcsixtwelve(FILE *fplog
, t_forcerec
*fr
, const gmx_mtop_t
*mtop
)
1011 const t_atoms
*atoms
, *atoms_tpi
;
1012 const t_blocka
*excl
;
1013 int mb
, nmol
, nmolc
, i
, j
, tpi
, tpj
, j1
, j2
, k
, nexcl
, q
;
1014 gmx_int64_t npair
, npair_ij
, tmpi
, tmpj
;
1015 double csix
, ctwelve
;
1016 int ntp
, *typecount
;
1019 real
*nbfp_comb
= nullptr;
1025 /* For LJ-PME, we want to correct for the difference between the
1026 * actual C6 values and the C6 values used by the LJ-PME based on
1027 * combination rules. */
1029 if (EVDW_PME(fr
->ic
->vdwtype
))
1031 nbfp_comb
= mk_nbfp_combination_rule(&mtop
->ffparams
,
1032 (fr
->ljpme_combination_rule
== eljpmeLB
) ? eCOMB_ARITHMETIC
: eCOMB_GEOMETRIC
);
1033 for (tpi
= 0; tpi
< ntp
; ++tpi
)
1035 for (tpj
= 0; tpj
< ntp
; ++tpj
)
1037 C6(nbfp_comb
, ntp
, tpi
, tpj
) =
1038 C6(nbfp
, ntp
, tpi
, tpj
) - C6(nbfp_comb
, ntp
, tpi
, tpj
);
1039 C12(nbfp_comb
, ntp
, tpi
, tpj
) = C12(nbfp
, ntp
, tpi
, tpj
);
1044 for (q
= 0; q
< (fr
->efep
== efepNO
? 1 : 2); q
++)
1052 /* Count the types so we avoid natoms^2 operations */
1053 snew(typecount
, ntp
);
1054 gmx_mtop_count_atomtypes(mtop
, q
, typecount
);
1056 for (tpi
= 0; tpi
< ntp
; tpi
++)
1058 for (tpj
= tpi
; tpj
< ntp
; tpj
++)
1060 tmpi
= typecount
[tpi
];
1061 tmpj
= typecount
[tpj
];
1064 npair_ij
= tmpi
*tmpj
;
1068 npair_ij
= tmpi
*(tmpi
- 1)/2;
1072 /* nbfp now includes the 6.0 derivative prefactor */
1073 csix
+= npair_ij
*BHAMC(nbfp
, ntp
, tpi
, tpj
)/6.0;
1077 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1078 csix
+= npair_ij
* C6(nbfp
, ntp
, tpi
, tpj
)/6.0;
1079 ctwelve
+= npair_ij
* C12(nbfp
, ntp
, tpi
, tpj
)/12.0;
1085 /* Subtract the excluded pairs.
1086 * The main reason for substracting exclusions is that in some cases
1087 * some combinations might never occur and the parameters could have
1088 * any value. These unused values should not influence the dispersion
1091 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1093 nmol
= mtop
->molblock
[mb
].nmol
;
1094 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
1095 excl
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].excls
;
1096 for (i
= 0; (i
< atoms
->nr
); i
++)
1100 tpi
= atoms
->atom
[i
].type
;
1104 tpi
= atoms
->atom
[i
].typeB
;
1106 j1
= excl
->index
[i
];
1107 j2
= excl
->index
[i
+1];
1108 for (j
= j1
; j
< j2
; j
++)
1115 tpj
= atoms
->atom
[k
].type
;
1119 tpj
= atoms
->atom
[k
].typeB
;
1123 /* nbfp now includes the 6.0 derivative prefactor */
1124 csix
-= nmol
*BHAMC(nbfp
, ntp
, tpi
, tpj
)/6.0;
1128 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1129 csix
-= nmol
*C6 (nbfp
, ntp
, tpi
, tpj
)/6.0;
1130 ctwelve
-= nmol
*C12(nbfp
, ntp
, tpi
, tpj
)/12.0;
1140 /* Only correct for the interaction of the test particle
1141 * with the rest of the system.
1144 &mtop
->moltype
[mtop
->molblock
[mtop
->nmolblock
-1].type
].atoms
;
1147 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1149 nmol
= mtop
->molblock
[mb
].nmol
;
1150 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
1151 for (j
= 0; j
< atoms
->nr
; j
++)
1154 /* Remove the interaction of the test charge group
1157 if (mb
== mtop
->nmolblock
-1)
1161 if (mb
== 0 && nmol
== 1)
1163 gmx_fatal(FARGS
, "Old format tpr with TPI, please generate a new tpr file");
1168 tpj
= atoms
->atom
[j
].type
;
1172 tpj
= atoms
->atom
[j
].typeB
;
1174 for (i
= 0; i
< fr
->n_tpi
; i
++)
1178 tpi
= atoms_tpi
->atom
[i
].type
;
1182 tpi
= atoms_tpi
->atom
[i
].typeB
;
1186 /* nbfp now includes the 6.0 derivative prefactor */
1187 csix
+= nmolc
*BHAMC(nbfp
, ntp
, tpi
, tpj
)/6.0;
1191 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1192 csix
+= nmolc
*C6 (nbfp
, ntp
, tpi
, tpj
)/6.0;
1193 ctwelve
+= nmolc
*C12(nbfp
, ntp
, tpi
, tpj
)/12.0;
1200 if (npair
- nexcl
<= 0 && fplog
)
1202 fprintf(fplog
, "\nWARNING: There are no atom pairs for dispersion correction\n\n");
1208 csix
/= npair
- nexcl
;
1209 ctwelve
/= npair
- nexcl
;
1213 fprintf(debug
, "Counted %d exclusions\n", nexcl
);
1214 fprintf(debug
, "Average C6 parameter is: %10g\n", (double)csix
);
1215 fprintf(debug
, "Average C12 parameter is: %10g\n", (double)ctwelve
);
1217 fr
->avcsix
[q
] = csix
;
1218 fr
->avctwelve
[q
] = ctwelve
;
1221 if (EVDW_PME(fr
->ic
->vdwtype
))
1226 if (fplog
!= nullptr)
1228 if (fr
->eDispCorr
== edispcAllEner
||
1229 fr
->eDispCorr
== edispcAllEnerPres
)
1231 fprintf(fplog
, "Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
1232 fr
->avcsix
[0], fr
->avctwelve
[0]);
1236 fprintf(fplog
, "Long Range LJ corr.: <C6> %10.4e\n", fr
->avcsix
[0]);
1242 static real
calcBuckinghamBMax(FILE *fplog
, const gmx_mtop_t
*mtop
)
1244 const t_atoms
*at1
, *at2
;
1245 int mt1
, mt2
, i
, j
, tpi
, tpj
, ntypes
;
1250 fprintf(fplog
, "Determining largest Buckingham b parameter for table\n");
1252 ntypes
= mtop
->ffparams
.atnr
;
1255 real bham_b_max
= 0;
1256 for (mt1
= 0; mt1
< mtop
->nmoltype
; mt1
++)
1258 at1
= &mtop
->moltype
[mt1
].atoms
;
1259 for (i
= 0; (i
< at1
->nr
); i
++)
1261 tpi
= at1
->atom
[i
].type
;
1264 gmx_fatal(FARGS
, "Atomtype[%d] = %d, maximum = %d", i
, tpi
, ntypes
);
1267 for (mt2
= mt1
; mt2
< mtop
->nmoltype
; mt2
++)
1269 at2
= &mtop
->moltype
[mt2
].atoms
;
1270 for (j
= 0; (j
< at2
->nr
); j
++)
1272 tpj
= at2
->atom
[j
].type
;
1275 gmx_fatal(FARGS
, "Atomtype[%d] = %d, maximum = %d", j
, tpj
, ntypes
);
1277 b
= mtop
->ffparams
.iparams
[tpi
*ntypes
+ tpj
].bham
.b
;
1282 if ((b
< bmin
) || (bmin
== -1))
1292 fprintf(fplog
, "Buckingham b parameters, min: %g, max: %g\n",
1299 static void make_nbf_tables(FILE *fp
,
1300 const interaction_const_t
*ic
, real rtab
,
1301 const char *tabfn
, char *eg1
, char *eg2
,
1307 if (tabfn
== nullptr)
1311 fprintf(debug
, "No table file name passed, can not read table, can not do non-bonded interactions\n");
1316 sprintf(buf
, "%s", tabfn
);
1319 /* Append the two energy group names */
1320 sprintf(buf
+ strlen(tabfn
) - strlen(ftp2ext(efXVG
)) - 1, "_%s_%s.%s",
1321 eg1
, eg2
, ftp2ext(efXVG
));
1323 nbl
->table_elec_vdw
= make_tables(fp
, ic
, buf
, rtab
, 0);
1324 /* Copy the contents of the table to separate coulomb and LJ tables too,
1325 * to improve cache performance.
1327 /* For performance reasons we want
1328 * the table data to be aligned to 16-byte. The pointers could be freed
1329 * but currently aren't.
1331 snew(nbl
->table_elec
, 1);
1332 nbl
->table_elec
->interaction
= GMX_TABLE_INTERACTION_ELEC
;
1333 nbl
->table_elec
->format
= nbl
->table_elec_vdw
->format
;
1334 nbl
->table_elec
->r
= nbl
->table_elec_vdw
->r
;
1335 nbl
->table_elec
->n
= nbl
->table_elec_vdw
->n
;
1336 nbl
->table_elec
->scale
= nbl
->table_elec_vdw
->scale
;
1337 nbl
->table_elec
->formatsize
= nbl
->table_elec_vdw
->formatsize
;
1338 nbl
->table_elec
->ninteractions
= 1;
1339 nbl
->table_elec
->stride
= nbl
->table_elec
->formatsize
* nbl
->table_elec
->ninteractions
;
1340 snew_aligned(nbl
->table_elec
->data
, nbl
->table_elec
->stride
*(nbl
->table_elec
->n
+1), 32);
1342 snew(nbl
->table_vdw
, 1);
1343 nbl
->table_vdw
->interaction
= GMX_TABLE_INTERACTION_VDWREP_VDWDISP
;
1344 nbl
->table_vdw
->format
= nbl
->table_elec_vdw
->format
;
1345 nbl
->table_vdw
->r
= nbl
->table_elec_vdw
->r
;
1346 nbl
->table_vdw
->n
= nbl
->table_elec_vdw
->n
;
1347 nbl
->table_vdw
->scale
= nbl
->table_elec_vdw
->scale
;
1348 nbl
->table_vdw
->formatsize
= nbl
->table_elec_vdw
->formatsize
;
1349 nbl
->table_vdw
->ninteractions
= 2;
1350 nbl
->table_vdw
->stride
= nbl
->table_vdw
->formatsize
* nbl
->table_vdw
->ninteractions
;
1351 snew_aligned(nbl
->table_vdw
->data
, nbl
->table_vdw
->stride
*(nbl
->table_vdw
->n
+1), 32);
1353 for (i
= 0; i
<= nbl
->table_elec_vdw
->n
; i
++)
1355 for (j
= 0; j
< 4; j
++)
1357 nbl
->table_elec
->data
[4*i
+j
] = nbl
->table_elec_vdw
->data
[12*i
+j
];
1359 for (j
= 0; j
< 8; j
++)
1361 nbl
->table_vdw
->data
[8*i
+j
] = nbl
->table_elec_vdw
->data
[12*i
+4+j
];
1366 /*!\brief If there's bonded interactions of type \c ftype1 or \c
1367 * ftype2 present in the topology, build an array of the number of
1368 * interactions present for each bonded interaction index found in the
1371 * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
1372 * valid type with that parameter.
1374 * \c count will be reallocated as necessary to fit the largest bonded
1375 * interaction index found, and its current size will be returned in
1376 * \c ncount. It will contain zero for every bonded interaction index
1377 * for which no interactions are present in the topology.
1379 static void count_tables(int ftype1
, int ftype2
, const gmx_mtop_t
*mtop
,
1380 int *ncount
, int **count
)
1382 const gmx_moltype_t
*molt
;
1384 int mt
, ftype
, stride
, i
, j
, tabnr
;
1386 // Loop over all moleculetypes
1387 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1389 molt
= &mtop
->moltype
[mt
];
1390 // Loop over all interaction types
1391 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1393 // If the current interaction type is one of the types whose tables we're trying to count...
1394 if (ftype
== ftype1
|| ftype
== ftype2
)
1396 il
= &molt
->ilist
[ftype
];
1397 stride
= 1 + NRAL(ftype
);
1398 // ... and there are actually some interactions for this type
1399 for (i
= 0; i
< il
->nr
; i
+= stride
)
1401 // Find out which table index the user wanted
1402 tabnr
= mtop
->ffparams
.iparams
[il
->iatoms
[i
]].tab
.table
;
1405 gmx_fatal(FARGS
, "A bonded table number is smaller than 0: %d\n", tabnr
);
1407 // Make room for this index in the data structure
1408 if (tabnr
>= *ncount
)
1410 srenew(*count
, tabnr
+1);
1411 for (j
= *ncount
; j
< tabnr
+1; j
++)
1417 // Record that this table index is used and must have a valid file
1425 /*!\brief If there's bonded interactions of flavour \c tabext and type
1426 * \c ftype1 or \c ftype2 present in the topology, seek them in the
1427 * list of filenames passed to mdrun, and make bonded tables from
1430 * \c ftype1 or \c ftype2 may be set to -1 to disable seeking for a
1431 * valid type with that parameter.
1433 * A fatal error occurs if no matching filename is found.
1435 static bondedtable_t
*make_bonded_tables(FILE *fplog
,
1436 int ftype1
, int ftype2
,
1437 const gmx_mtop_t
*mtop
,
1438 const t_filenm
*tabbfnm
,
1448 count_tables(ftype1
, ftype2
, mtop
, &ncount
, &count
);
1450 // Are there any relevant tabulated bond interactions?
1454 for (int i
= 0; i
< ncount
; i
++)
1456 // Do any interactions exist that requires this table?
1459 // This pattern enforces the current requirement that
1460 // table filenames end in a characteristic sequence
1461 // before the file type extension, and avoids table 13
1462 // being recognized and used for table 1.
1463 std::string patternToFind
= gmx::formatString("_%s%d.%s", tabext
, i
, ftp2ext(efXVG
));
1464 bool madeTable
= false;
1465 for (int j
= 0; j
< tabbfnm
->nfiles
&& !madeTable
; ++j
)
1467 std::string
filename(tabbfnm
->fns
[j
]);
1468 if (gmx::endsWith(filename
, patternToFind
))
1470 // Finally read the table from the file found
1471 tab
[i
] = make_bonded_table(fplog
, tabbfnm
->fns
[j
], NRAL(ftype1
)-2);
1477 bool isPlural
= (ftype2
!= -1);
1478 gmx_fatal(FARGS
, "Tabulated interaction of type '%s%s%s' with index %d cannot be used because no table file whose name matched '%s' was passed via the gmx mdrun -tableb command-line option.",
1479 interaction_function
[ftype1
].longname
,
1480 isPlural
? "' or '" : "",
1481 isPlural
? interaction_function
[ftype2
].longname
: "",
1483 patternToFind
.c_str());
1493 void forcerec_set_ranges(t_forcerec
*fr
,
1494 int ncg_home
, int ncg_force
,
1496 int natoms_force_constr
, int natoms_f_novirsum
)
1501 /* fr->ncg_force is unused in the standard code,
1502 * but it can be useful for modified code dealing with charge groups.
1504 fr
->ncg_force
= ncg_force
;
1505 fr
->natoms_force
= natoms_force
;
1506 fr
->natoms_force_constr
= natoms_force_constr
;
1508 if (fr
->natoms_force_constr
> fr
->nalloc_force
)
1510 fr
->nalloc_force
= over_alloc_dd(fr
->natoms_force_constr
);
1513 if (fr
->haveDirectVirialContributions
)
1515 fr
->forceBufferForDirectVirialContributions
->resize(natoms_f_novirsum
);
1518 if (fr
->ic
->cutoff_scheme
== ecutsVERLET
)
1520 fr
->forceBufferIntermediate
->resize(ncg_home
);
1524 static real
cutoff_inf(real cutoff
)
1528 cutoff
= GMX_CUTOFF_INF
;
1534 gmx_bool
can_use_allvsall(const t_inputrec
*ir
, gmx_bool bPrintNote
, t_commrec
*cr
, FILE *fp
)
1541 ir
->rcoulomb
== 0 &&
1543 ir
->ePBC
== epbcNONE
&&
1544 ir
->vdwtype
== evdwCUT
&&
1545 ir
->coulombtype
== eelCUT
&&
1546 ir
->efep
== efepNO
&&
1547 (ir
->implicit_solvent
== eisNO
||
1548 (ir
->implicit_solvent
== eisGBSA
&& (ir
->gb_algorithm
== egbSTILL
||
1549 ir
->gb_algorithm
== egbHCT
||
1550 ir
->gb_algorithm
== egbOBC
))) &&
1551 getenv("GMX_NO_ALLVSALL") == nullptr
1554 if (bAllvsAll
&& ir
->opts
.ngener
> 1)
1556 const char *note
= "NOTE: Can not use all-vs-all force loops, because there are multiple energy monitor groups; you might get significantly higher performance when using only a single energy monitor group.\n";
1562 fprintf(fp
, "\n%s\n", note
);
1568 if (bAllvsAll
&& fp
&& MASTER(cr
))
1570 fprintf(fp
, "\nUsing SIMD all-vs-all kernels.\n\n");
1577 gmx_bool
nbnxn_simd_supported(const gmx::MDLogger
&mdlog
,
1578 const t_inputrec
*ir
)
1580 if (ir
->vdwtype
== evdwPME
&& ir
->ljpme_combination_rule
== eljpmeLB
)
1582 /* LJ PME with LB combination rule does 7 mesh operations.
1583 * This so slow that we don't compile SIMD non-bonded kernels
1585 GMX_LOG(mdlog
.warning
).asParagraph().appendText("LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling back to plain C kernels");
1593 static void pick_nbnxn_kernel_cpu(const t_inputrec gmx_unused
*ir
,
1597 *kernel_type
= nbnxnk4x4_PlainC
;
1598 *ewald_excl
= ewaldexclTable
;
1602 #ifdef GMX_NBNXN_SIMD_4XN
1603 *kernel_type
= nbnxnk4xN_SIMD_4xN
;
1605 #ifdef GMX_NBNXN_SIMD_2XNN
1606 *kernel_type
= nbnxnk4xN_SIMD_2xNN
;
1609 #if defined GMX_NBNXN_SIMD_2XNN && defined GMX_NBNXN_SIMD_4XN
1610 /* We need to choose if we want 2x(N+N) or 4xN kernels.
1611 * Currently this is based on the SIMD acceleration choice,
1612 * but it might be better to decide this at runtime based on CPU.
1614 * 4xN calculates more (zero) interactions, but has less pair-search
1615 * work and much better kernel instruction scheduling.
1617 * Up till now we have only seen that on Intel Sandy/Ivy Bridge,
1618 * which doesn't have FMA, both the analytical and tabulated Ewald
1619 * kernels have similar pair rates for 4x8 and 2x(4+4), so we choose
1620 * 2x(4+4) because it results in significantly fewer pairs.
1621 * For RF, the raw pair rate of the 4x8 kernel is higher than 2x(4+4),
1622 * 10% with HT, 50% without HT. As we currently don't detect the actual
1623 * use of HT, use 4x8 to avoid a potential performance hit.
1624 * On Intel Haswell 4x8 is always faster.
1626 *kernel_type
= nbnxnk4xN_SIMD_4xN
;
1628 #if !GMX_SIMD_HAVE_FMA
1629 if (EEL_PME_EWALD(ir
->coulombtype
) ||
1630 EVDW_PME(ir
->vdwtype
))
1632 /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
1633 * There are enough instructions to make 2x(4+4) efficient.
1635 *kernel_type
= nbnxnk4xN_SIMD_2xNN
;
1638 #endif /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
1641 if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
1643 #ifdef GMX_NBNXN_SIMD_4XN
1644 *kernel_type
= nbnxnk4xN_SIMD_4xN
;
1646 gmx_fatal(FARGS
, "SIMD 4xN kernels requested, but GROMACS has been compiled without support for these kernels");
1649 if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
1651 #ifdef GMX_NBNXN_SIMD_2XNN
1652 *kernel_type
= nbnxnk4xN_SIMD_2xNN
;
1654 gmx_fatal(FARGS
, "SIMD 2x(N+N) kernels requested, but GROMACS has been compiled without support for these kernels");
1658 /* Analytical Ewald exclusion correction is only an option in
1660 * Since table lookup's don't parallelize with SIMD, analytical
1661 * will probably always be faster for a SIMD width of 8 or more.
1662 * With FMA analytical is sometimes faster for a width if 4 as well.
1663 * On BlueGene/Q, this is faster regardless of precision.
1664 * In single precision, this is faster on Bulldozer.
1665 * On Skylake table is faster in single and double. TODO: Test 5xxx series.
1667 #if ((GMX_SIMD_REAL_WIDTH >= 8 || (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)) \
1668 && !GMX_SIMD_X86_AVX_512) || GMX_SIMD_IBM_QPX
1669 *ewald_excl
= ewaldexclAnalytical
;
1671 if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
1673 *ewald_excl
= ewaldexclTable
;
1675 if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
1677 *ewald_excl
= ewaldexclAnalytical
;
1685 const char *lookup_nbnxn_kernel_name(int kernel_type
)
1687 const char *returnvalue
= nullptr;
1688 switch (kernel_type
)
1691 returnvalue
= "not set";
1693 case nbnxnk4x4_PlainC
:
1694 returnvalue
= "plain C";
1696 case nbnxnk4xN_SIMD_4xN
:
1697 case nbnxnk4xN_SIMD_2xNN
:
1699 returnvalue
= "SIMD";
1701 returnvalue
= "not available";
1704 case nbnxnk8x8x8_GPU
: returnvalue
= "GPU"; break;
1705 case nbnxnk8x8x8_PlainC
: returnvalue
= "plain C"; break;
1709 gmx_fatal(FARGS
, "Illegal kernel type selected");
1710 returnvalue
= nullptr;
1716 static void pick_nbnxn_kernel(FILE *fp
,
1717 const gmx::MDLogger
&mdlog
,
1718 gmx_bool use_simd_kernels
,
1720 EmulateGpuNonbonded emulateGpu
,
1721 const t_inputrec
*ir
,
1724 gmx_bool bDoNonbonded
)
1726 assert(kernel_type
);
1728 *kernel_type
= nbnxnkNotSet
;
1729 *ewald_excl
= ewaldexclTable
;
1731 if (emulateGpu
== EmulateGpuNonbonded::Yes
)
1733 *kernel_type
= nbnxnk8x8x8_PlainC
;
1737 GMX_LOG(mdlog
.warning
).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
1742 *kernel_type
= nbnxnk8x8x8_GPU
;
1745 if (*kernel_type
== nbnxnkNotSet
)
1747 if (use_simd_kernels
&&
1748 nbnxn_simd_supported(mdlog
, ir
))
1750 pick_nbnxn_kernel_cpu(ir
, kernel_type
, ewald_excl
);
1754 *kernel_type
= nbnxnk4x4_PlainC
;
1758 if (bDoNonbonded
&& fp
!= nullptr)
1760 fprintf(fp
, "\nUsing %s %dx%d non-bonded kernels\n\n",
1761 lookup_nbnxn_kernel_name(*kernel_type
),
1762 nbnxn_kernel_to_cluster_i_size(*kernel_type
),
1763 nbnxn_kernel_to_cluster_j_size(*kernel_type
));
1765 if (nbnxnk4x4_PlainC
== *kernel_type
||
1766 nbnxnk8x8x8_PlainC
== *kernel_type
)
1768 GMX_LOG(mdlog
.warning
).asParagraph().appendTextFormatted(
1769 "WARNING: Using the slow %s kernels. This should\n"
1770 "not happen during routine usage on supported platforms.",
1771 lookup_nbnxn_kernel_name(*kernel_type
));
1776 /*! \brief Print Coulomb Ewald citations and set ewald coefficients */
1777 static void initCoulombEwaldParameters(FILE *fp
, const t_inputrec
*ir
,
1778 bool systemHasNetCharge
,
1779 interaction_const_t
*ic
)
1781 if (!EEL_PME_EWALD(ir
->coulombtype
))
1788 fprintf(fp
, "Will do PME sum in reciprocal space for electrostatic interactions.\n");
1790 if (ir
->coulombtype
== eelP3M_AD
)
1792 please_cite(fp
, "Hockney1988");
1793 please_cite(fp
, "Ballenegger2012");
1797 please_cite(fp
, "Essmann95a");
1800 if (ir
->ewald_geometry
== eewg3DC
)
1804 fprintf(fp
, "Using the Ewald3DC correction for systems with a slab geometry%s.\n",
1805 systemHasNetCharge
? " and net charge" : "");
1807 please_cite(fp
, "In-Chul99a");
1808 if (systemHasNetCharge
)
1810 please_cite(fp
, "Ballenegger2009");
1815 ic
->ewaldcoeff_q
= calc_ewaldcoeff_q(ir
->rcoulomb
, ir
->ewald_rtol
);
1818 fprintf(fp
, "Using a Gaussian width (1/beta) of %g nm for Ewald\n",
1819 1/ic
->ewaldcoeff_q
);
1822 if (ic
->coulomb_modifier
== eintmodPOTSHIFT
)
1824 GMX_RELEASE_ASSERT(ic
->rcoulomb
!= 0, "Cutoff radius cannot be zero");
1825 ic
->sh_ewald
= std::erfc(ic
->ewaldcoeff_q
*ic
->rcoulomb
) / ic
->rcoulomb
;
1833 /*! \brief Print Van der Waals Ewald citations and set ewald coefficients */
1834 static void initVdwEwaldParameters(FILE *fp
, const t_inputrec
*ir
,
1835 interaction_const_t
*ic
)
1837 if (!EVDW_PME(ir
->vdwtype
))
1844 fprintf(fp
, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n");
1845 please_cite(fp
, "Essmann95a");
1847 ic
->ewaldcoeff_lj
= calc_ewaldcoeff_lj(ir
->rvdw
, ir
->ewald_rtol_lj
);
1850 fprintf(fp
, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n",
1851 1/ic
->ewaldcoeff_lj
);
1854 if (ic
->vdw_modifier
== eintmodPOTSHIFT
)
1856 real crc2
= gmx::square(ic
->ewaldcoeff_lj
*ic
->rvdw
);
1857 ic
->sh_lj_ewald
= (std::exp(-crc2
)*(1 + crc2
+ 0.5*crc2
*crc2
) - 1)/gmx::power6(ic
->rvdw
);
1861 ic
->sh_lj_ewald
= 0;
1865 gmx_bool
uses_simple_tables(int cutoff_scheme
,
1866 nonbonded_verlet_t
*nbv
,
1869 gmx_bool bUsesSimpleTables
= TRUE
;
1872 switch (cutoff_scheme
)
1875 bUsesSimpleTables
= TRUE
;
1878 assert(NULL
!= nbv
&& NULL
!= nbv
->grp
);
1879 grp_index
= (group
< 0) ? 0 : (nbv
->ngrp
- 1);
1880 bUsesSimpleTables
= nbnxn_kernel_pairlist_simple(nbv
->grp
[grp_index
].kernel_type
);
1883 gmx_incons("unimplemented");
1885 return bUsesSimpleTables
;
1888 static void init_ewald_f_table(interaction_const_t
*ic
,
1893 /* Get the Ewald table spacing based on Coulomb and/or LJ
1894 * Ewald coefficients and rtol.
1896 ic
->tabq_scale
= ewald_spline3_table_scale(ic
);
1898 if (ic
->cutoff_scheme
== ecutsVERLET
)
1900 maxr
= ic
->rcoulomb
;
1904 maxr
= std::max(ic
->rcoulomb
, rtab
);
1906 ic
->tabq_size
= static_cast<int>(maxr
*ic
->tabq_scale
) + 2;
1908 sfree_aligned(ic
->tabq_coul_FDV0
);
1909 sfree_aligned(ic
->tabq_coul_F
);
1910 sfree_aligned(ic
->tabq_coul_V
);
1912 sfree_aligned(ic
->tabq_vdw_FDV0
);
1913 sfree_aligned(ic
->tabq_vdw_F
);
1914 sfree_aligned(ic
->tabq_vdw_V
);
1916 if (EEL_PME_EWALD(ic
->eeltype
))
1918 /* Create the original table data in FDV0 */
1919 snew_aligned(ic
->tabq_coul_FDV0
, ic
->tabq_size
*4, 32);
1920 snew_aligned(ic
->tabq_coul_F
, ic
->tabq_size
, 32);
1921 snew_aligned(ic
->tabq_coul_V
, ic
->tabq_size
, 32);
1922 table_spline3_fill_ewald_lr(ic
->tabq_coul_F
, ic
->tabq_coul_V
, ic
->tabq_coul_FDV0
,
1923 ic
->tabq_size
, 1/ic
->tabq_scale
, ic
->ewaldcoeff_q
, v_q_ewald_lr
);
1926 if (EVDW_PME(ic
->vdwtype
))
1928 snew_aligned(ic
->tabq_vdw_FDV0
, ic
->tabq_size
*4, 32);
1929 snew_aligned(ic
->tabq_vdw_F
, ic
->tabq_size
, 32);
1930 snew_aligned(ic
->tabq_vdw_V
, ic
->tabq_size
, 32);
1931 table_spline3_fill_ewald_lr(ic
->tabq_vdw_F
, ic
->tabq_vdw_V
, ic
->tabq_vdw_FDV0
,
1932 ic
->tabq_size
, 1/ic
->tabq_scale
, ic
->ewaldcoeff_lj
, v_lj_ewald_lr
);
1936 void init_interaction_const_tables(FILE *fp
,
1937 interaction_const_t
*ic
,
1940 if (EEL_PME_EWALD(ic
->eeltype
) || EVDW_PME(ic
->vdwtype
))
1942 init_ewald_f_table(ic
, rtab
);
1946 fprintf(fp
, "Initialized non-bonded Ewald correction tables, spacing: %.2e size: %d\n\n",
1947 1/ic
->tabq_scale
, ic
->tabq_size
);
1952 static void clear_force_switch_constants(shift_consts_t
*sc
)
1959 static void force_switch_constants(real p
,
1963 /* Here we determine the coefficient for shifting the force to zero
1964 * between distance rsw and the cut-off rc.
1965 * For a potential of r^-p, we have force p*r^-(p+1).
1966 * But to save flops we absorb p in the coefficient.
1968 * force/p = r^-(p+1) + c2*r^2 + c3*r^3
1969 * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot
1971 sc
->c2
= ((p
+ 1)*rsw
- (p
+ 4)*rc
)/(pow(rc
, p
+ 2)*gmx::square(rc
- rsw
));
1972 sc
->c3
= -((p
+ 1)*rsw
- (p
+ 3)*rc
)/(pow(rc
, p
+ 2)*gmx::power3(rc
- rsw
));
1973 sc
->cpot
= -pow(rc
, -p
) + p
*sc
->c2
/3*gmx::power3(rc
- rsw
) + p
*sc
->c3
/4*gmx::power4(rc
- rsw
);
1976 static void potential_switch_constants(real rsw
, real rc
,
1977 switch_consts_t
*sc
)
1979 /* The switch function is 1 at rsw and 0 at rc.
1980 * The derivative and second derivate are zero at both ends.
1981 * rsw = max(r - r_switch, 0)
1982 * sw = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
1983 * dsw = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
1984 * force = force*dsw - potential*sw
1987 sc
->c3
= -10/gmx::power3(rc
- rsw
);
1988 sc
->c4
= 15/gmx::power4(rc
- rsw
);
1989 sc
->c5
= -6/gmx::power5(rc
- rsw
);
1992 /*! \brief Construct interaction constants
1994 * This data is used (particularly) by search and force code for
1995 * short-range interactions. Many of these are constant for the whole
1996 * simulation; some are constant only after PME tuning completes.
1999 init_interaction_const(FILE *fp
,
2000 interaction_const_t
**interaction_const
,
2001 const t_inputrec
*ir
,
2002 const gmx_mtop_t
*mtop
,
2003 bool systemHasNetCharge
)
2005 interaction_const_t
*ic
;
2009 ic
->cutoff_scheme
= ir
->cutoff_scheme
;
2011 /* Just allocate something so we can free it */
2012 snew_aligned(ic
->tabq_coul_FDV0
, 16, 32);
2013 snew_aligned(ic
->tabq_coul_F
, 16, 32);
2014 snew_aligned(ic
->tabq_coul_V
, 16, 32);
2017 ic
->vdwtype
= ir
->vdwtype
;
2018 ic
->vdw_modifier
= ir
->vdw_modifier
;
2019 ic
->reppow
= mtop
->ffparams
.reppow
;
2020 ic
->rvdw
= cutoff_inf(ir
->rvdw
);
2021 ic
->rvdw_switch
= ir
->rvdw_switch
;
2022 ic
->ljpme_comb_rule
= ir
->ljpme_combination_rule
;
2023 ic
->useBuckingham
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
2024 if (ic
->useBuckingham
)
2026 ic
->buckinghamBMax
= calcBuckinghamBMax(fp
, mtop
);
2029 initVdwEwaldParameters(fp
, ir
, ic
);
2031 clear_force_switch_constants(&ic
->dispersion_shift
);
2032 clear_force_switch_constants(&ic
->repulsion_shift
);
2034 switch (ic
->vdw_modifier
)
2036 case eintmodPOTSHIFT
:
2037 /* Only shift the potential, don't touch the force */
2038 ic
->dispersion_shift
.cpot
= -1.0/gmx::power6(ic
->rvdw
);
2039 ic
->repulsion_shift
.cpot
= -1.0/gmx::power12(ic
->rvdw
);
2041 case eintmodFORCESWITCH
:
2042 /* Switch the force, switch and shift the potential */
2043 force_switch_constants(6.0, ic
->rvdw_switch
, ic
->rvdw
,
2044 &ic
->dispersion_shift
);
2045 force_switch_constants(12.0, ic
->rvdw_switch
, ic
->rvdw
,
2046 &ic
->repulsion_shift
);
2048 case eintmodPOTSWITCH
:
2049 /* Switch the potential and force */
2050 potential_switch_constants(ic
->rvdw_switch
, ic
->rvdw
,
2054 case eintmodEXACTCUTOFF
:
2055 /* Nothing to do here */
2058 gmx_incons("unimplemented potential modifier");
2061 ic
->sh_invrc6
= -ic
->dispersion_shift
.cpot
;
2063 /* Electrostatics */
2064 ic
->eeltype
= ir
->coulombtype
;
2065 ic
->coulomb_modifier
= ir
->coulomb_modifier
;
2066 ic
->rcoulomb
= cutoff_inf(ir
->rcoulomb
);
2067 ic
->rcoulomb_switch
= ir
->rcoulomb_switch
;
2068 ic
->epsilon_r
= ir
->epsilon_r
;
2070 /* Set the Coulomb energy conversion factor */
2071 if (ic
->epsilon_r
!= 0)
2073 ic
->epsfac
= ONE_4PI_EPS0
/ic
->epsilon_r
;
2077 /* eps = 0 is infinite dieletric: no Coulomb interactions */
2081 /* Reaction-field */
2082 if (EEL_RF(ic
->eeltype
))
2084 ic
->epsilon_rf
= ir
->epsilon_rf
;
2085 /* Generalized reaction field parameters are updated every step */
2086 if (ic
->eeltype
!= eelGRF
)
2088 calc_rffac(fp
, ic
->eeltype
, ic
->epsilon_r
, ic
->epsilon_rf
,
2089 ic
->rcoulomb
, 0, 0, NULL
,
2090 &ic
->k_rf
, &ic
->c_rf
);
2093 if (ir
->cutoff_scheme
== ecutsGROUP
&& ic
->eeltype
== eelRF_ZERO
)
2095 /* grompp should have done this, but this scheme is obsolete */
2096 ic
->coulomb_modifier
= eintmodEXACTCUTOFF
;
2101 /* For plain cut-off we might use the reaction-field kernels */
2102 ic
->epsilon_rf
= ic
->epsilon_r
;
2104 if (ir
->coulomb_modifier
== eintmodPOTSHIFT
)
2106 ic
->c_rf
= 1/ic
->rcoulomb
;
2114 initCoulombEwaldParameters(fp
, ir
, systemHasNetCharge
, ic
);
2118 real dispersion_shift
;
2120 dispersion_shift
= ic
->dispersion_shift
.cpot
;
2121 if (EVDW_PME(ic
->vdwtype
))
2123 dispersion_shift
-= ic
->sh_lj_ewald
;
2125 fprintf(fp
, "Potential shift: LJ r^-12: %.3e r^-6: %.3e",
2126 ic
->repulsion_shift
.cpot
, dispersion_shift
);
2128 if (ic
->eeltype
== eelCUT
)
2130 fprintf(fp
, ", Coulomb %.e", -ic
->c_rf
);
2132 else if (EEL_PME(ic
->eeltype
))
2134 fprintf(fp
, ", Ewald %.3e", -ic
->sh_ewald
);
2139 *interaction_const
= ic
;
2142 /* TODO deviceInfo should be logically const, but currently
2143 * init_gpu modifies it to set up NVML support. This could
2144 * happen during the detection phase, and deviceInfo could
2145 * the become const. */
2146 static void init_nb_verlet(FILE *fp
,
2147 const gmx::MDLogger
&mdlog
,
2148 nonbonded_verlet_t
**nb_verlet
,
2149 gmx_bool bFEP_NonBonded
,
2150 const t_inputrec
*ir
,
2151 const t_forcerec
*fr
,
2152 const t_commrec
*cr
,
2153 gmx_device_info_t
*deviceInfo
,
2154 const gmx_mtop_t
*mtop
,
2157 nonbonded_verlet_t
*nbv
;
2160 nbnxn_alloc_t
*nb_alloc
;
2161 nbnxn_free_t
*nb_free
;
2163 nbv
= new nonbonded_verlet_t();
2165 nbv
->emulateGpu
= ((getenv("GMX_EMULATE_GPU") != nullptr) ? EmulateGpuNonbonded::Yes
: EmulateGpuNonbonded::No
);
2166 nbv
->bUseGPU
= deviceInfo
!= nullptr;
2168 GMX_RELEASE_ASSERT(!(nbv
->emulateGpu
== EmulateGpuNonbonded::Yes
&& nbv
->bUseGPU
), "When GPU emulation is active, there cannot be a GPU assignment");
2172 /* Use the assigned GPU. */
2173 init_gpu(mdlog
, cr
->nodeid
, deviceInfo
);
2177 nbv
->min_ci_balanced
= 0;
2179 nbv
->ngrp
= (DOMAINDECOMP(cr
) ? 2 : 1);
2180 for (int i
= 0; i
< nbv
->ngrp
; i
++)
2182 nbv
->grp
[i
].nbl_lists
.nnbl
= 0;
2183 nbv
->grp
[i
].kernel_type
= nbnxnkNotSet
;
2185 if (i
== 0) /* local */
2187 pick_nbnxn_kernel(fp
, mdlog
, fr
->use_simd_kernels
,
2188 nbv
->bUseGPU
, nbv
->emulateGpu
, ir
,
2189 &nbv
->grp
[i
].kernel_type
,
2190 &nbv
->grp
[i
].ewald_excl
,
2193 else /* non-local */
2195 /* Use the same kernel for local and non-local interactions */
2196 nbv
->grp
[i
].kernel_type
= nbv
->grp
[0].kernel_type
;
2197 nbv
->grp
[i
].ewald_excl
= nbv
->grp
[0].ewald_excl
;
2201 nbv
->listParams
= std::unique_ptr
<NbnxnListParameters
>(new NbnxnListParameters(ir
->rlist
));
2202 setupDynamicPairlistPruning(fp
, ir
, mtop
, box
, nbv
->bUseGPU
, fr
->ic
,
2203 nbv
->listParams
.get());
2205 nbnxn_init_search(&nbv
->nbs
,
2206 DOMAINDECOMP(cr
) ? &cr
->dd
->nc
: nullptr,
2207 DOMAINDECOMP(cr
) ? domdec_zones(cr
->dd
) : nullptr,
2209 gmx_omp_nthreads_get(emntPairsearch
));
2211 gpu_set_host_malloc_and_free(nbv
->grp
[0].kernel_type
== nbnxnk8x8x8_GPU
,
2212 &nb_alloc
, &nb_free
);
2214 for (int i
= 0; i
< nbv
->ngrp
; i
++)
2216 nbnxn_init_pairlist_set(&nbv
->grp
[i
].nbl_lists
,
2217 nbnxn_kernel_pairlist_simple(nbv
->grp
[i
].kernel_type
),
2218 /* 8x8x8 "non-simple" lists are ATM always combined */
2219 !nbnxn_kernel_pairlist_simple(nbv
->grp
[i
].kernel_type
),
2223 int enbnxninitcombrule
;
2224 if (fr
->ic
->vdwtype
== evdwCUT
&&
2225 (fr
->ic
->vdw_modifier
== eintmodNONE
||
2226 fr
->ic
->vdw_modifier
== eintmodPOTSHIFT
) &&
2227 getenv("GMX_NO_LJ_COMB_RULE") == nullptr)
2229 /* Plain LJ cut-off: we can optimize with combination rules */
2230 enbnxninitcombrule
= enbnxninitcombruleDETECT
;
2232 else if (fr
->ic
->vdwtype
== evdwPME
)
2234 /* LJ-PME: we need to use a combination rule for the grid */
2235 if (fr
->ljpme_combination_rule
== eljpmeGEOM
)
2237 enbnxninitcombrule
= enbnxninitcombruleGEOM
;
2241 enbnxninitcombrule
= enbnxninitcombruleLB
;
2246 /* We use a full combination matrix: no rule required */
2247 enbnxninitcombrule
= enbnxninitcombruleNONE
;
2251 bool bSimpleList
= nbnxn_kernel_pairlist_simple(nbv
->grp
[0].kernel_type
);
2252 nbnxn_atomdata_init(fp
,
2254 nbv
->grp
[0].kernel_type
,
2256 fr
->ntype
, fr
->nbfp
,
2258 bSimpleList
? gmx_omp_nthreads_get(emntNonbonded
) : 1,
2263 /* init the NxN GPU data; the last argument tells whether we'll have
2264 * both local and non-local NB calculation on GPU */
2265 nbnxn_gpu_init(&nbv
->gpu_nbv
,
2268 nbv
->listParams
.get(),
2273 /* With tMPI + GPUs some ranks may be sharing GPU(s) and therefore
2274 * also sharing texture references. To keep the code simple, we don't
2275 * treat texture references as shared resources, but this means that
2276 * the coulomb_tab and nbfp texture refs will get updated by multiple threads.
2277 * Hence, to ensure that the non-bonded kernels don't start before all
2278 * texture binding operations are finished, we need to wait for all ranks
2279 * to arrive here before continuing.
2281 * Note that we could omit this barrier if GPUs are not shared (or
2282 * texture objects are used), but as this is initialization code, there
2283 * is no point in complicating things.
2290 #endif /* GMX_THREAD_MPI */
2292 if ((env
= getenv("GMX_NB_MIN_CI")) != nullptr)
2296 nbv
->min_ci_balanced
= strtol(env
, &end
, 10);
2297 if (!end
|| (*end
!= 0) || nbv
->min_ci_balanced
< 0)
2299 gmx_fatal(FARGS
, "Invalid value passed in GMX_NB_MIN_CI=%s, non-negative integer required", env
);
2304 fprintf(debug
, "Neighbor-list balancing parameter: %d (passed as env. var.)\n",
2305 nbv
->min_ci_balanced
);
2310 nbv
->min_ci_balanced
= nbnxn_gpu_min_ci_balanced(nbv
->gpu_nbv
);
2313 fprintf(debug
, "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU multi-processors)\n",
2314 nbv
->min_ci_balanced
);
2323 gmx_bool
usingGpu(nonbonded_verlet_t
*nbv
)
2325 return nbv
!= nullptr && nbv
->bUseGPU
;
2328 void init_forcerec(FILE *fp
,
2329 const gmx::MDLogger
&mdlog
,
2332 const t_inputrec
*ir
,
2333 const gmx_mtop_t
*mtop
,
2334 const t_commrec
*cr
,
2338 const t_filenm
*tabbfnm
,
2339 gmx_device_info_t
*deviceInfo
,
2340 gmx_bool bNoSolvOpt
,
2343 int i
, m
, negp_pp
, negptable
, egi
, egj
;
2348 gmx_bool bGenericKernelOnly
;
2349 gmx_bool needGroupSchemeTables
, bSomeNormalNbListsAreInUse
;
2350 gmx_bool bFEP_NonBonded
;
2351 int *nm_ind
, egp_flags
;
2353 /* By default we turn SIMD kernels on, but it might be turned off further down... */
2354 fr
->use_simd_kernels
= TRUE
;
2356 fr
->bDomDec
= DOMAINDECOMP(cr
);
2358 if (check_box(ir
->ePBC
, box
))
2360 gmx_fatal(FARGS
, check_box(ir
->ePBC
, box
));
2363 /* Test particle insertion ? */
2366 /* Set to the size of the molecule to be inserted (the last one) */
2367 /* Because of old style topologies, we have to use the last cg
2368 * instead of the last molecule type.
2370 cgs
= &mtop
->moltype
[mtop
->molblock
[mtop
->nmolblock
-1].type
].cgs
;
2371 fr
->n_tpi
= cgs
->index
[cgs
->nr
] - cgs
->index
[cgs
->nr
-1];
2372 if (fr
->n_tpi
!= mtop
->mols
.index
[mtop
->mols
.nr
] - mtop
->mols
.index
[mtop
->mols
.nr
-1])
2374 gmx_fatal(FARGS
, "The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group.");
2382 if (ir
->coulombtype
== eelRF_NEC_UNSUPPORTED
)
2384 gmx_fatal(FARGS
, "%s electrostatics is no longer supported",
2385 eel_names
[ir
->coulombtype
]);
2390 gmx_fatal(FARGS
, "AdResS simulations are no longer supported");
2392 if (ir
->useTwinRange
)
2394 gmx_fatal(FARGS
, "Twin-range simulations are no longer supported");
2396 /* Copy the user determined parameters */
2397 fr
->userint1
= ir
->userint1
;
2398 fr
->userint2
= ir
->userint2
;
2399 fr
->userint3
= ir
->userint3
;
2400 fr
->userint4
= ir
->userint4
;
2401 fr
->userreal1
= ir
->userreal1
;
2402 fr
->userreal2
= ir
->userreal2
;
2403 fr
->userreal3
= ir
->userreal3
;
2404 fr
->userreal4
= ir
->userreal4
;
2407 fr
->fc_stepsize
= ir
->fc_stepsize
;
2410 fr
->efep
= ir
->efep
;
2411 fr
->sc_alphavdw
= ir
->fepvals
->sc_alpha
;
2412 if (ir
->fepvals
->bScCoul
)
2414 fr
->sc_alphacoul
= ir
->fepvals
->sc_alpha
;
2415 fr
->sc_sigma6_min
= gmx::power6(ir
->fepvals
->sc_sigma_min
);
2419 fr
->sc_alphacoul
= 0;
2420 fr
->sc_sigma6_min
= 0; /* only needed when bScCoul is on */
2422 fr
->sc_power
= ir
->fepvals
->sc_power
;
2423 fr
->sc_r_power
= ir
->fepvals
->sc_r_power
;
2424 fr
->sc_sigma6_def
= gmx::power6(ir
->fepvals
->sc_sigma
);
2426 env
= getenv("GMX_SCSIGMA_MIN");
2430 sscanf(env
, "%20lf", &dbl
);
2431 fr
->sc_sigma6_min
= gmx::power6(dbl
);
2434 fprintf(fp
, "Setting the minimum soft core sigma to %g nm\n", dbl
);
2438 fr
->bNonbonded
= TRUE
;
2439 if (getenv("GMX_NO_NONBONDED") != nullptr)
2441 /* turn off non-bonded calculations */
2442 fr
->bNonbonded
= FALSE
;
2443 GMX_LOG(mdlog
.warning
).asParagraph().appendText(
2444 "Found environment variable GMX_NO_NONBONDED.\n"
2445 "Disabling nonbonded calculations.");
2448 bGenericKernelOnly
= FALSE
;
2450 /* We now check in the NS code whether a particular combination of interactions
2451 * can be used with water optimization, and disable it if that is not the case.
2454 if (getenv("GMX_NB_GENERIC") != nullptr)
2459 "Found environment variable GMX_NB_GENERIC.\n"
2460 "Disabling all interaction-specific nonbonded kernels, will only\n"
2461 "use the slow generic ones in src/gmxlib/nonbonded/nb_generic.c\n\n");
2463 bGenericKernelOnly
= TRUE
;
2466 if (bGenericKernelOnly
== TRUE
)
2471 if ( (getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr) )
2473 fr
->use_simd_kernels
= FALSE
;
2477 "\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
2478 "Disabling the usage of any SIMD-specific non-bonded & bonded kernel routines\n"
2479 "(e.g. SSE2/SSE4.1/AVX).\n\n");
2483 fr
->bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
2485 /* Check if we can/should do all-vs-all kernels */
2486 fr
->bAllvsAll
= can_use_allvsall(ir
, FALSE
, nullptr, nullptr);
2487 fr
->AllvsAll_work
= nullptr;
2488 fr
->AllvsAll_workgb
= nullptr;
2490 /* All-vs-all kernels have not been implemented in 4.6 and later.
2491 * See Redmine #1249. */
2494 fr
->bAllvsAll
= FALSE
;
2498 "\nYour simulation settings would have triggered the efficient all-vs-all\n"
2499 "kernels in GROMACS 4.5, but these have not been implemented in GROMACS\n"
2500 "4.6 and 5.x. If performance is important, please use GROMACS 4.5.7\n"
2501 "or try cutoff-scheme = Verlet.\n\n");
2505 /* Neighbour searching stuff */
2506 fr
->cutoff_scheme
= ir
->cutoff_scheme
;
2507 fr
->bGrid
= (ir
->ns_type
== ensGRID
);
2508 fr
->ePBC
= ir
->ePBC
;
2510 if (fr
->cutoff_scheme
== ecutsGROUP
)
2512 const char *note
= "NOTE: This file uses the deprecated 'group' cutoff_scheme. This will be\n"
2513 "removed in a future release when 'verlet' supports all interaction forms.\n";
2517 fprintf(stderr
, "\n%s\n", note
);
2521 fprintf(fp
, "\n%s\n", note
);
2526 GMX_LOG(mdlog
.warning
).asParagraph()
2527 .appendText("There is no SIMD implementation of the group scheme kernels on "
2528 "BlueGene/Q. You will observe better performance from using the "
2529 "Verlet cut-off scheme.");
2533 /* Determine if we will do PBC for distances in bonded interactions */
2534 if (fr
->ePBC
== epbcNONE
)
2536 fr
->bMolPBC
= FALSE
;
2540 if (!DOMAINDECOMP(cr
))
2544 bSHAKE
= (ir
->eConstrAlg
== econtSHAKE
&&
2545 (gmx_mtop_ftype_count(mtop
, F_CONSTR
) > 0 ||
2546 gmx_mtop_ftype_count(mtop
, F_CONSTRNC
) > 0));
2548 /* The group cut-off scheme and SHAKE assume charge groups
2549 * are whole, but not using molpbc is faster in most cases.
2550 * With intermolecular interactions we need PBC for calculating
2551 * distances between atoms in different molecules.
2553 if ((fr
->cutoff_scheme
== ecutsGROUP
|| bSHAKE
) &&
2554 !mtop
->bIntermolecularInteractions
)
2556 fr
->bMolPBC
= ir
->bPeriodicMols
;
2558 if (bSHAKE
&& fr
->bMolPBC
)
2560 gmx_fatal(FARGS
, "SHAKE is not supported with periodic molecules");
2565 /* Not making molecules whole is faster in most cases,
2566 * but With orientation restraints we need whole molecules.
2568 fr
->bMolPBC
= (fcd
->orires
.nr
== 0);
2570 if (getenv("GMX_USE_GRAPH") != nullptr)
2572 fr
->bMolPBC
= FALSE
;
2575 GMX_LOG(mdlog
.warning
).asParagraph().appendText("GMX_USE_GRAPH is set, using the graph for bonded interactions");
2578 if (mtop
->bIntermolecularInteractions
)
2580 GMX_LOG(mdlog
.warning
).asParagraph().appendText("WARNING: Molecules linked by intermolecular interactions have to reside in the same periodic image, otherwise artifacts will occur!");
2584 GMX_RELEASE_ASSERT(fr
->bMolPBC
|| !mtop
->bIntermolecularInteractions
, "We need to use PBC within molecules with inter-molecular interactions");
2586 if (bSHAKE
&& fr
->bMolPBC
)
2588 gmx_fatal(FARGS
, "SHAKE is not properly supported with intermolecular interactions. For short simulations where linked molecules remain in the same periodic image, the environment variable GMX_USE_GRAPH can be used to override this check.\n");
2594 fr
->bMolPBC
= dd_bonded_molpbc(cr
->dd
, fr
->ePBC
);
2597 fr
->bGB
= (ir
->implicit_solvent
== eisGBSA
);
2599 fr
->rc_scaling
= ir
->refcoord_scaling
;
2600 copy_rvec(ir
->posres_com
, fr
->posres_com
);
2601 copy_rvec(ir
->posres_comB
, fr
->posres_comB
);
2602 fr
->rlist
= cutoff_inf(ir
->rlist
);
2603 fr
->ljpme_combination_rule
= ir
->ljpme_combination_rule
;
2605 /* This now calculates sum for q and c6*/
2606 bool systemHasNetCharge
= set_chargesum(fp
, fr
, mtop
);
2608 /* fr->ic is used both by verlet and group kernels (to some extent) now */
2609 init_interaction_const(fp
, &fr
->ic
, ir
, mtop
, systemHasNetCharge
);
2610 init_interaction_const_tables(fp
, fr
->ic
, ir
->rlist
+ ir
->tabext
);
2612 const interaction_const_t
*ic
= fr
->ic
;
2614 /* TODO: Replace this Ewald table or move it into interaction_const_t */
2615 if (ir
->coulombtype
== eelEWALD
)
2617 init_ewald_tab(&(fr
->ewald_table
), ir
, fp
);
2620 /* Electrostatics: Translate from interaction-setting-in-mdp-file to kernel interaction format */
2621 switch (ic
->eeltype
)
2624 fr
->nbkernel_elec_interaction
= (fr
->bGB
) ? GMX_NBKERNEL_ELEC_GENERALIZEDBORN
: GMX_NBKERNEL_ELEC_COULOMB
;
2629 fr
->nbkernel_elec_interaction
= GMX_NBKERNEL_ELEC_REACTIONFIELD
;
2633 fr
->nbkernel_elec_interaction
= GMX_NBKERNEL_ELEC_REACTIONFIELD
;
2634 GMX_RELEASE_ASSERT(ic
->coulomb_modifier
== eintmodEXACTCUTOFF
, "With the group scheme RF-zero needs the exact cut-off modifier");
2643 case eelPMEUSERSWITCH
:
2644 fr
->nbkernel_elec_interaction
= GMX_NBKERNEL_ELEC_CUBICSPLINETABLE
;
2650 fr
->nbkernel_elec_interaction
= GMX_NBKERNEL_ELEC_EWALD
;
2654 gmx_fatal(FARGS
, "Unsupported electrostatic interaction: %s", eel_names
[ic
->eeltype
]);
2658 /* Vdw: Translate from mdp settings to kernel format */
2659 switch (ic
->vdwtype
)
2664 fr
->nbkernel_vdw_interaction
= GMX_NBKERNEL_VDW_BUCKINGHAM
;
2668 fr
->nbkernel_vdw_interaction
= GMX_NBKERNEL_VDW_LENNARDJONES
;
2672 fr
->nbkernel_vdw_interaction
= GMX_NBKERNEL_VDW_LJEWALD
;
2678 case evdwENCADSHIFT
:
2679 fr
->nbkernel_vdw_interaction
= GMX_NBKERNEL_VDW_CUBICSPLINETABLE
;
2683 gmx_fatal(FARGS
, "Unsupported vdw interaction: %s", evdw_names
[ic
->vdwtype
]);
2687 if (ir
->cutoff_scheme
== ecutsGROUP
)
2689 fr
->bvdwtab
= ((ic
->vdwtype
!= evdwCUT
|| !gmx_within_tol(ic
->reppow
, 12.0, 10*GMX_DOUBLE_EPS
))
2690 && !EVDW_PME(ic
->vdwtype
));
2691 /* We have special kernels for standard Ewald and PME, but the pme-switch ones are tabulated above */
2692 fr
->bcoultab
= !(ic
->eeltype
== eelCUT
||
2693 ic
->eeltype
== eelEWALD
||
2694 ic
->eeltype
== eelPME
||
2695 ic
->eeltype
== eelP3M_AD
||
2696 ic
->eeltype
== eelRF
||
2697 ic
->eeltype
== eelRF_ZERO
);
2699 /* If the user absolutely wants different switch/shift settings for coul/vdw, it is likely
2700 * going to be faster to tabulate the interaction than calling the generic kernel.
2701 * However, if generic kernels have been requested we keep things analytically.
2703 if (fr
->nbkernel_elec_modifier
== eintmodPOTSWITCH
&&
2704 fr
->nbkernel_vdw_modifier
== eintmodPOTSWITCH
&&
2705 bGenericKernelOnly
== FALSE
)
2707 if ((ic
->rcoulomb_switch
!= ic
->rvdw_switch
) || (ic
->rcoulomb
!= ic
->rvdw
))
2709 fr
->bcoultab
= TRUE
;
2710 /* Once we tabulate electrostatics, we can use the switch function for LJ,
2711 * which would otherwise need two tables.
2715 else if ((fr
->nbkernel_elec_modifier
== eintmodPOTSHIFT
&& fr
->nbkernel_vdw_modifier
== eintmodPOTSHIFT
) ||
2716 ((fr
->nbkernel_elec_interaction
== GMX_NBKERNEL_ELEC_REACTIONFIELD
&&
2717 fr
->nbkernel_elec_modifier
== eintmodEXACTCUTOFF
&&
2718 (fr
->nbkernel_vdw_modifier
== eintmodPOTSWITCH
|| fr
->nbkernel_vdw_modifier
== eintmodPOTSHIFT
))))
2720 if ((ic
->rcoulomb
!= ic
->rvdw
) && (bGenericKernelOnly
== FALSE
))
2722 fr
->bcoultab
= TRUE
;
2726 if (fr
->nbkernel_elec_modifier
== eintmodFORCESWITCH
)
2728 fr
->bcoultab
= TRUE
;
2730 if (fr
->nbkernel_vdw_modifier
== eintmodFORCESWITCH
)
2735 if (getenv("GMX_REQUIRE_TABLES"))
2738 fr
->bcoultab
= TRUE
;
2743 fprintf(fp
, "Table routines are used for coulomb: %s\n",
2744 gmx::boolToString(fr
->bcoultab
));
2745 fprintf(fp
, "Table routines are used for vdw: %s\n",
2746 gmx::boolToString(fr
->bvdwtab
));
2749 if (fr
->bvdwtab
== TRUE
)
2751 fr
->nbkernel_vdw_interaction
= GMX_NBKERNEL_VDW_CUBICSPLINETABLE
;
2752 fr
->nbkernel_vdw_modifier
= eintmodNONE
;
2754 if (fr
->bcoultab
== TRUE
)
2756 fr
->nbkernel_elec_interaction
= GMX_NBKERNEL_ELEC_CUBICSPLINETABLE
;
2757 fr
->nbkernel_elec_modifier
= eintmodNONE
;
2761 if (ir
->cutoff_scheme
== ecutsVERLET
)
2763 if (!gmx_within_tol(ic
->reppow
, 12.0, 10*GMX_DOUBLE_EPS
))
2765 gmx_fatal(FARGS
, "Cut-off scheme %S only supports LJ repulsion power 12", ecutscheme_names
[ir
->cutoff_scheme
]);
2767 fr
->bvdwtab
= FALSE
;
2768 fr
->bcoultab
= FALSE
;
2771 /* 1-4 interaction electrostatics */
2772 fr
->fudgeQQ
= mtop
->ffparams
.fudgeQQ
;
2774 /* Parameters for generalized RF */
2778 if (ic
->eeltype
== eelGRF
)
2780 init_generalized_rf(fp
, mtop
, ir
, fr
);
2783 fr
->haveDirectVirialContributions
=
2784 (EEL_FULL(ic
->eeltype
) || EVDW_PME(ic
->vdwtype
) ||
2785 fr
->forceProviders
->hasForceProvider() ||
2786 gmx_mtop_ftype_count(mtop
, F_POSRES
) > 0 ||
2787 gmx_mtop_ftype_count(mtop
, F_FBPOSRES
) > 0 ||
2792 if (fr
->haveDirectVirialContributions
)
2794 fr
->forceBufferForDirectVirialContributions
= new std::vector
<gmx::RVec
>;
2797 fr
->forceBufferIntermediate
= new std::vector
<gmx::RVec
>; //TODO add proper conditionals
2799 if (fr
->cutoff_scheme
== ecutsGROUP
&&
2800 ncg_mtop(mtop
) > fr
->cg_nalloc
&& !DOMAINDECOMP(cr
))
2802 /* Count the total number of charge groups */
2803 fr
->cg_nalloc
= ncg_mtop(mtop
);
2804 srenew(fr
->cg_cm
, fr
->cg_nalloc
);
2806 if (fr
->shift_vec
== nullptr)
2808 snew(fr
->shift_vec
, SHIFTS
);
2811 if (fr
->fshift
== nullptr)
2813 snew(fr
->fshift
, SHIFTS
);
2816 if (fr
->nbfp
== nullptr)
2818 fr
->ntype
= mtop
->ffparams
.atnr
;
2819 fr
->nbfp
= mk_nbfp(&mtop
->ffparams
, fr
->bBHAM
);
2820 if (EVDW_PME(ic
->vdwtype
))
2822 fr
->ljpme_c6grid
= make_ljpme_c6grid(&mtop
->ffparams
, fr
);
2826 /* Copy the energy group exclusions */
2827 fr
->egp_flags
= ir
->opts
.egp_flags
;
2829 /* Van der Waals stuff */
2830 if ((ic
->vdwtype
!= evdwCUT
) && (ic
->vdwtype
!= evdwUSER
) && !fr
->bBHAM
)
2832 if (ic
->rvdw_switch
>= ic
->rvdw
)
2834 gmx_fatal(FARGS
, "rvdw_switch (%f) must be < rvdw (%f)",
2835 ic
->rvdw_switch
, ic
->rvdw
);
2839 fprintf(fp
, "Using %s Lennard-Jones, switch between %g and %g nm\n",
2840 (ic
->eeltype
== eelSWITCH
) ? "switched" : "shifted",
2841 ic
->rvdw_switch
, ic
->rvdw
);
2845 if (fr
->bBHAM
&& EVDW_PME(ic
->vdwtype
))
2847 gmx_fatal(FARGS
, "LJ PME not supported with Buckingham");
2850 if (fr
->bBHAM
&& (ic
->vdwtype
== evdwSHIFT
|| ic
->vdwtype
== evdwSWITCH
))
2852 gmx_fatal(FARGS
, "Switch/shift interaction not supported with Buckingham");
2855 if (fr
->bBHAM
&& fr
->cutoff_scheme
== ecutsVERLET
)
2857 gmx_fatal(FARGS
, "Verlet cutoff-scheme is not supported with Buckingham");
2862 fprintf(fp
, "Cut-off's: NS: %g Coulomb: %g %s: %g\n",
2863 fr
->rlist
, ic
->rcoulomb
, fr
->bBHAM
? "BHAM" : "LJ", ic
->rvdw
);
2866 fr
->eDispCorr
= ir
->eDispCorr
;
2867 fr
->numAtomsForDispersionCorrection
= mtop
->natoms
;
2868 if (ir
->eDispCorr
!= edispcNO
)
2870 set_avcsixtwelve(fp
, fr
, mtop
);
2873 fr
->gb_epsilon_solvent
= ir
->gb_epsilon_solvent
;
2875 /* Copy the GBSA data (radius, volume and surftens for each
2876 * atomtype) from the topology atomtype section to forcerec.
2878 snew(fr
->atype_radius
, fr
->ntype
);
2879 snew(fr
->atype_vol
, fr
->ntype
);
2880 snew(fr
->atype_surftens
, fr
->ntype
);
2881 snew(fr
->atype_gb_radius
, fr
->ntype
);
2882 snew(fr
->atype_S_hct
, fr
->ntype
);
2884 if (mtop
->atomtypes
.nr
> 0)
2886 for (i
= 0; i
< fr
->ntype
; i
++)
2888 fr
->atype_radius
[i
] = mtop
->atomtypes
.radius
[i
];
2890 for (i
= 0; i
< fr
->ntype
; i
++)
2892 fr
->atype_vol
[i
] = mtop
->atomtypes
.vol
[i
];
2894 for (i
= 0; i
< fr
->ntype
; i
++)
2896 fr
->atype_surftens
[i
] = mtop
->atomtypes
.surftens
[i
];
2898 for (i
= 0; i
< fr
->ntype
; i
++)
2900 fr
->atype_gb_radius
[i
] = mtop
->atomtypes
.gb_radius
[i
];
2902 for (i
= 0; i
< fr
->ntype
; i
++)
2904 fr
->atype_S_hct
[i
] = mtop
->atomtypes
.S_hct
[i
];
2908 /* Generate the GB table if needed */
2912 fr
->gbtabscale
= 2000;
2914 fr
->gbtabscale
= 500;
2918 fr
->gbtab
= make_gb_table(fr
);
2920 init_gb(&fr
->born
, fr
, ir
, mtop
, ir
->gb_algorithm
);
2922 /* Copy local gb data (for dd, this is done in dd_partition_system) */
2923 if (!DOMAINDECOMP(cr
))
2925 make_local_gb(cr
, fr
->born
, ir
->gb_algorithm
);
2929 /* Construct tables for the group scheme. A little unnecessary to
2930 * make both vdw and coul tables sometimes, but what the
2931 * heck. Note that both cutoff schemes construct Ewald tables in
2932 * init_interaction_const_tables. */
2933 needGroupSchemeTables
= (ir
->cutoff_scheme
== ecutsGROUP
&&
2934 (fr
->bcoultab
|| fr
->bvdwtab
));
2936 negp_pp
= ir
->opts
.ngener
- ir
->nwall
;
2938 if (!needGroupSchemeTables
)
2940 bSomeNormalNbListsAreInUse
= TRUE
;
2945 bSomeNormalNbListsAreInUse
= FALSE
;
2946 for (egi
= 0; egi
< negp_pp
; egi
++)
2948 for (egj
= egi
; egj
< negp_pp
; egj
++)
2950 egp_flags
= ir
->opts
.egp_flags
[GID(egi
, egj
, ir
->opts
.ngener
)];
2951 if (!(egp_flags
& EGP_EXCL
))
2953 if (egp_flags
& EGP_TABLE
)
2959 bSomeNormalNbListsAreInUse
= TRUE
;
2964 if (bSomeNormalNbListsAreInUse
)
2966 fr
->nnblists
= negptable
+ 1;
2970 fr
->nnblists
= negptable
;
2972 if (fr
->nnblists
> 1)
2974 snew(fr
->gid2nblists
, ir
->opts
.ngener
*ir
->opts
.ngener
);
2978 snew(fr
->nblists
, fr
->nnblists
);
2980 /* This code automatically gives table length tabext without cut-off's,
2981 * in that case grompp should already have checked that we do not need
2982 * normal tables and we only generate tables for 1-4 interactions.
2984 rtab
= ir
->rlist
+ ir
->tabext
;
2986 if (needGroupSchemeTables
)
2988 /* make tables for ordinary interactions */
2989 if (bSomeNormalNbListsAreInUse
)
2991 make_nbf_tables(fp
, ic
, rtab
, tabfn
, nullptr, nullptr, &fr
->nblists
[0]);
3000 /* Read the special tables for certain energy group pairs */
3001 nm_ind
= mtop
->groups
.grps
[egcENER
].nm_ind
;
3002 for (egi
= 0; egi
< negp_pp
; egi
++)
3004 for (egj
= egi
; egj
< negp_pp
; egj
++)
3006 egp_flags
= ir
->opts
.egp_flags
[GID(egi
, egj
, ir
->opts
.ngener
)];
3007 if ((egp_flags
& EGP_TABLE
) && !(egp_flags
& EGP_EXCL
))
3009 if (fr
->nnblists
> 1)
3011 fr
->gid2nblists
[GID(egi
, egj
, ir
->opts
.ngener
)] = m
;
3013 /* Read the table file with the two energy groups names appended */
3014 make_nbf_tables(fp
, ic
, rtab
, tabfn
,
3015 *mtop
->groups
.grpname
[nm_ind
[egi
]],
3016 *mtop
->groups
.grpname
[nm_ind
[egj
]],
3020 else if (fr
->nnblists
> 1)
3022 fr
->gid2nblists
[GID(egi
, egj
, ir
->opts
.ngener
)] = 0;
3029 /* Tables might not be used for the potential modifier
3030 * interactions per se, but we still need them to evaluate
3031 * switch/shift dispersion corrections in this case. */
3032 if (fr
->eDispCorr
!= edispcNO
)
3034 fr
->dispersionCorrectionTable
= makeDispersionCorrectionTable(fp
, ic
, rtab
, tabfn
);
3037 /* We want to use unmodified tables for 1-4 coulombic
3038 * interactions, so we must in general have an extra set of
3040 if (gmx_mtop_ftype_count(mtop
, F_LJ14
) > 0 ||
3041 gmx_mtop_ftype_count(mtop
, F_LJC14_Q
) > 0 ||
3042 gmx_mtop_ftype_count(mtop
, F_LJC_PAIRS_NB
) > 0)
3044 fr
->pairsTable
= make_tables(fp
, ic
, tabpfn
, rtab
,
3045 GMX_MAKETABLES_14ONLY
);
3049 fr
->nwall
= ir
->nwall
;
3050 if (ir
->nwall
&& ir
->wall_type
== ewtTABLE
)
3052 make_wall_tables(fp
, ir
, tabfn
, &mtop
->groups
, fr
);
3057 // Need to catch std::bad_alloc
3058 // TODO Don't need to catch this here, when merging with master branch
3061 fcd
->bondtab
= make_bonded_tables(fp
,
3062 F_TABBONDS
, F_TABBONDSNC
,
3063 mtop
, tabbfnm
, "b");
3064 fcd
->angletab
= make_bonded_tables(fp
,
3066 mtop
, tabbfnm
, "a");
3067 fcd
->dihtab
= make_bonded_tables(fp
,
3069 mtop
, tabbfnm
, "d");
3071 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3077 fprintf(debug
, "No fcdata or table file name passed, can not read table, can not do bonded interactions\n");
3081 /* QM/MM initialization if requested
3085 fprintf(stderr
, "QM/MM calculation requested.\n");
3088 fr
->bQMMM
= ir
->bQMMM
;
3089 fr
->qr
= mk_QMMMrec();
3091 /* Set all the static charge group info */
3092 fr
->cginfo_mb
= init_cginfo_mb(fp
, mtop
, fr
, bNoSolvOpt
,
3094 &fr
->bExcl_IntraCGAll_InterCGNone
);
3095 if (DOMAINDECOMP(cr
))
3097 fr
->cginfo
= nullptr;
3101 fr
->cginfo
= cginfo_expand(mtop
->nmolblock
, fr
->cginfo_mb
);
3104 if (!DOMAINDECOMP(cr
))
3106 forcerec_set_ranges(fr
, ncg_mtop(mtop
), ncg_mtop(mtop
),
3107 mtop
->natoms
, mtop
->natoms
, mtop
->natoms
);
3110 fr
->print_force
= print_force
;
3113 /* coarse load balancing vars */
3118 /* Initialize neighbor search */
3120 init_ns(fp
, cr
, fr
->ns
, fr
, mtop
);
3122 if (thisRankHasDuty(cr
, DUTY_PP
))
3124 gmx_nonbonded_setup(fr
, bGenericKernelOnly
);
3127 /* Initialize the thread working data for bonded interactions */
3128 init_bonded_threading(fp
, mtop
->groups
.grps
[egcENER
].nr
,
3129 &fr
->bonded_threading
);
3131 fr
->nthread_ewc
= gmx_omp_nthreads_get(emntBonded
);
3132 snew(fr
->ewc_t
, fr
->nthread_ewc
);
3134 if (fr
->cutoff_scheme
== ecutsVERLET
)
3136 // We checked the cut-offs in grompp, but double-check here.
3137 // We have PME+LJcutoff kernels for rcoulomb>rvdw.
3138 if (EEL_PME_EWALD(ir
->coulombtype
) && ir
->vdwtype
== eelCUT
)
3140 GMX_RELEASE_ASSERT(ir
->rcoulomb
>= ir
->rvdw
, "With Verlet lists and PME we should have rcoulomb>=rvdw");
3144 GMX_RELEASE_ASSERT(ir
->rcoulomb
== ir
->rvdw
, "With Verlet lists and no PME rcoulomb and rvdw should be identical");
3147 init_nb_verlet(fp
, mdlog
, &fr
->nbv
, bFEP_NonBonded
, ir
, fr
,
3152 if (ir
->eDispCorr
!= edispcNO
)
3154 calc_enervirdiff(fp
, ir
->eDispCorr
, fr
);
3158 /* Frees GPU memory and destroys the GPU context.
3160 * Note that this function needs to be called even if GPUs are not used
3161 * in this run because the PME ranks have no knowledge of whether GPUs
3162 * are used or not, but all ranks need to enter the barrier below.
3164 void free_gpu_resources(const t_forcerec
*fr
,
3165 const t_commrec
*cr
,
3166 const gmx_device_info_t
*deviceInfo
)
3168 gmx_bool bIsPPrankUsingGPU
;
3169 char gpu_err_str
[STRLEN
];
3171 bIsPPrankUsingGPU
= thisRankHasDuty(cr
, DUTY_PP
) && fr
&& fr
->nbv
&& fr
->nbv
->bUseGPU
;
3173 if (bIsPPrankUsingGPU
)
3175 /* free nbnxn data in GPU memory */
3176 nbnxn_gpu_free(fr
->nbv
->gpu_nbv
);
3177 /* stop the GPU profiler (only CUDA) */
3181 /* With tMPI we need to wait for all ranks to finish deallocation before
3182 * destroying the CUDA context in free_gpu() as some tMPI ranks may be sharing
3185 * This is not a concern in OpenCL where we use one context per rank which
3186 * is freed in nbnxn_gpu_free().
3188 * Note: it is safe to not call the barrier on the ranks which do not use GPU,
3189 * but it is easier and more futureproof to call it on the whole node.
3192 if (PAR(cr
) || MULTISIM(cr
))
3194 gmx_barrier_physical_node(cr
);
3196 #endif /* GMX_THREAD_MPI */
3198 if (bIsPPrankUsingGPU
)
3200 /* uninitialize GPU (by destroying the context) */
3201 if (!free_cuda_gpu(deviceInfo
, gpu_err_str
))
3203 gmx_warning("On rank %d failed to free GPU #%d: %s",
3204 cr
->nodeid
, get_current_cuda_gpu_device_id(), gpu_err_str
);