diff --git a/src/KOKKOS/fix_cmap_kokkos.cpp b/src/KOKKOS/fix_cmap_kokkos.cpp index e597c34334..ece0ece1db 100644 --- a/src/KOKKOS/fix_cmap_kokkos.cpp +++ b/src/KOKKOS/fix_cmap_kokkos.cpp @@ -56,6 +56,7 @@ FixCMAPKokkos::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) : datamask_modify = F_MASK; // allocate memory for CMAP data + memoryKK->create_kokkos(k_g_axis,g_axis,CMAPDIM,"cmap:g_axis"); memoryKK->create_kokkos(k_cmapgrid,cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:grid"); memoryKK->create_kokkos(k_d1cmapgrid,d1cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d1grid"); @@ -69,6 +70,7 @@ FixCMAPKokkos::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) : d_d12cmapgrid = k_d12cmapgrid.template view(); // read and setup CMAP data + read_grid_map(arg[3]); int i = 0; @@ -88,6 +90,7 @@ FixCMAPKokkos::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) : for( int i=0 ; i::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) : k_d1cmapgrid.template sync(); k_d2cmapgrid.template sync(); k_d12cmapgrid.template sync(); - } /* ---------------------------------------------------------------------- */ @@ -136,7 +138,6 @@ FixCMAPKokkos::~FixCMAPKokkos() memoryKK->destroy_kokkos(k_crossterm_atom5,crossterm_atom5); memoryKK->destroy_kokkos(d_crosstermlist); - } /* ---------------------------------------------------------------------- */ @@ -148,6 +149,7 @@ void FixCMAPKokkos::init() error->all(FLERR,"Cannot yet use respa with Kokkos"); // on KOKKOS, allocate enough for all crossterms on each GPU to avoid grow operation in device code + maxcrossterm = ncmap; memoryKK->create_kokkos(d_crosstermlist,maxcrossterm,CMAPMAX,"cmap:crosstermlist"); } @@ -159,7 +161,6 @@ void FixCMAPKokkos::init() template void FixCMAPKokkos::pre_neighbor() { - atomKK->sync(execution_space,X_MASK); d_x = atomKK->k_x.view(); int nlocal = atomKK->nlocal; @@ -179,14 +180,12 @@ void FixCMAPKokkos::pre_neighbor() copymode = 1; Kokkos::parallel_scan(Kokkos::RangePolicy(0,nlocal),*this,ncrosstermlist); copymode = 0; - } template KOKKOS_INLINE_FUNCTION void FixCMAPKokkos::operator()(TagFixCmapPreNeighbor, const int i, int &l_ncrosstermlist, const bool is_final ) const { - for( int m = 0; m < d_num_crossterm(i); m++) { int atom1 = AtomKokkos::map_kokkos(d_crossterm_atom1(i,m),map_style,k_map_array,k_map_hash); @@ -297,6 +296,7 @@ void FixCMAPKokkos::operator()(TagFixCmapPostForce, const int n, dou double vb45z = d_x(i4,2) - d_x(i5,2); // calculate normal vectors for planes that define the dihedral angles + double a1x = vb12y*vb23z - vb12z*vb23y; double a1y = vb12z*vb23x - vb12x*vb23z; double a1z = vb12x*vb23y - vb12y*vb23x; @@ -325,6 +325,7 @@ void FixCMAPKokkos::operator()(TagFixCmapPostForce, const int n, dou if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) return; // vectors needed to calculate the cross-term dihedral angles + double dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z; double dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z; double dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z; @@ -388,8 +389,8 @@ void FixCMAPKokkos::operator()(TagFixCmapPostForce, const int n, dou bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs,E,dEdPhi,dEdPsi); // sum up cmap energy contributions - // needed for compute_scalar() + double engfraction = 0.2 * E; if (i1 < nlocal) ecmapKK += engfraction; if (i2 < nlocal) ecmapKK += engfraction; @@ -479,6 +480,7 @@ void FixCMAPKokkos::grow_arrays(int nmax) k_crossterm_atom5.template sync(); // force reallocation on host + k_num_crossterm.template modify(); k_crossterm_type.template modify(); k_crossterm_atom1.template modify(); @@ -877,8 +879,6 @@ void FixCMAPKokkos::bc_interpol(double x1, double x2, int low1, int dEdPsi *= (180.0/MY_PI/CMAPDX); } - - /* ---------------------------------------------------------------------- return local index of atom J or any of its images that is closest to atom I if J is not a valid index like -1, just return it @@ -917,7 +917,6 @@ int FixCMAPKokkos::closest_image(const int i, int j) const return closest; } - namespace LAMMPS_NS { template class FixCMAPKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_limit_kokkos.cpp b/src/KOKKOS/fix_nve_limit_kokkos.cpp index 942ee41f3a..de77427e49 100644 --- a/src/KOKKOS/fix_nve_limit_kokkos.cpp +++ b/src/KOKKOS/fix_nve_limit_kokkos.cpp @@ -120,7 +120,6 @@ void FixNVELimitKokkos::initial_integrate(int /*vflag*/) ncount += d_ncount; atomKK->modified(execution_space, X_MASK | V_MASK ); - } /* ---------------------------------------------------------------------- */ @@ -190,7 +189,6 @@ void FixNVELimitKokkos::final_integrate() ncount += d_ncount; atomKK->modified(execution_space, V_MASK ); - } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/fix_recenter_kokkos.cpp b/src/KOKKOS/fix_recenter_kokkos.cpp index c3a840ff10..607f5ce8d9 100644 --- a/src/KOKKOS/fix_recenter_kokkos.cpp +++ b/src/KOKKOS/fix_recenter_kokkos.cpp @@ -50,7 +50,6 @@ FixRecenterKokkos::FixRecenterKokkos(LAMMPS *lmp, int narg, char **a template void FixRecenterKokkos::initial_integrate(int /*vflag*/) { - atomKK->sync(execution_space,datamask_read); int nlocal = atomKK->nlocal; if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; @@ -121,7 +120,6 @@ void FixRecenterKokkos::initial_integrate(int /*vflag*/) atomKK->modified(execution_space,datamask_modify); } - namespace LAMMPS_NS { template class FixRecenterKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_wall_region_kokkos.cpp b/src/KOKKOS/fix_wall_region_kokkos.cpp index 2eea884472..ab6f7186a1 100644 --- a/src/KOKKOS/fix_wall_region_kokkos.cpp +++ b/src/KOKKOS/fix_wall_region_kokkos.cpp @@ -91,8 +91,8 @@ void FixWallRegionKokkos::post_force(int vflag) // initilize ewall after region->prematch(), // so a dynamic region can access last timestep values - // energy intialize. - // eflag is used to track whether wall energies have been communicated. + // energy intialize + // eflag is used to track whether wall energies have been communicated eflag = 0; double result[10]; @@ -330,7 +330,6 @@ template KOKKOS_INLINE_FUNCTION void FixWallRegionKokkos::v_tally(value_type result, int i, double *v) const { - if (vflag_global) { result[4] += v[0]; result[5] += v[1]; @@ -348,7 +347,6 @@ void FixWallRegionKokkos::v_tally(value_type result, int i, double * Kokkos::atomic_add(&(d_vatom(i,4)),v[4]); Kokkos::atomic_add(&(d_vatom(i,5)),v[5]); } - } namespace LAMMPS_NS { diff --git a/src/KOKKOS/group_kokkos.cpp b/src/KOKKOS/group_kokkos.cpp index 01cf15f6c5..fb115eca0e 100644 --- a/src/KOKKOS/group_kokkos.cpp +++ b/src/KOKKOS/group_kokkos.cpp @@ -71,7 +71,6 @@ double GroupKokkos::mass(int igroup) return all; } - /* ---------------------------------------------------------------------- compute the center-of-mass coords of group of atoms masstotal = total mass diff --git a/src/KOKKOS/group_kokkos.h b/src/KOKKOS/group_kokkos.h index c8573b0d74..f62f192b84 100644 --- a/src/KOKKOS/group_kokkos.h +++ b/src/KOKKOS/group_kokkos.h @@ -23,7 +23,7 @@ template class GroupKokkos : public Group { public: GroupKokkos(class LAMMPS *); - double mass(int); // total mass of atoms in group + double mass(int); // total mass of atoms in group void xcm(int, double, double *); // center-of-mass coords of group }; diff --git a/src/KOKKOS/math_special_kokkos.h b/src/KOKKOS/math_special_kokkos.h index 12e04db1c0..1cc35e1969 100644 --- a/src/KOKKOS/math_special_kokkos.h +++ b/src/KOKKOS/math_special_kokkos.h @@ -271,7 +271,6 @@ namespace MathSpecialKokkos { return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; } - } // namespace MathSpecialKokkos } // namespace LAMMPS_NS diff --git a/src/KOKKOS/region_block_kokkos.h b/src/KOKKOS/region_block_kokkos.h index d0dcc12b88..052a6a4bcf 100644 --- a/src/KOKKOS/region_block_kokkos.h +++ b/src/KOKKOS/region_block_kokkos.h @@ -61,43 +61,44 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { KOKKOS_INLINE_FUNCTION int surface_kokkos(double x, double y, double z, double cutoff) -{ - int ncontact; - double xs, ys, zs; - double xnear[3], xorig[3]; + { + int ncontact; + double xs, ys, zs; + double xnear[3], xorig[3]; - if (dynamic) { - xorig[0] = x; xorig[1] = y; xorig[2] = z; - inverse_transform(x, y, z); - } - - xnear[0] = x; xnear[1] = y; xnear[2] = z; - - if (!openflag) { - if (interior) - ncontact = surface_interior_kokkos(xnear, cutoff); - else - ncontact = surface_exterior_kokkos(xnear, cutoff); - } else { - // one of surface_int/ext() will return 0 - // so no need to worry about offset of contact indices - ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); - } - - if (rotateflag && ncontact) { - for (int i = 0; i < ncontact; i++) { - xs = xnear[0] - d_contact[i].delx; - ys = xnear[1] - d_contact[i].dely; - zs = xnear[2] - d_contact[i].delz; - forward_transform(xs, ys, zs); - d_contact[i].delx = xorig[0] - xs; - d_contact[i].dely = xorig[1] - ys; - d_contact[i].delz = xorig[2] - zs; + if (dynamic) { + xorig[0] = x; xorig[1] = y; xorig[2] = z; + inverse_transform(x, y, z); } + + xnear[0] = x; xnear[1] = y; xnear[2] = z; + + if (!openflag) { + if (interior) + ncontact = surface_interior_kokkos(xnear, cutoff); + else + ncontact = surface_exterior_kokkos(xnear, cutoff); + } else { + // one of surface_int/ext() will return 0 + // so no need to worry about offset of contact indices + ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); + } + + if (rotateflag && ncontact) { + for (int i = 0; i < ncontact; i++) { + xs = xnear[0] - d_contact[i].delx; + ys = xnear[1] - d_contact[i].dely; + zs = xnear[2] - d_contact[i].delz; + forward_transform(xs, ys, zs); + d_contact[i].delx = xorig[0] - xs; + d_contact[i].dely = xorig[1] - ys; + d_contact[i].delz = xorig[2] - zs; + } + } + + return ncontact; } - return ncontact; -} Kokkos::View d_contact; private: @@ -106,317 +107,316 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { typename AT::t_x_array_randomread d_x; typename AT::t_int_1d_randomread d_mask; -KOKKOS_INLINE_FUNCTION -int surface_interior_kokkos(double *x, double cutoff) -{ - double delta; + KOKKOS_INLINE_FUNCTION + int surface_interior_kokkos(double *x, double cutoff) + { + double delta; - // x is exterior to block + // x is exterior to block - if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0; + if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0; - // x is interior to block or on its surface + // x is interior to block or on its surface - int n = 0; + int n = 0; - delta = x[0] - xlo; - if (delta < cutoff && !open_faces[0]) { - d_contact[n].r = delta; - d_contact[n].delx = delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 0; - n++; - } - delta = xhi - x[0]; - if (delta < cutoff && !open_faces[1]) { - d_contact[n].r = delta; - d_contact[n].delx = -delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 1; - n++; + delta = x[0] - xlo; + if (delta < cutoff && !open_faces[0]) { + d_contact[n].r = delta; + d_contact[n].delx = delta; + d_contact[n].dely = d_contact[n].delz = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 0; + n++; + } + delta = xhi - x[0]; + if (delta < cutoff && !open_faces[1]) { + d_contact[n].r = delta; + d_contact[n].delx = -delta; + d_contact[n].dely = d_contact[n].delz = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 1; + n++; + } + + delta = x[1] - ylo; + if (delta < cutoff && !open_faces[2]) { + d_contact[n].r = delta; + d_contact[n].dely = delta; + d_contact[n].delx = d_contact[n].delz = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 2; + n++; + } + delta = yhi - x[1]; + if (delta < cutoff && !open_faces[3]) { + d_contact[n].r = delta; + d_contact[n].dely = -delta; + d_contact[n].delx = d_contact[n].delz = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 3; + n++; + } + + delta = x[2] - zlo; + if (delta < cutoff && !open_faces[4]) { + d_contact[n].r = delta; + d_contact[n].delz = delta; + d_contact[n].delx = d_contact[n].dely = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 4; + n++; + } + delta = zhi - x[2]; + if (delta < cutoff && !open_faces[5]) { + d_contact[n].r = delta; + d_contact[n].delz = -delta; + d_contact[n].delx = d_contact[n].dely = 0.0; + d_contact[n].radius = 0; + d_contact[n].iwall = 5; + n++; + } + + return n; } - delta = x[1] - ylo; - if (delta < cutoff && !open_faces[2]) { - d_contact[n].r = delta; - d_contact[n].dely = delta; - d_contact[n].delx = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 2; - n++; - } - delta = yhi - x[1]; - if (delta < cutoff && !open_faces[3]) { - d_contact[n].r = delta; - d_contact[n].dely = -delta; - d_contact[n].delx = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 3; - n++; - } + KOKKOS_INLINE_FUNCTION + int surface_exterior_kokkos(double *x, double cutoff) + { + double xp, yp, zp; + double xc, yc, zc, dist, mindist; - delta = x[2] - zlo; - if (delta < cutoff && !open_faces[4]) { - d_contact[n].r = delta; - d_contact[n].delz = delta; - d_contact[n].delx = d_contact[n].dely = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 4; - n++; - } - delta = zhi - x[2]; - if (delta < cutoff && !open_faces[5]) { - d_contact[n].r = delta; - d_contact[n].delz = -delta; - d_contact[n].delx = d_contact[n].dely = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 5; - n++; - } + // x is far enough from block that there is no contact + // x is interior to block - return n; -} + if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff || + x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff) + return 0; + if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0; -KOKKOS_INLINE_FUNCTION -int surface_exterior_kokkos(double *x, double cutoff) -{ - double xp, yp, zp; - double xc, yc, zc, dist, mindist; + // x is exterior to block or on its surface + // xp,yp,zp = point on surface of block that x is closest to + // could be edge or corner pt of block + // do not add contact point if r >= cutoff - // x is far enough from block that there is no contact - // x is interior to block - - if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff || - x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff) - return 0; - if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0; - - // x is exterior to block or on its surface - // xp,yp,zp = point on surface of block that x is closest to - // could be edge or corner pt of block - // do not add contact point if r >= cutoff - - if (!openflag) { - if (x[0] < xlo) - xp = xlo; - else if (x[0] > xhi) - xp = xhi; - else - xp = x[0]; - if (x[1] < ylo) - yp = ylo; - else if (x[1] > yhi) - yp = yhi; - else - yp = x[1]; - if (x[2] < zlo) - zp = zlo; - else if (x[2] > zhi) - zp = zhi; - else - zp = x[2]; - } else { - mindist = MAXDOUBLEINT; - for (int i = 0; i < 6; i++) { - if (open_faces[i]) continue; - dist = find_closest_point(i, x, xc, yc, zc); - if (dist < mindist) { - xp = xc; - yp = yc; - zp = zc; - mindist = dist; + if (!openflag) { + if (x[0] < xlo) + xp = xlo; + else if (x[0] > xhi) + xp = xhi; + else + xp = x[0]; + if (x[1] < ylo) + yp = ylo; + else if (x[1] > yhi) + yp = yhi; + else + yp = x[1]; + if (x[2] < zlo) + zp = zlo; + else if (x[2] > zhi) + zp = zhi; + else + zp = x[2]; + } else { + mindist = MAXDOUBLEINT; + for (int i = 0; i < 6; i++) { + if (open_faces[i]) continue; + dist = find_closest_point(i, x, xc, yc, zc); + if (dist < mindist) { + xp = xc; + yp = yc; + zp = zc; + mindist = dist; + } } } + + add_contact(0, x, xp, yp, zp); + d_contact[0].iwall = 0; + if (d_contact[0].r < cutoff) return 1; + return 0; } - add_contact(0, x, xp, yp, zp); - d_contact[0].iwall = 0; - if (d_contact[0].r < cutoff) return 1; - return 0; -} - -KOKKOS_INLINE_FUNCTION -void add_contact(int n, double *x, double xp, double yp, double zp) -{ - double delx = x[0] - xp; - double dely = x[1] - yp; - double delz = x[2] - zp; - d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); - d_contact[n].radius = 0; - d_contact[n].delx = delx; - d_contact[n].dely = dely; - d_contact[n].delz = delz; -} - -KOKKOS_INLINE_FUNCTION -int k_inside(double x, double y, double z) const -{ - if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) - return 1; - return 0; -} - -KOKKOS_INLINE_FUNCTION -void forward_transform(double &x, double &y, double &z) const -{ - if (rotateflag) rotate(x, y, z, theta); - if (moveflag) { - x += dx; - y += dy; - z += dz; + KOKKOS_INLINE_FUNCTION + void add_contact(int n, double *x, double xp, double yp, double zp) + { + double delx = x[0] - xp; + double dely = x[1] - yp; + double delz = x[2] - zp; + d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); + d_contact[n].radius = 0; + d_contact[n].delx = delx; + d_contact[n].dely = dely; + d_contact[n].delz = delz; } -} -KOKKOS_INLINE_FUNCTION -void inverse_transform(double &x, double &y, double &z) const -{ - if (moveflag) { - x -= dx; - y -= dy; - z -= dz; - } - if (rotateflag) rotate(x,y,z,-theta); -} - -KOKKOS_INLINE_FUNCTION -void rotate(double &x, double &y, double &z, double angle) const -{ - double a[3],b[3],c[3],d[3],disp[3]; - - double sine = sin(angle); - double cosine = cos(angle); - d[0] = x - point[0]; - d[1] = y - point[1]; - d[2] = z - point[2]; - double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; - c[0] = x0dotr * runit[0]; - c[1] = x0dotr * runit[1]; - c[2] = x0dotr * runit[2]; - a[0] = d[0] - c[0]; - a[1] = d[1] - c[1]; - a[2] = d[2] - c[2]; - b[0] = runit[1]*a[2] - runit[2]*a[1]; - b[1] = runit[2]*a[0] - runit[0]*a[2]; - b[2] = runit[0]*a[1] - runit[1]*a[0]; - disp[0] = a[0]*cosine + b[0]*sine; - disp[1] = a[1]*cosine + b[1]*sine; - disp[2] = a[2]*cosine + b[2]*sine; - x = point[0] + c[0] + disp[0]; - y = point[1] + c[1] + disp[1]; - z = point[2] + c[2] + disp[2]; -} - -KOKKOS_INLINE_FUNCTION -void point_on_line_segment(double *a, double *b, double *c, double *d) -{ - double ba[3], ca[3]; - - sub3(b, a, ba); - sub3(c, a, ca); - double t = dot3(ca, ba) / dot3(ba, ba); - if (t <= 0.0) { - d[0] = a[0]; - d[1] = a[1]; - d[2] = a[2]; - } else if (t >= 1.0) { - d[0] = b[0]; - d[1] = b[1]; - d[2] = b[2]; - } else { - d[0] = a[0] + t * ba[0]; - d[1] = a[1] + t * ba[1]; - d[2] = a[2] + t * ba[2]; - } -} - -KOKKOS_INLINE_FUNCTION -double inside_face(double *xproj, int iface) -{ - if (iface < 2) { - if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; - } else if (iface < 4) { - if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo))) + KOKKOS_INLINE_FUNCTION + int k_inside(double x, double y, double z) const + { + if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) return 1; - } else { - if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1; + return 0; } - return 0; -} - -KOKKOS_INLINE_FUNCTION -double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) -{ - double dot, d2, d2min; - double xr[3], xproj[3], p[3]; - - xr[0] = x[0] - corners[i][0][0]; - xr[1] = x[1] - corners[i][0][1]; - xr[2] = x[2] - corners[i][0][2]; - dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2]; - xproj[0] = xr[0] - dot * face[i][0]; - xproj[1] = xr[1] - dot * face[i][1]; - xproj[2] = xr[2] - dot * face[i][2]; - - d2min = MAXDOUBLEINT; - - // check if point projects inside of face - - if (inside_face(xproj, i)) { - d2 = d2min = dot * dot; - xc = xproj[0] + corners[i][0][0]; - yc = xproj[1] + corners[i][0][1]; - zc = xproj[2] + corners[i][0][2]; - - // check each edge - - } else { - point_on_line_segment(corners[i][0], corners[i][1], x, p); - d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + - (p[2] - x[2]) * (p[2] - x[2]); - if (d2 < d2min) { - d2min = d2; - xc = p[0]; - yc = p[1]; - zc = p[2]; - } - - point_on_line_segment(corners[i][1], corners[i][2], x, p); - d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + - (p[2] - x[2]) * (p[2] - x[2]); - if (d2 < d2min) { - d2min = d2; - xc = p[0]; - yc = p[1]; - zc = p[2]; - } - - point_on_line_segment(corners[i][2], corners[i][3], x, p); - d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + - (p[2] - x[2]) * (p[2] - x[2]); - if (d2 < d2min) { - d2min = d2; - xc = p[0]; - yc = p[1]; - zc = p[2]; - } - - point_on_line_segment(corners[i][3], corners[i][0], x, p); - d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + - (p[2] - x[2]) * (p[2] - x[2]); - if (d2 < d2min) { - d2min = d2; - xc = p[0]; - yc = p[1]; - zc = p[2]; + KOKKOS_INLINE_FUNCTION + void forward_transform(double &x, double &y, double &z) const + { + if (rotateflag) rotate(x, y, z, theta); + if (moveflag) { + x += dx; + y += dy; + z += dz; } } - return d2min; -} + KOKKOS_INLINE_FUNCTION + void inverse_transform(double &x, double &y, double &z) const + { + if (moveflag) { + x -= dx; + y -= dy; + z -= dz; + } + if (rotateflag) rotate(x,y,z,-theta); + } + KOKKOS_INLINE_FUNCTION + void rotate(double &x, double &y, double &z, double angle) const + { + double a[3],b[3],c[3],d[3],disp[3]; + + double sine = sin(angle); + double cosine = cos(angle); + d[0] = x - point[0]; + d[1] = y - point[1]; + d[2] = z - point[2]; + double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2]; + c[0] = x0dotr * runit[0]; + c[1] = x0dotr * runit[1]; + c[2] = x0dotr * runit[2]; + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1]*a[2] - runit[2]*a[1]; + b[1] = runit[2]*a[0] - runit[0]*a[2]; + b[2] = runit[0]*a[1] - runit[1]*a[0]; + disp[0] = a[0]*cosine + b[0]*sine; + disp[1] = a[1]*cosine + b[1]*sine; + disp[2] = a[2]*cosine + b[2]*sine; + x = point[0] + c[0] + disp[0]; + y = point[1] + c[1] + disp[1]; + z = point[2] + c[2] + disp[2]; + } + + KOKKOS_INLINE_FUNCTION + void point_on_line_segment(double *a, double *b, double *c, double *d) + { + double ba[3], ca[3]; + + sub3(b, a, ba); + sub3(c, a, ca); + double t = dot3(ca, ba) / dot3(ba, ba); + if (t <= 0.0) { + d[0] = a[0]; + d[1] = a[1]; + d[2] = a[2]; + } else if (t >= 1.0) { + d[0] = b[0]; + d[1] = b[1]; + d[2] = b[2]; + } else { + d[0] = a[0] + t * ba[0]; + d[1] = a[1] + t * ba[1]; + d[2] = a[2] + t * ba[2]; + } + } + + KOKKOS_INLINE_FUNCTION + double inside_face(double *xproj, int iface) + { + if (iface < 2) { + if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; + } else if (iface < 4) { + if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo))) + return 1; + } else { + if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1; + } + + return 0; + } + + KOKKOS_INLINE_FUNCTION + double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) + { + double dot, d2, d2min; + double xr[3], xproj[3], p[3]; + + xr[0] = x[0] - corners[i][0][0]; + xr[1] = x[1] - corners[i][0][1]; + xr[2] = x[2] - corners[i][0][2]; + dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2]; + xproj[0] = xr[0] - dot * face[i][0]; + xproj[1] = xr[1] - dot * face[i][1]; + xproj[2] = xr[2] - dot * face[i][2]; + + d2min = MAXDOUBLEINT; + + // check if point projects inside of face + + if (inside_face(xproj, i)) { + d2 = d2min = dot * dot; + xc = xproj[0] + corners[i][0][0]; + yc = xproj[1] + corners[i][0][1]; + zc = xproj[2] + corners[i][0][2]; + + // check each edge + + } else { + point_on_line_segment(corners[i][0], corners[i][1], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); + if (d2 < d2min) { + d2min = d2; + xc = p[0]; + yc = p[1]; + zc = p[2]; + } + + point_on_line_segment(corners[i][1], corners[i][2], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); + if (d2 < d2min) { + d2min = d2; + xc = p[0]; + yc = p[1]; + zc = p[2]; + } + + point_on_line_segment(corners[i][2], corners[i][3], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); + if (d2 < d2min) { + d2min = d2; + xc = p[0]; + yc = p[1]; + zc = p[2]; + } + + point_on_line_segment(corners[i][3], corners[i][0], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); + if (d2 < d2min) { + d2min = d2; + xc = p[0]; + yc = p[1]; + zc = p[2]; + } + } + + return d2min; + } }; diff --git a/src/KOKKOS/region_sphere_kokkos.cpp b/src/KOKKOS/region_sphere_kokkos.cpp index b9a305d9fb..07275ee69e 100644 --- a/src/KOKKOS/region_sphere_kokkos.cpp +++ b/src/KOKKOS/region_sphere_kokkos.cpp @@ -75,7 +75,6 @@ void RegSphereKokkos::operator()(TagRegSphereMatchAll, const int &i) } } - /* ---------------------------------------------------------------------- */ namespace LAMMPS_NS { diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index d9859326c9..02116965b5 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -131,6 +131,7 @@ FixCMAP::~FixCMAP() if (copymode) return; // unregister callbacks to this fix from Atom class + atom->delete_callback(id,Atom::GROW); atom->delete_callback(id,Atom::RESTART); diff --git a/src/region.cpp b/src/region.cpp index 6bcbc4470a..7399b14adb 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -50,6 +50,7 @@ Region::Region(LAMMPS *lmp, int /*narg*/, char **arg) : Region::~Region() { if (copymode) return; + delete[] id; delete[] style; delete[] xstr; diff --git a/src/region_block.cpp b/src/region_block.cpp index 9376016843..36c38f517c 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -262,6 +262,7 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : RegBlock::~RegBlock() { if (copymode) return; + delete[] xlostr; delete[] xhistr; delete[] ylostr; diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index ea6e39d894..ec472c031c 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -102,6 +102,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::~RegSphere() { if (copymode) return; + delete[] xstr; delete[] ystr; delete[] zstr;