Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64843096
fix_langevin_kokkos.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, May 29, 21:18
Size
28 KB
Mime Type
text/x-c
Expires
Fri, May 31, 21:18 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17966597
Attached To
rLAMMPS lammps
fix_langevin_kokkos.cpp
View Options
/* ----------------------------------------------------------------------
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 <math.h>
#include <stdio.h>
#include <string.h>
#include "fix_langevin_kokkos.h"
#include "atom_masks.h"
#include "atom_kokkos.h"
#include "force.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "memory.h"
#include "group.h"
#include "random_mars.h"
#include "compute.h"
#include "comm.h"
#include "modify.h"
#include "input.h"
#include "variable.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NOBIAS,BIAS};
enum{CONSTANT,EQUAL,ATOM};
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **arg) :
FixLangevin(lmp, narg, arg),rand_pool(seed + comm->me)
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
int ntypes = atomKK->ntypes;
// allocate per-type arrays for force prefactors
memory->create_kokkos(k_gfactor1,gfactor1,ntypes+1,"langevin:gfactor1");
memory->create_kokkos(k_gfactor2,gfactor2,ntypes+1,"langevin:gfactor2");
memory->create_kokkos(k_ratio,ratio,ntypes+1,"langevin:ratio");
d_gfactor1 = k_gfactor1.template view<DeviceType>();
h_gfactor1 = k_gfactor1.template view<LMPHostType>();
d_gfactor2 = k_gfactor2.template view<DeviceType>();
h_gfactor2 = k_gfactor2.template view<LMPHostType>();
d_ratio = k_ratio.template view<DeviceType>();
h_ratio = k_ratio.template view<LMPHostType>();
// optional args
for (int i = 1; i <= ntypes; i++) ratio[i] = 1.0;
k_ratio.template modify<LMPHostType>();
if(gjfflag){
nvalues = 3;
grow_arrays(atomKK->nmax);
atom->add_callback(0);
// initialize franprev to zero
for (int i = 0; i < atomKK->nlocal; i++) {
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
}
k_franprev.template modify<LMPHostType>();
}
if(zeroflag){
k_fsumall = tdual_double_1d_3n("langevin:fsumall");
h_fsumall = k_fsumall.template view<LMPHostType>();
d_fsumall = k_fsumall.template view<DeviceType>();
}
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = V_MASK | F_MASK | MASK_MASK | RMASS_MASK | TYPE_MASK;
datamask_modify = F_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixLangevinKokkos<DeviceType>::~FixLangevinKokkos()
{
memory->destroy_kokkos(k_gfactor1,gfactor1);
memory->destroy_kokkos(k_gfactor2,gfactor2);
memory->destroy_kokkos(k_ratio,ratio);
memory->destroy_kokkos(k_flangevin,flangevin);
if(gjfflag) memory->destroy_kokkos(k_franprev,franprev);
memory->destroy_kokkos(k_tforce,tforce);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::init()
{
FixLangevin::init();
if(oflag)
error->all(FLERR,"Fix langevin omega is not yet implemented with kokkos");
if(ascale)
error->all(FLERR,"Fix langevin angmom is not yet implemented with kokkos");
// prefactors are modified in the init
k_gfactor1.template modify<LMPHostType>();
k_gfactor2.template modify<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::grow_arrays(int nmax)
{
memory->grow_kokkos(k_franprev,franprev,nmax,3,"langevin:franprev");
d_franprev = k_franprev.template view<DeviceType>();
h_franprev = k_franprev.template view<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::post_force(int vflag)
{
// sync the device views which might have been modified on host
atomKK->sync(execution_space,datamask_read);
rmass = atomKK->rmass;
f = atomKK->k_f.template view<DeviceType>();
v = atomKK->k_v.template view<DeviceType>();
type = atomKK->k_type.template view<DeviceType>();
mask = atomKK->k_mask.template view<DeviceType>();
k_gfactor1.template sync<DeviceType>();
k_gfactor2.template sync<DeviceType>();
k_ratio.template sync<DeviceType>();
if(gjfflag) k_franprev.template sync<DeviceType>();
boltz = force->boltz;
dt = update->dt;
mvv2e = force->mvv2e;
ftm2v = force->ftm2v;
fran_prop_const = sqrt(24.0*boltz/t_period/dt/mvv2e);
compute_target(); // modifies tforce vector, hence sync here
k_tforce.template sync<DeviceType>();
double fsum[3],fsumall[3];
bigint count;
int nlocal = atomKK->nlocal;
if (zeroflag) {
fsum[0] = fsum[1] = fsum[2] = 0.0;
count = group->count(igroup);
if (count == 0)
error->all(FLERR,"Cannot zero Langevin force of 0 atoms");
}
// reallocate flangevin if necessary
if (tallyflag) {
if (nlocal > maxatom1) {
memory->destroy_kokkos(k_flangevin,flangevin);
maxatom1 = atomKK->nmax;
memory->create_kokkos(k_flangevin,flangevin,maxatom1,3,"langevin:flangevin");
d_flangevin = k_flangevin.template view<DeviceType>();
h_flangevin = k_flangevin.template view<LMPHostType>();
}
}
// account for bias velocity
if(tbiasflag == BIAS){
atomKK->sync(temperature->execution_space,temperature->datamask_read);
temperature->compute_scalar();
temperature->remove_bias_all(); // modifies velocities
// if temeprature compute is kokkosized host-devcie comm won't be needed
atomKK->modified(temperature->execution_space,temperature->datamask_modify);
atomKK->sync(execution_space,temperature->datamask_modify);
}
// compute langevin force in parallel on the device
FSUM s_fsum;
if (tstyle == ATOM)
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else{
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else{
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,1,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,1,0,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,1,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,1,0,0,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,1,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,1,0,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,1,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,1,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,1,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,1,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,1,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (rmass)
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,0,1,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,0,1,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
else
if (zeroflag) {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,0,0,1> post_functor(this);
Kokkos::parallel_reduce(nlocal,post_functor,s_fsum);
}
else {
FixLangevinKokkosPostForceFunctor<DeviceType,0,0,0,0,0,0> post_functor(this);
Kokkos::parallel_for(nlocal,post_functor);
}
DeviceType::fence();
if(tbiasflag == BIAS){
atomKK->sync(temperature->execution_space,temperature->datamask_read);
temperature->restore_bias_all(); // modifies velocities
atomKK->modified(temperature->execution_space,temperature->datamask_modify);
atomKK->sync(execution_space,temperature->datamask_modify);
}
// set modify flags for the views modified in post_force functor
if (gjfflag) k_franprev.template modify<DeviceType>();
if (tallyflag) k_flangevin.template modify<DeviceType>();
// set total force to zero
if (zeroflag) {
fsum[0] = s_fsum.fx; fsum[1] = s_fsum.fy; fsum[2] = s_fsum.fz;
MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
h_fsumall(0) = fsumall[0]/count;
h_fsumall(1) = fsumall[1]/count;
h_fsumall(2) = fsumall[2]/count;
k_fsumall.template modify<LMPHostType>();
k_fsumall.template sync<DeviceType>();
// set total force zero in parallel on the device
FixLangevinKokkosZeroForceFunctor<DeviceType> zero_functor(this);
Kokkos::parallel_for(nlocal,zero_functor);
DeviceType::fence();
}
// f is modified by both post_force and zero_force functors
atomKK->modified(execution_space,datamask_modify);
// thermostat omega and angmom
// if (oflag) omega_thermostat();
// if (ascale) angmom_thermostat();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
KOKKOS_INLINE_FUNCTION
FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
{
FSUM fsum;
double fdrag[3],fran[3];
double gamma1,gamma2;
double fswap;
double tsqrt_t = tsqrt;
if (mask[i] & groupbit) {
rand_type rand_gen = rand_pool.get_state();
if(Tp_TSTYLEATOM) tsqrt_t = sqrt(d_tforce[i]);
if(Tp_RMASS){
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * fran_prop_const / ftm2v;
gamma1 *= 1.0/d_ratio[type[i]];
gamma2 *= 1.0/sqrt(d_ratio[type[i]]) * tsqrt_t;
} else {
gamma1 = d_gfactor1[type[i]];
gamma2 = d_gfactor2[type[i]] * tsqrt_t;
}
fran[0] = gamma2 * (rand_gen.drand() - 0.5); //(random->uniform()-0.5);
fran[1] = gamma2 * (rand_gen.drand() - 0.5); //(random->uniform()-0.5);
fran[2] = gamma2 * (rand_gen.drand() - 0.5); //(random->uniform()-0.5);
if(Tp_BIAS){
fdrag[0] = gamma1*v(i,0);
fdrag[1] = gamma1*v(i,1);
fdrag[2] = gamma1*v(i,2);
if (v(i,0) == 0.0) fran[0] = 0.0;
if (v(i,1) == 0.0) fran[1] = 0.0;
if (v(i,2) == 0.0) fran[2] = 0.0;
}else{
fdrag[0] = gamma1*v(i,0);
fdrag[1] = gamma1*v(i,1);
fdrag[2] = gamma1*v(i,2);
}
if (Tp_GJF) {
fswap = 0.5*(fran[0]+d_franprev(i,0));
d_franprev(i,0) = fran[0];
fran[0] = fswap;
fswap = 0.5*(fran[1]+d_franprev(i,1));
d_franprev(i,1) = fran[1];
fran[1] = fswap;
fswap = 0.5*(fran[2]+d_franprev(i,2));
d_franprev(i,2) = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f(i,0) *= gjffac;
f(i,1) *= gjffac;
f(i,2) *= gjffac;
}
f(i,0) += fdrag[0] + fran[0];
f(i,1) += fdrag[1] + fran[1];
f(i,2) += fdrag[2] + fran[2];
if (Tp_TALLY) {
d_flangevin(i,0) = fdrag[0] + fran[0];
d_flangevin(i,1) = fdrag[1] + fran[1];
d_flangevin(i,2) = fdrag[2] + fran[2];
}
if (Tp_ZERO) {
fsum.fx = fran[0];
fsum.fy = fran[1];
fsum.fz = fran[2];
}
rand_pool.free_state(rand_gen);
}
return fsum;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::zero_force_item(int i) const
{
if (mask[i] & groupbit) {
f(i,0) -= d_fsumall[0];
f(i,1) -= d_fsumall[1];
f(i,2) -= d_fsumall[2];
}
}
/* ----------------------------------------------------------------------
set current t_target and t_sqrt
------------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::compute_target()
{
atomKK->sync(Host, MASK_MASK);
mask = atomKK->k_mask.template view<DeviceType>();
int nlocal = atomKK->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
// if variable temp, evaluate variable, wrap with clear/add
// reallocate tforce array if necessary
if (tstyle == CONSTANT) {
t_target = t_start + delta * (t_stop-t_start);
tsqrt = sqrt(t_target);
} else {
modify->clearstep_compute();
if (tstyle == EQUAL) {
t_target = input->variable->compute_equal(tvar);
if (t_target < 0.0)
error->one(FLERR,"Fix langevin variable returned negative temperature");
tsqrt = sqrt(t_target);
} else {
if (atom->nmax > maxatom2) {
maxatom2 = atom->nmax;
memory->destroy_kokkos(k_tforce,tforce);
memory->create_kokkos(k_tforce,tforce,maxatom2,"langevin:tforce");
d_tforce = k_tforce.template view<DeviceType>();
h_tforce = k_tforce.template view<LMPHostType>();
}
input->variable->compute_atom(tvar,igroup,tforce,1,0); // tforce is modified on host
k_tforce.template modify<LMPHostType>();
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (h_tforce[i] < 0.0)
error->one(FLERR,
"Fix langevin variable returned negative temperature");
}
modify->addstep_compute(update->ntimestep + 1);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::reset_dt()
{
if (atomKK->mass) {
for (int i = 1; i <= atomKK->ntypes; i++) {
h_gfactor2[i] = sqrt(atomKK->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
h_gfactor2[i] *= 1.0/sqrt(h_ratio[i]);
}
k_gfactor2.template modify<LMPHostType>();
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
double FixLangevinKokkos<DeviceType>::compute_scalar()
{
if (!tallyflag || flangevin == NULL) return 0.0;
v = atomKK->k_v.template view<DeviceType>();
mask = atomKK->k_mask.template view<DeviceType>();
// capture the very first energy transfer to thermal reservoir
if (update->ntimestep == update->beginstep) {
energy_onestep = 0.0;
atomKK->sync(execution_space,V_MASK | MASK_MASK);
int nlocal = atomKK->nlocal;
k_flangevin.template sync<DeviceType>();
FixLangevinKokkosTallyEnergyFunctor<DeviceType> scalar_functor(this);
Kokkos::parallel_reduce(nlocal,scalar_functor,energy_onestep);
DeviceType::fence();
energy = 0.5*energy_onestep*update->dt;
}
// convert midstep energy back to previous fullstep energy
double energy_me = energy - 0.5*energy_onestep*update->dt;
double energy_all;
MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
return -energy_all;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
double FixLangevinKokkos<DeviceType>::compute_energy_item(int i) const
{
double energy;
if (mask[i] & groupbit)
energy = d_flangevin(i,0)*v(i,0) + d_flangevin(i,1)*v(i,1) +
d_flangevin(i,2)*v(i,2);
return energy;
}
/* ----------------------------------------------------------------------
tally energy transfer to thermal reservoir
------------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::end_of_step()
{
if (!tallyflag) return;
v = atomKK->k_v.template view<DeviceType>();
mask = atomKK->k_mask.template view<DeviceType>();
atomKK->sync(execution_space,V_MASK | MASK_MASK);
int nlocal = atomKK->nlocal;
energy_onestep = 0.0;
k_flangevin.template sync<DeviceType>();
FixLangevinKokkosTallyEnergyFunctor<DeviceType> tally_functor(this);
Kokkos::parallel_reduce(nlocal,tally_functor,energy_onestep);
DeviceType::fence();
energy += energy_onestep*update->dt;
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nvalues; m++)
h_franprev(j,m) = h_franprev(i,m);
k_franprev.template modify<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::cleanup_copy()
{
random = NULL;
tstr = NULL;
gfactor1 = NULL;
gfactor2 = NULL;
ratio = NULL;
id_temp = NULL;
flangevin = NULL;
tforce = NULL;
gjfflag = 0;
franprev = NULL;
id = style = NULL;
vatom = NULL;
}
namespace LAMMPS_NS {
template class FixLangevinKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class FixLangevinKokkos<LMPHostType>;
#endif
}
Event Timeline
Log In to Comment