1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: energy, enyinternal
|
---|
4 | !
|
---|
5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
6 | ! Shura Hayryan, Chin-Ku
|
---|
7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
8 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
9 | !
|
---|
10 | ! **************************************************************
|
---|
11 |
|
---|
12 | real*8 function energy()
|
---|
13 | ! ------------------------------------------
|
---|
14 | ! PURPOSE: calculate the *total* energy of the system
|
---|
15 | !
|
---|
16 | ! ientyp = 0 for ECEPP3, 1 for FLEX, 2 for Lund
|
---|
17 | ! 3 for ECEPP3 with Abagyan corrections
|
---|
18 | !
|
---|
19 | ! CALLS: enyflx,enyreg,enyshe,enysol,esolan,setvar, eninteract
|
---|
20 | !
|
---|
21 | ! -------------------------------------------------
|
---|
22 |
|
---|
23 | include 'INCL.H'
|
---|
24 | double precision teysm, teyel, teyvw, teyhb, teyvr
|
---|
25 | ! print *,'energy function with ientyp = ',ientyp
|
---|
26 | esm = 0.d0
|
---|
27 | teysm = 0.d0
|
---|
28 | teyel = 0.d0
|
---|
29 | teyvw = 0.d0
|
---|
30 | teyhb = 0.d0
|
---|
31 | teyvr = 0.d0
|
---|
32 | teysl = 0.d0
|
---|
33 |
|
---|
34 | do i = 1,ntlml
|
---|
35 | eysm=0
|
---|
36 | eyel=0
|
---|
37 | eyvr=0
|
---|
38 | eyhb=0
|
---|
39 | eyvw=0
|
---|
40 | eysl=0
|
---|
41 |
|
---|
42 | call setvar(i,vlvr) ! set variables & rebuild
|
---|
43 |
|
---|
44 | if (ientyp.eq.0.or.ientyp.eq.3) then
|
---|
45 | esm=esm+enyshe(i)
|
---|
46 | else if (ientyp.eq.1) then
|
---|
47 | esm=esm+enyflx(i)
|
---|
48 | else if (ientyp.eq.2) then
|
---|
49 | esm=enylun(i)
|
---|
50 | endif
|
---|
51 |
|
---|
52 | teysm = teysm + eysm
|
---|
53 | teyel = teyel + eyel
|
---|
54 | teyhb = teyhb + eyhb
|
---|
55 | teyvr = teyvr + eyvr
|
---|
56 | if (ientyp.eq.2) then
|
---|
57 | ! The Lund term stores the hydrophobicity energy in eysl
|
---|
58 | teysl = teysl + eysl
|
---|
59 | else
|
---|
60 | ! .. and the excluded volume term in eyvw, which is calculated once.
|
---|
61 | teyvw = teyvw + eyvw
|
---|
62 | endif
|
---|
63 |
|
---|
64 | if (ireg.eq.1) eyrg=enyreg(i)
|
---|
65 |
|
---|
66 | enddo
|
---|
67 |
|
---|
68 | if (ientyp.ne.2) then
|
---|
69 | ! Don't touch eysl if using Lund potential, as enylun stores
|
---|
70 | ! its hydrophobicity term there.
|
---|
71 | if (itysol.gt.0) then
|
---|
72 | esm=esm+enysol(0)
|
---|
73 | teysl = teysl+eysl
|
---|
74 | else if (itysol.lt.0) then
|
---|
75 | eysl = esolan(0)
|
---|
76 | teysl = teysl + eysl
|
---|
77 | esm = esm + eysl
|
---|
78 | else
|
---|
79 | eysl=0.d0
|
---|
80 | endif
|
---|
81 | else
|
---|
82 | ! Add excluded volume term and save it in eyvw
|
---|
83 | esm=esm+exvlun(0)
|
---|
84 | teyvw = teyvw+eyvw
|
---|
85 | endif
|
---|
86 |
|
---|
87 | ! The Abagyan entropic corrections depend on the area exposed to the
|
---|
88 | ! solvent for each residue. So, this term has to be evaluated after the
|
---|
89 | ! solvent term.
|
---|
90 | eyab=0.0
|
---|
91 | if (ientyp.eq.3) then
|
---|
92 | do i = 1,ntlml
|
---|
93 | eyab=eyab+eyabgn(i)
|
---|
94 | enddo
|
---|
95 | endif
|
---|
96 | esm=esm+eyab
|
---|
97 | ! Partial energies for the entire system. If you need the partial
|
---|
98 | ! energies for a single molecule call enyinternal.
|
---|
99 | eysm = teysm
|
---|
100 | eyel = teyel
|
---|
101 | eyvw = teyvw
|
---|
102 | eyhb = teyhb
|
---|
103 | eyvr = teyvr
|
---|
104 | eysl = teysl
|
---|
105 |
|
---|
106 | if (ientyp.ne.2) then
|
---|
107 | ! This is temporary. eninteract() does not yet know how to calculate
|
---|
108 | ! interactions using the Lund potential.
|
---|
109 | energy = esm + eninteract()
|
---|
110 | return
|
---|
111 | endif
|
---|
112 | energy=esm
|
---|
113 | return
|
---|
114 | end
|
---|
115 |
|
---|
116 | !c Calculates the internal energy for a single molecule.
|
---|
117 | ! All the partial energies are thus set to their values for molecule
|
---|
118 | ! nml.
|
---|
119 | !
|
---|
120 | ! @param nml the ID of the molecule
|
---|
121 | ! @return internal energy of a single molecule
|
---|
122 | !
|
---|
123 | ! @author Jan H. Meinke <j.meinke@fz-juelich.de>
|
---|
124 | real*8 function enyinternal(nml)
|
---|
125 |
|
---|
126 | !f2py intent(in) nml
|
---|
127 |
|
---|
128 | include 'INCL.H'
|
---|
129 | esm = 0.d0
|
---|
130 |
|
---|
131 | call setvar(nml,vlvr)
|
---|
132 |
|
---|
133 | if (ientyp.eq.0.or.ientyp.eq.3) then
|
---|
134 | esm=esm+enyshe(nml)
|
---|
135 | else if (ientyp.eq.1) then
|
---|
136 | esm=esm+enyflx(nml)
|
---|
137 | else if (ientyp.eq.2) then
|
---|
138 | esm=esm+enylun(nml)
|
---|
139 | endif
|
---|
140 |
|
---|
141 | if (ireg.eq.1) eyrg=enyreg(nml)
|
---|
142 |
|
---|
143 | if (ientyp.ne.2) then
|
---|
144 | if (itysol.gt.0) then
|
---|
145 | esm=esm+enysol(nml)
|
---|
146 | elseif (itysol.lt.0) then
|
---|
147 | esm=esm+esolan(nml)
|
---|
148 | else
|
---|
149 | eysl=0.d0
|
---|
150 | endif
|
---|
151 | else
|
---|
152 | esm=esm+exvlun(nml)
|
---|
153 | endif
|
---|
154 | ! The Abagyan entropic corrections depend on the area exposed to the
|
---|
155 | ! solvent for each residue. So, this term has to be evaluated after the
|
---|
156 | ! solvent term.
|
---|
157 | eyab=0.0
|
---|
158 | if (ientyp.eq.3) then
|
---|
159 | do i=1,ntlml
|
---|
160 | eyab=eyab + eyabgn(i)
|
---|
161 | enddo
|
---|
162 | endif
|
---|
163 | esm=esm+eyab
|
---|
164 |
|
---|
165 | enyinternal = esm
|
---|
166 | return
|
---|
167 | end
|
---|