H diffuse in Pd

[yanxon]

Hi @timfrol,

I hope your vacation went well.

I created this repo in an attempt to calculate the diffusivity constant.

Please check out the LAMMPS scripts.

[yanxon]

Hi @timfrol,

1. I run diffusion calculation as we discussed before. Here is the average of 100 md run with LAMMPS: https://github.com/yanxon/diffusion/blob/master/results/100run.png

2. I also wrote my own Python script to track the hydrogen position in PBC condition. I validate the my MSD result with LAMMPS and they are identical. Here is the example for 1 md run with my script: https://github.com/yanxon/diffusion/blob/master/results/1run.png

3. I tried to implement your suggestion by shifting the reference position (r(t=0)). I don't think this method is reliable. For example, I use the result from #2 and shift the reference point. Here is what I obtain: Indeed, it smoothens the figure #2. But I don't think it obeys Einstein's very consistently.

[timfrol]

@yanxon Howard, great stuff. My van leaves at 5:20, not sure if there is enough time to chat today. We can certainly do it tomorrow.

I assume that you compared your own MSQD with the one from lammps and they match perfectly. This is very good. I would like to look at your X, Y, Z positions plots from which you generated MSQD. How many times did you shift the reference points. Obviusly it is a question of statistics. You need to increase the averaging until you get the straight line.

We should definitely chat. Looks very promising. Tim

[timfrol]

@yanxon

your 1. looks pretty linear, which is nice. I would try improving 3 and compare the diffusion coefficients you get. Definitely fit 1. with a straight line. Here is an example from one of my papers.

[yanxon]

@timfrol,

We can definitely chat tomorrow.

I shift it 1000 times. Here is where I shift 4000 times:

I have the script here: https://github.com/yanxon/diffusion/blob/master/test_script/shift.py

The script is messy. But it works. I would like to clean up so that it can be read.

[timfrol]

@yanxon

ok, I think you got it. Do as many averaging you can with this data. I would shorten the trajectory to the time interval from 0 to 0.2. Make all the averaging you can and then supperimpose it on your plot from 1.

Even now the section (0,0.2) begins to look pretty linear. I am wondering if it matches 1.

[yanxon]

@timfrol,

• I fit the linear equation into #1 (the title is the diffusion coefficient): https://github.com/yanxon/diffusion/blob/master/results/fit.png I also checked (eyeballed) my diffusion coefficient with the one in literature. They are in the same order of magnitude.
• I also try to get the statistics in the time interval of [0, 0.2] by shifting the msd 8000 times. https://github.com/yanxon/diffusion/blob/master/results/result.png The D coefficient seems like twice as high.

[timfrol]

@yanxon

sounds good Howard, do the [0,0.2] interval and let's chat after that

[timfrol]

@yanxon

Just to clarify

in 3d =6Dt 2d: =4Dt 1d =2Dt

When you compare the diffusion coefficient, you have to find one for the potential you are using, not the experiment or some other model.

[yanxon]

@timfrol

In my case, I use 6Dt throughout.

[0,0.2] interval: https://github.com/yanxon/diffusion/blob/master/results/result.png

