Hello everyone, I am a student who has just started learning LAMMPS. I am currently studying the fix ttm/mod command and attempting to reproduce some existing simulations of one-dimensional metal ablation. I noticed the following description in the LAMMPS manual:
The current fix ttm/mod implementation allows TTM simulations with a vacuum. The vacuum region is defined as the grid cells with zero electronic temperature. The numerical scheme does not allow energy exchange with such cells. Since the material can expand to previously unoccupied region in some simulations, the vacuum border can be allowed to move. It is controlled by the surface_movement parameter in the init_file. If it is set to 1, then “vacuum” cells can be changed to “electron-filled” cells with the temperature T_e_min if atoms move into them (currently only implemented for the case of 1-dimensional motion of a flat surface normal to the X axis). The initial locations of the interfaces of the electron density to the vacuum can be set in the init_file via lsurface and rsurface parameters.
However, during my actual use of fix ttm/mod, I found that in the first timestep of the simulation, the output file for the electron temperature grid shows that all grid cells that were initially set to vacuum (i.e., 0 K) were assigned T_e_min. I have tried many times and the result is always the same.
So I checked the fix_ttm_mod.cpp file on GitHub and found the following code between lines 857–870:
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);
if (movsur == 1) {
if (T_electron[ixnode][iynode][iznode] == 0.0)
T_electron[ixnode][iynode][iznode] = electron_temperature_min;
}
}
fmt::print(fp, " {}", T_a);
}
I believe this code is what caused the issue I encountered.
Therefore, I would like to ask the community: if I want to implement the functionality described in the LAMMPS manual, should I modify the code as follows?
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);
if (movsur == 1) {
if (T_electron[t_surface_l][iynode][iznode] == 0.0) {
T_electron[t_surface_l][iynode][iznode] = electron_temperature_min;
}
}
}
fmt::print(fp, " {}", T_a);
}
This is just a preliminary thought from me as a beginner, and I would greatly appreciate any feedback or guidance from experienced researchers.