QUESTION about transfer a force filed from GULP to LAMMPS

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

thanks for the reply, Axel. i use a simple model (the model contains harmonic bond and angle term, couloum and l-j term), to test the consistency bewteen GULP and lammps, as showed in the picture 1. The atomic position are the same since it is dumped from GULP as final optimization geometry, but the per atom force are not. and when i check energy calculated by individual potential terms, it shows that the energy associated with bond term, angle term, and l-j term are matched between two codes, but the electrostatic energy differs about 0.05 ev, (which is the same situation when i work on the model that i really interested in ), i think the per atom force difference may come from the electrostatic term, but i have check the ewald summation accuracy and real space cutoff, both are the same with the parameters in GULP. any ideas about what should i check further? i attached my input script for lammps. thanks very much for your help!!!

And for that issue about improper energy which you pointed out “not scientific”, i realized that in my lammps data file, each SiO4 tetrahedral has 12 improper, 8 of which are redundant, so that the calculated improper energy becomes three time the real energy.( see below, the fourth column are Si atoms, each Si are included in 12 improper ) . i generate these impropers in moltemplate with the follwing command (with atom1 =Si, atom 2= O)

write_once(“Data Impropers By Type (class2_imp.py)”) {
@improper:type1 @atom:type1 @atom:type2 @atom:type2 @atom:type2 @bond:* @bond:* @bond:*
}

does anyone have a better way to define these class2_impropers?

Best Regards

Jiasen Guo

impropers

pastedImage.png

in.MZHB_SOD_minization (1.12 KB)

pciture1.png

2

I am CCing Stan who is our electrostatics expert. You are saying

that for the same atom config, same box, that GULP vs LAMMPS is

giving a different answer for only the electrostatics term? Please

provide the details for how GULP computes it, and the input script

lines for how you are asking LAMMPS to compute it. Also

the thermo output for the snapshot from both codes. Are you certain

that you are treating the exclusion terms the same in both

codes? In LAMMPS this is the special_bonds command for

excluding the Coulombic interactinos in 1/2, 1/3, 1/4 interactions.

Also in LAMMPS, the total Coulombic energy is the sum of short-range

Coul + long-range Coul

  • is that what you are comparing to GULP?

Steve

pastedImage.png

3

A question:
    Why are you using an improper interaction on SiO4? I was under
the impression that, in traditional ("class1") force-fields improper
interactions are only used to enforce planarity. SiO4 is definitely
not coplanar.
    Is this different for class2 force-fields?

Regarding the class2 impropers, you can solve your problem in two ways:

1) For each type of improper interaction, divide the "K" coefficients
by the number of redundant dihedrals generated by moltemplate (6 in
your case, I think). This is actually what a lot of molecule builders
do. (If you are interested why, see the discussion below)

2) Use different rules for generating improper interactions to reduce
or eliminate the redundancy for your particular molecule. There are
several ways to do this. Try downloading the attached file
("nbody_Impropers_Jcenter_swapIKL.py"), copying it to your
"moltemplate/src/nbody_alternate_symmetry/" subdirectory, and
modifying your .LT file:

write_once("Data Impropers By Type (nbody_Impropers_Jcenter_swapIKL.py)")
{
    @improper:type1 @atom:type1 @atom:type2 @atom:type2 @atom:type2
}
(You can omit the "@bond:* @bond:* @bond:*")
This will generate only 1 improper interaction for every central atom,
regardless of the order of (I,K,L)
If you are curious how this differs from the "class2_imp.py", look
inside the "moltemplate/src/nbody_alternate_symmetry/" subdirectory
and compare the scripts.

If that eliminates too many interactions, then try this instead:
write_once("Data Impropers By Type")
{
    @improper:type1 @atom:type1 @atom:type2 @atom:type2 @atom:type2
}
This uses the default rules for creating improper interactions (which
are defined in "moltemplate/src/nbody_Impropers.py")

I have not tested the script I just sent you.
Free to tinker with these scripts, and let me know what works.

---- Discussion: ----

SiO4 may be symmetric, but more generally, molecules will have
different atoms at each site. The formula for the class2 improper
interaction energy has two terms: Ei and Eaa:
http://lammps.sandia.gov/doc/improper_class2.html
The first term (Ei) is invariant with respect to changing the order of
any of the 3 atoms (I,K,L) surrounding the central atom (J). Unless
I'm really confused, it appears that second term (Eaa) is not:

Eaa = M1(θijk-θ1)(θkjl-θ3) + M2(θijk-θ1)(θijl-θ2) +M3(θijl-θ2)(θkjl-θ3)

...where J is the central atom and θijk is the angle between the J-I
bond and the J-K bond. Perhaps I'm dense, but it seemed that any
modification of the atom order could potentially change the energy Eaa
(unless the parameters M1=M2=M3, or θ1=θ2=θ3). So to make sure that
moltemplate works with ALL molecules, the "class2_imp.py" script you
are using generates new improper interactions for every possible
permutation of (J,K,L) atoms even if this means in many cases the
results are terribly redundant. (This redundancy is usually absorbed
into the coefficient for that interaction. I admit this is not very
efficient. In any case, picking that coefficient is the
responsibility of the force-field author or converter ... you, in this
case.)

I hope this helps.

Andrew

nbody_Impropers_Jcenter_swapIKL.py (1.62 KB)