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, 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.
39 #include "nb_free_energy.h"
45 #include "gromacs/gmxlib/nrnb.h"
46 #include "gromacs/gmxlib/nonbonded/nb_kernel.h"
47 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
48 #include "gromacs/math/functions.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdtypes/forcerec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/utility/fatalerror.h"
55 gmx_nb_free_energy_kernel(const t_nblist
* gmx_restrict nlist
,
56 rvec
* gmx_restrict xx
,
57 rvec
* gmx_restrict ff
,
58 t_forcerec
* gmx_restrict fr
,
59 const t_mdatoms
* gmx_restrict mdatoms
,
60 nb_kernel_data_t
* gmx_restrict kernel_data
,
61 t_nrnb
* gmx_restrict nrnb
)
67 int i
, n
, ii
, is3
, ii3
, k
, nj0
, nj1
, jnr
, j3
, ggid
;
69 real tx
, ty
, tz
, Fscal
;
70 double FscalC
[NSTATES
], FscalV
[NSTATES
]; /* Needs double for sc_power==48 */
71 double Vcoul
[NSTATES
], Vvdw
[NSTATES
]; /* Needs double for sc_power==48 */
72 real rinv6
, r
, rtC
, rtV
;
74 real qq
[NSTATES
], vctot
;
75 int ntiA
, ntiB
, tj
[NSTATES
];
76 real Vvdw6
, Vvdw12
, vvtot
;
77 real ix
, iy
, iz
, fix
, fiy
, fiz
;
78 real dx
, dy
, dz
, rsq
, rinv
;
79 real c6
[NSTATES
], c12
[NSTATES
], c6grid
;
80 real LFC
[NSTATES
], LFV
[NSTATES
], DLF
[NSTATES
];
81 double dvdl_coul
, dvdl_vdw
;
82 real lfac_coul
[NSTATES
], dlfac_coul
[NSTATES
], lfac_vdw
[NSTATES
], dlfac_vdw
[NSTATES
];
83 real sigma6
[NSTATES
], alpha_vdw_eff
, alpha_coul_eff
, sigma2_def
, sigma2_min
;
84 double rp
, rpm2
, rC
, rV
, rinvC
, rpinvC
, rinvV
, rpinvV
; /* Needs double for sc_power==48 */
85 real sigma2
[NSTATES
], sigma_pow
[NSTATES
];
86 int do_tab
, tab_elemsize
= 0;
87 int n0
, n1C
, n1V
, nnn
;
88 real Y
, F
, Fp
, Geps
, Heps2
, epsC
, eps2C
, epsV
, eps2V
, VV
, FF
;
99 const real
* shiftvec
;
102 const real
* VFtab
= NULL
;
105 real facel
, krf
, crf
;
106 const real
* chargeA
;
107 const real
* chargeB
;
108 real sigma6_min
, sigma6_def
, lam_power
, sc_r_power
;
109 real alpha_coul
, alpha_vdw
, lambda_coul
, lambda_vdw
;
110 real ewcljrsq
, ewclj
, ewclj2
, exponent
, poly
, vvdw_disp
, vvdw_rep
, sh_lj_ewald
;
112 const real
* nbfp
, *nbfp_grid
;
116 gmx_bool bDoForces
, bDoShiftForces
, bDoPotential
;
117 real rcoulomb
, rvdw
, sh_invrc6
;
118 gmx_bool bExactElecCutoff
, bExactVdwCutoff
, bExactCutoffAll
;
119 gmx_bool bEwald
, bEwaldLJ
;
121 const real
* tab_ewald_F_lj
= nullptr;
122 const real
* tab_ewald_V_lj
= nullptr;
123 real d
, d2
, sw
, dsw
, rinvcorr
;
124 real elec_swV3
, elec_swV4
, elec_swV5
, elec_swF2
, elec_swF3
, elec_swF4
;
125 real vdw_swV3
, vdw_swV4
, vdw_swV5
, vdw_swF2
, vdw_swF3
, vdw_swF4
;
126 gmx_bool bConvertEwaldToCoulomb
, bConvertLJEwaldToLJ6
;
127 gmx_bool bComputeVdwInteraction
, bComputeElecInteraction
;
128 const real
* ewtab
= nullptr;
130 real ewrt
, eweps
, ewtabscale
= 0, ewtabhalfspace
= 0, sh_ewald
= 0;
132 const real onetwelfth
= 1.0/12.0;
133 const real onesixth
= 1.0/6.0;
134 const real zero
= 0.0;
135 const real half
= 0.5;
136 const real one
= 1.0;
137 const real two
= 2.0;
138 const real six
= 6.0;
139 const real fourtyeight
= 48.0;
144 fshift
= fr
->fshift
[0];
148 jindex
= nlist
->jindex
;
150 icoul
= nlist
->ielec
;
152 shift
= nlist
->shift
;
155 shiftvec
= fr
->shift_vec
[0];
156 chargeA
= mdatoms
->chargeA
;
157 chargeB
= mdatoms
->chargeB
;
161 Vc
= kernel_data
->energygrp_elec
;
162 typeA
= mdatoms
->typeA
;
163 typeB
= mdatoms
->typeB
;
166 nbfp_grid
= fr
->ljpme_c6grid
;
167 Vv
= kernel_data
->energygrp_vdw
;
168 lambda_coul
= kernel_data
->lambda
[efptCOUL
];
169 lambda_vdw
= kernel_data
->lambda
[efptVDW
];
170 dvdl
= kernel_data
->dvdl
;
171 alpha_coul
= fr
->sc_alphacoul
;
172 alpha_vdw
= fr
->sc_alphavdw
;
173 lam_power
= fr
->sc_power
;
174 sc_r_power
= fr
->sc_r_power
;
175 sigma6_def
= fr
->sc_sigma6_def
;
176 sigma6_min
= fr
->sc_sigma6_min
;
177 bDoForces
= kernel_data
->flags
& GMX_NONBONDED_DO_FORCE
;
178 bDoShiftForces
= kernel_data
->flags
& GMX_NONBONDED_DO_SHIFTFORCE
;
179 bDoPotential
= kernel_data
->flags
& GMX_NONBONDED_DO_POTENTIAL
;
181 rcoulomb
= fr
->rcoulomb
;
183 sh_invrc6
= fr
->ic
->sh_invrc6
;
184 sh_lj_ewald
= fr
->ic
->sh_lj_ewald
;
185 ewclj
= fr
->ewaldcoeff_lj
;
186 ewclj2
= ewclj
*ewclj
;
187 ewclj6
= ewclj2
*ewclj2
*ewclj2
;
189 if (fr
->coulomb_modifier
== eintmodPOTSWITCH
)
191 d
= fr
->rcoulomb
-fr
->rcoulomb_switch
;
192 elec_swV3
= -10.0/(d
*d
*d
);
193 elec_swV4
= 15.0/(d
*d
*d
*d
);
194 elec_swV5
= -6.0/(d
*d
*d
*d
*d
);
195 elec_swF2
= -30.0/(d
*d
*d
);
196 elec_swF3
= 60.0/(d
*d
*d
*d
);
197 elec_swF4
= -30.0/(d
*d
*d
*d
*d
);
201 /* Avoid warnings from stupid compilers (looking at you, Clang!) */
202 elec_swV3
= elec_swV4
= elec_swV5
= elec_swF2
= elec_swF3
= elec_swF4
= 0.0;
205 if (fr
->vdw_modifier
== eintmodPOTSWITCH
)
207 d
= fr
->rvdw
-fr
->rvdw_switch
;
208 vdw_swV3
= -10.0/(d
*d
*d
);
209 vdw_swV4
= 15.0/(d
*d
*d
*d
);
210 vdw_swV5
= -6.0/(d
*d
*d
*d
*d
);
211 vdw_swF2
= -30.0/(d
*d
*d
);
212 vdw_swF3
= 60.0/(d
*d
*d
*d
);
213 vdw_swF4
= -30.0/(d
*d
*d
*d
*d
);
217 /* Avoid warnings from stupid compilers (looking at you, Clang!) */
218 vdw_swV3
= vdw_swV4
= vdw_swV5
= vdw_swF2
= vdw_swF3
= vdw_swF4
= 0.0;
221 if (fr
->cutoff_scheme
== ecutsVERLET
)
223 const interaction_const_t
*ic
;
226 if (EVDW_PME(ic
->vdwtype
))
228 ivdw
= GMX_NBKERNEL_VDW_LJEWALD
;
232 ivdw
= GMX_NBKERNEL_VDW_LENNARDJONES
;
235 if (ic
->eeltype
== eelCUT
|| EEL_RF(ic
->eeltype
))
237 icoul
= GMX_NBKERNEL_ELEC_REACTIONFIELD
;
239 else if (EEL_PME_EWALD(ic
->eeltype
))
241 icoul
= GMX_NBKERNEL_ELEC_EWALD
;
245 gmx_incons("Unsupported eeltype with Verlet and free-energy");
248 bExactElecCutoff
= TRUE
;
249 bExactVdwCutoff
= TRUE
;
253 bExactElecCutoff
= (fr
->coulomb_modifier
!= eintmodNONE
) || fr
->eeltype
== eelRF_ZERO
;
254 bExactVdwCutoff
= (fr
->vdw_modifier
!= eintmodNONE
);
257 bExactCutoffAll
= (bExactElecCutoff
&& bExactVdwCutoff
);
258 rcutoff_max2
= std::max(fr
->rcoulomb
, fr
->rvdw
);
259 rcutoff_max2
= rcutoff_max2
*rcutoff_max2
;
261 bEwald
= (icoul
== GMX_NBKERNEL_ELEC_EWALD
);
262 bEwaldLJ
= (ivdw
== GMX_NBKERNEL_VDW_LJEWALD
);
264 if (bEwald
|| bEwaldLJ
)
266 sh_ewald
= fr
->ic
->sh_ewald
;
267 ewtab
= fr
->ic
->tabq_coul_FDV0
;
268 ewtabscale
= fr
->ic
->tabq_scale
;
269 ewtabhalfspace
= half
/ewtabscale
;
270 tab_ewald_F_lj
= fr
->ic
->tabq_vdw_F
;
271 tab_ewald_V_lj
= fr
->ic
->tabq_vdw_V
;
274 /* For Ewald/PME interactions we cannot easily apply the soft-core component to
275 * reciprocal space. When we use vanilla (not switch/shift) Ewald interactions, we
276 * can apply the small trick of subtracting the _reciprocal_ space contribution
277 * in this kernel, and instead apply the free energy interaction to the 1/r
278 * (standard coulomb) interaction.
280 * However, we cannot use this approach for switch-modified since we would then
281 * effectively end up evaluating a significantly different interaction here compared to the
282 * normal (non-free-energy) kernels, either by applying a cutoff at a different
283 * position than what the user requested, or by switching different
284 * things (1/r rather than short-range Ewald). For these settings, we just
285 * use the traditional short-range Ewald interaction in that case.
287 bConvertEwaldToCoulomb
= (bEwald
&& (fr
->coulomb_modifier
!= eintmodPOTSWITCH
));
288 /* For now the below will always be true (since LJ-PME only works with Shift in Gromacs-5.0),
289 * but writing it this way means we stay in sync with coulomb, and it avoids future bugs.
291 bConvertLJEwaldToLJ6
= (bEwaldLJ
&& (fr
->vdw_modifier
!= eintmodPOTSWITCH
));
293 /* We currently don't implement exclusion correction, needed with the Verlet cut-off scheme, without conversion */
294 if (fr
->cutoff_scheme
== ecutsVERLET
&&
295 ((bEwald
&& !bConvertEwaldToCoulomb
) ||
296 (bEwaldLJ
&& !bConvertLJEwaldToLJ6
)))
298 gmx_incons("Unimplemented non-bonded setup");
301 /* fix compiler warnings */
309 /* Lambda factor for state A, 1-lambda*/
310 LFC
[STATE_A
] = one
- lambda_coul
;
311 LFV
[STATE_A
] = one
- lambda_vdw
;
313 /* Lambda factor for state B, lambda*/
314 LFC
[STATE_B
] = lambda_coul
;
315 LFV
[STATE_B
] = lambda_vdw
;
317 /*derivative of the lambda factor for state A and B */
321 for (i
= 0; i
< NSTATES
; i
++)
323 lfac_coul
[i
] = (lam_power
== 2 ? (1-LFC
[i
])*(1-LFC
[i
]) : (1-LFC
[i
]));
324 dlfac_coul
[i
] = DLF
[i
]*lam_power
/sc_r_power
*(lam_power
== 2 ? (1-LFC
[i
]) : 1);
325 lfac_vdw
[i
] = (lam_power
== 2 ? (1-LFV
[i
])*(1-LFV
[i
]) : (1-LFV
[i
]));
326 dlfac_vdw
[i
] = DLF
[i
]*lam_power
/sc_r_power
*(lam_power
== 2 ? (1-LFV
[i
]) : 1);
329 sigma2_def
= std::cbrt(sigma6_def
);
330 sigma2_min
= std::cbrt(sigma6_min
);
332 /* Ewald (not PME) table is special (icoul==enbcoulFEWALD) */
334 do_tab
= (icoul
== GMX_NBKERNEL_ELEC_CUBICSPLINETABLE
||
335 ivdw
== GMX_NBKERNEL_VDW_CUBICSPLINETABLE
);
338 tabscale
= kernel_data
->table_elec_vdw
->scale
;
339 VFtab
= kernel_data
->table_elec_vdw
->data
;
340 /* we always use the combined table here */
341 tab_elemsize
= kernel_data
->table_elec_vdw
->stride
;
344 for (n
= 0; (n
< nri
); n
++)
346 int npair_within_cutoff
;
348 npair_within_cutoff
= 0;
352 shY
= shiftvec
[is3
+1];
353 shZ
= shiftvec
[is3
+2];
361 iqA
= facel
*chargeA
[ii
];
362 iqB
= facel
*chargeB
[ii
];
363 ntiA
= 2*ntype
*typeA
[ii
];
364 ntiB
= 2*ntype
*typeB
[ii
];
371 for (k
= nj0
; (k
< nj1
); k
++)
378 rsq
= dx
*dx
+ dy
*dy
+ dz
*dz
;
380 if (bExactCutoffAll
&& rsq
>= rcutoff_max2
)
382 /* We save significant time by skipping all code below.
383 * Note that with soft-core interactions, the actual cut-off
384 * check might be different. But since the soft-core distance
385 * is always larger than r, checking on r here is safe.
389 npair_within_cutoff
++;
393 rinv
= gmx::invsqrt(rsq
);
398 /* The force at r=0 is zero, because of symmetry.
399 * But note that the potential is in general non-zero,
400 * since the soft-cored r will be non-zero.
406 if (sc_r_power
== six
)
408 rpm2
= rsq
*rsq
; /* r4 */
409 rp
= rpm2
*rsq
; /* r6 */
411 else if (sc_r_power
== fourtyeight
)
413 rp
= rsq
*rsq
*rsq
; /* r6 */
414 rp
= rp
*rp
; /* r12 */
415 rp
= rp
*rp
; /* r24 */
416 rp
= rp
*rp
; /* r48 */
417 rpm2
= rp
/rsq
; /* r46 */
421 rp
= std::pow(r
, sc_r_power
); /* not currently supported as input, but can handle it */
427 qq
[STATE_A
] = iqA
*chargeA
[jnr
];
428 qq
[STATE_B
] = iqB
*chargeB
[jnr
];
430 tj
[STATE_A
] = ntiA
+2*typeA
[jnr
];
431 tj
[STATE_B
] = ntiB
+2*typeB
[jnr
];
433 if (nlist
->excl_fep
== NULL
|| nlist
->excl_fep
[k
])
435 c6
[STATE_A
] = nbfp
[tj
[STATE_A
]];
436 c6
[STATE_B
] = nbfp
[tj
[STATE_B
]];
438 for (i
= 0; i
< NSTATES
; i
++)
440 c12
[i
] = nbfp
[tj
[i
]+1];
441 if ((c6
[i
] > 0) && (c12
[i
] > 0))
443 /* c12 is stored scaled with 12.0 and c6 is scaled with 6.0 - correct for this */
444 sigma6
[i
] = half
*c12
[i
]/c6
[i
];
445 sigma2
[i
] = std::cbrt(sigma6
[i
]);
446 /* should be able to get rid of cbrt call eventually. Will require agreement on
447 what data to store externally. Can't be fixed without larger scale changes, so not 4.6 */
448 if (sigma6
[i
] < sigma6_min
) /* for disappearing coul and vdw with soft core at the same time */
450 sigma6
[i
] = sigma6_min
;
451 sigma2
[i
] = sigma2_min
;
456 sigma6
[i
] = sigma6_def
;
457 sigma2
[i
] = sigma2_def
;
459 if (sc_r_power
== six
)
461 sigma_pow
[i
] = sigma6
[i
];
463 else if (sc_r_power
== fourtyeight
)
465 sigma_pow
[i
] = sigma6
[i
]*sigma6
[i
]; /* sigma^12 */
466 sigma_pow
[i
] = sigma_pow
[i
]*sigma_pow
[i
]; /* sigma^24 */
467 sigma_pow
[i
] = sigma_pow
[i
]*sigma_pow
[i
]; /* sigma^48 */
470 { /* not really supported as input, but in here for testing the general case*/
471 sigma_pow
[i
] = std::pow(sigma2
[i
], sc_r_power
/2);
475 /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
476 if ((c12
[STATE_A
] > 0) && (c12
[STATE_B
] > 0))
483 alpha_vdw_eff
= alpha_vdw
;
484 alpha_coul_eff
= alpha_coul
;
487 for (i
= 0; i
< NSTATES
; i
++)
494 /* Only spend time on A or B state if it is non-zero */
495 if ( (qq
[i
] != 0) || (c6
[i
] != 0) || (c12
[i
] != 0) )
497 /* this section has to be inside the loop because of the dependence on sigma_pow */
498 rpinvC
= one
/(alpha_coul_eff
*lfac_coul
[i
]*sigma_pow
[i
]+rp
);
499 rinvC
= std::pow(rpinvC
, one
/sc_r_power
);
502 rpinvV
= one
/(alpha_vdw_eff
*lfac_vdw
[i
]*sigma_pow
[i
]+rp
);
503 rinvV
= std::pow(rpinvV
, one
/sc_r_power
);
512 n1C
= tab_elemsize
*n0
;
518 n1V
= tab_elemsize
*n0
;
521 /* Only process the coulomb interactions if we have charges,
522 * and if we either include all entries in the list (no cutoff
523 * used in the kernel), or if we are within the cutoff.
525 bComputeElecInteraction
= !bExactElecCutoff
||
526 ( bConvertEwaldToCoulomb
&& r
< rcoulomb
) ||
527 (!bConvertEwaldToCoulomb
&& rC
< rcoulomb
);
529 if ( (qq
[i
] != 0) && bComputeElecInteraction
)
533 case GMX_NBKERNEL_ELEC_COULOMB
:
535 Vcoul
[i
] = qq
[i
]*rinvC
;
536 FscalC
[i
] = Vcoul
[i
];
537 /* The shift for the Coulomb potential is stored in
538 * the RF parameter c_rf, which is 0 without shift.
540 Vcoul
[i
] -= qq
[i
]*fr
->ic
->c_rf
;
543 case GMX_NBKERNEL_ELEC_REACTIONFIELD
:
545 Vcoul
[i
] = qq
[i
]*(rinvC
+ krf
*rC
*rC
-crf
);
546 FscalC
[i
] = qq
[i
]*(rinvC
- two
*krf
*rC
*rC
);
549 case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE
:
550 /* non-Ewald tabulated coulomb */
554 Geps
= epsC
*VFtab
[nnn
+2];
555 Heps2
= eps2C
*VFtab
[nnn
+3];
558 FF
= Fp
+Geps
+two
*Heps2
;
560 FscalC
[i
] = -qq
[i
]*tabscale
*FF
*rC
;
563 case GMX_NBKERNEL_ELEC_GENERALIZEDBORN
:
564 gmx_fatal(FARGS
, "Free energy and GB not implemented.\n");
567 case GMX_NBKERNEL_ELEC_EWALD
:
568 if (bConvertEwaldToCoulomb
)
570 /* Ewald FEP is done only on the 1/r part */
571 Vcoul
[i
] = qq
[i
]*(rinvC
-sh_ewald
);
572 FscalC
[i
] = qq
[i
]*rinvC
;
576 ewrt
= rC
*ewtabscale
;
577 ewitab
= static_cast<int>(ewrt
);
580 FscalC
[i
] = ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
581 rinvcorr
= rinvC
-sh_ewald
;
582 Vcoul
[i
] = qq
[i
]*(rinvcorr
-(ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+FscalC
[i
])));
583 FscalC
[i
] = qq
[i
]*(rinvC
-rC
*FscalC
[i
]);
587 case GMX_NBKERNEL_ELEC_NONE
:
593 gmx_incons("Invalid icoul in free energy kernel");
597 if (fr
->coulomb_modifier
== eintmodPOTSWITCH
)
599 d
= rC
-fr
->rcoulomb_switch
;
600 d
= (d
> zero
) ? d
: zero
;
602 sw
= one
+d2
*d
*(elec_swV3
+d
*(elec_swV4
+d
*elec_swV5
));
603 dsw
= d2
*(elec_swF2
+d
*(elec_swF3
+d
*elec_swF4
));
605 FscalC
[i
] = FscalC
[i
]*sw
- rC
*Vcoul
[i
]*dsw
;
608 FscalC
[i
] = (rC
< rcoulomb
) ? FscalC
[i
] : zero
;
609 Vcoul
[i
] = (rC
< rcoulomb
) ? Vcoul
[i
] : zero
;
613 /* Only process the VDW interactions if we have
614 * some non-zero parameters, and if we either
615 * include all entries in the list (no cutoff used
616 * in the kernel), or if we are within the cutoff.
618 bComputeVdwInteraction
= !bExactVdwCutoff
||
619 ( bConvertLJEwaldToLJ6
&& r
< rvdw
) ||
620 (!bConvertLJEwaldToLJ6
&& rV
< rvdw
);
621 if ((c6
[i
] != 0 || c12
[i
] != 0) && bComputeVdwInteraction
)
625 case GMX_NBKERNEL_VDW_LENNARDJONES
:
627 if (sc_r_power
== six
)
634 rinv6
= rinv6
*rinv6
*rinv6
;
637 Vvdw12
= c12
[i
]*rinv6
*rinv6
;
639 Vvdw
[i
] = ( (Vvdw12
- c12
[i
]*sh_invrc6
*sh_invrc6
)*onetwelfth
640 - (Vvdw6
- c6
[i
]*sh_invrc6
)*onesixth
);
641 FscalV
[i
] = Vvdw12
- Vvdw6
;
644 case GMX_NBKERNEL_VDW_BUCKINGHAM
:
645 gmx_fatal(FARGS
, "Buckingham free energy not supported.");
648 case GMX_NBKERNEL_VDW_CUBICSPLINETABLE
:
654 Geps
= epsV
*VFtab
[nnn
+2];
655 Heps2
= eps2V
*VFtab
[nnn
+3];
658 FF
= Fp
+Geps
+two
*Heps2
;
660 FscalV
[i
] -= c6
[i
]*tabscale
*FF
*rV
;
665 Geps
= epsV
*VFtab
[nnn
+6];
666 Heps2
= eps2V
*VFtab
[nnn
+7];
669 FF
= Fp
+Geps
+two
*Heps2
;
670 Vvdw
[i
] += c12
[i
]*VV
;
671 FscalV
[i
] -= c12
[i
]*tabscale
*FF
*rV
;
674 case GMX_NBKERNEL_VDW_LJEWALD
:
675 if (sc_r_power
== six
)
682 rinv6
= rinv6
*rinv6
*rinv6
;
684 c6grid
= nbfp_grid
[tj
[i
]];
686 if (bConvertLJEwaldToLJ6
)
690 Vvdw12
= c12
[i
]*rinv6
*rinv6
;
692 Vvdw
[i
] = ( (Vvdw12
- c12
[i
]*sh_invrc6
*sh_invrc6
)*onetwelfth
693 - (Vvdw6
- c6
[i
]*sh_invrc6
- c6grid
*sh_lj_ewald
)*onesixth
);
694 FscalV
[i
] = Vvdw12
- Vvdw6
;
699 ewcljrsq
= ewclj2
*rV
*rV
;
700 exponent
= std::exp(-ewcljrsq
);
701 poly
= exponent
*(one
+ ewcljrsq
+ ewcljrsq
*ewcljrsq
*half
);
702 vvdw_disp
= (c6
[i
]-c6grid
*(one
-poly
))*rinv6
;
703 vvdw_rep
= c12
[i
]*rinv6
*rinv6
;
704 FscalV
[i
] = vvdw_rep
- vvdw_disp
- c6grid
*onesixth
*exponent
*ewclj6
;
705 Vvdw
[i
] = (vvdw_rep
- c12
[i
]*sh_invrc6
*sh_invrc6
)*onetwelfth
- (vvdw_disp
- c6
[i
]*sh_invrc6
- c6grid
*sh_lj_ewald
)/six
;
709 case GMX_NBKERNEL_VDW_NONE
:
715 gmx_incons("Invalid ivdw in free energy kernel");
719 if (fr
->vdw_modifier
== eintmodPOTSWITCH
)
721 d
= rV
-fr
->rvdw_switch
;
722 d
= (d
> zero
) ? d
: zero
;
724 sw
= one
+d2
*d
*(vdw_swV3
+d
*(vdw_swV4
+d
*vdw_swV5
));
725 dsw
= d2
*(vdw_swF2
+d
*(vdw_swF3
+d
*vdw_swF4
));
727 FscalV
[i
] = FscalV
[i
]*sw
- rV
*Vvdw
[i
]*dsw
;
730 FscalV
[i
] = (rV
< rvdw
) ? FscalV
[i
] : zero
;
731 Vvdw
[i
] = (rV
< rvdw
) ? Vvdw
[i
] : zero
;
735 /* FscalC (and FscalV) now contain: dV/drC * rC
736 * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p
737 * Further down we first multiply by r^p-2 and then by
738 * the vector r, which in total gives: dV/drC * (r/rC)^1-p
745 /* Assemble A and B states */
746 for (i
= 0; i
< NSTATES
; i
++)
748 vctot
+= LFC
[i
]*Vcoul
[i
];
749 vvtot
+= LFV
[i
]*Vvdw
[i
];
751 Fscal
+= LFC
[i
]*FscalC
[i
]*rpm2
;
752 Fscal
+= LFV
[i
]*FscalV
[i
]*rpm2
;
754 dvdl_coul
+= Vcoul
[i
]*DLF
[i
] + LFC
[i
]*alpha_coul_eff
*dlfac_coul
[i
]*FscalC
[i
]*sigma_pow
[i
];
755 dvdl_vdw
+= Vvdw
[i
]*DLF
[i
] + LFV
[i
]*alpha_vdw_eff
*dlfac_vdw
[i
]*FscalV
[i
]*sigma_pow
[i
];
758 else if (icoul
== GMX_NBKERNEL_ELEC_REACTIONFIELD
)
760 /* For excluded pairs, which are only in this pair list when
761 * using the Verlet scheme, we don't use soft-core.
762 * The group scheme also doesn't soft-core for these.
763 * As there is no singularity, there is no need for soft-core.
773 for (i
= 0; i
< NSTATES
; i
++)
775 vctot
+= LFC
[i
]*qq
[i
]*VV
;
776 Fscal
+= LFC
[i
]*qq
[i
]*FF
;
777 dvdl_coul
+= DLF
[i
]*qq
[i
]*VV
;
781 if (bConvertEwaldToCoulomb
&& ( !bExactElecCutoff
|| r
< rcoulomb
) )
783 /* See comment in the preamble. When using Ewald interactions
784 * (unless we use a switch modifier) we subtract the reciprocal-space
785 * Ewald component here which made it possible to apply the free
786 * energy interaction to 1/r (vanilla coulomb short-range part)
787 * above. This gets us closer to the ideal case of applying
788 * the softcore to the entire electrostatic interaction,
789 * including the reciprocal-space component.
794 ewitab
= static_cast<int>(ewrt
);
797 f_lr
= ewtab
[ewitab
]+eweps
*ewtab
[ewitab
+1];
798 v_lr
= (ewtab
[ewitab
+2]-ewtabhalfspace
*eweps
*(ewtab
[ewitab
]+f_lr
));
801 /* Note that any possible Ewald shift has already been applied in
802 * the normal interaction part above.
807 /* If we get here, the i particle (ii) has itself (jnr)
808 * in its neighborlist. This can only happen with the Verlet
809 * scheme, and corresponds to a self-interaction that will
810 * occur twice. Scale it down by 50% to only include it once.
815 for (i
= 0; i
< NSTATES
; i
++)
817 vctot
-= LFC
[i
]*qq
[i
]*v_lr
;
818 Fscal
-= LFC
[i
]*qq
[i
]*f_lr
;
819 dvdl_coul
-= (DLF
[i
]*qq
[i
])*v_lr
;
823 if (bConvertLJEwaldToLJ6
&& (!bExactVdwCutoff
|| r
< rvdw
))
825 /* See comment in the preamble. When using LJ-Ewald interactions
826 * (unless we use a switch modifier) we subtract the reciprocal-space
827 * Ewald component here which made it possible to apply the free
828 * energy interaction to r^-6 (vanilla LJ6 short-range part)
829 * above. This gets us closer to the ideal case of applying
830 * the softcore to the entire VdW interaction,
831 * including the reciprocal-space component.
833 /* We could also use the analytical form here
834 * iso a table, but that can cause issues for
835 * r close to 0 for non-interacting pairs.
840 rs
= rsq
*rinv
*ewtabscale
;
841 ri
= static_cast<int>(rs
);
843 f_lr
= (1 - frac
)*tab_ewald_F_lj
[ri
] + frac
*tab_ewald_F_lj
[ri
+1];
844 /* TODO: Currently the Ewald LJ table does not contain
845 * the factor 1/6, we should add this.
848 VV
= (tab_ewald_V_lj
[ri
] - ewtabhalfspace
*frac
*(tab_ewald_F_lj
[ri
] + f_lr
))/six
;
852 /* If we get here, the i particle (ii) has itself (jnr)
853 * in its neighborlist. This can only happen with the Verlet
854 * scheme, and corresponds to a self-interaction that will
855 * occur twice. Scale it down by 50% to only include it once.
860 for (i
= 0; i
< NSTATES
; i
++)
862 c6grid
= nbfp_grid
[tj
[i
]];
863 vvtot
+= LFV
[i
]*c6grid
*VV
;
864 Fscal
+= LFV
[i
]*c6grid
*FF
;
865 dvdl_vdw
+= (DLF
[i
]*c6grid
)*VV
;
877 /* OpenMP atomics are expensive, but this kernels is also
878 * expensive, so we can take this hit, instead of using
879 * thread-local output buffers and extra reduction.
881 * All the OpenMP regions in this file are trivial and should
882 * not throw, so no need for try/catch.
893 /* The atomics below are expensive with many OpenMP threads.
894 * Here unperturbed i-particles will usually only have a few
895 * (perturbed) j-particles in the list. Thus with a buffered list
896 * we can skip a significant number of i-reductions with a check.
898 if (npair_within_cutoff
> 0)
914 fshift
[is3
+1] += fiy
;
916 fshift
[is3
+2] += fiz
;
930 dvdl
[efptCOUL
] += dvdl_coul
;
932 dvdl
[efptVDW
] += dvdl_vdw
;
934 /* Estimate flops, average for free energy stuff:
935 * 12 flops per outer iteration
936 * 150 flops per inner iteration
939 inc_nrnb(nrnb
, eNR_NBKERNEL_FREE_ENERGY
, nlist
->nri
*12 + nlist
->jindex
[n
]*150);