diff --git a/src/KOKKOS/compute_temp_kokkos.cpp b/src/KOKKOS/compute_temp_kokkos.cpp new file mode 100755 index 0000000000..bb5e2799c2 --- /dev/null +++ b/src/KOKKOS/compute_temp_kokkos.cpp @@ -0,0 +1,153 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "string.h" +#include "compute_temp_kokkos.h" +#include "atom_kokkos.h" +#include "update.h" +#include "force.h" +#include "domain.h" +#include "comm.h" +#include "group.h" +#include "error.h" +#include "atom_masks.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +template +ComputeTempKokkos::ComputeTempKokkos(LAMMPS *lmp, int narg, char **arg) : + ComputeTemp(lmp, narg, arg) +{ + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + + datamask_read = V_MASK | MASK_MASK | RMASS_MASK | TYPE_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +double ComputeTempKokkos::compute_scalar() +{ + atomKK->sync(execution_space,datamask_read); + + invoked_scalar = update->ntimestep; + + v = atomKK->k_v.view(); + rmass = atomKK->rmass; + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); + int nlocal = atom->nlocal; + + double t = 0.0; + CTEMP t_kk; + + copymode = 1; + if (rmass) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,nlocal),*this,t_kk); + else + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,nlocal),*this,t_kk); + DeviceType::fence(); + copymode = 0; + + t = t_kk.t0; // could make this more efficient + + MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + if (dynamic) dof_compute(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); + scalar *= tfactor; + return scalar; +} + +template +template +KOKKOS_INLINE_FUNCTION +void ComputeTempKokkos::operator()(TagComputeTempScalar, const int &i, CTEMP& t_kk) const { + if (RMASS) { + if (mask[i] & groupbit) + t_kk.t0 += (v(i,0)*v(i,0) + v(i,1)*v(i,1) + v(i,2)*v(i,2)) * rmass[i]; + } else { + if (mask[i] & groupbit) + t_kk.t0 += (v(i,0)*v(i,0) + v(i,1)*v(i,1) + v(i,2)*v(i,2)) * + mass[type[i]]; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeTempKokkos::compute_vector() +{ + atomKK->sync(execution_space,datamask_read); + + int i; + + invoked_vector = update->ntimestep; + + v = atomKK->k_v.view(); + rmass = atomKK->rmass; + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); + int nlocal = atom->nlocal; + + double t[6]; + for (i = 0; i < 6; i++) t[i] = 0.0; + CTEMP t_kk; + + copymode = 1; + if (rmass) + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,nlocal),*this,t_kk); + else + Kokkos::parallel_reduce(Kokkos::RangePolicy >(0,nlocal),*this,t_kk); + DeviceType::fence(); + copymode = 0; + + t[0] = t_kk.t0; + t[1] = t_kk.t1; + t[2] = t_kk.t2; + t[3] = t_kk.t3; + t[4] = t_kk.t4; + t[5] = t_kk.t5; + + MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; +} + +template +template +KOKKOS_INLINE_FUNCTION +void ComputeTempKokkos::operator()(TagComputeTempVector, const int &i, CTEMP& t_kk) const { + if (mask[i] & groupbit) { + F_FLOAT massone = 0.0; + if (RMASS) massone = rmass[i]; + else massone = mass[type[i]]; + t_kk.t0 += massone * v(i,0)*v(i,0); + t_kk.t1 += massone * v(i,1)*v(i,1); + t_kk.t2 += massone * v(i,2)*v(i,2); + t_kk.t3 += massone * v(i,0)*v(i,1); + t_kk.t4 += massone * v(i,0)*v(i,2); + t_kk.t5 += massone * v(i,1)*v(i,2); + } +} + +template class ComputeTempKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class ComputeTempKokkos; +#endif \ No newline at end of file diff --git a/src/KOKKOS/compute_temp_kokkos.h b/src/KOKKOS/compute_temp_kokkos.h new file mode 100755 index 0000000000..9cdb87c8b9 --- /dev/null +++ b/src/KOKKOS/compute_temp_kokkos.h @@ -0,0 +1,107 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(temp/kk,ComputeTempKokkos) +ComputeStyle(temp/kk/device,ComputeTempKokkos) +ComputeStyle(temp/kk/host,ComputeTempKokkos) + +#else + +#ifndef LMP_COMPUTE_TEMP_KOKKOS_H +#define LMP_COMPUTE_TEMP_KOKKOS_H + +#include "compute_temp.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + + struct s_CTEMP { + double t0, t1, t2, t3, t4, t5; + KOKKOS_INLINE_FUNCTION + s_CTEMP() { + t0 = t1 = t2 = t3 = t4 = t5 = 0.0; + } + KOKKOS_INLINE_FUNCTION + s_CTEMP& operator+=(const s_CTEMP &rhs){ + t0 += rhs.t0; + t1 += rhs.t1; + t2 += rhs.t2; + t3 += rhs.t3; + t4 += rhs.t4; + t5 += rhs.t5; + return *this; + } + + KOKKOS_INLINE_FUNCTION + volatile s_CTEMP& operator+=(const volatile s_CTEMP &rhs) volatile { + t0 += rhs.t0; + t1 += rhs.t1; + t2 += rhs.t2; + t3 += rhs.t3; + t4 += rhs.t4; + t5 += rhs.t5; + return *this; + } + }; + typedef s_CTEMP CTEMP; + +template +struct TagComputeTempScalar{}; + +template +struct TagComputeTempVector{}; + +template +class ComputeTempKokkos : public ComputeTemp { + public: + typedef DeviceType device_type; + typedef CTEMP value_type; + typedef ArrayTypes AT; + + ComputeTempKokkos(class LAMMPS *, int, char **); + virtual ~ComputeTempKokkos() {} + double compute_scalar(); + void compute_vector(); + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagComputeTempScalar, const int&, CTEMP&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagComputeTempVector, const int&, CTEMP&) const; + + protected: + typename ArrayTypes::t_v_array_randomread v; + double *rmass; + typename ArrayTypes::t_float_1d_randomread mass; + typename ArrayTypes::t_int_1d_randomread type; + typename ArrayTypes::t_int_1d_randomread mask; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/src/KOKKOS/fix_deform_kokkos.cpp b/src/KOKKOS/fix_deform_kokkos.cpp new file mode 100755 index 0000000000..7da7d301be --- /dev/null +++ b/src/KOKKOS/fix_deform_kokkos.cpp @@ -0,0 +1,374 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Pieter in 't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_deform_kokkos.h" +#include "atom_kokkos.h" +#include "update.h" +#include "comm.h" +#include "irregular.h" +#include "domain_kokkos.h" +#include "lattice.h" +#include "force.h" +#include "modify.h" +#include "math_const.h" +#include "kspace.h" +#include "input.h" +#include "variable.h" +#include "error.h" +#include "atom_masks.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE}; +enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; + +// same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp + +enum{NO_REMAP,X_REMAP,V_REMAP}; + +/* ---------------------------------------------------------------------- */ + +FixDeformKokkos::FixDeformKokkos(LAMMPS *lmp, int narg, char **arg) : FixDeform(lmp, narg, arg) +{ + domainKK = (DomainKokkos *) domain; + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- + box flipped on previous step + reset box tilts for flipped config and create new box in domain + image_flip() adjusts image flags due to box shape change induced by flip + remap() puts atoms outside the new box back into the new box + perform irregular on atoms in lamda coords to migrate atoms to new procs + important that image_flip comes before remap, since remap may change + image flags to new values, making eqs in doc of Domain:image_flip incorrect +------------------------------------------------------------------------- */ + +void FixDeformKokkos::pre_exchange() +{ + if (flip == 0) return; + + domain->yz = set[3].tilt_target = set[3].tilt_flip; + domain->xz = set[4].tilt_target = set[4].tilt_flip; + domain->xy = set[5].tilt_target = set[5].tilt_flip; + domain->set_global_box(); + domain->set_local_box(); + + domainKK->image_flip(flipxy,flipxz,flipyz); + + domainKK->remap_all(); + + domainKK->x2lamda(atom->nlocal); + atomKK->sync(Host,ALL_MASK); + irregular->migrate_atoms(); + atomKK->modified(Host,ALL_MASK); + domainKK->lamda2x(atom->nlocal); + + flip = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixDeformKokkos::end_of_step() +{ + int i; + + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + // wrap variable evaluations with clear/add + + if (varflag) modify->clearstep_compute(); + + // set new box size + // for NONE, target is current box size + // for TRATE, set target directly based on current time, also set h_rate + // for WIGGLE, set target directly based on current time, also set h_rate + // for VARIABLE, set target directly via variable eval, also set h_rate + // for others except VOLUME, target is linear value between start and stop + + for (i = 0; i < 3; i++) { + if (set[i].style == NONE) { + set[i].lo_target = domain->boxlo[i]; + set[i].hi_target = domain->boxhi[i]; + } else if (set[i].style == TRATE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); + h_rate[i] = set[i].rate * domain->h[i]; + h_ratelo[i] = -0.5*h_rate[i]; + } else if (set[i].style == WIGGLE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].lo_target = set[i].lo_start - + 0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); + set[i].hi_target = set[i].hi_start + + 0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); + h_rate[i] = TWOPI/set[i].tperiod * set[i].amplitude * + cos(TWOPI*delt/set[i].tperiod); + h_ratelo[i] = -0.5*h_rate[i]; + } else if (set[i].style == VARIABLE) { + double del = input->variable->compute_equal(set[i].hvar); + set[i].lo_target = set[i].lo_start - 0.5*del; + set[i].hi_target = set[i].hi_start + 0.5*del; + h_rate[i] = input->variable->compute_equal(set[i].hratevar); + h_ratelo[i] = -0.5*h_rate[i]; + } else if (set[i].style != VOLUME) { + set[i].lo_target = set[i].lo_start + + delta*(set[i].lo_stop - set[i].lo_start); + set[i].hi_target = set[i].hi_start + + delta*(set[i].hi_stop - set[i].hi_start); + } + } + + // set new box size for VOLUME dims that are linked to other dims + // NOTE: still need to set h_rate for these dims + + for (int i = 0; i < 3; i++) { + if (set[i].style != VOLUME) continue; + + if (set[i].substyle == ONE_FROM_ONE) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); + + } else if (set[i].substyle == ONE_FROM_TWO) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].dynamic2].hi_target - + set[set[i].dynamic2].lo_target)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].dynamic2].hi_target - + set[set[i].dynamic2].lo_target)); + + } else if (set[i].substyle == TWO_FROM_ONE) { + set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - + 0.5*sqrt(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start - + set[set[i].fixed].lo_start) * + (set[i].hi_start - set[i].lo_start)); + set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + + 0.5*sqrt(set[i].vol_start / + (set[set[i].dynamic1].hi_target - + set[set[i].dynamic1].lo_target) / + (set[set[i].fixed].hi_start - + set[set[i].fixed].lo_start) * + (set[i].hi_start - set[i].lo_start)); + } + } + + // for triclinic, set new box shape + // for NONE, target is current tilt + // for TRATE, set target directly based on current time. also set h_rate + // for WIGGLE, set target directly based on current time. also set h_rate + // for VARIABLE, set target directly via variable eval. also set h_rate + // for other styles, target is linear value between start and stop values + + if (triclinic) { + double *h = domain->h; + + for (i = 3; i < 6; i++) { + if (set[i].style == NONE) { + if (i == 5) set[i].tilt_target = domain->xy; + else if (i == 4) set[i].tilt_target = domain->xz; + else if (i == 3) set[i].tilt_target = domain->yz; + } else if (set[i].style == TRATE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].tilt_target = set[i].tilt_start * exp(set[i].rate*delt); + h_rate[i] = set[i].rate * domain->h[i]; + } else if (set[i].style == WIGGLE) { + double delt = (update->ntimestep - update->beginstep) * update->dt; + set[i].tilt_target = set[i].tilt_start + + set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); + h_rate[i] = TWOPI/set[i].tperiod * set[i].amplitude * + cos(TWOPI*delt/set[i].tperiod); + } else if (set[i].style == VARIABLE) { + double delta_tilt = input->variable->compute_equal(set[i].hvar); + set[i].tilt_target = set[i].tilt_start + delta_tilt; + h_rate[i] = input->variable->compute_equal(set[i].hratevar); + } else { + set[i].tilt_target = set[i].tilt_start + + delta*(set[i].tilt_stop - set[i].tilt_start); + } + + // tilt_target can be large positive or large negative value + // add/subtract box lengths until tilt_target is closest to current value + + int idenom; + if (i == 5) idenom = 0; + else if (i == 4) idenom = 0; + else if (i == 3) idenom = 1; + double denom = set[idenom].hi_target - set[idenom].lo_target; + + double current = h[i]/h[idenom]; + + while (set[i].tilt_target/denom - current > 0.0) + set[i].tilt_target -= denom; + while (set[i].tilt_target/denom - current < 0.0) + set[i].tilt_target += denom; + if (fabs(set[i].tilt_target/denom - 1.0 - current) < + fabs(set[i].tilt_target/denom - current)) + set[i].tilt_target -= denom; + } + } + + if (varflag) modify->addstep_compute(update->ntimestep + nevery); + + // if any tilt ratios exceed 0.5, set flip = 1 and compute new tilt values + // do not flip in x or y if non-periodic (can tilt but not flip) + // this is b/c the box length would be changed (dramatically) by flip + // if yz tilt exceeded, adjust C vector by one B vector + // if xz tilt exceeded, adjust C vector by one A vector + // if xy tilt exceeded, adjust B vector by one A vector + // check yz first since it may change xz, then xz check comes after + // flip is performed on next timestep, before reneighboring in pre-exchange() + + if (triclinic && flipflag) { + double xprd = set[0].hi_target - set[0].lo_target; + double yprd = set[1].hi_target - set[1].lo_target; + double xprdinv = 1.0 / xprd; + double yprdinv = 1.0 / yprd; + if (set[3].tilt_target*yprdinv < -0.5 || + set[3].tilt_target*yprdinv > 0.5 || + set[4].tilt_target*xprdinv < -0.5 || + set[4].tilt_target*xprdinv > 0.5 || + set[5].tilt_target*xprdinv < -0.5 || + set[5].tilt_target*xprdinv > 0.5) { + set[3].tilt_flip = set[3].tilt_target; + set[4].tilt_flip = set[4].tilt_target; + set[5].tilt_flip = set[5].tilt_target; + + flipxy = flipxz = flipyz = 0; + + if (domain->yperiodic) { + if (set[3].tilt_flip*yprdinv < -0.5) { + set[3].tilt_flip += yprd; + set[4].tilt_flip += set[5].tilt_flip; + flipyz = 1; + } else if (set[3].tilt_flip*yprdinv > 0.5) { + set[3].tilt_flip -= yprd; + set[4].tilt_flip -= set[5].tilt_flip; + flipyz = -1; + } + } + if (domain->xperiodic) { + if (set[4].tilt_flip*xprdinv < -0.5) { + set[4].tilt_flip += xprd; + flipxz = 1; + } + if (set[4].tilt_flip*xprdinv > 0.5) { + set[4].tilt_flip -= xprd; + flipxz = -1; + } + if (set[5].tilt_flip*xprdinv < -0.5) { + set[5].tilt_flip += xprd; + flipxy = 1; + } + if (set[5].tilt_flip*xprdinv > 0.5) { + set[5].tilt_flip -= xprd; + flipxy = -1; + } + } + + flip = 0; + if (flipxy || flipxz || flipyz) flip = 1; + if (flip) next_reneighbor = update->ntimestep + 1; + } + } + + // convert atoms and rigid bodies to lamda coords + + if (remapflag == X_REMAP) { + int nlocal = atom->nlocal; + + domainKK->x2lamda(nlocal); + //for (i = 0; i < nlocal; i++) + // if (mask[i] & groupbit) + // domain->x2lamda(x[i],x[i]); + + if (nrigid) + error->all(FLERR,"Cannot (yet) use rigid bodies with fix deform and Kokkos"); + //for (i = 0; i < nrigid; i++) + // modify->fix[rfix[i]]->deform(0); + } + + // reset global and local box to new size/shape + // only if deform fix is controlling the dimension + + if (set[0].style) { + domain->boxlo[0] = set[0].lo_target; + domain->boxhi[0] = set[0].hi_target; + } + if (set[1].style) { + domain->boxlo[1] = set[1].lo_target; + domain->boxhi[1] = set[1].hi_target; + } + if (set[2].style) { + domain->boxlo[2] = set[2].lo_target; + domain->boxhi[2] = set[2].hi_target; + } + if (triclinic) { + if (set[3].style) domain->yz = set[3].tilt_target; + if (set[4].style) domain->xz = set[4].tilt_target; + if (set[5].style) domain->xy = set[5].tilt_target; + } + + domain->set_global_box(); + domain->set_local_box(); + + // convert atoms and rigid bodies back to box coords + + if (remapflag == X_REMAP) { + int nlocal = atom->nlocal; + + domainKK->lamda2x(nlocal); + //for (i = 0; i < nlocal; i++) + // if (mask[i] & groupbit) + // domain->lamda2x(x[i],x[i]); + + //if (nrigid) + // for (i = 0; i < nrigid; i++) + // modify->fix[rfix[i]]->deform(1); + } + + // redo KSpace coeffs since box has changed + + if (kspace_flag) force->kspace->setup(); +} + diff --git a/src/KOKKOS/fix_deform_kokkos.h b/src/KOKKOS/fix_deform_kokkos.h new file mode 100755 index 0000000000..cd6507a92a --- /dev/null +++ b/src/KOKKOS/fix_deform_kokkos.h @@ -0,0 +1,104 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(deform/kk,FixDeformKokkos) + +#else + +#ifndef LMP_FIX_DEFORM_KOKKOS_H +#define LMP_FIX_DEFORM_KOKKOS_H + +#include "fix_deform.h" + +namespace LAMMPS_NS { + +class FixDeformKokkos : public FixDeform { + public: + + FixDeformKokkos(class LAMMPS *, int, char **); + virtual ~FixDeformKokkos() {} + void pre_exchange(); + void end_of_step(); + + private: + class DomainKokkos *domainKK; + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix deform tilt factors require triclinic box + +Cannot deform the tilt factors of a simulation box unless it +is a triclinic (non-orthogonal) box. + +E: Cannot use fix deform on a shrink-wrapped boundary + +The x, y, z options cannot be applied to shrink-wrapped +dimensions. + +E: Cannot use fix deform tilt on a shrink-wrapped 2nd dim + +This is because the shrink-wrapping will change the value +of the strain implied by the tilt factor. + +E: Fix deform volume setting is invalid + +Cannot use volume style unless other dimensions are being controlled. + +E: More than one fix deform + +Only one fix deform can be defined at a time. + +E: Variable name for fix deform does not exist + +Self-explantory. + +E: Variable for fix deform is invalid style + +The variable must be an equal-style variable. + +E: Final box dimension due to fix deform is < 0.0 + +Self-explanatory. + +E: Cannot use fix deform trate on a box with zero tilt + +The trate style alters the current strain. + +E: Fix deform cannot use yz variable with xy + +The yz setting cannot be a variable if xy deformation is also +specified. This is because LAMMPS cannot determine if the yz setting +will induce a box flip which would be invalid if xy is also changing. + +E: Fix deform is changing yz too much with xy + +When both yz and xy are changing, it induces changes in xz if the +box must flip from one tilt extreme to another. Thus it is not +allowed for yz to grow so much that a flip is induced. + +*/ diff --git a/src/KOKKOS/fix_nh_kokkos.cpp b/src/KOKKOS/fix_nh_kokkos.cpp new file mode 100755 index 0000000000..bdf1bd3917 --- /dev/null +++ b/src/KOKKOS/fix_nh_kokkos.cpp @@ -0,0 +1,738 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Stan Moore (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_nh_kokkos.h" +#include "math_extra.h" +#include "atom.h" +#include "force.h" +#include "group.h" +#include "comm.h" +#include "neighbor.h" +#include "irregular.h" +#include "modify.h" +#include "fix_deform.h" +#include "compute.h" +#include "kspace.h" +#include "update.h" +#include "respa.h" +#include "domain_kokkos.h" +#include "memory.h" +#include "error.h" +#include "atom_masks.h" +#include "atom_kokkos.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define DELTAFLIP 0.1 +#define TILTMAX 1.5 + +enum{NOBIAS,BIAS}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO,TRICLINIC}; + +/* ---------------------------------------------------------------------- + NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion + ---------------------------------------------------------------------- */ + +template +FixNHKokkos::FixNHKokkos(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg) +{ + domainKK = (DomainKokkos *) domain; + execution_space = ExecutionSpaceFromDevice::space; + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +FixNHKokkos::~FixNHKokkos() +{ + + +} + +/* ---------------------------------------------------------------------- */ + +template +void FixNHKokkos::init() +{ + FixNH::init(); + + atomKK->k_mass.modify(); + atomKK->k_mass.sync(); +} + +/* ---------------------------------------------------------------------- + compute T,P before integrator starts +------------------------------------------------------------------------- */ + +template +void FixNHKokkos::setup(int vflag) +{ + // t_target is needed by NPH and NPT in compute_scalar() + // If no thermostat or using fix nphug, + // t_target must be defined by other means. + + if (tstat_flag && strcmp(style,"nphug") != 0) { + compute_temp_target(); + } else if (pstat_flag) { + + // t0 = reference temperature for masses + // cannot be done in init() b/c temperature cannot be called there + // is b/c Modify::init() inits computes after fixes due to dof dependence + // guesstimate a unit-dependent t0 if actual T = 0.0 + // if it was read in from a restart file, leave it be + + if (t0 == 0.0) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + t0 = temperature->compute_scalar(); + if (t0 == 0.0) { + if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; + else t0 = 300.0; + } + } + t_target = t0; + } + + if (pstat_flag) compute_press_target(); + + atomKK->sync(temperature->execution_space,temperature->datamask_read); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + if (pstat_flag) { + //atomKK->sync(pressure->execution_space,pressure->datamask_read); + //atomKK->modified(pressure->execution_space,pressure->datamask_modify); + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + // masses and initial forces on thermostat variables + + if (tstat_flag) { + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) { + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - + boltz * t_target) / eta_mass[ich]; + } + } + + // masses and initial forces on barostat variables + + if (pstat_flag) { + double kt = boltz * t_target; + double nkt = atom->natoms * kt; + + for (int i = 0; i < 3; i++) + if (p_flag[i]) + omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + } + + // masses and initial forces on barostat thermostat variables + + if (mpchain) { + etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - + boltz * t_target) / etap_mass[ich]; + } + + } +} + +/* ---------------------------------------------------------------------- + 1st half of Verlet update +------------------------------------------------------------------------- */ + +template +void FixNHKokkos::initial_integrate(int vflag) +{ + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); + + // update eta_dot + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // need to recompute pressure to account for change in KE + // t_current is up-to-date, but compute_temperature is not + // compute appropriately coupled elements of mvv_current + + if (pstat_flag) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + //atomKK->sync(pressure->execution_space,pressure->datamask_read); + //atomKK->modified(pressure->execution_space,pressure->datamask_modify); + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) { + compute_press_target(); + nh_omega_dot(); + nh_v_press(); + } + + nve_v(); + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + nve_x(); + + // remap simulation box by 1/2 step + // redo KSpace coeffs since volume has changed + + if (pstat_flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + } +} + +/* ---------------------------------------------------------------------- + 2nd half of Verlet update +------------------------------------------------------------------------- */ + +template +void FixNHKokkos::final_integrate() +{ + nve_v(); + + // re-compute temp before nh_v_press() + // only needed for temperature computes with BIAS on reneighboring steps: + // b/c some biases store per-atom values (e.g. temp/profile) + // per-atom values are invalid if reneigh/comm occurred + // since temp->compute() in initial_integrate() + + if (which == BIAS && neighbor->ago == 0) + atomKK->sync(temperature->execution_space,temperature->datamask_read); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + t_current = temperature->compute_scalar(); + + if (pstat_flag) nh_v_press(); + + // compute new T,P after velocities rescaled by nh_v_press() + // compute appropriately coupled elements of mvv_current + + atomKK->sync(temperature->execution_space,temperature->datamask_read); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + if (pstat_flag) { + //atomKK->sync(pressure->execution_space,pressure->datamask_read); + //atomKK->modified(pressure->execution_space,pressure->datamask_modify); + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) nh_omega_dot(); + + // update eta_dot + // update eta_press_dot + + if (tstat_flag) nhc_temp_integrate(); + if (pstat_flag && mpchain) nhc_press_integrate(); +} + +/* ---------------------------------------------------------------------- + change box size + remap all atoms or dilate group atoms depending on allremap flag + if rigid bodies exist, scale rigid body centers-of-mass +------------------------------------------------------------------------- */ + +template +void FixNHKokkos::remap() +{ + int i; + double oldlo,oldhi; + double expfac; + + int nlocal = atom->nlocal; + double *h = domain->h; + + // omega is not used, except for book-keeping + + for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i]; + + // convert pertinent atoms and rigid bodies to lamda coords + + domainKK->x2lamda(nlocal); + //if (allremap) domainKK->x2lamda(nlocal); + //else { + // for (i = 0; i < nlocal; i++) + // if (mask[i] & dilate_group_bit) + // domain->x2lamda(x[i],x[i]); + //} + + if (nrigid) + error->all(FLERR,"Cannot (yet) use rigid bodies with fix nh and Kokkos"); + //for (i = 0; i < nrigid; i++) + // modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + // this operation corresponds to applying the + // translate and scale operations + // corresponding to the solution of the following ODE: + // + // h_dot = omega_dot * h + // + // where h_dot, omega_dot and h are all upper-triangular + // 3x3 tensors. In Voigt notation, the elements of the + // RHS product tensor are: + // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1] + // + // Ordering of operations preserves time symmetry. + + double dto2 = dto/2.0; + double dto4 = dto/4.0; + double dto8 = dto/8.0; + + // off-diagonal components, first half + + if (pstyle == TRICLINIC) { + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + if (p_flag[3]) { + expfac = exp(dto4*omega_dot[1]); + h[3] *= expfac; + h[3] += dto2*(omega_dot[3]*h[2]); + h[3] *= expfac; + } + + if (p_flag[5]) { + expfac = exp(dto4*omega_dot[0]); + h[5] *= expfac; + h[5] += dto2*(omega_dot[5]*h[1]); + h[5] *= expfac; + } + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + } + + // scale diagonal components + // scale tilt factors with cell, if set + + if (p_flag[0]) { + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + expfac = exp(dto*omega_dot[0]); + domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; + } + + if (p_flag[1]) { + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + expfac = exp(dto*omega_dot[1]); + domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; + if (scalexy) h[5] *= expfac; + } + + if (p_flag[2]) { + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + expfac = exp(dto*omega_dot[2]); + domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; + if (scalexz) h[4] *= expfac; + if (scaleyz) h[3] *= expfac; + } + + // off-diagonal components, second half + + if (pstyle == TRICLINIC) { + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + if (p_flag[3]) { + expfac = exp(dto4*omega_dot[1]); + h[3] *= expfac; + h[3] += dto2*(omega_dot[3]*h[2]); + h[3] *= expfac; + } + + if (p_flag[5]) { + expfac = exp(dto4*omega_dot[0]); + h[5] *= expfac; + h[5] += dto2*(omega_dot[5]*h[1]); + h[5] *= expfac; + } + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + } + + domain->yz = h[3]; + domain->xz = h[4]; + domain->xy = h[5]; + + // tilt factor to cell length ratio can not exceed TILTMAX in one step + + if (domain->yz < -TILTMAX*domain->yprd || + domain->yz > TILTMAX*domain->yprd || + domain->xz < -TILTMAX*domain->xprd || + domain->xz > TILTMAX*domain->xprd || + domain->xy < -TILTMAX*domain->xprd || + domain->xy > TILTMAX*domain->xprd) + error->all(FLERR,"Fix npt/nph has tilted box too far in one step - " + "periodic cell is too far from equilibrium state"); + + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + domainKK->lamda2x(nlocal); + //if (allremap) domainKK->lamda2x(nlocal); + //else { + // for (i = 0; i < nlocal; i++) + // if (mask[i] & dilate_group_bit) + // domain->lamda2x(x[i],x[i]); + //} + + //if (nrigid) + // for (i = 0; i < nrigid; i++) + // modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + perform half-step barostat scaling of velocities +-----------------------------------------------------------------------*/ + +template +void FixNHKokkos::nh_v_press() +{ + atomKK->sync(execution_space,V_MASK | MASK_MASK); + + v = atomKK->k_v.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2)); + factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2)); + factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2)); + + if (which == BIAS) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + temperature->remove_bias_all(); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + } + + copymode = 1; + if (pstyle == TRICLINIC) + Kokkos::parallel_for(Kokkos::RangePolicy >(0,nlocal),*this); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; + + atomKK->modified(execution_space,V_MASK); + + if (which == BIAS) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + temperature->restore_bias_all(); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + } + +} + +template +template +KOKKOS_INLINE_FUNCTION +void FixNHKokkos::operator()(TagFixNH_nh_v_press, const int &i) const { + if (mask[i] & groupbit) { + v(i,0) *= factor[0]; + v(i,1) *= factor[1]; + v(i,2) *= factor[2]; + if (TRICLINIC_FLAG) { + v(i,0) += -dthalf*(v(i,1)*omega_dot[5] + v(i,2)*omega_dot[4]); + v(i,1) += -dthalf*v(i,2)*omega_dot[3]; + } + v(i,0) *= factor[0]; + v(i,1) *= factor[1]; + v(i,2) *= factor[2]; + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of velocities +-----------------------------------------------------------------------*/ + +template +void FixNHKokkos::nve_v() +{ + atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | MASK_MASK | RMASS_MASK | TYPE_MASK); + atomKK->modified(execution_space,V_MASK); + + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + rmass = atomKK->rmass; + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + copymode = 1; + if (rmass) + Kokkos::parallel_for(Kokkos::RangePolicy >(0,nlocal),*this); + else + Kokkos::parallel_for(Kokkos::RangePolicy >(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; +} + +template +template +KOKKOS_INLINE_FUNCTION +void FixNHKokkos::operator()(TagFixNH_nve_v, const int &i) const { + if (RMASS) { + if (mask[i] & groupbit) { + const F_FLOAT dtfm = dtf / rmass[i]; + v(i,0) += dtfm*f(i,0); + v(i,1) += dtfm*f(i,1); + v(i,2) += dtfm*f(i,2); + } + } else { + if (mask[i] & groupbit) { + const F_FLOAT dtfm = dtf / mass[type[i]]; + v(i,0) += dtfm*f(i,0); + v(i,1) += dtfm*f(i,1); + v(i,2) += dtfm*f(i,2); + } + } +} + +/* ---------------------------------------------------------------------- + perform full-step update of positions +-----------------------------------------------------------------------*/ + +template +void FixNHKokkos::nve_x() +{ + atomKK->sync(execution_space,X_MASK | V_MASK | MASK_MASK); + atomKK->modified(execution_space,X_MASK); + + x = atomKK->k_x.view(); + v = atomKK->k_v.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + // x update by full step only for atoms in group + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; +} + +template +KOKKOS_INLINE_FUNCTION +void FixNHKokkos::operator()(TagFixNH_nve_x, const int &i) const { + if (mask[i] & groupbit) { + x(i,0) += dtv * v(i,0); + x(i,1) += dtv * v(i,1); + x(i,2) += dtv * v(i,2); + } +} + +/* ---------------------------------------------------------------------- + perform half-step thermostat scaling of velocities +-----------------------------------------------------------------------*/ + +template +void FixNHKokkos::nh_v_temp() +{ + atomKK->sync(execution_space,V_MASK | MASK_MASK); + + v = atomKK->k_v.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + if (which == BIAS) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + temperature->remove_bias_all(); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + } + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; + + atomKK->modified(execution_space,V_MASK); + + if (which == BIAS) { + atomKK->sync(temperature->execution_space,temperature->datamask_read); + temperature->restore_bias_all(); + atomKK->modified(temperature->execution_space,temperature->datamask_modify); + } +} + +template +KOKKOS_INLINE_FUNCTION +void FixNHKokkos::operator()(TagFixNH_nh_v_temp, const int &i) const { + if (mask[i] & groupbit) { + v(i,0) *= factor_eta; + v(i,1) *= factor_eta; + v(i,2) *= factor_eta; + } +} + +/* ---------------------------------------------------------------------- + if any tilt ratios exceed limits, set flip = 1 and compute new tilt values + do not flip in x or y if non-periodic (can tilt but not flip) + this is b/c the box length would be changed (dramatically) by flip + if yz tilt exceeded, adjust C vector by one B vector + if xz tilt exceeded, adjust C vector by one A vector + if xy tilt exceeded, adjust B vector by one A vector + check yz first since it may change xz, then xz check comes after + if any flip occurs, create new box in domain + image_flip() adjusts image flags due to box shape change induced by flip + remap() puts atoms outside the new box back into the new box + perform irregular on atoms in lamda coords to migrate atoms to new procs + important that image_flip comes before remap, since remap may change + image flags to new values, making eqs in doc of Domain:image_flip incorrect +------------------------------------------------------------------------- */ + +template +void FixNHKokkos::pre_exchange() +{ + double xprd = domain->xprd; + double yprd = domain->yprd; + + // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP + // this avoids immediate re-flipping due to tilt oscillations + + double xtiltmax = (0.5+DELTAFLIP)*xprd; + double ytiltmax = (0.5+DELTAFLIP)*yprd; + + int flipxy,flipxz,flipyz; + flipxy = flipxz = flipyz = 0; + + if (domain->yperiodic) { + if (domain->yz < -ytiltmax) { + domain->yz += yprd; + domain->xz += domain->xy; + flipyz = 1; + } else if (domain->yz >= ytiltmax) { + domain->yz -= yprd; + domain->xz -= domain->xy; + flipyz = -1; + } + } + + if (domain->xperiodic) { + if (domain->xz < -xtiltmax) { + domain->xz += xprd; + flipxz = 1; + } else if (domain->xz >= xtiltmax) { + domain->xz -= xprd; + flipxz = -1; + } + if (domain->xy < -xtiltmax) { + domain->xy += xprd; + flipxy = 1; + } else if (domain->xy >= xtiltmax) { + domain->xy -= xprd; + flipxy = -1; + } + } + + int flip = 0; + if (flipxy || flipxz || flipyz) flip = 1; + + if (flip) { + domain->set_global_box(); + domain->set_local_box(); + + domainKK->image_flip(flipxy,flipxz,flipyz); + + domainKK->remap_all(); + + domainKK->x2lamda(atom->nlocal); + atomKK->sync(Host,ALL_MASK); + irregular->migrate_atoms(); + atomKK->modified(Host,ALL_MASK); + domainKK->lamda2x(atom->nlocal); + } +} + +template class FixNHKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class FixNHKokkos; +#endif \ No newline at end of file diff --git a/src/KOKKOS/fix_nh_kokkos.h b/src/KOKKOS/fix_nh_kokkos.h new file mode 100755 index 0000000000..61f3918e72 --- /dev/null +++ b/src/KOKKOS/fix_nh_kokkos.h @@ -0,0 +1,86 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_NH_KOKKOS_H +#define LMP_FIX_NH_KOKKOS_H + +#include "fix_nh.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +template +struct TagFixNH_nh_v_press{}; + +template +struct TagFixNH_nve_v{}; + +struct TagFixNH_nve_x{}; + +struct TagFixNH_nh_v_temp{}; + +template +class FixNHKokkos : public FixNH { + public: + typedef DeviceType device_type; + + FixNHKokkos(class LAMMPS *, int, char **); + virtual ~FixNHKokkos(); + virtual void init(); + virtual void setup(int); + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual void pre_exchange(); + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagFixNH_nh_v_press, const int&) const; + + template + KOKKOS_INLINE_FUNCTION + void operator()(TagFixNH_nve_v, const int&) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagFixNH_nve_x, const int&) const; + + KOKKOS_INLINE_FUNCTION + void operator()(TagFixNH_nh_v_temp, const int&) const; + + protected: + virtual void remap(); + + virtual void nve_x(); // may be overwritten by child classes + virtual void nve_v(); + virtual void nh_v_press(); + virtual void nh_v_temp(); + + F_FLOAT factor[3]; + + class DomainKokkos *domainKK; + + typename ArrayTypes::t_x_array x; + typename ArrayTypes::t_v_array v; + typename ArrayTypes::t_f_array_const f; + double *rmass; + typename ArrayTypes::t_float_1d_randomread mass; + typename ArrayTypes::t_int_1d type; + typename ArrayTypes::t_int_1d mask; +}; + +} + +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/KOKKOS/fix_nph_kokkos.cpp b/src/KOKKOS/fix_nph_kokkos.cpp new file mode 100755 index 0000000000..1eddaad8ab --- /dev/null +++ b/src/KOKKOS/fix_nph_kokkos.cpp @@ -0,0 +1,74 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_nph_kokkos.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +template +FixNPHKokkos::FixNPHKokkos(LAMMPS *lmp, int narg, char **arg) : + FixNHKokkos(lmp, narg, arg) +{ + if (this->tstat_flag) + this->error->all(FLERR,"Temperature control can not be used with fix nph"); + if (!this->pstat_flag) + this->error->all(FLERR,"Pressure control must be used with fix nph"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(this->id) + 6; + this->id_temp = new char[n]; + strcpy(this->id_temp,this->id); + strcat(this->id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = this->id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/kk"; + + this->modify->add_compute(3,newarg); + delete [] newarg; + this->tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(this->id) + 7; + this->id_press = new char[n]; + strcpy(this->id_press,this->id); + strcat(this->id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = this->id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = this->id_temp; + this->modify->add_compute(4,newarg); + delete [] newarg; + this->pflag = 1; +} + +template class FixNPHKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class FixNPHKokkos; +#endif \ No newline at end of file diff --git a/src/KOKKOS/fix_nph_kokkos.h b/src/KOKKOS/fix_nph_kokkos.h new file mode 100755 index 0000000000..86f3dd0ee3 --- /dev/null +++ b/src/KOKKOS/fix_nph_kokkos.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(nph/kk,FixNPHKokkos) +FixStyle(nph/kk/device,FixNPHKokkos) +FixStyle(nph/kk/host,FixNPHKokkos) + +#else + +#ifndef LMP_FIX_NPH_KOKKOS_H +#define LMP_FIX_NPH_KOKKOS_H + +#include "fix_nh_kokkos.h" + +namespace LAMMPS_NS { + +template +class FixNPHKokkos : public FixNHKokkos { + public: + FixNPHKokkos(class LAMMPS *, int, char **); + ~FixNPHKokkos() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/KOKKOS/fix_npt_kokkos.cpp b/src/KOKKOS/fix_npt_kokkos.cpp new file mode 100755 index 0000000000..880bf0126f --- /dev/null +++ b/src/KOKKOS/fix_npt_kokkos.cpp @@ -0,0 +1,74 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_npt_kokkos.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +template +FixNPTKokkos::FixNPTKokkos(LAMMPS *lmp, int narg, char **arg) : + FixNHKokkos(lmp, narg, arg) +{ + if (!this->tstat_flag) + this->error->all(FLERR,"Temperature control must be used with fix npt"); + if (!this->pstat_flag) + this->error->all(FLERR,"Pressure control must be used with fix npt"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(this->id) + 6; + this->id_temp = new char[n]; + strcpy(this->id_temp,this->id); + strcat(this->id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = this->id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/kk"; + + this->modify->add_compute(3,newarg); + delete [] newarg; + this->tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(this->id) + 7; + this->id_press = new char[n]; + strcpy(this->id_press,this->id); + strcat(this->id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = this->id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = this->id_temp; + this->modify->add_compute(4,newarg); + delete [] newarg; + this->pflag = 1; +} + +template class FixNPTKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class FixNPTKokkos; +#endif diff --git a/src/KOKKOS/fix_npt_kokkos.h b/src/KOKKOS/fix_npt_kokkos.h new file mode 100755 index 0000000000..eddd34669b --- /dev/null +++ b/src/KOKKOS/fix_npt_kokkos.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(npt/kk,FixNPTKokkos) +FixStyle(npt/kk/device,FixNPTKokkos) +FixStyle(npt/kk/host,FixNPTKokkos) + +#else + +#ifndef LMP_FIX_NPT_KOKKOS_H +#define LMP_FIX_NPT_KOKKOS_H + +#include "fix_nh_kokkos.h" + +namespace LAMMPS_NS { + +template +class FixNPTKokkos : public FixNHKokkos { + public: + FixNPTKokkos(class LAMMPS *, int, char **); + ~FixNPTKokkos() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/KOKKOS/fix_nvt_kokkos.cpp b/src/KOKKOS/fix_nvt_kokkos.cpp new file mode 100755 index 0000000000..3045f1a4f4 --- /dev/null +++ b/src/KOKKOS/fix_nvt_kokkos.cpp @@ -0,0 +1,55 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_nvt_kokkos.h" +#include "group.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +template +FixNVTKokkos::FixNVTKokkos(LAMMPS *lmp, int narg, char **arg) : + FixNHKokkos(lmp, narg, arg) +{ + if (!this->tstat_flag) + this->error->all(FLERR,"Temperature control must be used with fix nvt"); + if (this->pstat_flag) + this->error->all(FLERR,"Pressure control can not be used with fix nvt"); + + // create a new compute temp style + // id = fix-ID + temp + + int n = strlen(this->id) + 6; + this->id_temp = new char[n]; + strcpy(this->id_temp,this->id); + strcat(this->id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = this->id_temp; + newarg[1] = this->group->names[this->igroup]; + newarg[2] = (char *) "temp/kk"; + + this->modify->add_compute(3,newarg); + delete [] newarg; + this->tflag = 1; +} + +template class FixNVTKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class FixNVTKokkos; +#endif \ No newline at end of file diff --git a/src/KOKKOS/fix_nvt_kokkos.h b/src/KOKKOS/fix_nvt_kokkos.h new file mode 100755 index 0000000000..5f6eb5c24d --- /dev/null +++ b/src/KOKKOS/fix_nvt_kokkos.h @@ -0,0 +1,44 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(nvt/kk,FixNVTKokkos) +FixStyle(nvt/kk/device,FixNVTKokkos) +FixStyle(nvt/kk/host,FixNVTKokkos) + +#else + + +#ifndef LMP_FIX_NVT_KOKKOS_H +#define LMP_FIX_NVT_KOKKOS_H + +#include "fix_nh_kokkos.h" + +namespace LAMMPS_NS { + +template +class FixNVTKokkos : public FixNHKokkos { + public: + FixNVTKokkos(class LAMMPS *, int, char **); + ~FixNVTKokkos() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/KOKKOS/fix_wall_reflect_kokkos.cpp b/src/KOKKOS/fix_wall_reflect_kokkos.cpp new file mode 100755 index 0000000000..2876b0836a --- /dev/null +++ b/src/KOKKOS/fix_wall_reflect_kokkos.cpp @@ -0,0 +1,111 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "fix_wall_reflect_kokkos.h" +#include "atom_kokkos.h" +#include "comm.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "lattice.h" +#include "input.h" +#include "variable.h" +#include "error.h" +#include "force.h" +#include "atom_masks.h" + + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5}; +enum{NONE=0,EDGE,CONSTANT,VARIABLE}; + + +/* ---------------------------------------------------------------------- */ + +template +FixWallReflectKokkos::FixWallReflectKokkos(LAMMPS *lmp, int narg, char **arg) : + FixWallReflect(lmp, narg, arg) +{ + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = X_MASK | V_MASK | MASK_MASK; + datamask_modify = X_MASK | V_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +void FixWallReflectKokkos::post_integrate() +{ + // coord = current position of wall + // evaluate variable if necessary, wrap with clear/add + + atomKK->sync(execution_space,datamask_read); + atomKK->modified(execution_space,datamask_modify); + + x = atomKK->k_x.view(); + v = atomKK->k_v.view(); + mask = atomKK->k_mask.view(); + int nlocal = atom->nlocal; + + + if (varflag) modify->clearstep_compute(); + + for (int m = 0; m < nwall; m++) { + if (wallstyle[m] == VARIABLE) { + coord = input->variable->compute_equal(varindex[m]); + if (wallwhich[m] < YLO) coord *= xscale; + else if (wallwhich[m] < ZLO) coord *= yscale; + else coord *= zscale; + } else coord = coord0[m]; + + dim = wallwhich[m] / 2; + side = wallwhich[m] % 2; + + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; + } + + if (varflag) modify->addstep_compute(update->ntimestep + 1); +} + +template +KOKKOS_INLINE_FUNCTION +void FixWallReflectKokkos::operator()(TagFixWallReflectPostIntegrate, const int &i) const { + if (mask[i] & groupbit) { + if (side == 0) { + if (x(i,dim) < coord) { + x(i,dim) = coord + (coord - x(i,dim)); + v(i,dim) = -v(i,dim); + } + } else { + if (x(i,dim) > coord) { + x(i,dim) = coord - (x(i,dim) - coord); + v(i,dim) = -v(i,dim); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +template class FixWallReflectKokkos; +#ifdef KOKKOS_HAVE_CUDA +template class FixWallReflectKokkos; +#endif diff --git a/src/KOKKOS/fix_wall_reflect_kokkos.h b/src/KOKKOS/fix_wall_reflect_kokkos.h new file mode 100755 index 0000000000..f5e28796fd --- /dev/null +++ b/src/KOKKOS/fix_wall_reflect_kokkos.h @@ -0,0 +1,62 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(wall/reflect/kk,FixWallReflectKokkos) +FixStyle(wall/reflect/kk/device,FixWallReflectKokkos) +FixStyle(wall/reflect/kk/host,FixWallReflectKokkos) + +#else + +#ifndef LMP_FIX_WALL_REFLECT_KOKKOS_H +#define LMP_FIX_WALL_REFLECT_KOKKOS_H + +#include "fix_wall_reflect.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +struct TagFixWallReflectPostIntegrate{}; + +template +class FixWallReflectKokkos : public FixWallReflect { + public: + typedef DeviceType device_type; + typedef ArrayTypes AT; + FixWallReflectKokkos(class LAMMPS *, int, char **); + void post_integrate(); + + KOKKOS_INLINE_FUNCTION + void operator()(TagFixWallReflectPostIntegrate, const int&) const; + + protected: + class AtomKokkos *atomKK; + + typename AT::t_x_array x; + typename AT::t_v_array v; + typename AT::t_int_1d_randomread mask; + + + int dim,side; + X_FLOAT coord; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/