Hii,
I am studying shock propagation in Ta. The temperature that it is calculating is very high (approx 5 times), which is not expected. Other values of pressure etc are absolutely fine. I feel the problem is that the calculation of kinetic energy is also adding centre of mass velocity and hence the temperature is high. Can I know how to solve this issue. The script is attached.
Regards,
Anuradha
Tashock.txt (5.22 KB)
Hii,
I am studying shock propagation in Ta. The temperature that it is
calculating is very high (approx 5 times), which is not expected. Other
values of pressure etc are absolutely fine. I feel the problem is that the
calculation of kinetic energy is also adding centre of mass velocity and
hence the temperature is high. Can I know how to solve this issue. The
script is attached.
So, I have done certain modifications so that the script calculates the temperature of a group of atoms, after subtracting out the center-of-mass velocity of the group using “temp/com” . the compute statement works but it shows an error in the statement where I am trying to save this temperature in a file called temp.${vpis} . the error is " Fix ave/chunk compute does not calculate per-atom values" .The modified script is :
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
variable a loop 9
variable Ts equal 300
variable alat equal 3.304
variable myseed equal 12345
variable xm equal 10
variable ym equal 10
variable zm equal 420 #LJ
variable dt equal 0.001
variable t_eq equal 10000
variable t_sh equal 50000
variable vpis index 0.1 0.3 0.5 1.0 1.3 1.8 2.2 2.5 3
variable Nevery equal 10
variable Nrepeat equal 5
variable Nfreq equal 500
variable dz equal 10
variable atomrate equal 5000
variable tdamp equal “v_dt100"
variable pdamp equal "v_dt1000”
variable Up equal “10*v_vpis”
timestep ${dt}
lattice bcc ${alat} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
#tmd 27g/NA 6e23 *4 atoms/cell/(4.05e-8cm)3 = 2.7 g/cm3
region sim_box block 0 {xm} 0 {ym} 0 ${zm} units lattice
create_box 2 sim_box
region atom_box block 0 {xm} 0 {ym} 0 ${zm}
create_atoms 1 region atom_box
group atom_box region atom_box
region piston block INF INF INF INF INF 4
region bulk block INF INF INF INF 4 INF
group piston region piston
group bulk region bulk
set group piston type 1
set group bulk type 2
mass 1 180.95
mass 2 180.95
#pair_style lj/cut 6.38
#pair_coeff * * 3.34 2.55
pair_style eam/alloy
pair_coeff * * Ta.eam.alloy Ta Ta
#mass 1 180.95
#mass 2 180.95
compute myCN bulk cna/atom 3.0
compute myKE bulk ke/atom
compute myPE bulk pe/atom
compute myCOM bulk com
compute peratom bulk stress/atom NULL
compute vz bulk property/atom vz
compute vorvol bulk voronoi/atom
compute tempcomm bulk temp/com
reset_timestep 0
velocity atom_box create {Ts} {myseed} rot yes dist gaussian
fix equilibration bulk npt temp {Ts} {Ts} {tdamp} iso 0 0 {pdamp} drag 1
variable eq1 equal “step”
variable eq2 equal “pxx/10000”
variable eq3 equal “pyy/10000”
variable eq4 equal “pzz/10000”
variable eq5 equal “lx”
variable eq6 equal “ly”
variable eq7 equal “lz”
variable eq8 equal “temp”
variable eq9 equal “etotal”
variable eq11 equal “density”
fix output1 bulk print 10 “{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} &
{eq7} {eq8} {eq9} {eq11}” file run.out.${vpis} screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal density
run ${t_eq}
unfix equilibration
unfix output1
change_box all boundary p p s
reset_timestep 0
fix 1 all nve
velocity piston set 0 0 v_Up sum no units box
fix 2 piston setforce 0.0 0.0 0.0
compute 3 bulk chunk/atom bin/1d z lower ${dz} units box
fix density_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 density/mass file denz.{vpis}
fix temp_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 c_tempcomm file temp.{vpis}
variable meanpress atom -(c_peratom[1]+c_peratom[2]+c_peratom[3])/3/(c_vorvol[1]+1e-99)
fix pressure_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 v_meanpress file pressure.{vpis}
fix velZ_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 c_vz file velocityZcomp.{vpis}
variable eq1 equal “step”
variable eq2 equal “pxx/10000”
variable eq3 equal “pyy/10000”
variable eq4 equal “pzz/10000”
variable eq5 equal “lx”
variable eq6 equal “ly”
variable eq7 equal “lz”
variable eq8 equal “temp”
variable eq9 equal “etotal”
variable eq10 equal “c_myCOM[3]”
variable eq11 equal “density”
fix shock bulk print 10 "{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} {eq7} {eq8} {eq9} {eq10} {eq11}" file run.{Ts}K.${vpis} screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal c_myCOM[3] density
#Use cfg(?) for AtomEye
dump 1 all cfg {atomrate} dump._*.{vpis}.cfg mass type xs ys zs id c_myPE c_myKE c_myCN
dump_modify 1 element Ta Ta
#dump 2 all image 1000 image.*.jpg type type axes yes 0.8 0.02 view 120 -90
run ${t_sh}
clear
next vpis
next a
jump Tashock.in
temp/com is a global scalar compute; fix ave/chunk is designed for per-atom computes. You need to use fix ave/time if you want lammps to time average global scalars.
Thanx for the help. As shown in the script below , rather than use of temp/com I used temp/chunk . . But now the problem is how to print the temperature in different chunk. I am able to print time averaged temperature using fix temperature_profile statement. But not the time averaged temperature profile in different bins.
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
variable a loop 9
variable Ts equal 300
variable alat equal 3.304
variable myseed equal 12345
variable xm equal 10
variable ym equal 10
variable zm equal 420 #LJ
variable dt equal 0.001
variable t_eq equal 10000
variable t_sh equal 50000
variable vpis index 0.1 0.3 0.5 1.0 1.3 1.8 2.2 2.5 3
variable Nevery equal 10
variable Nrepeat equal 5
variable Nfreq equal 500
variable dz equal 10
variable atomrate equal 5000
variable tdamp equal “v_dt100"
variable pdamp equal "v_dt1000”
variable Up equal “10*v_vpis”
timestep ${dt}
lattice bcc ${alat} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
#tmd 27g/NA 6e23 *4 atoms/cell/(4.05e-8cm)3 = 2.7 g/cm3
region sim_box block 0 {xm} 0 {ym} 0 ${zm} units lattice
create_box 2 sim_box
region atom_box block 0 {xm} 0 {ym} 0 ${zm}
create_atoms 1 region atom_box
group atom_box region atom_box
region piston block INF INF INF INF INF 4
region bulk block INF INF INF INF 4 INF
group piston region piston
group bulk region bulk
set group piston type 1
set group bulk type 2
mass 1 180.95
mass 2 180.95
#pair_style lj/cut 6.38
#pair_coeff * * 3.34 2.55
pair_style eam/alloy
pair_coeff * * Ta.eam.alloy Ta Ta
#mass 1 180.95
#mass 2 180.95
compute myCN bulk cna/atom 3.0
compute myKE bulk ke/atom
compute myPE bulk pe/atom
compute myCOM bulk com
compute peratom bulk stress/atom NULL
compute vz bulk property/atom vz
compute vorvol bulk voronoi/atom
reset_timestep 0
velocity atom_box create {Ts} {myseed} rot yes dist gaussian
fix equilibration bulk npt temp {Ts} {Ts} {tdamp} iso 0 0 {pdamp} drag 1
variable eq1 equal “step”
variable eq2 equal “pxx/10000”
variable eq3 equal “pyy/10000”
variable eq4 equal “pzz/10000”
variable eq5 equal “lx”
variable eq6 equal “ly”
variable eq7 equal “lz”
variable eq8 equal “temp”
variable eq9 equal “etotal”
variable eq11 equal “density”
fix output1 bulk print 10 “{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} &
{eq7} {eq8} {eq9} {eq11}” file run.out.${vpis} screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal density
run ${t_eq}
unfix equilibration
unfix output1
change_box all boundary p p s
reset_timestep 0
fix 1 all nve
velocity piston set 0 0 v_Up sum no units box
fix 2 piston setforce 0.0 0.0 0.0
compute 3 bulk chunk/atom bin/1d z lower ${dz} units box
fix density_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 density/mass file denz.{vpis}
compute tempcomm bulk temp/chunk 3 temp com yes
fix temp_profile bulk ave/time {Nevery} {Nrepeat} {Nfreq} c_tempcomm file temp.{vpis}
variable meanpress atom -(c_peratom[1]+c_peratom[2]+c_peratom[3])/3/(c_vorvol[1]+1e-99)
fix pressure_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 v_meanpress file pressure.{vpis}
fix velZ_profile bulk ave/chunk {Nevery} {Nrepeat} {Nfreq} 3 c_vz file velocityZcomp.{vpis}
variable eq1 equal “step”
variable eq2 equal “pxx/10000”
variable eq3 equal “pyy/10000”
variable eq4 equal “pzz/10000”
variable eq5 equal “lx”
variable eq6 equal “ly”
variable eq7 equal “lz”
variable eq8 equal “temp”
variable eq9 equal “etotal”
variable eq10 equal “c_myCOM[3]”
variable eq11 equal “density”
fix shock bulk print 10 "{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} {eq7} {eq8} {eq9} {eq10} {eq11}" file run.{Ts}K.${vpis} screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal c_myCOM[3] density
#Use cfg(?) for AtomEye
dump 1 all cfg {atomrate} dump._*.{vpis}.cfg mass type xs ys zs id c_myPE c_myKE c_myCN
dump_modify 1 element Ta Ta
#dump 2 all image 1000 image.*.jpg type type axes yes 0.8 0.02 view 120 -90
run ${t_sh}
clear
next vpis
next a
jump Tashock.in
Your script is 336 lines long. We can not possibly debug it for you. Instead of looking for a shortcuts to enlightment on the mailing list, try reading the documentation and doing your own testing. I suggest starting with http://lammps.sandia.gov/doc/fix_ave_chunk.html
