Hi,

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 ...

Steve

Arnab,

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.

Paul

Paul,

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.

Thanks

Arnab

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.

Steve

Arnab,

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

http://www.ks.uiuc.edu/Research/namd/doxygen/ComputeDihedrals_8C-source.html

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

Paul

Arnab,

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.

Paul

#!/usr/bin/python

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-k*IQ) - IR*k

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

f.close()

random = Random(12345)

for i in xrange(100):

d = int(360*random())
n = int(6*random())

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

Masses

1 1

Pair Coeffs

1 0.0 0.0 0.0 0.0

Atoms

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

Dihedrals

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],

r4[0],r4[1],r4[2],K,n,d)

print >>f,txt

f.close()

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)]

subtract(r1,r2,r12)

subtract(r2,r3,r23)

subtract(r3,r4,r34)

cross(r12,r23,A)

rAinv = 1/sqrt(dot(A,A))

cross(r23,r34,B)

rBinv = 1/sqrt(dot(B,B))

cross(r23,A,C)

rCinv = 1/sqrt(dot(C,C))

cos_phi = dot(A,B)*(rAinv*rBinv)

sin_phi = dot(C,B)*(rCinv*rBinv)

phi = atan2(sin_phi,cos_phi)

energy = K*(1 + cos(n*phi - float(d*PI/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)

Hi,

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).

Thanks

Arnab

Paul Crozier wrote: