Pair_style table buck/coul/long

LAMMPS LAMMPS Beginners
1

I am trying to perform molecular dynamics on a borophosphate system with the potentials from the Deng and Du article (https://ceramics.onlinelibrary.wiley.com/doi/10.1111/jace.16082). Using short-range buckingham and long-range coulombic. I got a lot of atom loss and range out errors. So I made the splice correction with a script, and I’m using the pair_style table. However, I continue to get these range out errors and loss of atoms, I wanted to understand what I can do. Below is the input and script.
test.py (1.5 KB)
////////////////////////////////////////////////////////////////////////////////////////////
pair_table_full.txt (323.5 KB)

Definir os tipos de átomos

units metal
dimension 3
boundary p p p
processors * * *

#create box
atom_style charge
region box block 0 34.19 0 34.19 0 34.19
create_box 5 box
create_atoms 1 random 700 12345 NULL overlap 1.7 maxtry 200000
create_atoms 2 random 1960 12345 NULL overlap 1.7 maxtry 200000
create_atoms 3 random 280 12345 NULL overlap 1.7 maxtry 200000
create_atoms 4 random 140 12345 NULL overlap 1.7 maxtry 200000
create_atoms 5 random 140 12345 NULL overlap 1.7 maxtry 200000

mass 1 10.811 # B
mass 2 15.9994 # O
mass 3 30.9738 # P
mass 4 39.0983 # K
mass 5 87.62 # Sr

atom charges

set type 1 charge 1.8000 # B
set type 2 charge -1.2000 # O
set type 3 charge 3.0000 # P
set type 4 charge 0.6000 # K
set type 5 charge 1.2000 # Sr
kspace_style ewald 1.0e-5

inter-atomic potential

pair_style hybrid/overlay coul/long 10 table linear 1000 ewald
pair_coeff * * coul/long
pair_coeff 2 2 table pair_table_full.txt O-O 10
pair_coeff 2 4 table pair_table_full.txt K-O 10
pair_coeff 2 5 table pair_table_full.txt Sr-O 10
pair_coeff 2 3 table pair_table_full.txt P-O 10
pair_coeff 1 2 table pair_table_full.txt B-O 10
neighbor 2.0 bin
neigh_modify every 1 check yes

outputs

thermo 1000
thermo_style custom step temp pe press lx density
fix 1 all nve/limit 0.01

Definindo o integrador

timestep 0.005
run 10000
unfix 1
timestep 0.005
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
minimize 1e-10 1e-10 1000 10000
unfix 1
write_data relaxed1_configuration.data

Relaxamento a 300 K em NVT

reset_timestep 0
fix 2 all nvt temp 300.0 300.0 100.0
run 10000 # Ajuste o número de passos conforme necessário
unfix 2

Gravação da configuração relaxada

write_data relaxed1_configuration.data

timestep 1

outputs

thermo 1000
thermo_style custom step temp pe press lx density

mixing

minimize 1.0e-8 1.0e-8 10000 10000
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
write_data final.dat
neighbor 2.0 bin
neigh_modify every 1 check yes
timestep 1

outputs

thermo 1000
thermo_style custom step temp pe press lx density

mixing

minimize 1.0e-8 1.0e-8 10000 10000
fix 1 all nvt temp 300 300 100
run 10000
unfix 1
write_data final.dat
////////////////////////////////////////////////////////////////////////////////////////////////////////

2

Hello,

Any reason you keep alternating between energy minimization and molecular dynamics steps ?

Simon

3

I was doing this to try to relax the system, as a suggestion to avoid such errors, but it keeps going wrong.

4

If the simulation is already running in molecular dynamics with a timestep of 1, it won’t help to perform an additional energy minimization unless the topology was modified (i.e. added/removed atoms).

Most likely here your simulation keep crashing because some pair_coeffs are badly defined or missing.

5

I understood. But I’ve already tested it in single steps, and I keep getting these errors. In the article, it mentions the Buckingham repulsion, which can be fixed with a splice, however I generated the table through the script, and the system still does not work. The interactions are B-O, P-O, Sr-O, K-O, O-O.

6

Did you define interaction between the other pairs ? I.e. is there anything that keep them from overlapping ?

7

any other interaction I define as “* *” for type 0 1 0

8

Please see this post about what to report (e.g. LAMMPS version and platform you are running on, where exactly the run stops) and how to correctly quote input files (otherwise the forum software will interpret special characters as markup, e.g. comments become headlines) and they can become hard to read.

What kind of a “splice correction” are you talking about? I have not heard of that term before.

Some general comment (pun intended): This is an international forum using the English language as a common ground to communicate. Thus it is counterproductive to use comments in your native(?) language when quoting files or providing them as attachments. You are withholding important information from people that do not know your language. English isn’t my native language either, but I learned very early on in my training as a researcher that it is a big advantage to use English for all my notes, comments, and even variable and function names so it would be easier to share and discuss with collaborators.

There are multiple issues here and worrying about the potential function is the least concern. Of course, if you still feel that there could be a problem then you first have to validate that your tables are correct. You would make your life much simpler by using tabulation tool provided in the LAMMPS distribution: 10. Auxiliary tools — LAMMPS documentation. This allows to determine the derivative numerically and thus validate if your force computation for the table and also the format is correct.
Even more so, I would first tabulate the unmodified Buckingham potential and compare it to the output of pair_write.

Some more comments about your input:

Here you are creating completely randomized positions. This is a very bad idea in general. There is no way that the expected borophosphate structure will form correctly. A big problem is that the Buckingham potential is not sufficiently repulsive at short distances to be used for anything but geometries that are already reasonably close to the final structure. A better approach would be to start from templates use molecule files and for cases where different elements show up at equivalent positions, first generate the geometry with just one element and then replace a fraction of them with the other according to the expected stoichiometry and then try to carefully melt and sinter the structure.

What did the people do in the paper that you are referencing to prepare their initial geometries?

Here you are using a timestep of 5 femtoseconds with is quite aggressive for a system with first row elements like boron and oxygen. More reasonable would be 2 fs and conservative (especially since you are starting from high potential energy with your random geometry) would be 1 fs.

Here you are switching to a 1 picosecond timestep. That is completely bonkers and will never work (see above).

In summary, I believe the major problem is your preparation of the initial geometry and your too large timestep is making matters worse. No tweak of the potential function can correct for that.

9

In a test I carried out, I completed a dynamic with other less energetic potentials, took the final structure and applied these Buckingham potentials, and still got the errors.

10

You are just disregarding all the advises, including the absolutely critical comment by @akohlmey on timestep.

11

When a simulation runs to completion it says nothing about its correctness. It can still be completely bogus or describing a system different from what you intend.

That said, I dare you to show us an input deck that can run a stable time integration of a condensed atomic/ionic system at a 1 picosecond time step.

12

Thank you for all the comments, I will rework my inputs.

13

What about answers to the questions that I had asked you?

The “contract” of a forum like this one between the people asking questions and the people responding is that the people asking “pay” for the responses by satisfying the curiosity of the people providing the responses. TANSTAAFL!