[gmx-users] Force field for graphene sheet: deformation at the boundaries

[Kathleen Kirchner]

Dear colleagues,

I came across a puzzling behaviour when simulating graphene sheets 
(periodic in xy).

When running a simulation of just one position restrained graphene sheet 
in a large vacuum slab with periodic_molecules=yes and pbc=xyz,  I found 
a static deformation of the boundaries of the graphene sheet.

On search for a more suitable force field I finally found the deletion 
of angle/dihedral information from the topology the best result.

Why is a large graphene sheet with the same force field information for 
all atoms, bonds, angles and dihedrals deformed only at the border of 
the periodic simulation box? In case of errors in the force field I 
would have suspected bents all over the graphene sheet.


Some information on the structure and topology preparation:

The initial structure has been generated with the help of of 
ase.structure 
(https://wiki.fysik.dtu.dk/ase/epydoc/ase.structure-module.html). A 3-4 
armchair nanoribbon resulted in a graphene sheet of 6 nm * 6 nm area.

Afterwards I run g_x2top to obtain the topology:

g_x2top -f Gl1.gro -o out.top -name Gl1 -pbc

with the addition of the line

C   opls_145   0.0   12.011   3   C   0.140   C   0.140   C   0.140

in gromacs4.5.5/share/gromacs/top/oplsaa.ff/atomname2type.n2t

and added position restraints in all directions for all carbon atoms 
(k_xyz=1000 kJ mol^-1 nm^-2).

Best regards
Kathleen

Kathleen Kirchner
PhD student
Max Planck Institute for Mathematics in the Sciences
(MPI f. Mathematik in den Naturwissenschaften)
Inselstr. 22-26, D04103 Leipzig
e-mail: kirchner at mis.mpg.de
web: http://www.mis.mpg.de/scicomp/CompPhysChem/
Tel +49 341 9959 725
Fax +49 341 9959 999