This fix computes 2 output quantities stored in a vector of +length 2, which can be accessed by various output +commands. The first quantity is +the total energy of the electronic subsystem. The second quantity +is the energy transferred from the electronic to the atomic subsystem +on that timestep. Note that the velocity verlet integrator applies the +fix ttm forces to the atomic subsystem as two half-step velocity +updates: one on the current timestep and one on the subsequent timestep. +Consequently, the change in the atomic subsystem energy is lagged by +half a timestep relative to the change in the electronic subsystem +energy. As a result of this, users may notice slight fluctuations in +the sum of the atomic and electronic subsystem energies reported at +the end of the timestep. +
+The vector values calculated by this fix are "extensive", +meaning they scale with the number of atoms in the simulation. +
IMPORTANT NOTE: The current implementation creates a copy of the electron grid that overlays the entire simulation domain, for each processor. Values on the grid are summed across all processors. Thus @@ -173,6 +190,8 @@
Default: none
+(Duffy) D M Duffy and A M Rutherford, J. Phys.: Condens. Matter, 19,
diff -Naur lammps-1Jun09/doc/fix_ttm.txt lammps-2Jun09/doc/fix_ttm.txt
--- lammps-1Jun09/doc/fix_ttm.txt 2009-04-06 15:03:10.000000000 -0600
+++ lammps-2Jun09/doc/fix_ttm.txt 2009-06-03 14:33:21.000000000 -0600
@@ -142,6 +142,23 @@
temperature controlled by another fix - e.g. "fix nvt"_fix_nvt.html or
"fix langevin"_fix_langevin.html.
+This fix computes 2 output quantities stored in a vector of
+length 2, which can be accessed by various "output
+commands"_Section_howto.html#4_15. The first quantity is
+the total energy of the electronic subsystem. The second quantity
+is the energy transferred from the electronic to the atomic subsystem
+on that timestep. Note that the velocity verlet integrator applies the
+fix ttm forces to the atomic subsystem as two half-step velocity
+updates: one on the current timestep and one on the subsequent timestep.
+Consequently, the change in the atomic subsystem energy is lagged by
+half a timestep relative to the change in the electronic subsystem
+energy. As a result of this, users may notice slight fluctuations in
+the sum of the atomic and electronic subsystem energies reported at
+the end of the timestep.
+
+The vector values calculated by this fix are "extensive",
+meaning they scale with the number of atoms in the simulation.
+
IMPORTANT NOTE: The current implementation creates a copy of the
electron grid that overlays the entire simulation domain, for each
processor. Values on the grid are summed across all processors. Thus
@@ -170,6 +187,8 @@
[Default:] none
+:line
+
:link(Duffy)
[(Duffy)] D M Duffy and A M Rutherford, J. Phys.: Condens. Matter, 19,
016207-016218 (2007).
diff -Naur lammps-1Jun09/src/fix_ttm.cpp lammps-2Jun09/src/fix_ttm.cpp
--- lammps-1Jun09/src/fix_ttm.cpp 2009-05-19 09:04:19.000000000 -0600
+++ lammps-2Jun09/src/fix_ttm.cpp 2009-06-03 14:33:16.000000000 -0600
@@ -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 @@
{
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 @@
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 @@
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 mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
+ mask |= END_OF_STEP;
return mask;
}
@@ -216,6 +207,11 @@
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::post_force(int vflag)
{
- update_electron_temperatures();
-
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@@ -266,6 +260,9 @@
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 @@
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::update_electron_temperatures()
+void FixTTM::end_of_step()
{
double **x = atom->x;
double **v = atom->v;
@@ -344,32 +345,14 @@
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 @@
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 @@
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