Please check this discussion on region intersect. I am attaching the revised script, which runs without errors. Your homework is to define a fix to compute the amount of energy added to the two Cu blocks over time: with your initial choice (and the NVE integrator for the heated region), the surfaces kind of evaporate at supersonic speed. Instead, you may want to control carefully the heat flux added to the system.
The new sample (cut into half to visualise the internal structure):
# 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 two electrodes
lattice fcc 3.615
region electrode block 100 200 100 150 230 330 units box
region workpiece block 100 200 100 150 100 200 units box
create_atoms 1 region electrode
create_atoms 1 region workpiece
group electrode region electrode
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
# Partition the electrodes into fixed and mobile regions.
region el_bulk block 105 195 104 145 230 327 units box
region el_surf block 100 200 100 150 230 235 units box
region el_thermo union 2 el_bulk el_surf
region wp_bulk block 105 195 104 145 104 200 units box
region wp_surf block 100 200 100 150 195 200 units box
region wp_thermo union 2 wp_bulk wp_surf
group el_thermo region el_thermo
group wp_thermo region wp_thermo
group el_side subtract electrode el_thermo
group wp_side subtract workpiece wp_thermo
group mobile union el_thermo wp_thermo
group side union el_side wp_side
set group side type 2
# The thermostating region of the electrode
region el_hwc block 105 195 104 145 250 327 units box
region wp_hwc block 105 195 104 145 104 180 units box
group el_hwc region el_hwc
group wp_hwc region wp_hwc
group hwc union el_hwc wp_hwc # nvt
# Newtonian layer
group el_ndc subtract el_thermo el_hwc # addforce
group wp_ndc subtract wp_thermo wp_hwc # addforce
group ndc union el_ndc wp_ndc # nve
# Thermodynamic output.
comm_style tiled
balance 1.0 rcb
compute mobile_temp mobile temp
thermo 500
thermo_style custom step temp pe ke etotal press pxx pyy pzz
thermo_modify temp mobile_temp flush yes
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 3000
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 * 10."
variable heat_mult2 equal "v_heat_active * 5."
# The first heat source
variable x1 equal 120
variable y1 equal 115
variable z1 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_x1)^2)*random(-1,1,37324)
variable fy atom v_Pm*exp(-v_k*(y-v_y1)^2)*random(-1,1,82396)
variable fz atom v_Pm*exp(-v_k*(z-v_z1)^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 100
# The second heat source
variable x2 equal 155
variable y2 equal 130
variable z2 equal 190
variable Pm2 equal ${heat_mult2}
variable fx2 atom v_Pm2*exp(-v_k*(x-v_x2)^2)*random(-1,1,24535)
variable fy2 atom v_Pm2*exp(-v_k*(y-v_y2)^2)*random(-1,1,61699)
variable fz2 atom v_Pm2*exp(-v_k*(z-v_z2)^2)*random(-1,1,12293)
fix heat2 wp_ndc addforce v_fx2 v_fy2 v_fz2 every 100
# 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