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