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.