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 | 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 |
|
---|
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)
|
---|
33 | ! Atom class of current atom
|
---|
34 | ity=ityat(iat)
|
---|
35 | ! Point charge at current atom
|
---|
36 | cqi=conv*cgat(iat)
|
---|
37 | ! Cartesian coordinates of current atom
|
---|
38 | xi=xat(iat)
|
---|
39 | yi=yat(iat)
|
---|
40 | zi=zat(iat)
|
---|
41 |
|
---|
42 | do jat = iatrs1(jres), iatrs2(jres)
|
---|
43 | ! Atom type of partner
|
---|
44 | jty=ityat(jat)
|
---|
45 | ! Differences in cartesian coordinates
|
---|
46 | xj=xat(jat)
|
---|
47 | yj=yat(jat)
|
---|
48 | zj=zat(jat)
|
---|
49 |
|
---|
50 | xij=xat(jat)-xi
|
---|
51 | yij=yat(jat)-yi
|
---|
52 | zij=zat(jat)-zi
|
---|
53 | ! Cartesian distance and higher powers
|
---|
54 | rij2=xij*xij+yij*yij+zij*zij
|
---|
55 | rij4=rij2*rij2
|
---|
56 | rij6=rij4*rij2
|
---|
57 | rij=sqrt(rij2)
|
---|
58 | ! Are we using a distance dependent dielectric constant?
|
---|
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
|
---|
65 | ! Coulomb interaction
|
---|
66 | eyeli=eyeli+cqi*cgat(jat)/(rij*ep)
|
---|
67 | ! If the two atoms cannot form a hydrogen bond use 6-12 Lennard-Jones potential
|
---|
68 | if (ihbty(ity,jty).eq.0) then
|
---|
69 | eyvwi=eyvwi+aij(ity,jty)/(rij6*rij6)
|
---|
70 | & -cij(ity,jty)/rij6
|
---|
71 | else
|
---|
72 | ! For hydrogen bonding use 10-12 Lennard-Jones potential
|
---|
73 | eyhbi=eyhbi+ahb(ity,jty)/(rij6*rij6)
|
---|
74 | & -chb(ity,jty)/(rij6*rij4)
|
---|
75 | endif
|
---|
76 | enddo
|
---|
77 | enddo
|
---|
78 | enddo
|
---|
79 | enddo
|
---|
80 | enddo
|
---|
81 | enddo
|
---|
82 |
|
---|
83 | eysmi = eyeli + eyvwi + eyhbi
|
---|
84 | eninteract = eysmi
|
---|
85 | return
|
---|
86 | end
|
---|