BOP potential for carbon

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS developers and users,

I am encountering some problems while using the BOP potential for carbon (ref: : X.W. Zhou, D.K. Ward, and M.E. Foster (2015), “An analytical bond-order potential for carbon”, Journal of Computational Chemistry, 36). I will try to be detailed so I apologize from now for the long email.

I am using Molecular Static simulations to sample different grain boundary structures in diamond cubic carbon. My simulations are quite “standard”:

  1. I define a bicrystal geometry with an input data file

  2. I rigidly translate one crystal with respect to the other and then perform energy minimization to relax the structure. Multiple translations are considered so the final output is a series of grain boundary structures

  3. I am interested only in the minimum energy structure so only one grain boundary structure is printed at the end of the simulation. However, I save the energy for all the configurations just to check that everything is ok.

I attach an example script (example.inp) to give you an idea of the simulations I am performing.

I am simulating diamond cubic carbon. At the beginning, I used the Tersoff-type potential developed by Erhart and Albe (ref:P. Erhart, and K. Albe (2005), “Analytical potential for atomistic simulations of silicon, carbon, and silicon carbide”, Physical Review B, 71 ). Everything worked fine.

I started using the BOP potential, same script, same simulation “protocol” (clearly, different initial configurations, the equilibrium lattice constant is different). Everything worked well except that I obtained some relaxed configurations showing atoms with a really-really-low/high per-atom energy (i.e. ~ 0.0 eV or ~ -3500 eV) i.e. nonphysical results (see attached an example “conf.dump”).

To understand where is the problem, I took a step back and checked the interatomic forces and per-atom energy of the corresponding un-relaxed configuration (i.e. the configuration obtained after translation before energy minimization). In the following, I will refer to this configuration as configuration A. I also consider a configuration which differs by a really small translation from the previous one (< 0.05 A) but which DOES NOT lead to non-physical per-atom energy values after minimization. In the following, I will refer to this configuration as configuration B.

In both the configurations, there are atoms which are really close and so developed high interatomic forces and energies (NB. I was aware that for some translations local atomic configurations with high energy would be created. However, this had no impact when using the Tersoff-type potential.). These atoms are the same for configuration A and B, because the two configurations differ by a rather small in-plane translation. I checked the interatomic forces and energy for these atoms in both the configurations. Results are reported in the attached pdf “results.pdf”.

For configuration A (the one leading to nonphysical results after energy minimization), interatomic forces are totally out-of-range and I suppose this is the reason of the nonphysical results after energy minimization. However, energy values are contained in a small range (strange).

For configuration B (which DOES NOT lead to nonphysical results after energy minimization, although it is almost identical to configuration A), interatomic forces have more contained values but the per-atom energy value is really high.

I repeated these calculations on the two configurations by using the Tersoff-type potential (I just scaled the interatomic distances in order to be coherent with the different equilibrium lattice constant). For the two configurations A and B, I checked the interatomic forces and the per-atom energy values for the same couple of atom. Results are reported in the attached pdf. With the Tersoff potential, in both cases we have high interatomic forces and per-atom energy values (which is ok, these atoms are really close to each other) but there is no “discontinuity” of these values between configuration A and B.

I have no idea if the results with the BOP potential are linked to its analytical form or to its numerical implementation. I would like to know your opinion.

Best regards,

Carolina Baruffi

Attached:

  • example.in: example LAMMPS script just to give you the idea of the context

  • conf.dump: example of a configuration with non-physical per atom energy values after minimization

  • files used for the test with the Tersoff potential (S79_initial.data + tersoff potential file + minimizing.inp)

  • files used for the test with the BOP potential (S79_initial_BOP.data + BOP potential file + minimizing_BOP.inp)

  • results.pdf: pdf with the results of the tests

PS. unfortunately I am not able to send you the material because my systems are quite big so, also if I compress the files, I am beyond 900 K of attachement. I send you part of the attachement, please tell me if you need more material. Best

results.pdf (340 KB)

example.inp (3.55 KB)

2

thanks for the detailed and specific description of your issue and the steps you did to confirm it.
unfortunately, there is little that i can do myself to investigate this further. my suspicion is, that this may be some issue with the bop pair style implementation.
please keep in mind, that the per-atom energy is usually only a diagnostic property and not used for minimization or MD, so chances are larger, that things have not been checked as thoroughly during implementation.

it would be extremely helpful, if you could construct a tiny test system that exhibits the same anomaly. and then you should check, if that anomaly is consistent when you use a different number of MPI ranks, i.e. run in parallel with different numbers of CPUs. if the per-atom energies change, then this would be a hint what too look for.

axel.

3

I am CCing the folks who implemented BOP for LAMMPS.

Steve

4

Thanks for the heads up. This is very strange indeed. As these types of tests are exactly some of the checks we had to ensure the implementation of the potential in lammps was correct. It would be extremely helpful if we could get a small system as axel described. Unfortunately, my time is very limited right now, but I will see if I can figure out the problem.
Don

5

Rather than plunge into something substantial, I would like to first make sure that we are using the correct potential version and the lammps command. Then can Carolina do a simple calculation (eg, fcc carbon) to make sure that you can reproduce the values in the paper? That is a good way to test if the simulation set up is correct.
Xiaowang

6

Goodmorning,
thank you all for your answers.

ABOUT THE MPI
My simulations have been running on clusters. The first thing I tested was running the minimization for the configuration leading to non-physical values of the per-atom energy on my computer, using just one processor. The problem persist (see results_2.pdf attached, enriched pdf containing results).

NB. to run the simulations on the clusters I modified the cut-off radius of COPPER because in my simulations there is NO copper. The cut-off radius of copper (3.736 A) was forcing me to use a ghost-atom cut-off radius > 3x3.736=11.21 A. By setting it lower (1.736 A rather than 3.736 A) than the one of diamond (2.37 A) I was able to use a ghost-cutoff radius of 7.20 A and speed-up my simulations. Just to be sure, I REPEAT the minimization on my computer with the original values of cut-off, no change in the results (see results_2.pdf)

In any case, all the other test have been performed by using the original cutoff radius.

ABOUT THE POTENTIAL FILE

I downloaded the potential file from NIST:
https://www.ctcms.nist.gov/potentials/
I checked some basic properties of diamond (equilibrium lattice constant, energy, etc) and the values are comparable to the one reported in the paper. I repeat the calculations, results in the pdf results_2.pdf (attached).

The pair_style command I am using is: pair_style bop

ATTACHMENT
Unfortunately, I am dealing with low angle csl so all the configuration files are quite large to be attached. However, I can send you the files to make the test, it takes 10 min in total to run the script and check the force/per-atom energy range in Ovito or whatever visualization program you are using (although I know that, when we are very busy, also 10 min are important…).

Attachements:

  • a python script. Just run it, you will create the initial configurations for the BOP and Tersoff tests
  • the potential files and the script to run the tests with the BOP and with the Tersoff-type potential. To generate the two configurations A and B, you just need to coment/uncomment 2 lines in the LAMMPS script (I put some comments, just search for: # LINES TO COMMENT-UNCOMMENT)
  • an ovito file to visualize configuration A and B from the BOP test.
  • the “enriched” results_2.pdf file

I tried to sent everything in a zip file but the mail is rejected.

Let me know if you need some more info.

Best,
Carolina

visualization.ovito (25.6 KB)

SiC_Erhart-Albe.tersoff (1.73 KB)

test_BOP.inp (2.09 KB)

test_EA.inp (2.06 KB)

results_2.pdf (360 KB)

ini_conf.py (1.54 KB)

CCu_v2.bop.table (476 KB)