The first problem is that you didn’t read the forum guidelines, but you can still do it and edit your first post accordingly.
The second problem is that you need to communicate your problem effectively, describing rigorously what problem(s) you encountered, providing context, etc. It is way too vague to state that “the simulation doesn’t work”, as no one is interested in doing the homework on your behalf.
Regarding your script, I have made the following changes:
- Simplified the definition of fixed and mobile regions,
- Used the atom types to distinguish between fixed and mobile atoms,
- Changed the position of the heat source (please check
z2), - Added a random component to
fix addforceto mimic the effect of a heat source and changed the Gaussian definition to my own taste, - Removed the
fix setforce, as it suffices to integrate the mobile atoms. - Thermalised the mobile atoms before applying the heat source,
- Used an NVE integrator for the atoms on which the force is applied,
Here is the revised script, with proper quotation:
# Initialization
units metal
dimension 3
boundary p p p
atom_style atomic
processors 3 2 2
# Define the simulation area
region whole block 0 300 0 300 0 500 units box
create_box 2 whole
# Create an electrode
lattice fcc 3.615
region electrode block 100 200 100 150 230 330 units box
create_atoms 1 region electrode
group electrode region electrode
# Create artifacts
lattice fcc 3.615
region workpiece block 100 200 100 150 100 200 units box
create_atoms 1 region workpiece
group workpiece region workpiece
# Define atomic type and mass
mass 1 63.546 # Cu
mass 2 63.546 # Cu
# Select the applicable force field
pair_style eam/alloy
pair_coeff * * Cu_zhou.eam.alloy Cu Cu
# Define the side boundary layer area of the electrode
region el_thermo block 105 195 105 145 230 325 units box
region el_side intersect 2 electrode el_thermo
region wp_thermo block 105 195 105 145 105 200 units box
region wp_side intersect 2 workpiece wp_thermo
region mobile union 2 el_thermo wp_thermo
region side union 2 el_side wp_side
group mobile region mobile
group side region side
set group side type 2
# The thermostating region of the electrode
region el_hwc block 105 195 105 145 250 325 units box
region wp_hwc block 105 195 105 145 105 180 units box
region hwc union 2 el_hwc wp_hwc # nvt
group hwc region hwc
# Newtonian layer
region el_ndc intersect 2 el_thermo el_hwc # addforce
region wp_ndc intersect 2 wp_thermo wp_hwc # addforce
region ndc union 2 el_ndc wp_ndc # nve
group wp_ndc region wp_ndc
group el_ndc region el_ndc
group ndc region ndc
# Thermodynamic output.
comm_style tiled
balance 1.0 rcb
compute mobile_temp mobile temp
thermo 3000
thermo_style custom step temp pe ke etotal press pxx pyy pzz
thermo_modify temp mobile_temp
dump 1 all custom 500 dump.Cu.xyz id type x y z
# First, relax the Cu blocks.
velocity all create 300 68748
fix NVT mobile nvt temp 300 300 0.1
run 10000
unfix NVT
# A variable that controls whether the heat source is activated or not, stop heating after 15000 steps
variable heat_active equal "step < 15000"
variable heat_mult equal "v_heat_active * 20.8"
variable heat_mult_2 equal "v_heat_active * 8.8"
# The first heat source
variable x0 equal 120
variable y0 equal 115
variable z0 equal 240 # Central position of heat source
variable Pm equal ${heat_mult} # Maximum heating intensity
variable k equal 0.03 # Parameters that control the range of Gaussian distribution
variable fx atom v_Pm*exp(-v_k*(x-v_x0)^2)*random(-1,1,37324)
variable fy atom v_Pm*exp(-v_k*(y-v_y0)^2)*random(-1,1,82396)
variable fz atom v_Pm*exp(-v_k*(z-v_z0)^2)*random(-1,1,72383)
# Applying the force of Gaussian distribution to the Newtonian layer
fix heat el_ndc addforce v_fx v_fy v_fz every 10
# The second heat source
variable x2 equal 155
variable y2 equal 130
variable z2 equal 190
variable Pm2 equal ${heat_mult_2}
variable fx_2 atom v_Pm2*exp(-v_k*(x-v_x2)^2)*random(-1,1,24535)
variable fy_2 atom v_Pm2*exp(-v_k*(y-v_y2)^2)*random(-1,1,61699)
variable fz_2 atom v_Pm2*exp(-v_k*(z-v_z2)^2)*random(-1,1,12293)
fix heat2 wp_ndc addforce v_fx_2 v_fy_2 v_fz_2 every 10
# Integration and temperature control on the electrodes' bulk.
fix NVT hwc nvt temp 300 300 0.1
# Integration only on the electrodes' surface.
fix NVE ndc nve
run 60000
WARNING: One or more atoms are time-integrated more than once. I am not sure why this happens, as the regions hwc and ndc should be defined as not to have any overlap. I leave the debugging to you (or any other curious user).