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”:
I define a bicrystal geometry with an input data file
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
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.
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)