ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add original TTM for testing

Browse Source
This commit is contained in:
Steve Plimpton
2021-08-19 16:55:57 -06:00
parent 95bae4d78c
commit f0a041799f
2 changed files with 859 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

704
src/EXTRA-FIX/fix_ttm_old.cpp Normal file
View File

@ -0,0 +1,704 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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: Paul Crozier (SNL)
Carolyn Phillips (University of Michigan)
------------------------------------------------------------------------- */
#include "fix_ttm_old.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "domain.h"
#include "respa.h"
#include "comm.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
FixTTMOld::FixTTMOld(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
random(nullptr), fp(nullptr), nsum(nullptr), nsum_all(nullptr),
gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), flangevin(nullptr),
T_electron(nullptr), T_electron_old(nullptr), sum_vsq(nullptr), sum_mass_vsq(nullptr),
sum_vsq_all(nullptr), sum_mass_vsq_all(nullptr), net_energy_transfer(nullptr),
net_energy_transfer_all(nullptr)
{
if (narg < 15) error->all(FLERR,"Illegal fix ttm command");
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extvector = 1;
nevery = 1;
restart_peratom = 1;
restart_global = 1;
seed = utils::inumeric(FLERR,arg[3],false,lmp);
electronic_specific_heat = utils::numeric(FLERR,arg[4],false,lmp);
electronic_density = utils::numeric(FLERR,arg[5],false,lmp);
electronic_thermal_conductivity = utils::numeric(FLERR,arg[6],false,lmp);
gamma_p = utils::numeric(FLERR,arg[7],false,lmp);
gamma_s = utils::numeric(FLERR,arg[8],false,lmp);
v_0 = utils::numeric(FLERR,arg[9],false,lmp);
nxnodes = utils::inumeric(FLERR,arg[10],false,lmp);
nynodes = utils::inumeric(FLERR,arg[11],false,lmp);
nznodes = utils::inumeric(FLERR,arg[12],false,lmp);
nfileevery = utils::inumeric(FLERR,arg[14],false,lmp);
if (nfileevery) {
if (narg != 16) error->all(FLERR,"Illegal fix ttm command");
if (comm->me == 0) {
fp = fopen(arg[15],"w");
if (fp == nullptr)
error->one(FLERR,"Cannot open output file {}: {}",
arg[15], utils::getsyserror());
}
}
// error check
if (seed <= 0)
error->all(FLERR,"Invalid random number seed in fix ttm command");
if (electronic_specific_heat <= 0.0)
error->all(FLERR,"Fix ttm electronic_specific_heat must be > 0.0");
if (electronic_density <= 0.0)
error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
if (electronic_thermal_conductivity < 0.0)
error->all(FLERR,"Fix ttm electronic_thermal_conductivity must be >= 0.0");
if (gamma_p <= 0.0) error->all(FLERR,"Fix ttm gamma_p must be > 0.0");
if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
error->all(FLERR,"Fix ttm number of nodes must be > 0");
v_0_sq = v_0*v_0;
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp,seed + comm->me);
// allocate per-type arrays for force prefactors
gfactor1 = new double[atom->ntypes+1];
gfactor2 = new double[atom->ntypes+1];
// allocate 3d grid variables
// check for allowed maxium number of total grid nodes
total_nnodes = (bigint)nxnodes * (bigint)nynodes * (bigint)nznodes;
if (total_nnodes > MAXSMALLINT)
error->all(FLERR,"Too many nodes in fix ttm");
memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum");
memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all");
memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq");
memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq");
memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all");
memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq_all");
memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old");
memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron");
memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer");
memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer_all");
flangevin = nullptr;
grow_arrays(atom->nmax);
// zero out the flangevin array
for (int i = 0; i < atom->nmax; i++) {
flangevin[i][0] = 0;
flangevin[i][1] = 0;
flangevin[i][2] = 0;
}
atom->add_callback(Atom::GROW);
atom->add_callback(Atom::RESTART);
// set initial electron temperatures from user input file
if (comm->me == 0) read_initial_electron_temperatures(arg[13]);
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
}
/* ---------------------------------------------------------------------- */
FixTTMOld::~FixTTMOld()
{
if (fp) fclose(fp);
delete random;
delete [] gfactor1;
delete [] gfactor2;
memory->destroy(nsum);
memory->destroy(nsum_all);
memory->destroy(sum_vsq);
memory->destroy(sum_mass_vsq);
memory->destroy(sum_vsq_all);
memory->destroy(sum_mass_vsq_all);
memory->destroy(T_electron_old);
memory->destroy(T_electron);
memory->destroy(flangevin);
memory->destroy(net_energy_transfer);
memory->destroy(net_energy_transfer_all);
}
/* ---------------------------------------------------------------------- */
int FixTTMOld::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::init()
{
if (domain->dimension == 2)
error->all(FLERR,"Cannot use fix ttm with 2d simulation");
if (domain->nonperiodic != 0)
error->all(FLERR,"Cannot use non-periodic boundares with fix ttm");
if (domain->triclinic)
error->all(FLERR,"Cannot use fix ttm with triclinic box");
// set force prefactors
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = - gamma_p / force->ftm2v;
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer_all[ixnode][iynode][iznode] = 0;
if (utils::strmatch(update->integrate_style,"^respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::setup(int vflag)
{
if (utils::strmatch(update->integrate_style,"^verlet")) {
post_force_setup(vflag);
} else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa_setup(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::post_force(int /*vflag*/)
{
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double gamma1,gamma2;
// apply damping and thermostat to all atoms in fix group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
if (T_electron[ixnode][iynode][iznode] < 0)
error->all(FLERR,"Electronic temperature dropped below zero");
double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]);
gamma1 = gfactor1[type[i]];
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::post_force_setup(int /*vflag*/)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply langevin forces that have been stored from previous run
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::post_force_respa_setup(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == nlevels_respa-1) post_force_setup(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::reset_dt()
{
for (int i = 1; i <= atom->ntypes; i++)
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTMOld::read_initial_electron_temperatures(const char *filename)
{
int ***T_initial_set;
memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
memset(&T_initial_set[0][0][0],0,total_nnodes*sizeof(int));
std::string name = utils::get_potential_file_path(filename);
if (name.empty())
error->one(FLERR,"Cannot open input file: {}",
filename);
FILE *fpr = fopen(name.c_str(),"r");
// read initial electron temperature values from file
char line[MAXLINE];
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == nullptr) break;
ValueTokenizer values(line);
if (values.has_next()) ixnode = values.next_int();
if (values.has_next()) iynode = values.next_int();
if (values.has_next()) iznode = values.next_int();
if (values.has_next()) T_tmp = values.next_double();
else error->one(FLERR,"Incorrect format in fix ttm input file");
// check correctness of input data
if ((ixnode < 0) || (ixnode >= nxnodes)
|| (iynode < 0) || (iynode >= nynodes)
|| (iznode < 0) || (iznode >= nznodes))
error->one(FLERR,"Fix ttm invalide node index in fix ttm input");
if (T_tmp < 0.0)
error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
fclose(fpr);
// check completeness of input data
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm");
memory->destroy(T_initial_set);
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::end_of_step()
{
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer[ixnode][iynode][iznode] = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
net_energy_transfer[ixnode][iynode][iznode] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
}
MPI_Allreduce(&net_energy_transfer[0][0][0],
&net_energy_transfer_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
// num_inner_timesteps = # of inner steps (thermal solves)
// required this MD step to maintain a stable explicit solve
int num_inner_timesteps = 1;
double inner_dt = update->dt;
double stability_criterion = 1.0 -
2.0*inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
if (stability_criterion < 0.0) {
inner_dt = 0.5*(electronic_specific_heat*electronic_density) /
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning(FLERR,"Too many inner timesteps in fix ttm");
}
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron_old[ixnode][iynode][iznode] =
T_electron[ixnode][iynode][iznode];
// compute new electron T profile
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
int right_xnode = ixnode + 1;
int right_ynode = iynode + 1;
int right_znode = iznode + 1;
if (right_xnode == nxnodes) right_xnode = 0;
if (right_ynode == nynodes) right_ynode = 0;
if (right_znode == nznodes) right_znode = 0;
int left_xnode = ixnode - 1;
int left_ynode = iynode - 1;
int left_znode = iznode - 1;
if (left_xnode == -1) left_xnode = nxnodes - 1;
if (left_ynode == -1) left_ynode = nynodes - 1;
if (left_znode == -1) left_znode = nznodes - 1;
T_electron[ixnode][iynode][iznode] =
T_electron_old[ixnode][iynode][iznode] +
inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity *
((T_electron_old[right_xnode][iynode][iznode] +
T_electron_old[left_xnode][iynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dx/dx +
(T_electron_old[ixnode][right_ynode][iznode] +
T_electron_old[ixnode][left_ynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dy/dy +
(T_electron_old[ixnode][iynode][right_znode] +
T_electron_old[ixnode][iynode][left_znode] -
2*T_electron_old[ixnode][iynode][iznode])/dz/dz) -
(net_energy_transfer_all[ixnode][iynode][iznode])/del_vol);
}
}
// output nodal temperatures for current timestep
if ((nfileevery) && !(update->ntimestep % nfileevery)) {
// compute atomic Ta for each grid point
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
nsum[ixnode][iynode][iznode] = 0;
nsum_all[ixnode][iynode][iznode] = 0;
sum_vsq[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq[ixnode][iynode][iznode] = 0.0;
sum_vsq_all[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0;
}
double massone;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
nsum[ixnode][iynode][iznode] += 1;
sum_vsq[ixnode][iynode][iznode] += vsq;
sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq;
}
MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes,
MPI_INT,MPI_SUM,world);
MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
if (comm->me == 0) {
fmt::print(fp,"{}",update->ntimestep);
double T_a;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
T_a = 0;
if (nsum_all[ixnode][iynode][iznode] > 0)
T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/
(3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e);
fmt::print(fp," {}",T_a);
}
fputs("\t",fp);
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
fmt::print(fp," {}",T_electron[ixnode][iynode][iznode]);
fputs("\n",fp);
}
}
}
/* ----------------------------------------------------------------------
memory usage of 3d grid
------------------------------------------------------------------------- */
double FixTTMOld::memory_usage()
{
double bytes = 0.0;
bytes += (double)5*total_nnodes * sizeof(int);
bytes += (double)14*total_nnodes * sizeof(double);
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixTTMOld::grow_arrays(int ngrow)
{
memory->grow(flangevin,ngrow,3,"TTM:flangevin");
}
/* ----------------------------------------------------------------------
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
double FixTTMOld::compute_vector(int n)
{
double e_energy = 0.0;
double transfer_energy = 0.0;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
e_energy +=
T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy +=
net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
}
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTTMOld::write_restart(FILE *fp)
{
double *rlist;
memory->create(rlist,nxnodes*nynodes*nznodes+1,"TTM:rlist");
int n = 0;
rlist[n++] = seed;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
rlist[n++] = T_electron[ixnode][iynode][iznode];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(rlist,sizeof(double),n,fp);
}
memory->destroy(rlist);
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTTMOld::restart(char *buf)
{
int n = 0;
double *rlist = (double *) buf;
// the seed must be changed from the initial seed
seed = static_cast<int> (0.5*rlist[n++]);
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron[ixnode][iynode][iznode] = rlist[n++];
delete random;
random = new RanMars(lmp,seed+comm->me);
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixTTMOld::pack_restart(int i, double *buf)
{
// pack buf[0] this way because other fixes unpack it
buf[0] = 4;
buf[1] = flangevin[i][0];
buf[2] = flangevin[i][1];
buf[3] = flangevin[i][2];
return 4;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixTTMOld::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
// unpack the Nth first values this way because other fixes pack them
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
flangevin[nlocal][0] = extra[nlocal][m++];
flangevin[nlocal][1] = extra[nlocal][m++];
flangevin[nlocal][2] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixTTMOld::maxsize_restart()
{
return 4;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixTTMOld::size_restart(int /*nlocal*/)
{
return 4;
}

155
src/EXTRA-FIX/fix_ttm_old.h Normal file
View File

@ -0,0 +1,155 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(ttm/old,FixTTMOld);
// clang-format on
#else
#ifndef LMP_FIX_TTM_OLD_H
#define LMP_FIX_TTM_OLD_H
#include "fix.h"
namespace LAMMPS_NS {
class FixTTMOld : public Fix {
public:
FixTTMOld(class LAMMPS *, int, char **);
~FixTTMOld();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void post_force_setup(int);
void post_force_respa_setup(int, int, int);
void end_of_step();
void reset_dt();
void write_restart(FILE *);
void restart(char *);
int pack_restart(int, double *);
void unpack_restart(int, int);
int size_restart(int);
int maxsize_restart();
double memory_usage();
void grow_arrays(int);
double compute_vector(int);
private:
int nfileevery;
int nlevels_respa;
int seed;
class RanMars *random;
FILE *fp;
int nxnodes, nynodes, nznodes;
bigint total_nnodes;
int ***nsum, ***nsum_all;
double *gfactor1, *gfactor2, *ratio, **flangevin;
double ***T_electron, ***T_electron_old;
double ***sum_vsq, ***sum_mass_vsq;
double ***sum_vsq_all, ***sum_mass_vsq_all;
double ***net_energy_transfer, ***net_energy_transfer_all;
double electronic_specific_heat, electronic_density;
double electronic_thermal_conductivity;
double gamma_p, gamma_s, v_0, v_0_sq;
void read_initial_electron_temperatures(const char *);
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
correct. If the file is a compressed file, also check that the gzip
executable can be found and run.
E: Cannot open fix ttm file %s
The output file for the fix ttm command cannot be opened. Check that
the path and name are correct.
E: Invalid random number seed in fix ttm command
Random number seed must be > 0.
E: Fix ttm electronic_specific_heat must be > 0.0
Self-explanatory.
E: Fix ttm electronic_density must be > 0.0
Self-explanatory.
E: Fix ttm electronic_thermal_conductivity must be >= 0.0
Self-explanatory.
E: Fix ttm gamma_p must be > 0.0
Self-explanatory.
E: Fix ttm gamma_s must be >= 0.0
Self-explanatory.
E: Fix ttm v_0 must be >= 0.0
Self-explanatory.
E: Fix ttm number of nodes must be > 0
Self-explanatory.
E: Cannot use fix ttm with 2d simulation
This is a current restriction of this fix due to the grid it creates.
E: Cannot use non-periodic boundares with fix ttm
This fix requires a fully periodic simulation box.
E: Cannot use fix ttm with triclinic box
This is a current restriction of this fix due to the grid it creates.
E: Electronic temperature dropped below zero
Something has gone wrong with the fix ttm electron temperature model.
E: Fix ttm electron temperatures must be > 0.0
Self-explanatory.
E: Initial temperatures not all set in fix ttm
Self-explanatory.
W: Too many inner timesteps in fix ttm
Self-explanatory.
*/