tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'github/master' into doc-intro-updates

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-26 05:29:53 -04:00
parent 1d33bd1264 e3af1544f1
commit 6fcd4c0e30
15 changed files with 1107 additions and 375 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_compute.rst
View File

@ -152,7 +152,7 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`temp/chunk <compute_temp_chunk>`
* :doc:`temp/com <compute_temp_com>`
* :doc:`temp/cs <compute_temp_cs>`
* :doc:`temp/deform <compute_temp_deform>`
* :doc:`temp/deform (k) <compute_temp_deform>`
* :doc:`temp/deform/eff <compute_temp_deform_eff>`
* :doc:`temp/drude <compute_temp_drude>`
* :doc:`temp/eff <compute_temp_eff>`

2
doc/src/Commands_fix.rst
View File

@ -148,7 +148,7 @@ OPT.
* :doc:`nvt/body <fix_nvt_body>`
* :doc:`nvt/eff <fix_nh_eff>`
* :doc:`nvt/manifold/rattle <fix_nvt_manifold_rattle>`
* :doc:`nvt/sllod (io) <fix_nvt_sllod>`
* :doc:`nvt/sllod (iko) <fix_nvt_sllod>`
* :doc:`nvt/sllod/eff <fix_nvt_sllod_eff>`
* :doc:`nvt/sphere (o) <fix_nvt_sphere>`
* :doc:`nvt/uef <fix_nh_uef>`

3
doc/src/compute_temp_deform.rst
View File

@ -1,8 +1,11 @@
.. index:: compute temp/deform
.. index:: compute temp/deform/kk
compute temp/deform command
===========================
Accelerator Variants: *temp/deform/kk*
Syntax
""""""

3
doc/src/fix_nvt_sllod.rst
View File

@ -1,11 +1,12 @@
.. index:: fix nvt/sllod
.. index:: fix nvt/sllod/intel
.. index:: fix nvt/sllod/omp
.. index:: fix nvt/sllod/kk
fix nvt/sllod command
=====================
Accelerator Variants: *nvt/sllod/intel*, *nvt/sllod/omp*
Accelerator Variants: *nvt/sllod/intel*, *nvt/sllod/omp*, *nvt/sllod/kk*
Syntax
""""""

178
lib/gpu/lal_lj_tip4p_long.cu
View File

