ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8206 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2012-06-01 17:20:38 +00:00
parent 9377bfd64b
commit 978a99e098
16 changed files with 37 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/MOLECULE/dihedral_charmm.cpp
View File

@ -170,7 +170,7 @@ void DihedralCharmm::compute(int eflag, int vflag)
m = multiplicity[type];
p = 1.0;
df1 = 0.0;
ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s;

2
src/MOLECULE/dihedral_harmonic.cpp
View File

@ -158,7 +158,7 @@ void DihedralHarmonic::compute(int eflag, int vflag)
m = multiplicity[type];
p = 1.0;
df1 = 0.0;
ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s;

5
src/USER-CG-CMM/angle_sdk.cpp
View File

@ -121,6 +121,8 @@ void AngleSDK::compute(int eflag, int vflag)
// so this has to be done here and not in the
// general non-bonded code.
f13 = e13 = delx3 = dely3 = delz3 = 0.0;
if (repflag) {
delx3 = x[i1][0] - x[i3][0];
@ -208,7 +210,8 @@ void AngleSDK::compute(int eflag, int vflag)
if (evflag) {
ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
ev_tally13(i1,i3,nlocal,newton_bond,e13,f13,delx3,dely3,delz3);
if (repflag)
ev_tally13(i1,i3,nlocal,newton_bond,e13,f13,delx3,dely3,delz3);
}
}
}

4
src/USER-MISC/angle_cosine_shift.cpp
View File

@ -104,10 +104,10 @@ void AngleCosineShift::compute(int eflag, int vflag)
if (s < SMALL) s = SMALL;
// force & energy
if (eflag) eangle = -k[type]-kcos*c-ksin*s;
kcos=kcost[type];
ksin=ksint[type];
if (eflag) eangle = -k[type]-kcos*c-ksin*s;
cps = c/s; // NOTE absorbed one c
a11 = (-kcos +ksin*cps )*c/ rsq1;

9
src/USER-OMP/angle_cosine_shift_omp.cpp
View File

@ -85,13 +85,14 @@ void AngleCosineShiftOMP::eval(int nfrom, int nto, ThrData * const thr)
{
int i1,i2,i3,n,type;
double delx1,dely1,delz1,delx2,dely2,delz2;
double eangle,f1[3],f3[3];
double rsq1,rsq2,r1,r2,c,s,cps,kcos,ksin,a11,a12,a22;
double f1[3],f3[3];
double rsq1,rsq2,r1,r2,c,s,cps,a11,a12,a22;
const double * const * const x = atom->x;
double * const * const f = thr->get_f();
const int * const * const anglelist = neighbor->anglelist;
const int nlocal = atom->nlocal;
double eangle = 0.0;
for (n = nfrom; n < nto; n++) {
i1 = anglelist[n][0];
@ -130,10 +131,10 @@ void AngleCosineShiftOMP::eval(int nfrom, int nto, ThrData * const thr)
if (s < SMALL) s = SMALL;
// force & energy
const double kcos=kcost[type];
const double ksin=ksint[type];
if (EFLAG) eangle = -k[type]-kcos*c-ksin*s;
kcos=kcost[type];
ksin=ksint[type];
cps = c/s; // NOTE absorbed one c
a11 = (-kcos +ksin*cps )*c/ rsq1;

2
src/USER-OMP/angle_dipole_omp.cpp
View File

@ -83,7 +83,7 @@ void AngleDipoleOMP::eval(int nfrom, int nto, ThrData * const thr)
int iRef,iDip,iDummy,n,type;
double delx,dely,delz;
double eangle,tangle;
double r,dr,cosGamma,deltaGamma,kdg,rmu;
double r,cosGamma,deltaGamma,kdg,rmu;
const double * const * const x = atom->x; // position vector
const double * const * const mu = atom->mu; // point-dipole components and moment magnitude

16
src/USER-OMP/angle_sdk_omp.cpp
View File

@ -140,6 +140,8 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr)
// so this has to be done here and not in the
// general non-bonded code.
f13 = e13 = delx3 = dely3 = delz3 = 0.0;
if (repflag) {
delx3 = x[i1][0] - x[i3][0];
@ -151,9 +153,6 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr)
const int type1 = atom->type[i1];
const int type3 = atom->type[i3];
f13=0.0;
e13=0.0;
if (rsq3 < rminsq[type1][type3]) {
const int ljt = lj_type[type1][type3];
const double r2inv = 1.0/rsq3;
@ -224,10 +223,11 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr)
f[i3][2] += f3[2] - f13*delz3;
}
if (EVFLAG) ev_tally_thr(this,i1,i2,i3,nlocal,NEWTON_BOND,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2,thr);
if (EVFLAG) ev_tally13_thr(this,i1,i3,nlocal,NEWTON_BOND,
e13,f13,delx3,dely3,delz3,thr);
if (EVFLAG) {
ev_tally_thr(this,i1,i2,i3,nlocal,NEWTON_BOND,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2,thr);
if (repflag) ev_tally13_thr(this,i1,i3,nlocal,NEWTON_BOND,
e13,f13,delx3,dely3,delz3,thr);
}
}
}

2
src/USER-OMP/dihedral_charmm_omp.cpp
View File

@ -195,7 +195,7 @@ void DihedralCharmmOMP::eval(int nfrom, int nto, ThrData * const thr)
m = multiplicity[type];
p = 1.0;
df1 = 0.0;
ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s;

2
src/USER-OMP/dihedral_harmonic_omp.cpp
View File

