Pair_style hybrid with EAM/alloy and comb3

Hi all. I am trying to simulate a Cu-Zr-O metallic glass using an EAM potential for interaction between Cu and Zr atoms and the comb3 potential for interactions with O. Please refer to the input script below. I am unsure if I am implementing the comb3 potential correctly since I haven’t used it before.

units metal # Angstrom, eV, ps, bar
boundary p p p # Periodic boundaries in all dimensions.
atom_style charge # Pair potentials.
newton on

#------------------------------------------------------------------

Input data

#------------------------------------------------------------------

read_data Zr50Cu50_50K_small

set type 1 type/fraction 3 0.02 4647
set type 2 type/fraction 3 0.02 5746

mass 1 91.2
mass 2 63.5
mass 3 16.0

#------------------------------------------------------------------

Define interatomic potential.

#------------------------------------------------------------------
pair_style hybrid comb3 polar_off eam/alloy

pair_coeff * * comb3 ffield.comb3 Zr Cu O
pair_coeff * * eam/alloy ZrCu.lammps.eam Zr Cu NULL

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

#------------------------------------------------------------------

Set up thermo.

#------------------------------------------------------------------
thermo 1000
thermo_style custom step temp lx ly lz press pxx pyy pzz pe etotal enthalpy vol density
thermo_modify norm yes flush yes

#------------------------------------------------------------------

Set up charge equilibration.

#------------------------------------------------------------------
fix 1 all qeq/comb 1 0.0001

#------------------------------------------------------------------

Minimize input structure.

#------------------------------------------------------------------
minimize 1.0e-4 1.0e-6 1000 10000

#------------------------------------------------------------------

Initialize velocity and equilibrate in NVE.

#------------------------------------------------------------------
fix 2 all nve
timestep 0.0005
velocity all create 50 4928459 mom yes rot yes dist gaussian
run 100000
unfix 2

Every time I run this script, I get non-numeric pressure and energies and lammps says there is a segmentation fault. Please refer to the output below.

LAMMPS (29 Oct 2020)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (1.0918818 1.0918818 1.0918818) to (53.308118 53.308118 53.308118)
4 by 4 by 5 MPI processor grid
reading atoms …
8192 atoms
read_data CPU = 0.015 seconds
Setting atom values …
85 settings made for type/fraction
Setting atom values …
92 settings made for type/fraction
PairComb3: polarization is off
Reading COMB3 potential file ffield.comb3 with DATE: 2015-09-18
element[1] = Zr, z = 1.83153
element[2] = Cu, z = 1.39726
element[3] = O , z = 1.37179
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 13
ghost atom cutoff = 13
binsize = 6.5, bins = 9 9 9
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair comb3, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(2) pair eam/alloy, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(3) neighbor class addition, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Setting up cg style minimization …
Unit style : metal
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 8.817 | 8.820 | 8.822 Mbytes
Step Temp Lx Ly Lz Press Pxx Pyy Pzz PotEng TotEng Enthalpy Volume Density
0 0 52.216236 52.216236 52.216236 -45557.926 -47746.027 -44066.905 -44860.844 -4.9107505 -4.9107505 -5.404925 142369.41 7.2503239
19 0 52.216236 52.216236 52.216236 nan nan nan nan nan nan nan 142369.41 7.2503239
Loop time of 29451 on 80 procs for 19 steps with 8192 atoms

99.9% CPU use with 80 MPI tasks x 1 OpenMP threads

Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
-4.91075047656185 -5.39175998754189 nan
Force two-norm initial, final = 187.24871 nan
Force max component initial, final = 15.388346 0.19915337
Final line search alpha, max atom move = 2.4977204e-12 4.9742944e-13
Iterations, force evaluations = 19 43

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

Pair | 53.381 | 64.55 | 76.818 | 81.0 | 0.22
Neigh | 0.032452 | 0.033827 | 0.035448 | 0.5 | 0.00
Comm | 77.592 | 95.593 | 108.63 | 68.1 | 0.32
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 29276 | 29291 | 29296 | 2.3 | 99.46
Other | | 0.004407 | | | 0.00

Nlocal: 102.400 ave 113 max 92 min
Histogram: 3 2 7 11 16 17 15 6 2 1
Nghost: 3093.35 ave 3121 max 3068 min
Histogram: 4 9 9 8 12 10 14 6 4 4
Neighs: 7268.34 ave 8102 max 6286 min
Histogram: 5 3 5 9 7 14 16 14 5 2
FullNghs: 2319.60 ave 5309 max 730 min
Histogram: 5 21 12 17 10 5 4 1 2 3

Total # of neighbors = 581467
Ave neighs/atom = 70.979858
Neighbor list builds = 2
Dangerous builds = 0
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 13
ghost atom cutoff = 13
binsize = 6.5, bins = 9 9 9
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair comb3, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(2) pair eam/alloy, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(3) neighbor class addition, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Setting up Verlet run …
Unit style : metal
Current step : 19
Time step : 0.0005
Per MPI rank memory allocation (min/avg/max) = 7.692 | 7.695 | 7.697 Mbytes
Step Temp Lx Ly Lz Press Pxx Pyy Pzz PotEng TotEng Enthalpy Volume Density
19 50 52.216236 52.216236 52.216236 nan nan nan nan nan nan nan 142369.41 7.2503239

===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 20 PID 58796 RUNNING AT hpc406
= KILLED BY SIGNAL: 11 (Segmentation fault)

I would really appreciate any insights on what could be going wrong and how I can fix it. Thanks!

What is wrong is the use of pair style hybrid. It is wrong it two ways: a) it doesn’t make sense in terms of the physics. Both of the potentials are manybody potentials and have different ways to include the manybody interactions. You cannot just say that you want one for some pairs and another for some other pairs. That is not how manybody interactions work, and there is no consistent way to enter these. Even if there would be a way, the result would be bogus.
And b) I don’t think that COMB (or COMB3) for your LAMMPS version is compatible with a hybrid potential at all.

The only way to model your system consistently would be to use COMB3 for all interactions.

Thank you Axel.