I am trying to set up LAMMPS simulation to relax twisted bilayer graphene using Drip potential.
As it says Drip only can be applied to interlayer interactions, I have used tersoff potential for intralayer.
This LAMMPS input script to relax twisted bilayer graphene structure
basic parameters
units metal
atom_style full
boundary p p p
newton on
processors 8 8 1
dimension 3
fix frelax all nve
fix controltemp all temp/rescale 100 10.0 10.0 10.0 1.0
run 1000
unfix controltemp
unfix frelax
undump Graphene
The error says “ERROR on proc 40: No enough neighbors to construct normal. Check the configuration to see whether atoms fly away. (src/USER-MISC/pair_drip.cpp:696)”
Can you please guide me on how to fix this problem and relax this system.
Thank you.
I think the problem is a combination of bad cell vectors and input geometry. Plus, the manual says to
assign different molecular IDs to atoms in different layers, [as] DRIP is implemented such that
only atoms with different molecular ID can interact with each other.
while in your DATA file the same molecule ID is used for both layers. You better start with a different initial guess, where the layers are perpendicular to a Cartesian axis (just to make life easier) and periodic boundaries are set to the actual graphene supercell.
Thank you very much for your reply.
I have generated these two twisted layers with a python code. Visualize it with Vesta. Import xyz coordinates and convert into lammps input file using atomsk.
Can you please explain whether I can use Vesta to arrange my struture perpendicular to a Cartesian axis.
Thank you so much
The crucial required change is to make certain that each layer has a different molecule ID.
Without that, there is no meaning in trying to realign.
The second thing you must do is to check, if the input will run with just 1 processor.
You are running a system with only 1324 atoms on 64 CPUs. that can create problems with the computation of the orientation of the layes in the interlayer potential, since each domain may own too small a piece of the layer.
I have changed the ID of the mol.
this is the output file after defining two layers.
LAMMPS (3 Mar 2020)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (0 0 0) to (83.5081 83.5081 83.5081)
8 by 8 by 1 MPI processor grid
reading atoms …
1324 atoms
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.0036819 secs
read_data CPU = 0.013231 secs
Reading potential file C.drip with DATE: 2019-04-19
Reading potential file SiC.tersoff with DATE: 2011-04-26
WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (src/min.cpp:190)
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 17.7
ghost atom cutoff = 17.7
binsize = 8.85, bins = 10 10 10
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair drip, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(2) pair tersoff, perpetual, copy from (1)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Setting up cg style minimization …
Unit style : metal
Current step : 0
ERROR on proc 40: No enough neighbors to construct normal. Check the configuration to see whether atoms fly away. (src/USER-MISC/pair_drip.cpp:696)
Last command: minimize 1e-4 1e-6 1000 100000
You are still running your input script on 64 CPUs. As suggested by Axel, comment out the processors command and run on a single core.
If the error remains, it could be due to the dangling bonds in the graphene layers.
To create you input, start from the crystallographic structure of graphite (tipically a CIF file) and replicate it with a program such as GDIS. Then make a supercell and export the coordinates in Cartesian format (XYZ or PDB). Make sure to use the correct PBC vectors.
LAMMPS (3 Mar 2020)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (0 0 0) to (83.5081 83.5081 83.5081)
1 by 1 by 1 MPI processor grid
reading atoms …
1324 atoms
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000568867 secs
read_data CPU = 0.0363438 secs
Reading potential file C.drip with DATE: 2019-04-19
Reading potential file SiC.tersoff with DATE: 2011-04-26
WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (src/min.cpp:190)
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 17.7
ghost atom cutoff = 17.7
binsize = 8.85, bins = 10 10 10
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair drip, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(2) pair tersoff, perpetual, copy from (1)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Setting up cg style minimization …
Unit style : metal
Current step : 0
ERROR on proc 0: No enough neighbors to construct normal. Check the configuration to see whether atoms fly away. (src/USER-MISC/pair_drip.cpp:696)
Last command: minimize 1e-4 1e-6 1000 100000
Then the problem is a bogus input structure. For instance, each plane is composed of 662 atoms, which is not a even division of the graphite structure.
I attach a DATA file that actually works. It is your job to learn crystallography and how to handle crystal structures. graphite.data (48.3 KB)
Thank you very much for all the feedbacks. I will learn how to generate a proper input structure. Meantime if possible please refer me some tutorial or notes which describe generate twisted bilayer graphene. (or any other twisted systems are also ok)
It is not clear to me, how this “twist” can be realized properly.
Have you tried running your system without the interlayer potential added?
Does it maintain the expected structure?
When I visualize your data file, I see two small graphene leaflets in a huge box. When I look at the graphene sheet and interlayer potential examples, I see graphene sheets with proper periodic continuation across boundaries. Without that periodicity, the sheets will likely deform and crumble or do a variety of other strange changes, once you add some kinetic energy.
But then again, with periodic boundaries, it is not possible to have an arbitrary in plane rotation of one sheet against the other. Only selected angles are possible and which ones are possible depends on the system size as the rotated system also needs to have perfect periodic continuation.
I agree what you are saying.
I used the following equation in the PRB 86 155449 (2012) to generate the twisted bilayer graphene system. Then convert in to lammps input file using atomsk software. I understand my there is some problem with periodic boundary
[quote=“akohlmey, post:11, topic:39978”]
Only selected angles are possible and which ones are possible depends on the system size as the rotated system also needs to have perfect periodic continuation.
condition as you were saying. box size is too large.
“Only selected angles are possible and which ones are possible depends on the system size as the rotated system also needs to have perfect periodic continuation.” I understand what you are saying. Can you please guide me how to make proper lammps input coordinates file for twisted bilayer graphene.
I have run the system with just tersoff potential. (without Drip) The minimization finished without giving errors. log.lammps (8.5 KB) tBLGm10.xyz (400.3 KB)
Sorry, I don’t have time to work this out and then teach you. It is easy enough for me to spot the inconsistencies, but figuring out a workflow is a different thing, especially when it is not my research project.
Key problem is that you seem to unconditionally trust some software to “do the right thing™” and in some way they do, but - for example - the box dimensions your data file has are not suitable for proper periodic boundary continuation.
I strongly suggest you develop your skills in stages. First you need to be able to reproduce a (non-twisted) data file like it is used in the “drip” example in the LAMMPS distribution. Then (and only then) you can move on to try and figure out how to build a suitable data file with allowed in plane rotations that are consistent with your needs. Please note that the problem with allowed rotation angles is equivalent to that of building surfaces that are not of the Miller indices [1 0 0], [0 1 0], or [0 0 1]. Those usually require using supercells with suitable multiples of the corresponding primitive cell in different directions.
Minimization is not a good measure since it is equivalent to observing the system at 0K. Nothing much will happen. Also, the fact that a simulation completes says nothing about its correctness. Calculations can complete and still be completely bogus.
I will learn how to generate a correct input data file with correct periodic box size.
As relaxation of the graphene system means energy minimization, how can one be sure whether the system is in relaxed state. I though by checking the minimized potential energy one can do it.
For instance with tersoff potential I got minimum potential energy per atom=-7.0162 eV
