jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2013-06-27 21:24:14 +00:00
parent 77866dd578
commit 34b028f059
2 changed files with 32 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/MOLECULE/pair_hbond_dreiding_lj.cpp
View File

@ -79,8 +79,8 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
{
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d;
double fi[3],fj[3],delr1[3],delr2[3];
double r2inv,r10inv;
double switch1,switch2;
@ -175,14 +175,20 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
powint(c,pm.ap-1)*s;
eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4);
force_switch=0.0;
if (rsq > pm.cut_innersq) {
switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
(pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) /
pm.denom_vdw;
switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
(rsq-pm.cut_innersq) / pm.denom_vdw;
force_kernel = force_kernel*switch1 + eng_lj*switch2;
eng_lj *= switch1;
force_kernel *= switch1;
force_angle *= switch1;
force_switch = eng_lj*switch2/rsq;
eng_lj *= switch1;
}
if (eflag) {
@ -193,6 +199,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
a = factor_hb*force_angle/s;
b = factor_hb*force_kernel;
d = factor_hb*force_switch;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
@ -205,12 +212,12 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
vz1 = a11*delr1[2] + a12*delr2[2];
vz2 = a22*delr2[2] + a12*delr1[2];
fi[0] = vx1 + b*delx;
fi[1] = vy1 + b*dely;
fi[2] = vz1 + b*delz;
fj[0] = vx2 - b*delx;
fj[1] = vy2 - b*dely;
fj[2] = vz2 - b*delz;
fi[0] = vx1 + b*delx + d*delx;
fi[1] = vy1 + b*dely + d*dely;
fi[2] = vz1 + b*delz + d*delz;
fj[0] = vx2 - b*delx - d*delx;
fj[1] = vy2 - b*dely - d*dely;
fj[2] = vz2 - b*delz - d*delz;
f[i][0] += fi[0];
f[i][1] += fi[1];

25
src/MOLECULE/pair_hbond_dreiding_morse.cpp
View File

@ -50,8 +50,8 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
{
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,evdwl,ehbond;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond;
double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double fi[3],fj[3],delr1[3],delr2[3];
double r,dr,dexp,eng_morse,switch1,switch2;
int *ilist,*jlist,*klist,*numneigh,**firstneigh;
@ -143,6 +143,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp);
force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap);
force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s;
force_switch = 0.0;
if (rsq > pm.cut_innersq) {
switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
@ -150,8 +151,11 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
pm.denom_vdw;
switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
(rsq-pm.cut_innersq) / pm.denom_vdw;
force_kernel = force_kernel*switch1 + eng_morse*switch2;
eng_morse *= switch1;
force_kernel *= switch1;
force_angle *= switch1;
force_switch = eng_morse*switch2/rsq;
eng_morse *= switch1;
}
if (eflag) {
@ -162,6 +166,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
a = factor_hb*force_angle/s;
b = factor_hb*force_kernel;
d = factor_hb*force_switch;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
@ -174,12 +179,12 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
vz1 = a11*delr1[2] + a12*delr2[2];
vz2 = a22*delr2[2] + a12*delr1[2];
fi[0] = vx1 + b*delx;
fi[1] = vy1 + b*dely;
fi[2] = vz1 + b*delz;
fj[0] = vx2 - b*delx;
fj[1] = vy2 - b*dely;
fj[2] = vz2 - b*delz;
fi[0] = vx1 + (b+d)*delx;
fi[1] = vy1 + (b+d)*dely;
fi[2] = vz1 + (b+d)*delz;
fj[0] = vx2 - (b+d)*delx;
fj[1] = vy2 - (b+d)*dely;
fj[2] = vz2 - (b+d)*delz;
f[i][0] += fi[0];
f[i][1] += fi[1];