LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Wrong temperature and pressure output because of wrong DOF - ellipsoidal particles
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Wrong temperature and pressure output because of wrong DOF - ellipsoidal particles


From: Thi Lan Huong Nguyen <a1672879_student@...7249...>
Date: Tue, 21 Nov 2017 01:38:58 +0000

Dear Lammps experts,

I have a question regarding to the DOF of ellipsoid and simulations relating to this. The lammps version I am using is from 16Feb16

So I am trying to do a test simulation to equilibrate ellipsoidal particles at 300K and 1atm. I tried to use a soft potential to avoid overlapping before equilibrating the system in npt. My ellipsoids are 3D particles with 3 different diameters.

As I understand from the Lammps documentations, fix nvt/asphere and fix npt/asphere allow thermostat to be applied to both translational and rotational motion, so as long as I use this for my system (regarding to the fact that my ellipsoids should have 6 dof) I don't have to use any fix_modify or compute_modify to adjust the DOF.

However, when I didn't adjust the DOF, the temperature output was doubled (~600K). Then when I tried to subtract 3 DOF from each particle the temperature was correct. Below this message is the script with modified DOF (by adding compute_modify to modify the temp in the fixes).

I do not understand why subtracting DOF helps correcting the temperature. Any idea for this?

I would appreciate any help for this problem.


#Input file for simulation with modified DOF

#----------Initialisation----------
clear
units         real
atom_style   ellipsoid
dimension    3

boundary     p p p

#----------Create atoms----------
lattice      sc 16                                
region    lower block 0 32 0 32 5 14                  

create_box    1 lower
create_atoms    1 region lower

#---------Atom Settings-----------
set         type 1 mass 1021.59
set         type 1 shape 32.9 15.9 3.5  
set         group all quat/random 18238

compute         q all property/atom quatw quati quatj quatk
compute      shape all property/atom shapex shapey shapez

dump      1 all custom 10 dump.soft.all &
         id type x y z c_q[1] c_q[2] c_q[3] c_q[4] &
          c_shape[1] c_shape[2] c_shape[3]  
 
######################---------Soft potential to prevent ovelapping---------
pair_style soft 33.0
pair_coeff 1 1 0.0 5.4


group         spheroid type 1
variable     dof equal count(spheroid)*3
velocity     all create 300 87287 loop geom

variable prefactor equal ramp(0,100)
fix 1 all nvt/asphere temp 300 300 50
fix 2 all adapt 1 pair soft a * * v_prefactor
compute_modify 1_temp extra ${dof} #this gives T=300, without it T=600

timestep 0.5

thermo 100
thermo_style custom step pe etotal vol press temp

write_restart restart.soft
run 500000

########################----------Real Force Fields ----------

pair_style gayberne 1.0 1.0 1.0 16.0
pair_coeff   1 1 0.1 3.5 2 2 10 0 0 0

neighbor     33.0 bin
neigh_modify every 5 delay 0 check yes


################------------- Run NPT ----------------------------
reset_timestep 0

thermo_style custom step epair etotal press vol temp density
thermo         100

timestep     0.1

#----------Visualisation----------

uncompute    q
uncompute    shape
compute         q all property/atom quatw quati quatj quatk
compute      shape all property/atom shapex shapey shapez

undump         1
dump         1 all custom 100 dump.equilibrium.ellipsoid &
         id type x y z c_q[1] c_q[2] c_q[3] c_q[4] &
          c_shape[1] c_shape[2] c_shape[3]
dump         trj all atom 100 dump.equilibrium.lammpstrj

unfix 1
unfix 2
fix          1 all npt/asphere temp 300 300 10 iso 1.0 1.0 100
compute_modify 1_temp extra ${dof}
run         100000