@ -181,7 +181,7 @@ void DihedralHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr)
m = multiplicity[type];
p = 1.0;
df1 = 0.0;
ddf1 = df1 = 0.0;
for (i = 0; i < m; i++) {
ddf1 = p*c - df1*s;

7
src/USER-OMP/dihedral_table_omp.cpp
View File

@ -82,8 +82,7 @@ void DihedralTableOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG, int NEWTON_BOND>
void DihedralTableOMP::eval(int nfrom, int nto, ThrData * const thr)
{
int i1,i2,i3,i4,i,m,n,type;
int i1,i2,i3,i4,n,type;
double edihedral,f1[3],f2[3],f3[3],f4[3];
const double * const * const x = atom->x;
@ -133,8 +132,8 @@ void DihedralTableOMP::eval(int nfrom, int nto, ThrData * const thr)
// n123 & n234: These two unit vectors are normal to the planes
// defined by atoms 1,2,3 and 2,3,4.
double n123[g_dim]; //n123=vb12 x vb23 / |vb12 x vb23| ("x" is cross product)
double n234[g_dim]; //n234=vb34 x vb23 / |vb34 x vb23| ("x" is cross product)
double n123[g_dim]; //n123=vb23 x vb12 / |vb23 x vb12| ("x" is cross product)
double n234[g_dim]; //n234=vb23 x vb34 / |vb23 x vb34| ("x" is cross product)
double proj12on23[g_dim];
// proj12on23[d] = (vb23[d]/|vb23|) * DotProduct(vb12,vb23)/|vb12|*|vb23|

10
src/USER-OMP/fix_omp.cpp
View File

@ -197,9 +197,9 @@ void FixOMP::init()
int check_hybrid;
last_pair_hybrid = NULL;
last_omp_style = NULL;
char *last_omp_name = NULL;
char *last_hybrid_name = NULL;
char *last_force_name = NULL;
const char *last_omp_name = NULL;
const char *last_hybrid_name = NULL;
const char *last_force_name = NULL;
// determine which is the last force style with OpenMP
// support as this is the one that has to reduce the forces
@ -211,7 +211,7 @@ void FixOMP::init()
(strcmp(force->name ## _style,"hybrid/overlay") == 0) ) \
check_hybrid=1; \
if (force->name->suffix_flag & Suffix::OMP) { \
last_force_name = #name; \
last_force_name = (const char *) #name; \
last_omp_name = force->name ## _style; \
last_omp_style = (void *) force->name; \
} \
@ -222,7 +222,7 @@ void FixOMP::init()
Class ## Hybrid *style = (Class ## Hybrid *) force->name; \
for (int i=0; i < style->nstyles; i++) { \
if (style->styles[i]->suffix_flag & Suffix::OMP) { \
last_force_name = #name; \
last_force_name = (const char *) #name; \
last_omp_name = style->keywords[i]; \
last_omp_style = style->styles[i]; \
} \

4
src/USER-OMP/pair_comb_omp.cpp
View File

@ -567,7 +567,7 @@ void PairCombOMP::Short_neigh_thr()
#pragma omp parallel default(none)
#endif
{
int nj,npntj,*neighptrj,itype,jtype;
int nj,npntj,*neighptrj;
int iparam_ij,*ilist,*jlist,*numneigh,**firstneigh;
int jnum,i,j,ii,jj;
double xtmp,ytmp,ztmp,rsq,delrj[3];
@ -597,7 +597,6 @@ void PairCombOMP::Short_neigh_thr()
if (iifrom < inum) {
for (ii = iifrom; ii < iito; ii++) {
i = ilist[ii];
itype = type[i];
#if defined(_OPENMP)
#pragma omp critical
@ -621,7 +620,6 @@ void PairCombOMP::Short_neigh_thr()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delrj[0] = xtmp - x[j][0];
delrj[1] = ytmp - x[j][1];

4
src/USER-OMP/pair_coul_wolf_omp.cpp
View File

@ -78,7 +78,7 @@ void PairCoulWolfOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
int i,j,ii,jj,jnum;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double rsq,forcecoul,factor_coul;
double prefactor;
@ -117,7 +117,6 @@ void PairCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp=fytmp=fztmp=0.0;
@ -135,7 +134,6 @@ void PairCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cut_coulsq) {
r = sqrt(rsq);

4
src/USER-OMP/pair_dpd_tstat_omp.cpp
View File

@ -106,13 +106,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairDPDTstatOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,randnum,factor_dpd;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
const double * const * const x = atom->x;
const double * const * const v = atom->v;
double * const * const f = thr->get_f();

3
src/USER-OMP/pair_gauss_omp.cpp
View File

@ -80,7 +80,7 @@ double PairGaussOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double r,rsq,r2inv,forcelj,factor_lj;
double rsq,r2inv,forcelj,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
int occ = 0;
@ -129,7 +129,6 @@ double PairGaussOMP::eval(int iifrom, int iito, ThrData * const thr)
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
forcelj = - 2.0*a[itype][jtype]*b[itype][jtype] * rsq *
exp(-b[itype][jtype]*rsq);
fpair = factor_lj*forcelj*r2inv;

3
src/USER-OMP/pair_yukawa_colloid_omp.cpp
View File

@ -76,7 +76,7 @@ void PairYukawaColloidOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,radi,radj;
double rsq,r,rinv,r2inv,screening,forceyukawa,factor;
double rsq,r,rinv,screening,forceyukawa,factor;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
@ -120,7 +120,6 @@ void PairYukawaColloidOMP::eval(int iifrom, int iito, ThrData * const thr)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
screening = exp(-kappa*(r-(radi+radj)));