How can I do this - Create an FCC Copper System with hydrogen atoms placed at random positions to observe diffusion

LAMMPS LAMMPS Mailing List Mirror
1

I am trying to generate a fcc copper box of dimensions 75 x 50 x 30 angstroms. In this system I want to manually add a few hydrogen atoms at random positions. Then I would like to thermalize the system at a particular temperature, say 450 K. Eventually, after the system has thermalized to that particular temperature, I would like to calculate the MSD and thereby obtain the diffusion coefficient. Until now what I have done are the followings:

1) Generated a copper box using the following script:

units metal

dimension 3

boundary p p p

atom_style atomic

atom_modify map array

lattice fcc 3.61

region box block 0 75 0 50 0 30 units box

create_box 1 box

lattice fcc 3.61 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

create_atoms 1 box

mass 1 63.5460

write_data Data.R1

print “Done”

2) In the generated Data.R1 file, the first few lines were:

9996 atoms

2 atom types

0.0000000000000000e+00 7.5000000000000000e+01 xlo xhi

0.0000000000000000e+00 5.0000000000000000e+01 ylo yhi

0.0000000000000000e+00 3.0000000000000000e+01 zlo zhi

Masses

1 63.5460

Atoms # atomic

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0

2 1 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

So, in order to incorporate Hydrogen atoms for a test run, I added a new mass, the mass of hydrogen in the following manner:

Masses

1 63.5460

2 1.00794

and also, to place a hydrogen atom at the second atoms position, I replaced the atom type with 2 instead of 1 in the second row second column of the following lines

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0

2 2 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

My first question is if this is a useable/correct method of adding hydrogen atoms in this system? If not, in what ways can I add hydrogen atoms in the copper system?

Next what I did was thermalized the system (I think) at 450 K using the following code snippet:

units metal

dimension 3

boundary p p p

atom_style atomic

atom_modify map array sort 0 0

-------------Create atoms and box-------------------------

read_data Data.R1

region e1 block INF INF INF INF INF INF

group e1 dynamic all region e1

-------------Define Interatomic Potential-------------------------

pair_style bop

pair_coeff * * CuH.bop.table Cu H

comm_modify cutoff 11.4

thermo_modify lost ignore

velocity all create 450 25674 dist gaussian

fix 1 all npt temp 450 450 0.1 iso 0.0 0.0 10

fix 2 all momentum 1 linear 1 1 1

minimize 1.0e-8 1.0e-8 1000000 1000000

timestep 0.001

reset_timestep 0

compute 1 e1 temp

thermo 20000

dump 2 all atom 20000 copper_box_300K2.lammpstrj

dump 3 all custom 100000 Copper2_box_300k_final.dump id type x y z vx vy vz

thermo_modify lost ignore

restart 1000000 poly.restart

run 10000000

unfix 1

write_data Data25.R1

write_restart Copper2_box_300k.restart

print “Done”

IS THIS CORRECT?

2

I am trying to generate a fcc copper box of dimensions 75 x 50 x 30 angstroms. In this system I want to manually add a few hydrogen atoms at random positions. Then I would like to thermalize the system at a particular temperature, say 450 K. Eventually, after the system has thermalized to that particular temperature, I would like to calculate the MSD and thereby obtain the diffusion coefficient. Until now what I have done are the followings:

1) Generated a copper box using the following script:

units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
lattice fcc 3.61
region box block 0 75 0 50 0 30 units box
create_box 1 box
lattice fcc 3.61 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
mass 1 63.5460

write_data Data.R1
print “Done”

2) In the generated Data.R1 file, the first few lines were:

9996 atoms
2 atom types

0.0000000000000000e+00 7.5000000000000000e+01 xlo xhi
0.0000000000000000e+00 5.0000000000000000e+01 ylo yhi
0.0000000000000000e+00 3.0000000000000000e+01 zlo zhi

Masses

1 63.5460

Atoms # atomic

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
2 1 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

So, in order to incorporate Hydrogen atoms for a test run, I added a new mass, the mass of hydrogen in the following manner:

Masses

1 63.5460
2 1.00794

and also, to place a hydrogen atom at the second atoms position, I replaced the atom type with 2 instead of 1 in the second row second column of the following lines

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
2 2 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

My first question is if this is a useable/correct method of adding hydrogen atoms in this system?

