#pragma omp atomic
inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150);
}
-
-real
-nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, real alpha_vdw,
- real tabscale, real *vftab,
- real qqA, real c6A, real c12A, real qqB, real c6B, real c12B,
- real LFC[2], real LFV[2], real DLF[2],
- real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], real dlfac_vdw[2],
- real sigma6_def, real sigma6_min, real sigma2_def, real sigma2_min,
- real *velectot, real *vvdwtot, real *dvdl)
-{
- real r, rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
- real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2], sigma_powm2[2];
- real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
- real rpinv, r_coul, r_vdw, velecsum, vvdwsum;
- real fscal_vdw[2], fscal_elec[2];
- real velec[2], vvdw[2];
- int i, ntab;
-
- const real half = 0.5;
- const real one = 1.0;
- const real two = 2.0;
- const real six = 6.0;
- const real fourtyeight = 48.0;
-
- qq[0] = qqA;
- qq[1] = qqB;
- c6[0] = c6A;
- c6[1] = c6B;
- c12[0] = c12A;
- c12[1] = c12B;
-
- if (sc_r_power == six)
- {
- rpm2 = r2*r2; /* r4 */
- rp = rpm2*r2; /* r6 */
- }
- else if (sc_r_power == fourtyeight)
- {
- rp = r2*r2*r2; /* r6 */
- rp = rp*rp; /* r12 */
- rp = rp*rp; /* r24 */
- rp = rp*rp; /* r48 */
- rpm2 = rp/r2; /* r46 */
- }
- else
- {
- rp = pow(r2, half*sc_r_power); /* not currently supported as input, but can handle it */
- rpm2 = rp/r2;
- }
-
- /* Loop over state A(0) and B(1) */
- for (i = 0; i < 2; i++)
- {
- if ((c6[i] > 0) && (c12[i] > 0))
- {
- /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
- * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
- */
- sigma6[i] = half*c12[i]/c6[i];
- sigma2[i] = pow(half*c12[i]/c6[i], 1.0/3.0);
- /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on
- what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */
- if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */
- {
- sigma6[i] = sigma6_min;
- sigma2[i] = sigma2_min;
- }
- }
- else
- {
- sigma6[i] = sigma6_def;
- sigma2[i] = sigma2_def;
- }
- if (sc_r_power == six)
- {
- sigma_pow[i] = sigma6[i];
- sigma_powm2[i] = sigma6[i]/sigma2[i];
- }
- else if (sc_r_power == fourtyeight)
- {
- sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */
- sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
- sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
- sigma_powm2[i] = sigma_pow[i]/sigma2[i];
- }
- else
- { /* not really supported as input, but in here for testing the general case*/
- sigma_pow[i] = pow(sigma2[i], sc_r_power/2);
- sigma_powm2[i] = sigma_pow[i]/(sigma2[i]);
- }
- }
-
- /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
- if ((c12[0] > 0) && (c12[1] > 0))
- {
- alpha_vdw_eff = 0;
- alpha_coul_eff = 0;
- }
- else
- {
- alpha_vdw_eff = alpha_vdw;
- alpha_coul_eff = alpha_coul;
- }
-
- /* Loop over A and B states again */
- for (i = 0; i < 2; i++)
- {
- fscal_elec[i] = 0;
- fscal_vdw[i] = 0;
- velec[i] = 0;
- vvdw[i] = 0;
-
- /* Only spend time on A or B state if it is non-zero */
- if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
- {
- /* Coulomb */
- rpinv = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
- r_coul = pow(rpinv, -one/sc_r_power);
-
- /* Electrostatics table lookup data */
- rtab = r_coul*tabscale;
- ntab = rtab;
- eps = rtab-ntab;
- eps2 = eps*eps;
- ntab = 12*ntab;
- /* Electrostatics */
- Y = vftab[ntab];
- F = vftab[ntab+1];
- Geps = eps*vftab[ntab+2];
- Heps2 = eps2*vftab[ntab+3];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+two*Heps2;
- velec[i] = qq[i]*VV;
- fscal_elec[i] = -qq[i]*FF*r_coul*rpinv*tabscale;
-
- /* Vdw */
- rpinv = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
- r_vdw = pow(rpinv, -one/sc_r_power);
- /* Vdw table lookup data */
- rtab = r_vdw*tabscale;
- ntab = rtab;
- eps = rtab-ntab;
- eps2 = eps*eps;
- ntab = 12*ntab;
- /* Dispersion */
- Y = vftab[ntab+4];
- F = vftab[ntab+5];
- Geps = eps*vftab[ntab+6];
- Heps2 = eps2*vftab[ntab+7];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+two*Heps2;
- vvdw[i] = c6[i]*VV;
- fscal_vdw[i] = -c6[i]*FF;
-
- /* Repulsion */
- Y = vftab[ntab+8];
- F = vftab[ntab+9];
- Geps = eps*vftab[ntab+10];
- Heps2 = eps2*vftab[ntab+11];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+two*Heps2;
- vvdw[i] += c12[i]*VV;
- fscal_vdw[i] -= c12[i]*FF;
- fscal_vdw[i] *= r_vdw*rpinv*tabscale;
- }
- }
- /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
- /* Assemble A and B states */
- velecsum = 0;
- vvdwsum = 0;
- dvdl_coul = 0;
- dvdl_vdw = 0;
- fscal = 0;
- for (i = 0; i < 2; i++)
- {
- velecsum += LFC[i]*velec[i];
- vvdwsum += LFV[i]*vvdw[i];
-
- fscal += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2;
-
- dvdl_coul += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i];
- dvdl_vdw += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i];
- }
-
- dvdl[efptCOUL] += dvdl_coul;
- dvdl[efptVDW] += dvdl_vdw;
-
- *velectot = velecsum;
- *vvdwtot = vvdwsum;
-
- return fscal;
-}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff