- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Eric Shamay <eric.shamay.gmail.com>

Date: Fri, 11 Jun 2010 17:13:56 -0700

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

Date: Fri, 11 Jun 2010 17:13:56 -0700

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

Custom Search