Yes.
(Note: If you are using atom_style full, then you would modify the
number in the 3rd column, not the 2nd column. But this format is
correct for atom_style atomic)

If not, in what ways can I add hydrogen atoms in the copper system?

What you are doing should work. And I think it's a good idea that you
are learning how to create and edit the DATA file yourself. If at
some later point you want to use external tools to help you build it,
here are a list of them:
https://lammps.sandia.gov/prepost.html
…and some examples for creating crystalline materials:
http://moltemplate.org/doc/moltemplate_manual.pdf#page=38
https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---part-2
(These tools are usually only necessary if you have molecules and
bonds somewhere in your system.)

Next what I did was thermalized the system (I think) at 450 K using the following code snippet:

units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array sort 0 0

# -------------Create atoms and box-------------------------

read_data Data.R1
region e1 block INF INF INF INF INF INF
group e1 dynamic all region e1

# -------------Define Interatomic Potential-------------------------

pair_style bop
pair_coeff * * CuH.bop.table Cu H
comm_modify cutoff 11.4
thermo_modify lost ignore
velocity all create 450 25674 dist gaussian

fix 1 all npt temp 450 450 0.1 iso 0.0 0.0 10
fix 2 all momentum 1 linear 1 1 1

minimize 1.0e-8 1.0e-8 1000000 1000000
timestep 0.001
reset_timestep 0
compute 1 e1 temp

thermo 20000
dump 2 all atom 20000 copper_box_300K2.lammpstrj
dump 3 all custom 100000 Copper2_box_300k_final.dump id type x y z vx vy vz
thermo_modify lost ignore
restart 1000000 poly.restart
run 10000000

unfix 1

write_data Data25.R1
write_restart Copper2_box_300k.restart

print “Done”

IS THIS CORRECT?

The best way to do this is to run the simulation yourself.
Once you eventually are able to run the simulation without crashing,
then it is time to visualize the system.

In LAMMPS use the "dump" command and use VMD or OVITO to read the file
and display the trajectory.
For your system, I recommend OVITO. If your system has a lot of
molecules with complicated bond topology, I recommend VMD.

For VMD, at least, I posted some instructions here:

https://github.com/jewettaij/moltemplate/blob/master/examples/all_atom/force_field_explicit_parameters/aluminum_crystal_strain/README_visualize.txt

Good luck.

Andrew

3

Thanks for your response.

I am little confused about the first answer though - “(Note: If you are using atom_style full, then you would modify the number in the 3rd column, not the 2nd column. But this format is
correct for atom_style atomic)”

The 3rd column in the following lines is the bold one:

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
2 2 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

as in the x - coordinate. Why would I be changing the number on that column to change the atom type?
Also, in these 2 lines of index, atom type and coordinates, what do the last 3 columns represent? (0 0 0)

Also, I did run the thermalization code snippet and it ran fine without crashing. So if I were to calculate MSD using the newly generated data file after a 10 ns thermalization simulation as mentioned above, would the following code snippet suffice?

clear

log copper_box_300k.log

units metal

dimension 3

boundary p p p

atom_style atomic

atom_modify map array sort 0 0

read_data Data25.R1

pair_style bop

pair_coeff * * CuH.bop.table Cu H

comm_modify cutoff 11.4

mass 1 63.5460

mass 2 1.00794

timestep 0.001

thermo 200000

fix 1 all nvt temp 600 600 0.1

fix 3 all momentum 1 linear 1 1 1

thermo_style custom step temp press vol enthalpy ke pe etotal lx ly lz density atoms

group Cu type 1

group H type 2

compute 2 H msd

fix 4 H ave/time 20000 5 100000 c_2 file H.txt mode vector

thermo_style custom step temp press vol enthalpy ke pe etotal lx ly lz density atoms

thermo_modify lost warn

dump 2 all atom 200000 copper_box_300K.lammpstrj

dump 3 all custom 200000 copper_box_300k_final.dump id type x y z vx vy vz

restart 100000 poly.restart

write_data Data.R2

write_restart copper_box_300k.restart

run 1000000 upto

4

I ran the following code snippet in order to calculate MSD of two hydrogen atoms that I placed manually in an FCC copper system.

units metal

dimension 3

boundary p p p
atom_style atomic

atom_modify map array sort 0 0
read_data Data25.R1

