Mesoporous silica simulation

Hello,

I am trying to simulate MCM-41 (a mesoporous silica) using LAMMPS but am encountering some problems. I have two structures of MCM-41 one that I have prepared and the other is an optimized MCM-41 structure available from Piero Ugliengo. Both structures have been optimized using GULP as well as the Ugliengo model using DFT. I’m using the potentials for hydroxylated silica from Schroeder but have also tried the FFSiOH potentials. I have run using 1 processor however the same errors are present when using more. However in all cases the simulation runs for a few (10-40) time steps before giving errors. I’ve visualized the system in VMD which shows that certain atoms fly off from the structure giving errors such as: Lost atoms, bond atoms missing or angle atoms missing. The output also shows very large temperatures but whether these cause the atoms to fly off or if they as a result of the large movement im not sure.

Step Temp E_pair E_mol TotEng Press
0 0 44577.249 6302.8849 50880.134 -207049.6
10 16755.925 36989.993 6770.228 56857.257 -217308.19
20 8795890.6 -4481.2057 3528100.5 10398805 16563180
ERROR on proc 0: Bond atoms 25 381 missing on proc 0 at step 30 (…/neigh_bond.cpp:65)

I’m pretty sure the structures and potentials are ok as I’ve previously run the same structures with the same potentials using DL_POLY which run fine.

I’ve attached an input file to see if im missing some LAMMPS specific commands. I’m using LAMMPS version 7 Dec 2015 version.

Thanks for any help.

Daniel

units metal

dimension 3

boundary p p p

atom_style hybrid angle charge

atom_modify sort 0 0

read_data mcm.lmpdat

group Siatoms type 1

group Hatoms type 2

group Oatoms type 3

group Ohatoms type 4

set group Hatoms charge 0.426

set group Oatoms charge -2

set group Siatoms charge 4

set group Ohatoms charge -1.426

kspace_style ewald 1e-6

pair_style hybrid buck/coul/long 10

pair_coeff 1 3 buck/coul/long 1283.9 0.3205 10.6615 #O Si

pair_coeff 1 4 buck/coul/long 983.5566 0.32052 10.66158 #Si Oh

pair_coeff 3 3 buck/coul/long 22764 0.149 27.88 # O O

pair_coeff 4 4 buck/coul/long 22764 0.149 27.88 # Oh Oh

pair_coeff 2 3 buck/coul/long 396.97 0.25 0 # H O

pair_coeff 3 4 buck/coul/long 22764 0.149 13.94 # O Oh

pair_coeff 2 4 buck/coul/long 312 0.25 0 #H Oh

#pair_coeff 2 4 morse 7.05 3.1749 0.942

pair_coeff 1 2 none # H Si

pair_coeff 2 2 none # H H

pair_coeff 1 1 none # Si Si

pair_modify table 0

bond_style hybrid harmonic morse

bond_coeff 1 harmonic 25 1.626

bond_coeff 2 morse 7.05 3.1749 0.942

angle_style harmonic

angle_coeff 1 2.1 109.47

angle_coeff 2 2.1 109.47

timestep 0.001

thermo 10

dump positions all atom 10 history.lammpstrj

fix 1 all nve

run 10000

Hello,

I am trying to simulate MCM-41 (a mesoporous silica) using LAMMPS but am encountering some problems. I have two structures of MCM-41 one that I have prepared and the other is an optimized MCM-41 structure available from Piero Ugliengo. Both structures have been optimized using GULP as well as the Ugliengo model using DFT. I’m using the potentials for hydroxylated silica from Schroeder but have also tried the FFSiOH potentials. I have run using 1 processor however the same errors are present when using more. However in all cases the simulation runs for a few (10-40) time steps before giving errors. I’ve visualized the system in VMD which shows that certain atoms fly off from the structure giving errors such as: Lost atoms, bond atoms missing or angle atoms missing. The output also shows very large temperatures but whether these cause the atoms to fly off or if they as a result of the large movement im not sure.

Step Temp E_pair E_mol TotEng Press
0 0 44577.249 6302.8849 50880.134 -207049.6
10 16755.925 36989.993 6770.228 56857.257 -217308.19
20 8795890.6 -4481.2057 3528100.5 10398805 16563180
ERROR on proc 0: Bond atoms 25 381 missing on proc 0 at step 30 (…/neigh_bond.cpp:65)

I’m pretty sure the structures and potentials are ok as I’ve previously run the same structures with the same potentials using DL_POLY which run fine.

But that doesn’t mean, there cannot be any mistakes in translating the geometry to lammps or that there are typos or mistakes converting potential parameters.

Axel

It also looks as if your totaleng is growing by 2x in 20 steps. So you are

not conserving energy with NVE and a 1 fmsec timestep. I would turn

off parts of your potential (start as simple as possible) and build up

to see where things are going wrong.

Steve