diff --git a/doc/html/Manual.html b/doc/html/Manual.html index 5c7ec2f31d..24875c671b 100644 --- a/doc/html/Manual.html +++ b/doc/html/Manual.html @@ -134,8 +134,8 @@

LAMMPS Documentation

-
-

27 Aug 2016 version

+
+

7 Sep 2016 version

Version info:

diff --git a/doc/html/Section_packages.html b/doc/html/Section_packages.html index f882ee234b..fedd3cc1ae 100644 --- a/doc/html/Section_packages.html +++ b/doc/html/Section_packages.html @@ -493,7 +493,7 @@ name links to a sub-section below with more details.

python lib/python -REAX +REAX ReaxFF potential Aidan Thompson (Sandia) pair_style reax diff --git a/doc/html/compute_rdf.html b/doc/html/compute_rdf.html index c8d486553f..ef1d1d8cfc 100644 --- a/doc/html/compute_rdf.html +++ b/doc/html/compute_rdf.html @@ -227,7 +227,7 @@ atoms of type jtypeN.

to a file is to use the fix ave/time command, for example:

compute myRDF all rdf 50
-fix 1 all ave/time 100 1 100 c_myRDF file tmp.rdf mode vector
+fix 1 all ave/time 100 1 100 c_myRDF[*] file tmp.rdf mode vector

Output info:

diff --git a/doc/html/fix_lb_pc.html b/doc/html/fix_lb_pc.html index 8a5281c6c3..4fe4f72ced 100644 --- a/doc/html/fix_lb_pc.html +++ b/doc/html/fix_lb_pc.html @@ -146,7 +146,7 @@

Description

Update the positions and velocities of the individual particles described by group-ID, experiencing velocity-dependent hydrodynamic -forces, using the integration algorithm described in Mackay et al.. This integration algorithm should only be used if a +forces, using the integration algorithm described in Mackay et al.. This integration algorithm should only be used if a user-specified value for the force-coupling constant used in fix lb/fluid has been set; do not use this integration algorithm if the force coupling constant has been set by default.

diff --git a/doc/html/pair_thole.html b/doc/html/pair_thole.html index b2742247f9..e3a5f8b9cb 100644 --- a/doc/html/pair_thole.html +++ b/doc/html/pair_thole.html @@ -189,8 +189,8 @@ short distances by a function

\frac{s_{ij} r_{ij} }{2} \right) \exp \left( - s_{ij} r_{ij} \right) \end{equation}\]

This function results from an adaptation to point charges -(Noskov) of the dipole screening scheme originally proposed -by Thole. The scaling coefficient \(s_{ij}\) is determined +(Noskov) of the dipole screening scheme originally proposed +by Thole. The scaling coefficient \(s_{ij}\) is determined by the polarizability of the atoms, \(\alpha_i\), and by a Thole damping parameter \(a\). This Thole damping parameter usually takes a value of 2.6, but in certain force fields the value can depend upon diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt index adf6b5c47e..b7edc5ea99 100644 --- a/doc/src/Manual.txt +++ b/doc/src/Manual.txt @@ -1,7 +1,7 @@ LAMMPS Users Manual - + @@ -21,7 +21,7 @@

