| 1 | ! *********************************************************************
|
|---|
| 2 | ! This file contains eninteract
|
|---|
| 3 | !
|
|---|
| 4 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 5 | ! Jan H. Meinke, Sandipan Mohanty
|
|---|
| 6 | !
|
|---|
| 7 | !
|
|---|
| 8 | ! Description: Calculates the interaction energy between molecules
|
|---|
| 9 | ! The function assumes that all molecules are up-to-date. If in doubt
|
|---|
| 10 | ! call energy first.
|
|---|
| 11 | ! The energy function is based on the ECEPP/3 dataset.
|
|---|
| 12 | !
|
|---|
| 13 | ! TODO: Intermolecular interaction energy for FLEX and ECEPP/2
|
|---|
| 14 | real*8 function eninteract()
|
|---|
| 15 |
|
|---|
| 16 | include 'INCL.H'
|
|---|
| 17 |
|
|---|
| 18 | eysmi=0.0
|
|---|
| 19 | eyeli=0.0
|
|---|
| 20 | eyvwi=0.0
|
|---|
| 21 | eyhbi=0.0
|
|---|
| 22 |
|
|---|
| 23 | do iml = 1, ntlml
|
|---|
| 24 | do jml = iml + 1, ntlml
|
|---|
| 25 | do ires= irsml1(iml), irsml2(iml)
|
|---|
| 26 | do jres= irsml1(jml), irsml2(jml)
|
|---|
| 27 | do iat = iatrs1(ires), iatrs2(ires)
|
|---|
| 28 | ! Atom class of current atom
|
|---|
| 29 | ity=ityat(iat)
|
|---|
| 30 | ! Point charge at current atom
|
|---|
| 31 | cqi=conv*cgat(iat)
|
|---|
| 32 | ! Cartesian coordinates of current atom
|
|---|
| 33 | xi=xat(iat)
|
|---|
| 34 | yi=yat(iat)
|
|---|
| 35 | zi=zat(iat)
|
|---|
| 36 |
|
|---|
| 37 | do jat = iatrs1(jres), iatrs2(jres)
|
|---|
| 38 | ! Atom type of partner
|
|---|
| 39 | jty=ityat(jat)
|
|---|
| 40 | ! Differences in cartesian coordinates
|
|---|
| 41 | xj=xat(jat)
|
|---|
| 42 | yj=yat(jat)
|
|---|
| 43 | zj=zat(jat)
|
|---|
| 44 |
|
|---|
| 45 | xij=xat(jat)-xi
|
|---|
| 46 | yij=yat(jat)-yi
|
|---|
| 47 | zij=zat(jat)-zi
|
|---|
| 48 | ! Cartesian distance and higher powers
|
|---|
| 49 | rij2=xij*xij+yij*yij+zij*zij
|
|---|
| 50 | rij4=rij2*rij2
|
|---|
| 51 | rij6=rij4*rij2
|
|---|
| 52 | rij=sqrt(rij2)
|
|---|
| 53 | ! Are we using a distance dependent dielectric constant?
|
|---|
| 54 | if(epsd) then
|
|---|
| 55 | sr=slp*rij
|
|---|
| 56 | ep=plt-(sr*sr+2.0*sr+2.0)*(plt-1.0)*exp(-sr)/2.0
|
|---|
| 57 | else
|
|---|
| 58 | ep = 1.0d0
|
|---|
| 59 | end if
|
|---|
| 60 | ! Coulomb interaction
|
|---|
| 61 | eyeli=eyeli+cqi*cgat(jat)/(rij*ep)
|
|---|
| 62 | ! If the two atoms cannot form a hydrogen bond use 6-12 Lennard-Jones potential
|
|---|
| 63 | if (ihbty(ity,jty).eq.0) then
|
|---|
| 64 | eyvwi=eyvwi+aij(ity,jty)/(rij6*rij6)
|
|---|
| 65 | & -cij(ity,jty)/rij6
|
|---|
| 66 | else
|
|---|
| 67 | ! For hydrogen bonding use 10-12 Lennard-Jones potential
|
|---|
| 68 | eyhbi=eyhbi+ahb(ity,jty)/(rij6*rij6)
|
|---|
| 69 | & -chb(ity,jty)/(rij6*rij4)
|
|---|
| 70 | endif
|
|---|
| 71 | enddo
|
|---|
| 72 | enddo
|
|---|
| 73 | enddo
|
|---|
| 74 | enddo
|
|---|
| 75 | enddo
|
|---|
| 76 | enddo
|
|---|
| 77 |
|
|---|
| 78 | eysmi = eyeli + eyvwi + eyhbi
|
|---|
| 79 | eninteract = eysmi
|
|---|
| 80 | return
|
|---|
| 81 | end |
|---|
