Hello everyone!
I have a problem about “fix aveforce” command. Under units metal, I try to apply pressure on my upper boundary and make it move down. In general, the compression process cannot continue after the boundary moves down a certain distance, but it keeps moving downward and completely deforms the entire structure. Below is my code.
boundary p p s
dimension 3
atom_style atomic
timestep 0.001
neighbor 2.0 bin
neigh_modify delay 1 check yes
#modeling
region box block -20.2 20.2 -20.2 20.2 0 120.8 units box
create_box 6 box
region fix block INF INF INF INF 0 10 units box
region thermo1 block INF INF INF INF 10 20 units box
region al block INF INF INF INF 20 60.4 units box
region cu block INF INF INF INF 60.4 100.8 units box
region thermo2 block INF INF INF INF 100.8 110.8 units box
region upper block INF INF INF INF 110.4 120.8 units box
lattice fcc 4.05
create_atoms 1 region fix
create_atoms 2 region thermo1
create_atoms 3 region al
lattice fcc 3.61
create_atoms 4 region cu
create_atoms 5 region thermo2
create_atoms 6 region upper
mass 1 27
mass 2 27
mass 3 27
mass 4 64
mass 5 64
mass 6 64
group fix region fix
group thermo1 region thermo1
group al region al
group cu region cu
group thermo2 region thermo2
group upper region upper
group work type 2 3 4 5
group boundary type 1 6
group up type 4 5
group low type 2 3
pair_style eam/alloy
pair_coeff * * AlCu.eam.alloy Al Al Al Cu Cu Cu
minimize 1.0e-5 1.0e-6 1000 1000
thermo 1000
thermo_style custom step temp pe ke etotal press
dump process1 all atom 1000 5mp.xyz
dump_modify process1 scale no first yes
#computation
compute thermo1temp thermo1 temp
compute cutemp cu temp
compute altemp al temp
compute thermo2temp thermo2 temp
compute s1 al stress/atom NULL
compute s11 al reduce sum c_s1[3]
variable press1 equal -(c_s11)/vol
compute s2 cu stress/atom NULL
compute s22 cu reduce sum c_s2[3]
variable press2 equal -(c_s22)/vol
compute s3 upper stress/atom NULL
compute s33 upper reduce sum c_s3[3]
variable press3 equal -(c_s33)/vol
variable dis equal xcm(move,z)
variable fzz equal fcm(move,z)
#rescale temperature at 300K
fix 1 all nve
fix boundary_f boundary setforce 0 0 0
compute uptemp up temp
velocity up create 300 5 temp uptemp units box
compute lowtemp low temp
velocity low create 300 5 temp lowtemp units box
fix reuptemp up temp/rescale 10 300 300 10.0 1
fix_modify reuptemp temp uptemp
fix relowtemp low temp/rescale 10 300 300 10.0 1
fix_modify relowtemp temp lowtemp
thermo_style custom step temp press etotal c_uptemp c_lowtemp c_cutemp c_altemp v_press1 v_press2 v_press3
run 10000
unfix reuptemp
unfix relowtemp
unfix 1
#applying pressure
variable p equal 10 #Mpa
variable s equal lx*ly*1.0e-20
variable N equal count(upper)
variable ev equal 624150913 #unit conversion newton to ev/Angstom
variable extforce equal ```1.0e6*($p*$s*${ev})/$N```
fix externalforce upper aveforce 0 0 ```$(-v_extforce)```
fix upper_nve upper nve
fix 1 work nvt temp 300 300 0.1
thermo_style custom step temp press etotal c_thermo2temp c_thermo1temp c_cutemp c_altemp v_press1 v_press2 v_press3 v_dis v_disp
run 200000
I don't know why the result is not correct, looking forward to your reply, thank you!
Best wishes
Li
Your input, as it is, contains some mistakes. Like, the move group is not defined. The least you can do is to simplify your input and perhaps use a smaller system for debugging purposes. Also, please use the syntax highlight to format the code.
Have a look at the lines that I have commented out just to visualise your system:
units metal
boundary p p s
dimension 3
atom_style atomic
timestep 0.001
neighbor 2.0 bin
neigh_modify delay 1 check yes
#modeling
region box block -20.2 20.2 -20.2 20.2 0 120.8 units box
create_box 6 box
region fix block INF INF INF INF 0 10 units box
region thermo1 block INF INF INF INF 10 20 units box
region al block INF INF INF INF 20 60.4 units box
region cu block INF INF INF INF 60.4 100.8 units box
region thermo2 block INF INF INF INF 100.8 110.8 units box
region upper block INF INF INF INF 110.4 120.8 units box
lattice fcc 4.05
create_atoms 1 region fix
create_atoms 2 region thermo1
create_atoms 3 region al
lattice fcc 3.61
create_atoms 4 region cu
create_atoms 5 region thermo2
create_atoms 6 region upper
mass 1 27
mass 2 27
mass 3 27
mass 4 64
mass 5 64
mass 6 64
group fix region fix
group thermo1 region thermo1
group al region al
group cu region cu
group thermo2 region thermo2
group upper region upper
group work type 2 3 4 5
group boundary type 1 6
group up type 4 5
group low type 2 3
pair_style eam/alloy
pair_coeff * * .local/share/lammps/potentials/AlCu.eam.alloy Al Al Al Cu Cu Cu
minimize 1.0e-5 1.0e-6 1000 1000
thermo 1000
thermo_style custom step temp pe ke etotal press
dump process1 all atom 1000 5mp.xyz
dump_modify process1 scale no first yes
#computation
compute thermo1temp thermo1 temp
compute cutemp cu temp
compute altemp al temp
compute thermo2temp thermo2 temp
compute s1 al stress/atom NULL
compute s11 al reduce sum c_s1[3]
variable press1 equal -(c_s11)/vol
compute s2 cu stress/atom NULL
compute s22 cu reduce sum c_s2[3]
variable press2 equal -(c_s22)/vol
#compute s3 move stress/atom NULL
#compute s33 move reduce sum c_s3[3]
#variable press3 equal -(c_s33)/vol
variable dis equal xcm(move,z)
variable fzz equal fcm(move,z)
#rescale temperature at 300K
fix 1 all nve
fix boundary_f boundary setforce 0 0 0
compute uptemp up temp
velocity up create 300 5 temp uptemp units box
compute lowtemp low temp
velocity low create 300 5 temp lowtemp units box
fix reuptemp up temp/rescale 10 300 300 10.0 1
fix_modify reuptemp temp uptemp
fix relowtemp low temp/rescale 10 300 300 10.0 1
fix_modify relowtemp temp lowtemp
#thermo_style custom step temp press etotal c_uptemp c_lowtemp c_cutemp c_altemp v_press1 v_press2 v_press3
thermo_style custom step temp press etotal c_uptemp c_lowtemp c_cutemp c_altemp v_press1 v_press2
run 0
Hi Dr. Roscioni,
Thank you for your reply. I don’t know how should I copy my code from VS Code to the web page keeping the syntax highlight and I try to make my code clearly. I am so sorry that it brings some troubles to you. I originally used region move instead of region upper. When I copy my code to the website, I used region upper to make it easier for people to understand my model. And it’s already a relatively small model. I reduced the pressure applied but the upper layer still keeps going down. I don’t know why and I hope you can answer my questions.
Best Regards
Li
The text editor says it supports the markdown syntax. Please see here for a quick overview: Markdown Cheat Sheet | Markdown Guide
To show code blocks you need to bracket them with triple backquotes (```).
Hi Li,
you should do the effort of posting a working input script. Even after the changes, there were missing lines and undefined variables (eg v_disp).
Anyway, a working example of your script is attached for your convenience. From what I can see, the simulation is doing exactly what you asked for. The problem you observe is due to the use of fix aveforce which keeps pushing the wall along -z, squeezing the sample unphysically. If you want the process to stop when the pressure reaches a certain threshold, you should loop over a short run and set a suitable exit criterium. See my example below.
Finally, formatting your input is just a button away:
Here is the modified input script. Run it with the command lmp -in myscript.in.
# Use a loop and an exit criterium.
variable pmax equal 100000
units metal
boundary p p s
dimension 3
atom_style atomic
timestep 0.001
neighbor 2.0 bin
neigh_modify delay 1 check yes
#modeling
region box block -20.2 20.2 -20.2 20.2 0 120.8 units box
create_box 6 box
region fix block INF INF INF INF 0 10 units box
region thermo1 block INF INF INF INF 10 20 units box
region al block INF INF INF INF 20 60.4 units box
region cu block INF INF INF INF 60.4 100.8 units box
region thermo2 block INF INF INF INF 100.8 110.8 units box
region upper block INF INF INF INF 110.4 120.8 units box
lattice fcc 4.05
create_atoms 1 region fix
create_atoms 2 region thermo1
create_atoms 3 region al
lattice fcc 3.61
create_atoms 4 region cu
create_atoms 5 region thermo2
create_atoms 6 region upper
mass 1 27
mass 2 27
mass 3 27
mass 4 64
mass 5 64
mass 6 64
group fix region fix
group thermo1 region thermo1
group al region al
group cu region cu
group thermo2 region thermo2
group upper region upper
group work type 2 3 4 5
group boundary type 1 6
group up type 4 5
group low type 2 3
pair_style eam/alloy
pair_coeff * * /home/otello/.local/share/lammps/potentials/AlCu.eam.alloy Al Al Al Cu Cu Cu
minimize 1.0e-5 1.0e-6 1000 1000
thermo 100
thermo_style custom step temp pe ke etotal press
dump process1 all atom 100 5mp_2.dump
dump_modify process1 scale no first yes
#computation
compute thermo1temp thermo1 temp
compute cutemp cu temp
compute altemp al temp
compute thermo2temp thermo2 temp
compute s1 al stress/atom NULL
compute s11 al reduce sum c_s1[3]
variable press1 equal -(c_s11)/vol
compute s2 cu stress/atom NULL
compute s22 cu reduce sum c_s2[3]
variable press2 equal -(c_s22)/vol
compute s3 upper stress/atom NULL
compute s33 upper reduce sum c_s3[3]
variable press3 equal -(c_s33)/vol
variable dis equal xcm(upper,z)
variable fzz equal fcm(upper,z)
#rescale temperature at 300K
fix 1 all nve
fix boundary_f boundary setforce 0 0 0
compute uptemp up temp
velocity up create 300 5 temp uptemp units box
compute lowtemp low temp
velocity low create 300 5 temp lowtemp units box
fix reuptemp up temp/rescale 10 300 300 10.0 1
fix_modify reuptemp temp uptemp
fix relowtemp low temp/rescale 10 300 300 10.0 1
fix_modify relowtemp temp lowtemp
thermo_style custom step etotal temp c_uptemp c_lowtemp c_cutemp c_altemp press v_press1 v_press2 v_press3 cpu
run 10000
unfix reuptemp
unfix relowtemp
unfix 1
#applying pressure
variable p equal 10 #Mpa
variable s equal lx*ly*1.0e-20
variable N equal count(upper)
variable ev equal 624150913 #unit conversion newton to ev/Angstom
variable extforce equal 1.0e6*($p*$s*${ev})/$N
fix externalforce upper aveforce 0 0 $(-v_extforce)
fix upper_nve upper nve
fix 1 work nvt temp 300 300 0.1
thermo_style custom step etotal temp c_thermo2temp c_thermo1temp c_cutemp c_altemp press v_press1 v_press2 v_press3 v_dis cpu
# Loop the simulation until the exit criterium is met.
label compress
variable ptest equal c_thermo_press
run 1000
# Repeat until the threshold is met.
if "${ptest} >= ${pmax}" then "jump SELF break"
variable ptest delete
jump SELF compress
# Exit the loop and end this simulation.
label break
print "Max pressure reached!"
Dear Dr. Roscioni,
Thank you for your reply! I will try it according to your instructions.
Best Regards
Li
Hi Dr. Roscioni,
I still have the question. The ideal situation should be that when my sample is deformed to a certain extent, the applied pressure cannot continue to squeeze the sample unless the pressure is increased. So why the wall is kept pushed after using fix aveforce. After all a structure cannot be compressed unlimitedly. Thank you!
Best Regards
Li
I see no other explanation other than the applied force is too strong.
You can exploit the loop by measuring the deformation (your job to figure that out) and start decreasing the applied force once you reach the desired conditions.
If the external force matches some experimental parameters, then the EAM potential is not adequate to describe the sample under the specified conditions.
I hope this helps.
Best,
Otello
Thank your Dr. Roscioni !
