I have a need to perform a modified calculation for water bond energies. I
plan to model systems with multiple residue types of which one will be a
water (flexible water as per Ferguson 1995). In that model the bond energy
includes a cubic term: Kl * Kcub * (l-l0)^3.
The problem seems simple enough on paper, so I decided to dive into the code
and track down where the bond calculation happens. I'm not very Fortran
savvy (all my coding is in C++, python, bash, etc.) and so I'm seeking help
in this endeavor.
Thus far I have found what looks like the proper calculation for the bond
energy occurs in the file: src/sander/ene.f
The relevant code is as follows (this is Amber 9.0):
108 do jn = 1,maxlen
109 rij0 = xij(jn)*xij(jn)+yij(jn)*yij(jn)+zij(jn)*zij(jn)
110 rij(jn) = sqrt(rij0)
111 end do
112
113 ! ----- CALCULATION OF THE ENERGY AND DER -----
114
115 do jn = 1,maxlen
116 ic = icb(jn+ist)
117 rij0 = rij(jn)
118 da = rij0-req(ic)
119 ! for rms deviation from ideal
bonds:
120 ebdev = ebdev + da*da
121 df = rk(ic)*da
122 eaw(jn) = df*da
123 if(idecomp == 1 .or. idecomp == 2) then
124 ii = (ib(jn+ist) + 3)/3
125 jj = (jb(jn+ist) + 3)/3
126 call decpair(4,ii,jj,eaw(jn))
127 end if
108-111: calculate the bond length.
118: find the difference between the actual bond length and the equilibrium
length
121: Calculate... df? It's looks like the bond constant x the bond length
difference. (noted later below)
122: eaw appears to hold the bond energy term. eaw = the above energy term,
namely ' bond constant x (R - Req)^2 '
123-127: do bond decomposition (I'm not interested in this right now)
Okay, good. So the eaw variable is our bond energy term. My gut reaction is
that I can simply add to eaw another term such as:
! --- Add in a cubic bond term ---
eaw(jn) = eaw(jn) + rc(ic)*df*da*da
where rc(ic) would correspond to a new cubic-bond-constant that I'd have to
read in from somewhere (no clue where variables are read in right now). This
new added term is equal to the cubic term but with another cubic-constant
multiplier.
There are still a few issues here:
1) Where else do I need to modify behavior? It looks like further down in
ene.f:
155 dfw(jn) = (df+df)/rij0
156 end do
157
158 ! ----- CALCULATION OF THE FORCE -----
159
160 do jn = 1,maxlen
161 i3 = ib(jn+ist)
162 j3 = jb(jn+ist)
163 df = dfw(jn)
164 xa = df*xij(jn)
165 ya = df*yij(jn)
166 za = df*zij(jn)
167 f(i3+1) = f(i3+1)-xa
168 f(i3+2) = f(i3+2)-ya
169 f(i3+3) = f(i3+3)-za
170 f(j3+1) = f(j3+1)+xa
171 f(j3+2) = f(j3+2)+ya
172 f(j3+3) = f(j3+3)+za
155: We again use df in calculating the force on the bond (?) which
translates into a force on the atoms. How should the added cubic term above
affect the df value to keep things sane?
2) What needs to be modified to allow for a cubic bond constant term to be
supplied. Ideally this would be in the frcmod file. We could define a BOND
section such that we could do this:
BOND
OW-HW 320.0 1.000 50.0 SPC/E and POL3 water
this way we could keep the previous format, but add in a third field that is
read as the cubic term.
3) How do I specify that *ONLY* certain bond types are to be handled with
the cubic term, while leaving all the others alone? I only want water O-H
bonds processed, and all my other residues to skip this calculation or
modification to the bond energy.
I hope someone out there has hacked the Amber code enough to be familiar
with these types of modifications. They seem simple enough to me in theory,
I'm hoping that someone will chime in with some advice.
Thanks!
~Eric
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Jun 11 2010 - 17:30:03 PDT
