fix_msst.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jul 17, 06:44

fix_msst.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 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)
	------------------------------------------------------------------------- */

	#include <string.h>
	#include <stdlib.h>
	#include <math.h>
	#include "fix_msst.h"
	#include "atom.h"
	#include "force.h"
	#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 "thermo.h"
	#include "memory.h"
	#include "error.h"

	using namespace LAMMPS_NS;
	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)
	{
	if (narg < 4) error->all(FLERR,"Illegal fix msst command");

	restart_global = 1;
	box_change_size = 1;
	time_integrate = 1;
	scalar_flag = 1;
	vector_flag = 1;
	size_vector = 4;
	global_freq = 1;
	extscalar = 1;
	extvector = 0;

	// set defaults

	velocity = 0.0;
	dilation[0] = dilation[1] = dilation[2] = 1.0;
	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;
	direction = 2;
	p0_set = 0;
	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");

	velocity = force->numeric(FLERR,arg[4]);
	if (velocity < 0) error->all(FLERR,"Illegal fix msst command");

	// 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 += 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 += 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]);
	p0_set = 1;
	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 += 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 += 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 += 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");
	if (direction == 0) fprintf(screen," Shock in x direction\n");
	else if (direction == 1) fprintf(screen," Shock in y direction\n");
	else if (direction == 2) fprintf(screen," Shock in z direction\n");
	fprintf(screen," Cell mass-like parameter qmass "
	"(units of mass^2/length^4) = %12.5e\n", qmass);
	fprintf(screen," Shock velocity = %12.5e\n", velocity);
	fprintf(screen," Artificial viscosity "
	"(units of mass/length/time) = %12.5e\n", mu);

	if (p0_set)
	fprintf(screen," Initial pressure specified to be %12.5e\n", p0);
	else fprintf(screen," Initial pressure calculated on first step\n");

	if (v0_set)
	fprintf(screen," Initial volume specified to be %12.5e\n", v0);
	else fprintf(screen," Initial volume calculated on first step\n");

	if (e0_set)
	fprintf(screen," Initial energy specified to be %12.5e\n", e0);
	else fprintf(screen," Initial energy calculated on first step\n");
	}
	if (logfile) {
	fprintf(logfile,"MSST parameters:\n");
	if (direction == 0) fprintf(logfile," Shock in x direction\n");
	else if (direction == 1) fprintf(logfile," Shock in y direction\n");
	else if (direction == 2) fprintf(logfile," Shock in z direction\n");
	fprintf(logfile," Cell mass-like parameter qmass "
	"(units of mass^2/length^4) = %12.5e\n", qmass);
	fprintf(logfile," Shock velocity = %12.5e\n", velocity);
	fprintf(logfile," Artificial viscosity "
	"(units of mass/length/time) = %12.5e\n", mu);

	if (p0_set)
	fprintf(logfile," Initial pressure specified to be %12.5e\n", p0);
	else fprintf(logfile," Initial pressure calculated on first step\n");

	if (v0_set)
	fprintf(logfile," Initial volume specified to be %12.5e\n", v0);
	else fprintf(logfile," Initial volume calculated on first step\n");

	if (e0_set)
	fprintf(logfile," Initial energy specified to be %12.5e\n", e0);
	else fprintf(logfile," Initial energy calculated on first step\n");
	}
	}

	// check for periodicity in controlled dimensions

	if (domain->nonperiodic) error->all(FLERR,"Fix msst requires a periodic box");

	// 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) + 10;
	id_temp = new char[n];
	strcpy(id_temp,id);
	strcat(id_temp,"MSST_temp");

	char *newarg = new char[3];
	newarg[0] = id_temp;
	newarg[1] = (char *) "all";
	newarg[2] = (char *) "temp";
	modify->add_compute(3,newarg);
	delete [] newarg;
	tflag = 1;

	// 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) + 11;
	id_press = new char[n];
	strcpy(id_press,id);
	strcat(id_press,"MSST_press");

	newarg = new char*[4];
	newarg[0] = id_press;
	newarg[1] = (char *) "all";
	newarg[2] = (char *) "pressure";
	newarg[3] = id_temp;
	modify->add_compute(4,newarg);
	delete [] newarg;
	pflag = 1;

	// create a new potential energy compute
	// id = fix-ID + "MSST_pe", compute group = all

	n = strlen(id) + 8;
	id_pe = new char[n];
	strcpy(id_pe,id);
	strcat(id_pe,"MSST_pe");

	newarg = new char*[3];
	newarg[0] = id_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;
	nrigid = 0;
	rfix = NULL;

	maxold = -1;
	old_velocity = NULL;
	}

	/* ---------------------------------------------------------------------- */

	FixMSST::~FixMSST()
	{
	delete [] rfix;

	// delete temperature and pressure if fix created them

	if (tflag) modify->delete_compute(id_temp);
	if (pflag) modify->delete_compute(id_press);
	if (peflag) modify->delete_compute(id_pe);

	delete [] id_temp;
	delete [] id_press;
	delete [] id_pe;

	memory->destroy(old_velocity);
	}

	/* ---------------------------------------------------------------------- */

	int FixMSST::setmask()
	{
	int mask = 0;
	mask \|= INITIAL_INTEGRATE;
	mask \|= FINAL_INTEGRATE;
	mask \|= THERMO_ENERGY;
	return mask;
	}

	/* ---------------------------------------------------------------------- */

	void FixMSST::init()
	{
	if (atom->mass == NULL)
	error->all(FLERR,"Cannot use fix msst without per-type mass defined");

	// set compute ptrs

	int itemp = modify->find_compute(id_temp);
	int ipress = modify->find_compute(id_press);
	int ipe = modify->find_compute(id_pe);
	if (itemp < 0 \|\| ipress < 0\|\| ipe < 0)
	error->all(FLERR,"Could not find fix msst compute ID");
	if (modify->compute[itemp]->tempflag == 0)
	error->all(FLERR,"Fix msst compute ID does not compute temperature");
	if (modify->compute[ipress]->pressflag == 0)
	error->all(FLERR,"Fix msst compute ID does not compute pressure");
	if (modify->compute[ipe]->peflag == 0)
	error->all(FLERR,"Fix msst compute ID does not compute potential energy");

	temperature = modify->compute[itemp];
	pressure = modify->compute[ipress];
	pe = modify->compute[ipe];

	dtv = update->dt;
	dtf = 0.5 * update->dt * force->ftm2v;
	dthalf = 0.5 * update->dt;

	boltz = force->boltz;
	nktv2p = force->nktv2p;
	mvv2e = force->mvv2e;

	double mass = 0.0;
	for (int i = 0; i < atom->nlocal; i++) mass += atom->mass[atom->type[i]];
	MPI_Allreduce(&mass,&total_mass,1,MPI_DOUBLE,MPI_SUM,world);

	if (force->kspace) kspace_flag = 1;
	else kspace_flag = 0;

	// detect if any fix rigid exist so rigid bodies move when box is dilated
	// rfix[] = indices to each fix rigid

	delete [] rfix;
	nrigid = 0;
	rfix = NULL;

	for (int i = 0; i < modify->nfix; i++)
	if (strcmp(modify->fix[i]->style,"rigid") == 0 \|\|
	strcmp(modify->fix[i]->style,"poems") == 0) nrigid++;
	if (nrigid) {
	rfix = new int[nrigid];
	nrigid = 0;
	for (int i = 0; i < modify->nfix; i++)
	if (strcmp(modify->fix[i]->style,"rigid") == 0 \|\|
	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");
	}
	}

	/* ----------------------------------------------------------------------
	compute T,P before integrator starts
	------------------------------------------------------------------------- */

	void FixMSST::setup(int vflag)
	{
	lagrangian_position = 0.0;

	temperature->compute_vector();
	pressure->compute_vector();
	couple();
	velocity_sum = compute_vsum();

	if ( v0_set == 0 ) {
	v0 = compute_vol();
	v0_set = 1;
	if (comm->me == 0) {
	if ( screen ) fprintf(screen,"Fix MSST v0 = %12.5e\n", v0);
	if ( logfile ) fprintf(logfile,"Fix MSST v0 = %12.5e\n", v0);
	}
	}

	if ( p0_set == 0 ) {
	p0 = p_current[direction];
	p0_set = 1;

	if ( comm->me == 0 ) {
	if ( screen ) fprintf(screen,"Fix MSST p0 = %12.5e\n", p0);
	if ( logfile ) fprintf(logfile,"Fix MSST p0 = %12.5e\n", p0);
	}
	}

	if ( e0_set == 0 ) {
	e0 = compute_etotal();
	e0_set = 1;

	if ( comm->me == 0 ) {
	if ( screen ) fprintf(screen,"Fix MSST e0 = to be %12.5e\n",e0);
	if ( logfile ) fprintf(logfile,"Fix MSST e0 = to be %12.5e\n",e0);
	}

	}

	temperature->compute_vector();
	double *ke_tensor = temperature->vector;
	double ke_temp = ke_tensor[0]+ke_tensor[1]+ke_tensor[2];
	if (ke_temp > 0.0 && tscale > 0.0 ) {

	// transfer energy from atom velocities to cell volume motion
	// to bias initial compression

	double **v = atom->v;
	int *mask = atom->mask;
	double sqrt_initial_temperature_scaling = sqrt(1.0-tscale);

	double fac1 = tscaletotal_mass/qmasske_temp/force->mvv2e;

	omega[direction]=-1*sqrt(fac1);
	double fac2 = omega[direction]/v0;

	if ( comm->me == 0 && tscale != 1.0) {
	if ( screen )
	fprintf(screen,"Fix MSST initial strain rate of %12.5e established "
	"by reducing temperature by factor of %12.5e\n",
	fac2,tscale);
	if ( logfile )
	fprintf(logfile,"Fix MSST initial strain rate of %12.5e established "
	"by reducing temperature by factor of %12.5e\n",
	fac2,tscale);
	}
	for (int i = 0; i < atom->nlocal; i++) {
	if (mask[i] & groupbit) {
	for (int k = 0; k < 3; k++ ) {
	v[i][k]*=sqrt_initial_temperature_scaling;
	}
	}
	}
	}

	// trigger virial computation on next timestep

	pe->addstep(update->ntimestep+1);
	pressure->addstep(update->ntimestep+1);
	}

	/* ----------------------------------------------------------------------
	1st half of Verlet update
	------------------------------------------------------------------------- */

	void FixMSST::initial_integrate(int vflag)
	{
	int i,k;
	double p_msst; // MSST driving pressure
	double vol;

	int nlocal = atom->nlocal;
	int *mask = atom->mask;
	double **v = atom->v;
	double **f = atom->f;
	double *mass = atom->mass;
	int *type = atom->type;
	double **x = atom->x;
	int sd = direction;

	// realloc old_velocity if necessary

	if (nlocal > maxold) {
	memory->destroy(old_velocity);
	maxold = atom->nmax;
	memory->create(old_velocity,maxold,3,"msst:old_velocity");
	}

	// for DFTB, extract TS_dftb from fix external
	// must convert energy to mv^2 units

	if (dftb) {
	const double TS_dftb = fix_external->compute_vector(0);
	const double TS = force->ftm2v*TS_dftb;
	// update S_elec terms and compute TS_dot via finite differences
	S_elec_2 = S_elec_1;
	S_elec_1 = S_elec;
	const double Temp = temperature->compute_scalar();
	S_elec = TS/Temp;
	TS_dot = Temp(3.0S_elec-4.0S_elec_1+S_elec_2)/(2.0update->dt);
	TS_int += (update->dt*TS_dot);
	if (update->ntimestep == 1) T0S0 = TS;
	}

	// compute new pressure and volume

	temperature->compute_vector();
	pressure->compute_vector();
	couple();
	vol = compute_vol();

	// 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

	p_msst = nktv2p * mvv2e * velocity * velocity * total_mass *
	( v0 - vol)/( v0 * v0);
	double A = total_mass * ( p_current[sd] - p0 - p_msst ) /
	(qmass * nktv2p * mvv2e);
	double B = total_mass * mu / ( qmass * vol );

	// prevent blow-up of the volume

	if (vol > v0 && A > 0.0) A = -A;

	// 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 )
	* exp(-B * dthalf);
	} else {
	omega[sd] = omega[sd] + (A - B * omega[sd]) * dthalf +
	0.5 * (B * B * omega[sd] - A * B ) * dthalf * dthalf;
	}

	// propagate velocity sum 1/2 step by
	// temporarily propagating the velocities

	velocity_sum = compute_vsum();

	if (dftb) {
	for (i = 0; i < nlocal; i++) {
	if (mask[i] & groupbit) {
	for ( k = 0; k < 3; k++ ) {
	const double C = f[i][k] * force->ftm2v / mass[type[i]];
	const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
	const double escale_term = force->ftm2vbeta(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 ) {
	const 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++ ) {
	const 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 -= 2.0 * omega[sd] / vol;
	}
	if ( fabs(dthalf * D) > 1.0e-06 ) {
	const 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

	for (i = 0; i < nlocal; i++) {
	if (mask[i] & groupbit) {
	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

	if (dftb) {
	for (i = 0; i < nlocal; i++) {
	if (mask[i] & groupbit) {
	for ( k = 0; k < 3; k++ ) {
	const double C = f[i][k] * force->ftm2v / mass[type[i]];
	const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
	const double escale_term = force->ftm2vbeta(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 ) {
	const 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++ ) {
	const 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 -= 2.0 * omega[sd] / vol;
	}
	if ( fabs(dthalf * D) > 1.0e-06 ) {
	const 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

	double vol1 = vol + omega[sd] * dthalf;

	// rescale positions and change box size

	dilation[sd] = vol1/vol;
	remap(0);

	// propagate particle positions 1 time step

	for (i = 0; i < nlocal; i++) {
	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];
	}
	}

	// propagate the volume 1/2 step

	double vol2 = vol1 + omega[sd] * dthalf;

	// rescale positions and change box size

	dilation[sd] = vol2/vol1;
	remap(0);

	if (kspace_flag) force->kspace->setup();
	}

	/* ----------------------------------------------------------------------
	2nd half of Verlet update
	------------------------------------------------------------------------- */

	void FixMSST::final_integrate()
	{
	int i;
	double p_msst; // MSST driving pressure

	// v update only for atoms in MSST group

	double **v = atom->v;
	double **f = atom->f;
	double *mass = atom->mass;
	int *type = atom->type;
	int *mask = atom->mask;
	int nlocal = atom->nlocal;
	double vol = compute_vol();

	int sd = direction;

	// compute etot + extra terms for conserved quantity

	double e_scale = compute_etotal() + compute_scalar();

	// for DFTB, extract TS_dftb from fix external
	// must convert energy to mv^2 units

	if (dftb) {
	const double TS_dftb = fix_external->compute_vector(0);
	const double TS = force->ftm2v*TS_dftb;
	S_elec_2 = S_elec_1;
	S_elec_1 = S_elec;
	const double Temp = temperature->compute_scalar();
	// update S_elec terms and compute TS_dot via finite differences
	S_elec = TS/Temp;
	TS_dot = Temp(3.0S_elec-4.0S_elec_1+S_elec_2)/(2.0update->dt);
	TS_int += (update->dt*TS_dot);
	if (update->ntimestep == 1) T0S0 = TS;
	}

	// 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++ ) {
	const double C = f[i][k] * force->ftm2v / mass[type[i]];
	const double TS_term = TS_dot/(mass[type[i]]*velocity_sum);
	const double escale_term = force->ftm2vbeta(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 ) {
	const 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++ ) {
	const 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 -= 2.0 * omega[sd] / vol;
	}
	if ( fabs(dthalf * D) > 1.0e-06 ) {
	const 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

	temperature->compute_vector();

	pressure->compute_vector();
	couple();
	velocity_sum = compute_vsum();
	vol = compute_vol();

	// 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 );
	double A = total_mass * ( p_current[sd] - p0 - p_msst ) /
	( qmass * nktv2p * mvv2e );
	const double B = total_mass * mu / ( qmass * vol );

	// prevent blow-up of the volume

	if ( vol > v0 && A > 0.0 ) A = -A;

	// 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 ) * exp(-B * dthalf);
	} else {
	omega[sd] = omega[sd] + (A - B * omega[sd]) * dthalf +
	0.5 * (B * B * omega[sd] - A * B ) * dthalf * dthalf;
	}

	// calculate Lagrangian position of computational cell

	lagrangian_position -= velocityvol/v0update->dt;

	// trigger energy and virial computation on next timestep

	pe->addstep(update->ntimestep+1);
	pressure->addstep(update->ntimestep+1);
	}

	/* ---------------------------------------------------------------------- */

	void FixMSST::couple()
	{
	double *tensor = pressure->vector;

	p_current[0] = tensor[0];
	p_current[1] = tensor[1];
	p_current[2] = tensor[2];
	}

	/* ----------------------------------------------------------------------
	change box size
	remap owned or owned+ghost atoms depending on flag
	if rigid bodies exist, scale rigid body centers-of-mass
	------------------------------------------------------------------------- */

	void FixMSST::remap(int flag)
	{
	int i,n;
	double oldlo,oldhi,ctr;

	double **v = atom->v;
	if (flag) n = atom->nlocal + atom->nghost;
	else n = atom->nlocal;

	// convert pertinent atoms and rigid bodies to lamda coords

	domain->x2lamda(n);

	if (nrigid)
	for (i = 0; i < nrigid; i++)
	modify->fix[rfix[i]]->deform(0);

	// reset global and local box to new size/shape

	for (i = 0; i < 3; i++) {
	if ( direction == i ) {
	oldlo = domain->boxlo[i];
	oldhi = domain->boxhi[i];
	ctr = 0.5 * (oldlo + oldhi);
	domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
	domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
	}
	}

	domain->set_global_box();
	domain->set_local_box();

	// convert pertinent atoms and rigid bodies back to box coords

	domain->lamda2x(n);

	if (nrigid)
	for (i = 0; i < nrigid; i++)
	modify->fix[rfix[i]]->deform(1);

	for (i = 0; i < n; i++) {
	v[i][direction] = v[i][direction] *
	dilation[direction];
	}
	}

	/* ----------------------------------------------------------------------
	pack entire state of Fix into one write
	------------------------------------------------------------------------- */

	void FixMSST::write_restart(FILE *fp)
	{
	int n = 0;
	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);
	fwrite(&list,sizeof(double),n,fp);
	}
	}

	/* ----------------------------------------------------------------------
	use state info from restart file to restart the Fix
	------------------------------------------------------------------------- */

	void FixMSST::restart(char *buf)
	{
	int n = 0;
	double list = (double ) buf;
	omega[direction] = list[n++];
	e0 = list[n++];
	v0 = 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;
	}

	/* ---------------------------------------------------------------------- */

	int FixMSST::modify_param(int narg, char **arg)
	{
	if (strcmp(arg[0],"temp") == 0) {
	if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
	if (tflag) {
	modify->delete_compute(id_temp);
	tflag = 0;
	}
	delete [] id_temp;
	int n = strlen(arg[1]) + 1;
	id_temp = new char[n];
	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");
	temperature = modify->compute[icompute];

	if (temperature->tempflag == 0)
	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");

	return 2;

	} else if (strcmp(arg[0],"press") == 0) {
	if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
	if (pflag) {
	modify->delete_compute(id_press);
	pflag = 0;
	}
	delete [] id_press;
	int n = strlen(arg[1]) + 1;
	id_press = new char[n];
	strcpy(id_press,arg[1]);

	int icompute = modify->find_compute(id_press);
	if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
	pressure = modify->compute[icompute];

	if (pressure->pressflag == 0)
	error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
	return 2;
	}
	return 0;
	}

	/* ---------------------------------------------------------------------- */

	double FixMSST::compute_scalar()
	{
	// compute new pressure and volume

	temperature->compute_vector();
	pressure->compute_vector();
	couple();

	double volume = compute_vol();

	double energy = 0.0;
	int i;

	i = direction;
	energy = qmass * omega[i] * omega[i] / (2.0 * total_mass) * mvv2e;
	energy -= 0.5 * total_mass * velocity * velocity *
	(1.0 - volume/ v0) *
	(1.0 - volume/ v0) * mvv2e;
	energy -= p0 * ( v0 - volume ) / nktv2p;

	// subtract off precomputed TS_int integral value
	// TS_int = 0 for non DFTB calculations

	if (dftb) energy -= TS_int;

	return energy;
	}

	/* ----------------------------------------------------------------------
	return a single element from the following vector,
	[dhug,dray,lgr_vel,lgr_pos]
	------------------------------------------------------------------------- */

	double FixMSST::compute_vector(int n)
	{
	if (n == 0) {
	return compute_hugoniot();
	} else if (n == 1) {
	return compute_rayleigh();
	} else if (n == 2) {
	return compute_lagrangian_speed();
	} else if (n == 3) {
	return compute_lagrangian_position();
	}
	return 0.0;
	}

	/* ----------------------------------------------------------------------
	Computes the deviation of the current point
	from the Hugoniot in Kelvin for the MSST
	------------------------------------------------------------------------- */

	double FixMSST::compute_hugoniot()
	{
	double v, e, p;
	double dhugo;

	e = compute_etotal();

	temperature->compute_vector();
	pressure->compute_vector();
	p = pressure->vector[direction];

	v = compute_vol();

	dhugo = (0.5 * (p + p0 ) * ( v0 - v)) /
	force->nktv2p + e0 - e;
	dhugo /= temperature->dof * force->boltz;

	return dhugo;
	}

	/* ----------------------------------------------------------------------
	Computes the deviation of the current point from the Rayleigh
	in pressure units for the MSST
	------------------------------------------------------------------------- */

	double FixMSST::compute_rayleigh()
	{
	double v, p;
	double drayleigh;

	temperature->compute_vector();
	pressure->compute_vector();
	p = pressure->vector[direction];

	v = compute_vol();

	drayleigh = p - p0 -
	total_mass * velocity * velocity * force->mvv2e *
	(1.0 - v / v0 ) * force->nktv2p / v0;

	return drayleigh;
	}

	/* ----------------------------------------------------------------------
	Computes the speed of the MSST computational cell in the
	unshocked material rest-frame
	------------------------------------------------------------------------- */

	double FixMSST::compute_lagrangian_speed()
	{
	double v = compute_vol();
	return velocity*(1.0-v/v0);
	}

	/* ----------------------------------------------------------------------
	Computes the distance behind the
	shock front of the MSST computational cell.
	------------------------------------------------------------------------- */

	double FixMSST::compute_lagrangian_position()
	{
	return lagrangian_position;
	}

	/* ---------------------------------------------------------------------- */

	double FixMSST::compute_etotal()
	{
	double epot,ekin,etot;
	epot = pe->compute_scalar();
	if (thermo_energy) epot -= compute_scalar();
	ekin = temperature->compute_scalar();
	ekin = 0.5 temperature->dof * force->boltz;
	etot = epot+ekin;
	return etot;
	}

	/* ---------------------------------------------------------------------- */

	double FixMSST::compute_vol()
	{
	if (domain->dimension == 3)
	return domain->xprd * domain->yprd * domain->zprd;
	else
	return domain->xprd * domain->yprd;
	}

	/* ---------------------------------------------------------------------- */

	double FixMSST::compute_vsum()
	{
	double vsum;

	double **v = atom->v;
	int *mask = atom->mask;
	int nlocal = atom->nlocal;

	double t = 0.0;

	for (int i = 0; i < nlocal; i++) {
	if (mask[i] & groupbit) {
	t += (v[i][0]v[i][0] + v[i][1]v[i][1] + v[i][2]*v[i][2]) ;
	}
	}

	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 = 3atom->nmax sizeof(double);
	return bytes;
	}

fix_msst.cppNo OneTemporaryActions

File Metadata

fix_msst.cppView Options

Event Timeline

fix_msst.cpp
No OneTemporary
Actions

fix_msst.cpp
View Options