Dear All,
Since I saw this warning on the heat/flux page
“The compute heat/flux has been reported to produce unphysical values for angle, dihedral, improper and constraint force contributions when used with compute stress/atom, as discussed in (Surblys2019), (Boone) and (Surblys2021). You are strongly advised to use compute centroid/stress/atom, which has been implemented specifically for such cases.”
I tried using the compute centroid/stress/atom command. However, I’m using a hybrid/overlay pair_style with lj/cut and eam.
I get the following error:
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:231)
ERROR: Pair style does not support compute centroid/stress/atom (src/compute_centroid_stress_atom.cpp:130)
However, in the compute centroid/stress/atom command page under restrictions it says it does not support pair styles with many-body interactions but EAM is an exception.
Since EAM is an exception why do I get this error?
Do you have any idea for a work around this problem? Thank you!
Below is my input script. I have lammps 29Sep2021 installed.
Please find the log file attached.
------------------------ INITIALIZATION ----------------------------
units real
dimension 3
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
variable latparam equal 5.0553
variable T equal 298
variable V equal vol
variable dt equal 0.04
variable p equal 20000 # correlation length
variable s equal 1 # sample interval
variable d equal $p*s # dump interval
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {kCal2J}*{kCal2J}/{fs2s}/${A2m}
----------------------- ATOM DEFINITION ----------------------------
lattice fcc ${latparam}
region box block -50 50 -50 50 -50 50 units box
create_box 5 box bond/types 9 angle/types 14 extra/bond/per/atom 9 extra/angle/per/atom 14 extra/special/per/atom 10
region liquid_1 sphere 0 0 0 35 side out units box
region liquid_0 block -50 50 -50 50 -50 50 units box
region liquid intersect 2 liquid_1 liquid_0 units box
molecule ethylene_glycol lammps.data
create_atoms 0 region liquid mol ethylene_glycol 52874 ratio 0.041346233 74637
region particle_1 sphere 0 0 0 35 units box
#create_atoms 5 region particle_1 #ratio 0.363636364 74637
create_atoms 5 random 540 489571 particle_1
#read_data argoncopper1per.data
group solid type 5
group liquid type 1 2 3 4
mass 1 12.011 #C
mass 2 15.999 #O
mass 3 1.008 #H1
mass 4 1.008 #H2
mass 5 63.546 #Cu
------------------------ FORCE FIELDS ------------------------------
pair_style hybrid/overlay lj/cut 8 eam
#pair_style eam/fs
pair_coeff 5 5 eam Cu_u3.eam
#pair_coef type_i type_j epsilon sigma rcut
#pair_coeff * * lj/cut 0.030 2.5 8 #all
pair_coeff 1 1 lj/cut 0.066 3.5 8 #C-C
pair_coeff 2 2 lj/cut 0.170 3.12 8 #O-O
pair_coeff 3 3 lj/cut 0.030 2.5 8 #H1-H1
pair_coeff 4 4 lj/cut 0 0 8 #H2-H2
pair_coeff 5 5 lj/cut 4.72 2.616 8 #Cu-Cu
pair_coeff 1 2 lj/cut 0.01122 2.5 8 #C-O
pair_coeff 1 3 lj/cut 0.030 2.5 8 #C-H1
pair_coeff 1 4 lj/cut 0 0 8 #C-H2
pair_coeff 1 5 lj/cut 0.31152 2.5 8 #C-Cu
pair_coeff 2 3 lj/cut 0.170 3.12 8 #O-H1
pair_coeff 2 4 lj/cut 0 0 8 #O-H2
pair_coeff 2 5 lj/cut 0.06387 2.7172 8 #O-Cu
pair_coeff 3 4 lj/cut 0 0 8 #H1-H2
pair_coeff 3 5 lj/cut 0.03396 1.335 8 #H1-Cu
pair_coeff 4 5 lj/cut 0 0 8 #H2-Cu
pair_modify pair lj/cut mix geometric
bond_coeff 1 268.0000 1.5290
bond_coeff 2 320.0000 1.4100
bond_coeff 3 320.0000 1.4100
bond_coeff 4 340.0000 1.0900
bond_coeff 5 340.0000 1.0900
bond_coeff 6 340.0000 1.0900
bond_coeff 7 340.0000 1.0900
bond_coeff 8 553.0000 0.9450
bond_coeff 9 553.0000 0.9450
angle_coeff 1 50.000 109.500
angle_coeff 2 50.000 109.500
angle_coeff 3 37.500 110.700
angle_coeff 4 37.500 110.700
angle_coeff 5 37.500 110.700
angle_coeff 6 37.500 110.700
angle_coeff 7 55.000 108.500
angle_coeff 8 55.000 108.500
angle_coeff 9 35.000 109.500
angle_coeff 10 33.000 107.800
angle_coeff 11 35.000 109.500
angle_coeff 12 35.000 109.500
angle_coeff 13 35.000 109.500
angle_coeff 14 33.000 107.800
---------------------Minimization---------------------------
reset_timestep 0
min_style cg
minimize 1.0e-13 1.0e-12 1000 10000
print “Finished Minimizing”
---------------------EQUILIBRATION--------------------------
reset_timestep 0
timestep ${dt}
velocity all create 298.0 97563
fix 1 solid nvt temp 298 298 (100.0*dt) tchain 1
fix 2 liquid nvt temp 298 298 (100.0*dt) tchain 1
compute 1 solid temp
compute 2 liquid temp
Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp density c_1 c_2
#dump 1 all atom 1000 dump.cunp.lammpstrj
Run for at least 100 picosecond (assuming 0.04 fs timestep)
run 2500000
unfix 1
unfix 2
---------------------Production-----------------------------
reset_timestep 0
timestep ${dt}
velocity all create 298.0 54664 mom yes rot no
fix 1 all npt temp 298.0 298.0 (100.0*dt) iso 1 1 100
#fix 2 liquid npt temp 298.0 298.0 (100.0*dt) iso 1 1 100
Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp density c_1 c_2
run 100000
unfix 1
#unfix 2
------------Thermal_Conductivity_Calculation----------------
reset_timestep 0
timestep ${dt}
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all centroid/stress/atom NULL virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix 1 all nve
fix JJ all ave/correlate $s $p d &
c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running
variable scale equal {convert}/${kB}/$T/$T/$V*s*{dt}
variable k11 equal trap(f_JJ[3]){scale}
variable k22 equal trap(f_JJ[4])*{scale}
variable k33 equal trap(f_JJ[5])${scale}
thermo $d
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 temp density c_1 c_2
#dump TC_dump all custom 1000 TC_dump.dump id type x y z
run 100000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print “average conductivity: $k[W/mK] @ T K, {ndens} /A^3”
Best regards,
Yunes Salman
log.lammps (436 KB)