[lammps-users] Dihedral style

       I was looking into the code of dihedral_charmm.cpp and trying to figure out the calculation. A for loop starting from Line no 176 apparently looks missing something. Inside the loop as pasted below there seems no involvement of the variable 'i' and certainly seems erroneous. I appreciate any kind of input in this.

Thanks - we'll look into the code snippet. It does seem odd.
The multiharmonic formula was a bad link - I'll patch it in a minute ...



The variable i is just a looping index. The code looks strange, but seems to work fine. Have you found otherwise? Do you have a test case where it fails?

I did find one minor issue with the docs. This page: http://lammps.sandia.gov/doc/dihedral_style_charmm.html

should have

E = K[1 + cos(n*phi - d)]

instead of

E = K[1 + cos(n*phi + d)]

to agree with the code and with the cited MacKerell paper.


          I haven't tested the code. I am actually trying to extend the code in order to express the dihedral energy as a sum of two terms of the same form. I did miss the typo in the documentation though.


Paul Crozier wrote:

I just didn't look at the loop carefully enough. The loop is
over the multiplicity of the dihedral angle and it is accumulating
the quantity df1, resulting in a final value of p.



You might want to look at how they do this calc in NAMD:


Their source code for this calculation may be easier for you to follow.



I wrote a little python script to double check that LAMMPS is indeed computing the charmm dihedral energies correctly. I’ll paste it below and attach it. It picks random positions for 4 atoms, a random value for n, and a random value for d. Then it creates a LAMMPS input file and data file, and runs LAMMPS a single timestep to get the dihedral energy. Then it computes the dihedral energy on its own, in a way similar to what NAMD does, and compares it to what LAMMPS reported.

If you run it, you’ll see that the two very different approaches give answers that are essentially identical. This proves that LAMMPS is doing the calculation right, even if the code looks weird and more obfuscated than what NAMD has. The advantage is that the way LAMMPS does the calculation is probably a bit faster.



import commands, re
from math import acos, asin, atan, atan2, cos, sin, sqrt

PI = 4.0*atan(1.0)

v1 = [[] for i in xrange(3)]
v2 = [[] for i in xrange(3)]
v3 = [[] for i in xrange(3)]

def subtract(v1,v2,v3):
v3[0] = v2[0] - v1[0]
v3[1] = v2[1] - v1[1]
v3[2] = v2[2] - v1[2]

def dot(v1,v2):
sum = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
return sum

def cross(v1,v2,v3):
v3[0] = v1[1]*v2[2] - v1[2]*v2[1]
v3[1] = v1[2]*v2[0] - v1[0]*v2[2]
v3[2] = v1[0]*v2[1] - v1[1]*v2[0]

def norm(v1):
mag = sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2])
v1[0] /= mag
v1[1] /= mag
v1[2] /= mag


random number generator class

IM = 2147483647
AM = 1.0/IM
IA = 16807
IQ = 127773
IR = 2836

class Random:
def init(self,seed):
self.seed = seed

def call(self):
k = self.seed/IQ
self.seed = IA*(self.seed-kIQ) - IRk
if self.seed < 0: self.seed += IM
return AM*self.seed


def seed(new_seed):
random.seed = new_seed


f = open(“dih.in”,“w”)
text = ‘’‘units real
atom_style full
dihedral_style charmm
pair_style lj/charmm/coul/charmm 8 10
read_data dih.data
thermo_style custom edihed
run 0’’’
print >>f,text

random = Random(12345)

for i in xrange(100):

d = int(360random())
n = int(6
K = 1.0

r1 = [random(),random(),random()]
r2 = [random(),random(),random()]
r3 = [random(),random(),random()]
r4 = [random(),random(),random()]

f = open(“dih.data”,“w”)
text = ‘’’

4 atoms
0 bonds
0 angles
1 dihedrals
0 impropers

1 atom types
0 bond types
0 angle types
1 dihedral types
0 improper types

0 100 xlo xhi
0 100 ylo yhi
0 100 zlo zhi


1 1

Pair Coeffs

1 0.0 0.0 0.0 0.0


1 1 1 0.0 %f %f %f
2 1 1 0.0 %f %f %f
3 1 1 0.0 %f %f %f
4 1 1 0.0 %f %f %f

Dihedral Coeffs

1 %d %d %d 1


1 1 1 2 3 4

txt = text % (r1[0],r1[1],r1[2],r2[0],r2[1],r2[2],r3[0],r3[1],r3[2],

print >>f,txt

outtext = commands.getoutput(“lmp < dih.in”)

words = outtext.split()

r12 = [[] for i in xrange(3)]
r23 = [[] for i in xrange(3)]
r34 = [[] for i in xrange(3)]
A = [[] for i in xrange(3)]
B = [[] for i in xrange(3)]
C = [[] for i in xrange(3)]

rAinv = 1/sqrt(dot(A,A))
rBinv = 1/sqrt(dot(B,B))
rCinv = 1/sqrt(dot(C,C))

cos_phi = dot(A,B)(rAinvrBinv)
sin_phi = dot(C,B)(rCinvrBinv)
phi = atan2(sin_phi,cos_phi)

energy = K*(1 + cos(nphi - float(dPI/180.0)))

lammps_energy = float(words[80])
ediff = energy - lammps_energy
print "n = ", n, " d = ", d, " phi = ", phi, " energy = ", energy, ", error = ", ediff
if (abs(ediff) > 1e-4): print “Discrepancy found!”

dih.py (3.34 KB)

       Thanks for all the effort. Everything helped a lot and i believe i have understood the portion of the code i wanted to. As Steve pointed out the accumulation of 'p' and 'df1' value in the loop it basically means as i understand that the loop calculates the value of Cos (nPhi) and Sin(nPhi). I guess it is faster than directly calculating Cos (nPhi) and Sin(nPhi).


Paul Crozier wrote: