Questions

Dear all,
I am trying to implement to apply cyclic shear on a polycrystalline materials. I have used
fix NVT/SLLOD and fix deform command. The box deformation rate is time dependent.
The system is sheared by 1000 MD steps along +x direction and then under sheared condition another 1000 MD steps are simulated with the same NVT/SLLOD (Is it right to simulate the system under NVT/SLLOD when the box is not deforming anymore?). Then with the same rate the system comes back to the original position and a cycle ends with another 1000 MD steps under unsheared condition.

First question:
When MSD of the system is calculated it comes out that at a certain low strain rate (consequently lower maximum strain) MSD is becoming larger as compared to the higher strain rate. It sounds quiet strange to me. If you make some comment it would be good.

Second question:
In some cases, I found the translational drift of the centre of mass. To remove this problem I used fix momentum command. How much logical it is? What is the reason of COM movement?

Third question:
When I calculate the velocity of the top and the bottom layer (which is defined by region command) without using the fix momentum command it shows that the bottom layer has a small velocity along -x direction whereas one expects that its velocity should fluctuate around zero.
However, if I sum up the speed of the top and lower layer it matches with box deformation rate.
Why does the lower layer move with a certain velocity (-x direction)? How can I fix the problem?

The input file is shown below.

units lj
atom_style atomic
dimension 2

boundary p p p

read_data system.data

variable DT equal 0.005
timestep ${DT}

region upper block INF INF 104.409 INF INF INF
region lower block INF INF INF 1.0 INF INF

pair_style lj/cut 3.5 #LJ potential with cut off 3.5

pair_coeff 1 1 1.0 1.0 3.0 #First 1 1 says that the interaction between "1" type of atoms; Epsilon=1.0; Sigma=1.0; Cutoff=3.0
pair_coeff 2 2 1.0 1.8 5.4

pair_coeff 1 2 1.0 1.4 4.2 #First 1 2 says that the interaction between "1" and "2" type of atoms; Epsilon=1.0; Sigma=1.4; Cutoff=4.2;

pair_modify shift yes
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes

variable SS equal 1000
variable S1 equal 500
variable S2 equal 500
variable S3 equal 500
variable S4 equal 500
variable S equal "v_S1+v_S2+v_SS”
variable SR equal "v_S+v_S3"

variable f equal 0.1

reset_timestep 0

variable u equal 0.0
variable rate equal "v_f*(step*dt)" # deformation rate for the first 500 MD steps
variable displace equal "(v_f*((step*dt)*(step*dt)))/2" # Displacement of the box during the deformations
fix 2 all deform 1 xy variable v_displace v_rate remap v # Command to deform the box and v_rate and

fix 1 all nvt/sllod temp 0.01 0.01 0.1 tchain 3
compute myTemp all temp/deform
fix_modify 1 temp myTemp

compute MSD all msd com yes
fix ave_MSD all ave/time 100 1 100 c_MSD[1] c_MSD[2] c_MSD[3] c_MSD[4] file MSD.data

compute ave_vel all reduce ave vx vy vz
compute sum_vel_up all reduce/region upper ave vx vy vz
compute sum_vel_low all reduce/region lower ave vx vy vz
fix ave_vel1 all ave/time 1 1 1 c_sum_vel_up[1] c_sum_vel_up[2] c_sum_vel_low[1] c_sum_vel_low[2] v_rate c_ave_vel[1] file vel.dump

run ${S1} #by this steps box deformation rate reaches the maximum

variable u equal "v_f*v_S1*dt" #initial velocity when the box deformation rate decreases
variable rate equal "v_u-v_f*((step-v_S1)*dt)" # box deformation rate decreases
variable displace equal "(v_u*((step-v_S1)*dt)-(v_f*(((step-v_S1)*dt)*((step-v_S1)*dt)))/2)" # displacement

run ${S2} #by this MD steps period box deformation rate comes back to zero

Sorry, this script is too complicated to easily understand. The drop in MSD at higher strain rates might be caused by slip occurring near the walls. The best way to understand what is happening in your simulation is to vizualize a small system.

Also, in order to achieve a smoother oscillation, you could try using fix deform xy wiggle.

I agree that one can use fix deform xy wiggle for smoother oscillation. However, I need to relax the system under sheared condition which can not be done by “fix deform xy wiggle”. Therefore, I need to do “fix deform xy variable”. The way what I use to have the smoother shearing is as follows: First 500 MD steps the box is deforming with the rate ft where f is a constant and t is calculated as (stepdt). In the next 500 MD steps, box deformation rate slowly decreases to zero. Then under the sheared condition the system is asked to stay for 1000 MD steps where fix deform and fix NVT/SLLOD commands are still working.
Input script for those steps are here.

In the same way the system comes back to the original position.

Now, At the end of every cycle a translational drift of the COM is observed. To fix the problem I used “fix momentum linear” command. Consequently, the particles close to the upper and the lower part of the box start to move with the same speed and opposite direction as fix momentum command keeps linear momentum fixed by rescaling the velocity. Therefore, the particles are not exactly following the box deformation rate. When I checked the movie it looks that the particles close to the lower part of the box moves from one side to another due to the periodic boundary condition. However when unwrapped coordinates are checked it shows nicely that the upper particles are moved +x direction and lower particles are moved -x direction. Do you think that it is wrong as the particles are not following the box deformation rate due to the velocity rescaling by fix momentum command?

With best regards,
Pritam

Then under the sheared condition the system is asked to stay for 1000 MD steps where fix deform and fix
NVT/SLLOD commands are still working.

I don’t know what this means. If the box is stationary, then fix deform can (should) be off.

E.g. use the unfix command, then run for 1000 steps.

Likewise if the system is not deforming you should not be using fix nvt/sllod.

Unfix it and use fix nvt or fix nve for the 1000 step run.

Steve

Dear Steve,
Thank you for the suggestion. I agree that one should not use NVT/SLLOD when the box is stationary even in the deformed condition. So, I change the input script and now it is

variable f equal 0.1

variable u equal 0.0
variable rate equal "v_f*(stepdt)"
variable displace equal "(v_f
((stepdt)(step*dt)))/2”

fix 2 all deform 1 xy variable v_displace v_rate remap v
fix 1 all nvt/sllod temp 0.01 0.01 0.1 tchain 3

compute myTemp all temp/deform
fix_modify 1 temp myTemp

run 500

variable u equal “v_f500dt”

variable rate equal “v_u-v_f*((step-500)dt)"
variable displace equal "(v_u
((step-500)dt)-(v_f(((step-500)dt)((step-500)*dt)))/2)”
run 500

unfix 1
unfix 2
uncompute myTemp

fix NVT all nvt temp 0.01 0.01 0.1 tchain 3
compute myTemp all temp
run 1000

The system still has the translational drift after every cycle either in the +x or in the -x direction. Do you have any suggestion regarding this problem?
The system is consist of two different types of particles. They have same mass but different sizes. The number of bigger particles is 25 whereas the smaller is 10,000. The system size is 105.409*105.409.
Is there any limitation of minimum tilt factor one can use by fix deform command?

With best regards,
Pritam

If you shear/drive a system in one direction,

then turn off the shear, it doesn’t mean it will

stop drifting.

Steve

Thank you very much. Do you think one can use fix momentum command to keep COM fixed by moving particles close to upper and lower boundaries in opposite directions though the lower part of the box remains fixed?

With best regards,
Pritam

Thank you very much. Do you think one can use fix momentum command to keep
COM fixed by moving particles close to upper and lower boundaries in
opposite directions though the lower part of the box remains fixed?

the most straightforward approach to keep the COM of a group of
particles from drifting would be to use fix spring with the tether
option.
whatever method you are using, you need to make certain that it
doesn't taint the results of your simulation more than the error due
to other manipulations that you are already performing.

axel.

Thank you very much. To apply strain on the system, I have used fix NVT/SLLOD and fix deform. I have defined the upper layer by region command to compute the average velocity of the particles close to the upper part of the box. I was expecting that upper layer should follow the box deformation rate and lower region which is close to the bottom of the box should be fixed. However, I find some deviation from the expectation. Here, I have attached the figure. If you make some comment it would be nice.

Best regards,
Pritam

vel_f_0-1.pdf (39.4 KB)

There are many issues with how to use fix deform and fix nvt/sllod. Like
how fast you shear, how big your box is, whether the system is a liquid or solid,
You say polycrystalline, so I assume a solid. Nvt/sllod is really designed for
liquids, which flow with the deforming box, which is not the case for solids.
Hence the requirement that you use fix deform remap v with nvt/sllod. See
the discussion of the remap keyword on the fix deform doc page.

I suggest you first focus on getteing your system to do the right thing while the box is deforming
before worrying about what happens after you have deformend it.

Steve

Thank you for the reply.

One more thing I am trying to understand that if I use fix NVT instead of NVT/SLLOD, then same type of behaviour is observed. However, if I use fix langevin, the system shows a nice behaviour, i.e. the velocity of lower layer is fluctuating around zero and top layer also follows the box deformation rate. Is it like the problem comes from the thermostat?

It will be extremely good if you suggest some other ways to perform shear in the bulk (i.e. with periodic boundary condition)? Because I found fix deform is very convenient to use with thermostat when we need to shear in the bulk. My system is prepared by cooling down the mixture of Lennard Jones particles where the final structure has some defects.

With best regards,
Pritam