LAMMPS Documentation :c,h3 -27 Aug 2016 version :c,h4 +7 Sep 2016 version :c,h4 Version info: :h4 diff --git a/src/KOKKOS/pair_eam_kokkos.cpp b/src/KOKKOS/pair_eam_kokkos.cpp index 37a26145e0..d91da280ac 100755 --- a/src/KOKKOS/pair_eam_kokkos.cpp +++ b/src/KOKKOS/pair_eam_kokkos.cpp @@ -340,81 +340,64 @@ void PairEAMKokkos::array2spline() rdr = 1.0/dr; rdrho = 1.0/drho; - tdual_ffloat4 k_frho_spline_a = tdual_ffloat4("pair:frho_a",nfrho,nrho+1); - tdual_ffloat4 k_rhor_spline_a = tdual_ffloat4("pair:rhor_a",nrhor,nr+1); - tdual_ffloat4 k_z2r_spline_a = tdual_ffloat4("pair:z2r_a",nz2r,nr+1); - tdual_ffloat4 k_frho_spline_b = tdual_ffloat4("pair:frho_b",nfrho,nrho+1); - tdual_ffloat4 k_rhor_spline_b = tdual_ffloat4("pair:rhor_b",nrhor,nr+1); - tdual_ffloat4 k_z2r_spline_b = tdual_ffloat4("pair:z2r_b",nz2r,nr+1); + tdual_ffloat_2d_n7 k_frho_spline = tdual_ffloat_2d_n7("pair:frho",nfrho,nrho+1); + tdual_ffloat_2d_n7 k_rhor_spline = tdual_ffloat_2d_n7("pair:rhor",nrhor,nr+1); + tdual_ffloat_2d_n7 k_z2r_spline = tdual_ffloat_2d_n7("pair:z2r",nz2r,nr+1); - t_host_ffloat4 h_frho_spline_a = k_frho_spline_a.h_view; - t_host_ffloat4 h_rhor_spline_a = k_rhor_spline_a.h_view; - t_host_ffloat4 h_z2r_spline_a = k_z2r_spline_a.h_view; - t_host_ffloat4 h_frho_spline_b = k_frho_spline_b.h_view; - t_host_ffloat4 h_rhor_spline_b = k_rhor_spline_b.h_view; - t_host_ffloat4 h_z2r_spline_b = k_z2r_spline_b.h_view; + t_host_ffloat_2d_n7 h_frho_spline = k_frho_spline.h_view; + t_host_ffloat_2d_n7 h_rhor_spline = k_rhor_spline.h_view; + t_host_ffloat_2d_n7 h_z2r_spline = k_z2r_spline.h_view; for (int i = 0; i < nfrho; i++) - interpolate(nrho,drho,frho[i],h_frho_spline_a,h_frho_spline_b,i); - k_frho_spline_a.template modify(); - k_frho_spline_a.template sync(); - k_frho_spline_b.template modify(); - k_frho_spline_b.template sync(); + interpolate(nrho,drho,frho[i],h_frho_spline,i); + k_frho_spline.template modify(); + k_frho_spline.template sync(); for (int i = 0; i < nrhor; i++) - interpolate(nr,dr,rhor[i],h_rhor_spline_a,h_rhor_spline_b,i); - k_rhor_spline_a.template modify(); - k_rhor_spline_a.template sync(); - k_rhor_spline_b.template modify(); - k_rhor_spline_b.template sync(); + interpolate(nr,dr,rhor[i],h_rhor_spline,i); + k_rhor_spline.template modify(); + k_rhor_spline.template sync(); for (int i = 0; i < nz2r; i++) - interpolate(nr,dr,z2r[i],h_z2r_spline_a,h_z2r_spline_b,i); - k_z2r_spline_a.template modify(); - k_z2r_spline_a.template sync(); - k_z2r_spline_b.template modify(); - k_z2r_spline_b.template sync(); - - d_frho_spline_a = k_frho_spline_a.d_view; - d_rhor_spline_a = k_rhor_spline_a.d_view; - d_z2r_spline_a = k_z2r_spline_a.d_view; - d_frho_spline_b = k_frho_spline_b.d_view; - d_rhor_spline_b = k_rhor_spline_b.d_view; - d_z2r_spline_b = k_z2r_spline_b.d_view; - + interpolate(nr,dr,z2r[i],h_z2r_spline,i); + k_z2r_spline.template modify(); + k_z2r_spline.template sync(); + d_frho_spline = k_frho_spline.d_view; + d_rhor_spline = k_rhor_spline.d_view; + d_z2r_spline = k_z2r_spline.d_view; } /* ---------------------------------------------------------------------- */ template -void PairEAMKokkos::interpolate(int n, double delta, double *f, t_host_ffloat4 h_spline_a, t_host_ffloat4 h_spline_b, int i) +void PairEAMKokkos::interpolate(int n, double delta, double *f, t_host_ffloat_2d_n7 h_spline, int i) { - for (int m = 1; m <= n; m++) h_spline_b(i,m).w = f[m]; + for (int m = 1; m <= n; m++) h_spline(i,m,6) = f[m]; - h_spline_b(i,1).z = h_spline_b(i,2).w - h_spline_b(i,1).w; - h_spline_b(i,2).z = 0.5 * (h_spline_b(i,3).w-h_spline_b(i,1).w); - h_spline_b(i,n-1).z = 0.5 * (h_spline_b(i,n).w-h_spline_b(i,n-2).w); - h_spline_b(i,n).z = h_spline_b(i,n).w - h_spline_b(i,n-1).w; + h_spline(i,1,5) = h_spline(i,2,6) - h_spline(i,1,6); + h_spline(i,2,5) = 0.5 * (h_spline(i,3,6)-h_spline(i,1,6)); + h_spline(i,n-1,5) = 0.5 * (h_spline(i,n,6)-h_spline(i,n-2,6)); + h_spline(i,n,5) = h_spline(i,n,6) - h_spline(i,n-1,6); for (int m = 3; m <= n-2; m++) - h_spline_b(i,m).z = ((h_spline_b(i,m-2).w-h_spline_b(i,m+2).w) + - 8.0*(h_spline_b(i,m+1).w-h_spline_b(i,m-1).w)) / 12.0; + h_spline(i,m,5) = ((h_spline(i,m-2,6)-h_spline(i,m+2,6)) + + 8.0*(h_spline(i,m+1,6)-h_spline(i,m-1,6))) / 12.0; for (int m = 1; m <= n-1; m++) { - h_spline_b(i,m).y = 3.0*(h_spline_b(i,m+1).w-h_spline_b(i,m).w) - - 2.0*h_spline_b(i,m).z - h_spline_b(i,m+1).z; - h_spline_b(i,m).x = h_spline_b(i,m).z + h_spline_b(i,m+1).z - - 2.0*(h_spline_b(i,m+1).w-h_spline_b(i,m).w); + h_spline(i,m,4) = 3.0*(h_spline(i,m+1,6)-h_spline(i,m,6)) - + 2.0*h_spline(i,m,5) - h_spline(i,m+1,5); + h_spline(i,m,3) = h_spline(i,m,5) + h_spline(i,m+1,5) - + 2.0*(h_spline(i,m+1,6)-h_spline(i,m,6)); } - h_spline_b(i,n).y = 0.0; - h_spline_b(i,n).x = 0.0; + h_spline(i,n,4) = 0.0; + h_spline(i,n,3) = 0.0; for (int m = 1; m <= n; m++) { - h_spline_a(i,m).z = h_spline_b(i,m).z/delta; - h_spline_a(i,m).y = 2.0*h_spline_b(i,m).y/delta; - h_spline_a(i,m).x = 3.0*h_spline_b(i,m).x/delta; + h_spline(i,m,2) = h_spline(i,m,5)/delta; + h_spline(i,m,1) = 2.0*h_spline(i,m,4)/delta; + h_spline(i,m,0) = 3.0*h_spline(i,m,3)/delta; } } @@ -558,13 +541,12 @@ void PairEAMKokkos::operator()(TagPairEAMKernelA::operator()(TagPairEAMKernelB, const int & p -= m; p = MIN(p,1.0); const int d_type2frho_i = d_type2frho[itype]; - const F_FLOAT4 frho = d_frho_spline_a(d_type2frho_i,m); - d_fp[i] = (frho.x*p + frho.y)*p + frho.z; + d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2); if (EFLAG) { - const F_FLOAT4 frho_b = d_frho_spline_b(d_type2frho_i,m); - F_FLOAT phi = ((frho_b.x*p + frho_b.y)*p + frho_b.z)*p + frho_b.w; + F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p + + d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6); if (d_rho[i] > rhomax) phi += d_fp[i] * (d_rho[i]-rhomax); if (eflag_global) ev.evdwl += phi; if (eflag_atom) d_eatom[i] += phi; @@ -651,8 +632,8 @@ void PairEAMKokkos::operator()(TagPairEAMKernelAB, const int p -= m; p = MIN(p,1.0); const int d_type2rhor_ji = d_type2rhor(jtype,itype); - const F_FLOAT4 rhor = d_rhor_spline_b(d_type2rhor_ji,m); - rhotmp += ((rhor.x*p + rhor.y)*p + rhor.z)*p + rhor.w; + rhotmp += ((d_rhor_spline(d_type2rhor_ji,m,3)*p + d_rhor_spline(d_type2rhor_ji,m,4))*p + + d_rhor_spline(d_type2rhor_ji,m,5))*p + d_rhor_spline(d_type2rhor_ji,m,6); } } @@ -669,11 +650,10 @@ void PairEAMKokkos::operator()(TagPairEAMKernelAB, const int p -= m; p = MIN(p,1.0); const int d_type2frho_i = d_type2frho[itype]; - const F_FLOAT4 frho = d_frho_spline_a(d_type2frho_i,m); - d_fp[i] = (frho.x*p + frho.y)*p + frho.z; + d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2); if (EFLAG) { - const F_FLOAT4 frho_b = d_frho_spline_b(d_type2frho_i,m); - F_FLOAT phi = ((frho_b.x*p + frho_b.y)*p + frho_b.z)*p + frho_b.w; + F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p + + d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6); if (d_rho[i] > rhomax) phi += d_fp[i] * (d_rho[i]-rhomax); if (eflag_global) ev.evdwl += phi; if (eflag_atom) d_eatom[i] += phi; @@ -742,18 +722,16 @@ void PairEAMKokkos::operator()(TagPairEAMKernelC; #ifdef KOKKOS_HAVE_CUDA template class PairEAMKokkos; #endif -} - +} \ No newline at end of file diff --git a/src/KOKKOS/pair_eam_kokkos.h b/src/KOKKOS/pair_eam_kokkos.h index b3d33bb37f..89d1d61b6b 100755 --- a/src/KOKKOS/pair_eam_kokkos.h +++ b/src/KOKKOS/pair_eam_kokkos.h @@ -136,14 +136,14 @@ class PairEAMKokkos : public PairEAM { DAT::t_int_2d_randomread d_type2rhor; DAT::t_int_2d_randomread d_type2z2r; - typedef Kokkos::DualView tdual_ffloat4; - typedef typename tdual_ffloat4::t_dev_const_randomread t_ffloat4_randomread; - typedef typename tdual_ffloat4::t_host t_host_ffloat4; + typedef Kokkos::DualView tdual_ffloat_2d_n7; + typedef typename tdual_ffloat_2d_n7::t_dev_const_randomread t_ffloat_2d_n7_randomread; + typedef typename tdual_ffloat_2d_n7::t_host t_host_ffloat_2d_n7; - t_ffloat4_randomread d_frho_spline_a, d_frho_spline_b; - t_ffloat4_randomread d_rhor_spline_a, d_rhor_spline_b; - t_ffloat4_randomread d_z2r_spline_a, d_z2r_spline_b; - void interpolate(int, double, double *, t_host_ffloat4, t_host_ffloat4, int); + t_ffloat_2d_n7_randomread d_frho_spline; + t_ffloat_2d_n7_randomread d_rhor_spline; + t_ffloat_2d_n7_randomread d_z2r_spline; + void interpolate(int, double, double *, t_host_ffloat_2d_n7, int); virtual void file2array(); void array2spline(); diff --git a/src/compute_omega_chunk.cpp b/src/compute_omega_chunk.cpp index ca4af7931f..f2a8394964 100644 --- a/src/compute_omega_chunk.cpp +++ b/src/compute_omega_chunk.cpp @@ -18,11 +18,14 @@ #include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" +#include "math_extra.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; +#define EPSILON 1.0e-6 + /* ---------------------------------------------------------------------- */ ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) : @@ -192,49 +195,97 @@ void ComputeOmegaChunk::compute_array() MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk, MPI_DOUBLE,MPI_SUM,world); - // compute omega for each chunk from L = Iw, inverting I to solve for w + // compute omega for each chunk - double ione[3][3],inverse[3][3]; + double determinant,invdeterminant; + double idiag[3],ex[3],ey[3],ez[3],cross[3]; + double ione[3][3],inverse[3][3],evectors[3][3]; + double *iall,*mall; for (m = 0; m < nchunk; m++) { - ione[0][0] = inertiaall[m][0]; - ione[1][1] = inertiaall[m][1]; - ione[2][2] = inertiaall[m][2]; - ione[0][1] = inertiaall[m][3]; - ione[1][2] = inertiaall[m][4]; - ione[0][2] = inertiaall[m][5]; - ione[1][0] = ione[0][1]; - ione[2][1] = ione[1][2]; - ione[2][0] = ione[0][2]; - inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1]; - inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]); - inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1]; + // determinant = triple product of rows of inertia matrix - inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]); - inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0]; - inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]); + iall = &inertiaall[m][0]; + determinant = iall[0] * (iall[1]*iall[2] - iall[4]*iall[4]) + + iall[3] * (iall[4]*iall[5] - iall[3]*iall[2]) + + iall[5] * (iall[3]*iall[4] - iall[1]*iall[5]); - inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0]; - inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]); - inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0]; + ione[0][0] = iall[0]; + ione[1][1] = iall[1]; + ione[2][2] = iall[2]; + ione[0][1] = ione[1][0] = iall[3]; + ione[1][2] = ione[2][1] = iall[4]; + ione[0][2] = ione[2][0] = iall[5]; - double determinant = ione[0][0]*ione[1][1]*ione[2][2] + - ione[0][1]*ione[1][2]*ione[2][0] + ione[0][2]*ione[1][0]*ione[2][1] - - ione[0][0]*ione[1][2]*ione[2][1] - ione[0][1]*ione[1][0]*ione[2][2] - - ione[2][0]*ione[1][1]*ione[0][2]; + // non-singular I matrix + // use L = Iw, inverting I to solve for w - if (determinant > 0.0) + if (determinant > EPSILON) { + inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1]; + inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]); + inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1]; + + inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]); + inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0]; + inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]); + + inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0]; + inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]); + inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0]; + + invdeterminant = 1.0/determinant; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) - inverse[i][j] /= determinant; + inverse[i][j] *= invdeterminant; + + mall = &angmomall[m][0]; + omega[m][0] = inverse[0][0]*mall[0] + inverse[0][1]*mall[1] + + inverse[0][2]*mall[2]; + omega[m][1] = inverse[1][0]*mall[0] + inverse[1][1]*mall[1] + + inverse[1][2]*mall[2]; + omega[m][2] = inverse[2][0]*mall[0] + inverse[2][1]*mall[1] + + inverse[2][2]*mall[2]; - omega[m][0] = inverse[0][0]*angmom[m][0] + inverse[0][1]*angmom[m][1] + - inverse[0][2]*angmom[m][2]; - omega[m][1] = inverse[1][0]*angmom[m][0] + inverse[1][1]*angmom[m][1] + - inverse[1][2]*angmom[m][2]; - omega[m][2] = inverse[2][0]*angmom[m][0] + inverse[2][1]*angmom[m][1] + - inverse[2][2]*angmom[i][2]; + // handle each (nearly) singular I matrix + // due to 2-atom chunk or linear molecule + // use jacobi() and angmom_to_omega() to calculate valid omega + + } else { + int ierror = MathExtra::jacobi(ione,idiag,evectors); + if (ierror) error->all(FLERR, + "Insufficient Jacobi rotations for omega/chunk"); + + ex[0] = evectors[0][0]; + ex[1] = evectors[1][0]; + ex[2] = evectors[2][0]; + ey[0] = evectors[0][1]; + ey[1] = evectors[1][1]; + ey[2] = evectors[2][1]; + ez[0] = evectors[0][2]; + ez[1] = evectors[1][2]; + ez[2] = evectors[2][2]; + + // enforce 3 evectors as a right-handed coordinate system + // flip 3rd vector if needed + + MathExtra::cross3(ex,ey,cross); + if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez); + + // if any principal moment < scaled EPSILON, set to 0.0 + + double max; + max = MAX(idiag[0],idiag[1]); + max = MAX(max,idiag[2]); + + if (idiag[0] < EPSILON*max) idiag[0] = 0.0; + if (idiag[1] < EPSILON*max) idiag[1] = 0.0; + if (idiag[2] < EPSILON*max) idiag[2] = 0.0; + + // calculate omega using diagonalized inertia matrix + + MathExtra::angmom_to_omega(&angmomall[m][0],ex,ey,ez,idiag,&omega[m][0]); + } } } diff --git a/src/group.cpp b/src/group.cpp index 5b572b2070..da0e94fc11 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -29,6 +29,7 @@ #include "input.h" #include "variable.h" #include "dump.h" +#include "math_extra.h" #include "memory.h" #include "error.h" @@ -37,6 +38,7 @@ using namespace LAMMPS_NS; #define MAX_GROUP 32 +#define EPSILON 1.0e-6 enum{TYPE,MOLECULE,ID}; enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN}; @@ -1662,42 +1664,47 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) } /* ---------------------------------------------------------------------- - compute angular velocity omega from L = Iw, inverting I to solve for w + compute angular velocity omega from L and I really not a group/region operation, but L,I were computed for a group/region + diagonalize I instead of inverting it, to allow for a singular matrix ------------------------------------------------------------------------- */ void Group::omega(double *angmom, double inertia[3][3], double *w) { - double inverse[3][3]; + double idiag[3],ex[3],ey[3],ez[3],cross[3]; + double evectors[3][3]; + + int ierror = MathExtra::jacobi(inertia,idiag,evectors); + if (ierror) error->all(FLERR, + "Insufficient Jacobi rotations for group::omega"); - inverse[0][0] = inertia[1][1]*inertia[2][2] - inertia[1][2]*inertia[2][1]; - inverse[0][1] = -(inertia[0][1]*inertia[2][2] - inertia[0][2]*inertia[2][1]); - inverse[0][2] = inertia[0][1]*inertia[1][2] - inertia[0][2]*inertia[1][1]; - - inverse[1][0] = -(inertia[1][0]*inertia[2][2] - inertia[1][2]*inertia[2][0]); - inverse[1][1] = inertia[0][0]*inertia[2][2] - inertia[0][2]*inertia[2][0]; - inverse[1][2] = -(inertia[0][0]*inertia[1][2] - inertia[0][2]*inertia[1][0]); - - inverse[2][0] = inertia[1][0]*inertia[2][1] - inertia[1][1]*inertia[2][0]; - inverse[2][1] = -(inertia[0][0]*inertia[2][1] - inertia[0][1]*inertia[2][0]); - inverse[2][2] = inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0]; - - double determinant = inertia[0][0]*inertia[1][1]*inertia[2][2] + - inertia[0][1]*inertia[1][2]*inertia[2][0] + - inertia[0][2]*inertia[1][0]*inertia[2][1] - - inertia[0][0]*inertia[1][2]*inertia[2][1] - - inertia[0][1]*inertia[1][0]*inertia[2][2] - - inertia[2][0]*inertia[1][1]*inertia[0][2]; - - if (determinant > 0.0) - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) - inverse[i][j] /= determinant; - - w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] + - inverse[0][2]*angmom[2]; - w[1] = inverse[1][0]*angmom[0] + inverse[1][1]*angmom[1] + - inverse[1][2]*angmom[2]; - w[2] = inverse[2][0]*angmom[0] + inverse[2][1]*angmom[1] + - inverse[2][2]*angmom[2]; + ex[0] = evectors[0][0]; + ex[1] = evectors[1][0]; + ex[2] = evectors[2][0]; + ey[0] = evectors[0][1]; + ey[1] = evectors[1][1]; + ey[2] = evectors[2][1]; + ez[0] = evectors[0][2]; + ez[1] = evectors[1][2]; + ez[2] = evectors[2][2]; + + // enforce 3 evectors as a right-handed coordinate system + // flip 3rd vector if needed + + MathExtra::cross3(ex,ey,cross); + if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez); + + // if any principal moment < scaled EPSILON, set to 0.0 + + double max; + max = MAX(idiag[0],idiag[1]); + max = MAX(max,idiag[2]); + + if (idiag[0] < EPSILON*max) idiag[0] = 0.0; + if (idiag[1] < EPSILON*max) idiag[1] = 0.0; + if (idiag[2] < EPSILON*max) idiag[2] = 0.0; + + // calculate omega using diagonalized inertia matrix + + MathExtra::angmom_to_omega(angmom,ex,ey,ez,idiag,w); } diff --git a/src/version.h b/src/version.h index 9725014aa4..006495e9f1 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "27 Aug 2016" +#define LAMMPS_VERSION "7 Sep 2016"