Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81051599
fix_nh_intel.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
Tue, Sep 3, 23:37
Size
13 KB
Mime Type
text/x-c
Expires
Thu, Sep 5, 23:37 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20332313
Attached To
rLAMMPS lammps
fix_nh_intel.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: W. Michael Brown (Intel)
------------------------------------------------------------------------- */
#include "fix_nh_intel.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "update.h"
#include <cstring>
#include <math.h>
#include <stdio.h>
using namespace LAMMPS_NS;
using namespace FixConst;
#define TILTMAX 1.5
enum{NOBIAS,BIAS};
enum{ISO,ANISO,TRICLINIC};
typedef struct { double x,y,z; } dbl3_t;
/* ----------------------------------------------------------------------
NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
---------------------------------------------------------------------- */
FixNHIntel::FixNHIntel(LAMMPS *lmp, int narg, char **arg) :
FixNH(lmp, narg, arg)
{
_dtfm = 0;
_nlocal3 = 0;
_nlocal_max = 0;
}
/* ---------------------------------------------------------------------- */
FixNHIntel::~FixNHIntel()
{
}
/* ----------------------------------------------------------------------
change box size
remap all atoms or dilate group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixNHIntel::remap()
{
double oldlo,oldhi;
double expfac;
dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
int *mask = atom->mask;
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
const double hi0 = domain->h_inv[0];
const double hi1 = domain->h_inv[1];
const double hi2 = domain->h_inv[2];
const double hi3 = domain->h_inv[3];
const double hi4 = domain->h_inv[4];
const double hi5 = domain->h_inv[5];
const double b0 = domain->boxlo[0];
const double b1 = domain->boxlo[1];
const double b2 = domain->boxlo[2];
if (allremap) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
const double d0 = x[i].x - b0;
const double d1 = x[i].y - b1;
const double d2 = x[i].z - b2;
x[i].x = hi0*d0 + hi5*d1 + hi4*d2;
x[i].y = hi1*d1 + hi3*d2;
x[i].z = hi2*d2;
}
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
if (mask[i] & dilate_group_bit) {
const double d0 = x[i].x - b0;
const double d1 = x[i].y - b1;
const double d2 = x[i].z - b2;
x[i].x = hi0*d0 + hi5*d1 + hi4*d2;
x[i].y = hi1*d1 + hi3*d2;
x[i].z = hi2*d2;
}
}
}
if (nrigid)
for (int 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
const double h0 = domain->h[0];
const double h1 = domain->h[1];
const double h2 = domain->h[2];
const double h3 = domain->h[3];
const double h4 = domain->h[4];
const double h5 = domain->h[5];
const double nb0 = domain->boxlo[0];
const double nb1 = domain->boxlo[1];
const double nb2 = domain->boxlo[2];
if (allremap) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
x[i].x = h0*x[i].x + h5*x[i].y + h4*x[i].z + nb0;
x[i].y = h1*x[i].y + h3*x[i].z + nb1;
x[i].z = h2*x[i].z + nb2;
}
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
if (mask[i] & dilate_group_bit) {
x[i].x = h0*x[i].x + h5*x[i].y + h4*x[i].z + nb0;
x[i].y = h1*x[i].y + h3*x[i].z + nb1;
x[i].z = h2*x[i].z + nb2;
}
}
}
if (nrigid)
for (int i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ---------------------------------------------------------------------- */
void FixNHIntel::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
dto = dthalf;
// If using respa, then remap is performed in innermost level
if (strstr(update->integrate_style,"respa"))
dto = 0.5*step_respa[0];
if (pstat_flag)
pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
if (tstat_flag)
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
const int * const mask = atom->mask;
const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst :
atom->nlocal;
if (nlocal > _nlocal_max) {
if (_nlocal_max) memory->destroy(_dtfm);
_nlocal_max = static_cast<int>(1.20 * nlocal);
memory->create(_dtfm, _nlocal_max * 3, "fix_nve_intel:dtfm");
}
_nlocal3 = nlocal * 3;
if (igroup == 0) {
if (atom->rmass) {
const double * const rmass = atom->rmass;
int n = 0;
for (int i = 0; i < nlocal; i++) {
_dtfm[n++] = dtf / rmass[i];
_dtfm[n++] = dtf / rmass[i];
_dtfm[n++] = dtf / rmass[i];
}
} else {
const double * const mass = atom->mass;
const int * const type = atom->type;
int n = 0;
for (int i = 0; i < nlocal; i++) {
_dtfm[n++] = dtf / mass[type[i]];
_dtfm[n++] = dtf / mass[type[i]];
_dtfm[n++] = dtf / mass[type[i]];
}
}
} else {
if (atom->rmass) {
const double * const rmass = atom->rmass;
int n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
_dtfm[n++] = dtf / rmass[i];
_dtfm[n++] = dtf / rmass[i];
_dtfm[n++] = dtf / rmass[i];
} else {
_dtfm[n++] = 0.0;
_dtfm[n++] = 0.0;
_dtfm[n++] = 0.0;
}
} else {
const double * const mass = atom->mass;
const int * const type = atom->type;
int n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
_dtfm[n++] = dtf / mass[type[i]];
_dtfm[n++] = dtf / mass[type[i]];
_dtfm[n++] = dtf / mass[type[i]];
} else {
_dtfm[n++] = 0.0;
_dtfm[n++] = 0.0;
_dtfm[n++] = 0.0;
}
}
}
}
/* ----------------------------------------------------------------------
perform half-step barostat scaling of velocities
-----------------------------------------------------------------------*/
void FixNHIntel::nh_v_press()
{
if (pstyle == TRICLINIC || which == BIAS) {
FixNH::nh_v_press();
return;
}
dbl3_t * _noalias const v = (dbl3_t *)atom->v[0];
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double f0 = exp(-dt4*(omega_dot[0]+mtk_term2));
double f1 = exp(-dt4*(omega_dot[1]+mtk_term2));
double f2 = exp(-dt4*(omega_dot[2]+mtk_term2));
f0 *= f0;
f1 *= f1;
f2 *= f2;
if (igroup == 0) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
v[i].x *= f0;
v[i].y *= f1;
v[i].z *= f2;
}
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i].x *= f0;
v[i].y *= f1;
v[i].z *= f2;
}
}
}
}
/* ----------------------------------------------------------------------
perform half-step update of velocities
-----------------------------------------------------------------------*/
void FixNHIntel::nve_v()
{
if (neighbor->ago == 0) reset_dt();
double * _noalias const v = atom->v[0];
const double * _noalias const f = atom->f[0];
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < _nlocal3; i++)
v[i] += _dtfm[i] * f[i];
}
/* ----------------------------------------------------------------------
perform full-step update of positions
-----------------------------------------------------------------------*/
void FixNHIntel::nve_x()
{
double * _noalias const x = atom->x[0];
double * _noalias const v = atom->v[0];
const double * _noalias const f = atom->f[0];
// x update by full step only for atoms in group
if (igroup == 0) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < _nlocal3; i++)
x[i] += dtv * v[i];
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < _nlocal3; i++) {
if (_dtfm[i] != 0.0)
x[i] += dtv * v[i];
}
}
}
/* ----------------------------------------------------------------------
perform half-step thermostat scaling of velocities
-----------------------------------------------------------------------*/
void FixNHIntel::nh_v_temp()
{
if (which == BIAS) {
FixNH::nh_v_temp();
return;
}
double * _noalias const v = atom->v[0];
if (igroup == 0) {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < _nlocal3; i++)
v[i] *= factor_eta;
} else {
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for (int i = 0; i < _nlocal3; i++) {
if (_dtfm[i] != 0.0)
v[i] *= factor_eta;
}
}
}
double FixNHIntel::memory_usage()
{
return FixNH::memory_usage() + _nlocal_max * 3 * sizeof(double);
}
Event Timeline
Log In to Comment