* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
- * Copyright (c) 2011 Christoph Junghans, Sebastian Fritsch
- * Copyright (c) 2012, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
- * Copyright (c) 2012, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
+ * Copyright (c) 2011 Christoph Junghans, Sebastian Fritsch.
+ * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
*
* GROMACS is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
-#include <math.h>
-
-#include "types/simple.h"
-#include "vec.h"
-#include "typedefs.h"
#include "nb_generic_adress.h"
-#include "nrnb.h"
-#include "nonbonded.h"
-#include "nb_kernel.h"
+#include <math.h>
+
+#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
+#include "gromacs/legacyheaders/nonbonded.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/simple.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/fatalerror.h"
#define ALMOST_ZERO 1e-30
#define ALMOST_ONE 1-(1e-30)
void
gmx_nb_generic_adress_kernel(t_nblist * nlist,
- rvec * xx,
- rvec * ff,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb)
+ rvec * xx,
+ rvec * ff,
+ t_forcerec * fr,
+ t_mdatoms * mdatoms,
+ nb_kernel_data_t * kernel_data,
+ t_nrnb * nrnb)
{
- int nri,ntype,table_nelements,ielec,ivdw;
- real facel,gbtabscale;
- int n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid,nnn,n0;
- real shX,shY,shZ;
- real fscal,felec,fvdw,velec,vvdw,tx,ty,tz;
+ int nri, ntype, table_nelements, ielec, ivdw;
+ real facel, gbtabscale;
+ int n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid, nnn, n0;
+ real shX, shY, shZ;
+ real fscal, felec, fvdw, velec, vvdw, tx, ty, tz;
real rinvsq;
real iq;
- real qq,vctot;
- int nti,nvdwparam;
+ real qq, vctot;
+ int nti, nvdwparam;
int tj;
- real rt,r,eps,eps2,Y,F,Geps,Heps2,VV,FF,Fp,fijD,fijR;
+ real rt, r, eps, eps2, Y, F, Geps, Heps2, VV, FF, Fp, fijD, fijR;
real rinvsix;
real vvdwtot;
- real vvdw_rep,vvdw_disp;
- real ix,iy,iz,fix,fiy,fiz;
- real jx,jy,jz;
- real dx,dy,dz,rsq,rinv;
- real c6,c12,cexp1,cexp2,br;
+ real vvdw_rep, vvdw_disp;
+ real ix, iy, iz, fix, fiy, fiz;
+ real jx, jy, jz;
+ real dx, dy, dz, rsq, rinv;
+ real c6, c12, cexp1, cexp2, br;
real * charge;
real * shiftvec;
real * vdwparam;
real * x;
real * f;
int ewitab;
- real ewtabscale,eweps,sh_ewald,ewrt,ewtabhalfspace;
+ real ewtabscale, eweps, sh_ewald, ewrt, ewtabhalfspace;
real * ewtab;
- real rcoulomb2,rvdw,rvdw2,sh_invrc6;
- real rcutoff,rcutoff2;
- real rswitch_elec,rswitch_vdw,d,d2,sw,dsw,rinvcorr;
- real elec_swV3,elec_swV4,elec_swV5,elec_swF2,elec_swF3,elec_swF4;
- real vdw_swV3,vdw_swV4,vdw_swV5,vdw_swF2,vdw_swF3,vdw_swF4;
- gmx_bool bExactElecCutoff,bExactVdwCutoff,bExactCutoff;
-
- real * wf;
- real weight_cg1;
- real weight_cg2;
- real weight_product;
- real hybscal; /* the multiplicator to the force for hybrid interactions*/
- real force_cap;
- gmx_bool bCG;
- int egp_nr;
+ real rcoulomb2, rvdw, rvdw2, sh_dispersion, sh_repulsion;
+ real rcutoff, rcutoff2;
+ real rswitch_elec, rswitch_vdw, d, d2, sw, dsw, rinvcorr;
+ real elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
+ real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
+ gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoff;
+
+ real * wf;
+ real weight_cg1;
+ real weight_cg2;
+ real weight_product;
+ real hybscal; /* the multiplicator to the force for hybrid interactions*/
+ real force_cap;
+ gmx_bool bCG;
+ int egp_nr;
wf = mdatoms->wf;
rcoulomb2 = fr->rcoulomb*fr->rcoulomb;
rvdw = fr->rvdw;
rvdw2 = rvdw*rvdw;
- sh_invrc6 = fr->ic->sh_invrc6;
+ sh_dispersion = fr->ic->dispersion_shift.cpot;
+ sh_repulsion = fr->ic->repulsion_shift.cpot;
- if(fr->coulomb_modifier==eintmodPOTSWITCH)
+ if (fr->coulomb_modifier == eintmodPOTSWITCH)
{
d = fr->rcoulomb-fr->rcoulomb_switch;
elec_swV3 = -10.0/(d*d*d);
/* Avoid warnings from stupid compilers (looking at you, Clang!) */
elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
}
- if(fr->vdw_modifier==eintmodPOTSWITCH)
+ if (fr->vdw_modifier == eintmodPOTSWITCH)
{
d = fr->rvdw-fr->rvdw_switch;
vdw_swV3 = -10.0/(d*d*d);
vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
}
- bExactElecCutoff = (fr->coulomb_modifier!=eintmodNONE) || fr->eeltype == eelRF_ZERO;
- bExactVdwCutoff = (fr->vdw_modifier!=eintmodNONE);
+ bExactElecCutoff = (fr->coulomb_modifier != eintmodNONE) || fr->eeltype == eelRF_ZERO;
+ bExactVdwCutoff = (fr->vdw_modifier != eintmodNONE);
bExactCutoff = bExactElecCutoff || bExactVdwCutoff;
- if(bExactCutoff)
+ if (bExactCutoff)
{
rcutoff = ( fr->rcoulomb > fr->rvdw ) ? fr->rcoulomb : fr->rvdw;
rcutoff2 = rcutoff*rcutoff;
eps2 = 0.0;
/* 3 VdW parameters for buckingham, otherwise 2 */
- nvdwparam = (ivdw==GMX_NBKERNEL_VDW_BUCKINGHAM) ? 3 : 2;
+ nvdwparam = (ivdw == GMX_NBKERNEL_VDW_BUCKINGHAM) ? 3 : 2;
table_nelements = 12;
charge = mdatoms->chargeA;
vdwparam = fr->nbfp;
ntype = fr->ntype;
- for(n=0; (n<nlist->nri); n++)
+ for (n = 0; (n < nlist->nri); n++)
{
is3 = 3*nlist->shift[n];
- shX = shiftvec[is3];
+ shX = shiftvec[is3];
shY = shiftvec[is3+1];
shZ = shiftvec[is3+2];
- nj0 = nlist->jindex[n];
- nj1 = nlist->jindex[n+1];
- ii = nlist->iinr[n];
- ii3 = 3*ii;
+ nj0 = nlist->jindex[n];
+ nj1 = nlist->jindex[n+1];
+ ii = nlist->iinr[n];
+ ii3 = 3*ii;
ix = shX + x[ii3+0];
iy = shY + x[ii3+1];
iz = shZ + x[ii3+2];
iq = facel*charge[ii];
nti = nvdwparam*ntype*type[ii];
- vctot = 0;
+ vctot = 0;
vvdwtot = 0;
fix = 0;
fiy = 0;
- fiz = 0;
+ fiz = 0;
/* We need to find out if this i atom is part of an
all-atom or CG energy group */
- egp_nr=mdatoms->cENER[ii];
- bCG= !fr->adress_group_explicit[egp_nr];
-
+ egp_nr = mdatoms->cENER[ii];
+ bCG = !fr->adress_group_explicit[egp_nr];
+
weight_cg1 = wf[ii];
- if ((!bCG) && weight_cg1 < ALMOST_ZERO){
- continue;
- }
+ if ((!bCG) && weight_cg1 < ALMOST_ZERO)
+ {
+ continue;
+ }
- for(k=nj0; (k<nj1); k++)
+ for (k = nj0; (k < nj1); k++)
{
- jnr = nlist->jjnr[k];
+ jnr = nlist->jjnr[k];
weight_cg2 = wf[jnr];
weight_product = weight_cg1*weight_cg2;
if (weight_product < ALMOST_ZERO)
- {
- /* if it's a explicit loop, skip this atom */
+ {
+ /* if it's a explicit loop, skip this atom */
if (!bCG)
{
continue;
}
else /* if it's a coarse grained loop, include this atom */
{
- hybscal = 1.0;
+ hybscal = 1.0;
}
}
else if (weight_product >= ALMOST_ONE)
{
-
- /* if it's a explicit loop, include this atom */
- if(!bCG)
+
+ /* if it's a explicit loop, include this atom */
+ if (!bCG)
{
- hybscal = 1.0;
- }
+ hybscal = 1.0;
+ }
else /* if it's a coarse grained loop, skip this atom */
{
continue;
{
hybscal = weight_product;
- if(bCG)
+ if (bCG)
{
hybscal = 1.0 - hybscal;
}
}
- j3 = 3*jnr;
- jx = x[j3+0];
- jy = x[j3+1];
- jz = x[j3+2];
- dx = ix - jx;
- dy = iy - jy;
- dz = iz - jz;
+ j3 = 3*jnr;
+ jx = x[j3+0];
+ jy = x[j3+1];
+ jz = x[j3+2];
+ dx = ix - jx;
+ dy = iy - jy;
+ dz = iz - jz;
rsq = dx*dx+dy*dy+dz*dz;
rinv = gmx_invsqrt(rsq);
- rinvsq = rinv*rinv;
+ rinvsq = rinv*rinv;
felec = 0;
fvdw = 0;
velec = 0;
vvdw = 0;
- if(bExactCutoff && rsq>rcutoff2)
+ if (bExactCutoff && rsq > rcutoff2)
{
continue;
}
- if(ielec==GMX_NBKERNEL_ELEC_CUBICSPLINETABLE || ivdw==GMX_NBKERNEL_VDW_CUBICSPLINETABLE)
+ if (ielec == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE || ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE)
{
r = rsq*rinv;
rt = r*tabscale;
}
/* Coulomb interaction. ielec==0 means no interaction */
- if(ielec!=GMX_NBKERNEL_ELEC_NONE)
+ if (ielec != GMX_NBKERNEL_ELEC_NONE)
{
qq = iq*charge[jnr];
- switch(ielec)
+ switch (ielec)
{
case GMX_NBKERNEL_ELEC_NONE:
break;
case GMX_NBKERNEL_ELEC_GENERALIZEDBORN:
/* GB */
- gmx_fatal(FARGS,"Death & horror! GB generic interaction not implemented.\n");
+ gmx_fatal(FARGS, "Death & horror! GB generic interaction not implemented.\n");
break;
case GMX_NBKERNEL_ELEC_EWALD:
eweps = ewrt-ewitab;
ewitab = 4*ewitab;
felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
- rinvcorr = (fr->coulomb_modifier==eintmodPOTSHIFT) ? rinv-fr->ic->sh_ewald : rinv;
+ rinvcorr = (fr->coulomb_modifier == eintmodPOTSHIFT) ? rinv-fr->ic->sh_ewald : rinv;
velec = qq*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
felec = qq*rinv*(rinvsq-felec);
break;
default:
- gmx_fatal(FARGS,"Death & horror! No generic coulomb interaction for ielec=%d.\n",ielec);
+ gmx_fatal(FARGS, "Death & horror! No generic coulomb interaction for ielec=%d.\n", ielec);
break;
}
- if(fr->coulomb_modifier==eintmodPOTSWITCH)
+ if (fr->coulomb_modifier == eintmodPOTSWITCH)
{
d = rsq*rinv-fr->rcoulomb_switch;
- d = (d>0.0) ? d : 0.0;
+ d = (d > 0.0) ? d : 0.0;
d2 = d*d;
sw = 1.0+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
dsw = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
/* Once we have used velec to update felec we can modify velec too */
velec *= sw;
}
- if(bExactElecCutoff)
+ if (bExactElecCutoff)
{
- felec = (rsq<=rcoulomb2) ? felec : 0.0;
- velec = (rsq<=rcoulomb2) ? velec : 0.0;
+ felec = (rsq <= rcoulomb2) ? felec : 0.0;
+ velec = (rsq <= rcoulomb2) ? velec : 0.0;
}
vctot += velec;
} /* End of coulomb interactions */
/* VdW interaction. ivdw==0 means no interaction */
- if(ivdw!=GMX_NBKERNEL_VDW_NONE)
+ if (ivdw != GMX_NBKERNEL_VDW_NONE)
{
tj = nti+nvdwparam*type[jnr];
- switch(ivdw)
+ switch (ivdw)
{
case GMX_NBKERNEL_VDW_NONE:
break;
vvdw_disp = c6*rinvsix;
vvdw_rep = c12*rinvsix*rinvsix;
fvdw = (vvdw_rep-vvdw_disp)*rinvsq;
- if(fr->vdw_modifier==eintmodPOTSHIFT)
+ if (fr->vdw_modifier == eintmodPOTSHIFT)
{
- vvdw = (vvdw_rep-c12*sh_invrc6*sh_invrc6)*(1.0/12.0)-(vvdw_disp-c6*sh_invrc6)*(1.0/6.0);
+ vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion)/6.0;
}
else
{
br = cexp2*rsq*rinv;
vvdw_rep = cexp1*exp(-br);
fvdw = (br*vvdw_rep-vvdw_disp)*rinvsq;
- if(fr->vdw_modifier==eintmodPOTSHIFT)
+ if (fr->vdw_modifier == eintmodPOTSHIFT)
{
- vvdw = (vvdw_rep-cexp1*exp(-cexp2*rvdw))-(vvdw_disp-c6*sh_invrc6)/6.0;
+ vvdw = (vvdw_rep-cexp1*exp(-cexp2*rvdw)) - (vvdw_disp + c6*sh_dispersion)/6.0;
}
else
{
break;
default:
- gmx_fatal(FARGS,"Death & horror! No generic VdW interaction for ivdw=%d.\n",ivdw);
+ gmx_fatal(FARGS, "Death & horror! No generic VdW interaction for ivdw=%d.\n", ivdw);
break;
}
- if(fr->vdw_modifier==eintmodPOTSWITCH)
+ if (fr->vdw_modifier == eintmodPOTSWITCH)
{
d = rsq*rinv-fr->rvdw_switch;
- d = (d>0.0) ? d : 0.0;
+ d = (d > 0.0) ? d : 0.0;
d2 = d*d;
sw = 1.0+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
fvdw = fvdw*sw - rinv*vvdw*dsw;
vvdw *= sw;
}
- if(bExactVdwCutoff)
+ if (bExactVdwCutoff)
{
- fvdw = (rsq<=rvdw2) ? fvdw : 0.0;
- vvdw = (rsq<=rvdw2) ? vvdw : 0.0;
+ fvdw = (rsq <= rvdw2) ? fvdw : 0.0;
+ vvdw = (rsq <= rvdw2) ? vvdw : 0.0;
}
vvdwtot += vvdw;
} /* end VdW interactions */
fscal = felec+fvdw;
- if(!bCG && force_cap>0 && (fabs(fscal)> force_cap))
+ if (!bCG && force_cap > 0 && (fabs(fscal) > force_cap))
{
- fscal=force_cap*fscal/fabs(fscal);
+ fscal = force_cap*fscal/fabs(fscal);
}
fscal *= hybscal;
velecgrp[ggid] += vctot;
vvdwgrp[ggid] += vvdwtot;
}
- /* Estimate flops, average for generic kernel:
- * 12 flops per outer iteration
- * 50 flops per inner iteration
+ /* Estimate flops, average for generic adress kernel:
+ * 14 flops per outer iteration
+ * 54 flops per inner iteration
*/
- inc_nrnb(nrnb,eNR_NBKERNEL_GENERIC,nlist->nri*12 + nlist->jindex[n]*50);
+ inc_nrnb(nrnb, eNR_NBKERNEL_GENERIC_ADRESS, nlist->nri*14 + nlist->jindex[n]*54);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff