jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #562 from lammps/fix-external 

additional fix external hooks for calling programs
This commit is contained in:
sjplimp 2017-07-05 14:46:10 -06:00 committed by GitHub
parent d0a397d6cb 4de9cec1b6
commit be8360ac4b
6 changed files with 431 additions and 197 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
doc/src/fix_msst.txt
View File

@ -17,19 +17,22 @@ msst = style name of this fix :l
dir = {x} or {y} or {z} :l
shockvel = shock velocity (strictly positive, distance/time units) :l
zero or more keyword value pairs may be appended :l
keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} :l
keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} or {beta} or {dftb} :l
{q} value = cell mass-like parameter (mass^2/distance^4 units)
{mu} value = artificial viscosity (mass/length/time units)
{p0} value = initial pressure in the shock equations (pressure units)
{v0} value = initial simulation cell volume in the shock equations (distance^3 units)
{e0} value = initial total energy (energy units)
{tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0) :pre
{tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0)
{dftb} value = {yes} or {no} for whether using MSST in conjunction with DFTB+
{beta} value = scale factor on energy contribution of DFTB+ :pre
:ule
[Examples:]
fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5
fix 2 all msst z 50.0 q 1.0e4 mu 1.0e4 v0 4.3419e+03 p0 3.7797e+03 e0 -9.72360e+02 tscale 0.01 :pre
fix 2 all msst z 50.0 q 1.0e4 mu 1.0e4 v0 4.3419e+03 p0 3.7797e+03 e0 -9.72360e+02 tscale 0.01
fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5 dftb yes beta 0.5 :pre
[Description:]
@ -58,11 +61,11 @@ oscillations have physical significance in some cases. The optional
symmetry to equilibrate to the shock Hugoniot and Rayleigh line more
rapidly in such cases.
{tscale} is a factor between 0 and 1 that determines what fraction of
thermal kinetic energy is converted to compressive strain kinetic
energy at the start of the simulation. Setting this parameter to a
non-zero value may assist in compression at the start of simulations
where it is slow to occur.
The keyword {tscale} is a factor between 0 and 1 that determines what
fraction of thermal kinetic energy is converted to compressive strain
kinetic energy at the start of the simulation. Setting this parameter
to a non-zero value may assist in compression at the start of
simulations where it is slow to occur.
If keywords {e0}, {p0},or {v0} are not supplied, these quantities will
be calculated on the first step, after the energy specified by
@ -77,17 +80,40 @@ For all pressure styles, the simulation box stays orthogonal in shape.
Parrinello-Rahman boundary conditions (tilted box) are supported by
LAMMPS, but are not implemented for MSST.
This fix computes a temperature and pressure each timestep. To do
this, the fix creates its own computes of style "temp" and "pressure",
as if these commands had been issued:
This fix computes a temperature and pressure and potential energy each
timestep. To do this, the fix creates its own computes of style "temp"
"pressure", and "pe", as if these commands had been issued:
compute fix-ID_temp group-ID temp
compute fix-ID_press group-ID pressure fix-ID_temp :pre
compute fix-ID_MSST_temp all temp
compute fix-ID_MSST_press all pressure fix-ID_MSST_temp :pre
compute fix-ID_MSST_pe all pe :pre
See the "compute temp"_compute_temp.html and "compute
pressure"_compute_pressure.html commands for details. Note that the
IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID
+ underscore + "press". The group for the new computes is "all".
IDs of the new computes are the fix-ID + "_MSST_temp" or "_MSST_press"
or "_MSST_pe". The group for the new computes is "all".
:line
The {dftb} and {beta} keywords are to allow this fix to be used when
LAMMPS is being driven by DFTB+, a density-functional tight-binding
code.
If the keyword {dftb} is used with a value of {yes}, then the MSST
equations are altered to account for an energy contribution compute by
DFTB+. In this case, you must define a "fix
external"_fix_external.html command in your input script, which is
used to callback to DFTB+ during the LAMMPS timestepping. DFTB+ will
communicate its info to LAMMPS via that fix.
The keyword {beta} is a scale factor on the DFTB+ energy contribution.
The value of {beta} must be between 0.0 and 1.0 inclusive. A value of
0.0 means no contribution, a value of 1.0 means a full contribution.
(July 2017) More information about these keywords and the use of
LAMMPS with DFTB+ will be added to the LAMMMPS documention soon.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
@ -149,8 +175,9 @@ all.
[Default:]
The keyword defaults are q = 10, mu = 0, tscale = 0.01. p0, v0, and e0
are calculated on the first step.
The keyword defaults are q = 10, mu = 0, tscale = 0.01, dftb = no,
beta = 0.0. Note that p0, v0, and e0 are calculated on the first
timestep.
:line

454
src/SHOCK/fix_msst.cpp
View File

@ -13,8 +13,8 @@
/* ----------------------------------------------------------------------
Contributing authors: Laurence Fried (LLNL), Evan Reed (LLNL, Stanford)
implementation of the Multi-Scale Shock Method
See Reed, Fried, Joannopoulos, Phys Rev Lett, 90, 235503 (2003)
implementation of the Multi-Scale Shock Method
see Reed, Fried, Joannopoulos, Phys Rev Lett, 90, 235503 (2003)
------------------------------------------------------------------------- */
#include <string.h>
@ -26,13 +26,15 @@
#include "comm.h"
#include "output.h"
#include "modify.h"
#include "fix_external.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
#include "thermo.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -42,7 +44,7 @@ using namespace FixConst;
FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), old_velocity(NULL), rfix(NULL),
id_temp(NULL), id_press(NULL), id_pe(NULL), temperature(NULL),
pressure(NULL), pe(NULL), atoms_allocated(0)
pressure(NULL), pe(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix msst command");
@ -63,6 +65,11 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
p0 = 0.0;
v0 = 1.0;
e0 = 0.0;
TS_int = 0;
T0S0 = 0.0;
S_elec = 0.0;
S_elec_1 = 0.0;
S_elec_2 = 0.0;
qmass = 1.0e1;
mu = 0.0;
@ -71,48 +78,67 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
v0_set = 0;
e0_set = 0;
tscale = 0.01;
dftb = 0;
beta = 0.0;
if ( strcmp(arg[3],"x") == 0 )
direction = 0;
else if ( strcmp(arg[3],"y") == 0 )
direction = 1;
else if ( strcmp(arg[3],"z") == 0 )
direction = 2;
else {
error->all(FLERR,"Illegal fix msst command");
}
if (strcmp(arg[3],"x") == 0) direction = 0;
else if (strcmp(arg[3],"y") == 0) direction = 1;
else if (strcmp(arg[3],"z") == 0) direction = 2;
else error->all(FLERR,"Illegal fix msst command");
velocity = force->numeric(FLERR,arg[4]);
if ( velocity < 0 )
error->all(FLERR,"Illegal fix msst command");
if (velocity < 0) error->all(FLERR,"Illegal fix msst command");
for ( int iarg = 5; iarg < narg; iarg++ ) {
if ( strcmp(arg[iarg],"q") == 0 ) {
// optional args
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg],"q") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
qmass = force->numeric(FLERR,arg[iarg+1]);
iarg++;
} else if ( strcmp(arg[iarg],"mu") == 0 ) {
iarg += 2;
} else if (strcmp(arg[iarg],"mu") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
mu = force->numeric(FLERR,arg[iarg+1]);
iarg++;
} else if ( strcmp(arg[iarg],"p0") == 0 ) {
iarg += 2;
} else if (strcmp(arg[iarg],"p0") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
p0 = force->numeric(FLERR,arg[iarg+1]);
iarg++;
p0_set = 1;
} else if ( strcmp(arg[iarg],"v0") == 0 ) {
iarg += 2;
} else if (strcmp(arg[iarg],"v0") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
v0 = force->numeric(FLERR,arg[iarg+1]);
v0_set = 1;
iarg++;
} else if ( strcmp(arg[iarg],"e0") == 0 ) {
iarg += 2;
} else if (strcmp(arg[iarg],"e0") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
e0 = force->numeric(FLERR,arg[iarg+1]);
e0_set = 1;
iarg++;
} else if ( strcmp(arg[iarg],"tscale") == 0 ) {
iarg += 2;
} else if (strcmp(arg[iarg],"tscale") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
tscale = force->numeric(FLERR,arg[iarg+1]);
if (tscale < 0.0 || tscale > 1.0)
error->all(FLERR,"Fix msst tscale must satisfy 0 <= tscale < 1");
iarg++;
iarg += 2;
} else if (strcmp(arg[iarg],"dftb") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
if (strcmp(arg[iarg+1],"yes") == 0) dftb = 1;
else if (strcmp(arg[iarg+1],"yes") == 0) dftb = 0;
else error->all(FLERR,"Illegal fix msst command");
iarg += 2;
} else if (strcmp(arg[iarg],"beta") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command");
beta = force->numeric(FLERR,arg[iarg+1]);
if (beta < 0.0 || beta > 1.0)
error->all(FLERR,"Illegal fix msst command");
iarg += 2;
} else error->all(FLERR,"Illegal fix msst command");
}
// output MSST info
if (comm->me == 0) {
if (screen) {
fprintf(screen,"MSST parameters:\n");
@ -166,15 +192,15 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
if (domain->nonperiodic) error->all(FLERR,"Fix msst requires a periodic box");
// create a new compute temp style
// id = fix-ID + temp
// create a new temperature compute
// id = fix-ID + "MSST_temp"
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
int n = strlen(id) + 10;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
strcat(id_temp,"MSST_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
@ -184,14 +210,14 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
delete [] newarg;
tflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// create a new pressure compute
// id = fix-ID + "MSST_press", compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
n = strlen(id) + 11;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
strcat(id_press,"MSST_press");
newarg = new char*[4];
newarg[0] = id_press;
@ -202,33 +228,30 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) :
delete [] newarg;
pflag = 1;
// create a new compute potential energy compute
// create a new potential energy compute
// id = fix-ID + "MSST_pe", compute group = all
n = strlen(id) + 4;
n = strlen(id) + 8;
id_pe = new char[n];
strcpy(id_pe,id);
strcat(id_pe,"_pe");
strcat(id_pe,"MSST_pe");
newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char*)"all";
newarg[2] = (char*)"pe";
newarg[1] = (char*) "all";
newarg[2] = (char*) "pe";
modify->add_compute(3,newarg);
delete [] newarg;
peflag = 1;
// initialize the time derivative of the volume.
omega[0] = omega[1] = omega[2] = 0.0;
// initialize the time derivative of the volume
omega[0] = omega[1] = omega[2] = 0.0;
nrigid = 0;
rfix = NULL;
old_velocity = new double* [atom->nlocal];
for ( int j = 0; j < atom->nlocal; j++ ) {
old_velocity[j] = new double [3];
}
atoms_allocated = atom->nlocal;
maxold = -1;
old_velocity = NULL;
}
/* ---------------------------------------------------------------------- */
@ -247,11 +270,7 @@ FixMSST::~FixMSST()
delete [] id_press;
delete [] id_pe;
for ( int j = 0; j < atoms_allocated; j++ ) {
delete [] old_velocity[j];
}
delete [] old_velocity;
memory->destroy(old_velocity);
}
/* ---------------------------------------------------------------------- */
@ -323,6 +342,15 @@ void FixMSST::init()
strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i;
}
// find fix external being used to drive LAMMPS from DFTB+
if (dftb) {
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"external") == 0)
fix_external = (FixExternal *) modify->fix[i];
if (fix_external == NULL)
error->all(FLERR,"Fix msst dftb cannot be used w/out fix external");
}
}
/* ----------------------------------------------------------------------
@ -415,10 +443,10 @@ void FixMSST::setup(int vflag)
void FixMSST::initial_integrate(int vflag)
{
int sd;
double p_msst; // MSST driving pressure.
int i, k;
double vol;
int i,k;
double p_msst; // MSST driving pressure
double vol,TS,TS_term,escale_term;
int nlocal = atom->nlocal;
int *mask = atom->mask;
double **v = atom->v;
@ -426,21 +454,51 @@ void FixMSST::initial_integrate(int vflag)
double *mass = atom->mass;
int *type = atom->type;
double **x = atom->x;
int sd = direction;
// check to see if old_velocity is correctly allocated
// realloc old_velocity if necessary
check_alloc(nlocal);
if (nlocal > maxold) {
memory->destroy(old_velocity);
maxold = atom->nmax;
memory->create(old_velocity,maxold,3,"msst:old_velocity");
}
sd = direction;
// for DFTB, extract TS_dftb from fix external
// must convert energy to mv^2 units
if (dftb) {
TS_dftb = fix_external->compute_vector(0);
TS = force->ftm2v*TS_dftb;
if (update->ntimestep == 1) T0S0 = TS;
} else {
TS = 0.0;
T0S0 = 0.0;
}
// compute new pressure and volume
// compute new pressure and volume.
temperature->compute_vector();
pressure->compute_vector();
couple();
vol = compute_vol();
// update S_elec terms and compute TS_dot via finite differences
S_elec_2 = S_elec_1;
S_elec_1 = S_elec;
double Temp = temperature->compute_scalar();
S_elec = TS/Temp;
TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt);
TS_int += (update->dt*TS_dot);
//TS_int += (update->dt*TS_dot)/total_mass;
// compute etot + extra terms for conserved quantity
double e_scale = compute_etotal() + compute_scalar();
// propagate the time derivative of
// the volume 1/2 step at fixed vol, r, rdot.
// the volume 1/2 step at fixed vol, r, rdot
p_msst = nktv2p * mvv2e * velocity * velocity * total_mass *
( v0 - vol)/( v0 * v0);
@ -448,13 +506,11 @@ void FixMSST::initial_integrate(int vflag)
(qmass * nktv2p * mvv2e);
double B = total_mass * mu / ( qmass * vol );
// prevent blow-up of the volume.
// prevent blow-up of the volume
if ( vol > v0 && A > 0.0 ) {
A = -A;
}
if (vol > v0 && A > 0.0) A = -A;
// use taylor expansion to avoid singularity at B == 0.
// use Taylor expansion to avoid singularity at B = 0
if ( B * dthalf > 1.0e-06 ) {
omega[sd] = ( omega[sd] + A * ( exp(B * dthalf) - 1.0 ) / B )
@ -465,73 +521,124 @@ void FixMSST::initial_integrate(int vflag)
}
// propagate velocity sum 1/2 step by
// temporarily propagating the velocities.
// temporarily propagating the velocities
velocity_sum = compute_vsum();
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
old_velocity[i][k] = v[i][k];
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
if (dftb) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
TS_term = TS_dot/(mass[type[i]]*velocity_sum);
escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
D += escale_term - TS_term;
old_velocity[i][k] = v[i][k];
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
old_velocity[i][k] = v[i][k];
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
}
}
}
velocity_sum = compute_vsum();
// reset the velocities.
// reset the velocities
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
v[i][k] = old_velocity[i][k];
}
v[i][0] = old_velocity[i][0];
v[i][1] = old_velocity[i][1];
v[i][2] = old_velocity[i][2];
}
}
// propagate velocities 1/2 step using the new velocity sum.
// propagate velocities 1/2 step using the new velocity sum
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
if (dftb) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
TS_term = TS_dot/(mass[type[i]]*velocity_sum);
escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
D += escale_term - TS_term;
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
}
}
}
// propagate the volume 1/2 step.
// propagate the volume 1/2 step
double vol1 = vol + omega[sd] * dthalf;
// rescale positions and change box size.
// rescale positions and change box size
dilation[sd] = vol1/vol;
remap(0);
// propagate particle positions 1 time step.
// propagate particle positions 1 time step
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -541,11 +648,11 @@ void FixMSST::initial_integrate(int vflag)
}
}
// propagate the volume 1/2 step.
// propagate the volume 1/2 step
double vol2 = vol1 + omega[sd] * dthalf;
// rescale positions and change box size.
// rescale positions and change box size
dilation[sd] = vol2/vol1;
remap(0);
@ -560,6 +667,8 @@ void FixMSST::initial_integrate(int vflag)
void FixMSST::final_integrate()
{
int i;
double p_msst; // MSST driving pressure
double TS,TS_term,escale_term;
// v update only for atoms in MSST group
@ -570,32 +679,68 @@ void FixMSST::final_integrate()
int *mask = atom->mask;
int nlocal = atom->nlocal;
double vol = compute_vol();
double p_msst;
int sd = direction;
// propagate particle velocities 1/2 step.
// for DFTB, extract TS_dftb from fix external
// must convert energy to mv^2 units
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( int k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
if (dftb) {
TS_dftb = fix_external->compute_vector(0);
TS = force->ftm2v*TS_dftb;
} else TS = 0.0;
// compute etot + extra terms for conserved quantity
double e_scale = compute_etotal() + compute_scalar();
// propagate particle velocities 1/2 step
if (dftb) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( int k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
TS_term = TS_dot/(mass[type[i]]*velocity_sum);
escale_term = force->ftm2v*beta*(e0-e_scale) /
(mass[type[i]]*velocity_sum);
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
D += escale_term - TS_term;
if ( k == direction ) D -= 2.0 * omega[sd] / vol;
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for ( int k = 0; k < 3; k++ ) {
double C = f[i][k] * force->ftm2v / mass[type[i]];
double D = mu * omega[sd] * omega[sd] /
(velocity_sum * mass[type[i]] * vol );
if ( k == direction ) {
D = D - 2.0 * omega[sd] / vol;
}
if ( fabs(dthalf * D) > 1.0e-06 ) {
double expd = exp(D * dthalf);
v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D;
} else {
v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf +
0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf;
}
}
}
}
}
// compute new pressure and volume.
// compute new pressure and volume
temperature->compute_vector();
@ -604,7 +749,7 @@ void FixMSST::final_integrate()
velocity_sum = compute_vsum();
vol = compute_vol();
// propagate the time derivative of the volume 1/2 step at fixed V, r, rdot.
// propagate the time derivative of the volume 1/2 step at fixed V, r, rdot
p_msst = nktv2p * mvv2e * velocity * velocity * total_mass *
( v0 - vol )/( v0 * v0 );
@ -612,11 +757,9 @@ void FixMSST::final_integrate()
( qmass * nktv2p * mvv2e );
double B = total_mass * mu / ( qmass * vol );
// prevent blow-up of the volume.
// prevent blow-up of the volume
if ( vol > v0 && A > 0.0 ) {
A = -A;
}
if ( vol > v0 && A > 0.0 ) A = -A;
// use taylor expansion to avoid singularity at B == 0.
@ -632,8 +775,9 @@ void FixMSST::final_integrate()
lagrangian_position -= velocity*vol/v0*update->dt;
// trigger virial computation on next timestep
// trigger energy and virial computation on next timestep
pe->addstep(update->ntimestep+1);
pressure->addstep(update->ntimestep+1);
}
@ -707,11 +851,12 @@ void FixMSST::remap(int flag)
void FixMSST::write_restart(FILE *fp)
{
int n = 0;
double list[4];
double list[5];
list[n++] = omega[direction];
list[n++] = e0;
list[n++] = v0;
list[n++] = p0;
list[n++] = TS_int;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
@ -730,7 +875,9 @@ void FixMSST::restart(char *buf)
omega[direction] = list[n++];
e0 = list[n++];
v0 = list[n++];
p0 = list[n++];
p0 = list[n++];
TS_int = list[n++];
tscale = 0.0; // set tscale to zero for restart
p0_set = 1;
v0_set = 1;
e0_set = 1;
@ -752,11 +899,13 @@ int FixMSST::modify_param(int narg, char **arg)
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID");
if (icompute < 0)
error->all(FLERR,"Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,"Fix_modify temperature ID does not compute temperature");
error->all(FLERR,"Fix_modify temperature ID does not "
"compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning(FLERR,"Temperature for MSST is not for group all");
@ -806,6 +955,10 @@ double FixMSST::compute_scalar()
(1.0 - volume/ v0) * mvv2e;
energy -= p0 * ( v0 - volume ) / nktv2p;
// subtract off precomputed TS_int integral value
energy -= TS_int;
return energy;
}
@ -920,25 +1073,7 @@ double FixMSST::compute_vol()
return domain->xprd * domain->yprd;
}
/* ----------------------------------------------------------------------
Checks to see if the allocated size of old_velocity is >= n
The number of local atoms can change during a parallel run.
------------------------------------------------------------------------- */
void FixMSST::check_alloc(int n)
{
if ( atoms_allocated < n ) {
for ( int j = 0; j < atoms_allocated; j++ ) {
delete [] old_velocity[j];
}
delete [] old_velocity;
old_velocity = new double* [n];
for ( int j = 0; j < n; j++ )
old_velocity[j] = new double [3];
atoms_allocated = n;
}
}
/* ---------------------------------------------------------------------- */
double FixMSST::compute_vsum()
{
@ -959,3 +1094,14 @@ double FixMSST::compute_vsum()
MPI_Allreduce(&t,&vsum,1,MPI_DOUBLE,MPI_SUM,world);
return vsum;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixMSST::memory_usage()
{
double bytes = 3*atom->nmax * sizeof(double);
return bytes;
}

53
src/SHOCK/fix_msst.h
View File

@ -10,10 +10,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* Implementation of the Multi-Scale Shock Method.
See Reed, Fried, Joannopoulos, Phys. Rev. Lett., 90, 235503(2003).
Implementation by Laurence Fried, LLNL, 4/2007.
*/
#ifdef FIX_CLASS
@ -42,55 +38,68 @@ class FixMSST : public Fix {
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
double memory_usage();
private:
double dtv,dtf,dthalf; // Full and half step sizes.
double boltz,nktv2p, mvv2e; // Boltzmann factor and unit conversions.
double total_mass; // Mass of the computational cell.
double dtv,dtf,dthalf; // Full and half step sizes
double boltz,nktv2p, mvv2e; // Boltzmann factor and unit conversions
double total_mass; // Mass of the computational cell
double omega[3]; // Time derivative of the volume.
double omega[3]; // Time derivative of the volume
double p_current[3],dilation[3];
double qmass; // Effective cell mass.
double mu; // Effective cell viscosity.
double qmass; // Effective cell mass
double mu; // Effective cell viscosity
double tscale; // Converts thermal energy to compressive
// strain ke at simulation start
int dftb; // flag for use with DFTB+
double velocity_sum; // Sum of the velocities squared.
double velocity_sum; // Sum of the velocities squared
double damping; // Damping function for TS force term at
// small volume difference (v0 - vol)
double T0S0; // Initial TS term for DFTB+ simulations
double S_elec,S_elec_1,S_elec_2; // time history of electron entropy
// for DFTB+ simulaitons
double TS_dot; // time derivative of TS term for
// DFTB+ simulations
double **old_velocity; // Saved velocities.
double **old_velocity; // Saved velocities
int kspace_flag; // 1 if KSpace invoked, 0 if not
int nrigid; // number of rigid fixes
int *rfix; // indices of rigid fixes
char *id_temp,*id_press; // Strings with identifiers of
char *id_pe; // created computes.
char *id_pe; // created computes
class Compute *temperature; // Computes created to evaluate
class Compute *pressure; // thermodynamic quantities.
class Compute *pressure; // thermodynamic quantities
class Compute *pe;
int tflag,pflag,vsflag,peflag; // Flags to keep track of computes that
// were created.
// were created
// shock initial conditions.
// shock initial conditions
double e0; // Initial energy
double v0; // Initial volume
double p0; // Initial pressure
double velocity; // Velocity of the shock.
double velocity; // Velocity of the shock
double lagrangian_position; // Lagrangian location of computational cell
int direction; // Direction of shock
int p0_set; // Is pressure set.
int v0_set; // Is volume set.
int e0_set; // Is energy set.
int p0_set; // Is pressure set
int v0_set; // Is volume set
int e0_set; // Is energy set
double TS_int; // Needed for conserved quantity
// with thermal electronic excitations
double beta; // Energy conservation scaling factor
int atoms_allocated; // The number of allocated atoms in old_velocity.
int maxold; // allocated size of old_velocity
double TS_dftb; // value needed from DFTB+ via fix external
class FixExternal *fix_external; // ptr to fix external
// functions
void couple();
void remap(int);
void check_alloc(int n);
double compute_etotal();
double compute_vol();
double compute_hugoniot();

55
src/fix_external.cpp
View File

@ -30,7 +30,7 @@ enum{PF_CALLBACK,PF_ARRAY};
FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
fexternal(NULL)
fexternal(NULL), caller_vector(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix external command");
@ -62,6 +62,11 @@ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
user_energy = 0.0;
// optional vector of values provided by caller
// vector_flag and size_vector are setup via set_vector_length()
caller_vector = NULL;
}
/* ---------------------------------------------------------------------- */
@ -73,6 +78,7 @@ FixExternal::~FixExternal()
atom->delete_callback(id,0);
memory->destroy(fexternal);
delete [] caller_vector;
}
/* ---------------------------------------------------------------------- */
@ -167,10 +173,15 @@ void FixExternal::min_post_force(int vflag)
post_force(vflag);
}
// ----------------------------------------------------------------------
// "set" methods caller can invoke directly
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
caller invokes this method to set its contribution to global energy
do not just return if eflag_global is not set
input script could access this quantity via compute_scalar()
unlike other energy/virial set methods:
do not just return if eflag_global is not set
b/c input script could access this quantity via compute_scalar()
even if eflag is not set on a particular timestep
------------------------------------------------------------------------- */
@ -220,6 +231,34 @@ void FixExternal::set_virial_peratom(double **caller_virial)
vatom[i][j] = caller_virial[i][j];
}
/* ----------------------------------------------------------------------
caller invokes this method to set length of vector of values
assume all vector values are extensive, could make this an option
------------------------------------------------------------------------- */
void FixExternal::set_vector_length(int n)
{
delete [] caller_vector;
vector_flag = 1;
size_vector = n;
extvector = 1;
caller_vector = new double[n];
}
/* ----------------------------------------------------------------------
caller invokes this method to set Index value in vector
index ranges from 1 to N inclusive
------------------------------------------------------------------------- */
void FixExternal::set_vector(int index, double value)
{
if (index >= size_vector)
error->all(FLERR,"Invalid set_vector index in fix external");
caller_vector[index-1] = value;
}
/* ----------------------------------------------------------------------
potential energy of added force
up to user to set it via set_energy()
@ -230,6 +269,16 @@ double FixExternal::compute_scalar()
return user_energy;
}
/* ----------------------------------------------------------------------
arbitrary value computed by caller
up to user to set it via set_vector()
------------------------------------------------------------------------- */
double FixExternal::compute_vector(int n)
{
return caller_vector[n];
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */

4
src/fix_external.h
View File

@ -39,11 +39,14 @@ class FixExternal : public Fix {
void post_force(int);
void min_post_force(int);
double compute_scalar();
double compute_vector(int);
void set_energy_global(double);
void set_virial_global(double *);
void set_energy_peratom(double *);
void set_virial_peratom(double **);
void set_vector_length(int);
void set_vector(int,double);
double memory_usage();
void grow_arrays(int);
@ -59,6 +62,7 @@ class FixExternal : public Fix {
FnPtr callback;
void *ptr_caller;
double user_energy;
double *caller_vector;
};
}

1
src/modify.cpp
View File

@ -1066,7 +1066,6 @@ int Modify::check_rigid_region_overlap(int groupbit, Region *reg)
int Modify::check_rigid_list_overlap(int *select)
{
const int * const mask = atom->mask;
const int nlocal = atom->nlocal;
int dim;