Compression in NPT ensemble with ReaxFF potential

LAMMPS LAMMPS General Discussion
, ,
1

Hello everyone
As I am new to LAMMPS and molecular dynamics, I am trying to regenerate some data reported in literature with LAMMPS. Specifically I am working with the paper titled ‘An Atomistic Carbide-Derived Carbon Model Generated Using ReaxFF-Based Quenched Molecular Dynamics’ (doi:
C | Free Full-Text | An Atomistic Carbide-Derived Carbon Model Generated Using ReaxFF-Based Quenched Molecular Dynamics).
The first step reported in the paper involves quenching a structure in NVT ensemble. I was able to do that and produced results in close agreement with the paper. The LAMMPS input script is as follows:

units real # for reaxff forcefiled
atom_style charge # for reaxff forcefiled
boundary p p p # 3d periodic boundary

region box_ block 0 74.88 0 74.88 0 74.88 # box dimension 7.488 nm
create_box 1 box_
region whole block INF INF INF INF INF INF units box
create_atoms 1 random 20000 539816 NULL # 20000 atoms are placed at random inside the box
group 1 region whole

neighbor 2 bin
neigh_modify every 1 delay 0 check yes

mass 1 12.0107 # mass of carbon atom
pair_style reax/c NULL # interatomic potential
pair_coeff * * ffield.reax-2013.C C # coefficient for carbon atom

fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c # for reaxff

thermo 10
minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

velocity all create 3500 729014 mom yes rot yes dist gaussian
timestep 0.5 # fs
dump 1 all custom 1000 C.* id xs ys zs fx fy fz
thermo 100
thermo_style custom step time temp density press pe ke etotal

restart 10000 file.*.restart
restart 500 file1.restart file2.restart

fix 2 all nvt temp 3500 3000 10 # tdmap = 10 fs (unit = real = fs)

run 100000 # 50 ps

This same input script is run for three different quench steps. The density, cell sizes and ring size distribution for all simulation matches to that of reference paper.

1kps
1kps717×578 12.3 KB

10kps
10kps718×573 7.72 KB

100kps
100kps740×569 7.26 KB

The next step is compressing the previously quenched structure in NPT ensemble. As stated in paper, 'Isothermal temperature control and isobaric pressure control at 3000 K and 20,000 atm were used with damping constants of 10 fs and 100 fs,respectively. I cannot reproduce the data reported in this step.
In the paper it is reported that the final cell size and density for QMD-10c structure will be 6.489 nm and 1.467 g/cc after 250 ps of simulation time. In my simulation, after 52.5 ps the cell size and density appear to be 5.804 nm and 2.0398 g/cc.
For simulating the NPT step, the restart file generated previously was read in and following command line was used:

fix 2 all npt temp 3000 3000 10 iso 20000 20000 100 # tdmap = 10 fs (unit = real = fs) pdamp = 100fs

I cannot figure out how is it possible that the results from NVT ensemble matched but not the results from NPT ensemble? Any help or suggestion is much appreciated.

Regards
Koushik

1 Like
2

Hi Koushik,
Have you tried to use aniso instead of iso for the pressure control in fix npt? I am not sure what your system exactly is, but iso is usually only appropriate for fluids.
Best,
Simon

2 Likes
3

Sir
The initial structure is fluid-like and structure after quenching is solid. As stated in paper “…major structural changes do not occur after the system has condensed from a fluid to solid state, so additional quenching below 3000 K was not considered.”
I will try to run the simulation with aniso. Many thanks for your suggetion.
Sincerely
Koushik

4

Sir
I tried using aniso for pressure control. But there’s no improvement in result. After 12 ps, the density becomes 1.75 g/cc which is much larger than the final density 1.467 g/cc. Moreover the cell dimensions are also smaller than the desired values. The final values in paper are taken after 250 ps of simulation. Here is the input script

units real # for reaxff forcefiled
atom_style charge # for reaxff forcefiled
boundary p p p # 3d periodic boundary

read_restart file.100000.restart

neighbor 2 bin
neigh_modify every 1 delay 0 check yes

mass 1 12.0107 # mass of carbon atom
pair_style reax/c NULL # interatomic potential
pair_coeff * * ffield.reax-2013.C C # coefficient for carbon atom

fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c # for reaxff

timestep 0.5 # fs
dump 1 all custom 1000 C.* id xs ys zs fx fy fz
thermo 100
thermo_style custom step time temp press density pe ke etotal lx ly lz

fix 2 all npt temp 3000 3000 10 aniso 20000 20000 100 # tdmap = 10 fs (unit = real = fs) pdamp = 100fs

run 500000 # 250 ps

I am wondering if my units, atom_style, neighbor settings and other settings are appropriate for this simulation. I would be grateful if you could let me know your thoughts.

Sincerely
Koushik Sarkar

5

Sir
The parameters at the beginning of the NPT simulation as taken from the log file is as follows:

Step Time Temp Press Density PotEng KinEng TotEng Lx Ly Lz
100000 0 3000.0056 -7139.7631 0.95005745 -3034307.1 178839.99 -2855467.1 74.88 74.88 74.88

