Dear Lammps users,
Hi, I am trying to calculate thermal conductivity of bulk silicon using Green-Kubo approach. I was able to get the conductivity of ~200 W/m-K as setting up the structure within the input file (lattice, region, creeate_box, etc…) with sw potential. Now I am trying to use ‘read_data’ command with the data file containing silicon atom positions in 6 x 6 x 6 unit cells, to see its isotopic effect.
However, I keep getting ‘nan’ for Temp, Epair, TotEng, … in the output file, and the simulation ended up with an error message " Invalid thermo keyword in variable formula"
I wonder why this happened with ‘read_data’ command despite the successful simulation with the material structure defined within the input file.
I am attaching the input file, data file, and output file. Please let me know if you can see that I am wrong on doing this.
Thank you in advance,
%-------------input file-------------------------
units metal
variable T equal 300
variable V equal vol
variable dt equal 0.001
variable p equal 10000 # correlation length
variable s equal 5 # sample interval
variable d equal $p*$s # dump interval
convert from LAMMPS real units to SI
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.602e-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal {eV2J}*{eV2J}/{ps2s}/{A2m}
setup problem
dimension 3
boundary p p p
#lattice diamond 5.43 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
#region box block 0 8 0 8 0 8
#create_box 1 box
#create_atoms 1 box
atom_style atomic
read_data siData.meam
#mass 1 28.0855
pair_style sw
pair_coeff * * Si.sw Si
timestep ${dt}
thermo $d
equilibration and thermalization
velocity all create $T 4928459 mom yes rot yes dist gaussian
fix NPT all npt temp $T $T 1 x 0.0 0.0 10 y 0.0 0.0 10 z 0.0 0.0 10
run 2000000
thermal conductivity calculation, switch to NVE if desired
unfix NPT
fix NVE all nve
reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom 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 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_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run 1000000
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”
%-------------data file-------------------------
Position data for Bulk Silicon
1963 atoms
1 atom types
0 32.5800 xlo xhi
0 32.5800 ylo yhi
0 32.5800 zlo zhi
Masses
1 28.0855
Atoms
1 1 0 0 0
2 1 5.43000 0 0
3 1 5.43000 5.43000 0
4 1 5.43000 0 5.43000
5 1 5.43000 5.43000 5.43000
6 1 0 5.43000 0
7 1 0 5.43000 5.43000
%-------------output file-------------------------
LAMMPS (10 Jan 2014)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (0 0 0) to (32.58 32.58 32.58)
3 by 4 by 4 MPI processor grid
1963 atoms
Setting up run …
Memory usage per processor = 2.09671 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -nan 0 -nan -nan 34582.25
50000 -nan -nan 0 -nan -nan -nan
100000 -nan -nan 0 -nan -nan -nan
150000 -nan -nan 0 -nan -nan -nan
200000 -nan -nan 0 -nan -nan -nan
250000 -nan -nan 0 -nan -nan -nan
300000 -nan -nan 0 -nan -nan -nan
%----------------ends----------------------------------