git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14065 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp
2015-09-24 20:35:13 +00:00
parent 546186739c
commit d2ee16c936
14 changed files with 2068 additions and 0 deletions

View File

@ -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<class DeviceType>
ComputeTempKokkos<DeviceType>::ComputeTempKokkos(LAMMPS *lmp, int narg, char **arg) :
ComputeTemp(lmp, narg, arg)
{
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = V_MASK | MASK_MASK | RMASS_MASK | TYPE_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
double ComputeTempKokkos<DeviceType>::compute_scalar()
{
atomKK->sync(execution_space,datamask_read);
invoked_scalar = update->ntimestep;
v = atomKK->k_v.view<DeviceType>();
rmass = atomKK->rmass;
mass = atomKK->k_mass.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atom->nlocal;
double t = 0.0;
CTEMP t_kk;
copymode = 1;
if (rmass)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempScalar<1> >(0,nlocal),*this,t_kk);
else
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempScalar<0> >(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<class DeviceType>
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void ComputeTempKokkos<DeviceType>::operator()(TagComputeTempScalar<RMASS>, 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<class DeviceType>
void ComputeTempKokkos<DeviceType>::compute_vector()
{
atomKK->sync(execution_space,datamask_read);
int i;
invoked_vector = update->ntimestep;
v = atomKK->k_v.view<DeviceType>();
rmass = atomKK->rmass;
mass = atomKK->k_mass.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
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<DeviceType, TagComputeTempVector<1> >(0,nlocal),*this,t_kk);
else
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempVector<0> >(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<class DeviceType>
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void ComputeTempKokkos<DeviceType>::operator()(TagComputeTempVector<RMASS>, 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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class ComputeTempKokkos<LMPHostType>;
#endif

107
src/KOKKOS/compute_temp_kokkos.h Executable file
View File

@ -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<LMPDeviceType>)
ComputeStyle(temp/kk/device,ComputeTempKokkos<LMPDeviceType>)
ComputeStyle(temp/kk/host,ComputeTempKokkos<LMPHostType>)
#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<int RMASS>
struct TagComputeTempScalar{};
template<int RMASS>
struct TagComputeTempVector{};
template<class DeviceType>
class ComputeTempKokkos : public ComputeTemp {
public:
typedef DeviceType device_type;
typedef CTEMP value_type;
typedef ArrayTypes<DeviceType> AT;
ComputeTempKokkos(class LAMMPS *, int, char **);
virtual ~ComputeTempKokkos() {}
double compute_scalar();
void compute_vector();
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempScalar<RMASS>, const int&, CTEMP&) const;
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempVector<RMASS>, const int&, CTEMP&) const;
protected:
typename ArrayTypes<DeviceType>::t_v_array_randomread v;
double *rmass;
typename ArrayTypes<DeviceType>::t_float_1d_randomread mass;
typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
typename ArrayTypes<DeviceType>::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.
*/

374
src/KOKKOS/fix_deform_kokkos.cpp Executable file
View File

@ -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();
}

104
src/KOKKOS/fix_deform_kokkos.h Executable file
View File

@ -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.
*/

738
src/KOKKOS/fix_nh_kokkos.cpp Executable file
View File

@ -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<class DeviceType>
FixNHKokkos<DeviceType>::FixNHKokkos(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg)
{
domainKK = (DomainKokkos *) domain;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixNHKokkos<DeviceType>::~FixNHKokkos()
{
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixNHKokkos<DeviceType>::init()
{
FixNH::init();
atomKK->k_mass.modify<LMPHostType>();
atomKK->k_mass.sync<LMPDeviceType>();
}
/* ----------------------------------------------------------------------
compute T,P before integrator starts
------------------------------------------------------------------------- */
template<class DeviceType>
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::nh_v_press()
{
atomKK->sync(execution_space,V_MASK | MASK_MASK);
v = atomKK->k_v.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
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<DeviceType, TagFixNH_nh_v_press<1> >(0,nlocal),*this);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNH_nh_v_press<0> >(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<class DeviceType>
template<int TRICLINIC_FLAG>
KOKKOS_INLINE_FUNCTION
void FixNHKokkos<DeviceType>::operator()(TagFixNH_nh_v_press<TRICLINIC_FLAG>, 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<class DeviceType>
void FixNHKokkos<DeviceType>::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<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
rmass = atomKK->rmass;
mass = atomKK->k_mass.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atomKK->nlocal;
if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
copymode = 1;
if (rmass)
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNH_nve_v<1> >(0,nlocal),*this);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNH_nve_v<0> >(0,nlocal),*this);
DeviceType::fence();
copymode = 0;
}
template<class DeviceType>
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void FixNHKokkos<DeviceType>::operator()(TagFixNH_nve_v<RMASS>, 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<class DeviceType>
void FixNHKokkos<DeviceType>::nve_x()
{
atomKK->sync(execution_space,X_MASK | V_MASK | MASK_MASK);
atomKK->modified(execution_space,X_MASK);
x = atomKK->k_x.view<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
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<DeviceType, TagFixNH_nve_x>(0,nlocal),*this);
DeviceType::fence();
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::nh_v_temp()
{
atomKK->sync(execution_space,V_MASK | MASK_MASK);
v = atomKK->k_v.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
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<DeviceType, TagFixNH_nh_v_temp>(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<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixNHKokkos<DeviceType>::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<class DeviceType>
void FixNHKokkos<DeviceType>::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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixNHKokkos<LMPHostType>;
#endif

86
src/KOKKOS/fix_nh_kokkos.h Executable file
View File

@ -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<int TRICLINIC_FLAG>
struct TagFixNH_nh_v_press{};
template<int RMASS>
struct TagFixNH_nve_v{};
struct TagFixNH_nve_x{};
struct TagFixNH_nh_v_temp{};
template<class DeviceType>
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<int TRICLINIC_FLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagFixNH_nh_v_press<TRICLINIC_FLAG>, const int&) const;
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void operator()(TagFixNH_nve_v<RMASS>, 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<DeviceType>::t_x_array x;
typename ArrayTypes<DeviceType>::t_v_array v;
typename ArrayTypes<DeviceType>::t_f_array_const f;
double *rmass;
typename ArrayTypes<DeviceType>::t_float_1d_randomread mass;
typename ArrayTypes<DeviceType>::t_int_1d type;
typename ArrayTypes<DeviceType>::t_int_1d mask;
};
}
#endif
/* ERROR/WARNING messages:
*/

74
src/KOKKOS/fix_nph_kokkos.cpp Executable file
View File

@ -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<class DeviceType>
FixNPHKokkos<DeviceType>::FixNPHKokkos(LAMMPS *lmp, int narg, char **arg) :
FixNHKokkos<DeviceType>(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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixNPHKokkos<LMPHostType>;
#endif

43
src/KOKKOS/fix_nph_kokkos.h Executable file
View File

@ -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<LMPDeviceType>)
FixStyle(nph/kk/device,FixNPHKokkos<LMPDeviceType>)
FixStyle(nph/kk/host,FixNPHKokkos<LMPHostType>)
#else
#ifndef LMP_FIX_NPH_KOKKOS_H
#define LMP_FIX_NPH_KOKKOS_H
#include "fix_nh_kokkos.h"
namespace LAMMPS_NS {
template<class DeviceType>
class FixNPHKokkos : public FixNHKokkos<DeviceType> {
public:
FixNPHKokkos(class LAMMPS *, int, char **);
~FixNPHKokkos() {}
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/

74
src/KOKKOS/fix_npt_kokkos.cpp Executable file
View File

@ -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<class DeviceType>
FixNPTKokkos<DeviceType>::FixNPTKokkos(LAMMPS *lmp, int narg, char **arg) :
FixNHKokkos<DeviceType>(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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixNPTKokkos<LMPHostType>;
#endif

43
src/KOKKOS/fix_npt_kokkos.h Executable file
View File

@ -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<LMPDeviceType>)
FixStyle(npt/kk/device,FixNPTKokkos<LMPDeviceType>)
FixStyle(npt/kk/host,FixNPTKokkos<LMPHostType>)
#else
#ifndef LMP_FIX_NPT_KOKKOS_H
#define LMP_FIX_NPT_KOKKOS_H
#include "fix_nh_kokkos.h"
namespace LAMMPS_NS {
template<class DeviceType>
class FixNPTKokkos : public FixNHKokkos<DeviceType> {
public:
FixNPTKokkos(class LAMMPS *, int, char **);
~FixNPTKokkos() {}
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/

55
src/KOKKOS/fix_nvt_kokkos.cpp Executable file
View File

@ -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<class DeviceType>
FixNVTKokkos<DeviceType>::FixNVTKokkos(LAMMPS *lmp, int narg, char **arg) :
FixNHKokkos<DeviceType>(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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixNVTKokkos<LMPHostType>;
#endif

44
src/KOKKOS/fix_nvt_kokkos.h Executable file
View File

@ -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<LMPDeviceType>)
FixStyle(nvt/kk/device,FixNVTKokkos<LMPDeviceType>)
FixStyle(nvt/kk/host,FixNVTKokkos<LMPHostType>)
#else
#ifndef LMP_FIX_NVT_KOKKOS_H
#define LMP_FIX_NVT_KOKKOS_H
#include "fix_nh_kokkos.h"
namespace LAMMPS_NS {
template<class DeviceType>
class FixNVTKokkos : public FixNHKokkos<DeviceType> {
public:
FixNVTKokkos(class LAMMPS *, int, char **);
~FixNVTKokkos() {}
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/

View File

@ -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<class DeviceType>
FixWallReflectKokkos<DeviceType>::FixWallReflectKokkos(LAMMPS *lmp, int narg, char **arg) :
FixWallReflect(lmp, narg, arg)
{
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | V_MASK | MASK_MASK;
datamask_modify = X_MASK | V_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixWallReflectKokkos<DeviceType>::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<DeviceType>();
v = atomKK->k_v.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
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<DeviceType, TagFixWallReflectPostIntegrate>(0,nlocal),*this);
DeviceType::fence();
copymode = 0;
}
if (varflag) modify->addstep_compute(update->ntimestep + 1);
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixWallReflectKokkos<DeviceType>::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<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixWallReflectKokkos<LMPHostType>;
#endif

View File

@ -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<LMPDeviceType>)
FixStyle(wall/reflect/kk/device,FixWallReflectKokkos<LMPDeviceType>)
FixStyle(wall/reflect/kk/host,FixWallReflectKokkos<LMPHostType>)
#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 DeviceType>
class FixWallReflectKokkos : public FixWallReflect {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> 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:
*/