@ -150,7 +150,8 @@ __kernel void k_lj_tip4p_long_distrib(const __global numtyp4 *restrict x_,
engv[inum*engv_iter + i] += vM.z * (acctyp)0.5 * alpha;
}
}
} else {
}
if (itype == typeO) {
fM = ansO[i];
int iH1 = hneigh[i*4 ];
int iH2 = hneigh[i*4+1];
@ -212,7 +213,8 @@ __kernel void k_lj_tip4p_reneigh(const __global numtyp4 *restrict x_,
hneigh[i*4+1] = iH2;
hneigh[i*4+2] = -1;
}
} else {
}
if (itype == typeH) {
if (hneigh[i*4+2] != -1) {
int iI, iH;
iI = atom_mapping(map,tag[i] - 1);
@ -316,12 +318,13 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
int non_local_oxy = 0;
int iH1, iH2, iO;
if(itype == typeO) {
if (itype == typeO) {
iO = i;
iH1 = hneigh[i*4 ];
iH2 = hneigh[i*4+1];
x1 = m[iO];
} else {
}
if (itype == typeH) {
iO = hneigh[i *4 ];
iH1 = hneigh[iO*4 ];
iH2 = hneigh[iO*4+1];
@ -400,16 +403,16 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
prefactor *= qqrd2e*qtmp/r;
numtyp force_coul = r2inv*prefactor * (_erfc + EWALD_F*grij*expm2 - factor_coul);
if (itype == typeH) {
f.x += delx * force_coul;
f.y += dely * force_coul;
f.z += delz * force_coul;
f.w += 0;
} else {
if (itype == typeO) {
fO.x += delx * force_coul;
fO.y += dely * force_coul;
fO.z += delz * force_coul;
fO.w += 0;
} else {
f.x += delx * force_coul;
f.y += dely * force_coul;
f.z += delz * force_coul;
f.w += 0;
}
if (EVFLAG && eflag) {
e_coul += prefactor*(_erfc-factor_coul);
@ -419,15 +422,33 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
fd.x = delx*force_coul;
fd.y = dely*force_coul;
fd.z = delz*force_coul;
if (itype == typeH) {
if (jtype == typeH) {
virial[0] += delx*fd.x;
virial[1] += dely*fd.y;
virial[2] += delz*fd.z;
virial[3] += delx*fd.y;
virial[4] += delx*fd.z;
virial[5] += dely*fd.z;
} else {
if (itype == typeO) {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdi, vdj;
numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
numtyp4 xO; fetch4(xO,iO,pos_tex);
vdi.x = xO.x*cO + xH1.x*cH + xH2.x*cH;
vdi.y = xO.y*cO + xH1.y*cH + xH2.y*cH;
vdi.z = xO.z*cO + xH1.z*cH + xH2.z*cH;
//vdi.w = vdi.w;
if (jtype == typeO) {
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
numtyp4 xjO; fetch4(xjO,jO,pos_tex);
vdj.x = xjO.x*cO + xjH1.x*cH + xjH2.x*cH;
vdj.y = xjO.y*cO + xjH1.y*cH + xjH2.y*cH;
vdj.z = xjO.z*cO + xjH1.z*cH + xjH2.z*cH;
//vdj.w = vdj.w;
} else vdj = jx;
vO[0] += 0.5*(vdi.x - vdj.x)*fd.x;
vO[1] += 0.5*(vdi.y - vdj.y)*fd.y;
vO[2] += 0.5*(vdi.z - vdj.z)*fd.z;
vO[3] += 0.5*(vdi.x - vdj.x)*fd.y;
vO[4] += 0.5*(vdi.x - vdj.x)*fd.z;
vO[5] += 0.5*(vdi.y - vdj.y)*fd.z;
} else {
if (jtype == typeO) {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdj;
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
@ -443,32 +464,14 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
virial[3] += (ix.x - vdj.x)*fd.y;
virial[4] += (ix.x - vdj.x)*fd.z;
virial[5] += (ix.y - vdj.y)*fd.z;
} else {
virial[0] += delx*fd.x;
virial[1] += dely*fd.y;
virial[2] += delz*fd.z;
virial[3] += delx*fd.y;
virial[4] += delx*fd.z;
virial[5] += dely*fd.z;
}
} else {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdi, vdj;
numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
numtyp4 xO; fetch4(xO,iO,pos_tex);
vdi.x = xO.x*cO + xH1.x*cH + xH2.x*cH;
vdi.y = xO.y*cO + xH1.y*cH + xH2.y*cH;
vdi.z = xO.z*cO + xH1.z*cH + xH2.z*cH;
//vdi.w = vdi.w;
if (jtype != typeH) {
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
numtyp4 xjO; fetch4(xjO,jO,pos_tex);
vdj.x = xjO.x*cO + xjH1.x*cH + xjH2.x*cH;
vdj.y = xjO.y*cO + xjH1.y*cH + xjH2.y*cH;
vdj.z = xjO.z*cO + xjH1.z*cH + xjH2.z*cH;
//vdj.w = vdj.w;
} else vdj = jx;
vO[0] += 0.5*(vdi.x - vdj.x)*fd.x;
vO[1] += 0.5*(vdi.y - vdj.y)*fd.y;
vO[2] += 0.5*(vdi.z - vdj.z)*fd.z;
vO[3] += 0.5*(vdi.x - vdj.x)*fd.y;
vO[4] += 0.5*(vdi.x - vdj.x)*fd.z;
vO[5] += 0.5*(vdi.y - vdj.y)*fd.z;
}
}
}
@ -633,7 +636,7 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
}
__syncthreads();
if (ii<inum) {
if (ii<inum) {
int i, numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -646,12 +649,13 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
int non_local_oxy = 0;
int iH1, iH2, iO;
if(itype == typeO) {
if (itype == typeO) {
iO = i;
iH1 = hneigh[i*4 ];
iH2 = hneigh[i*4+1];
x1 = m[iO];
} else {
}
if (itype == typeH) {
iO = hneigh[i *4 ];
iH1 = hneigh[iO*4 ];
iH2 = hneigh[iO*4+1];
@ -730,16 +734,16 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
prefactor *= qqrd2e*qtmp/r;
numtyp force_coul = r2inv*prefactor * (_erfc + EWALD_F*grij*expm2 - factor_coul);
if (itype == typeH) {
f.x += delx * force_coul;
f.y += dely * force_coul;
f.z += delz * force_coul;
f.w += 0;
} else {
if (itype == typeO) {
fO.x += delx * force_coul;
fO.y += dely * force_coul;
fO.z += delz * force_coul;
fO.w += 0;
} else {
f.x += delx * force_coul;
f.y += dely * force_coul;
f.z += delz * force_coul;
f.w += 0;
}
if (EVFLAG && eflag) {
e_coul += prefactor*(_erfc-factor_coul);
@ -749,15 +753,33 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
fd.x = delx*force_coul;
fd.y = dely*force_coul;
fd.z = delz*force_coul;
if (itype == typeH) {
if (jtype == typeH) {
virial[0] += delx*fd.x;
virial[1] += dely*fd.y;
virial[2] += delz*fd.z;
virial[3] += delx*fd.y;
virial[4] += delx*fd.z;
virial[5] += dely*fd.z;
} else {
if (itype == typeO) {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdi, vdj;
numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
numtyp4 xO; fetch4(xO,iO,pos_tex);
vdi.x = xO.x*cO + xH1.x*cH + xH2.x*cH;
vdi.y = xO.y*cO + xH1.y*cH + xH2.y*cH;
vdi.z = xO.z*cO + xH1.z*cH + xH2.z*cH;
//vdi.w = vdi.w;
if (jtype == typeO) {
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
numtyp4 xjO; fetch4(xjO,jO,pos_tex);
vdj.x = xjO.x*cO + xjH1.x*cH + xjH2.x*cH;
vdj.y = xjO.y*cO + xjH1.y*cH + xjH2.y*cH;
vdj.z = xjO.z*cO + xjH1.z*cH + xjH2.z*cH;
//vdj.w = vdj.w;
} else vdj = jx;
vO[0] += 0.5*(vdi.x - vdj.x)*fd.x;
vO[1] += 0.5*(vdi.y - vdj.y)*fd.y;
vO[2] += 0.5*(vdi.z - vdj.z)*fd.z;
vO[3] += 0.5*(vdi.x - vdj.x)*fd.y;
vO[4] += 0.5*(vdi.x - vdj.x)*fd.z;
vO[5] += 0.5*(vdi.y - vdj.y)*fd.z;
} else {
if (jtype == typeO) {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdj;
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
@ -773,32 +795,14 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
virial[3] += (ix.x - vdj.x)*fd.y;
virial[4] += (ix.x - vdj.x)*fd.z;
virial[5] += (ix.y - vdj.y)*fd.z;
} else {
virial[0] += delx*fd.x;
virial[1] += dely*fd.y;
virial[2] += delz*fd.z;
virial[3] += delx*fd.y;
virial[4] += delx*fd.z;
virial[5] += dely*fd.z;
}
} else {
numtyp cO = 1 - alpha, cH = 0.5*alpha;
numtyp4 vdi, vdj;
numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
numtyp4 xO; fetch4(xO,iO,pos_tex);
vdi.x = xO.x*cO + xH1.x*cH + xH2.x*cH;
vdi.y = xO.y*cO + xH1.y*cH + xH2.y*cH;
vdi.z = xO.z*cO + xH1.z*cH + xH2.z*cH;
//vdi.w = vdi.w;
if (jtype != typeH) {
numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
numtyp4 xjO; fetch4(xjO,jO,pos_tex);
vdj.x = xjO.x*cO + xjH1.x*cH + xjH2.x*cH;
vdj.y = xjO.y*cO + xjH1.y*cH + xjH2.y*cH;
vdj.z = xjO.z*cO + xjH1.z*cH + xjH2.z*cH;
//vdj.w = vdj.w;
} else vdj = jx;
vO[0] += 0.5*(vdi.x - vdj.x)*fd.x;
vO[1] += 0.5*(vdi.y - vdj.y)*fd.y;
vO[2] += 0.5*(vdi.z - vdj.z)*fd.z;
vO[3] += 0.5*(vdi.x - vdj.x)*fd.y;
vO[4] += 0.5*(vdi.x - vdj.x)*fd.z;
vO[5] += 0.5*(vdi.y - vdj.y)*fd.z;
}
}
}

6
src/KOKKOS/Install.sh
View File

@ -43,7 +43,7 @@ elif (test $mode = 0) then
fi
fi
# list of files with optional dependcies
# list of files with optional dependencies
action angle_charmm_kokkos.cpp angle_charmm.cpp
action angle_charmm_kokkos.h angle_charmm.h
@ -94,6 +94,8 @@ action compute_orientorder_atom_kokkos.cpp
action compute_orientorder_atom_kokkos.h
action compute_temp_kokkos.cpp
action compute_temp_kokkos.h
action compute_temp_deform_kokkos.cpp
action compute_temp_deform_kokkos.h
action dihedral_charmm_kokkos.cpp dihedral_charmm.cpp
action dihedral_charmm_kokkos.h dihedral_charmm.h
action dihedral_class2_kokkos.cpp dihedral_class2.cpp
@ -135,6 +137,8 @@ action fix_nve_sphere_kokkos.cpp
action fix_nve_sphere_kokkos.h
action fix_nvt_kokkos.cpp
action fix_nvt_kokkos.h
action fix_nvt_sllod_kokkos.cpp
action fix_nvt_sllod_kokkos.h
action fix_property_atom_kokkos.cpp
action fix_property_atom_kokkos.h
action fix_qeq_reaxff_kokkos.cpp fix_qeq_reaxff.cpp

270
src/KOKKOS/compute_temp_deform_kokkos.cpp Normal file
View File

@ -0,0 +1,270 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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 authors: Emily Kahl (Uni. of QLD, e.kahl@uq.edu.au)
------------------------------------------------------------------------- */
#include "compute_temp_deform_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "comm.h"
#include "domain_kokkos.h"
#include "error.h"
#include "force.h"
#include "memory_kokkos.h"
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
template<class DeviceType>
ComputeTempDeformKokkos<DeviceType>::ComputeTempDeformKokkos(LAMMPS *lmp, int narg, char **arg) :
ComputeTempDeform(lmp, narg, arg)
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
domainKK = (DomainKokkos *) domain;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = V_MASK | MASK_MASK | RMASS_MASK | TYPE_MASK;
datamask_modify = EMPTY_MASK;
maxbias = 0;
}
template<class DeviceType>
ComputeTempDeformKokkos<DeviceType>::~ComputeTempDeformKokkos()
{
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
double ComputeTempDeformKokkos<DeviceType>::compute_scalar()
{
atomKK->sync(execution_space,datamask_read);
atomKK->k_mass.sync<DeviceType>();
invoked_scalar = update->ntimestep;
v = atomKK->k_v.view<DeviceType>();
x = atomKK->k_x.view<DeviceType>();
if (atomKK->rmass)
rmass = atomKK->k_rmass.view<DeviceType>();
else
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;
domainKK->x2lamda(nlocal);
h_rate = domainKK->h_rate;
h_ratelo = domainKK->h_ratelo;
copymode = 1;
if (atomKK->rmass)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformScalar<1> >(0,nlocal),*this,t_kk);
else
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformScalar<0> >(0,nlocal),*this,t_kk);
copymode = 0;
domainKK->lamda2x(nlocal);
t = t_kk.t0;
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 ComputeTempDeformKokkos<DeviceType>::operator()(TagComputeTempDeformScalar<RMASS>, const int &i, CTEMP& t_kk) const {
double vstream[3],vthermal[3];
vstream[0] = h_rate[0]*x(i,0) + h_rate[5]*x(i,1) + h_rate[4]*x(i,2) + h_ratelo[0];
vstream[1] = h_rate[1]*x(i,1) + h_rate[3]*x(i,2) + h_ratelo[1];
vstream[2] = h_rate[2]*x(i,2) + h_ratelo[2];
vthermal[0] = v(i,0) - vstream[0];
vthermal[1] = v(i,1) - vstream[1];
vthermal[2] = v(i,2) - vstream[2];
if (RMASS) {
if (mask[i] & groupbit)
t_kk.t0 += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i];
} else {
if (mask[i] & groupbit)
t_kk.t0 += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]];
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void ComputeTempDeformKokkos<DeviceType>::compute_vector()
{
atomKK->sync(execution_space,datamask_read);
int i;
invoked_vector = update->ntimestep;
v = atomKK->k_v.view<DeviceType>();
x = atomKK->k_x.view<DeviceType>();
if (atomKK->rmass)
rmass = atomKK->k_rmass.view<DeviceType>();
else
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;
domainKK->x2lamda(nlocal);
h_rate = domainKK->h_rate;
h_ratelo = domainKK->h_ratelo;
copymode = 1;
if (atomKK->rmass)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformVector<1> >(0,nlocal),*this,t_kk);
else
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformVector<0> >(0,nlocal),*this,t_kk);
copymode = 0;
domainKK->lamda2x(nlocal);
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 ComputeTempDeformKokkos<DeviceType>::operator()(TagComputeTempDeformVector<RMASS>, const int &i, CTEMP& t_kk) const {
double vstream[3],vthermal[3];
vstream[0] = h_rate[0]*x(i,0) + h_rate[5]*x(i,1) + h_rate[4]*x(i,2) + h_ratelo[0];
vstream[1] = h_rate[1]*x(i,1) + h_rate[3]*x(i,2) + h_ratelo[1];
vstream[2] = h_rate[2]*x(i,2) + h_ratelo[2];
vthermal[0] = v(i,0) - vstream[0];
vthermal[1] = v(i,1) - vstream[1];
vthermal[2] = v(i,2) - vstream[2];
if (mask[i] & groupbit) {
F_FLOAT massone = 0.0;
if (RMASS) massone = rmass[i];
else massone = mass[type[i]];
t_kk.t0 += massone * vthermal[0]*vthermal[0];
t_kk.t1 += massone * vthermal[1]*vthermal[1];
t_kk.t2 += massone * vthermal[2]*vthermal[2];
t_kk.t3 += massone * vthermal[0]*vthermal[1];
t_kk.t4 += massone * vthermal[0]*vthermal[2];
t_kk.t5 += massone * vthermal[1]*vthermal[2];
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void ComputeTempDeformKokkos<DeviceType>::remove_bias_all()
{
atomKK->sync(execution_space,datamask_read);
v = atomKK->k_v.view<DeviceType>();
x = atomKK->k_x.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atom->nlocal;
if (atom->nmax > maxbias) {
maxbias = atom->nmax;
vbiasall = typename ArrayTypes<DeviceType>::t_v_array("temp/deform/kk:vbiasall", maxbias);
}
domainKK->x2lamda(nlocal);
h_rate = domain->h_rate;
h_ratelo = domain->h_ratelo;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformRemoveBias >(0,nlocal),*this);
copymode = 0;
domainKK->lamda2x(nlocal);
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void ComputeTempDeformKokkos<DeviceType>::operator()(TagComputeTempDeformRemoveBias, const int &i) const {
if (mask[i] & groupbit) {
vbiasall(i,0) = h_rate[0]*x(i,0) + h_rate[5]*x(i,1) + h_rate[4]*x(i,2) + h_ratelo[0];
vbiasall(i,1) = h_rate[1]*x(i,1) + h_rate[3]*x(i,2) + h_ratelo[1];
vbiasall(i,2) = h_rate[2]*x(i,2) + h_ratelo[2];
v(i,0) -= vbiasall(i,0);
v(i,1) -= vbiasall(i,1);
v(i,2) -= vbiasall(i,2);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void ComputeTempDeformKokkos<DeviceType>::restore_bias_all()
{
atomKK->sync(execution_space,datamask_read);
v = atomKK->k_v.view<DeviceType>();
x = atomKK->k_x.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atom->nlocal;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagComputeTempDeformRestoreBias >(0,nlocal),*this);
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void ComputeTempDeformKokkos<DeviceType>::operator()(TagComputeTempDeformRestoreBias, const int &i) const {
if (mask[i] & groupbit) {
v(i,0) += vbiasall(i,0);
v(i,1) += vbiasall(i,1);
v(i,2) += vbiasall(i,2);
}
}
namespace LAMMPS_NS {
template class ComputeTempDeformKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class ComputeTempDeformKokkos<LMPHostType>;
#endif
}

126
src/KOKKOS/compute_temp_deform_kokkos.h Normal file
View File

@ -0,0 +1,126 @@
// clang-format off
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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
// clang-format off
ComputeStyle(temp/deform/kk,ComputeTempDeformKokkos<LMPDeviceType>);
ComputeStyle(temp/deform/kk/device,ComputeTempDeformKokkos<LMPDeviceType>);
ComputeStyle(temp/deform/kk/host,ComputeTempDeformKokkos<LMPHostType>);
// clang-format on
#else
#ifndef LMP_COMPUTE_TEMP_DEFORM_KOKKOS_H
#define LMP_COMPUTE_TEMP_DEFORM_KOKKOS_H
#include "compute_temp_deform.h"
#include "kokkos_few.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
template<int RMASS>
struct TagComputeTempDeformScalar{};
template<int RMASS>
struct TagComputeTempDeformVector{};
struct TagComputeTempDeformRemoveBias{};
struct TagComputeTempDeformRestoreBias{};
template<class DeviceType>
class ComputeTempDeformKokkos: public ComputeTempDeform {
public:
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
void 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;
}
};
typedef s_CTEMP CTEMP;
typedef CTEMP value_type;
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
ComputeTempDeformKokkos(class LAMMPS *, int, char **);
~ComputeTempDeformKokkos();
double compute_scalar();
void compute_vector();
void remove_bias_all();
void restore_bias_all();
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempDeformScalar<RMASS>, const int&, CTEMP&) const;
template<int RMASS>
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempDeformVector<RMASS>, const int&, CTEMP&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempDeformRemoveBias, const int &i) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagComputeTempDeformRestoreBias, const int &i) const;
protected:
typename ArrayTypes<DeviceType>::t_x_array_randomread x;
typename ArrayTypes<DeviceType>::t_v_array v;
typename ArrayTypes<DeviceType>::t_v_array vbiasall;
typename ArrayTypes<DeviceType>::t_float_1d_randomread 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;
class DomainKokkos *domainKK;
Few<double, 6> h_rate, h_ratelo;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature compute degrees of freedom < 0
This should not happen if you are calculating the temperature
on a valid set of atoms.
*/

21
src/KOKKOS/compute_temp_kokkos.h
View File

@ -28,6 +28,16 @@ ComputeStyle(temp/kk/host,ComputeTempKokkos<LMPHostType>);
namespace LAMMPS_NS {
template<int RMASS>
struct TagComputeTempScalar{};
template<int RMASS>
struct TagComputeTempVector{};
template<class DeviceType>
class ComputeTempKokkos : public ComputeTemp {
public:
struct s_CTEMP {
double t0, t1, t2, t3, t4, t5;
KOKKOS_INLINE_FUNCTION
@ -55,17 +65,8 @@ namespace LAMMPS_NS {
t5 += rhs.t5;
}
};
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;

158
src/KOKKOS/fix_nvt_sllod_kokkos.cpp Normal file
View File

@ -0,0 +1,158 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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 authors: Emily Kahl (Uni. of QLD, e.kahl@uq.edu.au)
------------------------------------------------------------------------- */
#include "fix_nvt_sllod_kokkos.h"
#include "atom.h"
#include "atom.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "fix_deform_kokkos.h"
#include "group.h"
#include "kokkos_few.h"
#include "math_extra.h"
#include "memory_kokkos.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
template<class DeviceType>
FixNVTSllodKokkos<DeviceType>::FixNVTSllodKokkos(LAMMPS *lmp, int narg, char **arg) :
FixNHKokkos<DeviceType>(lmp, narg, arg)
{
atomKK = (AtomKokkos *) this->atom;
this->kokkosable = 1;
this->domainKK = (DomainKokkos *) this->domain;
if (!this->tstat_flag)
this->error->all(FLERR,"Temperature control must be used with fix nvt/kk");
if (this->pstat_flag)
this->error->all(FLERR,"Pressure control can not be used with fix nvt/kk");
if (this->mtchain_default_flag) this->mtchain = 1;
this->id_temp = utils::strdup(std::string(this->id)+"_temp");
this->modify->add_compute(fmt::format("{} all temp/deform/kk",this->id_temp));
this->tcomputeflag = 1;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixNVTSllodKokkos<DeviceType>::init()
{
FixNHKokkos<DeviceType>::init();
vdelu = typename ArrayTypes<DeviceType>::t_v_array("nvt/sllod/kk:vdelu", atomKK->nlocal);
if (!this->temperature->tempbias)
this->error->all(FLERR,"Temperature for fix nvt/sllod does not have a bias");
nondeformbias = 0;
if (utils::strmatch(this->temperature->style,"^temp/deform")) nondeformbias = 1;
// check fix deform remap settings
int i;
for (i = 0; i < this->modify->nfix; i++)
if (utils::strmatch(this->modify->fix[i]->style,"^deform")) {
if (((FixDeform *) this->modify->fix[i])->remapflag != Domain::V_REMAP)
this->error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform remap option");
break;
}
if (i == this->modify->nfix)
this->error->all(FLERR,"Using fix nvt/sllod with no fix deform defined");
}
/* ----------------------------------------------------------------------
perform half-step scaling of velocities
-----------------------------------------------------------------------*/
template<class DeviceType>
void FixNVTSllodKokkos<DeviceType>::nh_v_temp()
{
// remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
// thermostat thermal velocity only
// vdelu = SLLOD correction = Hrate*Hinv*vthermal
// for non temp/deform BIAS:
// calculate temperature since some computes require temp
// computed on current nlocal atoms to remove bias
if (nondeformbias){
atomKK->sync(this->temperature->execution_space,this->temperature->datamask_read);
this->temperature->compute_scalar();
atomKK->modified(this->temperature->execution_space,this->temperature->datamask_modify);
}
v = atomKK->k_v.view<DeviceType>();
mask = atomKK->k_mask.view<DeviceType>();
int nlocal = atomKK->nlocal;
if (this->igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
double h_two[6];
MathExtra::multiply_shape_shape(this->domain->h_rate,this->domain->h_inv,h_two);
d_h_two = Few<double, 6>(h_two);
this->copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNVTSllod_temp1>(0,nlocal),*this);
this->copymode = 0;
this->temperature->remove_bias_all();
this->copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNVTSllod_temp2>(0,nlocal),*this);
this->copymode = 0;
this->temperature->restore_bias_all();
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixNVTSllodKokkos<DeviceType>::operator()(TagFixNVTSllod_temp1, const int &i) const {
if (mask[i] & this->groupbit) {
vdelu(i,0) = d_h_two[0]*v(i,0) + d_h_two[5]*v(i,1) + d_h_two[4]*v(i,2);
vdelu(i,1) = d_h_two[1]*v(i,1) + d_h_two[3]*v(i,2);
vdelu(i,2) = d_h_two[2]*v(i,2);
}
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixNVTSllodKokkos<DeviceType>::operator()(TagFixNVTSllod_temp2, const int &i) const {
if (mask[i] & this->groupbit) {
v(i,0) = v(i,0)*this->factor_eta - this->dthalf*vdelu(i,0);
v(i,1) = v(i,1)*this->factor_eta - this->dthalf*vdelu(i,1);
v(i,2) = v(i,2)*this->factor_eta - this->dthalf*vdelu(i,2);
}
}
namespace LAMMPS_NS {
template class FixNVTSllodKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class FixNVTSllodKokkos<LMPHostType>;
#endif
}

97
src/KOKKOS/fix_nvt_sllod_kokkos.h Normal file
View File

@ -0,0 +1,97 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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
// clang-format off
FixStyle(nvt/sllod/kk,FixNVTSllodKokkos<LMPDeviceType>);
FixStyle(nvt/sllod/kk/device,FixNVTSllodKokkos<LMPDeviceType>);
FixStyle(nvt/sllod/kk/host,FixNVTSllodKokkos<LMPHostType>);
// clang-format on
#else
#ifndef LMP_FIX_NVT_SLLOD_KOKKOS_H
#define LMP_FIX_NVT_SLLOD_KOKKOS_H
#include "fix_nh_kokkos.h"
#include "kokkos_few.h"
#include "kokkos_type.h"
// clang-format off
namespace LAMMPS_NS {
struct TagFixNVTSllod_temp1{};
struct TagFixNVTSllod_temp2{};
template<class DeviceType>
class FixNVTSllodKokkos : public FixNHKokkos<DeviceType> {
public:
FixNVTSllodKokkos(class LAMMPS *, int, char **);
~FixNVTSllodKokkos() {}
void init();
KOKKOS_INLINE_FUNCTION
void operator()(TagFixNVTSllod_temp1, const int& i) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagFixNVTSllod_temp2, const int& i) const;
private:
int nondeformbias;
void nh_v_temp();
protected:
typename ArrayTypes<DeviceType>::t_x_array x;
typename ArrayTypes<DeviceType>::t_v_array v;
typename ArrayTypes<DeviceType>::t_v_array vdelu;
typename ArrayTypes<DeviceType>::t_f_array_const f;
typename ArrayTypes<DeviceType>::t_float_1d rmass;
typename ArrayTypes<DeviceType>::t_float_1d mass;
typename ArrayTypes<DeviceType>::t_int_1d type;
typename ArrayTypes<DeviceType>::t_int_1d mask;
Few<double, 6> d_h_two;
class DomainKokkos *domainKK;
class AtomKokkos *atomKK;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix nvt/sllod
Self-explanatory.
E: Pressure control can not be used with fix nvt/sllod
Self-explanatory.
E: Temperature for fix nvt/sllod does not have a bias
The specified compute must compute temperature with a bias.
E: Using fix nvt/sllod with inconsistent fix deform remap option
Fix nvt/sllod requires that deforming atoms have a velocity profile
provided by "remap v" as a fix deform option.
E: Using fix nvt/sllod with no fix deform defined
Self-explanatory.
*/

455
src/KOKKOS/pppm_kokkos.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -31,6 +30,7 @@
#include "memory_kokkos.h"
#include "pair.h"
#include "remap_kokkos.h"
#include "kokkos_few.h"
#include <cmath>
@ -67,7 +67,7 @@ PPPMKokkos<DeviceType>::PPPMKokkos(LAMMPS *lmp) : PPPM(lmp)
pppmflag = 1;
group_group_enable = 0;
triclinic_support = 0;
triclinic_support = 1;
nfactors = 3;
//factors = new int[nfactors];
@ -455,80 +455,95 @@ void PPPMKokkos<DeviceType>::operator()(TagPPPM_setup4, const int &n) const
template<class DeviceType>
void PPPMKokkos<DeviceType>::setup_triclinic()
{
// int i,j,k,n;
// double *prd;
//
// // volume-dependent factors
// // adjust z dimension for 2d slab PPPM
// // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0
//
// prd = domain->prd;
//
// double xprd = prd[0];
// double yprd = prd[1];
// double zprd = prd[2];
// double zprd_slab = zprd*slab_volfactor;
// volume = xprd * yprd * zprd_slab;
//
// // use lamda (0-1) coordinates
//
// delxinv = nx_pppm;
// delyinv = ny_pppm;
// delzinv = nz_pppm;
// delvolinv = delxinv*delyinv*delzinv/volume;
//
// // d_fkx,d_fky,d_fkz for my FFT grid pts
//
// double per_i,per_j,per_k;
//
// n = 0;
// for (k = nzlo_fft; k <= nzhi_fft; k++) { // parallel_for
// per_k = k - nz_pppm*(2*k/nz_pppm);
// for (j = nylo_fft; j <= nyhi_fft; j++) {
// per_j = j - ny_pppm*(2*j/ny_pppm);
// for (i = nxlo_fft; i <= nxhi_fft; i++) {
// per_i = i - nx_pppm*(2*i/nx_pppm);
//
// double unitk_lamda[3];
// unitk_lamda[0] = 2.0*MY_PI*per_i;
// unitk_lamda[1] = 2.0*MY_PI*per_j;
// unitk_lamda[2] = 2.0*MY_PI*per_k;
// x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
// d_fkx[n] = unitk_lamda[0];
// d_fky[n] = unitk_lamda[1];
// d_fkz[n] = unitk_lamda[2];
// n++;
// }
// }
// }
//
// // virial coefficients
//
// double sqk,vterm;
//
// for (n = 0; n < nfft; n++) { // parallel_for
// sqk = d_fkx[n]*d_fkx[n] + d_fky[n]*d_fky[n] + d_fkz[n]*d_fkz[n];
// if (sqk == 0.0) {
// d_vg(n,0) = 0.0;
// d_vg(n,1) = 0.0;
// d_vg(n,2) = 0.0;
// d_vg(n,3) = 0.0;
// d_vg(n,4) = 0.0;
// d_vg(n,5) = 0.0;
// } else {
// vterm = -2.0 * (1.0/sqk + 0.25/(g_ewald*g_ewald));
// d_vg(n,0) = 1.0 + vterm*d_fkx[n]*d_fkx[n];
// d_vg(n,1) = 1.0 + vterm*d_fky[n]*d_fky[n];
// d_vg(n,2) = 1.0 + vterm*d_fkz[n]*d_fkz[n];
// d_vg(n,3) = vterm*d_fkx[n]*d_fky[n];
// d_vg(n,4) = vterm*d_fkx[n]*d_fkz[n];
// d_vg(n,5) = vterm*d_fky[n]*d_fkz[n];
// }
// }
//
// compute_gf_ik_triclinic();
double *prd;
// volume-dependent factors
// adjust z dimension for 2d slab PPPM
// z dimension for 3d PPPM is zprd since slab_volfactor = 1.0
prd = domain->prd;
// Update simulation box parameters
h = Few<double, 6>(domain->h);
h_inv = Few<double, 6>(domain->h_inv);
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
double zprd_slab = zprd*slab_volfactor;
volume = xprd * yprd * zprd_slab;
// use lamda (0-1) coordinates
delxinv = nx_pppm;
delyinv = ny_pppm;
delzinv = nz_pppm;
delvolinv = delxinv*delyinv*delzinv/volume;
numz_fft = nzhi_fft-nzlo_fft + 1;
numy_fft = nyhi_fft-nylo_fft + 1;
numx_fft = nxhi_fft-nxlo_fft + 1;
const int inum_fft = numz_fft*numy_fft*numx_fft;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_setup_triclinic1>(0,inum_fft),*this);
copymode = 0;
// virial coefficients
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_setup_triclinic2>(0,nfft),*this);
copymode = 0;
compute_gf_ik_triclinic();
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_setup_triclinic1, const int &n) const
{
int k = n/(numy_fft*numx_fft);
int j = (n - k*numy_fft*numx_fft) / numx_fft;
int i = n - k*numy_fft*numx_fft - j*numx_fft;
k += nzlo_fft;
j += nylo_fft;
i += nxlo_fft;
double per_k = k - nz_pppm*(2*k/nz_pppm);
double per_j = j - ny_pppm*(2*j/ny_pppm);
double per_i = i - nx_pppm*(2*i/nx_pppm);
double unitk_lamda[3];
unitk_lamda[0] = 2.0*MY_PI*per_i;
unitk_lamda[1] = 2.0*MY_PI*per_j;
unitk_lamda[2] = 2.0*MY_PI*per_k;
x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
d_fkx[n] = unitk_lamda[0];
d_fky[n] = unitk_lamda[1];
d_fkz[n] = unitk_lamda[2];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_setup_triclinic2, const int &n) const
{
const double sqk = d_fkx[n]*d_fkx[n] + d_fky[n]*d_fky[n] + d_fkz[n]*d_fkz[n];
if (sqk == 0.0) {
d_vg(n,0) = 0.0;
d_vg(n,1) = 0.0;
d_vg(n,2) = 0.0;
d_vg(n,3) = 0.0;
d_vg(n,4) = 0.0;
d_vg(n,5) = 0.0;
} else {
const double vterm = -2.0 * (1.0/sqk + 0.25/(g_ewald*g_ewald));
d_vg(n,0) = 1.0 + vterm*d_fkx[n]*d_fkx[n];
d_vg(n,1) = 1.0 + vterm*d_fky[n]*d_fky[n];
d_vg(n,2) = 1.0 + vterm*d_fkz[n]*d_fkz[n];
d_vg(n,3) = vterm*d_fkx[n]*d_fky[n];
d_vg(n,4) = vterm*d_fkx[n]*d_fkz[n];
d_vg(n,5) = vterm*d_fky[n]*d_fkz[n];
}
}
/* ----------------------------------------------------------------------
reset local grid arrays and communication stencils
called by fix balance b/c it changed sizes of processor sub-domains
@ -1004,7 +1019,7 @@ void PPPMKokkos<DeviceType>::set_grid_global()
tmp[0] = nx_pppm;
tmp[1] = ny_pppm;
tmp[2] = nz_pppm;
x2lamdaT(&tmp[0],&tmp[0]);
KSpace::x2lamdaT(&tmp[0],&tmp[0]);
h_x = 1.0/tmp[0];
h_y = 1.0/tmp[1];
h_z = 1.0/tmp[2];
@ -1232,7 +1247,7 @@ void PPPMKokkos<DeviceType>::compute_gf_ik()
numz_fft = nzhi_fft-nzlo_fft + 1;
numy_fft = nyhi_fft-nylo_fft + 1;
numx_fft = nxhi_fft-nxlo_fft + 1;
const int inum_fft = numz_fft*numy_fft*numx_fft;
const int inum_fft = numx_fft*numy_fft*numz_fft;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_compute_gf_ik>(0,inum_fft),*this);
@ -1311,6 +1326,10 @@ void PPPMKokkos<DeviceType>::compute_gf_ik_triclinic()
nby = static_cast<int> (tmp[1]);
nbz = static_cast<int> (tmp[2]);
// Update the local copy of the domain box tilt
h = Few<double, 6>(domain->h);
h_inv = Few<double, 6>(domain->h_inv);
twoorder = 2*order;
copymode = 1;
@ -1320,71 +1339,71 @@ void PPPMKokkos<DeviceType>::compute_gf_ik_triclinic()
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_compute_gf_ik_triclinic, const int &/*m*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_compute_gf_ik_triclinic, const int &m) const
{
//int n = (m - nzlo_fft)*(nyhi_fft+1 - nylo_fft)*(nxhi_fft+1 - nxlo_fft);
//
//const int mper = m - nz_pppm*(2*m/nz_pppm);
//const double snz = square(sin(MY_PI*mper/nz_pppm));
//
//for (int l = nylo_fft; l <= nyhi_fft; l++) {
// const int lper = l - ny_pppm*(2*l/ny_pppm);
// const double sny = square(sin(MY_PI*lper/ny_pppm));
//
// for (int k = nxlo_fft; k <= nxhi_fft; k++) {
// const int kper = k - nx_pppm*(2*k/nx_pppm);
// const double snx = square(sin(MY_PI*kper/nx_pppm));
//
// double unitk_lamda[3];
// unitk_lamda[0] = 2.0*MY_PI*kper;
// unitk_lamda[1] = 2.0*MY_PI*lper;
// unitk_lamda[2] = 2.0*MY_PI*mper;
// x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
//
// const double sqk = square(unitk_lamda[0]) + square(unitk_lamda[1]) + square(unitk_lamda[2]);
//
// if (sqk != 0.0) {
// const double numerator = 12.5663706/sqk;
// const double denominator = gf_denom(snx,sny,snz);
// double sum1 = 0.0;
//
// for (int nx = -nbx; nx <= nbx; nx++) {
// const double argx = MY_PI*kper/nx_pppm + MY_PI*nx;
// const double wx = powsinxx(argx,twoorder);
//
// for (int ny = -nby; ny <= nby; ny++) {
// const double argy = MY_PI*lper/ny_pppm + MY_PI*ny;
// const double wy = powsinxx(argy,twoorder);
//
// for (int nz = -nbz; nz <= nbz; nz++) {
// const double argz = MY_PI*mper/nz_pppm + MY_PI*nz;
// const double wz = powsinxx(argz,twoorder);
//
// double b[3];
// b[0] = 2.0*MY_PI*nx_pppm*nx;
// b[1] = 2.0*MY_PI*ny_pppm*ny;
// b[2] = 2.0*MY_PI*nz_pppm*nz;
// x2lamdaT(&b[0],&b[0]);
//
// const double qx = unitk_lamda[0]+b[0];
// const double sx = exp(-0.25*square(qx/g_ewald));
//
// const double qy = unitk_lamda[1]+b[1];
// const double sy = exp(-0.25*square(qy/g_ewald));
//
// const double qz = unitk_lamda[2]+b[2];
// const double sz = exp(-0.25*square(qz/g_ewald));
//
// const double dot1 = unitk_lamda[0]*qx + unitk_lamda[1]*qy + unitk_lamda[2]*qz;
// const double dot2 = qx*qx+qy*qy+qz*qz;
// sum1 += (dot1/dot2) * sx*sy*sz * wx*wy*wz;
// }
// }
// }
// d_greensfn[n++] = numerator*sum1/denominator;
// } else d_greensfn[n++] = 0.0;
// }
//}
int n = (m - nzlo_fft)*(nyhi_fft+1 - nylo_fft)*(nxhi_fft+1 - nxlo_fft);
const int mper = m - nz_pppm*(2*m/nz_pppm);
const double snz = square(sin(MY_PI*mper/nz_pppm));
for (int l = nylo_fft; l <= nyhi_fft; l++) {
const int lper = l - ny_pppm*(2*l/ny_pppm);
const double sny = square(sin(MY_PI*lper/ny_pppm));
for (int k = nxlo_fft; k <= nxhi_fft; k++) {
const int kper = k - nx_pppm*(2*k/nx_pppm);
const double snx = square(sin(MY_PI*kper/nx_pppm));
double unitk_lamda[3];
unitk_lamda[0] = 2.0*MY_PI*kper;
unitk_lamda[1] = 2.0*MY_PI*lper;
unitk_lamda[2] = 2.0*MY_PI*mper;
x2lamdaT(&unitk_lamda[0],&unitk_lamda[0]);
const double sqk = square(unitk_lamda[0]) + square(unitk_lamda[1]) + square(unitk_lamda[2]);
if (sqk != 0.0) {
const double numerator = 12.5663706/sqk;
const double denominator = gf_denom(snx,sny,snz);
double sum1 = 0.0;
for (int nx = -nbx; nx <= nbx; nx++) {
const double argx = MY_PI*kper/nx_pppm + MY_PI*nx;
const double wx = powsinxx(argx,twoorder);
for (int ny = -nby; ny <= nby; ny++) {
const double argy = MY_PI*lper/ny_pppm + MY_PI*ny;
const double wy = powsinxx(argy,twoorder);
for (int nz = -nbz; nz <= nbz; nz++) {
const double argz = MY_PI*mper/nz_pppm + MY_PI*nz;
const double wz = powsinxx(argz,twoorder);
double b[3];
b[0] = 2.0*MY_PI*nx_pppm*nx;
b[1] = 2.0*MY_PI*ny_pppm*ny;
b[2] = 2.0*MY_PI*nz_pppm*nz;
x2lamdaT(&b[0],&b[0]);
const double qx = unitk_lamda[0]+b[0];
const double sx = exp(-0.25*square(qx/g_ewald));
const double qy = unitk_lamda[1]+b[1];
const double sy = exp(-0.25*square(qy/g_ewald));
const double qz = unitk_lamda[2]+b[2];
const double sz = exp(-0.25*square(qz/g_ewald));
const double dot1 = unitk_lamda[0]*qx + unitk_lamda[1]*qy + unitk_lamda[2]*qz;
const double dot2 = qx*qx+qy*qy+qz*qz;
sum1 += (dot1/dot2) * sx*sy*sz * wx*wy*wz;
}
}
}
d_greensfn[n++] = numerator*sum1/denominator;
} else d_greensfn[n++] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
@ -1867,107 +1886,122 @@ void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik10, const int &ii) con
template<class DeviceType>
void PPPMKokkos<DeviceType>::poisson_ik_triclinic()
{
// int i,j,k,n;
//
// // compute gradients of V(r) in each of 3 dims by transforming ik*V(k)
// // FFT leaves data in 3d brick decomposition
// // copy it into inner portion of vdx,vdy,vdz arrays
//
// // x direction gradient
//
// n = 0;
// for (i = 0; i < nfft; i++) { // parallel_for1
// d_work2[n] = -d_fkx[i]*d_work1[n+1];
// d_work2[n+1] = d_fkx[i]*d_work1[n];
// n += 2;
// }
//
// fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
//
// n = 0;
// for (k = nzlo_in-nzlo_out; k <= nzhi_in-nzlo_out; k++) // parallel_for2
//
//
// // y direction gradient
//
// n = 0;
// for (i = 0; i < nfft; i++) { // parallel_for3
// d_work2[n] = -d_fky[i]*d_work1[n+1];
// d_work2[n+1] = d_fky[i]*d_work1[n];
// n += 2;
// }
//
// fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
//
// n = 0;
// for (k = nzlo_in-nzlo_out; k <= nzhi_in-nzlo_out; k++) // parallel_for4
// for (j = nylo_in-nylo_out; j <= nyhi_in-nylo_out; j++)
// for (i = nxlo_in-nxlo_out; i <= nxhi_in-nxlo_out; i++) {
// d_vdy_brick(k,j,i) = d_work2[n];
// n += 2;
// }
//
// // z direction gradient
//
// n = 0;
// for (i = 0; i < nfft; i++) { // parallel_for5
// d_work2[n] = -d_fkz[i]*d_work1[n+1];
// d_work2[n+1] = d_fkz[i]*d_work1[n];
// n += 2;
// }
//
// fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
//
// n = 0;
// for (k = nzlo_in-nzlo_out; k <= nzhi_in-nzlo_out; k++) // parallel_for6
//
// compute gradients of V(r) in each of 3 dims by transforming ik*V(k)
// FFT leaves data in 3d brick decomposition
// copy it into inner portion of vdx,vdy,vdz arrays
// x direction gradient
// merge three outer loops into one for better threading
numz_fft = nzhi_fft-nzlo_fft + 1;
numy_fft = nyhi_fft-nylo_fft + 1;
numx_fft = nxhi_fft-nxlo_fft + 1;
const int inum_fft = numz_fft*numy_fft*numx_fft;
numz_inout = (nzhi_in-nzlo_out)-(nzlo_in-nzlo_out) + 1;
numy_inout = (nyhi_in-nylo_out)-(nylo_in-nylo_out) + 1;
numx_inout = (nxhi_in-nxlo_out)-(nxlo_in-nxlo_out) + 1;
const int inum_inout = numz_inout*numy_inout*numx_inout;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic1>(0,nfft),*this);
copymode = 0;
fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic2>(0,inum_inout),*this);
copymode = 0;
// y direction gradient
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic3>(0,nfft),*this);
copymode = 0;
fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic4>(0,inum_inout),*this);
copymode = 0;
// z direction gradient
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic5>(0,nfft),*this);
copymode = 0;
fft2->compute(d_work2,d_work2,FFT3dKokkos<DeviceType>::BACKWARD);
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPPPM_poisson_ik_triclinic6>(0,inum_inout),*this);
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic1, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic1, const int &ii) const
{
d_work2[2*ii] = -d_fkx[ii]*d_work1[2*ii+1];
d_work2[2*ii+1] = d_fkx[ii]*d_work1[2*ii];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic2, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic2, const int &ii) const
{
// for (j = nylo_in-nylo_out; j <= nyhi_in-nylo_out; j++)
// for (i = nxlo_in-nxlo_out; i <= nxhi_in-nxlo_out; i++) {
// d_vdx_brick(k,j,i) = d_work2[n];
// n += 2;
// }
const int n = ii*2;
int k = ii/(numy_inout*numx_inout);
int j = (ii - k*numy_inout*numx_inout) / numx_inout;
int i = ii - k*numy_inout*numx_inout - j*numx_inout;
k += nzlo_in-nzlo_out;
j += nylo_in-nylo_out;
i += nxlo_in-nxlo_out;
d_vdx_brick(k,j,i) = d_work2[n];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic3, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic3, const int &ii) const
{
// int n = (k - (nzlo_in-nzlo_out))*((nyhi_in-nylo_out) - (nylo_in-nylo_out) + 1)*((nxhi_in-nxlo_out) - (nxlo_in-nxlo_out) + 1)*2;
d_work2[2*ii] = -d_fky[ii]*d_work1[2*ii+1];
d_work2[2*ii+1] = d_fky[ii]*d_work1[2*ii];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic4, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic4, const int &ii) const
{
// int n = (k - (nzlo_in-nzlo_out))*((nyhi_in-nylo_out) - (nylo_in-nylo_out) + 1)*((nxhi_in-nxlo_out) - (nxlo_in-nxlo_out) + 1)*2;
//
const int n = ii*2;
int k = ii/(numy_inout*numx_inout);
int j = (ii - k*numy_inout*numx_inout) / numx_inout;
int i = ii - k*numy_inout*numx_inout - j*numx_inout;
k += nzlo_in-nzlo_out;
j += nylo_in-nylo_out;
i += nxlo_in-nxlo_out;
d_vdy_brick(k,j,i) = d_work2[n];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic5, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic5, const int &ii) const
{
// int n = (k - (nzlo_in-nzlo_out))*((nyhi_in-nylo_out) - (nylo_in-nylo_out) + 1)*((nxhi_in-nxlo_out) - (nxlo_in-nxlo_out) + 1)*2;
//
d_work2[2*ii] = -d_fkz[ii]*d_work1[2*ii+1];
d_work2[2*ii+1] = d_fkz[ii]*d_work1[2*ii];
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic6, const int &/*k*/) const
void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic6, const int &ii) const
{
// int n = (k - (nzlo_in-nzlo_out))*((nyhi_in-nylo_out) - (nylo_in-nylo_out) + 1)*((nxhi_in-nxlo_out) - (nxlo_in-nxlo_out) + 1)*2;
//
@ -1976,6 +2010,14 @@ void PPPMKokkos<DeviceType>::operator()(TagPPPM_poisson_ik_triclinic6, const int
// d_vdz_brick(k,j,i) = d_work2[n];
// n += 2;
// }
const int n = ii*2;
int k = ii/(numy_inout*numx_inout);
int j = (ii - k*numy_inout*numx_inout) / numx_inout;
int i = ii - k*numy_inout*numx_inout - j*numx_inout;
k += nzlo_in-nzlo_out;
j += nylo_in-nylo_out;
i += nxlo_in-nxlo_out;
d_vdz_brick(k,j,i) = d_work2[n];
}
/* ----------------------------------------------------------------------
@ -2899,3 +2941,4 @@ template class PPPMKokkos<LMPDeviceType>;
template class PPPMKokkos<LMPHostType>;
#endif
}

29
src/KOKKOS/pppm_kokkos.h
View File

@ -29,6 +29,9 @@ KSpaceStyle(pppm/kk/host,PPPMKokkos<LMPHostType>);
#include "kokkos_base_fft.h"
#include "fftdata_kokkos.h"
#include "kokkos_type.h"
#include "kokkos_few.h"
// clang-format off
// fix up FFT defines for KOKKOS with CUDA
@ -55,6 +58,8 @@ struct TagPPPM_setup1{};
struct TagPPPM_setup2{};
struct TagPPPM_setup3{};
struct TagPPPM_setup4{};
struct TagPPPM_setup_triclinic1{};
struct TagPPPM_setup_triclinic2{};
struct TagPPPM_compute_gf_ik{};
struct TagPPPM_compute_gf_ik_triclinic{};
struct TagPPPM_self1{};
@ -138,6 +143,12 @@ class PPPMKokkos : public PPPM, public KokkosBaseFFT {
KOKKOS_INLINE_FUNCTION
void operator()(TagPPPM_setup4, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPPPM_setup_triclinic1, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPPPM_setup_triclinic2, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPPPM_compute_gf_ik, const int&) const;
@ -309,6 +320,23 @@ class PPPMKokkos : public PPPM, public KokkosBaseFFT {
int numx_out,numy_out,numz_out;
int ix,iy,nlocal;
// Local copies of the domain box tilt etc.
Few<double,6> h, h_inv;
KOKKOS_INLINE_FUNCTION
void x2lamdaT(double* v, double* lamda) const
{
double lamda_tmp[3];
lamda_tmp[0] = h_inv[0]*v[0];
lamda_tmp[1] = h_inv[5]*v[0] + h_inv[1]*v[1];
lamda_tmp[2] = h_inv[4]*v[0] + h_inv[3]*v[1] + h_inv[2]*v[2];
lamda[0] = lamda_tmp[0];
lamda[1] = lamda_tmp[1];
lamda[2] = lamda_tmp[2];
}
int nx,ny,nz;
typename AT::t_int_1d_um d_list_index;
typename FFT_AT::t_FFT_SCALAR_1d_um d_buf;
@ -596,3 +624,4 @@ accuracy. This error should not occur for typical problems. Please
send an email to the developers.
*/

131
src/compute_temp_deform.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -18,27 +17,25 @@
#include "compute_temp_deform.h"
#include <cstring>
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "fix_deform.h"
#include "force.h"
#include "group.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute temp/deform command");
if (narg != 3) error->all(FLERR, "Illegal compute temp/deform command");
scalar_flag = vector_flag = 1;
size_vector = 6;
@ -56,8 +53,10 @@ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) :
ComputeTempDeform::~ComputeTempDeform()
{
memory->destroy(vbiasall);
delete [] vector;
if (!copymode) {
memory->destroy(vbiasall);
delete[] vector;
}
}
/* ---------------------------------------------------------------------- */
@ -69,16 +68,14 @@ void ComputeTempDeform::init()
// check fix deform remap settings
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
if (((FixDeform *) modify->fix[i])->remapflag == Domain::X_REMAP &&
comm->me == 0)
error->warning(FLERR,"Using compute temp/deform with inconsistent "
"fix deform remap option");
if (utils::strmatch(modify->fix[i]->style, "^deform")) {
if (((FixDeform *) modify->fix[i])->remapflag == Domain::X_REMAP && comm->me == 0)
error->warning(FLERR,
"Using compute temp/deform with inconsistent fix deform remap option");
break;
}
if (i == modify->nfix && comm->me == 0)
error->warning(FLERR,
"Using compute temp/deform with no fix deform defined");
error->warning(FLERR, "Using compute temp/deform with no fix deform defined");
}
/* ---------------------------------------------------------------------- */
@ -98,15 +95,17 @@ void ComputeTempDeform::dof_compute()
natoms_temp = group->count(igroup);
dof = domain->dimension * natoms_temp;
dof -= extra_dof + fix_dof;
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
if (dof > 0)
tfactor = force->mvv2e / (dof * force->boltz);
else
tfactor = 0.0;
}
/* ---------------------------------------------------------------------- */
double ComputeTempDeform::compute_scalar()
{
double lamda[3],vstream[3],vthermal[3];
double lamda[3], vstream[3], vthermal[3];
invoked_scalar = update->ntimestep;
@ -129,26 +128,25 @@ double ComputeTempDeform::compute_scalar()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
domain->x2lamda(x[i],lamda);
vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
domain->x2lamda(x[i], lamda);
vstream[0] = h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2] * lamda[2] + h_ratelo[2];
vthermal[0] = v[i][0] - vstream[0];
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (rmass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * rmass[i];
t += (vthermal[0] * vthermal[0] + vthermal[1] * vthermal[1] + vthermal[2] * vthermal[2]) *
rmass[i];
else
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
t += (vthermal[0] * vthermal[0] + vthermal[1] * vthermal[1] + vthermal[2] * vthermal[2]) *
mass[type[i]];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
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");
error->all(FLERR, "Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}
@ -157,7 +155,7 @@ double ComputeTempDeform::compute_scalar()
void ComputeTempDeform::compute_vector()
{
double lamda[3],vstream[3],vthermal[3];
double lamda[3], vstream[3], vthermal[3];
invoked_vector = update->ntimestep;
@ -172,31 +170,32 @@ void ComputeTempDeform::compute_vector()
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
double massone,t[6];
double massone, t[6];
for (int i = 0; i < 6; i++) t[i] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
domain->x2lamda(x[i],lamda);
vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
domain->x2lamda(x[i], lamda);
vstream[0] = h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2] * lamda[2] + h_ratelo[2];
vthermal[0] = v[i][0] - vstream[0];
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
t[1] += massone * vthermal[1]*vthermal[1];
t[2] += massone * vthermal[2]*vthermal[2];
t[3] += massone * vthermal[0]*vthermal[1];
t[4] += massone * vthermal[0]*vthermal[2];
t[5] += massone * vthermal[1]*vthermal[2];
if (rmass)
massone = rmass[i];
else
massone = mass[type[i]];
t[0] += massone * vthermal[0] * vthermal[0];
t[1] += massone * vthermal[1] * vthermal[1];
t[2] += massone * vthermal[2] * vthermal[2];
t[3] += massone * vthermal[0] * vthermal[1];
t[4] += massone * vthermal[0] * vthermal[2];
t[5] += massone * vthermal[1] * vthermal[2];
}
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world);
for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
@ -210,11 +209,10 @@ void ComputeTempDeform::remove_bias(int i, double *v)
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
domain->x2lamda(atom->x[i],lamda);
vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2];
domain->x2lamda(atom->x[i], lamda);
vbias[0] = h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
vbias[1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
vbias[2] = h_rate[2] * lamda[2] + h_ratelo[2];
v[0] -= vbias[0];
v[1] -= vbias[1];
v[2] -= vbias[2];
@ -230,11 +228,10 @@ void ComputeTempDeform::remove_bias_thr(int i, double *v, double *b)
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
domain->x2lamda(atom->x[i],lamda);
b[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
b[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
b[2] = h_rate[2]*lamda[2] + h_ratelo[2];
domain->x2lamda(atom->x[i], lamda);
b[0] = h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
b[1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
b[2] = h_rate[2] * lamda[2] + h_ratelo[2];
v[0] -= b[0];
v[1] -= b[1];
v[2] -= b[2];
@ -253,7 +250,7 @@ void ComputeTempDeform::remove_bias_all()
if (atom->nmax > maxbias) {
memory->destroy(vbiasall);
maxbias = atom->nmax;
memory->create(vbiasall,maxbias,3,"temp/deform:vbiasall");
memory->create(vbiasall, maxbias, 3, "temp/deform:vbiasall");
}
double lamda[3];
@ -262,11 +259,11 @@ void ComputeTempDeform::remove_bias_all()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
domain->x2lamda(atom->x[i],lamda);
vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2];
domain->x2lamda(atom->x[i], lamda);
vbiasall[i][0] =
h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
vbiasall[i][1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
vbiasall[i][2] = h_rate[2] * lamda[2] + h_ratelo[2];
v[i][0] -= vbiasall[i][0];
v[i][1] -= vbiasall[i][1];
v[i][2] -= vbiasall[i][2];
@ -320,6 +317,6 @@ void ComputeTempDeform::restore_bias_all()
double ComputeTempDeform::memory_usage()
{
double bytes = 3*maxbias * sizeof(double);
double bytes = 3 * maxbias * sizeof(double);
return bytes;
}

1
unittest/force-styles/tests/mol-pair-lj_cut_tip4p_long.yaml
View File

@ -2,7 +2,6 @@
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:08:49 2021
epsilon: 1e-13
skip_tests: gpu
prerequisites: ! |
atom full
pair lj/cut/tip4p/long