Dear Lammps experts,

I am simulating impact of a particle on a similar type substrate with LAMMPS . my problem is about calculating the particle temperature as a moving object actually.

First the initial structure of the particle and its energy minimize at the room temperature(300K) with fix nve and fix temp/rescale and then the particle is thrown toward the substrate. I applied the compute temp/com in order to calculate the particle temperature as a moving object as per the Lammps document.

These are the problems:

1- Exactly upon throwing the particle ( after unfix temp/rescale) the particle temperature decreases about 150K which is not expected. I have no idea why temperature decreases TOO MUCH upon unfixing. How about minimizing the particle at the room temperature longer???

2- and then the temperature increases until the end of the simulation. However, it is expected that the temp. only increases at the impact moment and then decreases to a steady or constant temp at the end of the simulation.

So, any comment could be highly appreciated.

below you can find my input file:

#Phase 1 ------------------------------------------Simulation main setup----------------------------------

dimension 3

units real

atom_style charge

boundary p p p

variable radius equal 50 # makes 100 Ang or 10nm particle

variable n equal 50 # makes 81*3.786 length=306

variable a1 equal 3.786

variable a2 equal 3.786

variable a3 equal 9.514

variable boxx equal “v_n*v_a1"
variable boxy equal "v_n*v_a2”

variable diameter equal “2*v_radius"
variable x0 equal "0.5*v_boxx”

variable y0 equal “0.5

*v_boxy"*

variable subs_thick equal "3v_a3” #v_a3 equal 9.514

variable subs_thick equal "3

variable subslow_thick equal “1

*v_a3"*

variable z_gap equal "0.5v_radius”

variable z_gap equal "0.5

variable distance equal “v_subs_thick+v_subslow_thick+v_z_gap+v_radius”

variable boxz equal “v_distance+v_diameter+v_radius”

#simulation box

region box block 0 {boxx} 0 {boxy} 0 ${boxz} units box

create_box 2 box #2 is number of atoms

# it was very difficult to finally find these positions

lattice custom 1 a1 3.786 0.00000 0.00000 a2 0.0000 3.786 0.00000 a3 0 0 9.514 &

basis 0.0000 0.2500 0.3750 &

basis 0.0000 0.7500 0.6250 &

basis 0.5000 0.7500 0.8750 &

basis 0.5000 0.2500 0.1250 &

basis 0.0000 0.0000 0.1700 &

basis 0.0000 0.7500 0.4200 &

basis 0.5000 0.2500 0.3300 &

basis 0.5000 0.7500 0.0800 &

basis 0.500 0.5000 0.6700 &

basis 0.500 0.2500 0.9200 &

basis 0.000 0.7500 0.8300 &

basis 0.000 0.2500 0.5800

mass 1 47.86000

mass 2 15.99940

#particle

region particle sphere {x0} {y0} {distance} {radius} units box

create_atoms 2 region particle &

basis 1 1 &

basis 2 1 &

basis 3 1 &

basis 4 1 &

basis 5 2 &

basis 6 2 &

basis 7 2 &

basis 8 2 &

basis 9 2 &

basis 10 2 &

basis 11 2 &

basis 12 2

group particle region particle

set type 1 charge 2.196

set type 2 charge -1.098

#substrate

region substrate block 0 {boxx} 0 {boxy} {subslow_thick} {subs_thick} units box

create_atoms 2 region substrate &

basis 1 1 &

basis 2 1 &

basis 3 1 &

basis 4 1 &

basis 5 2 &

basis 6 2 &

basis 7 2 &

basis 8 2 &

basis 9 2 &

basis 10 2 &

basis 11 2 &

basis 12 2

group substrate region substrate

set type 1 charge 2.196

set type 2 charge -1.098

#group model union particle substrate

#lower_substrate

region lower_substrate block 0 {boxx} 0 {boxy} 0 ${subslow_thick} units box

create_atoms 2 region lower_substrate &

basis 1 1 &

basis 2 1 &

basis 3 1 &

basis 4 1 &

basis 5 2 &

basis 6 2 &

basis 7 2 &

basis 8 2 &

basis 9 2 &

basis 10 2 &

basis 11 2 &

basis 12 2

group lower_substrate region lower_substrate

set type 1 charge 2.196

set type 2 charge -1.098

group model union particle substrate

#–Phase 2----------------------------------------Buckingham Potential-----------------------------------------------

pair_style buck/coul/long 15

pair_coeff 1 1 717647.40 0.154 121.067

pair_coeff 1 2 391049.10 0.194 290.331

pair_coeff 2 2 271716.30 0.234 696.888

neighbor 2.0 bin # skin distance for real units is by default 2.0

neigh_modify every 1 delay 0 check yes

kspace_style pppm 0.0001

neigh_modify delay 5

#—Phase 3----------------------------------------Compute—how can I monitor Vparticle!!! ------------------------

compute temp_substrate all temp/region substrate

compute temp_particle_com particle temp/com

compute temp_particle all temp/region particle

#----Phase 4-------------------------------------Initial Equilibration at 300K ----------------------------------------

reset_timestep 0

timestep 1.0 # or 2

velocity all create 300 12345 mom yes rot no

fix 1 lower_substrate setforce 0.0 0.0 0.0

fix 2 particle nve

fix 3 substrate nve

fix 4 substrate temp/rescale 100 300.0 300.0 10.0 1.0

fix 5 particle temp/rescale 100 300.0 300.0 10.0 1.0

thermo 100

thermo_modify lost ignore flush yes

thermo_style custom step c_temp_particle c_temp_particle_com c_temp_substrate

dump 1 all xyz 100 dump1000.txt

run 7000

unfix 5

#----Phase 5---------------------------------------Particle Impact at the room temperature -------------------

velocity particle set 0 0 -0.010 units box

thermo 100

thermo_style custom step c_temp_particle c_temp_particle_com

run 50000

Best regards

Bahman Daneshian