diff --git a/doc/src/fix_msst.txt b/doc/src/fix_msst.txt index bd1edd805a..43f35d6880 100644 --- a/doc/src/fix_msst.txt +++ b/doc/src/fix_msst.txt @@ -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 diff --git a/src/SHOCK/fix_msst.cpp b/src/SHOCK/fix_msst.cpp index c514a9f086..c710bdb77f 100644 --- a/src/SHOCK/fix_msst.cpp +++ b/src/SHOCK/fix_msst.cpp @@ -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 @@ -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 = 0; + 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(1); + 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(1); + 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; +} + diff --git a/src/SHOCK/fix_msst.h b/src/SHOCK/fix_msst.h index e44c34dde3..450256bcab 100644 --- a/src/SHOCK/fix_msst.h +++ b/src/SHOCK/fix_msst.h @@ -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(); diff --git a/src/fix_external.cpp b/src/fix_external.cpp index 40e538c6c6..85523b912c 100644 --- a/src/fix_external.cpp +++ b/src/fix_external.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/fix_external.h b/src/fix_external.h index b5f6d7ae15..1e9e78365f 100644 --- a/src/fix_external.h +++ b/src/fix_external.h @@ -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; }; }