I am unable to get right equation of motion for SW potential

LAMMPS LAMMPS Beginners
1

I am trying to get equation of state of silicon in sw potential for T =0.05 and density =0.445 but I am unable to get right pressure at these density could anyone tell me what’s wrong in my script.


#variable
variable   sigma equal 2.0951
variable   epsilon equal 2.1683

variable   tempScale equal 2.5173e4
variable   densityScale equal 5.0571
variable   T equal 0.05
variable   t equal  $(v_T*v_tempScale)
variable   dt equal 0.005
variable   rho equal 0.489
variable   a equal (8.0*$(v_sigma)*$(v_sigma)*$(v_sigma)/${rho})^(1.0/3.0)



#main script

units 	metal 
atom_style	atomic
boundary	p p p




lattice         custom $a               &
                a1 1.0 0.0 0.0          &
                a2 0.0 1.0 0.0          &
                a3 0.0 0.0 1.0          &
                basis 0.0 0.0 0.0       &
                basis 0.0 0.5 0.5       &
                basis 0.5 0.0 0.5       &
                basis 0.5 0.5 0.0       &
                basis 0.25 0.25 0.25    &
                basis 0.25 0.75 0.75    &
                basis 0.75 0.25 0.75    &
                basis 0.75 0.75 0.25

region          myreg block     0 4 &
                                0 4 &
                                0 4

create_box      8 myreg
create_atoms    1 region myreg &
                basis 1 1  &
                basis 2 2  &
                basis 3 3  &
                basis 4 4  &
                basis 5 5  &
                basis 6 6  &
                basis 7 7  &
                basis 8 8
                
mass     *    28.0855

velocity   all create $t 42342 loop geom

pair_style	sw
pair_coeff    * * Si.sw Si Si Si Si Si Si Si Si

neighbor	2.0 bin
neigh_modify  every 1 delay 0 check no

timestep	${dt}

thermo		10
thermo_style	custom step density ke pe etotal vol temp press

fix	1 all nvt temp $t $t 0.1



#output variable

compute   myTemp all temp
compute    Press all pressure myTemp


variable  curT equal temp
variable  curStep equal step
variable  curV equal vol
variable  curH  equal enthalpy
variable  curKE equal ke
variable  curPE equal pe
variable  curE  equal etotal
variable     curP equal (c_Press[1]+c_Press[2]+c_Press[3])/3.0
variable  curD equal density

fix     2 all print 1 "${curStep} ${curD} ${curKE} ${curPE} ${curE} ${curV} ${curH} ${curT} ${curP}" screen no file Si_rho${rho}_T0.05_metal.dat title "#Step Density KE PE TotalE Vol Enthalpy Temp Press"

dump  1 all custom 10 Si_rho${rho}_T0.05_.dump id type xu yu zu x y z vx vy vz
log    Si_rho${rho}_T0.05_.log

run  100000

Thanks

2

Hi @pankaj,

your question makes no sense.

  • On the one hand you apparently ask question using reduced units with T=0.05 and \rho=0.445 in you question but your script launch simulates a temperature of \approx{}1258\ \mathrm{K} and a density of \approx{}2.48\ \mathrm{g\cdot{}cm^{-3}} using `metal units and Stillinger-Weber potential. These value is quite dense for Silicon at high temperature so what’s wrong here? Why use reduced units in the first place? What did you expect as an output and what is the issue here? Which version of LAMMPS are you using? There are no information to help you in any meaningful way.
  • On the other hand “getting the equation of motion” makes no sense. LAMMPS integrates equation of motion at the given condition using velocity verlet algorithm (or respa) and you get the output. The pressure is what it is then. You should get first a proper understanding of what the hamiltoniant used in the NVT ensemble is and what the MD script you provided is doing.

But overall your script looks very complicated. I think this shows a lack of proper understanding on how LAMMPS and MD work in the first place. You might get more insightful advices from direct colleagues or advisors and having a look at the relevant literature in the first place.

3

yes ,right I am trying to plot Pressure vs dinsity plot for silicon at 1258K but I am not getting the right slope.And I dont know whats the error?

4

Your script is a constant volume (or constant density) simulation. I don’t know how you can expect to get any slope on density variation out of this.

5

i am checking pressure at different densities then plotting it as average pressure and density .

6

Does that mean you are trying to get data for fitting to an equation of state?
If yes, how do you determine what is “correct”?

7

I have reference data from literature where they have shown that in silicon if you plot pressure vs density you will get a non linear curve but i am getting linear curve.

8

Is this data for the specific Stillinger-Weber potential parameterization that you are using?

9

Its using same Si.sw files default values. In my code I my doing NVT simulation at different densities and taking output.

10

These are a lot of details you didn’t provide initially and are very different from your initial question.

“Equation of motion” and “equation of state” are two very different things.

11

oh sorry about that it was a mistake.

12

If you have problems reproducing data from a paper, then this is either a topic for discussion with your adviser/supervisor or with the author('s) of the publication. It is not really a LAMMPS issue.

I am going to “close” this topic (and the corresponding previous one on the same issue) since there is nothing left to do here from the LAMMPS side of things.

1 Like