diff -Naur lammps-17Jan08/src/ASPHERE/pair_gayberne.cpp lammps-18Jan08/src/ASPHERE/pair_gayberne.cpp --- lammps-17Jan08/src/ASPHERE/pair_gayberne.cpp 2007-12-06 17:11:14.000000000 -0700 +++ lammps-18Jan08/src/ASPHERE/pair_gayberne.cpp 2008-01-17 10:15:25.000000000 -0700 @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + cetain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -35,7 +35,7 @@ using namespace LAMMPS_NS; -enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_ELLIPSE}; +enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE}; /* ---------------------------------------------------------------------- */ @@ -105,11 +105,13 @@ i = ilist[ii]; itype = type[i]; - MathExtra::quat_to_mat_trans(quat[i],a1); - MathExtra::diag_times3(well[itype],a1,temp); - MathExtra::transpose_times3(a1,temp,b1); - MathExtra::diag_times3(shape[itype],a1,temp); - MathExtra::transpose_times3(a1,temp,g1); + if (form[itype][itype] == ELLIPSE_ELLIPSE) { + MathExtra::quat_to_mat_trans(quat[i],a1); + MathExtra::diag_times3(well[itype],a1,temp); + MathExtra::transpose_times3(a1,temp,b1); + MathExtra::diag_times3(shape[itype],a1,temp); + MathExtra::transpose_times3(a1,temp,g1); + } jlist = firstneigh[i]; jnum = numneigh[i]; @@ -151,6 +153,20 @@ rtor[0] = rtor[1] = rtor[2] = 0.0; break; + case SPHERE_ELLIPSE: + MathExtra::quat_to_mat_trans(quat[j],a2); + MathExtra::diag_times3(well[jtype],a2,temp); + MathExtra::transpose_times3(a2,temp,b2); + MathExtra::diag_times3(shape[jtype],a2,temp); + MathExtra::transpose_times3(a2,temp,g2); + one_eng = gayberne_lj(j,i,a2,b2,g2,r12,rsq,fforce,rtor); + break; + + case ELLIPSE_SPHERE: + one_eng = gayberne_lj(i,j,a1,b1,g1,r12,rsq,fforce,ttor); + rtor[0] = rtor[1] = rtor[2] = 0.0; + break; + default: MathExtra::quat_to_mat_trans(quat[j],a2); MathExtra::diag_times3(well[jtype],a2,temp); @@ -292,14 +308,14 @@ well[i][0] = pow(eia_one,-1.0/mu); well[i][1] = pow(eib_one,-1.0/mu); well[i][2] = pow(eic_one,-1.0/mu); - if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0) setwell[i] = 2; + if (eia_one == eib_one && eib_one == eic_one) setwell[i] = 2; else setwell[i] = 1; } if (eja_one != 0.0 || ejb_one != 0.0 || ejc_one != 0.0) { well[j][0] = pow(eja_one,-1.0/mu); well[j][1] = pow(ejb_one,-1.0/mu); well[j][2] = pow(ejc_one,-1.0/mu); - if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0) setwell[j] = 2; + if (eja_one == ejb_one && ejb_one == ejc_one) setwell[j] = 2; else setwell[j] = 1; } setflag[i][j] = 1; @@ -372,11 +388,17 @@ atom->shape[j][1] != atom->shape[j][2]) jshape = 1; if (setwell[j] == 1) jshape = 1; - if (ishape == 0 && jshape == 0) form[i][j] = SPHERE_SPHERE; - else if (ishape == 0 || jshape == 0) form[i][j] = SPHERE_ELLIPSE; - else form[i][j] = ELLIPSE_ELLIPSE; + if (ishape == 0 && jshape == 0) + form[i][i] = form[j][j] = form[i][j] = form[j][i] = SPHERE_SPHERE; + else if (ishape == 0) { + form[i][i] = SPHERE_SPHERE; form[j][j] = ELLIPSE_ELLIPSE; + form[i][j] = SPHERE_ELLIPSE; form[j][i] = ELLIPSE_SPHERE; + } else if (jshape == 0) { + form[j][j] = SPHERE_SPHERE; form[i][i] = ELLIPSE_ELLIPSE; + form[j][i] = SPHERE_ELLIPSE; form[i][j] = ELLIPSE_SPHERE; + } else + form[i][i] = form[j][j] = form[i][j] = form[j][i] = ELLIPSE_ELLIPSE; - form[j][i] = form[i][j]; epsilon[j][i] = epsilon[i][j]; sigma[j][i] = sigma[i][j]; cut[j][i] = cut[i][j]; @@ -661,6 +683,156 @@ } /* ---------------------------------------------------------------------- + compute analytic energy, force (fforce), and torque (ttor) + between ellipsoid and lj particle +------------------------------------------------------------------------- */ + +double PairGayBerne::gayberne_lj(const int i,const int j,double a1[3][3], + double b1[3][3],double g1[3][3], + double *r12,const double rsq,double *fforce, + double *ttor) +{ + double tempv[3], tempv2[3]; + double temp[3][3]; + double temp1,temp2,temp3; + + int *type = atom->type; + int newton_pair = force->newton_pair; + int nlocal = atom->nlocal; + + double r12hat[3]; + MathExtra::normalize3(r12,r12hat); + double r = sqrt(rsq); + + // compute distance of closest approach + + double g12[3][3]; + g12[0][0] = g1[0][0]+shape[type[j]][0]; + g12[1][1] = g1[1][1]+shape[type[j]][0]; + g12[2][2] = g1[2][2]+shape[type[j]][0]; + g12[0][1] = g1[0][1]; g12[1][0] = g1[1][0]; + g12[0][2] = g1[0][2]; g12[2][0] = g1[2][0]; + g12[1][2] = g1[1][2]; g12[2][1] = g1[2][1]; + double kappa[3]; + MathExtra::mldivide3(g12,r12,kappa,error); + + // tempv = G12^-1*r12hat + + tempv[0] = kappa[0]/r; + tempv[1] = kappa[1]/r; + tempv[2] = kappa[2]/r; + double sigma12 = MathExtra::dot3(r12hat,tempv); + sigma12 = pow(0.5*sigma12,-0.5); + double h12 = r-sigma12; + + // energy + // compute u_r + + double varrho = sigma[type[i]][type[j]]/(h12+gamma*sigma[type[i]][type[j]]); + double varrho6 = pow(varrho,6.0); + double varrho12 = varrho6*varrho6; + double u_r = 4.0*epsilon[type[i]][type[j]]*(varrho12-varrho6); + + // compute eta_12 + + double eta = 2.0*lshape[type[i]]*lshape[type[j]]; + double det_g12 = MathExtra::det3(g12); + eta = pow(eta/det_g12,upsilon); + + // compute chi_12 + + double b12[3][3]; + double iota[3]; + b12[0][0] = b1[0][0] + well[type[j]][0]; + b12[1][1] = b1[1][1] + well[type[j]][0]; + b12[2][2] = b1[2][2] + well[type[j]][0]; + b12[0][1] = b1[0][1]; b12[1][0] = b1[1][0]; + b12[0][2] = b1[0][2]; b12[2][0] = b1[2][0]; + b12[1][2] = b1[1][2]; b12[2][1] = b1[2][1]; + MathExtra::mldivide3(b12,r12,iota,error); + + // tempv = G12^-1*r12hat + + tempv[0] = iota[0]/r; + tempv[1] = iota[1]/r; + tempv[2] = iota[2]/r; + double chi = MathExtra::dot3(r12hat,tempv); + chi = pow(chi*2.0,mu); + + // force + // compute dUr/dr + + temp1 = (2.0*varrho12*varrho-varrho6*varrho)/sigma[type[i]][type[j]]; + temp1 = temp1*24.0*epsilon[type[i]][type[j]]; + double u_slj = temp1*pow(sigma12,3.0)/2.0; + double dUr[3]; + temp2 = MathExtra::dot3(kappa,r12hat); + double uslj_rsq = u_slj/rsq; + dUr[0] = temp1*r12hat[0]+uslj_rsq*(kappa[0]-temp2*r12hat[0]); + dUr[1] = temp1*r12hat[1]+uslj_rsq*(kappa[1]-temp2*r12hat[1]); + dUr[2] = temp1*r12hat[2]+uslj_rsq*(kappa[2]-temp2*r12hat[2]); + + // compute dChi_12/dr + + double dchi[3]; + temp1 = MathExtra::dot3(iota,r12hat); + temp2 = -4.0/rsq*mu*pow(chi,(mu-1.0)/mu); + dchi[0] = temp2*(iota[0]-temp1*r12hat[0]); + dchi[1] = temp2*(iota[1]-temp1*r12hat[1]); + dchi[2] = temp2*(iota[2]-temp1*r12hat[2]); + + temp1 = -eta*u_r; + temp2 = eta*chi; + fforce[0] = temp1*dchi[0]-temp2*dUr[0]; + fforce[1] = temp1*dchi[1]-temp2*dUr[1]; + fforce[2] = temp1*dchi[2]-temp2*dUr[2]; + + // torque for particle 1 and 2 + // compute dUr + + tempv[0] = -uslj_rsq*kappa[0]; + tempv[1] = -uslj_rsq*kappa[1]; + tempv[2] = -uslj_rsq*kappa[2]; + MathExtra::row_times3(kappa,g1,tempv2); + MathExtra::cross3(tempv,tempv2,dUr); + + // compute d_chi + + MathExtra::row_times3(iota,b1,tempv); + MathExtra::cross3(tempv,iota,dchi); + temp1 = -4.0/rsq; + dchi[0] *= temp1; + dchi[1] *= temp1; + dchi[2] *= temp1; + + // compute d_eta + + double deta[3]; + deta[0] = deta[1] = deta[2] = 0.0; + compute_eta_torque(g12,a1,shape[type[i]],temp); + temp1 = -eta*upsilon; + for (int m = 0; m < 3; m++) { + for (int y = 0; y < 3; y++) tempv[y] = temp1*temp[m][y]; + MathExtra::cross3(a1[m],tempv,tempv2); + deta[0] += tempv2[0]; + deta[1] += tempv2[1]; + deta[2] += tempv2[2]; + } + + // torque + + temp1 = u_r*eta; + temp2 = u_r*chi; + temp3 = chi*eta; + + ttor[0] = (temp1*dchi[0]+temp2*deta[0]+temp3*dUr[0]) * -1.0; + ttor[1] = (temp1*dchi[1]+temp2*deta[1]+temp3*dUr[1]) * -1.0; + ttor[2] = (temp1*dchi[2]+temp2*deta[2]+temp3*dUr[2]) * -1.0; + + return temp1*chi; +} + +/* ---------------------------------------------------------------------- torque contribution from eta computes trace in the last doc equation for the torque derivative code comes from symbolic solver dump diff -Naur lammps-17Jan08/src/ASPHERE/pair_gayberne.h lammps-18Jan08/src/ASPHERE/pair_gayberne.h --- lammps-17Jan08/src/ASPHERE/pair_gayberne.h 2007-10-03 10:13:59.000000000 -0600 +++ lammps-18Jan08/src/ASPHERE/pair_gayberne.h 2008-01-17 10:15:25.000000000 -0700 @@ -53,6 +53,9 @@ double g1[3][3], double g2[3][3], double *r12, const double rsq, double *fforce, double *ttor, double *rtor); + double gayberne_lj(const int i, const int j, double a1[3][3], + double b1[3][3],double g1[3][3],double *r12, + const double rsq, double *fforce, double *ttor); void compute_eta_torque(double m[3][3], double m2[3][3], double *s, double ans[3][3]); }; diff -Naur lammps-17Jan08/tools/pymol_asphere/src/cartesian.cpp lammps-18Jan08/tools/pymol_asphere/src/cartesian.cpp --- lammps-17Jan08/tools/pymol_asphere/src/cartesian.cpp 2007-06-21 16:53:58.000000000 -0600 +++ lammps-18Jan08/tools/pymol_asphere/src/cartesian.cpp 2008-01-17 10:11:02.000000000 -0700 @@ -170,12 +170,6 @@ } template -numtyp TwoD::dist(const TwoD &two) const { - TwoD diff=*this-two; - return numtyp(sqrt(double(diff.dot(diff)))); -} - -template numtyp TwoD::dist2(const TwoD &two) const { TwoD diff=*this-two; return diff.dot(diff); @@ -294,26 +288,6 @@ } template -numtyp & ThreeD::operator[](unsigned i) { - switch(i) { - case X: return x; - case Y: return y; - case Z: return z; - } - return x; -} - -template -numtyp ThreeD::operator[](unsigned i) const { - switch(i) { - case X: return x; - case Y: return y; - case Z: return z; - } - return x; -} - -template void ThreeD::operator = (const ThreeD &two) { #ifdef NANCHECK //assert(a::not_nan(two.x) && a::not_nan(two.y) && a::not_nan(two.z)); @@ -396,11 +370,6 @@ return (x*two.x+y*two.y+z*two.z); } -template -void ThreeD::operator *= (const numtyp &two) { - x*=two; y*=two; z*=two; -} - // Move coordinates into array template void ThreeD::to_array(numtyp *array) { @@ -500,26 +469,11 @@ } template -numtyp ThreeD::dist(const ThreeD &two) { - return (*this-two).norm(); -} - -template numtyp ThreeD::dist2(const ThreeD &two) { ThreeD diff=*this-two; return diff.dot(diff); } -// For normalizing a vector -template -void ThreeD::normalize() { - numtyp temp=norm(); - #ifdef NANCHECK - assert(temp!=0); - #endif - *this/=temp; -} - // Return unit vector template ThreeD ThreeD::unit() { diff -Naur lammps-17Jan08/tools/pymol_asphere/src/cartesian.h lammps-18Jan08/tools/pymol_asphere/src/cartesian.h --- lammps-17Jan08/tools/pymol_asphere/src/cartesian.h 2007-06-21 16:53:58.000000000 -0600 +++ lammps-18Jan08/tools/pymol_asphere/src/cartesian.h 2008-01-17 10:11:02.000000000 -0700 @@ -204,8 +204,23 @@ friend ThreeD operator* <>(const numtyp, const ThreeD &two); friend ThreeD operator/ <>(const numtyp, const ThreeD &two); - numtyp &operator[](unsigned i); - numtyp operator[](unsigned i) const; + inline numtyp &operator[](unsigned i) { + switch(i) { + case X: return x; + case Y: return y; + case Z: return z; + } + return x; + } + + inline numtyp operator[](unsigned i) const { + switch(i) { + case X: return x; + case Y: return y; + case Z: return z; + } + return x; + } bool operator == (const ThreeD &two) const; bool operator != (const ThreeD &two) const; @@ -222,7 +237,10 @@ void operator += (const ThreeD &two); void operator -= (const numtyp &two); void operator -= (const ThreeD &two); - void operator *= (const numtyp &two); + inline void operator *= (const numtyp &two) { + x*=two; y*=two; z*=two; + } + void operator /= (const numtyp &two); /// Move coordinates into array @@ -253,11 +271,20 @@ /// Magnitude of vector numtyp hypot() const; /// Distance between two points - numtyp dist(const ThreeD &two); + inline numtyp dist(const ThreeD &two) { + return (*this-two).norm(); + } /// Distance squared between two points numtyp dist2(const ThreeD &two); /// Converts \b *this to the unit vector - void normalize(); + inline void normalize() { + numtyp temp=norm(); + #ifdef NANCHECK + assert(temp!=0); + #endif + *this/=temp; + } + /// Return the unit vector of \b *this ThreeD unit(); diff -Naur lammps-17Jan08/tools/pymol_asphere/src/glsurface.cpp lammps-18Jan08/tools/pymol_asphere/src/glsurface.cpp --- lammps-17Jan08/tools/pymol_asphere/src/glsurface.cpp 2007-10-22 15:36:28.000000000 -0600 +++ lammps-18Jan08/tools/pymol_asphere/src/glsurface.cpp 2008-01-17 10:11:02.000000000 -0700 @@ -333,7 +333,8 @@ writepymolheader(out); for (i=0; i