OPLS-aa + Tersoff for fullerene(C60)-derivative PCBM

Jason · July 15, 2013, 9:33pm

hi lammps users,

I am modeling the PCBM system, which is a C60 (atom-type 1, atom-id 1-60) and a 28 atom “satellite” molecule (atom-type 2-11, atom-id 61-88), see attached pcbm.png.

The atom-type 1, atom-id 1-60 I want to model using the Tersoff potential.

The atom-type 1, atom-id 1,2 are bonded to atoms in the “satellite” molecule through harmonic bonds:

16 7 61 1
17 7 61 2

3-body angle harmonic bonds:

22 10 1 61 70
23 10 2 61 70
24 11 1 61 2

47 19 1 61 70
48 19 2 61 70
49 20 1 61 62
50 20 2 61 62

55 25 2 1 61
56 25 6 1 61
57 25 9 1 61
58 25 1 2 61
59 25 3 2 61
60 25 12 2 61

ad 4-body dihedral bonds:

1 20 1 61 70 71
2 20 2 61 70 71
1 21 1 61 70 71
2 21 2 61 70 71
1 22 1 61 70 71
2 22 2 61 70 71
1 23 1 61 70 80
2 23 1 61 70 81
3 23 2 61 70 80
4 23 2 61 70 81
1 24 2 1 61 62
2 24 6 1 61 62
3 24 9 1 61 62
4 24 1 2 61 62
5 24 3 2 61 62
6 24 12 2 61 62
1 25 2 1 61 70
2 25 6 1 61 70
3 25 9 1 61 70
4 25 1 2 61 70
5 25 3 2 61 70
6 25 12 2 61 70

I see that the special_bonds command does not support Tersoff:

http://lammps.sandia.gov/doc/special_bonds.html

"

Is there anyway to enforce the Tersoff interaction between all atom-type 1 while also retaining the bonding term with the satellite molecule without significantly changing my input structure format?

When I run the attached structure and script the C60 atoms “explode”. I’m guessing this because either:

a) the Tersoff potential is not being applied to atom-type 1, atom-id 1-60 ecause these atoms are part of molecule-id 1.
or
b) the non-bonded lj/cut and/or coul/wolf are turned on between atom-type 1, atom-id 1-60.

I still need the non-bonded (lj/cut and coul/wolf) interactions to apply to intra-C60 atoms, so that I could simulate crystalline C60 for example.

If I cannot keep a similar input structure, it is not obvious to me how to set the potentials up properly.

I apologize if the solution for this is obvious with more careful reading of the documentation.

Thanks,

Jason

C.tersoff (734 Bytes)

fort.26 (13.4 KB)

lmp.in.anneal (5.89 KB)

sjplimp · July 16, 2013, 1:27pm

Why do you think it is a good idea to model C60
with permanent bonds and angles as well as
a Tersoff potential? What is the Tersoff potential
doing if the structure is already bonded together?
Is this a model from the literature or your own
creation?

Steve

Jason · July 16, 2013, 2:02pm

Hi Steve,

I am modeling the C60 Carbons (C1 or atom-type 1) using Tersoff because I don't have a program which can identify the unique bonded harmonic, angle harmonic, and dihedral terms for the 5- and 6-member rings. As Axel told me "this is a non-trivial task". Tersoff handles this automatically.

It may not be obvious from the input files, but the only permanent-bonded interactions that the C1 atoms experience is with the first few atoms in the satellite molecule (which I pointed out in the first message).

So, this potential is my own creation, but one which I have successfully optimised and performed lattice dynamics on using GULP, see:

https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/gulp.gin
https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/olps.lib

I guess I just need to know if this is the correct approach in LAMMPS, or if I have to modify my input, perhaps with one more atom-type to identify the 2 C1 atoms which are permanently bonded. I was hoping to avoid this as it complicates my transfer from GULP.

Thanks for the reply,

Jason

sjplimp · July 17, 2013, 2:11pm

ok - now I understand your geometry better.

Is there anyway to enforce the Tersoff interaction between all atom-type 1 while also retaining the bonding term with the satellite >molecule without significantly changing my input structure format?

yes, you should be able to use pair hybrid with Tersoff
only applied to type 1, then apply other pair styles
to the other types or to interactions bewteen
type 1 and N (N>1). Bonds/angles will simply
be computed in addition to the pair styles.
So as long as there are no bonds between two type 1
atoms, that should be fined.

The pair hybrid doc page has a long
section on how to do those kinds of things with
many-body potentials. It starts with this comment:

IMPORTANT: Several of the potentials defined via the pair_style
command in LAMMPS are really many-body potentials, such as Tersoff,
AIREBO, MEAM, ReaxFF, etc. The way to think about using these
potentials in a hybrid setting is as follows.

Steve

Jason · July 24, 2013, 6:43am

ok - now I understand your geometry better.

>Is there anyway to enforce the Tersoff interaction between all atom-type 1 while also retaining the bonding term with the satellite >molecule without significantly changing my input structure format?

yes, you should be able to use pair hybrid with Tersoff
only applied to type 1, then apply other pair styles
to the other types or to interactions bewteen
type 1 and N (N>1). Bonds/angles will simply
be computed in addition to the pair styles.
So as long as there are no bonds between two type 1
atoms, that should be fined.

So I believe I need to have Tersoff and bond terms acting between atom-id=1,2 (atom-type=1), both atom's are of type 1.

I have read through the pair hybrid doc page, and I think I have the correct input.

I need 2 of the Tersoff-modeled C's (atom type=1) to be bonded together with harmonic angle and dihedral terms. This is necessary to connect my satellite molecule to the C60, see :

https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_duda.png

The contents of isolated.tar.gz are:

https://github.com/jasonlarkin/disorder/tree/master/pcbm/pcbm_pdb/isolated/anneal

and

https://github.com/jasonlarkin/disorder/tree/master/pcbm/pcbm_pdb/isolated/anneal_nobond

For anneal, I include bonds between atom-id=1,2 (atom-type=1). The result is unstable with large negative energy:

https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/anneal/log.lammps

For anneal_nobond, I remove bonds between atom-id=1,2 (atom-type=1). The result is unstable with large positive energy:

https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/anneal_nobond/log.lammps

For both anneal and anneal_nobond, GULP produces a reasonably stable molecule:

https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/anneal/freq.gout
https://github.com/jasonlarkin/disorder/blob/master/pcbm/pcbm_pdb/isolated/anneal_nobond/freq.gout

I can also confirm that GULP MD produces reasonably stable structures for both.

isolated.tar.gz (104 KB)