FYI, the this paper (https://www.sciencedirect.com/science/article/pii/S1359646218300824) also using the same Potential. Therefore, I can compare.

[timfrol]

@yanxon Howard,

the plot looks good. Can you overlay it with the plot you got by averaging lammps MSQD over 100 separate simulations? Do you get the same number for Ds?

[timfrol]

@yanxon

Howard,

it is good that you are looking at this paper. A good test of the methodology would be to take their system size and concentration and repeat the calculation to make sure you get the same number.

Tim

[yanxon]

@timfrol

I put the D coefficient in the title of the plot. Sorry if it's too small.

[yanxon]

@timfrol

How do I tell LAMMPS that I want Pd10H4?

[yanxon]

@timfrol

I think figure it out. So, tell if I'm wrong. In the paper, they use 2048 Pd atoms. I just have to input ~819 hydrogen atoms. Is this correct?

[yanxon]

@timfrol

1. I ran a NPT calculation just like the paper suggested, and used the same parameters, number of atoms, etc. The figure is the result of PdHx (x=0.4) at 600 K. D = 1853.2 A^2/ns. It's close to the result in the paper which is ~2000 A^2/ns. I only ran the MD simulation for 0.15 ns, whereas Zhou et al. ran it for 2.2 ns. I think my method is correct.

2. I also ran NVT on the same system. But, they yield very different result on the diffusion coefficients, almost 1 order of magnitude. Can you please explain why does this happen?

[yanxon]

@timfrol

I think the difference can be explained by pressure. With NVT simulation, the pressure is at about 20 GPa. Whereas, the NPT simulation is forcing the pressure of the system to be at 0 GPa.

[timfrol]

@yanxon

very good Howard. You are becoming a master of MD. That tells me that you did not scale the simulation block correctly. Ideally, your scaled block should have had the same dimensions as the averaged dimensions in NPT. 20Gpa is a lot, you must have been off by a lot.

Did you even confirm that MSQD reported by LAMMPS gives you the same answer as the D you get from the positions and averaging by shifting the reference time?

[timfrol]

@yanxon

"Did you even confirm"

was a typo

I meant to write did you ever confirmed

[yanxon]

@timfrol

I concluded that they are not the same.

[timfrol]

@yanxon

Howard, if you tested your MSQD from atomic positions against MSQD from LAMMPS and you get identical result, there are only two possibilities for the discrepancy 1. Something is not right about how you do the shifting and the averaging. 2. If 1. is correct then it is a question of statistics.

The blue line is taken from a single trajectory in a single simulation. If you take a different simulation you should get a very different 'blue' line. By averaging them you will get the red line.

The shifting the reference averaging is one of the effective methods to improve your statistics. We should figure this out before we proceed to GB simulations.

[yanxon]

@timfrol

I just finished generating this graph: At a first glance, I think the average blue curve is populating towards the bottom of the red curve. But the average of the blue curve matches the red curve:

[timfrol]

@yanxon

Ok, Howard,

it looks pretty good to me and I feel you understood how this averaging by shifting the reference point as well as by running MD in parallel could improve the statistics for the D calculations. I suggest that you try running larger and longer jobs on borax, so we produce the data very quickly. At which point I would calculate D for a couple of other Ts to compare with the Arrhenius plot of D from the paper.

GBs after that and you can write your awesome paper.

[timfrol]

@yanxon

btw, once you run some simulations on borax I would like to repeat/change your runs just to make sure all the parameters and the equilibration looks good. We can exhange the files on lc by using give yanxon1 file, take frolov2 file commands.

[yanxon]

@timfrol

1. I will try to do this. I will use the same parameters as mentioned in the paper and reproduce figure 1(a).

2. I install lamps on borax and quartz

3. Do you want to meet with Brandon this Friday? I will send out an email to figure out time.

[timfrol]

@yanxon Howard, I have been having meetings each day in the past month. I need to do some writing and I would like to spend some time on the benchmark script. I would prefer not to have a meeting for a while. It sounds like we are making progress, I would focus on work for now. Feel free to have a meeting without me.

1. and 2. sounds great. You should have parallel lammps and the submission script for borax from last year right?

[timfrol]

@yanxon just gave you an expample of a simulation for borax. There should be a binary and the submission script

[timfrol]

@yanxon

just type on borax

take frolov2

[yanxon]

@timfrol

Sounds good. I think I will just postpone the meeting until I get back from the conference. I have nothing really to tell Brandon except I'm making progress on my MD training.

I do have the submission script for borax.

Thank you, @timfrol.

I wanna be done with this grain boundary project and move on with helping you and Qiang developing the benchmark.

[timfrol]

@yanxon

great!

[yanxon]

@timfrol

Can you please resend Cu01.eam.alloy? I accidentally delete it.

[yanxon]

@timfrol

Also, this script gives me error in the burn initial step part. As you suggested, we should not consider the first nanosecond of the calculation.

# ---------- Initialize Simulation --------------------- 
variable alat equal 3.914295757
variable xdim equal ${alat}*8
variable T equal 600
variable T2 equal ${T}
variable tstep equal 0.001
variable runtime equal 10000000
variable Srandomseed equal 2129

clear
units metal
dimension 3 
boundary p p p
atom_style atomic

# ---------- Create Atomistic Structure --------------------- 
lattice fcc ${alat}
region whole block 0 ${xdim} 0 ${xdim} 0 ${xdim} units box
create_box 2 whole
create_atoms 1 region whole
create_atoms 2 random 819 ${Srandomseed} NULL

group 1 type 1
group 2 type 2

# ---------- Define Interatomic Potential --------------------- 
pair_style  eam/alloy
pair_coeff * * PdH_zhou.eam.alloy Pd H

timestep ${tstep}

# ---------- set random velocity -----------
velocity all create ${T2} ${Srandomseed}

# Burn initial steps
fix 1 all npt temp ${T} ${T} 10 iso 0 0 10
run 10000
unfix 1

# ----- perform MD -------
fix 1 all npt temp ${T} ${T} 1 iso 0 0 1
compute msd 2 msd

#----- print ------
thermo 100 
thermo_style custom step temp c_msd[4] press
thermo_modify lost ignore flush yes

dump 1 2 custom 100 hydrogen/*.chkpt id type mass x y z

run ${runtime}

I just give you the potential on Borax.

[timfrol]

@yanxon

Howard,

I forgot how to run multiple simulations from one script, but you can search it on lammps manual. you may need to split the runs by using run 0 command or reset run. I am sorry I cannot remmenber.

[timfrol]

@yanxon

Howard,

you should use commands like

reset_timestep 0 thermo 0 run 0

to split the 2 simulations.

On Borax, I gave you my lammps script that builds a triple junction. Doesnt matter what it is , but I run multiple simulations in that script. Take a look.

Tim

[yanxon]

@timfrol

You're right. Due to the atom escaping, the Arrhenius plot makes no sense. I will construct the PdH system by hand. If you know happen to know a shortcut to do this, I will really appreciate it.

[timfrol]

@yanxon

Howard, I dont think there is a shortcut with lammp, but you can write a simple script, it would be good to have it for future simulations of GBs.

[yanxon]

@timfrol

I generated the PdH_x (x=0.4) with Python. I'm running the same input script with PdH_x, and I obtain this error:

ERROR on proc 0: Neighbor list overflow, boost neigh_modify one (../npair_half_bin_atomonly_newton.cpp:114)
Last command: run ${runtime}

Do you have any idea what's wrong? Btw, I give you the files in Borax.

I tried several method suggested from LAMMPS email such as using neigh_modify etc. Nothing seems to work.

[timfrol]

@yanxon Neighbor list overflow probably means you have too many atoms in one spot. Try building PdH containing less atoms and check that it works.

[yanxon]

@timfrol

It's running now. The dump files coordinates are in fractional coordinates. But, read_data command read the dump file in the format of cartesian coordinates.

I have to convert the fractional coordinates to cartesian coordinates. Now, it's working.

[yanxon]

@timfrol

I have a question. When you run npt simulation and the pressure fluctuates from -3000 to 3000 bars. Is this normal?

[timfrol]

@yanxon

Howard, first of all the fluctuations depend on the number of atoms. The more atoms you have smaller fluctuations. What is your barostat constant? there is a discussion on lammps website regarding what it should be.

[yanxon]

@timfrol

I use 1 as the Pdamp. I tried 10, 100, 1000. It gives the same result. I have 2867 Pd and H atoms.

[timfrol]

@yanxon

what if you try 0.1? Also this constant becomes less important for large systems. Try using about 500 to 1000 atoms per core. On borax you have 36 cores, so you can afford bigger systems, this way you get better statistics as well.

[yanxon]

@timfrol

From what I read in the lammps documentation, the small Pdamp means more fluctuation. That's why I increase it. But, I'm trying 0.1, and it's the same.

[timfrol]

@yanxon

I think it is ok to keep it at 1. Good that you get the same

[yanxon]

@timfrol

I generated the Arrhenius plot:

I felt like this is wrong. Here is the 300K MD simulation:

400K:

500K:

600K:

[timfrol]

@yanxon Excellent work Howard! The Arrhenius plot looks good. What you have found is sometimes called the ballistic region. This is when an atom moves in a straight path without hitting anything. This is usually a shorter displacement which is different from a random walk, as such it is not expected to follow the Einstein relation. This ballistic region becomes less noticable at high T when the random walk is substantial.

To make your low T D coefficients more accurate, people usually skip the ballistic region and only fit the large MSQDs.

Very good work, I am glad you carefully checking your data.

I think we will be running the GB simulations in no time!

[timfrol]

@yanxon

Howard, I am having second thoughts about the ballistic range. The displacement is quite large. Is it possible that your system did not equilibrate at lower temperatures?

How long is your equilibration stage?

Tim

[yanxon]

@timfrol

The temperature equilibrate just fine. The equilibration stage is 0.1 ns. The pressure fluctuates a lot.

[timfrol]

@yanxon

it is not the temperature or pressure only, the structure needs to equilibrate as well. In the past I ran 1 or 2 ns of equilibration. Lets us test the longer equilibration time at low T and see if this funny slope goes away.

[yanxon]

@timfrol

I will try that.