I am confused in Getting charge per atom on the edge of graphene sheet

LAMMPS LAMMPS Beginners
,
1

Hello All,
I tried to calculate charge per atom using ReaxFF Potential where it reads data from Graphene.dat file


1 atom types

0 149.1 xlo xhi
0 88.5424 ylo yhi
-10 10 zlo zhi

Masses

1 12.0107

Atoms

1 1 0 0 0 0
2 1 0 0.71 1.22976 0
3 1 0 2.13 1.22976 0
4 1 0 2.84 0 0
5 1 0 0 2.45951 0
6 1 0 0.71 3.68927 0
7 1 0 2.13 3.68927 0
8 1 0 2.84 2.45951 0`

It’s perfectly regular and continued across periodic boundaries with ATOM ID, MOL ID , ATOM TYPE, CHARGE X Y Z I had to add charge column because it was giving me a bad formatting error.

This is my reaxff input file

units          real
dimension      3
boundary       p p f
atom_style     charge

#---------------ATOM DEFINITION------------------------------
read_data      graphene.dat 

#---------------FORCE FIELDS---------------------------------
pair_style     reax/c NULL
pair_coeff     * * ffield.reax C

#---------------CHARGE EQUILIBRATION------------------------- 
fix            qeq all qeq/reax 1 0.0 0.0 1e-6 reax/c

#--------------------SETTINGS--------------------------------
timestep       0.001 

#---------------COMPUTES AND VARIABLES-----------------------
region         left block 0.0 5.0 0.0 5.0 -0.5 0.5
group          left_region region left

compute        charge left_region property/atom q

compute        q1 left_region reduce ave c_charge
variable       qtot equal count(left_region)*c_q1

thermo_style   custom step c_q1 v_qtot pe
#---------------TIME INTEGRATION-----------------------------
fix            nve all nve

thermo         100
run            1000

print "All done"

fix qeq all qeq/reax 1 0.0 0.0 1e-6 reax/c

I tried to keep equilibration as low as possible (it can be changed to high value)

time step is 0.001 ( was mentioned in documentation)

I grouped only left part because I wanted calculate charge behaviour on the edge ( a research paper mentioned that at centre its nearly zero and on the edge it’s around o.13e)

region         left block 0.0 5.0 0.0 5.0 -0.5 0.5
group          left_region region left

compute        charge left_region property/atom q

compute        q1 left_region reduce ave c_charge
variable       qtot equal count(left_region)*c_q1

thermo_style   custom step c_q1 v_qtot pe
#---------------TIME INTEGRATION-----------------------------
fix            nve all nve

thermo         100
run            1000

print "All done"

I directly asked to compute charge/atom which it did but results were like this

Step          c_q1             v_qtot                  PotEng
         0  -4.6018744e-14 -6.9028117e-13 -967733.14
       100   5.884182e-15   8.826273e-14  -967733.14
       200  -2.8199665e-14 -4.2299497e-13 -967733.14
       300  -6.0396133e-14 -9.0594199e-13 -967733.14
       400  -7.6827433e-14 -1.1524115e-12 -967733.14
       500  -4.0245585e-14 -6.0368377e-13 -967733.14
       600   1.4988011e-15  2.2482016e-14 -967733.14
       700   5.2347016e-14  7.8520523e-13 -967733.14
       800  -1.3822277e-14 -2.0733415e-13 -967733.14
       900   2.8754776e-14  4.3132165e-13 -967733.14
      1000  -1.6042723e-14 -2.4064084e-13 -967733.14

The issues i have with this is it doesn’t give me constant value of q/atom ( to draw some conclusion )

I tried changing my code by using airebo’s Potential ( its mostly utilized for hydrocarbon )

units          real
dimension      3
boundary       p p f
atom_style     charge


#---------------ATOM DEFINITION------------------------------
read_data      graphene.dat

#---------------FORCE FIELDS---------------------------------

pair_style 	airebo 3.0
pair_coeff     * * CH.airebo C


#--------------------SETTINGS-------------------------------------
timestep       0.001

#---------------COMPUTES AND VARIABLES-----------------------
group          type1 type 1
compute        charge1 type1 property/atom q

compute        q1 type1 reduce ave c_charge1
variable       qtot equal count(type1)*c_q1



thermo_style   custom step c_q1 v_qtot pe
fix            nve all nve

thermo       100
run          1000

print "All done"

butt it gave me initial charge value that was given to it in data file

         0   0              0             -37339.066
       100   0              0             -37339.066
       200   0              0             -37339.066
       300   0              0             -37339.066

I Don’t know how to proceed further, should I be using hybrid potential ? Feel free to correct my code since am new to this things.

2

These numbers are essentially zero and just representative of the numerical noise for double precision floating point numbers that will come out of the charge equilibration process that is part of the ReaxFF model. Note that you only converge to a finite precision.

Given the structure and geometry of your system, it makes no sense to determine partial charges, since there is nothing that would induce a polarization and thus it is expected that those charges will be zero (or very close to zero).

The AIREBO potential (as can be seen from the publications describing it) does not include any explicit charges so it will neither read, use, or modify the per-atom charges you set. You can set those charges to random values and the result of the force calculation will be the same.

This is not a problem of the LAMMPS code or your input, but a problem of understanding the physics of the models you are using. Using a hybrid potential (how? which?) will likely only make the mess worse. At any point this kind of discussion is not a subject for this forum, but more a topic you need to discuss with your adviser or tutor that is advising/training your in doing simulations and is familiar with your research (and cares about the results you produce).

3

Even if change the geometry of my graphene(I made a triangular hole in the middle) sheet it gives the same distorted numerical noise.

Initially my instructor told me to prove that symmetric graphene is neutral and doesn’t contain any partial charge but my main work is with non-symmetric graphene where I would analyse its partial charge but I am struggling to get proper potential which could get me non zero charge value. Apart from the potential, Is my other code is well enough to get me what I want ?

PS: My data file contains only carbon it doesn’t contain hydrogen at edge will it bother graphene properties

4

This forum is not an input file review & correction service.

It seems to me that you don’t understand how ReaxFF and charge equilibration works and thus are having the wrong expectations. For as long as you have a system with the same element everywhere, where should the polarization come from?

At any rate, what you are asking is not a question about LAMMPS, but about the science of your simulations and thus you need to discuss with your adviser, not us.