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