The pressure of the structure generated previously from NVT quenching appears to be -7139.7631 atm. My question is it possible to have a negative pressure? Is it possible to set the pressure to zero without altering the structure before carrying out the actual NPT compression?

6

Yes. A negative pressure indicates, that the system wants to shrink.

The (instantaneous) pressure of the system is a consequence of the geometry and the force field. Thus you cannot set it.

1 Like
7

Why performing anisotropic NPT from an already solidified structure? The whole point of using aniso NPT is to allow the solid structure to adopt the shape it wants…
Simon

1 Like
8

According to the reference paper, the NPT compression is carried out to change the pore size distribution of the solid structure so as to conform to the experimental pore size distribution. As stated in the paper, “We also investigate the role of a novel post-quenching compression step as a means to adjust pore size distribution.”

9

Before wasting too much time on fine tuning the protocol details, I would suggest running through every step up to the NPT run, for a second set of random initial positions (using a different seed).

If you get two very different results up to that point then you already know you will get different results from the NPT. (You already have some deviation from the paper results for the one set of initial conditions you’ve tried.)

I had a quick glance through the paper you linked and they do not mention doing statistically independent replicates. They do not appear to have error bars or multiple-trajectory result comparisons. So, I don’t know if their protocol is repeatable. I can predict that it’s not very repeatable from theory*, but I certainly can’t prove it one way or another. You should check if the protocol even gives the same results twice before worrying that it doesn’t give the same results as the paper.

*Any structure as amorphous as the one they produced would be very heterogeneous at the nanoscale. Why would one nanoscale block have a pressure response that is at all similar to another very different nanoscale block?

3 Likes
10

I have carried out the NVT ensemble simulation for quench rate of 100 K/ps and 10 K/ps with three different sets of random initial position and velocity. The ring size distributions obtained from these simulation appear to have close agreement with each other and to the data reported in the paper. The figures are as follows:
image
image

I have not carried out the simulations for 1 K/ps as it takes a long time to complete.
I think these results indicates the data can be repeated. Please let me know your thoughts on this issue.

Also I would like to point out that I have found the github profile of the corresponding author and input script for the NVT ensemble quenching was given there. The data regenerated here is obtained from running the input script taken from github on my machine with different seeds. But I could not find the input script for NPT compression on his github profile.

11

Now you can check if the NPT part is repeatable!

1 Like
12

I have carried out the subsequent NPT compression on the quenched structure generated with different random number generator seeds. The compression step is supposed to run for 250 ps but for conserving time I ended each simulation at 50 ps. I have worked with three different quenched structures with anisotropic pressure control. I have also simulated the first quenched structure in isotropic pressure control. The results are summarized in the following table:

image
image1134×316 13.3 KB

From the above data, I think the results appears to be repeatable.

13

Now you have three options:

  1. Use your current NPT final structures, and note the discrepancy when you write your final paper.
  2. If you are confident that the final density from the paper is correct (from experiments), you can simply specify the density as your simulation target. First use fix deform to change your box to its target size, and then equilibrate it at fixed size, both using the NVT integrator over a timescale equivalent to the NPT equilibration duration you used.
  3. Contact the original paper authors (since they have a GitHub repo, you can probably easily reach them with a GitHub issue) and obtain their script and LAMMPS version to compare with yours.

I cannot tell you which of the three options is the best, without knowing a lot more about the field of interest, your remaining computational time available, your level of confidence in experimental or other “ground truths” – that is, these are questions you need to answer either yourself as a researcher, or with advice from a supervisor or local expert who can give you in-depth, face-to-face advice.

I am not being dismissive – you have done important work understanding your system, and you should be proud of yourself (I’ve learned something too!). But that work has led you somewhere where you really need advice from experts, not strangers on the Internet :slight_smile:

1 Like
14

Thanks for your detailed response.
I have been trying to reach the author for over a month but have been unable to reach him via email. Recently I was able to reach him via twitter. He gave me a new mail address and said he would get back to me as soon as possible.
I am simply trying to figure out the right input script to compress a solid structure in NPT ensemble. I wish to carry out similar procedure on other structures and possibly with other potential functions.
Thanks again for your encouraging words. Unfortunately I have no other resources to take help from at this moment. Even if it does not solve my problem altogether, I get to explore new ideas by asking questions on the internet. I hope for your continued support.

15

I was able to get the input script for NVT quenching from the github profile of the author. I have noticed that in the input script a fix is applied as follows:

fix 6 all reax/c/bonds 1000 bonds.quench

This fix produces the bonds.quench file. I could not find reax/c/bonds fix in LAMMPS documentation but I found reaxff/bonds. According to the documentation, ‘Write out the bond information computed by the ReaxFF potential specified by pair_style reaxff…’ . I am wodering if this bond information is necessary for the subsequent NPT compression.
Is it possible to define a LAMMPS readable input structure file with bonds between atoms using the file generated from reax/c/bonds fix and the dump file?
Should a structure with bonds between atoms behave significantly different from a structure without bonds, if simulated with same pair_style potential?

16

The reax/c/bonds style was renamed to reaxff/bonds

1 Like
17

No. It is for analysis only.

Yes, but adding explicit bonds to a ReaxFF simulation would be a big mistake.

1 Like