I would like to ask some question regarding the source code of Sander module. I hope somebody having experience in the development of the package could give an answer.
I understand that the dihedral angle energy is calculated in the subroutine ‘ephi’ in the file ‘ene.f’ under the ‘src/sander’ directory (Amber9). The 1-4 interactions are also calculated in that subroutine. My question is related to the part of the calculation of the 1-4 interactions. To calculate the 1-4 interaction, one should first calculate the distance between 1^st and 4^th atoms. I show the relevant code snippet from the “subroutine ephi” adding the line number counted from the first line of ‘ene.f’:
793 do jn = 1,maxlen
794 i3 = ip(jn+ist)
795 k3t = kp(jn+ist)
796 l3t = lp(jn+ist)
797 ic0 = icp(jn+ist)
798 idumi = sign(1,k3t)
799 iduml = sign(1,l3t)
800 kdiv = (2+idumi+iduml)/4
801 l3 = iabs(l3t)
802 fmuln = kdiv*fmn(ic0)
803
804 ii = (i3+3)/3
805 jj = (l3+3)/3
806 ia1 = iac(ii)
807 ia2 = iac(jj)
808 ibig = max0(ia1,ia2)
809 isml = min0(ia1,ia2)
810 ic = ibig*(ibig-1)/2+isml
811
812 ! ----- CALCULATE THE 14-EEL ENERGY -----
813
814 r2 = fmuln/ct(jn)
815 r1 = sqrt(r2)
“ct(jn)” at line 814 is the square of the 1-4 distance. The inverse of the square of the 1-4 distance is calculated at line 814 and it has a factor given as ‘fmuln’. ‘fmuln’ is calculated at line 802, where “kdiv” becomes zero if the 3^rd and 4^th atom indexes are negative. My understanding is that, when the four atoms specifying a dihedral angle are not sequentially connected, the 3^rd and 4^th atom indexes are negative in the prmtop file, thus making “kdiv” zero.
But there are cases where only 3th atom index is negative in the DIHEDRALS_WITHOUT_HYDROGEN section and DIHEDRALS_INC_HYDROGEN section of the prmtop file. I understand that these are cases where the dihedral potential term consists of more than one cosine functions. In these cases, the same four atom indexes are repeated as many as the number of the terms in the dihedral potential expression and from the second set the 3rd atom index is negative. Is this understanding correct?
In this situation I think that, when only 3^rd atom index is negative in the prmtop file, the value of the variable ‘fmn(ic0)’ should be zero in the line 802 so that fmuln be zero. (because “kdiv” in line 800 is not zero but 1/2) So my main question is how this happens. I found that the values of the variable ‘fmn’ are set at the subroutine “dihpar” in the file ‘set.f’. In the subroutine “dihpar”, the value of “fmn(i)” is set to zero only when “pn(i)” is negative. I understand that the variable “pn(i)” is the periodicity of the cosine function used to express the dihedral angle energy term. In the prmtop file I have used in the normal simulations of peptides, I couldn’t find any case when the value of “pn” parameter is negative. It seems to me that the values of “pn” are always positive in the prmtop file made from the LEAP program. So, what I would like to know is when and how the value of the variable “fmn” become zero. If !
there is misunderstanding in my part, please let me know.
------------------------------------------------------------------------
NAVER :: Korea's No.1 search portal
www.naver.com
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Oct 22 2006 - 06:07:27 PDT
