diff --git a/src/fix_ttm.cpp b/src/fix_ttm.cpp index 5dbc8915a8..3fdbd7109d 100644 --- a/src/fix_ttm.cpp +++ b/src/fix_ttm.cpp @@ -12,7 +12,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Paul Crozier (SNL) + Contributing authors: Paul Crozier (SNL) + Carolyn Phillips (University of Michigan) ------------------------------------------------------------------------- */ #include "mpi.h" @@ -45,6 +46,12 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) : { if (narg < 15) error->all("Illegal fix ttm command"); + vector_flag = 1; + size_vector = 2; + scalar_vector_freq = 1; + extvector = 1; + nevery = 1; + int seed = atoi(arg[3]); electronic_specific_heat = atof(arg[4]); electronic_density = atof(arg[5]); @@ -109,45 +116,35 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) : total_nnodes = nxnodes*nynodes*nznodes; nsum = memory->create_3d_int_array(nxnodes,nynodes,nznodes,"ttm:nsum"); - nsum_prime = memory->create_3d_int_array(nxnodes,nynodes,nznodes, - "ttm:nsum_prime"); nsum_all = memory->create_3d_int_array(nxnodes,nynodes,nznodes, - "ttm:nsum_all"); - nsum_prime_all = memory->create_3d_int_array(nxnodes,nynodes,nznodes, - "ttm:nsum_prime_all"); + "ttm:nsum_all"); T_initial_set = memory->create_3d_int_array(nxnodes,nynodes,nznodes, - "ttm:T_initial_set"); + "ttm:T_initial_set"); sum_vsq = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_vsq"); - sum_vsq_prime = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_vsq_prime"); + "ttm:sum_vsq"); sum_mass_vsq = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_mass_vsq"); - sum_mass_vsq_prime = - memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_mass_vsq_prime"); + "ttm:sum_mass_vsq"); sum_vsq_all = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_vsq_all"); - sum_vsq_prime_all = - memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_vsq_prime_all"); + "ttm:sum_vsq_all"); sum_mass_vsq_all = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_mass_vsq_all"); - sum_mass_vsq_prime_all = - memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:sum_mass_vsq_prime_all"); + "ttm:sum_mass_vsq_all"); T_electron_old = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:T_electron_old"); + "ttm:T_electron_old"); T_electron = memory->create_3d_double_array(nxnodes,nynodes,nznodes,"ttm:T_electron"); - T_a = memory->create_3d_double_array(nxnodes,nynodes,nznodes,"ttm:T_a"); - T_a_prime = memory->create_3d_double_array(nxnodes,nynodes,nznodes, - "ttm:T_a_prime"); - g_s = memory->create_3d_double_array(nxnodes,nynodes,nznodes,"ttm:g_s"); - g_p = memory->create_3d_double_array(nxnodes,nynodes,nznodes,"ttm:g_p"); + net_energy_transfer = + memory->create_3d_double_array(nxnodes,nynodes,nznodes, + "TTM:net_energy_transfer"); + net_energy_transfer_all = + memory->create_3d_double_array(nxnodes,nynodes,nznodes, + "TTM:net_energy_transfer_all"); + + flangevin = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); // set initial electron temperatures from user input file @@ -167,24 +164,17 @@ FixTTM::~FixTTM() delete [] gfactor2; memory->destroy_3d_int_array(nsum); - memory->destroy_3d_int_array(nsum_prime); memory->destroy_3d_int_array(nsum_all); - memory->destroy_3d_int_array(nsum_prime_all); memory->destroy_3d_int_array(T_initial_set); memory->destroy_3d_double_array(sum_vsq); - memory->destroy_3d_double_array(sum_vsq_prime); memory->destroy_3d_double_array(sum_mass_vsq); - memory->destroy_3d_double_array(sum_mass_vsq_prime); memory->destroy_3d_double_array(sum_vsq_all); - memory->destroy_3d_double_array(sum_vsq_prime_all); memory->destroy_3d_double_array(sum_mass_vsq_all); - memory->destroy_3d_double_array(sum_mass_vsq_prime_all); memory->destroy_3d_double_array(T_electron_old); memory->destroy_3d_double_array(T_electron); - memory->destroy_3d_double_array(T_a); - memory->destroy_3d_double_array(T_a_prime); - memory->destroy_3d_double_array(g_s); - memory->destroy_3d_double_array(g_p); + memory->destroy_2d_double_array(flangevin); + memory->destroy_3d_double_array(net_energy_transfer); + memory->destroy_3d_double_array(net_energy_transfer_all); } /* ---------------------------------------------------------------------- */ @@ -194,6 +184,7 @@ int FixTTM::setmask() int mask = 0; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; + mask |= END_OF_STEP; return mask; } @@ -216,6 +207,11 @@ void FixTTM::init() 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 (strcmp(update->integrate_style,"respa") == 0) nlevels_respa = ((Respa *) update->integrate)->nlevels; } @@ -237,8 +233,6 @@ void FixTTM::setup(int vflag) void FixTTM::post_force(int vflag) { - update_electron_temperatures(); - double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -266,6 +260,9 @@ void FixTTM::post_force(int vflag) while (iynode < 0) iynode += nynodes; while (iznode < 0) iznode += nznodes; + if (T_electron[ixnode][iynode][iznode] < 0) + error->all("Electronic temperature dropped to below zero"); + double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]); gamma1 = gfactor1[type[i]]; @@ -273,9 +270,13 @@ void FixTTM::post_force(int vflag) if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p; gamma2 = gfactor2[type[i]] * tsqrt; - f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5); - f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5); - f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5); + 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]; } } } @@ -335,7 +336,7 @@ void FixTTM::read_initial_electron_temperatures() /* ---------------------------------------------------------------------- */ -void FixTTM::update_electron_temperatures() +void FixTTM::end_of_step() { double **x = atom->x; double **v = atom->v; @@ -344,32 +345,14 @@ void FixTTM::update_electron_temperatures() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - - double massone; - - // compute atomic Ta, and Ta' (for high v atoms) 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_prime[ixnode][iynode][iznode] = 0; - nsum_all[ixnode][iynode][iznode] = 0; - nsum_prime_all[ixnode][iynode][iznode] = 0; - sum_vsq[ixnode][iynode][iznode] = 0.0; - sum_vsq_prime[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq_prime[ixnode][iynode][iznode] = 0.0; - sum_vsq_all[ixnode][iynode][iznode] = 0.0; - sum_vsq_prime_all[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq_prime_all[ixnode][iynode][iznode] = 0.0; - } - + 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) { - 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; @@ -382,65 +365,20 @@ void FixTTM::update_electron_temperatures() 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; - if (vsq > v_0_sq) { - nsum_prime[ixnode][iynode][iznode] += 1; - sum_vsq_prime[ixnode][iynode][iznode] += vsq; - sum_mass_vsq_prime[ixnode][iynode][iznode] += massone*vsq; - } + 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; - MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes, - MPI_INT,MPI_SUM,world); - MPI_Allreduce(&nsum_prime[0][0][0],&nsum_prime_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_vsq_prime[0][0][0],&sum_vsq_prime_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); - MPI_Allreduce(&sum_mass_vsq_prime[0][0][0],&sum_mass_vsq_prime_all[0][0][0], - total_nnodes,MPI_DOUBLE,MPI_SUM,world); - - double max_g_p = 0.0; - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - if (nsum_all[ixnode][iynode][iznode] > 0) { - T_a[ixnode][iynode][iznode] = - sum_mass_vsq_all[ixnode][iynode][iznode]/ - (3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e); - double g_p_tmp = force->mvv2e*gamma_p*sum_vsq_all[ixnode][iynode][iznode]/ - T_a[ixnode][iynode][iznode]/del_vol; - max_g_p = MAX(max_g_p,g_p_tmp); - g_p[ixnode][iynode][iznode] = g_p_tmp; - } else { - T_a[ixnode][iynode][iznode] = 0; - g_p[ixnode][iynode][iznode] = 0; - } - if (nsum_prime_all[ixnode][iynode][iznode] > 0) { - T_a_prime[ixnode][iynode][iznode] = - sum_mass_vsq_prime_all[ixnode][iynode][iznode]/ - (3.0*force->boltz*nsum_prime_all[ixnode][iynode][iznode] / - force->mvv2e); - g_s[ixnode][iynode][iznode] = - force->mvv2e*gamma_s*sum_vsq_prime_all[ixnode][iynode][iznode]/ - T_a_prime[ixnode][iynode][iznode]/del_vol; - } else { - T_a_prime[ixnode][iynode][iznode] = 0; - g_s[ixnode][iynode][iznode] = 0; - } - } - // num_inner_timesteps = # of inner steps (thermal solves) // required this MD step to maintain a stable explicit solve @@ -448,12 +386,10 @@ void FixTTM::update_electron_temperatures() 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) + - 0.5*max_g_p); + (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) + - 0.5*max_g_p); + (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) @@ -467,7 +403,7 @@ void FixTTM::update_electron_temperatures() 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]; + T_electron[ixnode][iynode][iznode]; // compute new electron T profile @@ -487,33 +423,81 @@ void FixTTM::update_electron_temperatures() 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) + - g_p[ixnode][iynode][iznode] * - (T_a[ixnode][iynode][iznode] - - T_electron_old[ixnode][iynode][iznode]) + - g_s[ixnode][iynode][iznode]*T_a_prime[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 (!(update->ntimestep % nfileevery) && (me == 0)) { + if ((nfileevery) && !(update->ntimestep % nfileevery) && (me == 0)) { fprintf(fp,"%d ",update->ntimestep); + + // 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++) - fprintf(fp,"%f ",T_a[ixnode][iynode][iznode]); + 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); + + 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); + fprintf(fp,"%f ",T_a); + } + fprintf(fp,"\t"); for (int ixnode = 0; ixnode < nxnodes; ixnode++) for (int iynode = 0; iynode < nynodes; iynode++) @@ -534,3 +518,42 @@ double FixTTM::memory_usage() bytes += 14*total_nnodes * sizeof(double); return bytes; } + +/* ---------------------------------------------------------------------- */ + +void FixTTM::grow_arrays(int ngrow) +{ + + flangevin = memory->grow_2d_double_array(flangevin,ngrow,3,"TTM:flangevin"); + +} + +/* ---------------------------------------------------------------------- + return the energy of the electronic subsystem or the net_energy transfer + between the subsystems +------------------------------------------------------------------------- */ + +double FixTTM::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; +} diff --git a/src/fix_ttm.h b/src/fix_ttm.h index a223c3138c..64b6b7d3fa 100644 --- a/src/fix_ttm.h +++ b/src/fix_ttm.h @@ -27,8 +27,11 @@ class FixTTM : public Fix { void setup(int); void post_force(int); void post_force_respa(int, int, int); + void end_of_step(); void reset_dt(); double memory_usage(); + void grow_arrays(int); + double compute_vector(int); private: int me; @@ -37,21 +40,19 @@ class FixTTM : public Fix { class RanMars *random; FILE *fp,*fpr; int nxnodes,nynodes,nznodes,total_nnodes; - int ***nsum,***nsum_prime; - int ***nsum_all,***nsum_prime_all,***T_initial_set; + int ***nsum; + int ***nsum_all,***T_initial_set; double *gfactor1,*gfactor2,*ratio; - double ***T_electron,***T_a,***T_a_prime,***g_p,***g_s; - double ***sum_vsq,***sum_vsq_prime; - double ***sum_mass_vsq,***sum_mass_vsq_prime; - double ***sum_vsq_all,***sum_vsq_prime_all; - double ***sum_mass_vsq_all,***sum_mass_vsq_prime_all; - double ***T_electron_old; + double **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(); - void update_electron_temperatures(); }; }