| [bd2278d] | 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
|
|---|
| [e40e335] | 14 | real*8 function eninteract()
|
|---|
| 15 |
|
|---|
| 16 | include 'INCL.H'
|
|---|
| [32289cd] | 17 |
|
|---|
| 18 | double precision cqi, xi, yi, zi, xj, yj, zj, xij, yij, zij, rij2
|
|---|
| 19 | double precision rij4, rij6, rij, sr, ep
|
|---|
| 20 |
|
|---|
| 21 | integer iml, jml, ires, jres, iat, ity, jat, jty
|
|---|
| 22 |
|
|---|
| [e40e335] | 23 | eysmi=0.0
|
|---|
| 24 | eyeli=0.0
|
|---|
| 25 | eyvwi=0.0
|
|---|
| 26 | eyhbi=0.0
|
|---|
| 27 |
|
|---|
| 28 | do iml = 1, ntlml
|
|---|
| 29 | do jml = iml + 1, ntlml
|
|---|
| 30 | do ires= irsml1(iml), irsml2(iml)
|
|---|
| 31 | do jres= irsml1(jml), irsml2(jml)
|
|---|
| 32 | do iat = iatrs1(ires), iatrs2(ires)
|
|---|
| [bd2278d] | 33 | ! Atom class of current atom
|
|---|
| [e40e335] | 34 | ity=ityat(iat)
|
|---|
| [bd2278d] | 35 | ! Point charge at current atom
|
|---|
| [e40e335] | 36 | cqi=conv*cgat(iat)
|
|---|
| [bd2278d] | 37 | ! Cartesian coordinates of current atom
|
|---|
| [e40e335] | 38 | xi=xat(iat)
|
|---|
| 39 | yi=yat(iat)
|
|---|
| 40 | zi=zat(iat)
|
|---|
| 41 |
|
|---|
| 42 | do jat = iatrs1(jres), iatrs2(jres)
|
|---|
| [bd2278d] | 43 | ! Atom type of partner
|
|---|
| [e40e335] | 44 | jty=ityat(jat)
|
|---|
| [bd2278d] | 45 | ! Differences in cartesian coordinates
|
|---|
| [e40e335] | 46 | xj=xat(jat)
|
|---|
| 47 | yj=yat(jat)
|
|---|
| 48 | zj=zat(jat)
|
|---|
| [32289cd] | 49 |
|
|---|
| [e40e335] | 50 | xij=xat(jat)-xi
|
|---|
| 51 | yij=yat(jat)-yi
|
|---|
| 52 | zij=zat(jat)-zi
|
|---|
| [bd2278d] | 53 | ! Cartesian distance and higher powers
|
|---|
| [e40e335] | 54 | rij2=xij*xij+yij*yij+zij*zij
|
|---|
| 55 | rij4=rij2*rij2
|
|---|
| 56 | rij6=rij4*rij2
|
|---|
| 57 | rij=sqrt(rij2)
|
|---|
| [bd2278d] | 58 | ! Are we using a distance dependent dielectric constant?
|
|---|
| [e40e335] | 59 | if(epsd) then
|
|---|
| 60 | sr=slp*rij
|
|---|
| 61 | ep=plt-(sr*sr+2.0*sr+2.0)*(plt-1.0)*exp(-sr)/2.0
|
|---|
| 62 | else
|
|---|
| 63 | ep = 1.0d0
|
|---|
| 64 | end if
|
|---|
| [bd2278d] | 65 | ! Coulomb interaction
|
|---|
| [e40e335] | 66 | eyeli=eyeli+cqi*cgat(jat)/(rij*ep)
|
|---|
| [bd2278d] | 67 | ! If the two atoms cannot form a hydrogen bond use 6-12 Lennard-Jones potential
|
|---|
| [e40e335] | 68 | if (ihbty(ity,jty).eq.0) then
|
|---|
| 69 | eyvwi=eyvwi+aij(ity,jty)/(rij6*rij6)
|
|---|
| [bd2278d] | 70 | & -cij(ity,jty)/rij6
|
|---|
| [e40e335] | 71 | else
|
|---|
| [bd2278d] | 72 | ! For hydrogen bonding use 10-12 Lennard-Jones potential
|
|---|
| [e40e335] | 73 | eyhbi=eyhbi+ahb(ity,jty)/(rij6*rij6)
|
|---|
| [bd2278d] | 74 | & -chb(ity,jty)/(rij6*rij4)
|
|---|
| [e40e335] | 75 | endif
|
|---|
| [32289cd] | 76 | enddo
|
|---|
| [e40e335] | 77 | enddo
|
|---|
| 78 | enddo
|
|---|
| 79 | enddo
|
|---|
| 80 | enddo
|
|---|
| 81 | enddo
|
|---|
| 82 |
|
|---|
| 83 | eysmi = eyeli + eyvwi + eyhbi
|
|---|
| 84 | eninteract = eysmi
|
|---|
| [32289cd] | 85 | return
|
|---|
| 86 | end
|
|---|
