--- /home/sjplimp/oldlammps/doc/dihedral_style.txt 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/doc/dihedral_style.txt 2005-12-07 10:07:23.000000000 -0700 @@ -44,7 +44,7 @@ torsions that contain the j-k bond in an i-j-k-l torsion. :l Some force fields let {n} be positive or negative which corresponds to -{d} = 1 or -1. :ule,l +{d} = 1 or -1 for the harmonic style. :ule,l A style of {none} means dihedral forces are not computed, even if dihedrals are defined. --- /home/sjplimp/oldlammps/doc/dihedral_coeff.txt 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/doc/dihedral_coeff.txt 2005-12-07 10:07:23.000000000 -0700 @@ -66,12 +66,12 @@ For style {charmm}, specify 4 coefficients: K (energy) -n (1,2,3,4,6) -d (0 or 180 degrees) +n (integer >= 0) +d (integer value of degrees) weighting factor (0.0 to 1.0) :ul The weighting factor is applied to pairwise interaction between the -1st and 4th atoms in the dihedral. +1st and 4th atoms in the dihedral. :line @@ -151,7 +151,7 @@ K (energy) d (+1 or -1) -n (1,2,3,4,6) :ul +n (integer >= 0) :ul :line --- /home/sjplimp/oldlammps/doc/dihedral_style.html 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/doc/dihedral_style.html 2005-12-07 10:07:29.000000000 -0700 @@ -47,7 +47,7 @@ torsions that contain the j-k bond in an i-j-k-l torsion.
  • Some force fields let n be positive or negative which corresponds to -d = 1 or -1. +d = 1 or -1 for the harmonic style.

    A style of none means dihedral forces are not computed, even if dihedrals are defined. --- /home/sjplimp/oldlammps/doc/dihedral_coeff.html 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/doc/dihedral_coeff.html 2005-12-07 10:07:29.000000000 -0700 @@ -69,12 +69,12 @@

    For style charmm, specify 4 coefficients:

    The weighting factor is applied to pairwise interaction between the -1st and 4th atoms in the dihedral. +1st and 4th atoms in the dihedral.

    @@ -154,7 +154,7 @@

    --- /home/sjplimp/oldlammps/src/MOLECULE/dihedral_charmm.cpp 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/src/MOLECULE/dihedral_charmm.cpp 2005-12-07 10:10:35.000000000 -0700 @@ -46,7 +46,9 @@ memory->sfree(setflag); memory->sfree(k); memory->sfree(multiplicity); - memory->sfree(sign); + memory->sfree(shift); + memory->sfree(cos_shift); + memory->sfree(sin_shift); memory->sfree(weight); } } @@ -55,15 +57,14 @@ void DihedralCharmm::compute(int eflag, int vflag) { - int m,n,i1,i2,i3,i4,type,factor; + int i,m,n,i1,i2,i3,i4,type,factor; double rfactor; double vb1x,vb1y,vb1z,vb2x,vb2y; - double vb2z,vb2xm,vb2ym,vb2zm,vb3x,vb3y,vb3z,sb1; - double sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; - double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; - double c2mag,sc1,sc2,s1,s12,c,p,pd,rc2,a,a11,a22; - double a33,a12,a13,a23,sx1,sx2,sx12,sy1,sy2,sy12; - double sz1,sz2,sz12,s2; + double vb2z,vb2xm,vb2ym,vb2zm,vb3x,vb3y,vb3z; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb,dfx,dfy,dfz,dgx,dgy,dgz; + double dhx,dhy,dhz,dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; + double c,s,p,sx1,sx2,sx12,sy1,sy2,sy12,sz1,sz2,sz12,s2; int itype,jtype; double delx,dely,delz,rsq,r2inv,r6inv; double fforce,forcecoul,forcelj,phicoul,philj; @@ -125,49 +126,26 @@ vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; domain->minimum_image(&vb3x,&vb3y,&vb3z); + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + rginv = 1.0/rg; + ra2inv = 1.0/rasq; + rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); - // c0 calculation - - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); - - rb1 = sqrt(sb1); - rb3 = sqrt(sb3); - - c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; - - // 1st and 2nd angle - - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; - b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; - b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; - b3mag = sqrt(b3mag2); - - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); - c1mag = ctmp * r12c1; - - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); - c2mag = ctmp * r12c2; - - // cos and sin of 2 angles and final c - - sc1 = sqrt(1.0 - c1mag*c1mag); - if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; - - sc2 = sqrt(1.0 - c2mag*c2mag); - if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; - - s1 = sc1 * sc1; - s2 = sc2 * sc2; - s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (ax*bx + ay*by + az*bz)*rabinv; + s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); // error check @@ -191,68 +169,58 @@ if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - - // force & energy - // p = 1 + cos(n*phi) for d = 1 - // p = 1 - cos(n*phi) for d = -1 - // pd = dp/dc / 2 - + m = multiplicity[type]; - - if (m == 2) { - p = 2.0*c*c; - pd = 2.0*c; - } else if (m == 3) { - rc2 = c*c; - p = (4.0*rc2-3.0)*c + 1.0; - pd = 6.0*rc2 - 1.5; - } else if (m == 4) { - rc2 = c*c; - p = 8.0*(rc2-1)*rc2 + 2.0; - pd = (16.0*rc2-8.0)*c; - } else if (m == 6) { - rc2 = c*c; - p = ((32.0*rc2-48.0)*rc2 + 18.0)*rc2; - pd = (96.0*(rc2-1.0)*rc2 + 18.0)*c; - } else if (m == 1) { - p = c + 1.0; - pd = 0.5; - } else if (m == 5) { - rc2 = c*c; - p = ((16.0*rc2-20.0)*rc2 + 5.0)*c + 1.0; - pd = (40.0*rc2-30.0)*rc2 + 2.5; - } else if (m == 0) { - p = 2.0; - pd = 0.0; - } + p = 1.0; + df1 = 0.0; - if (sign[type] == 180) { - p = 2.0 - p; - pd = -pd; + for (i = 0; i < m; i++) { + ddf1 = p*c - df1*s; + df1 = p*s + df1*c; + p = ddf1; } - if (eflag) energy += rfactor * k[type] * p; - - a = 2.0 * k[type] * pd; - c = c * a; - s12 = s12 * a; - a11 = (-c*sb1*s1); - a22 = sb2*(2.0*c0*s12 - c*(s1+s2)); - a33 = (-c*sb3*s2); - a12 = r12c1*(c1mag*c*s1 + c2mag*s12); - a13 = rb1*rb3*s12; - a23 = r12c2*(-c2mag*c*s2 - c1mag*s12); - - sx1 = a11*vb1x + a12*vb2x + a13*vb3x; - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sx12 = a13*vb1x + a23*vb2x + a33*vb3x; - sy1 = a11*vb1y + a12*vb2y + a13*vb3y; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sy12 = a13*vb1y + a23*vb2y + a33*vb3y; - sz1 = a11*vb1z + a12*vb2z + a13*vb3z; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; - sz12 = a13*vb1z + a23*vb2z + a33*vb3z; + p = p*cos_shift[type] + df1*sin_shift[type]; + df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; + df1 *= -m; + p += 1.0; + + if (m == 0) { + p = 1.0 + cos_shift[type]; + df1 = 0.0; + } + if (eflag) energy += rfactor * k[type] * p; + + fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; + hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; + fga = fg*ra2inv*rginv; + hgb = hg*rb2inv*rginv; + gaa = -ra2inv*rg; + gbb = rb2inv*rg; + + dtfx = gaa*ax; + dtfy = gaa*ay; + dtfz = gaa*az; + dtgx = fga*ax - hgb*bx; + dtgy = fga*ay - hgb*by; + dtgz = fga*az - hgb*bz; + dthx = gbb*bx; + dthy = gbb*by; + dthz = gbb*bz; + + df = k[type] * df1; + + sx1 = df*dtfx; + sy1 = df*dtfy; + sz1 = df*dtfz; + sx2 = -df*dtgx; + sy2 = -df*dtgy; + sz2 = -df*dtgz; + sx12 = df*dthx; + sy12 = df*dthy; + sz12 = df*dthz; + // apply force to each of 4 atoms if (newton_bond || i1 < nlocal) { @@ -350,9 +318,14 @@ int n = atom->ndihedraltypes; k = (double *) memory->smalloc((n+1)*sizeof(double),"dihedral:k"); - multiplicity = (int *) memory->smalloc((n+1)*sizeof(double), - "dihedral:multiplicity"); - sign = (int *) memory->smalloc((n+1)*sizeof(double),"dihedral:sign"); + multiplicity = (int *) + memory->smalloc((n+1)*sizeof(double),"dihedral:multiplicity"); + shift = (int *) + memory->smalloc((n+1)*sizeof(double),"dihedral:shift"); + cos_shift = (double *) + memory->smalloc((n+1)*sizeof(double),"dihedral:cos_shift"); + sin_shift = (double *) + memory->smalloc((n+1)*sizeof(double),"dihedral:sin_shift"); weight = (double *) memory->smalloc((n+1)*sizeof(double),"dihedral:weight"); setflag = (int *) memory->smalloc((n+1)*sizeof(int),"dihedral:setflag"); @@ -371,23 +344,29 @@ int ilo,ihi; force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi); - + + // require integer values of shift for backwards compatibility + // arbitrary phase angle shift could be allowed, but would break + // backwards compatibility and is probably not needed + double k_one = atof(arg[1]); int multiplicity_one = atoi(arg[2]); - int sign_one = atoi(arg[3]); + int shift_one = atoi(arg[3]); double weight_one = atof(arg[4]); - if (multiplicity_one < 1 || multiplicity_one > 6) - error->all("Incorrect args for dihedral coefficients"); - if (sign_one != 0 && sign_one != 180) - error->all("Incorrect args for dihedral coefficients"); + if (multiplicity_one < 0) + error->all("Incorrect multiplicity arg for dihedral coefficients"); if (weight_one < 0.0 || weight_one > 1.0) - error->all("Incorrect args for dihedral coefficients"); + error->all("Incorrect weight arg for dihedral coefficients"); + double PI = 4.0*atan(1.0); + int count = 0; for (int i = ilo; i <= ihi; i++) { k[i] = k_one; - sign[i] = sign_one; + shift[i] = shift_one; + cos_shift[i] = cos(PI*shift_one/180.0); + sin_shift[i] = sin(PI*shift_one/180.0); multiplicity[i] = multiplicity_one; weight[i] = weight_one; setflag[i] = 1; @@ -440,7 +419,7 @@ { fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp); fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); - fwrite(&sign[1],sizeof(int),atom->ndihedraltypes,fp); + fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp); fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp); } @@ -455,13 +434,18 @@ if (comm->me == 0) { fread(&k[1],sizeof(double),atom->ndihedraltypes,fp); fread(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp); - fread(&sign[1],sizeof(int),atom->ndihedraltypes,fp); + fread(&shift[1],sizeof(int),atom->ndihedraltypes,fp); fread(&weight[1],sizeof(double),atom->ndihedraltypes,fp); } MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); - MPI_Bcast(&sign[1],atom->ndihedraltypes,MPI_INT,0,world); + MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); - for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; + double PI = 4.0*atan(1.0); + for (int i = 1; i <= atom->ndihedraltypes; i++) { + setflag[i] = 1; + cos_shift[i] = cos(PI*shift[i]/180.0); + sin_shift[i] = sin(PI*shift[i]/180.0); + } } --- /home/sjplimp/oldlammps/src/MOLECULE/dihedral_charmm.h 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/src/MOLECULE/dihedral_charmm.h 2005-12-07 10:10:35.000000000 -0700 @@ -28,8 +28,8 @@ void read_restart(FILE *); private: - double *k,*weight; - int *multiplicity,*sign; + double *k,*weight,*cos_shift,*sin_shift; + int *multiplicity,*shift; double **lj14_1,**lj14_2,**lj14_3,**lj14_4; void allocate(); --- /home/sjplimp/oldlammps/src/MOLECULE/dihedral_harmonic.cpp 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/src/MOLECULE/dihedral_harmonic.cpp 2005-12-07 10:10:35.000000000 -0700 @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + // harmonic dihedral potential #include "mpi.h" @@ -18,9 +22,9 @@ #include "string.h" #include "dihedral_harmonic.h" #include "atom.h" +#include "comm.h" #include "neighbor.h" #include "domain.h" -#include "comm.h" #include "force.h" #include "update.h" #include "memory.h" @@ -40,6 +44,8 @@ memory->sfree(k); memory->sfree(sign); memory->sfree(multiplicity); + memory->sfree(cos_shift); + memory->sfree(sin_shift); } } @@ -47,15 +53,14 @@ void DihedralHarmonic::compute(int eflag, int vflag) { - int m,n,i1,i2,i3,i4,type,factor; + int i,m,n,i1,i2,i3,i4,type,factor; double rfactor; double vb1x,vb1y,vb1z,vb2x,vb2y; - double vb2z,vb2xm,vb2ym,vb2zm,vb3x,vb3y,vb3z,sb1; - double sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; - double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; - double c2mag,sc1,sc2,s1,s12,c,p,pd,rc2,a,a11,a22; - double a33,a12,a13,a23,sx1,sx2,sx12,sy1,sy2,sy12; - double sz1,sz2,sz12,s2; + double vb2z,vb2xm,vb2ym,vb2zm,vb3x,vb3y,vb3z; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb,dfx,dfy,dfz,dgx,dgy,dgz; + double dhx,dhy,dhz,dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; + double c,s,p,sx1,sx2,sx12,sy1,sy2,sy12,sz1,sz2,sz12,s2; energy = 0.0; if (vflag) for (n = 0; n < 6; n++) virial[n] = 0.0; @@ -110,49 +115,26 @@ vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; domain->minimum_image(&vb3x,&vb3y,&vb3z); + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + rginv = 1.0/rg; + ra2inv = 1.0/rasq; + rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); - // c0 calculation - - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); - - rb1 = sqrt(sb1); - rb3 = sqrt(sb3); - - c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; - - // 1st and 2nd angle - - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; - b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; - b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; - b3mag = sqrt(b3mag2); - - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); - c1mag = ctmp * r12c1; - - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); - c2mag = ctmp * r12c2; - - // cos and sin of 2 angles and final c - - sc1 = sqrt(1.0 - c1mag*c1mag); - if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; - - sc2 = sqrt(1.0 - c2mag*c2mag); - if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; - - s1 = sc1 * sc1; - s2 = sc2 * sc2; - s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (ax*bx + ay*by + az*bz)*rabinv; + s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); // error check @@ -176,68 +158,58 @@ if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - - // force & energy - // p = 1 + cos(n*phi) for d = 1 - // p = 1 - cos(n*phi) for d = -1 - // pd = dp/dc / 2 - + m = multiplicity[type]; - - if (m == 2) { - p = 2.0*c*c; - pd = 2.0*c; - } else if (m == 3) { - rc2 = c*c; - p = (4.0*rc2-3.0)*c + 1.0; - pd = 6.0*rc2 - 1.5; - } else if (m == 4) { - rc2 = c*c; - p = 8.0*(rc2-1)*rc2 + 2.0; - pd = (16.0*rc2-8.0)*c; - } else if (m == 6) { - rc2 = c*c; - p = ((32.0*rc2-48.0)*rc2 + 18.0)*rc2; - pd = (96.0*(rc2-1.0)*rc2 + 18.0)*c; - } else if (m == 1) { - p = c + 1.0; - pd = 0.5; - } else if (m == 5) { - rc2 = c*c; - p = ((16.0*rc2-20.0)*rc2 + 5.0)*c + 1.0; - pd = (40.0*rc2-30.0)*rc2 + 2.5; - } else if (m == 0) { - p = 2.0; - pd = 0.0; - } + p = 1.0; + df1 = 0.0; - if (sign[type] == -1) { - p = 2.0 - p; - pd = -pd; + for (i = 0; i < m; i++) { + ddf1 = p*c - df1*s; + df1 = p*s + df1*c; + p = ddf1; } - if (eflag) energy += rfactor * k[type] * p; - - a = 2.0 * k[type] * pd; - c = c * a; - s12 = s12 * a; - a11 = (-c*sb1*s1); - a22 = sb2*(2.0*c0*s12 - c*(s1+s2)); - a33 = (-c*sb3*s2); - a12 = r12c1*(c1mag*c*s1 + c2mag*s12); - a13 = rb1*rb3*s12; - a23 = r12c2*(-c2mag*c*s2 - c1mag*s12); - - sx1 = a11*vb1x + a12*vb2x + a13*vb3x; - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sx12 = a13*vb1x + a23*vb2x + a33*vb3x; - sy1 = a11*vb1y + a12*vb2y + a13*vb3y; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sy12 = a13*vb1y + a23*vb2y + a33*vb3y; - sz1 = a11*vb1z + a12*vb2z + a13*vb3z; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; - sz12 = a13*vb1z + a23*vb2z + a33*vb3z; + p = p*cos_shift[type] + df1*sin_shift[type]; + df1 = df1*cos_shift[type] - ddf1*sin_shift[type]; + df1 *= -m; + p += 1.0; + + if (m == 0) { + p = 1.0 + cos_shift[type]; + df1 = 0.0; + } + if (eflag) energy += rfactor * k[type] * p; + + fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; + hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; + fga = fg*ra2inv*rginv; + hgb = hg*rb2inv*rginv; + gaa = -ra2inv*rg; + gbb = rb2inv*rg; + + dtfx = gaa*ax; + dtfy = gaa*ay; + dtfz = gaa*az; + dtgx = fga*ax - hgb*bx; + dtgy = fga*ay - hgb*by; + dtgz = fga*az - hgb*bz; + dthx = gbb*bx; + dthy = gbb*by; + dthz = gbb*bz; + + df = k[type] * df1; + + sx1 = df*dtfx; + sy1 = df*dtfy; + sz1 = df*dtfz; + sx2 = -df*dtgx; + sy2 = -df*dtgy; + sz2 = -df*dtgz; + sx12 = df*dthx; + sy12 = df*dthy; + sz12 = df*dthz; + // apply force to each of 4 atoms if (newton_bond || i1 < nlocal) { @@ -288,6 +260,10 @@ sign = (int *) memory->smalloc((n+1)*sizeof(double),"dihedral:sign"); multiplicity = (int *) memory->smalloc((n+1)*sizeof(double),"dihedral:multiplicity"); + cos_shift = (double *) + memory->smalloc((n+1)*sizeof(double),"dihedral:cos_shift"); + sin_shift = (double *) + memory->smalloc((n+1)*sizeof(double),"dihedral:sin_shift"); setflag = (int *) memory->smalloc((n+1)*sizeof(int),"dihedral:setflag"); for (int i = 1; i < n; i++) setflag[i] = 0; @@ -310,15 +286,26 @@ int sign_one = atoi(arg[2]); int multiplicity_one = atoi(arg[3]); + // require sign = +/- 1 for backwards compatibility + // arbitrary phase angle shift could be allowed, but would break + // backwards compatibility and is probably not needed + if (sign_one != -1 && sign_one != 1) - error->all("Incorrect args for dihedral coefficients"); - if (multiplicity_one < 0 || multiplicity_one > 6) - error->all("Incorrect args for dihedral coefficients"); - + error->all("Incorrect sign arg for dihedral coefficients"); + if (multiplicity_one < 0) + error->all("Incorrect multiplicity arg for dihedral coefficients"); + int count = 0; for (int i = ilo; i <= ihi; i++) { k[i] = k_one; sign[i] = sign_one; + if (sign[i] == 1) { + cos_shift[i] = 1; + sin_shift[i] = 0; + } else { + cos_shift[i] = -1; + sin_shift[i] = 0; + } multiplicity[i] = multiplicity_one; setflag[i] = 1; count++; @@ -354,6 +341,15 @@ MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world); MPI_Bcast(&sign[1],atom->ndihedraltypes,MPI_INT,0,world); MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world); - - for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1; + + for (int i = 1; i <= atom->ndihedraltypes; i++) { + setflag[i] = 1; + if (sign[i] == 1) { + cos_shift[i] = 1; + sin_shift[i] = 0; + } else { + cos_shift[i] = -1; + sin_shift[i] = 0; + } + } } --- /home/sjplimp/oldlammps/src/MOLECULE/dihedral_harmonic.h 2005-11-09 15:08:07.000000000 -0700 +++ /home/sjplimp/lammps/src/MOLECULE/dihedral_harmonic.h 2005-12-07 10:10:35.000000000 -0700 @@ -27,7 +27,7 @@ void read_restart(FILE *); private: - double *k; + double *k,*cos_shift,*sin_shift; int *sign,*multiplicity; void allocate();