QUESTION about transfer a force filed from GULP to LAMMPS

[_Guo_Jiasen_MU-Stude]

Dear lammps users,

I am trying to transfer a force field described in GULP to lammps, and i run MD simulation with the same parameters in both the two package, using the same post pocessing tool, but the result spectrum is very different, see the following two spectrum, the one on the right is from GULP, and the other one from lammps. i really don’t know what will be the reason for the difference. someone please help me. thanks very much!!!

Best Regards

Jiasen Guo

the force field specification in GULP follows,

harmonic bond
Si O 33.57535587 1.61

three bond
Si O O 4.8678025239 109.47
O Si Si 0.748892696 149.8

bcross bond
Si O O 1.1857467686 1.61 1.61 0 1.8 0 1.8 0 3.6
O Si Si 1.5601931166 1.61 1.61 0 1.8 0 1.8 0 3.6
bacross bond
Si O O 0.8113004207 0.0000 1.61 1.61 109.47

xangleangle bond
Si O O O -0.6864849713 -0.6864849713 -0.6864849713 109.47 109.47 109.47

I attach the GULP doc for specific terms, hoping it helps you understand the way GULP defines force field.

bond cross term

bond angle cross term

angle cross term

the force field specification in lammps follows,

bond_style harmonic
bond_coeff 1 16.787678 1.61

angle_style class2
angle_coeff 1 109.47000 2.4339013 0.00000 0.0000
angle_coeff 1 bb 1.1857467686 1.61 1.61
angle_coeff 1 ba 0.8113004207 0.0000 1.61 1.61
angle_coeff 2 149.80000 0.3744463 0.00000 0.0000
angle_coeff 2 bb 1.5601931166 1.61 1.61
angle_coeff 2 ba 0.000000000 0.0000 1.61 1.61

improper_style class2
improper_coeff 1 0.0 100
improper_coeff 1 aa -0.6864849713 -0.6864849713 -0.6864849713 109.47 109.47 109.47

pair_style coul/long 12.0
pair_coeff * *
kspace_style ewald 1.0e-4
special_bonds coul 1.0 1.0 1.0

details from lammps doc.

angle_style class2

improper_style class2

[sjplimp]

No one is in a better position to “debug” the difference between

the 2 codes than you. I suggest you take a single snapshot

and compare the energy/pressure that LAMMPS vs GULP produce.

For LAMMPS, there are various thermo keywords that correspond

to individual terms in the energy expression, at nearly the

level of detail of your eqs. I assume there is something similar

in GULP. So compare the energy outputs term by term and

figure out where the difference(s) are.

Steve

[athomps]

You are wasting your time comparing spectral densities from the two codes. First you need to convince yourself that both codes are running the same interatomic potential. The first step is to pick a simple potential that is standard in both codes, such as Stillinger-Weber silicon, and verify that you can build a crystal in both codes, run for zero timesteps, print out the forces and obtain a numerical match to very high precision. The reason for doing this is that there are countless opportunities for discrepancies to occur even in a very simple comparison, and you will eliminate them quicker if you know the codes are supposed to agree. You should then repeat this process with your new potential. It is likely that you have made some mistakes in your definition of the GULP potential in LAMMPS. Only after you have eliminated those mistakes should you start comparing results from simulations. Even there, it is possible to make mistakes, due to differences in units, etc. Only after you can get simple results to match from both codes should you attempt to compare the spectral densities.

Aidan