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