X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxlib%2Fnonbonded%2Fnb_free_energy.c;h=b4a0200fd1f57be757749770ef1f175a32c784a3;hb=a1fbbb07241495656429dbb560e93769368609a4;hp=6b5d6d89dcae8b49ee779c473893c399e7fd8bce;hpb=19d3c2e5d0c401eb59010960d11a18b6ba2c54c6;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c index 6b5d6d89dc..b4a0200fd1 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c @@ -935,198 +935,3 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, #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; -}