Temperature gradient problem

LAMMPS LAMMPS General Discussion
1

dear all,
I am trying to calculate the thermal conductivity of CO2 using NEMD by employing langevin but after simulation, I was getting temperatures far less and far higher than the temperatures I set. Can anybody help me solve this? Below is my input script

units		real	  
atom_style	full      

timestep	1.0		

boundary        p p p

molecule co2 co2.mol

region box block 0 30 0 30 0 200
region box1 block 0 10 0 10 0 20
region box2 block 0 10 20 30 0 20
region box3 block 0 10 20 30 180 200
region box4 block 0 10 0 10 180 200 
region box5 block 20 30 0 10 0 20
region box6 block 20 30 0 10 180 200
region box7 block 20 30 20 30 0 20
region box8 block 20 30 20 30 180 200
create_box 2 box &
	bond/types 1 &
	angle/types 1 &
	extra/bond/per/atom 2 &
	extra/angle/per/atom 1 &
	extra/special/per/atom 2
	
create_atoms 0 random 20 456415 box1 mol co2 554756
create_atoms 0 random 20 456415 box2 mol co2 554756
create_atoms 0 random 20 456415 box3 mol co2 554756
create_atoms 0 random 20 456415 box4 mol co2 554756
create_atoms 0 random 20 456415 box5 mol co2 554756
create_atoms 0 random 20 456415 box6 mol co2 554756
create_atoms 0 random 20 456415 box7 mol co2 554756
create_atoms 0 random 20 456415 box8 mol co2 554756

mass 1 12.01070	
mass 2 15.99990

pair_style	lj/cut/coul/long 14.0 14.0  
pair_modify     tail yes			    
bond_style      harmonic
angle_style     harmonic
kspace_style	pppm 1.0e-4             
#kspace_style	ewald 1.0e-4            
                                   
group           carbon type 1
group           oxygen type 2

pair_coeff      1 1 0.053654828567 2.8000       
pair_coeff	2 2 0.156990053954 3.05000          
pair_coeff	1 2 0.091778398502 2.9250          

bond_style	harmonic   			     
angle_style	harmonic			     
dihedral_style	none				 
improper_style	none				 

bond_coeff	1 5000.00 1.160                     
angle_coeff	1 500.0 180.00			   
		  				    

#Neighbor List Parameters
neighbor        2.0 bin			     
neigh_modify    delay 0 every 1 check yes
                                                                                         
thermos 10    
min_style sd
minimize	1.0e-4 1.0e-6 1000 10000

reset_timestep 0

velocity all create 300 4932567 mom yes rot yes dist gaussian


region          cold block INF INF INF INF 0 10
region          hot block  INF INF INF INF 100 110
compute         Thot all temp/region hot
compute         Tcold all temp/region cold
variable        tdiff equal c_Thot-c_Tcold

variable        kinetic_energy equal ke
variable        potential_energy equal pe
variable        pressure equal press   
variable        Temperature equal temp   
variable        Density equal density    
variable        total_energy equal etotal

# 1st equilibration run

fix             1 all nvt temp 300 300 100 
thermos          100
thermo_style    custom step temp density epair ke etotal press 
run            300000

velocity        all scale 300
unfix 1

# 2nd equilibration run

fix             1 all nve
fix             hot all langevin 320 320 100 59804 tally yes
fix             cold all langevin 280 280 100 287859 tally yes
fix_modify      hot temp Thot
fix_modify      cold temp Tcold
thermo 1000
thermo_style    custom step temp press density epair ke etotal 
dump 1 all custom 100 dump1.lammpstrj id type x y z vx vy vz
run 3000000

undump 1

compute         ke all ke/atom
variable        temp atom c_ke*1000/(1.5*1.98721)

fix             hot all langevin 320 320 100 59804 tally yes
fix             cold all langevin 280 280 100 287859 tally yes
fix_modify      hot temp Thot
fix_modify      cold temp Tcold


fix             ave all ave/time 10 100 1000 v_tdiff ave running
fix             3 all ave/time 10 100 1000 v_Temperature v_kinetic_energy v_potential_energy v_total_energy v_pressure v_Density file ave.dens_18.0molL30.out format %.8g
thermo_style    custom step temp c_Thot c_Tcold v_tdiff f_hot f_cold f_ave press density epair ke etotal vol lx ly lz
#thermo_modify   colname c_Thot Temp_hot colname c_Tcold Temp_cold &
                colname f_hot E_hot colname f_cold E_cold &
                colname v_tdiff dTemp_step colname f_ave dTemp

compute         layers all chunk/atom bin/1d z lower 0.05 units reduced
fix             2 all ave/chunk 1 600000 600000 layers v_temp file profile.langevin

variable start_time equal time

variable kappa equal (0.5*(abs(f_hot)+abs(f_cold))/(time-${start_time})/(lx*ly)/2.0)*(lz/2.0)/f_ave
thermo  10000
run             6000000
print           "Running average thermal conductivity: $(v_kappa:%.10f)"
Summary

This text will be hidden

2

There are multiple issues here:

  • Your posted input has syntax errors and thus cannot even run without modifications. So how can you have tested this?
  • Your system is TINY with less than 500 atoms. How can you expect to get meaningful and converged results without large fluctuations from such a tiny system?
  • Have you checked whether your method does even work (in theory) for diluted/gas phase systems? Many methods for material properties require condensed systems.
  • Have you checked your system and simulation settings for whether they can properly conserve energy after being properly equilbrated?

Overall, it very much looks like you need to get your simulation fundamentals done properly before even thinking about applying some advanced simulation methods and also need to properly study and understand the methods you want to apply.

3

I was only trying to see if my input script was correct. so can you please show me the syntax errors

4

LAMMPS can do this by itself as it will print an error and stop on syntax errors.
How to correctly use LAMMPS is described in the documentation.

5

The input scripts did work so I tried a real simulation. During equilibration in NVT the average pressure was very close to the expected pressure, but I realized that after combining Langevin command with NVE, the average pressures were far far smaller than the expected pressure.
is there any explanation for that?

7

The quoted input cannot be executed, since it contains the command “thermos” which does not exist in LAMMPS.

This is a discussion about the science of your simulation and thus something you need to discuss with your adviser/tutor/mentor or similar; i.e. somebody that knows and cares about your research. As I already stated, I would not expect meaningful results at all from your input.

At any rate, this kind of discussion is off-topic for the LAMMPS forum categories. They are not LAMMPS issues.