diff --git a/src/EXTRA-FIX/fix_ttm_old.cpp b/src/EXTRA-FIX/fix_ttm_old.cpp new file mode 100644 index 0000000000..80db9e3fc2 --- /dev/null +++ b/src/EXTRA-FIX/fix_ttm_old.cpp @@ -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 +#include +#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(xscale*nxnodes); + int iynode = static_cast(yscale*nynodes); + int iznode = static_cast(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(xscale*nxnodes); + int iynode = static_cast(yscale*nynodes); + int iznode = static_cast(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(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(xscale*nxnodes); + int iynode = static_cast(yscale*nynodes); + int iznode = static_cast(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 (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 (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; +} diff --git a/src/EXTRA-FIX/fix_ttm_old.h b/src/EXTRA-FIX/fix_ttm_old.h new file mode 100644 index 0000000000..078e7d4cb7 --- /dev/null +++ b/src/EXTRA-FIX/fix_ttm_old.h @@ -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. + +*/