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 |
---|