[lammps-users] Fix in variable not computed at compatible time (../variable.cpp:4128)

Dear All,

I have lammps 29Sep21 installed.

I’m trying to calculate the thermal conductivity of my nanofluid using Green-Kubo method.

When I run my script I get this error “ERROR: Fix in variable not computed at compatible time (…/variable.cpp:4128)”.

I understand that this error has to do with frequency of output mismatch.

In my input script this is has to do with the dt and s variables.

If I do dt=10 and s=1 my simulation runs and I get rid from this error. However, it screws up my simulation because my time step needs to be 0.04

Do you know a work around or how I can fix this without changing my time step?

My input script looks like this:

------------------------ 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 300
variable V equal vol
variable dt equal 0.04
variable p equal 200 # 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 ----------------------------

region whole block -50 50 -50 50 -50 50 units lattice
create_box 4 whole bond/types 9 angle/types 14 extra/bond/per/atom 9 extra/angle/per/atom 14 extra/special/per/atom 10
lattice fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

molecule ethylene_glycol lammps.data
create_atoms 0 region whole mol ethylene_glycol 52874 ratio 0.035 74637

region nanoparticle sphere 0 0 0 3.5 units lattice
lattice fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 4 region nanoparticle

mass 1 12.011 #O
mass 2 15.999 #C
mass 3 1.008 #H
mass 4 63.546 #Cu

------------------------ FORCE FIELDS ------------------------------

#pair_style hybrid eam/fs lj/cut 2.5
pair_style eam/fs
pair_coeff * * Mendelev_Cu2_2012.eam.fs NULL NULL NULL Cu

pair_style lj/cut 2.5
#pair_coef type_i type_j epsilon sigma rcut
pair_coeff * * 0.030 2.5 2.5 #all
pair_coeff 2 3 0.030 2.5 2.5 #H-C
pair_coeff 2 2 0.066 3.5 2.5 #C-C
pair_coeff 1 3 0.170 3.07 2.5 #O-H
pair_coeff 3 4 0.03396 1.335 2.5 #H-Cu
pair_coeff 1 4 0.06387 2.7172 2.5 #O-Cu

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

info system
print “Finished INFO”

---------------------Minimization---------------------------

reset_timestep 0
minimize 1.0e-13 1.0e-12 1000 10000
print “Finished Minimizing”

---------------------EQUILIBRATION--------------------------

reset_timestep 0
timestep ${dt}

velocity all create 300.0 97563 mom yes rot no
fix 1 all nvt temp 300.0 300.0 $(100.0*dt) tchain 1

Set thermo output

thermo 1
thermo_style custom step lx ly lz press pxx pyy pzz pe temp

run 10
unfix 1

---------------------Production-----------------------------

reset_timestep 0
timestep ${dt}

velocity all create 300.0 54664 mom yes rot no
fix 1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0

Set thermo output

thermo 1
thermo_style custom step lx ly lz press pxx pyy pzz pe temp

run 10
unfix 1

------------Thermal_Conductivity_Calculation----------------

variable dt equal 10
reset_timestep 0

compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all 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 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 1
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33

Run for at least 100 picosecond (assuming 0.04 fs timestep)

dump 1 all custom 1 dump.comp_*.dump mass type x y z
run 10
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”

Outlook-glxkiqb5.png

Outlook-glxkiqb5.png

Well, debugging this should be straightforward.
You just look at where you have values that are derived from the length of the time step and adjust them so you get the values you need.
LAMMPS helps you with that as it will echo to the log file exactly how variables in the input are expanded.

The rest is up to you.

Axel.

Outlook-glxkiqb5.png

Outlook-glxkiqb5.png