running lammps code for CaCO3

Dear All

I have some problem for running lammps code for CaCO3 (calcite). I find position for each atom in that box and also I used pair potential parameter ibut when I run program I see some errors. I want to find thermal conductivity from green kubo formalism. I have used this code for alloys and I did not have any problem (but I used read data for alloys that I don’t have here for ionic material.

units metal
atom_style full

variable a equal 4.990
variable b equal 4.990
variable c equal 17.060

variable 1_3 equal 1.0/3.0
variable 2_3 equal 1.0/3.0
variable 5_6 equal 5.0/6.0
variable 1_6 equal 1.0/6.0
variable 1_4 equal 1.0/4.0
variable 3_4 equal 3.0/4.0
variable 1_12 equal 1.0/12.0
variable 5_12 equal 5.0/12.0
variable 7_12 equal 7.0/12.0
variable 11_12 equal 11.0/12.0

variable s equal 0.257

lattice custom 1.0 a1 $a 0 0 a2 0 $b 0 a3 0 0 c basis 0 0 0 basis {2_3} {1_3} {1_3}…

dimension 3
boundary p p p

region box block 0 8 0 8 0 8

create_box 3 box

create_atoms 1 box basis 1 1 basis…

Set masses

variable MA equal 40.078 # Ca
variable MB equal 12.0107 # C
variable MC equal 15.9994 # O

mass 1 {MA} mass 2 {MB}
mass 3 ${MC}

Create groups

group Ca type 1 #A for cations
#group CO3 type 2 3 #C for anions
group C type 2 #A for cations
group O type 3 #A for cations

Set charges

set group Ca charge 2

set group C charge +1.343539
set group O charge -1.114513

variable V equal vol

variable dt equal 0.0015

variable T equal 100

variable P equal 0.0

variable kB equal 8.6173423e-5 # [eV/K] Boltzman constant

variable eV2J equal 1.60217646e-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal {eV2J}/{A2m}/${ps2s}

variable NBR equal 1000000 # barostating run
variable NER equal 1000000 # equilibrization run
variable NPR equal 2000000 # production run

variable NBRave equal round({NBR}/2) variable NERave equal round({NER}/2)
variable NPRave equal round(${NPR}/2)

variable NCL equal 1000 # correlation length
variable NSI equal 1 # sample interval
variable NDI equal {NCL}*{NSI} # dump interval

Buckingham Potential

pair_style buck/coul/long 10.0
pair_coeff 1 3 1550 .297 0.0

pair_style buck/coul/long 15.0
pair_coeff 3 3 16732 0.213 3.47

kspace_style ewald 1.0e-4

Morse Potential

pair_style morse 2.5
pair_coeff 2 3 4.71 3.8 1.18

Three Body Potential

angle_style hybrid harmonic
angle_coeff 1 harmonic 1.69 120.0
ERROR: Numeric index is out of bounds (force.cpp:666)

angle_style hybrid harmonic
angle_coeff 1 harmonic 1.69 120.0
ERROR: Numeric index is out of bounds (force.cpp:666)

Since you didn;t read a data file you allocated no space

for angles to be defined. You can do this with “angle”

optional keywords added to the create_box command.

But how are you going to define what atoms are in

what angle interactions (I,J,K) w/out a data file?


Dear Steve,

Thanks for your attention.

Could you please help me, have can I create read data file?


There are many example data files in the examples dir.

There are many builder tools mentioned on the PrePost web

page of the LAMMPS site.


Dear Steve,

thanks, But I can not understand how many bonds, angles and dihedrals are in one crystal of CACO3 (include 6 Ca, 6 C and 18 O)?

could you please help me?



I can’t help you enumerate them. Maybe some pre-processing tool,

i.e. a builder, can. See the PrePost doc page on the web site

for ideas.