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-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
f.close()
random = Random(12345)
for i in xrange(100):
d = int(360random())
n = int(6random())
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)(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)
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: