diff --git a/src/USER-MISC/angle_fourier.cpp b/src/USER-MISC/angle_fourier.cpp index 3f62b059ea..beb3cc023b 100644 --- a/src/USER-MISC/angle_fourier.cpp +++ b/src/USER-MISC/angle_fourier.cpp @@ -200,7 +200,7 @@ double AngleFourier::equilibrium_angle(int i) double ret=MY_PI; if ( C2[i] != 0.0 ) { ret = (C1[i]/4.0/C2[i]); - if ( abs(ret) <= 1.0 ) ret = acos(-ret); + if ( fabs(ret) <= 1.0 ) ret = acos(-ret); } return ret; } diff --git a/src/USER-MISC/angle_fourier_simple.cpp b/src/USER-MISC/angle_fourier_simple.cpp index 63a972436c..2513050ae4 100644 --- a/src/USER-MISC/angle_fourier_simple.cpp +++ b/src/USER-MISC/angle_fourier_simple.cpp @@ -113,8 +113,7 @@ void AngleFourierSimple::compute(int eflag, int vflag) // handle sin(n th)/sin(th) singulatiries - if ( abs(c)-1.0 > 0.0001 ) - { + if ( fabs(c)-1.0 > 0.0001 ) { a = k[type]*C[type]*N[type]*sin(nth)/sin(th); } else { if ( c >= 0.0 ) { diff --git a/src/USER-MISC/dihedral_nharmonic.cpp b/src/USER-MISC/dihedral_nharmonic.cpp index c6cd4f7664..b5512de128 100644 --- a/src/USER-MISC/dihedral_nharmonic.cpp +++ b/src/USER-MISC/dihedral_nharmonic.cpp @@ -60,7 +60,7 @@ void DihedralNHarmonic::compute(int eflag, int vflag) double edihedral,f1[3],f2[3],f3[3],f4[3]; double sb1,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,a,a11,a22; + double c2mag,sc1,sc2,s1,s12,c,p,pd,a11,a22; double a33,a12,a13,a23,sx2,sy2,sz2; double s2,sin2; @@ -177,20 +177,19 @@ void DihedralNHarmonic::compute(int eflag, int vflag) // force & energy // p = sum (i=1,n) a_i * c**(i-1) // pd = dp/dc - p = this->a[type][0]; - pd = this->a[type][1]; + p = a[type][0]; + pd = a[type][1]; for (int i = 1; i < nterms[type]-1; i++) { - p += c * this->a[type][i]; - pd += c * static_cast(i+1) * this->a[type][i+1]; + p += c * a[type][i]; + pd += c * static_cast(i+1) * a[type][i+1]; c *= c; } - p += c * this->a[type][nterms[type]-1]; + p += c * a[type][nterms[type]-1]; if (eflag) edihedral = p; - a = pd; - c = c * a; - s12 = s12 * a; + c = c * pd; + s12 = s12 * pd; a11 = c*sb1*s1; a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); a33 = c*sb3*s2; diff --git a/src/USER-MISC/dihedral_quadratic.cpp b/src/USER-MISC/dihedral_quadratic.cpp index e617be70bd..59a5dc1991 100644 --- a/src/USER-MISC/dihedral_quadratic.cpp +++ b/src/USER-MISC/dihedral_quadratic.cpp @@ -54,14 +54,14 @@ DihedralQuadratic::~DihedralQuadratic() void DihedralQuadratic::compute(int eflag, int vflag) { - int i1,i2,i3,i4,i,m,n,type; + int i1,i2,i3,i4,n,type; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double edihedral,f1[3],f2[3],f3[3],f4[3]; double sb1,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,a,a11,a22; double a33,a12,a13,a23,sx2,sy2,sz2; - double s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + double s2,cx,cy,cz,cmag,dx,phi,si,sin2; edihedral = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -86,26 +86,22 @@ void DihedralQuadratic::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation @@ -179,7 +175,7 @@ void DihedralQuadratic::compute(int eflag, int vflag) me,x[i4][0],x[i4][1],x[i4][2]); } } - + if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; diff --git a/src/USER-MISC/dihedral_table.cpp b/src/USER-MISC/dihedral_table.cpp index 7465779c0d..493efcd4a3 100644 --- a/src/USER-MISC/dihedral_table.cpp +++ b/src/USER-MISC/dihedral_table.cpp @@ -23,10 +23,7 @@ #include #include #include -#define NDEBUG //(disable assert()) -#include #include -using namespace std; #include "atom.h" #include "comm.h" @@ -38,17 +35,371 @@ using namespace std; #include "error.h" #include "dihedral_table.h" -// Additional header files available for numerical debugging: -//#undef NDEBUG -//#define DIH_DEBUG_NUM -//#ifdef DIH_DEBUG_NUM -//#include "dihedral_table_dbg.h" //in "supporting_files/debug_numerical/" -//#endif -// and pointless posturing: -//#include "dihedral_table_nd_mod.h" //in "supporting_files/nd/" +#include "math_const.h" +#include "math_extra.h" +using namespace std; using namespace LAMMPS_NS; -using namespace DIHEDRAL_TABLE_NS; +using namespace MathConst; +using namespace MathExtra; + +// ------------------------------------------------------------------------ +// The following auxiliary functions were left out of the +// DihedralTable class either because they require template parameters, +// or because they have nothing to do with dihedral angles. +// ------------------------------------------------------------------------ + +// ------------------------------------------------------------------- +// --------- The function was stolen verbatim from the --------- +// --------- GNU Scientific Library (GSL, version 1.15) --------- +// ------------------------------------------------------------------- + +/* Author: Gerard Jungman */ +/* for description of method see [Engeln-Mullges + Uhlig, p. 96] + * + * diag[0] offdiag[0] 0 ..... offdiag[N-1] + * offdiag[0] diag[1] offdiag[1] ..... + * 0 offdiag[1] diag[2] + * 0 0 offdiag[2] ..... + * ... ... + * offdiag[N-1] ... + * + */ +// -- (A non-symmetric version of this function is also available.) -- + +enum { //GSL status return codes. + GSL_FAILURE = -1, + GSL_SUCCESS = 0, + GSL_ENOMEM = 8, + GSL_EZERODIV = 12, + GSL_EBADLEN = 19 +}; + + +static int solve_cyc_tridiag( const double diag[], size_t d_stride, + const double offdiag[], size_t o_stride, + const double b[], size_t b_stride, + double x[], size_t x_stride, + size_t N) +{ + int status = GSL_SUCCESS; + double * delta = (double *) malloc (N * sizeof (double)); + double * gamma = (double *) malloc (N * sizeof (double)); + double * alpha = (double *) malloc (N * sizeof (double)); + double * c = (double *) malloc (N * sizeof (double)); + double * z = (double *) malloc (N * sizeof (double)); + + if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0) { + cerr << "Internal Cyclic Spline Error: failed to allocate working space\n"; + exit(GSL_ENOMEM); + } + else + { + size_t i, j; + double sum = 0.0; + + /* factor */ + + if (N == 1) + { + x[0] = b[0] / diag[0]; + return GSL_SUCCESS; + } + + alpha[0] = diag[0]; + gamma[0] = offdiag[0] / alpha[0]; + delta[0] = offdiag[o_stride * (N-1)] / alpha[0]; + + if (alpha[0] == 0) { + status = GSL_EZERODIV; + } + + for (i = 1; i < N - 2; i++) + { + alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1]; + gamma[i] = offdiag[o_stride * i] / alpha[i]; + delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i]; + if (alpha[i] == 0) { + status = GSL_EZERODIV; + } + } + + for (i = 0; i < N - 2; i++) + { + sum += alpha[i] * delta[i] * delta[i]; + } + + alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3]; + + gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2]; + + alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2]; + + /* update */ + z[0] = b[0]; + for (i = 1; i < N - 1; i++) + { + z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1]; + } + sum = 0.0; + for (i = 0; i < N - 2; i++) + { + sum += delta[i] * z[i]; + } + z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2]; + for (i = 0; i < N; i++) + { + c[i] = z[i] / alpha[i]; + } + + /* backsubstitution */ + x[x_stride * (N - 1)] = c[N - 1]; + x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)]; + if (N >= 3) + { + for (i = N - 3, j = 0; j <= N - 3; j++, i--) + { + x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)]; + } + } + } + + if (z != 0) + free (z); + if (c != 0) + free (c); + if (alpha != 0) + free (alpha); + if (gamma != 0) + free (gamma); + if (delta != 0) + free (delta); + + if (status == GSL_EZERODIV) { + cerr <<"Internal Cyclic Spline Error: Matrix must be positive definite.\n"; + exit(GSL_EZERODIV); + } + + return status; +} //solve_cyc_tridiag() + +/* ---------------------------------------------------------------------- + spline and splint routines modified from Numerical Recipes +------------------------------------------------------------------------- */ + +static void cyc_spline(double const *xa, + double const *ya, + int n, + double period, + double *y2a) +{ + + double *diag = new double[n]; + double *offdiag = new double[n]; + double *rhs = new double[n]; + double xa_im1, xa_ip1; + + // In the cyclic case, there are n equations with n unknows. + // The for loop sets up the equations we need to solve. + // Later we invoke the GSL tridiagonal matrix solver to solve them. + + for(int i=0; i < n; i++) { + + // I have to lookup xa[i+1] and xa[i-1]. This gets tricky because of + // periodic boundary conditions. We handle that now. + int im1 = i-1; + if (im1<0) { + im1 += n; + xa_im1 = xa[im1] - period; + } + else + xa_im1 = xa[im1]; + + int ip1 = i+1; + if (ip1>=n) { + ip1 -= n; + xa_ip1 = xa[ip1] + period; + } + else + xa_ip1 = xa[ip1]; + + // Recall that we want to find the y2a[] parameters (there are n of them). + // To solve for them, we have a linear equation with n unknowns + // (in the cyclic case that is). For details, the non-cyclic case is + // explained in equation 3.3.7 in Numerical Recipes in C, p. 115. + diag[i] = (xa_ip1 - xa_im1) / 3.0; + offdiag[i] = (xa_ip1 - xa[i]) / 6.0; + rhs[i] = ((ya[ip1] - ya[i]) / (xa_ip1 - xa[i])) - + ((ya[i] - ya[im1]) / (xa[i] - xa_im1)); + } + + // Ordinarily we would have to invert a large matrix (with potentially + // thousands of rows and columns). However because this matix happens + // to be tridiagonal (and cyclic), we can use the following cheap method: + solve_cyc_tridiag(diag, 1, + offdiag, 1, + rhs, 1, + y2a, 1, + n); + + delete [] diag; + delete [] offdiag; + delete [] rhs; + +} // cyc_spline() + +/* ---------------------------------------------------------------------- */ + +// cyc_splint(): Evaluates a spline at position x, with n control +// points located at xa[], ya[], and with parameters y2a[] +// The xa[] must be monotonically increasing and their +// range should not exceed period (ie xa[n-1] < xa[0] + period). +// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] +// "period" is typically 2*PI. +static double cyc_splint(double const *xa, + double const *ya, + double const *y2a, + int n, + double period, + double x) +{ + int klo = -1; + int khi = n; + int k; + double xlo = xa[n-1] - period; + double xhi = xa[0] + period; + while (khi-klo > 1) { + k = (khi+klo) >> 1; //(k=(khi+klo)/2) + if (xa[k] > x) { + khi = k; + xhi = xa[k]; + } + else { + klo = k; + xlo = xa[k]; + } + } + + if (khi == n) khi = 0; + if (klo ==-1) klo = n-1; + + double h = xhi-xlo; + double a = (xhi-x) / h; + double b = (x-xlo) / h; + double y = a*ya[klo] + b*ya[khi] + + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; + + return y; + +} // cyc_splint() + +// cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x, +// with n control points at xa[], ya[], with parameters y2a[]. +// The xa[] must be monotonically increasing and their +// range should not exceed period (ie xa[n-1] < xa[0] + period). +// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] +// "period" is typically 2*PI. +static double cyc_splintD(double const *xa, + double const *ya, + double const *y2a, + int n, + double period, + double x) +{ + int klo = -1; + int khi = n; // (not n-1) + int k; + double xlo = xa[n-1] - period; + double xhi = xa[0] + period; + while (khi-klo > 1) { + k = (khi+klo) >> 1; //(k=(khi+klo)/2) + if (xa[k] > x) { + khi = k; + xhi = xa[k]; + } + else { + klo = k; + xlo = xa[k]; + } + } + + if (khi == n) khi = 0; + if (klo ==-1) klo = n-1; + + double yhi = ya[khi]; + double ylo = ya[klo]; + double h = xhi-xlo; + double g = yhi-ylo; + double a = (xhi-x) / h; + double b = (x-xlo) / h; + // Formula below taken from equation 3.3.5 of "numerical recipes in c" + // "yD" = the derivative of y + double yD = g/h - ( (3.0*a*a-1.0)*y2a[klo] - (3.0*b*b-1.0)*y2a[khi] ) * h/6.0; + // For rerefence: y = a*ylo + b*yhi + + // ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; + + return yD; + +} // cyc_splintD() + +// -------------------------------------------- +// ------- Calculate the dihedral angle ------- +// -------------------------------------------- +static const int g_dim=3; + +static double Phi(double const *x1, //array holding x,y,z coords atom 1 + double const *x2, // : : : : 2 + double const *x3, // : : : : 3 + double const *x4, // : : : : 4 + Domain *domain, //<-periodic boundary information + // The following arrays are of doubles with g_dim elements. + // (g_dim is a constant known at compile time, usually 3). + // Their contents is calculated by this function. + // Space for these vectors must be allocated in advance. + // (This is not hidden internally because these vectors + // may be needed outside the function, later on.) + double *vb12, // will store x2-x1 + double *vb23, // will store x3-x2 + double *vb34, // will store x4-x3 + double *n123, // will store normal to plane x1,x2,x3 + double *n234) // will store normal to plane x2,x3,x4 +{ + + for (int d=0; d < g_dim; ++d) { + vb12[d] = x2[d] - x1[d]; // 1st bond + vb23[d] = x3[d] - x2[d]; // 2nd bond + vb34[d] = x4[d] - x3[d]; // 3rd bond + } + + //Consider periodic boundary conditions: + domain->minimum_image(vb12[0],vb12[1],vb12[2]); + domain->minimum_image(vb23[0],vb23[1],vb23[2]); + domain->minimum_image(vb34[0],vb34[1],vb34[2]); + + //--- Compute the normal to the planes formed by atoms 1,2,3 and 2,3,4 --- + + cross3(vb23, vb12, n123); // <- n123=vb23 x vb12 + cross3(vb23, vb34, n234); // <- n234=vb23 x vb34 + + norm3(n123); + norm3(n234); + + double cos_phi = -dot3(n123, n234); + + if (cos_phi > 1.0) + cos_phi = 1.0; + else if (cos_phi < -1.0) + cos_phi = -1.0; + + double phi = acos(cos_phi); + + if (dot3(n123, vb34) > 0.0) { + phi = -phi; //(Note: Negative dihedral angles are possible only in 3-D.) + phi += MY_2PI; //<- This insure phi is always in the range 0 to 2*PI + } + return phi; +} // DihedralTable::Phi() + /* ---------------------------------------------------------------------- */ @@ -56,6 +407,7 @@ DihedralTable::DihedralTable(LAMMPS *lmp) : Dihedral(lmp) { ntables = 0; tables = NULL; + checkU_fname = checkF_fname = NULL; } /* ---------------------------------------------------------------------- */ @@ -64,6 +416,8 @@ DihedralTable::~DihedralTable() { for (int m = 0; m < ntables; m++) free_table(&tables[m]); memory->sfree(tables); + memory->sfree(checkU_fname); + memory->sfree(checkF_fname); if (allocated) { memory->destroy(setflag); @@ -133,9 +487,9 @@ void DihedralTable::compute(int eflag, int vflag) 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| + // proj12on23[d] = (vb23[d]/|vb23|) * dot3(vb12,vb23)/|vb12|*|vb23| double proj34on23[g_dim]; - // proj34on23[d] = (vb34[d]/|vb23|) * DotProduct(vb34,vb23)/|vb34|*|vb23| + // proj34on23[d] = (vb34[d]/|vb23|) * dot3(vb34,vb23)/|vb34|*|vb23| double perp12on23[g_dim]; // perp12on23[d] = v12[d] - proj12on23[d] double perp34on23[g_dim]; @@ -179,9 +533,9 @@ void DihedralTable::compute(int eflag, int vflag) double dphi_dx3[g_dim]; // d x[i1][d] double dphi_dx4[g_dim]; //where d=0,1,2 corresponds to x,y,z (if g_dim==3) - double dot123 = DotProduct(vb12, vb23); - double dot234 = DotProduct(vb23, vb34); - double L23sqr = DotProduct(vb23, vb23); + double dot123 = dot3(vb12, vb23); + double dot234 = dot3(vb23, vb34); + double L23sqr = dot3(vb23, vb23); double L23 = sqrt(L23sqr); // (central bond length) double inv_L23sqr = 0.0; double inv_L23 = 0.0; @@ -207,9 +561,9 @@ void DihedralTable::compute(int eflag, int vflag) // These two gradients point in the direction of n123 and n234, // and are scaled by the distances of atoms 1 and 4 from the central axis. // Distance of atom 1 to central axis: - double perp12on23_len = sqrt(DotProduct(perp12on23, perp12on23)); + double perp12on23_len = sqrt(dot3(perp12on23, perp12on23)); // Distance of atom 4 to central axis: - double perp34on23_len = sqrt(DotProduct(perp34on23, perp34on23)); + double perp34on23_len = sqrt(dot3(perp34on23, perp34on23)); double inv_perp12on23 = 0.0; if (perp12on23_len != 0.0) inv_perp12on23 = 1.0 / perp12on23_len; @@ -265,38 +619,12 @@ void DihedralTable::compute(int eflag, int vflag) dphi_dx3[d] = dphi123_dx3_coef*dphi_dx1[d] + dphi234_dx3_coef*dphi_dx4[d]; } - - - - #ifdef DIH_DEBUG_NUM - // ----- Numerical test? ----- - - cerr << " -- testing gradient for dihedral (n="<= TWOPI) { + if ((phihi - philo) >= MY_2PI) { string err_msg; err_msg = string("Dihedral table angle range must be < 2*PI radians (") + string(arg[2]) + string(")."); @@ -493,10 +821,10 @@ void DihedralTable::coeff(int narg, char **arg) // convert phi from degrees to radians if (tb->use_degrees) { for (int i=0; i < tb->ninput; i++) { - tb->phifile[i] *= PI/180.0; + tb->phifile[i] *= MY_PI/180.0; // I assume that if angles are in degrees, then the forces (f=dU/dphi) // are specified with "phi" in radians as well. - tb->ffile[i] *= 180.0/PI; + tb->ffile[i] *= 180.0/MY_PI; } } @@ -515,7 +843,7 @@ void DihedralTable::coeff(int narg, char **arg) for (int i=0; i < tb->ninput; i++) { double phi = tb->phifile[i]; // Add a multiple of 2*PI to phi until it lies in the range [0, 2*PI). - phi -= TWOPI * floor(phi/TWOPI); + phi -= MY_2PI * floor(phi/MY_2PI); phifile_tmp[i] = phi; efile_tmp[i] = tb->efile[i]; ffile_tmp[i] = tb->ffile[i]; @@ -523,7 +851,6 @@ void DihedralTable::coeff(int narg, char **arg) //There should only be at most one discontinuity, because we have //insured that the data was sorted before imaging, and because the //range of angle values does not exceed 2*PI. - assert(i_discontinuity == tb->ninput); //<-should only happen once i_discontinuity = i; } } @@ -533,17 +860,14 @@ void DihedralTable::coeff(int narg, char **arg) tb->phifile[I] = phifile_tmp[i]; tb->efile[I] = efile_tmp[i]; tb->ffile[I] = ffile_tmp[i]; - assert((I==0) || (tb->phifile[I] > tb->phifile[I-1])); I++; } for (int i = 0; i < i_discontinuity; i++) { tb->phifile[I] = phifile_tmp[i]; tb->efile[I] = efile_tmp[i]; tb->ffile[I] = ffile_tmp[i]; - assert((I==0) || (tb->phifile[I] > tb->phifile[I-1])); I++; } - assert(I == tb->ninput); } // spline read-in and compute r,e,f vectors within table @@ -559,10 +883,10 @@ void DihedralTable::coeff(int narg, char **arg) ofstream checkU_file; checkU_file.open(checkU_fname, ios::out); for (int i=0; i < tablength; i++) { - double phi = i*TWOPI/tablength; + double phi = i*MY_2PI/tablength; double u = tb->e[i]; if (tb->use_degrees) - phi *= 180.0/PI; + phi *= 180.0/MY_PI; checkU_file << phi << " " << u << "\n"; } checkU_file.close(); @@ -573,7 +897,7 @@ void DihedralTable::coeff(int narg, char **arg) checkF_file.open(checkF_fname, ios::out); for (int i=0; i < tablength; i++) { - double phi = i*TWOPI/tablength; + double phi = i*MY_2PI/tablength; double f; if ((tabstyle == SPLINE) && (tb->f_unspecified)) { double dU_dphi = @@ -585,17 +909,17 @@ void DihedralTable::coeff(int narg, char **arg) // -cyc_splintD() routine to calculate the derivative of the // energy spline, using the energy data (tb->e[]). // To be nice and report something, I do the same thing here.) - cyc_splintD(tb->phi, tb->e, tb->e2, tablength, TWOPI,phi); + cyc_splintD(tb->phi, tb->e, tb->e2, tablength, MY_2PI,phi); f = -dU_dphi; } else // Otherwise we calculated the tb->f[] array. Report its contents. f = tb->f[i]; if (tb->use_degrees) { - phi *= 180.0/PI; + phi *= 180.0/MY_PI; // If the user wants degree angle units, we should convert our // internal force tables (in energy/radians) to (energy/degrees) - f *= PI/180.0; + f *= MY_PI/180.0; } checkF_file << phi << " " << f << "\n"; } @@ -680,7 +1004,7 @@ void DihedralTable::free_table(Table *tb) /* ---------------------------------------------------------------------- read table file, only called by proc 0 ------------------------------------------------------------------------- */ - +static const int MAXLINE=2048; void DihedralTable::read_table(Table *tb, char *file, char *keyword) { char line[MAXLINE]; @@ -776,10 +1100,10 @@ void DihedralTable::spline_table(Table *tb) memory->create(tb->e2file,tb->ninput,"dihedral:e2file"); memory->create(tb->f2file,tb->ninput,"dihedral:f2file"); - cyc_spline(tb->phifile, tb->efile, tb->ninput, TWOPI, tb->e2file); + cyc_spline(tb->phifile, tb->efile, tb->ninput, MY_2PI, tb->e2file); if (! tb->f_unspecified) { - cyc_spline(tb->phifile, tb->ffile, tb->ninput, TWOPI, tb->f2file); + cyc_spline(tb->phifile, tb->ffile, tb->ninput, MY_2PI, tb->f2file); } // CHECK to help make sure the user calculated forces in a way @@ -801,14 +1125,14 @@ void DihedralTable::spline_table(Table *tb) int im1 = i-1; if (im1 < 0) { im1 += tb->ninput; - phi_im1 = tb->phifile[im1] - TWOPI; + phi_im1 = tb->phifile[im1] - MY_2PI; } else phi_im1 = tb->phifile[im1]; int ip1 = i+1; if (ip1 >= tb->ninput) { ip1 -= tb->ninput; - phi_ip1 = tb->phifile[ip1] + TWOPI; + phi_ip1 = tb->phifile[ip1] + MY_2PI; } else phi_ip1 = tb->phifile[ip1]; @@ -836,7 +1160,7 @@ void DihedralTable::spline_table(Table *tb) double f = -dU_dphi; // Alternately, we could use spline interpolation instead: // double f = - splintD(tb->phifile, tb->efile, tb->e2file, - // tb->ninput, TWOPI, tb->phifile[i]); + // tb->ninput, MY_2PI, tb->phifile[i]); // This is the way I originally did it, but I trust // the ugly simple linear way above better. // Recall this entire block of code doess not calculate @@ -867,7 +1191,7 @@ void DihedralTable::spline_table(Table *tb) void DihedralTable::compute_table(Table *tb) { //delta = table spacing in dihedral angle for tablength (cyclic) bins - tb->delta = TWOPI / tablength; + tb->delta = MY_2PI / tablength; tb->invdelta = 1.0/tb->delta; tb->deltasq6 = tb->delta*tb->delta / 6.0; @@ -889,9 +1213,9 @@ void DihedralTable::compute_table(Table *tb) for (int i = 0; i < tablength; i++) { double phi = i*tb->delta; tb->phi[i] = phi; - tb->e[i]= cyc_splint(tb->phifile,tb->efile,tb->e2file,tb->ninput,TWOPI,phi); + tb->e[i]= cyc_splint(tb->phifile,tb->efile,tb->e2file,tb->ninput,MY_2PI,phi); if (! tb->f_unspecified) - tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,TWOPI,phi); + tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi); } if (tabstyle == LINEAR) { @@ -918,9 +1242,9 @@ void DihedralTable::compute_table(Table *tb) } } // if (tabstyle == LINEAR) - cyc_spline(tb->phi, tb->e, tablength, TWOPI, tb->e2); + cyc_spline(tb->phi, tb->e, tablength, MY_2PI, tb->e2); if (! tb->f_unspecified) - cyc_spline(tb->phi, tb->f, tablength, TWOPI, tb->f2); + cyc_spline(tb->phi, tb->f, tablength, MY_2PI, tb->f2); } @@ -936,8 +1260,6 @@ void DihedralTable::param_extract(Table *tb, char *line) tb->ninput = 0; tb->f_unspecified = false; //default tb->use_degrees = true; //default - strcpy(checkU_fname, ""); - strcpy(checkF_fname, ""); char *word = strtok(line," \t\n\r\f"); while (word) { @@ -956,11 +1278,15 @@ void DihedralTable::param_extract(Table *tb, char *line) } else if (strcmp(word,"CHECKU") == 0) { word = strtok(NULL," \t\n\r\f"); - strncpy(checkU_fname, word, MAXLINE); + memory->sfree(checkU_fname); + memory->create(checkU_fname,strlen(word),"dihedral_table:checkU"); + strcpy(checkU_fname, word); } else if (strcmp(word,"CHECKF") == 0) { word = strtok(NULL," \t\n\r\f"); - strncpy(checkF_fname, word, MAXLINE); + memory->sfree(checkF_fname); + memory->create(checkF_fname,strlen(word),"dihedral_table:checkF"); + strcpy(checkF_fname, word); } // COMMENTING OUT: equilibrium angles are not supported //else if (strcmp(word,"EQ") == 0) { @@ -1012,316 +1338,3 @@ void DihedralTable::bcast_table(Table *tb) } -namespace LAMMPS_NS { -namespace DIHEDRAL_TABLE_NS { - - -// ------------------------------------------------------------------- -// --------- The function was stolen verbatim from the --------- -// --------- GNU Scientific Library (GSL, version 1.15) --------- -// ------------------------------------------------------------------- - -/* Author: Gerard Jungman */ -/* for description of method see [Engeln-Mullges + Uhlig, p. 96] - * - * diag[0] offdiag[0] 0 ..... offdiag[N-1] - * offdiag[0] diag[1] offdiag[1] ..... - * 0 offdiag[1] diag[2] - * 0 0 offdiag[2] ..... - * ... ... - * offdiag[N-1] ... - * - */ -// -- (A non-symmetric version of this function is also available.) -- - -enum { //GSL status return codes. - GSL_FAILURE = -1, - GSL_SUCCESS = 0, - GSL_ENOMEM = 8, - GSL_EZERODIV = 12, - GSL_EBADLEN = 19 -}; - - -int -solve_cyc_tridiag( - const double diag[], size_t d_stride, - const double offdiag[], size_t o_stride, - const double b[], size_t b_stride, - double x[], size_t x_stride, - size_t N) -{ - int status = GSL_SUCCESS; - double * delta = (double *) malloc (N * sizeof (double)); - double * gamma = (double *) malloc (N * sizeof (double)); - double * alpha = (double *) malloc (N * sizeof (double)); - double * c = (double *) malloc (N * sizeof (double)); - double * z = (double *) malloc (N * sizeof (double)); - - if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0) { - cerr << "Internal Cyclic Spline Error: failed to allocate working space\n"; - exit(GSL_ENOMEM); - } - else - { - size_t i, j; - double sum = 0.0; - - /* factor */ - - if (N == 1) - { - x[0] = b[0] / diag[0]; - return GSL_SUCCESS; - } - - alpha[0] = diag[0]; - gamma[0] = offdiag[0] / alpha[0]; - delta[0] = offdiag[o_stride * (N-1)] / alpha[0]; - - if (alpha[0] == 0) { - status = GSL_EZERODIV; - } - - for (i = 1; i < N - 2; i++) - { - alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1]; - gamma[i] = offdiag[o_stride * i] / alpha[i]; - delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i]; - if (alpha[i] == 0) { - status = GSL_EZERODIV; - } - } - - for (i = 0; i < N - 2; i++) - { - sum += alpha[i] * delta[i] * delta[i]; - } - - alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3]; - - gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2]; - - alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2]; - - /* update */ - z[0] = b[0]; - for (i = 1; i < N - 1; i++) - { - z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1]; - } - sum = 0.0; - for (i = 0; i < N - 2; i++) - { - sum += delta[i] * z[i]; - } - z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2]; - for (i = 0; i < N; i++) - { - c[i] = z[i] / alpha[i]; - } - - /* backsubstitution */ - x[x_stride * (N - 1)] = c[N - 1]; - x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)]; - if (N >= 3) - { - for (i = N - 3, j = 0; j <= N - 3; j++, i--) - { - x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)]; - } - } - } - - if (z != 0) - free (z); - if (c != 0) - free (c); - if (alpha != 0) - free (alpha); - if (gamma != 0) - free (gamma); - if (delta != 0) - free (delta); - - if (status == GSL_EZERODIV) { - cerr <<"Internal Cyclic Spline Error: Matrix must be positive definite.\n"; - exit(GSL_EZERODIV); - } - - return status; -} //solve_cyc_tridiag() - - - -/* ---------------------------------------------------------------------- - spline and splint routines modified from Numerical Recipes -------------------------------------------------------------------------- */ - -void cyc_spline(double const *xa, - double const *ya, - int n, - double period, - double *y2a) -{ - - double *diag = new double[n]; - double *offdiag = new double[n]; - double *rhs = new double[n]; - double xa_im1, xa_ip1; - - // In the cyclic case, there are n equations with n unknows. - // The for loop sets up the equations we need to solve. - // Later we invoke the GSL tridiagonal matrix solver to solve them. - - for(int i=0; i < n; i++) { - - // I have to lookup xa[i+1] and xa[i-1]. This gets tricky because of - // periodic boundary conditions. We handle that now. - int im1 = i-1; - if (im1<0) { - im1 += n; - xa_im1 = xa[im1] - period; - } - else - xa_im1 = xa[im1]; - - int ip1 = i+1; - if (ip1>=n) { - ip1 -= n; - xa_ip1 = xa[ip1] + period; - } - else - xa_ip1 = xa[ip1]; - - // Recall that we want to find the y2a[] parameters (there are n of them). - // To solve for them, we have a linear equation with n unknowns - // (in the cyclic case that is). For details, the non-cyclic case is - // explained in equation 3.3.7 in Numerical Recipes in C, p. 115. - diag[i] = (xa_ip1 - xa_im1) / 3.0; - offdiag[i] = (xa_ip1 - xa[i]) / 6.0; - rhs[i] = ((ya[ip1] - ya[i]) / (xa_ip1 - xa[i])) - - ((ya[i] - ya[im1]) / (xa[i] - xa_im1)); - } - - // Ordinarily we would have to invert a large matrix (with potentially - // thousands of rows and columns). However because this matix happens - // to be tridiagonal (and cyclic), we can use the following cheap method: - solve_cyc_tridiag(diag, 1, - offdiag, 1, - rhs, 1, - y2a, 1, - n); - - delete [] diag; - delete [] offdiag; - delete [] rhs; - -} // cyc_spline() - - - -/* ---------------------------------------------------------------------- */ - -// cyc_splint(): Evaluates a spline at position x, with n control -// points located at xa[], ya[], and with parameters y2a[] -// The xa[] must be monotonically increasing and their -// range should not exceed period (ie xa[n-1] < xa[0] + period). -// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] -// "period" is typically 2*PI. -double cyc_splint(double const *xa, - double const *ya, - double const *y2a, - int n, - double period, - double x) -{ - int klo = -1; - int khi = n; - int k; - double xlo = xa[n-1] - period; - double xhi = xa[0] + period; - assert((xlo <= x) && (x <= xhi)); - while (khi-klo > 1) { - k = (khi+klo) >> 1; //(k=(khi+klo)/2) - if (xa[k] > x) { - khi = k; - xhi = xa[k]; - } - else { - klo = k; - xlo = xa[k]; - } - } - assert((xlo <= x) && (x <= xhi)); - - if (khi == n) khi = 0; - if (klo ==-1) klo = n-1; - - double h = xhi-xlo; - double a = (xhi-x) / h; - double b = (x-xlo) / h; - double y = a*ya[klo] + b*ya[khi] + - ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; - - return y; - -} // cyc_splint() - - - -// cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x, -// with n control points at xa[], ya[], with parameters y2a[]. -// The xa[] must be monotonically increasing and their -// range should not exceed period (ie xa[n-1] < xa[0] + period). -// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] -// "period" is typically 2*PI. -double cyc_splintD(double const *xa, - double const *ya, - double const *y2a, - int n, - double period, - double x) -{ - int klo = -1; - int khi = n; // (not n-1) - int k; - double xlo = xa[n-1] - period; - double xhi = xa[0] + period; - assert((xlo <= x) && (x <= xhi)); - while (khi-klo > 1) { - k = (khi+klo) >> 1; //(k=(khi+klo)/2) - if (xa[k] > x) { - khi = k; - xhi = xa[k]; - } - else { - klo = k; - xlo = xa[k]; - } - } - assert((xlo <= x) && (x <= xhi)); - - if (khi == n) khi = 0; - if (klo ==-1) klo = n-1; - - double yhi = ya[khi]; - double ylo = ya[klo]; - double h = xhi-xlo; - double g = yhi-ylo; - double a = (xhi-x) / h; - double b = (x-xlo) / h; - // Formula below taken from equation 3.3.5 of "numerical recipes in c" - // "yD" = the derivative of y - double yD = g/h - ( (3.0*a*a-1.0)*y2a[klo] - (3.0*b*b-1.0)*y2a[khi] ) * h/6.0; - // For rerefence: y = a*ylo + b*yhi + - // ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; - - return yD; - -} // cyc_splintD() - - - -} //namespace DIHEDRAL_TABLE_NS -} //namespace LAMMPS_NS diff --git a/src/USER-MISC/dihedral_table.h b/src/USER-MISC/dihedral_table.h index 3d8de3d61f..08bd21ba83 100644 --- a/src/USER-MISC/dihedral_table.h +++ b/src/USER-MISC/dihedral_table.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -23,10 +23,6 @@ DihedralStyle(table,DihedralTable) #ifndef LMP_DIHEDRAL_TABLE_H #define LMP_DIHEDRAL_TABLE_H - -#include -#include -#include "domain.h" #include "dihedral.h" namespace LAMMPS_NS { @@ -45,6 +41,8 @@ class DihedralTable : public Dihedral { protected: int tabstyle,tablength; // double *phi0; <- equilibrium angles not supported + char *checkU_fname; + char *checkF_fname; struct Table { int ninput; @@ -147,169 +145,8 @@ class DihedralTable : public Dihedral { } } // u_lookup() - - // Pre-allocated strings to store file names for debugging splines. (One day - // I would really like to rewrite everything and use C++ strings instead.) - static const int MAXLINE=2048; - char checkU_fname[MAXLINE]; - char checkF_fname[MAXLINE]; - }; //class DihedralTable -// ------------------------------------------------------------------------ -// The following auxiliary functions were left out of the -// DihedralTable class either because they require template parameters, -// or because they have nothing to do with dihedral angles. -// ------------------------------------------------------------------------ - -namespace DIHEDRAL_TABLE_NS { - -static const double PI = 3.1415926535897931; -static const double TWOPI = 6.2831853071795862; - -// Determine the array of "y2" parameters of a cyclic spline from its control -// points at positions x[] and y[]. (spline() must be invoked before splint()) -// The x[] positions should be sorted in order and not exceed period. -void cyc_spline(double const *xa, double const *ya, int n, - double period, double *y2a); - -// Evaluate a cyclic spline at position x with n control points at xa[], ya[], -// (The y2a array must be calculated using cyc_spline() above in advance.) -// x (and all the xa[] positions) should lie in the range from 0 to period. -// (Typically period = 2*PI, but this is optional.) -double cyc_splint(double const *xa, double const *ya, double const *y2a, - int n, double period, double x); - -// Evaluate the deriviative of a cyclic spline at position x: -double cyc_splintD(double const *xa, double const *ya, double const *y2a, - int n, double period, double x); - -// ----------------------------------------------------------- -// ---- some simple vector operations are defined below. ---- -// ----------------------------------------------------------- - -// --- g_dim --- As elsewhere in the LAMMPS code, coordinates here are -// represented as entries in an array, not as named variables "x" "y" "z". -// (I like this style.) In this spirit, the vector operations here are -// defined for vectors of arbitrary size. For this to work, the number -// of dimensions, "g_dim", must be known at compile time: -const int g_dim = 3; -// In LAMMPS at least, this constant is always 3, and is only used inside -// the dihedral code here. (It should not conflict with 2-D simulations.) -// Note: Compiler optimizations should eliminate any performance overhead -// associated with loops like "for (int i=0; i -inline _Real -DotProduct(_Real const *A, _Real const *B) -{ - _Real AdotB = 0.0; - for (int d=0; d < g_dim; ++d) - AdotB += A[d]*B[d]; - return AdotB; -} - -// Normalize() divides the components of the vector "v" by it's length. -// Normalize() silently ignores divide-by-zero errors but does not -// crash. (If "v" has length 0, then we replace v with the unit vector in -// an arbitrary direction,(1,0,...).) -// It returns the length of v (useful for checking if the operation succeeded). -template -inline _Real -Normalize(_Real *v) -{ - _Real length = sqrt(DotProduct(v,v)); - if (length != 0.0) - { - _Real one_over_length = 1.0 / length; - for (int d=0; d < g_dim; ++d) - v[d] *= one_over_length; - } - else { - v[0] = 1.0; - for (int d=1; d < g_dim; ++d) - v[d] = 0.0; - } - return length; -} - - -// CrossProduct(A,B,dest) computes the cross-product (A x B) -// and stores the result in "dest". -template -inline void -CrossProduct(_Real const *A, _Real const *B, _Real *dest) -{ - dest[0] = A[1]*B[2] - A[2]*B[1]; - dest[1] = A[2]*B[0] - A[0]*B[2]; - dest[2] = A[0]*B[1] - A[1]*B[0]; -} - - -// -------------------------------------------- -// ------- Calculate the dihedral angle ------- -// -------------------------------------------- - -inline double Phi(double const *x1, //array holding x,y,z coords atom 1 - double const *x2, // : : : : 2 - double const *x3, // : : : : 3 - double const *x4, // : : : : 4 - Domain *domain, //<-periodic boundary information - // The following arrays are of doubles with g_dim elements. - // (g_dim is a constant known at compile time, usually 3). - // Their contents is calculated by this function. - // Space for these vectors must be allocated in advance. - // (This is not hidden internally because these vectors - // may be needed outside the function, later on.) - double *vb12, // will store x2-x1 - double *vb23, // will store x3-x2 - double *vb34, // will store x4-x3 - double *n123, // will store normal to plane x1,x2,x3 - double *n234) // will store normal to plane x2,x3,x4 -{ - - for (int d=0; d < g_dim; ++d) { - vb12[d] = x2[d] - x1[d]; // 1st bond - vb23[d] = x3[d] - x2[d]; // 2nd bond - vb34[d] = x4[d] - x3[d]; // 3rd bond - } - - //Consider periodic boundary conditions: - domain->minimum_image(vb12[0],vb12[1],vb12[2]); - domain->minimum_image(vb23[0],vb23[1],vb23[2]); - domain->minimum_image(vb34[0],vb34[1],vb34[2]); - - //--- Compute the normal to the planes formed by atoms 1,2,3 and 2,3,4 --- - - CrossProduct(vb23, vb12, n123); // <- n123=vb23 x vb12 - CrossProduct(vb23, vb34, n234); // <- n234=vb23 x vb34 - - Normalize(n123); - Normalize(n234); - - double cos_phi = -DotProduct(n123, n234); - - if (cos_phi > 1.0) - cos_phi = 1.0; - else if (cos_phi < -1.0) - cos_phi = -1.0; - - double phi = acos(cos_phi); - - if (DotProduct(n123, vb34) > 0.0) { - phi = -phi; //(Note: Negative dihedral angles are possible only in 3-D.) - phi += TWOPI; //<- This insure phi is always in the range 0 to 2*PI - } - return phi; -} // DihedralTable::Phi() - - -} // namespace DIHEDRAL_TABLE_NS - } // namespace LAMMPS_NS diff --git a/src/USER-MISC/improper_ring.cpp b/src/USER-MISC/improper_ring.cpp index fbbe9dd75b..031821df82 100644 --- a/src/USER-MISC/improper_ring.cpp +++ b/src/USER-MISC/improper_ring.cpp @@ -48,11 +48,13 @@ #include "force.h" #include "update.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; #define TOLERANCE 0.05 #define SMALL 0.001 @@ -169,7 +171,7 @@ void ImproperRing::compute(int eflag, int vflag) /* Append the current angle to the sum of angle differences. */ angle_summer += (bend_angle[icomb] - chi[type]); } - if (eflag) eimproper = (1.0/6.0) *k[type] * pow(angle_summer,6.0); + if (eflag) eimproper = (1.0/6.0) *k[type] * powint(angle_summer,6); /* printf("The tags: %d-%d-%d-%d, of type %d .\n",atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4],type); // printf("The coordinates of the first: %f, %f, %f.\n", x[i1][0], x[i1][1], x[i1][2]); @@ -183,7 +185,7 @@ void ImproperRing::compute(int eflag, int vflag) /* Force calculation acting on all atoms. Calculate the derivatives of the potential. */ - angfac = k[type] * pow(angle_summer,5.0); + angfac = k[type] * powint(angle_summer,5); f1[0] = 0.0; f1[1] = 0.0; f1[2] = 0.0; f3[0] = 0.0; f3[1] = 0.0; f3[2] = 0.0; diff --git a/src/USER-MISC/pair_dipole_sf.cpp b/src/USER-MISC/pair_dipole_sf.cpp index 61534e6ec3..34c7532898 100644 --- a/src/USER-MISC/pair_dipole_sf.cpp +++ b/src/USER-MISC/pair_dipole_sf.cpp @@ -259,8 +259,9 @@ void PairDipoleSF::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { - ecoul = qtmp * q[j] * rinv * - pow((1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])),2.0); + ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])); + ecoul *= ecoul; + ecoul *= qtmp * q[j] * rinv; if (mu[i][3] > 0.0 && mu[j][3] > 0.0) ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); if (mu[i][3] > 0.0 && q[j] != 0.0)