Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85427441
fix_nh_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
Sun, Sep 29, 03:50
Size
21 KB
Mime Type
text/x-c
Expires
Tue, Oct 1, 03:50 (2 d)
Engine
blob
Format
Raw Data
Handle
21179781
Attached To
rLAMMPS lammps
fix_nh_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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
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
)
{
kokkosable
=
1
;
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
()
{
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
);
}
}
namespace
LAMMPS_NS
{
template
class
FixNHKokkos
<
LMPDeviceType
>
;
#ifdef KOKKOS_HAVE_CUDA
template
class
FixNHKokkos
<
LMPHostType
>
;
#endif
}
Event Timeline
Log In to Comment