pair_style bop

pair_coeff * * CuH.bop.table Cu H

comm_modify cutoff 11.4

mass 1 63.5460

mass 2 1.00794

timestep 0.001

thermo 200000

fix 1 all nvt temp 600 600 0.1

fix 3 all momentum 1 linear 1 1 1
thermo_style custom step temp press vol enthalpy ke pe etotal lx ly lz density atoms
group Cu type 1

group H type 2

compute 2 H msd
fix 4 H ave/time 20000 5 100000 c_2 file H.txt mode vector

thermo_style custom step temp press vol enthalpy ke pe etotal lx ly lz density atoms

thermo_modify lost warn

dump 2 all atom 200000 copper_box_300K.lammpstrj

dump 3 all custom 200000 copper_box_300k_final.dump id type x y z vx vy vz

restart 100000 poly.restart

write_data Data.R2

write_restart copper_box_300k.restart

run 1000000 upto

The output H.txt file shows me the following data:

Time-averaged data for fix 4

TimeStep Number-of-rows

Row c_2

100000 4

1 1.48716

2 0.850136

3 0.255434

4 2.59273

200000 4

1 2.71356

2 0.428951

3 1.96094

4 5.10345

300000 4

1 13.7711

2 0.68944

3 4.2548

4 18.7153

400000 4

1 14.5508

2 0.413935

3 4.48876

4 19.4535

500000 4

1 12.8003

2 0.643756

3 5.38209

4 18.8262

600000 4

1 14.661

2 0.449091

3 5.24004

4 20.3502

700000 4

1 14.1748

2 0.434456

3 5.04747

4 19.6567

800000 4

1 15.5679

2 1.06356

3 5.82025

4 22.4517

900000 4

1 20.0385

2 7.14621

3 2.75025

4 29.935

1000000 4

1 7.53698

2 3.22376

3 3.24066

4 14.0014

Could you please help me understand the data? This data was produced based on the fix 4 H ave/time command. All I need to find is the MSD of the 2 hydrogen atoms in the system at every or multiple timesteps upto 1ns. Is there a simpler command that I can use to get the MSDs at different timesteps upto 1ns?

5

Thanks for your response.

I am little confused about the first answer though - "(Note: If you are using atom_style full, then you would modify the number in the 3rd column, not the 2nd column. But this format is
correct for atom_style atomic)"

I'm sorry I wrote that. I did not mean to confuse you. (I also
wanted help other users who use "atom_style full".) You are using
"atom_style atomic", so my comment was not relevant to you.

The 3rd column in the following lines is the bold one:

1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0 0 0
2 2 1.8049999999999999e+00 1.8049999999999999e+00 0.0000000000000000e+00 0 0 0

as in the x - coordinate. Why would I be changing the number on that column to change the atom type?

It doesn't. In your case, you would modify the 2nd column. (You did
it correctly in your first email.)

Also, in these 2 lines of index, atom type and coordinates, what do the last 3 columns represent? (0 0 0)

Those are "image flags". Most simulations in LAMMPS use periodic
boundary conditions. Suppose that at the beginning of your simulation
a particle is located in the middle of the box. It has indices 0,0,0.
Later on in the simulation, it might drift forward in the +x direction
out of the box and into the next (periodic) box. Now the particle has
indices 1,0,0 and the "x" position will be "wrapped" back to the
interval xlo...xhi (approximately; in your case xlo=0 and
xhi=7.5000e+01).

Also, I did run the thermalization code snippet and it ran fine without crashing. So if I were to calculate MSD using the newly generated data file after a 10 ns thermalization simulation as mentioned above, would the following code snippet suffice?

I don't know. I don't use "fix msd". (or "fix ave/time". For me it
is easier to write a short script with a for-loop which can parse the
trajectory file created by the "dump" command, and calculate things I
need, like MSD.) I'm sorry I cannot help you with this.

cheers

Andrew
P.S. I mostly simulate molecules. When calculating the MSD of
molecules (due to internal changes in the molecule's shape, not
rotation or drifting), I've been using this code:
https://github.com/jewettaij/superpose3d (in python)
https://github.com/jewettaij/superpose3d_cpp (in c++)
But since your simulations do not contain molecules, I don't think
this that code is relevant to you.
(...So please ignore this. Perhaps I should not have mentioned it.)