multi-harmonic potential for dihedral type

Hi,

Just checking if anyone has extended the multi-harmonic dihedral code in Lammps to include the 6th power of cos (Angle) for implementating the
Ryckaert-Bellemans cosine expansion form? The current one in Lammps lacks the 6th power term.

The energy missing from the total dihedral potential is about 1/3 of what it should be, based on the coefficients given in Ryckaert-Bellemans with a dihedral angle of ~127 degrees for a simple 4-carbon chain.

Thanks

william

Not that I'm aware but if you extend the existing potential
please send it to us and we will include it in the distro.

Steve

Sigh. I thought about doing this, but I was too lazy to write a
DihedralXXXX::coeff() function to read in the parameters.

    --- details ---
I wanted to do it for a general expansion with sines and cosines and
phi0. (Not an expansion that stops at the 4th or 6th or Nth term.)
For this reason, one has to be careful because I think that the
strings in the "char **arg" argument of the coeff() function have
restricted length (80 characters?). This would probably force us to
read the coefficients of the expansion from an external file.

Hi Andrew and Steve,

Thanks for your quick replies.

I guess I will try out the dihedral table and later create that potential as an exercise.

Just a small question for using the dihedral table, how fine should the interval be? I have thoughts of using 0.01 degree intervals but that seems to be an overkill. Any advices?

Thanks

william.

0.01 degrees sounds like plenty of resolution

Steve

Hello Steve and Andrew,

Thanks for your advices.

I was thinking if energy conservation (5 to 6 sig. fig.) might be a problem with insufficient resolution due to linear interpolation. Spline fits although more accurate might have some “wandering” effect near the interval limits but i guess it should not be too much a problem.

I am trying it out. I have also modified the multi-harmonic code to take in the 6th coefficient but i am still testing it.

Thanks

william

Forgive me if I posted this reply twice. (The first time did not go
through, it seems.)

I was thinking if energy conservation (5 to 6 sig. fig.) might be a problem
with insufficient resolution due to linear interpolation. Spline fits
although more accurate might have some "wandering" effect near the interval
limits but i guess it should not be too much a problem.

Oh I see.

That's an really good point. I missed that.

If you are worried about this, then I recommend that you use the "NOF"
option with "spline". The "NOF" option does not use separate splines
to calculate the energy and the forces. It uses the same spline (with
the same control points) and calculates the derivative of the spline
function to get the forces. (The two relevant functions are
"cyc_splintD()" and "cyc_splint()" in "dihedral_table.cpp"). This
way, the force at every angle is guaranteed to be the exact derivative
of the potential. (So you never have to worry about inconsistency
between the forces and energies, regardless of the table resolution.)
Of course, you still have to use a table with sufficiently high enough
resolution to approximate the multi-harmonic function you want to use.
(For this purpose, you can probably get away with a table with a few
hundred entries, but feel free to experiment and benchmark.)

I am trying it out. I have also modified the multi-harmonic code to take in
the 6th coefficient but i am still testing it.

Great.

I just realised that users can define an arbitrary Fourier expansion
(Ryckaert-Bellemans expansion) to approximate any dihedral potential
by defining multiple dihedral interactions between the same 4 atoms
(21, 22, 23, and 23 in the example below). (If you use something like
"dihedral_style charmm", then the resulting function does not have to
be an even function.) Example:

-- data file --
Data Dihedrals

1 1 21 22 23 24
1 2 21 22 23 24
1 3 21 22 23 24
1 4 21 22 23 24
1 5 21 22 23 24
1 6 21 22 23 24
:

# -- Input Script (init section) --
dihedral_style charmm

:

# -- Input Script (settings section) --
# dihType K n d w(set 0)

dihedral_coeff 1 -1.0 1 -90 0
dihedral_coeff 2 0.5 2 180 0
dihedral_coeff 3 -0.333 3 -270 0
dihedral_coeff 4 0.25 4 360 0
dihedral_coeff 5 -0.2 5 -540 0
dihedral_coeff 6 0.167 6 720 0
:

http://lammps.sandia.gov/doc/dihedral_charmm.html
http://en.wikipedia.org/wiki/Sawtooth_wave (I hope)

The only drawback I can think of is that the dihedral calculations
would be slower. (LAMMPS would have to recompute the same overhead
multiple times). But this is not usually the rate-limiting step, so
it should not matter much. I hope this helps. (Please let me know if
I'm wrong.)
Cheers

Andrew

(Oops I meant to write
dihedral_coeff 5 -0.2 5 -450 0
dihedral_coeff 6 0.167 6 540 0

Not that it matters. (otherwise somebody here will correct me I bet)
Anyway you get the idea.)
cheers

Interesting - if yo uwant to write up a brief paragraph
or 2 for the approipriate dihedral style doc page, I'll include
it as a note to users.

Steve