1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ###############################################################
9c     ##                                                           ##
10c     ##  subroutine analysis  --  energy components and analysis  ##
11c     ##                                                           ##
12c     ###############################################################
13c
14c
15c     "analysis" calls the series of routines needed to calculate
16c     the potential energy and perform energy partitioning analysis
17c     in terms of type of interaction or atom number
18c
19c
20      subroutine analysis (energy)
21      use analyz
22      use atoms
23      use energi
24      use group
25      use inter
26      use iounit
27      use limits
28      use potent
29      use vdwpot
30      implicit none
31      integer i
32      real*8 energy
33      real*8 cutoff
34c
35c
36c     zero out each of the potential energy components
37c
38      eb = 0.0d0
39      ea = 0.0d0
40      eba = 0.0d0
41      eub = 0.0d0
42      eaa = 0.0d0
43      eopb = 0.0d0
44      eopd = 0.0d0
45      eid = 0.0d0
46      eit = 0.0d0
47      et = 0.0d0
48      ept = 0.0d0
49      ebt = 0.0d0
50      eat = 0.0d0
51      ett = 0.0d0
52      ev = 0.0d0
53      er = 0.0d0
54      edsp = 0.0d0
55      ec = 0.0d0
56      ecd = 0.0d0
57      ed = 0.0d0
58      em = 0.0d0
59      ep = 0.0d0
60      ect = 0.0d0
61      erxf = 0.0d0
62      es = 0.0d0
63      elf = 0.0d0
64      eg = 0.0d0
65      ex = 0.0d0
66c
67c     perform dynamic allocation of some global arrays
68c
69      if (allocated(aesum)) then
70         if (size(aesum) .lt. n) then
71            deallocate (aesum)
72            deallocate (aeb)
73            deallocate (aea)
74            deallocate (aeba)
75            deallocate (aeub)
76            deallocate (aeaa)
77            deallocate (aeopb)
78            deallocate (aeopd)
79            deallocate (aeid)
80            deallocate (aeit)
81            deallocate (aet)
82            deallocate (aept)
83            deallocate (aebt)
84            deallocate (aeat)
85            deallocate (aett)
86            deallocate (aev)
87            deallocate (aer)
88            deallocate (aedsp)
89            deallocate (aec)
90            deallocate (aecd)
91            deallocate (aed)
92            deallocate (aem)
93            deallocate (aep)
94            deallocate (aect)
95            deallocate (aerxf)
96            deallocate (aes)
97            deallocate (aelf)
98            deallocate (aeg)
99            deallocate (aex)
100         end if
101      end if
102      if (.not. allocated(aesum)) then
103         allocate (aesum(n))
104         allocate (aeb(n))
105         allocate (aea(n))
106         allocate (aeba(n))
107         allocate (aeub(n))
108         allocate (aeaa(n))
109         allocate (aeopb(n))
110         allocate (aeopd(n))
111         allocate (aeid(n))
112         allocate (aeit(n))
113         allocate (aet(n))
114         allocate (aept(n))
115         allocate (aebt(n))
116         allocate (aeat(n))
117         allocate (aett(n))
118         allocate (aev(n))
119         allocate (aer(n))
120         allocate (aedsp(n))
121         allocate (aec(n))
122         allocate (aecd(n))
123         allocate (aed(n))
124         allocate (aem(n))
125         allocate (aep(n))
126         allocate (aect(n))
127         allocate (aerxf(n))
128         allocate (aes(n))
129         allocate (aelf(n))
130         allocate (aeg(n))
131         allocate (aex(n))
132      end if
133c
134c     zero out energy partitioning components for each atom
135c
136      do i = 1, n
137         aeb(i) = 0.0d0
138         aea(i) = 0.0d0
139         aeba(i) = 0.0d0
140         aeub(i) = 0.0d0
141         aeaa(i) = 0.0d0
142         aeopb(i) = 0.0d0
143         aeopd(i) = 0.0d0
144         aeid(i) = 0.0d0
145         aeit(i) = 0.0d0
146         aet(i) = 0.0d0
147         aept(i) = 0.0d0
148         aebt(i) = 0.0d0
149         aeat(i) = 0.0d0
150         aett(i) = 0.0d0
151         aev(i) = 0.0d0
152         aer(i) = 0.0d0
153         aedsp(i) = 0.0d0
154         aec(i) = 0.0d0
155         aecd(i) = 0.0d0
156         aed(i) = 0.0d0
157         aem(i) = 0.0d0
158         aep(i) = 0.0d0
159         aect(i) = 0.0d0
160         aerxf(i) = 0.0d0
161         aes(i) = 0.0d0
162         aelf(i) = 0.0d0
163         aeg(i) = 0.0d0
164         aex(i) = 0.0d0
165      end do
166c
167c     zero out the total intermolecular energy
168c
169      einter = 0.0d0
170c
171c     remove any previous use of the replicates method
172c
173      cutoff = 0.0d0
174      call replica (cutoff)
175c
176c     update the pairwise interaction neighbor lists
177c
178      if (use_list)  call nblist
179c
180c     many implicit solvation models require Born radii
181c
182      if (use_born)  call born
183c
184c     alter partial charges and multipoles for charge flux
185c
186      if (use_chgflx)  call alterchg
187c
188c     modify bond and torsion constants for pisystem
189c
190      if (use_orbit)  call picalc
191c
192c     call the local geometry energy component routines
193c
194      if (use_bond)  call ebond3
195      if (use_angle)  call eangle3
196      if (use_strbnd)  call estrbnd3
197      if (use_urey)  call eurey3
198      if (use_angang)  call eangang3
199      if (use_opbend)  call eopbend3
200      if (use_opdist)  call eopdist3
201      if (use_improp)  call eimprop3
202      if (use_imptor)  call eimptor3
203      if (use_tors)  call etors3
204      if (use_pitors)  call epitors3
205      if (use_strtor)  call estrtor3
206      if (use_angtor)  call eangtor3
207      if (use_tortor)  call etortor3
208c
209c     call the electrostatic energy component routines
210c
211      if (use_charge)  call echarge3
212      if (use_chgdpl)  call echgdpl3
213      if (use_dipole)  call edipole3
214      if (use_mpole)  call empole3
215      if (use_polar)  call epolar3
216      if (use_chgtrn)  call echgtrn3
217      if (use_rxnfld)  call erxnfld3
218c
219c     call the van der Waals energy component routines
220c
221      if (use_vdw) then
222         if (vdwtyp .eq. 'LENNARD-JONES')  call elj3
223         if (vdwtyp .eq. 'BUCKINGHAM')  call ebuck3
224         if (vdwtyp .eq. 'MM3-HBOND')  call emm3hb3
225         if (vdwtyp .eq. 'BUFFERED-14-7')  call ehal3
226         if (vdwtyp .eq. 'GAUSSIAN')  call egauss3
227      end if
228      if (use_repuls)  call erepel3
229      if (use_disp)  call edisp3
230c
231c     call any miscellaneous energy component routines
232c
233      if (use_solv)  call esolv3
234      if (use_metal)  call emetal3
235      if (use_geom)  call egeom3
236      if (use_extra)  call extra3
237c
238c     sum up to give the total potential energy
239c
240      esum = eb + ea + eba + eub + eaa + eopb + eopd + eid + eit
241     &          + et + ept + ebt + eat + ett + ev + er + edsp
242     &          + ec + ecd + ed + em + ep + ect + erxf + es + elf
243     &          + eg + ex
244      energy = esum
245c
246c     sum up to give the total potential energy per atom
247c
248      do i = 1, n
249         aesum(i) = aeb(i) + aea(i) + aeba(i) + aeub(i) + aeaa(i)
250     &                 + aeopb(i) + aeopd(i) + aeid(i) + aeit(i)
251     &                 + aet(i) + aept(i) + aebt(i) + aeat(i) + aett(i)
252     &                 + aev(i) + aer(i) + aedsp(i) + aec(i) + aecd(i)
253     &                 + aed(i) + aem(i) + aep(i) + aect(i) + aerxf(i)
254     &                 + aes(i) + aelf(i) + aeg(i) + aex(i)
255      end do
256c
257c     check for an illegal value for the total energy
258c
259c     if (isnan(esum)) then
260      if (esum .ne. esum) then
261         write (iout,10)
262   10    format (/,' ANALYSIS  --  Illegal Value for the Total',
263     &              ' Potential Energy')
264         call fatal
265      end if
266      return
267      end
268