real Vvdw6, Vvdw12, vvtot;
real ix, iy, iz, fix, fiy, fiz;
real dx, dy, dz, rsq, rinv;
- real c6[NSTATES], c12[NSTATES];
+ real c6[NSTATES], c12[NSTATES], c6grid[NSTATES];
real LFC[NSTATES], LFV[NSTATES], DLF[NSTATES];
double dvdl_coul, dvdl_vdw;
real lfac_coul[NSTATES], dlfac_coul[NSTATES], lfac_vdw[NSTATES], dlfac_vdw[NSTATES];
const real * chargeA;
const real * chargeB;
real sigma6_min, sigma6_def, lam_power, sc_power, sc_r_power;
- real alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc;
- const real * nbfp;
+ real alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc_lj;
+ const real * nbfp, *nbfp_grid;
real * dvdl;
real * Vv;
real * Vc;
real rcutoff, rcutoff2, rswitch, d, d2, swV3, swV4, swV5, swF2, swF3, swF4, sw, dsw, rinvcorr;
const real * tab_ewald_F;
const real * tab_ewald_V;
+ const real * tab_ewald_F_lj;
+ const real * tab_ewald_V_lj;
real tab_ewald_scale, tab_ewald_halfsp;
x = xx[0];
facel = fr->epsfac;
krf = fr->k_rf;
crf = fr->c_rf;
- ewc = fr->ewaldcoeff_q;
+ ewc_lj = fr->ewaldcoeff_lj;
Vc = kernel_data->energygrp_elec;
typeA = mdatoms->typeA;
typeB = mdatoms->typeB;
ntype = fr->ntype;
nbfp = fr->nbfp;
+ nbfp_grid = fr->ljpme_c6grid;
Vv = kernel_data->energygrp_vdw;
lambda_coul = kernel_data->lambda[efptCOUL];
lambda_vdw = kernel_data->lambda[efptVDW];
tab_ewald_F = fr->ic->tabq_coul_F;
tab_ewald_V = fr->ic->tabq_coul_V;
tab_ewald_scale = fr->ic->tabq_scale;
+ tab_ewald_F_lj = fr->ic->tabq_vdw_F;
+ tab_ewald_V_lj = fr->ic->tabq_vdw_V;
tab_ewald_halfsp = 0.5/tab_ewald_scale;
if (fr->coulomb_modifier == eintmodPOTSWITCH || fr->vdw_modifier == eintmodPOTSWITCH)
const interaction_const_t *ic;
ic = fr->ic;
-
- ivdw = GMX_NBKERNEL_VDW_LENNARDJONES;
+ if (EVDW_PME(ic->vdwtype))
+ {
+ ivdw = GMX_NBKERNEL_VDW_LJEWALD;
+ }
+ else
+ {
+ ivdw = GMX_NBKERNEL_VDW_LENNARDJONES;
+ }
if (ic->eeltype == eelCUT || EEL_RF(ic->eeltype))
{
qq[STATE_A] = iqA*chargeA[jnr];
qq[STATE_B] = iqB*chargeB[jnr];
+ tj[STATE_A] = ntiA+2*typeA[jnr];
+ tj[STATE_B] = ntiB+2*typeB[jnr];
+
+ if (ivdw == GMX_NBKERNEL_VDW_LJEWALD)
+ {
+ c6grid[STATE_A] = nbfp_grid[tj[STATE_A]];
+ c6grid[STATE_B] = nbfp_grid[tj[STATE_B]];
+ }
+
if (nlist->excl_fep == NULL || nlist->excl_fep[k])
{
- tj[STATE_A] = ntiA+2*typeA[jnr];
- tj[STATE_B] = ntiB+2*typeB[jnr];
+ c6[STATE_A] = nbfp[tj[STATE_A]];
+ c6[STATE_B] = nbfp[tj[STATE_B]];
for (i = 0; i < NSTATES; i++)
{
-
- c6[i] = nbfp[tj[i]];
c12[i] = nbfp[tj[i]+1];
if ((c6[i] > 0) && (c12[i] > 0))
{
}
if ((c6[i] != 0 || c12[i] != 0) &&
- !(bExactVdwCutoff && rV >= rvdw))
+ !(bExactVdwCutoff &&
+ ((ivdw != GMX_NBKERNEL_VDW_LJEWALD && rV >= rvdw) ||
+ (ivdw == GMX_NBKERNEL_VDW_LJEWALD && r >= rvdw))))
{
switch (ivdw)
{
case GMX_NBKERNEL_VDW_LENNARDJONES:
+ case GMX_NBKERNEL_VDW_LJEWALD:
/* cutoff LJ */
if (sc_r_power == 6.0)
{
}
}
+ if (ivdw == GMX_NBKERNEL_VDW_LJEWALD &&
+ !(bExactVdwCutoff && r >= rvdw))
+ {
+ real rs, frac, f_lr;
+ int ri;
+
+ rs = rsq*rinv*tab_ewald_scale;
+ ri = (int)rs;
+ frac = rs - ri;
+ f_lr = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1];
+ FF = f_lr*rinv;
+ VV = tab_ewald_V_lj[ri] - tab_ewald_halfsp*frac*(tab_ewald_F_lj[ri] + f_lr);
+ for (i = 0; i < NSTATES; i++)
+ {
+ vvtot += LFV[i]*c6grid[i]*VV*(1.0/6.0);
+ Fscal += LFV[i]*c6grid[i]*FF*(1.0/6.0);
+ dvdl_vdw += (DLF[i]*c6grid[i])*VV*(1.0/6.0);
+ }
+
+ }
+
if (bDoForces)
{
tx = Fscal*dx;
int nthread; /* The number of threads doing PME on our rank */
gmx_bool bPPnode; /* Node also does particle-particle forces */
-
gmx_bool bFEP; /* Compute Free energy contribution */
gmx_bool bFEP_q;
gmx_bool bFEP_lj;
-
int nkx, nky, nkz; /* Grid dimensions */
-
gmx_bool bP3M; /* Do P3M: optimize the influence function */
int pme_order;
real epsilon_r;
}
}
-/* TODO: when adding free-energy support, sigmaB will no longer be
- * unused */
int gmx_pme_do(gmx_pme_t pme,
int start, int homenr,
rvec x[], rvec f[],
real *chargeA, real *chargeB,
real *c6A, real *c6B,
- real *sigmaA, real gmx_unused *sigmaB,
+ real *sigmaA, real *sigmaB,
matrix box, t_commrec *cr,
int maxshift_x, int maxshift_y,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
real *dvdlambda_q, real *dvdlambda_lj,
int flags)
{
- int d, i, j, ntot, npme, grid_index, max_grid_index;
+ int d, i, j, k, ntot, npme, grid_index, max_grid_index;
int nx, ny, nz;
int n_d, local_ny;
pme_atomcomm_t *atc = NULL;
t_complex * cfftgrid;
int thread;
gmx_bool bFirst, bDoSplines;
+ int fep_state;
+ int fep_states_lj = pme->bFEP_lj ? 2 : 1;
const gmx_bool bCalcEnerVir = flags & GMX_PME_CALC_ENER_VIR;
const gmx_bool bCalcF = flags & GMX_PME_CALC_F;
if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB)
{
- real *local_c6 = NULL, *local_sigma = NULL;
- if (pme->nnodes == 1)
- {
- if (pme->lb_buf1 == NULL)
- {
- pme->lb_buf_nalloc = pme->atc[0].n;
- snew(pme->lb_buf1, pme->lb_buf_nalloc);
- }
- pme->atc[0].coefficient = pme->lb_buf1;
- local_c6 = c6A;
- local_sigma = sigmaA;
- }
- else
+ /* Loop over A- and B-state if we are doing FEP */
+ for (fep_state = 0; fep_state < fep_states_lj; ++fep_state)
{
- atc = &pme->atc[0];
-
- wallcycle_start(wcycle, ewcPME_REDISTXF);
-
- do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, c6A);
- if (pme->lb_buf_nalloc < atc->n)
- {
- pme->lb_buf_nalloc = atc->nalloc;
- srenew(pme->lb_buf1, pme->lb_buf_nalloc);
- srenew(pme->lb_buf2, pme->lb_buf_nalloc);
- }
- local_c6 = pme->lb_buf1;
- for (i = 0; i < atc->n; ++i)
- {
- local_c6[i] = atc->coefficient[i];
- }
- where();
-
- do_redist_pos_coeffs(pme, cr, start, homenr, FALSE, x, sigmaA);
- local_sigma = pme->lb_buf2;
- for (i = 0; i < atc->n; ++i)
+ real *local_c6 = NULL, *local_sigma = NULL, *RedistC6 = NULL, *RedistSigma = NULL;
+ if (pme->nnodes == 1)
{
- local_sigma[i] = atc->coefficient[i];
+ if (pme->lb_buf1 == NULL)
+ {
+ pme->lb_buf_nalloc = pme->atc[0].n;
+ snew(pme->lb_buf1, pme->lb_buf_nalloc);
+ }
+ pme->atc[0].coefficient = pme->lb_buf1;
+ switch (fep_state)
+ {
+ case 0:
+ local_c6 = c6A;
+ local_sigma = sigmaA;
+ break;
+ case 1:
+ local_c6 = c6B;
+ local_sigma = sigmaB;
+ break;
+ default:
+ gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
+ }
}
- where();
-
- wallcycle_stop(wcycle, ewcPME_REDISTXF);
- }
- calc_initial_lb_coeffs(pme, local_c6, local_sigma);
-
- /*Seven terms in LJ-PME with LB, grid_index < 2 reserved for electrostatics*/
- for (grid_index = 2; grid_index < 9; ++grid_index)
- {
- /* Unpack structure */
- pmegrid = &pme->pmegrid[grid_index];
- fftgrid = pme->fftgrid[grid_index];
- cfftgrid = pme->cfftgrid[grid_index];
- pfft_setup = pme->pfft_setup[grid_index];
- calc_next_lb_coeffs(pme, local_sigma);
- grid = pmegrid->grid.grid;
- where();
-
- if (flags & GMX_PME_SPREAD)
+ else
{
- wallcycle_start(wcycle, ewcPME_SPREADGATHER);
- /* Spread the c6 on a grid */
- spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
+ atc = &pme->atc[0];
+ switch (fep_state)
+ {
+ case 0:
+ RedistC6 = c6A;
+ RedistSigma = sigmaA;
+ break;
+ case 1:
+ RedistC6 = c6B;
+ RedistSigma = sigmaB;
+ break;
+ default:
+ gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
+ }
+ wallcycle_start(wcycle, ewcPME_REDISTXF);
- if (bFirst)
+ do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, RedistC6);
+ if (pme->lb_buf_nalloc < atc->n)
{
- inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
+ pme->lb_buf_nalloc = atc->nalloc;
+ srenew(pme->lb_buf1, pme->lb_buf_nalloc);
+ srenew(pme->lb_buf2, pme->lb_buf_nalloc);
}
- inc_nrnb(nrnb, eNR_SPREADBSP,
- pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
- if (pme->nthread == 1)
+ local_c6 = pme->lb_buf1;
+ for (i = 0; i < atc->n; ++i)
{
- wrap_periodic_pmegrid(pme, grid);
+ local_c6[i] = atc->coefficient[i];
+ }
+ where();
- /* sum contributions to local grid from other nodes */
-#ifdef GMX_MPI
- if (pme->nnodes > 1)
- {
- gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
- where();
- }
-#endif
- copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
+ do_redist_pos_coeffs(pme, cr, start, homenr, FALSE, x, RedistSigma);
+ local_sigma = pme->lb_buf2;
+ for (i = 0; i < atc->n; ++i)
+ {
+ local_sigma[i] = atc->coefficient[i];
}
+ where();
- wallcycle_stop(wcycle, ewcPME_SPREADGATHER);
+ wallcycle_stop(wcycle, ewcPME_REDISTXF);
}
+ calc_initial_lb_coeffs(pme, local_c6, local_sigma);
- /*Here we start a large thread parallel region*/
-#pragma omp parallel num_threads(pme->nthread) private(thread)
+ /*Seven terms in LJ-PME with LB, grid_index < 2 reserved for electrostatics*/
+ for (grid_index = 2; grid_index < 9; ++grid_index)
{
- thread = gmx_omp_get_thread_num();
- if (flags & GMX_PME_SOLVE)
+ /* Unpack structure */
+ pmegrid = &pme->pmegrid[grid_index];
+ fftgrid = pme->fftgrid[grid_index];
+ cfftgrid = pme->cfftgrid[grid_index];
+ pfft_setup = pme->pfft_setup[grid_index];
+ calc_next_lb_coeffs(pme, local_sigma);
+ grid = pmegrid->grid.grid;
+ where();
+
+ if (flags & GMX_PME_SPREAD)
{
- /* do 3d-fft */
- if (thread == 0)
+ wallcycle_start(wcycle, ewcPME_SPREADGATHER);
+ /* Spread the c6 on a grid */
+ spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
+
+ if (bFirst)
{
- wallcycle_start(wcycle, ewcPME_FFT);
+ inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
}
- gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
- thread, wcycle);
- if (thread == 0)
+ inc_nrnb(nrnb, eNR_SPREADBSP,
+ pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
+ if (pme->nthread == 1)
{
- wallcycle_stop(wcycle, ewcPME_FFT);
+ wrap_periodic_pmegrid(pme, grid);
+ /* sum contributions to local grid from other nodes */
+#ifdef GMX_MPI
+ if (pme->nnodes > 1)
+ {
+ gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
+ where();
+ }
+#endif
+ copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
}
- where();
+ wallcycle_stop(wcycle, ewcPME_SPREADGATHER);
}
- }
- bFirst = FALSE;
- }
- if (flags & GMX_PME_SOLVE)
- {
- int loop_count;
- /* solve in k-space for our local cells */
+ /*Here we start a large thread parallel region*/
#pragma omp parallel num_threads(pme->nthread) private(thread)
- {
- thread = gmx_omp_get_thread_num();
- if (thread == 0)
{
- wallcycle_start(wcycle, ewcLJPME);
- }
+ thread = gmx_omp_get_thread_num();
+ if (flags & GMX_PME_SOLVE)
+ {
+ /* do 3d-fft */
+ if (thread == 0)
+ {
+ wallcycle_start(wcycle, ewcPME_FFT);
+ }
- loop_count =
- solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE, ewaldcoeff_lj,
- box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
- bCalcEnerVir,
- pme->nthread, thread);
- if (thread == 0)
- {
- wallcycle_stop(wcycle, ewcLJPME);
- where();
- inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
+ gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
+ thread, wcycle);
+ if (thread == 0)
+ {
+ wallcycle_stop(wcycle, ewcPME_FFT);
+ }
+ where();
+ }
}
+ bFirst = FALSE;
}
- }
-
- if (bCalcEnerVir)
- {
- /* This should only be called on the master thread and
- * after the threads have synchronized.
- */
- get_pme_ener_vir_lj(pme, pme->nthread, &energy_AB[2], vir_AB[2]);
- }
-
- if (bCalcF)
- {
- bFirst = !(flags & GMX_PME_DO_COULOMB);
- calc_initial_lb_coeffs(pme, local_c6, local_sigma);
- for (grid_index = 8; grid_index >= 2; --grid_index)
+ if (flags & GMX_PME_SOLVE)
{
- /* Unpack structure */
- pmegrid = &pme->pmegrid[grid_index];
- fftgrid = pme->fftgrid[grid_index];
- cfftgrid = pme->cfftgrid[grid_index];
- pfft_setup = pme->pfft_setup[grid_index];
- grid = pmegrid->grid.grid;
- calc_next_lb_coeffs(pme, local_sigma);
- where();
+ int loop_count;
+ /* solve in k-space for our local cells */
#pragma omp parallel num_threads(pme->nthread) private(thread)
{
thread = gmx_omp_get_thread_num();
- /* do 3d-invfft */
if (thread == 0)
{
- where();
- wallcycle_start(wcycle, ewcPME_FFT);
+ wallcycle_start(wcycle, ewcLJPME);
}
- gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
- thread, wcycle);
+ loop_count =
+ solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE, ewaldcoeff_lj,
+ box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
+ bCalcEnerVir,
+ pme->nthread, thread);
if (thread == 0)
{
- wallcycle_stop(wcycle, ewcPME_FFT);
-
+ wallcycle_stop(wcycle, ewcLJPME);
where();
+ inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
+ }
+ }
+ }
+
+ if (bCalcEnerVir)
+ {
+ /* This should only be called on the master thread and
+ * after the threads have synchronized.
+ */
+ get_pme_ener_vir_lj(pme, pme->nthread, &energy_AB[2+fep_state], vir_AB[2+fep_state]);
+ }
- if (pme->nodeid == 0)
+ if (bCalcF)
+ {
+ bFirst = !(flags & GMX_PME_DO_COULOMB);
+ calc_initial_lb_coeffs(pme, local_c6, local_sigma);
+ for (grid_index = 8; grid_index >= 2; --grid_index)
+ {
+ /* Unpack structure */
+ pmegrid = &pme->pmegrid[grid_index];
+ fftgrid = pme->fftgrid[grid_index];
+ cfftgrid = pme->cfftgrid[grid_index];
+ pfft_setup = pme->pfft_setup[grid_index];
+ grid = pmegrid->grid.grid;
+ calc_next_lb_coeffs(pme, local_sigma);
+ where();
+#pragma omp parallel num_threads(pme->nthread) private(thread)
+ {
+ thread = gmx_omp_get_thread_num();
+ /* do 3d-invfft */
+ if (thread == 0)
{
- ntot = pme->nkx*pme->nky*pme->nkz;
- npme = ntot*log((real)ntot)/log(2.0);
- inc_nrnb(nrnb, eNR_FFT, 2*npme);
+ where();
+ wallcycle_start(wcycle, ewcPME_FFT);
+ }
+
+ gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
+ thread, wcycle);
+ if (thread == 0)
+ {
+ wallcycle_stop(wcycle, ewcPME_FFT);
+
+ where();
+
+ if (pme->nodeid == 0)
+ {
+ ntot = pme->nkx*pme->nky*pme->nkz;
+ npme = ntot*log((real)ntot)/log(2.0);
+ inc_nrnb(nrnb, eNR_FFT, 2*npme);
+ }
+ wallcycle_start(wcycle, ewcPME_SPREADGATHER);
}
- wallcycle_start(wcycle, ewcPME_SPREADGATHER);
- }
- copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
+ copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
- } /*#pragma omp parallel*/
+ } /*#pragma omp parallel*/
- /* distribute local grid to all nodes */
+ /* distribute local grid to all nodes */
#ifdef GMX_MPI
- if (pme->nnodes > 1)
- {
- gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
- }
+ if (pme->nnodes > 1)
+ {
+ gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
+ }
#endif
- where();
+ where();
- unwrap_periodic_pmegrid(pme, grid);
+ unwrap_periodic_pmegrid(pme, grid);
- /* interpolate forces for our local atoms */
- where();
- bClearF = (bFirst && PAR(cr));
- scale = pme->bFEP ? (grid_index < 9 ? 1.0-lambda_lj : lambda_lj) : 1.0;
- scale *= lb_scale_factor[grid_index-2];
+ /* interpolate forces for our local atoms */
+ where();
+ bClearF = (bFirst && PAR(cr));
+ scale = pme->bFEP ? (fep_state < 1 ? 1.0-lambda_lj : lambda_lj) : 1.0;
+ scale *= lb_scale_factor[grid_index-2];
#pragma omp parallel for num_threads(pme->nthread) schedule(static)
- for (thread = 0; thread < pme->nthread; thread++)
- {
- gather_f_bsplines(pme, grid, bClearF, &pme->atc[0],
- &pme->atc[0].spline[thread],
- scale);
- }
- where();
+ for (thread = 0; thread < pme->nthread; thread++)
+ {
+ gather_f_bsplines(pme, grid, bClearF, &pme->atc[0],
+ &pme->atc[0].spline[thread],
+ scale);
+ }
+ where();
- inc_nrnb(nrnb, eNR_GATHERFBSP,
- pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
- wallcycle_stop(wcycle, ewcPME_SPREADGATHER);
+ inc_nrnb(nrnb, eNR_GATHERFBSP,
+ pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
+ wallcycle_stop(wcycle, ewcPME_SPREADGATHER);
- bFirst = FALSE;
- } /*for (grid_index = 8; grid_index >= 2; --grid_index)*/
- } /*if (bCalcF)*/
- } /*if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB) */
+ bFirst = FALSE;
+ } /* for (grid_index = 8; grid_index >= 2; --grid_index) */
+ } /* if (bCalcF) */
+ } /* for (fep_state = 0; fep_state < fep_states_lj; ++fep_state) */
+ } /* if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB) */
if (bCalcF && pme->nnodes > 1)
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff