| 1 | ! *******************************************************************
|
|---|
| 2 | ! SMMP version of Anders Irback's force field, to be called the Lund
|
|---|
| 3 | ! force field. This file contains the function enylun, which in turn
|
|---|
| 4 | ! calls all the terms in the energy function. The terms Bias (ebias),
|
|---|
| 5 | ! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the
|
|---|
| 6 | ! Excluded volume (eexvol and eloexv) are also implemented in this
|
|---|
| 7 | ! file.
|
|---|
| 8 | !
|
|---|
| 9 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 10 | ! Jan H. Meinke, Sandipan Mohanty
|
|---|
| 11 | !
|
|---|
| 12 | subroutine init_lundff
|
|---|
| 13 | include 'INCL.H'
|
|---|
| 14 | include 'incl_lund.h'
|
|---|
| 15 | character mynm*4
|
|---|
| 16 | logical prlvr
|
|---|
| 17 |
|
|---|
| 18 | print *,'initializing Lund forcefield'
|
|---|
| 19 | ! Some parameters in the Lund force field.
|
|---|
| 20 | ! The correspondence between internal energy scale and kcal/mol
|
|---|
| 21 | eunit=1.3315
|
|---|
| 22 | ! Bias
|
|---|
| 23 | kbias=100.0*eunit
|
|---|
| 24 | ! print *,'Bias'
|
|---|
| 25 | ! Hydrogen bonds
|
|---|
| 26 | epshb1=3.1*eunit
|
|---|
| 27 | epshb2=2.0*eunit
|
|---|
| 28 | sighb=2.0
|
|---|
| 29 | cthb=4.5
|
|---|
| 30 | cthb2=cthb*cthb
|
|---|
| 31 | powa=0.5
|
|---|
| 32 | powb=0.5
|
|---|
| 33 | blhb=-30.0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
|
|---|
| 34 | alhb=-(5*((sighb/cthb)**12)-6*((sighb/cthb)**10))-blhb*cthb2
|
|---|
| 35 | sighb2=sighb*sighb
|
|---|
| 36 | cdon=1.0
|
|---|
| 37 | cacc=(1.0/1.23)**powb
|
|---|
| 38 | csacc=(1.0/1.25)**powb
|
|---|
| 39 | ! print *,'Hydrogen bonds'
|
|---|
| 40 | ! Hydrophobicity
|
|---|
| 41 | ! print *,'Hydrophobicity with nhptyp = ',nhptyp
|
|---|
| 42 |
|
|---|
| 43 | hpstrg(1)=0.0*eunit
|
|---|
| 44 | hpstrg(2)=0.1*eunit
|
|---|
| 45 | hpstrg(3)=0.1*eunit
|
|---|
| 46 | hpstrg(4)=0.1*eunit
|
|---|
| 47 | hpstrg(5)=0.9*eunit
|
|---|
| 48 | hpstrg(6)=2.8*eunit
|
|---|
| 49 | hpstrg(7)=0.1*eunit
|
|---|
| 50 | hpstrg(8)=2.8*eunit
|
|---|
| 51 | hpstrg(9)=3.2*eunit
|
|---|
| 52 |
|
|---|
| 53 | do i=1,mxrs
|
|---|
| 54 | do j=1,6
|
|---|
| 55 | ihpat(i,j)=0
|
|---|
| 56 | enddo
|
|---|
| 57 | nhpat(i)=0
|
|---|
| 58 | enddo
|
|---|
| 59 | do i=irsml1(1),irsml2(ntlml)
|
|---|
| 60 | mynm=seq(i)
|
|---|
| 61 | call tolost(mynm)
|
|---|
| 62 | if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
|
|---|
| 63 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
|---|
| 64 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
|---|
| 65 | prlvr=.true. ! residue i is a proline variant
|
|---|
| 66 | else
|
|---|
| 67 | prlvr=.false.
|
|---|
| 68 | endif
|
|---|
| 69 |
|
|---|
| 70 | if (mynm.eq.'ala') then
|
|---|
| 71 | nhpat(i)=1
|
|---|
| 72 | ihpat(i,1)=iCa(i)+2
|
|---|
| 73 | else if (mynm.eq.'val') then
|
|---|
| 74 | nhpat(i)=3
|
|---|
| 75 | ihpat(i,1)=iCa(i)+2
|
|---|
| 76 | ihpat(i,2)=iCa(i)+4
|
|---|
| 77 | ihpat(i,3)=iCa(i)+8
|
|---|
| 78 | else if (prlvr) then
|
|---|
| 79 | nhpat(i)=3
|
|---|
| 80 | ihpat(i,1)=iCa(i)+2
|
|---|
| 81 | ihpat(i,2)=iCa(i)+5
|
|---|
| 82 | ihpat(i,3)=iCa(i)+8
|
|---|
| 83 | else if (mynm.eq.'leu') then
|
|---|
| 84 | nhpat(i)=4
|
|---|
| 85 | ihpat(i,1)=iCa(i)+2
|
|---|
| 86 | ihpat(i,2)=iCa(i)+5
|
|---|
| 87 | ihpat(i,3)=iCa(i)+7
|
|---|
| 88 | ihpat(i,3)=iCa(i)+11
|
|---|
| 89 | else if (mynm.eq.'ile') then
|
|---|
| 90 | nhpat(i)=4
|
|---|
| 91 | ihpat(i,1)=iCa(i)+2
|
|---|
| 92 | ihpat(i,2)=iCa(i)+4
|
|---|
| 93 | ihpat(i,3)=iCa(i)+8
|
|---|
| 94 | ihpat(i,3)=iCa(i)+11
|
|---|
| 95 | else if (mynm.eq.'met') then
|
|---|
| 96 | nhpat(i)=4
|
|---|
| 97 | ihpat(i,1)=iCa(i)+2
|
|---|
| 98 | ihpat(i,2)=iCa(i)+5
|
|---|
| 99 | ihpat(i,3)=iCa(i)+8
|
|---|
| 100 | ihpat(i,3)=iCa(i)+9
|
|---|
| 101 | else if (mynm.eq.'phe') then
|
|---|
| 102 | nhpat(i)=6
|
|---|
| 103 | ihpat(i,1)=iCa(i)+5
|
|---|
| 104 | ihpat(i,2)=iCa(i)+6
|
|---|
| 105 | ihpat(i,3)=iCa(i)+8
|
|---|
| 106 | ihpat(i,3)=iCa(i)+10
|
|---|
| 107 | ihpat(i,3)=iCa(i)+12
|
|---|
| 108 | ihpat(i,3)=iCa(i)+14
|
|---|
| 109 | else if (mynm.eq.'tyr') then
|
|---|
| 110 | nhpat(i)=6
|
|---|
| 111 | ihpat(i,1)=iCa(i)+5
|
|---|
| 112 | ihpat(i,2)=iCa(i)+6
|
|---|
| 113 | ihpat(i,3)=iCa(i)+8
|
|---|
| 114 | ihpat(i,3)=iCa(i)+10
|
|---|
| 115 | ihpat(i,3)=iCa(i)+13
|
|---|
| 116 | ihpat(i,3)=iCa(i)+15
|
|---|
| 117 | else if (mynm.eq.'trp') then
|
|---|
| 118 | nhpat(i)=6
|
|---|
| 119 | ihpat(i,1)=iCa(i)+10
|
|---|
| 120 | ihpat(i,2)=iCa(i)+11
|
|---|
| 121 | ihpat(i,3)=iCa(i)+13
|
|---|
| 122 | ihpat(i,3)=iCa(i)+15
|
|---|
| 123 | ihpat(i,3)=iCa(i)+17
|
|---|
| 124 | ihpat(i,3)=iCa(i)+19
|
|---|
| 125 | endif
|
|---|
| 126 | enddo
|
|---|
| 127 | ! print *,'Hydrophobicity'
|
|---|
| 128 |
|
|---|
| 129 | ! Excluded volume and local pair excluded volume terms
|
|---|
| 130 | exvk=0.1*eunit
|
|---|
| 131 | exvcut=4.3
|
|---|
| 132 | exvcut2=exvcut*exvcut
|
|---|
| 133 |
|
|---|
| 134 | sigsa(1)=1.0 ! hydrogen
|
|---|
| 135 | sigsa(2)=1.0 ! hydrogen
|
|---|
| 136 | sigsa(3)=1.0 ! hydrogen
|
|---|
| 137 | sigsa(4)=1.0 ! hydrogen
|
|---|
| 138 | sigsa(5)=1.0 ! hydrogen
|
|---|
| 139 | sigsa(6)=1.0 ! hydrogen
|
|---|
| 140 | sigsa(7)=1.75 ! carbon
|
|---|
| 141 | sigsa(8)=1.75 ! carbon
|
|---|
| 142 | sigsa(9)=1.75 ! carbon
|
|---|
| 143 | sigsa(10)=1.42 ! oxygen
|
|---|
| 144 | sigsa(11)=1.42 ! oxygen
|
|---|
| 145 | sigsa(12)=1.42 ! oxygen
|
|---|
| 146 | sigsa(13)=1.55 ! nitrogen
|
|---|
| 147 | sigsa(14)=1.55 ! nitrogen
|
|---|
| 148 | sigsa(15)=1.55 ! nitrogen
|
|---|
| 149 | sigsa(16)=1.77 ! sulfur
|
|---|
| 150 | sigsa(17)=1.0 ! hydrogen
|
|---|
| 151 | sigsa(18)=1.75 ! carbon
|
|---|
| 152 |
|
|---|
| 153 | do i=1,mxtyat
|
|---|
| 154 | do j=1,mxtyat
|
|---|
| 155 | sig2lcp(i,j)=(sigsa(i)+sigsa(j))**2
|
|---|
| 156 | asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0)
|
|---|
| 157 | bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0)/exvcut2
|
|---|
| 158 | enddo
|
|---|
| 159 | enddo
|
|---|
| 160 | ! print *,'Local pair excluded volume constants'
|
|---|
| 161 |
|
|---|
| 162 | exvlam=0.75
|
|---|
| 163 | exvcutg=exvcut*exvlam
|
|---|
| 164 | exvcutg2=exvcutg*exvcutg
|
|---|
| 165 |
|
|---|
| 166 | do i=1,mxtyat
|
|---|
| 167 | do j=1,mxtyat
|
|---|
| 168 | sig2exv(i,j)=(exvlam*(sigsa(i)+sigsa(j)))**2
|
|---|
| 169 | asaexv(i,j)=-7*((sig2exv(i,j)/exvcutg2)**6.0)
|
|---|
| 170 | bsaexv(i,j)=6*((sig2exv(i,j)/exvcutg2)**6.0)/exvcutg2
|
|---|
| 171 | enddo
|
|---|
| 172 | enddo
|
|---|
| 173 | ! print *,'General excluded volume constants'
|
|---|
| 174 |
|
|---|
| 175 | ! Initialization of the connections matrix matcon(i,j). The index
|
|---|
| 176 | ! i runs from -mxconr to +mxconr, and j from 1 to mxat.
|
|---|
| 177 | ! matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed
|
|---|
| 178 | ! = 2, if atoms i1 and i2 are separated by 3 covalent
|
|---|
| 179 | ! bonds and their distance can change
|
|---|
| 180 | ! = 1, for all other pairs
|
|---|
| 181 | ! if abs(i2-i1) > mxconr, the atoms are assumed to be separated by
|
|---|
| 182 | ! many bonds, and with no restriction on their distances. On a protein
|
|---|
| 183 | ! molecule made of natural amino acids, atoms with indices separated
|
|---|
| 184 | ! by more than 35 can not be connected by three covalent bonds.
|
|---|
| 185 |
|
|---|
| 186 | do i=1,mxat
|
|---|
| 187 | do j=-mxconr,mxconr
|
|---|
| 188 | matcon(j,i)=1
|
|---|
| 189 | enddo
|
|---|
| 190 | matcon(0,i)=0
|
|---|
| 191 | enddo
|
|---|
| 192 | ! continued...
|
|---|
| 193 | do iml=1,ntlml
|
|---|
| 194 | do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
|
|---|
| 195 | do j=1,nbdat(iat1)
|
|---|
| 196 | iat2=ibdat(j,iat1)
|
|---|
| 197 | matcon(iat2-iat1,iat1)=0 ! directly bonded atoms
|
|---|
| 198 | matcon(iat1-iat2,iat2)=0 !
|
|---|
| 199 | do k=1,nbdat(iat2)
|
|---|
| 200 | iat3=ibdat(k,iat2)
|
|---|
| 201 | matcon(iat3-iat1,iat1)=0 !
|
|---|
| 202 | matcon(iat1-iat3,iat3)=0 ! 2 covalent bonds
|
|---|
| 203 | do l=1,nbdat(iat3)
|
|---|
| 204 | iat4=ibdat(l,iat3)
|
|---|
| 205 | if (iat4.ne.iat2) then
|
|---|
| 206 | matcon(iat4-iat1,iat1)=2 !
|
|---|
| 207 | matcon(iat1-iat4,iat4)=2 ! 3 covalent bonds
|
|---|
| 208 | endif
|
|---|
| 209 | do m=1,nbdat(iat2)
|
|---|
| 210 | iat3p=ibdat(m,iat2)
|
|---|
| 211 | if (iat3p.ne.iat3) then
|
|---|
| 212 | matcon(iat4-iat3p,iat3p)=2 ! 3 covalent bonds
|
|---|
| 213 | matcon(iat3p-iat4,iat4)=2 !
|
|---|
| 214 | endif
|
|---|
| 215 | enddo
|
|---|
| 216 | enddo
|
|---|
| 217 | enddo
|
|---|
| 218 | do k=1,nbdat(iat1)
|
|---|
| 219 | iat3=ibdat(k,iat1)
|
|---|
| 220 | matcon(iat3-iat2,iat2)=0 ! also 2 bonds iat2-iat1-iat3
|
|---|
| 221 | matcon(iat2-iat3,iat3)=0 !
|
|---|
| 222 | enddo
|
|---|
| 223 | enddo
|
|---|
| 224 | enddo
|
|---|
| 225 |
|
|---|
| 226 | ! print *,'going to initialize connections for first residue'
|
|---|
| 227 | ! print *,'iN,iCa,iC =',iN(irsml1(iml)),
|
|---|
| 228 | ! # iCa(irsml1(iml)),iC(irsml1(iml))
|
|---|
| 229 | do iat1=iN(irsml1(iml))+1,iCa(irsml1(iml))-1
|
|---|
| 230 | ! print *,'connections for iat1 = ',iat1
|
|---|
| 231 | matcon(iat1-iN(irsml1(iml)),iN(irsml1(iml)))=0
|
|---|
| 232 | matcon(iN(irsml1(iml))-iat1,iat1)=0
|
|---|
| 233 | matcon(iat1-iCa(irsml1(iml)),iCa(irsml1(iml)))=0
|
|---|
| 234 | matcon(iCa(irsml1(iml))-iat1,iat1)=0
|
|---|
| 235 |
|
|---|
| 236 | matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
|
|---|
| 237 | matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
|
|---|
| 238 | matcon(iat1-iCa(irsml1(iml))-2,iCa(irsml1(iml))+2)=2
|
|---|
| 239 | matcon(iCa(irsml1(iml))+2-iat1,iat1)=2
|
|---|
| 240 | matcon(iat1-iC(irsml1(iml)),iC(irsml1(iml)))=2
|
|---|
| 241 | matcon(iC(irsml1(iml))-iat1,iat1)=2
|
|---|
| 242 | enddo
|
|---|
| 243 |
|
|---|
| 244 | ! Below: for certain residues, some atoms separated by 3 or more bonds
|
|---|
| 245 | ! do not change distance. So, the connection matrix term for such pairs
|
|---|
| 246 | ! should be zero.
|
|---|
| 247 |
|
|---|
| 248 | do irs=irsml1(iml),irsml2(iml)
|
|---|
| 249 | if (irs.eq.irsml1(iml)) then
|
|---|
| 250 | iatoff=1
|
|---|
| 251 | else
|
|---|
| 252 | iatoff=0
|
|---|
| 253 | endif
|
|---|
| 254 | if (irs.eq.irsml2(iml)) then
|
|---|
| 255 | iatmrg=2
|
|---|
| 256 | else
|
|---|
| 257 | iatmrg=0
|
|---|
| 258 | endif
|
|---|
| 259 | mynm=seq(irs)
|
|---|
| 260 | call tolost(mynm)
|
|---|
| 261 | if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
|
|---|
| 262 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
|---|
| 263 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
|---|
| 264 | prlvr=.true. ! residue i is a proline variant
|
|---|
| 265 | else
|
|---|
| 266 | prlvr=.false.
|
|---|
| 267 | endif
|
|---|
| 268 |
|
|---|
| 269 | if ((mynm.eq.'arg').or.(mynm.eq.'arg+')) then
|
|---|
| 270 | do iat1=iatoff+iatrs1(irs)+13,iatrs2(irs)-2-iatmrg
|
|---|
| 271 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 272 | matcon(iat2-iat1,iat1)=0
|
|---|
| 273 | matcon(iat1-iat2,iat2)=0
|
|---|
| 274 | enddo
|
|---|
| 275 | enddo
|
|---|
| 276 | else if ((mynm.eq.'his').or.(mynm.eq.'hise')
|
|---|
| 277 | & .or.(mynm.eq.'hisd').or.(mynm.eq.'his+')) then
|
|---|
| 278 | do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
|
|---|
| 279 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 280 | matcon(iat2-iat1,iat1)=0
|
|---|
| 281 | matcon(iat1-iat2,iat2)=0
|
|---|
| 282 | enddo
|
|---|
| 283 | enddo
|
|---|
| 284 | else if (mynm.eq.'phe') then
|
|---|
| 285 | do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
|
|---|
| 286 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 287 | matcon(iat2-iat1,iat1)=0
|
|---|
| 288 | matcon(iat1-iat2,iat2)=0
|
|---|
| 289 | enddo
|
|---|
| 290 | enddo
|
|---|
| 291 | else if (mynm.eq.'tyr') then
|
|---|
| 292 | do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
|
|---|
| 293 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 294 | if (iat2.ne.iatrs1(irs)+15) then
|
|---|
| 295 | matcon(iat2-iat1,iat1)=0
|
|---|
| 296 | matcon(iat1-iat2,iat2)=0
|
|---|
| 297 | endif
|
|---|
| 298 | enddo
|
|---|
| 299 | enddo
|
|---|
| 300 | else if (mynm.eq.'trp') then
|
|---|
| 301 | do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
|
|---|
| 302 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 303 | matcon(iat2-iat1,iat1)=0
|
|---|
| 304 | matcon(iat1-iat2,iat2)=0
|
|---|
| 305 | enddo
|
|---|
| 306 | enddo
|
|---|
| 307 | else if (prlvr) then
|
|---|
| 308 | ! Proline. Many more distances are fixed because of the fixed
|
|---|
| 309 | ! phi angle
|
|---|
| 310 | do iat1=iatoff+iatrs1(irs),iatrs2(irs)-2-iatmrg
|
|---|
| 311 | do iat2=iat1+1,iatrs2(irs)-2-iatmrg
|
|---|
| 312 | matcon(iat2-iat1,iat1)=0
|
|---|
| 313 | matcon(iat1-iat2,iat2)=0
|
|---|
| 314 | enddo
|
|---|
| 315 | enddo
|
|---|
| 316 | ! distances to the C' atom of the previous residue are also fixed
|
|---|
| 317 | if (irs.ne.irsml1(iml)) then
|
|---|
| 318 | iat1=iowat(iatrs1(irs))
|
|---|
| 319 | do iat2=iatrs1(irs),iatrs2(irs)-1-iatmrg
|
|---|
| 320 | matcon(iat2-iat1,iat1)=0
|
|---|
| 321 | matcon(iat1-iat2,iat2)=0
|
|---|
| 322 | enddo
|
|---|
| 323 | endif
|
|---|
| 324 | endif
|
|---|
| 325 | enddo
|
|---|
| 326 | enddo
|
|---|
| 327 | ! finished initializing matrix conmat
|
|---|
| 328 | ! print *,'Connections matrix'
|
|---|
| 329 |
|
|---|
| 330 | ! Local pair excluded volume
|
|---|
| 331 | do i=1,mxml
|
|---|
| 332 | ilpst(i)=1
|
|---|
| 333 | ilpnd(i)=0
|
|---|
| 334 | enddo
|
|---|
| 335 | do i=1,50*mxrs
|
|---|
| 336 | lcp1(i)=0
|
|---|
| 337 | lcp2(i)=0
|
|---|
| 338 | enddo
|
|---|
| 339 | ilp=0
|
|---|
| 340 | do iml=1,ntlml
|
|---|
| 341 | do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
|
|---|
| 342 | do iat2=iat1+1,iatrs2(irsml2(iml))
|
|---|
| 343 | if ((iat2-iat1.le.mxconr).and.
|
|---|
| 344 | & matcon(iat2-iat1,iat1).eq.2) then
|
|---|
| 345 | ilp=ilp+1
|
|---|
| 346 | lcp1(ilp)=iat1
|
|---|
| 347 | lcp2(ilp)=iat2
|
|---|
| 348 | endif
|
|---|
| 349 | enddo
|
|---|
| 350 | enddo
|
|---|
| 351 |
|
|---|
| 352 | ilpnd(iml)=ilp
|
|---|
| 353 | if (iml.lt.ntlml) then
|
|---|
| 354 | ilpst(iml+1)=ilp+1
|
|---|
| 355 | endif
|
|---|
| 356 | ! print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
|
|---|
| 357 | ! print *,'local pair list'
|
|---|
| 358 | do lci=ilpst(iml),ilpnd(iml)
|
|---|
| 359 | iat1=lcp1(lci)
|
|---|
| 360 | iat2=lcp2(lci)
|
|---|
| 361 | ! print *,lci,iat1,iat2,matcon(iat2-iat1,iat1)
|
|---|
| 362 | enddo
|
|---|
| 363 | enddo
|
|---|
| 364 |
|
|---|
| 365 | print *,'finished initializing Lund force field'
|
|---|
| 366 | end
|
|---|
| 367 |
|
|---|
| 368 | integer function ihptype(iaa)
|
|---|
| 369 | include 'INCL.H'
|
|---|
| 370 | character mynm*4
|
|---|
| 371 | mynm=seq(iaa)
|
|---|
| 372 | ityp=-1
|
|---|
| 373 | call tolost(mynm)
|
|---|
| 374 | if (mynm.eq.'ala') then
|
|---|
| 375 | ityp=1
|
|---|
| 376 | else if ((mynm.eq.'val').or.(mynm.eq.'leu').or.(mynm.eq.'ile')
|
|---|
| 377 | & .or.(mynm.eq.'met').or.(mynm.eq.'pro').or.(mynm.eq.'cpro')
|
|---|
| 378 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
|---|
| 379 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
|---|
| 380 | ityp=2
|
|---|
| 381 | else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
|
|---|
| 382 | & then
|
|---|
| 383 | ityp=3
|
|---|
| 384 | endif
|
|---|
| 385 | ihptype=ityp
|
|---|
| 386 | return
|
|---|
| 387 | end
|
|---|
| 388 |
|
|---|
| 389 | real*8 function ebiasrs(irsd)
|
|---|
| 390 | include 'INCL.H'
|
|---|
| 391 | include 'incl_lund.h'
|
|---|
| 392 | dimension q1(2),q2(2)
|
|---|
| 393 | data q1/-0.2,0.2/
|
|---|
| 394 | data q2/0.42,-0.42/
|
|---|
| 395 |
|
|---|
| 396 | et=0.0
|
|---|
| 397 | do i=0,1
|
|---|
| 398 | iat1=iN(irsd)+i
|
|---|
| 399 | do j=0,1
|
|---|
| 400 | iat2=iC(irsd)+j
|
|---|
| 401 | xij=xat(iat1)-xat(iat2)
|
|---|
| 402 | yij=yat(iat1)-yat(iat2)
|
|---|
| 403 | zij=zat(iat1)-zat(iat2)
|
|---|
| 404 | et=et+q1(i+1)*q2(j+1)/sqrt(xij*xij+yij*yij+zij*zij)
|
|---|
| 405 | enddo
|
|---|
| 406 | enddo
|
|---|
| 407 | ebiasrs=kbias*et
|
|---|
| 408 | return
|
|---|
| 409 | end
|
|---|
| 410 | ! Evaluates backbone backbone hydrogen bond strength for residues
|
|---|
| 411 | ! i and j, taking the donor from residue i and acceptor from residue j
|
|---|
| 412 | real*8 function ehbmmrs(i,j)
|
|---|
| 413 | include 'INCL.H'
|
|---|
| 414 | include 'incl_lund.h'
|
|---|
| 415 | double precision r2,r4,r6,dx,dy,dz,ca,cb
|
|---|
| 416 | integer d1,d2,a1,a2
|
|---|
| 417 | d1=iN(i)
|
|---|
| 418 | d2=d1+1
|
|---|
| 419 | a1=iC(j)+1 ! dipoles are numbered from -ve to +ve charge
|
|---|
| 420 | a2=iC(j) ! so, acceptor 1 is the Oxygen, and a2, the carbon
|
|---|
| 421 | dx=xat(a1)-xat(d2)
|
|---|
| 422 | dy=yat(a1)-yat(d2)
|
|---|
| 423 | dz=zat(a1)-zat(d2)
|
|---|
| 424 | r2=dx*dx+dy*dy+dz*dz
|
|---|
| 425 | if (r2.gt.cthb2) then
|
|---|
| 426 | ! print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2
|
|---|
| 427 | ! print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb
|
|---|
| 428 | ehbmmrs=0
|
|---|
| 429 | return
|
|---|
| 430 | endif
|
|---|
| 431 | ca=(xat(d2)-xat(d1))*dx+(yat(d2)-yat(d1))*dy+(zat(d2)-zat(d1))*dz
|
|---|
| 432 | cb=(xat(a2)-xat(a1))*dx+(yat(a2)-yat(a1))*dy+(zat(a2)-zat(a1))*dz
|
|---|
| 433 | if (powa.gt.0.and.ca.le.0) then
|
|---|
| 434 | ! print *,'hbmm, returning 0 because of angle a'
|
|---|
| 435 | ehbmmrs=0
|
|---|
| 436 | return
|
|---|
| 437 | endif
|
|---|
| 438 | if (powb.gt.0.and.cb.le.0) then
|
|---|
| 439 | ! print *,'hbmm, returning 0 because of angle b'
|
|---|
| 440 | ehbmmrs=0
|
|---|
| 441 | return
|
|---|
| 442 | endif
|
|---|
| 443 | r6=sighb2/r2
|
|---|
| 444 | r4=r6*r6
|
|---|
| 445 | r6=r6*r4
|
|---|
| 446 | evlu=((ca*ca/r2)**(0.5*powa))*((cb*cb/r2)**(0.5*powb))
|
|---|
| 447 | evlu=evlu*(r6*(5*r6-6*r4)+alhb+blhb*r2)
|
|---|
| 448 | ! print *,'found hbmm contribution ',evlu
|
|---|
| 449 | ehbmmrs=epshb1*evlu
|
|---|
| 450 | return
|
|---|
| 451 | end
|
|---|
| 452 | real*8 function enylun(nml)
|
|---|
| 453 | ! nml = 1 .. ntlml. No provision exists to handle out of range values
|
|---|
| 454 | ! for nml inside this function.
|
|---|
| 455 | include 'INCL.H'
|
|---|
| 456 | include 'incl_lund.h'
|
|---|
| 457 | character mynm*4
|
|---|
| 458 | logical prlvr
|
|---|
| 459 | eyhb=0.0 ! backbone-backbone and sidechain-backbone HB
|
|---|
| 460 | eyel=0.0 ! Bias, or local electrostatic term
|
|---|
| 461 | eysl=0.0 ! Hydrophobicity, replaces the solvent term of SMMP
|
|---|
| 462 | eyvr=0.0 ! Local pair excluded volume, in a sense a variable potential
|
|---|
| 463 | eyvw=0.0 ! atom-atom repulsion, excluded volume
|
|---|
| 464 | ! atom-atom repulsion is calculated on a system wide basis, instead of
|
|---|
| 465 | ! molecule by molecule for efficiency. Look into function exvlun.
|
|---|
| 466 |
|
|---|
| 467 | istres=irsml1(nml)
|
|---|
| 468 | indres=irsml2(nml)
|
|---|
| 469 |
|
|---|
| 470 | ! First, all terms that can be calculated on a residue by residue basis
|
|---|
| 471 | do i=istres,indres
|
|---|
| 472 | mynm=seq(i)
|
|---|
| 473 | call tolost(mynm)
|
|---|
| 474 | if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
|
|---|
| 475 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
|---|
| 476 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
|---|
| 477 | prlvr=.true. ! residue i is a proline variant
|
|---|
| 478 | else
|
|---|
| 479 | prlvr=.false.
|
|---|
| 480 | endif
|
|---|
| 481 |
|
|---|
| 482 | ! Bias, or local electrostatic term. Excluded from the list are
|
|---|
| 483 | ! residues at the ends of the chain, glycine and all proline variants
|
|---|
| 484 | if ((i.ne.istres).and.(i.ne.indres).and.
|
|---|
| 485 | & .not.prlvr.and.mynm.ne.'gly') then
|
|---|
| 486 | eyel=eyel+ebiasrs(i)
|
|---|
| 487 | endif
|
|---|
| 488 | ! Backbone--backbone hydrogen bonds
|
|---|
| 489 | shbm1=1.0
|
|---|
| 490 | shbm2=1.0
|
|---|
| 491 | if ((i.eq.istres).or.(i.eq.indres)) shbm1=0.5
|
|---|
| 492 | ! Residue i contributes the donor, and j, the acceptor, so both i and
|
|---|
| 493 | ! j run over the whole set of amino acids.
|
|---|
| 494 | ! No terms for residue i, if it is a proline variant.
|
|---|
| 495 | if (.not.prlvr) then
|
|---|
| 496 | do j=istres,indres
|
|---|
| 497 | if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
|
|---|
| 498 | etmp=ehbmmrs(i,j)
|
|---|
| 499 | eyhb=eyhb+shbm1*shbm2*etmp
|
|---|
| 500 | enddo
|
|---|
| 501 | endif
|
|---|
| 502 | ! Hydrophobicity, only if residue i is hydrophobic to start with
|
|---|
| 503 | ihpi=ihptype(i)
|
|---|
| 504 | if (ihpi.ge.0) then
|
|---|
| 505 | ! Unlike hydrogen bonds, the hydrophobicity potential is symmetric
|
|---|
| 506 | ! in i and j. So, the loop for j runs from i+1 to the end.
|
|---|
| 507 |
|
|---|
| 508 | do j=i+1,indres
|
|---|
| 509 | ihpj=ihptype(j)
|
|---|
| 510 | if (ihpj.ge.0) then
|
|---|
| 511 | etmp=ehp(i,j,ihpi,ihpj)
|
|---|
| 512 | if (j.eq.(i+1)) etmp=0
|
|---|
| 513 | if (j.eq.(i+2)) etmp=0.5*etmp
|
|---|
| 514 | eysl=eysl+etmp
|
|---|
| 515 | endif
|
|---|
| 516 | enddo
|
|---|
| 517 | endif
|
|---|
| 518 | enddo
|
|---|
| 519 |
|
|---|
| 520 | ! Terms that are not calculated residue by residue ...
|
|---|
| 521 |
|
|---|
| 522 | ! Local pair or third-neighbour excluded volume
|
|---|
| 523 | ! Numerically this is normally the term with largest positive
|
|---|
| 524 | ! contribution to the energy in an equilibrated stystem.
|
|---|
| 525 |
|
|---|
| 526 | i1=ilpst(nml)
|
|---|
| 527 | i2=ilpnd(nml)
|
|---|
| 528 | etmp=0.0
|
|---|
| 529 | do i=i1,i2
|
|---|
| 530 | etmp1=0.0
|
|---|
| 531 | iat1=lcp1(i)
|
|---|
| 532 | iat2=lcp2(i)
|
|---|
| 533 | iatt1=ityat(iat1)
|
|---|
| 534 | iatt2=ityat(iat2)
|
|---|
| 535 | xij=xat(iat1)-xat(iat2)
|
|---|
| 536 | yij=yat(iat1)-yat(iat2)
|
|---|
| 537 | zij=zat(iat1)-zat(iat2)
|
|---|
| 538 | r2=(xij*xij + yij*yij + zij*zij)
|
|---|
| 539 | if (r2.le.exvcut2) then
|
|---|
| 540 | r6=sig2lcp(iatt1,iatt2)/r2
|
|---|
| 541 | r6=r6*r6*r6
|
|---|
| 542 | etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2)
|
|---|
| 543 | if (etmp1.ge.500) then
|
|---|
| 544 | print *,'local pair contribution ',iat1,iat2,etmp1
|
|---|
| 545 | endif
|
|---|
| 546 | etmp=etmp+etmp1
|
|---|
| 547 | endif
|
|---|
| 548 | ! print *,'pair : ',iat1,iat2,' contribution ',etmp1
|
|---|
| 549 | ! print *,exvcut2,r2
|
|---|
| 550 | enddo
|
|---|
| 551 | eyvr=exvk*etmp
|
|---|
| 552 | eysm=eyel+eyhb+eyvr+eysl
|
|---|
| 553 | enylun=eysm
|
|---|
| 554 | end
|
|---|
| 555 | real*8 function eninlun()
|
|---|
| 556 | include 'INCL.H'
|
|---|
| 557 | eysmi=0.0
|
|---|
| 558 | eyeli=0.0
|
|---|
| 559 | eyvwi=0.0
|
|---|
| 560 | eyhbi=0.0
|
|---|
| 561 | eninlun=0.0
|
|---|
| 562 | return
|
|---|
| 563 | end
|
|---|
| 564 | real*8 function ehp(i1,i2,ihp1,ihp2)
|
|---|
| 565 | include 'INCL.H'
|
|---|
| 566 | include 'incl_lund.h'
|
|---|
| 567 | dimension r2min(12)
|
|---|
| 568 |
|
|---|
| 569 | b2=20.25
|
|---|
| 570 | a2=12.25
|
|---|
| 571 | ! ihp1=ihptype(i1)
|
|---|
| 572 | ! ihp2=ihptype(i2)
|
|---|
| 573 | if ((ihp1.le.0).or.(ihp2.le.0)) then
|
|---|
| 574 | ehp=0.0
|
|---|
| 575 | return
|
|---|
| 576 | endif
|
|---|
| 577 | ni=nhpat(i1)
|
|---|
| 578 | nj=nhpat(i2)
|
|---|
| 579 | do i=1,ni+nj
|
|---|
| 580 | r2min(i)=100000000 ! quite far
|
|---|
| 581 | enddo
|
|---|
| 582 | do i=1,ni
|
|---|
| 583 | k=ihpat(i1,i)
|
|---|
| 584 | do j=1,nj
|
|---|
| 585 | l=ihpat(i2,j)
|
|---|
| 586 | xij=xat(k)-xat(l)
|
|---|
| 587 | yij=yat(k)-yat(l)
|
|---|
| 588 | zij=zat(k)-zat(l)
|
|---|
| 589 | dtmp=(xij*xij + yij*yij + zij*zij)
|
|---|
| 590 | if (dtmp.le.r2min(i)) r2min(i)=dtmp
|
|---|
| 591 | if (dtmp.le.r2min(ni+j)) r2min(ni+j)=dtmp
|
|---|
| 592 | enddo
|
|---|
| 593 | enddo
|
|---|
| 594 | sum=0
|
|---|
| 595 | do i=1,ni+nj
|
|---|
| 596 | if (r2min(i).le.b2) then
|
|---|
| 597 | if (r2min(i).lt.a2) then
|
|---|
| 598 | sum=sum+1
|
|---|
| 599 | else
|
|---|
| 600 | sum=sum+(b2-r2min(i))/(b2-a2)
|
|---|
| 601 | endif
|
|---|
| 602 | endif
|
|---|
| 603 | enddo
|
|---|
| 604 | ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*sum/(ni+nj)
|
|---|
| 605 | return
|
|---|
| 606 | end
|
|---|
| 607 |
|
|---|
| 608 | real*8 function exvlun(nml)
|
|---|
| 609 | include 'INCL.H'
|
|---|
| 610 | include 'incl_lund.h'
|
|---|
| 611 | ! For multi-chain systems it makes little sense to split the calculation
|
|---|
| 612 | ! of this term into an 'interaction part' and a contribution from
|
|---|
| 613 | ! individual molecules. So, normally this should always be called with
|
|---|
| 614 | ! argument nml=0. Only for diagnostic reasons, you might want to find
|
|---|
| 615 | ! the contribution from one molecule in a multi-chain system assuming
|
|---|
| 616 | ! there was no other molecule.
|
|---|
| 617 | dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
|
|---|
| 618 | dimension icell(mxat)
|
|---|
| 619 | if (nml.eq.0) then
|
|---|
| 620 | istat=iatrs1(irsml1(1))
|
|---|
| 621 | indat=iatrs2(irsml2(ntlml))
|
|---|
| 622 | else
|
|---|
| 623 | istat=iatrs1(irsml1(nml))
|
|---|
| 624 | indat=iatrs2(irsml2(nml))
|
|---|
| 625 | endif
|
|---|
| 626 |
|
|---|
| 627 | eyvw=0.0
|
|---|
| 628 | ! The beginning part of this implementation is very similar to the
|
|---|
| 629 | ! assignment of cells to the atoms during calculation of solvent
|
|---|
| 630 | ! accessible surface area. So, much of that part is similar. But
|
|---|
| 631 | ! unlike the accessible surface calculations, this term is symmetric
|
|---|
| 632 | ! in any two participating atoms. So, the part after the assignment
|
|---|
| 633 | ! of cells differs even in the structure of the code.
|
|---|
| 634 |
|
|---|
| 635 | do i=1,mxcell
|
|---|
| 636 | incell(i)=0
|
|---|
| 637 | enddo
|
|---|
| 638 | ! print *,'evaluating general excluded volume :',istat,',',indat
|
|---|
| 639 | ! Find minimal containing box
|
|---|
| 640 | xmin=xat(istat)
|
|---|
| 641 | ymin=yat(istat)
|
|---|
| 642 | zmin=zat(istat)
|
|---|
| 643 | xmax=xmin
|
|---|
| 644 | ymax=ymin
|
|---|
| 645 | zmax=zmin
|
|---|
| 646 |
|
|---|
| 647 | do i=istat,indat
|
|---|
| 648 | if (xat(i).le.xmin) then
|
|---|
| 649 | xmin=xat(i)
|
|---|
| 650 | else if (xat(i).ge.xmax) then
|
|---|
| 651 | xmax=xat(i)
|
|---|
| 652 | endif
|
|---|
| 653 | if (yat(i).le.ymin) then
|
|---|
| 654 | ymin=yat(i)
|
|---|
| 655 | else if (yat(i).ge.ymax) then
|
|---|
| 656 | ymax=yat(i)
|
|---|
| 657 | endif
|
|---|
| 658 | if (zat(i).le.zmin) then
|
|---|
| 659 | zmin=zat(i)
|
|---|
| 660 | else if (zat(i).ge.zmax) then
|
|---|
| 661 | zmax=zat(i)
|
|---|
| 662 | endif
|
|---|
| 663 | enddo
|
|---|
| 664 |
|
|---|
| 665 | sizex=xmax-xmin
|
|---|
| 666 | sizey=ymax-ymin
|
|---|
| 667 | sizez=zmax-zmin
|
|---|
| 668 | ! Number of cells along each directions that fit into the box.
|
|---|
| 669 | ndx=int(sizex/exvcutg)+1
|
|---|
| 670 | ndy=int(sizey/exvcutg)+1
|
|---|
| 671 | ndz=int(sizez/exvcutg)+1
|
|---|
| 672 |
|
|---|
| 673 | nxy=ndx*ndy
|
|---|
| 674 | ncell=nxy*ndz
|
|---|
| 675 | ! print *,'Number of cells along x,y,z = ',ndx,',',ndy,',',ndz
|
|---|
| 676 | if (ncell.ge.mxcell) then
|
|---|
| 677 | print *,'exvlun> required number of cells',ncell,
|
|---|
| 678 | & ' exceeded the limit ',mxcell
|
|---|
| 679 | print *,'recompile with a higher mxcell.'
|
|---|
| 680 | stop
|
|---|
| 681 | endif
|
|---|
| 682 | ! Expand box to contain an integral number of cells along each direction
|
|---|
| 683 | shiftx=(dble(ndx)*exvcutg-sizex)/2.0
|
|---|
| 684 | shifty=(dble(ndy)*exvcutg-sizey)/2.0
|
|---|
| 685 | shiftz=(dble(ndz)*exvcutg-sizez)/2.0
|
|---|
| 686 | xmin=xmin-shiftx
|
|---|
| 687 | ymin=ymin-shifty
|
|---|
| 688 | zmin=zmin-shiftz
|
|---|
| 689 | xmax=xmax+shiftx
|
|---|
| 690 | ymax=ymax+shifty
|
|---|
| 691 | zmax=zmax+shiftz
|
|---|
| 692 |
|
|---|
| 693 | ! Set occupied cells to zero. Note that the maximum number of occupied
|
|---|
| 694 | ! cells is the number of atoms in the system.
|
|---|
| 695 | nocccl=0
|
|---|
| 696 | do i=1,mxat
|
|---|
| 697 | locccl(i)=0
|
|---|
| 698 | enddo
|
|---|
| 699 |
|
|---|
| 700 | ! Put atoms in cells
|
|---|
| 701 | do j=istat,indat
|
|---|
| 702 | mx=min(int(max((xat(j)-xmin)/exvcutg,0.0d0)),ndx-1)
|
|---|
| 703 | my=min(int(max((yat(j)-ymin)/exvcutg,0.0d0)),ndy-1)
|
|---|
| 704 | mz=min(int(max((zat(j)-zmin)/exvcutg,0.0d0)),ndz-1)
|
|---|
| 705 | icellj=mx+my*ndx+mz*nxy+1
|
|---|
| 706 | icell(j)=icellj
|
|---|
| 707 | if (icellj.gt.mxcell) then
|
|---|
| 708 | print *,'exvlun> bad cell index ',icellj,' for atom ',j
|
|---|
| 709 | stop
|
|---|
| 710 | else
|
|---|
| 711 | if (incell(icellj).eq.0) then
|
|---|
| 712 | ! previously unoccupied cell
|
|---|
| 713 | nocccl=nocccl+1
|
|---|
| 714 | locccl(nocccl)=icellj
|
|---|
| 715 | endif
|
|---|
| 716 | incell(icellj)=incell(icellj)+1
|
|---|
| 717 | endif
|
|---|
| 718 | enddo
|
|---|
| 719 | ! print *,'finished assigning cells. nocccl = ',nocccl
|
|---|
| 720 | ! Cummulative occupancy of i'th cell
|
|---|
| 721 | do i=1,ncell
|
|---|
| 722 | incell(i+1)=incell(i+1)+incell(i)
|
|---|
| 723 | enddo
|
|---|
| 724 | ! print *,'finished making cumulative cell sums'
|
|---|
| 725 | ! Sorting atoms by their cell index
|
|---|
| 726 | do i=istat,indat
|
|---|
| 727 | j=icell(i)
|
|---|
| 728 | jj=incell(j)
|
|---|
| 729 | isort(jj)=i
|
|---|
| 730 | incell(j)=jj-1
|
|---|
| 731 | enddo
|
|---|
| 732 | ! print *,'sorted atoms by cell index'
|
|---|
| 733 | etmp=0.0
|
|---|
| 734 | do icl=1,nocccl
|
|---|
| 735 | ! loop through occupied cells
|
|---|
| 736 | lcell=locccl(icl)
|
|---|
| 737 | ix=mod(lcell-1,ndx)
|
|---|
| 738 | iy=(mod(lcell-1,nxy)-ix)/ndx
|
|---|
| 739 | iz=(lcell-1-ix-ndx*iy)/nxy
|
|---|
| 740 | ! print *,'icl=',icl,'absolute index of cell = ',lcell
|
|---|
| 741 | ! print *,'iz,iy,ix = ',iz,iy,ix
|
|---|
| 742 | ! find all atoms in current cell and all its forward-going neighbours
|
|---|
| 743 | nex=min(ix+1,ndx-1)
|
|---|
| 744 | ney=min(iy+1,ndy-1)
|
|---|
| 745 | nez=min(iz+1,ndz-1)
|
|---|
| 746 | nsame=0
|
|---|
| 747 | nngbr=0
|
|---|
| 748 | do jx=ix,nex
|
|---|
| 749 | do jy=iy,ney
|
|---|
| 750 | do jz=iz,nez
|
|---|
| 751 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 752 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 753 | ! count the total number of neighbours
|
|---|
| 754 | nngbr=nngbr+1
|
|---|
| 755 | if (jx.eq.ix.and.jy.eq.iy.and.jz.eq.iz) then
|
|---|
| 756 | ! count how many neighbours are from the same cell
|
|---|
| 757 | nsame=nsame+1
|
|---|
| 758 | endif
|
|---|
| 759 | ngbr(nngbr)=isort(ii)
|
|---|
| 760 | enddo
|
|---|
| 761 | enddo
|
|---|
| 762 | enddo
|
|---|
| 763 | enddo
|
|---|
| 764 | ! A few more cells need to be searched, so that we cover 13 of the 26
|
|---|
| 765 | ! neighbouring cells.
|
|---|
| 766 | ! 1
|
|---|
| 767 | jx=ix+1
|
|---|
| 768 | jy=iy
|
|---|
| 769 | jz=iz-1
|
|---|
| 770 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 771 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 772 | nngbr=nngbr+1
|
|---|
| 773 | ngbr(nngbr)=isort(ii)
|
|---|
| 774 | enddo
|
|---|
| 775 | ! 2
|
|---|
| 776 | jx=ix
|
|---|
| 777 | jy=iy-1
|
|---|
| 778 | jz=iz+1
|
|---|
| 779 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 780 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 781 | nngbr=nngbr+1
|
|---|
| 782 | ngbr(nngbr)=isort(ii)
|
|---|
| 783 | enddo
|
|---|
| 784 | ! 3
|
|---|
| 785 | jx=ix-1
|
|---|
| 786 | jy=iy+1
|
|---|
| 787 | jz=iz
|
|---|
| 788 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 789 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 790 | nngbr=nngbr+1
|
|---|
| 791 | ngbr(nngbr)=isort(ii)
|
|---|
| 792 | enddo
|
|---|
| 793 | ! 4
|
|---|
| 794 | jx=ix+1
|
|---|
| 795 | jy=iy+1
|
|---|
| 796 | jz=iz-1
|
|---|
| 797 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 798 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 799 | nngbr=nngbr+1
|
|---|
| 800 | ngbr(nngbr)=isort(ii)
|
|---|
| 801 | enddo
|
|---|
| 802 | ! 5
|
|---|
| 803 | jx=ix+1
|
|---|
| 804 | jy=iy-1
|
|---|
| 805 | jz=iz+1
|
|---|
| 806 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 807 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 808 | nngbr=nngbr+1
|
|---|
| 809 | ngbr(nngbr)=isort(ii)
|
|---|
| 810 | enddo
|
|---|
| 811 | ! 6
|
|---|
| 812 | jx=ix+1
|
|---|
| 813 | jy=iy-1
|
|---|
| 814 | jz=iz-1
|
|---|
| 815 | jcl=jx+ndx*jy+nxy*jz+1
|
|---|
| 816 | do ii=incell(jcl)+1,incell(jcl+1)
|
|---|
| 817 | nngbr=nngbr+1
|
|---|
| 818 | ngbr(nngbr)=isort(ii)
|
|---|
| 819 | enddo
|
|---|
| 820 |
|
|---|
| 821 | ! print *,'atoms in same cell ',nsame
|
|---|
| 822 | ! print *,'atoms in neighbouring cells ',nngbr
|
|---|
| 823 | do i1=1,nsame
|
|---|
| 824 | ! Over all atoms from the original cell
|
|---|
| 825 | iat1=ngbr(i1)
|
|---|
| 826 | do i2=i1,nngbr
|
|---|
| 827 | ! Over all atoms in the original+neighbouring cells
|
|---|
| 828 | iat2=ngbr(i2)
|
|---|
| 829 | xij=xat(iat1)-xat(iat2)
|
|---|
| 830 | yij=yat(iat1)-yat(iat2)
|
|---|
| 831 | zij=zat(iat1)-zat(iat2)
|
|---|
| 832 | r2=(xij*xij+yij*yij+zij*zij)
|
|---|
| 833 |
|
|---|
| 834 | if (r2.le.exvcutg2) then
|
|---|
| 835 | if (abs(iat2-iat1).gt.mxconr.or.
|
|---|
| 836 | & matcon(iat2-iat1,iat1).eq.1) then
|
|---|
| 837 | iatt1=ityat(iat1)
|
|---|
| 838 | iatt2=ityat(iat2)
|
|---|
| 839 | r6=sig2exv(iatt1,iatt2)/r2
|
|---|
| 840 | r6=r6*r6*r6
|
|---|
| 841 | etmp1=r6*r6+asaexv(iatt1,iatt2)
|
|---|
| 842 | & +bsaexv(iatt1,iatt2)*r2
|
|---|
| 843 | etmp=etmp+etmp1
|
|---|
| 844 | ! if (etmp1.ge.2000) then
|
|---|
| 845 | ! print *,'contribution ',iat1,iat2,etmp1
|
|---|
| 846 | ! call outpdb(1,'EXAMPLES/clash.pdb')
|
|---|
| 847 | ! stop
|
|---|
| 848 | ! endif
|
|---|
| 849 | endif
|
|---|
| 850 | endif
|
|---|
| 851 | enddo
|
|---|
| 852 | enddo
|
|---|
| 853 | enddo
|
|---|
| 854 | ! irs=1
|
|---|
| 855 | ! do iat=iatrs1(irs),iatrs2(irs)
|
|---|
| 856 | ! do j=-mxconr,mxconr
|
|---|
| 857 | ! print *,iat,j,':',matcon(j,iat)
|
|---|
| 858 | ! enddo
|
|---|
| 859 | ! enddo
|
|---|
| 860 | ! irs=irsml2(1)
|
|---|
| 861 | ! do iat=iatrs1(irs),iatrs2(irs)
|
|---|
| 862 | ! do j=-mxconr,mxconr
|
|---|
| 863 | ! print *,iat,j,':',matcon(j,iat)
|
|---|
| 864 | ! enddo
|
|---|
| 865 | ! enddo
|
|---|
| 866 |
|
|---|
| 867 | eyvw=exvk*etmp
|
|---|
| 868 | exvlun=eyvw
|
|---|
| 869 | return
|
|---|
| 870 | end
|
|---|
| 871 |
|
|---|
| 872 | real*8 function exvbrfc()
|
|---|
| 873 | ! Brute force excluded volume evaluation
|
|---|
| 874 | include 'INCL.H'
|
|---|
| 875 | include 'incl_lund.h'
|
|---|
| 876 |
|
|---|
| 877 | etmp=0.0
|
|---|
| 878 | etmp1=0.0
|
|---|
| 879 | print *,'max connection radius is ',mxconr
|
|---|
| 880 | do iat1=iatrs1(irsml1(1)),iatrs2(irsml2(ntlml))
|
|---|
| 881 | do iat2=iat1+1,iatrs2(irsml2(ntlml))
|
|---|
| 882 | xij=xat(iat1)-xat(iat2)
|
|---|
| 883 | yij=yat(iat1)-yat(iat2)
|
|---|
| 884 | zij=zat(iat1)-zat(iat2)
|
|---|
| 885 | r2=(xij*xij+yij*yij+zij*zij)
|
|---|
| 886 |
|
|---|
| 887 | if (r2.le.exvcutg2) then
|
|---|
| 888 | if (abs(iat2-iat1).gt.mxconr.or.
|
|---|
| 889 | & matcon(iat2-iat1,iat1).eq.1) then
|
|---|
| 890 | iatt1=ityat(iat1)
|
|---|
| 891 | iatt2=ityat(iat2)
|
|---|
| 892 | r6=sig2exv(iatt1,iatt2)/r2
|
|---|
| 893 | r6=r6*r6*r6
|
|---|
| 894 | etmp1=r6*r6+asaexv(iatt1,iatt2)
|
|---|
| 895 | & +bsaexv(iatt1,iatt2)*r2
|
|---|
| 896 | etmp=etmp+etmp1
|
|---|
| 897 | if (iat1.eq.43.and.iat2.eq.785) then
|
|---|
| 898 | print *,'contribution ',iat1,iat2,etmp1
|
|---|
| 899 | print *,'r2 = ',r2
|
|---|
| 900 | print *,'atom types : ',iatt1,iatt2
|
|---|
| 901 | print *,'sig2exv = ',sig2exv(iatt1,iatt2)
|
|---|
| 902 | print *,'asa = ',asaexv(iatt1,iatt2)
|
|---|
| 903 | print *,'bsa = ',bsaexv(iatt1,iatt2)
|
|---|
| 904 |
|
|---|
| 905 | ! call outpdb(1,'EXAMPLES/clash.pdb')
|
|---|
| 906 | ! stop
|
|---|
| 907 | endif
|
|---|
| 908 | else
|
|---|
| 909 | ! print *,'atoms ', iat1,' and ',iat2,' were close',
|
|---|
| 910 | ! # 'but matcon is ',matcon(iat2-iat1,iat1)
|
|---|
| 911 | endif
|
|---|
| 912 | endif
|
|---|
| 913 | enddo
|
|---|
| 914 | enddo
|
|---|
| 915 | exvbrfc=etmp*exvk
|
|---|
| 916 | return
|
|---|
| 917 | end
|
|---|
