[lammps-users] different results from two kinds of method for calculating shear modulus of silicon

Hi, all

I am calculating the shear modulus of silicon using two different methods, which yield two different results. In this two method, the silicon structure, potential parameter and temperature are all the same. The Silicon I use is 848 unit cells .

The first method is using " fix deform erate" command as follows

boundary p p p

fix 2 all nvt 300.0 300.0 10.0
fix 3 all deform 1 xz erate 1e-3 remap v

The shear stress is pxz, the shear strain rate is 1e-3. So from the plot I can fit the slope to the shear modulus, which is around 67 Gpa

The second method is using the method in the example “shear”. You can find the input file in the example folder. Several key commands are as follows

boundary s p s

region lower block INF INF INF INF INF 3 units box
region upper block INF INF INF INF 40 INF units box
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary

velocity upper set 0.5 0 0 units box
velocity mobile ramp vx 0 0.5 z 3 40 sum yes units box

I get 30 Gpa.

I do not know why these two methods give different results. I would really appreciate it if you can give me some suggestions.



Impressive that you got such great agreement using these very different approaches. Note that the latter approach includes some boundary effects. You’ll probably see the answers converge in the limit of bigger and longer simulations.

This is more a simulation methodology question than a LAMMPS question. It would probably be useful for you to consult papers and books that deal with computing shear using MD.