1 | real*8 function enyshe_simple()
|
---|
2 |
|
---|
3 | include "INCL.H"
|
---|
4 |
|
---|
5 | integer i,j,numberOfAtoms,indexOfMolecule, i14
|
---|
6 | integer indexOfFirstResidue, indexOfLastResidue
|
---|
7 | integer indexOfFirstAtom, indexOfLastAtom, atomTypeI, indexOfAtom
|
---|
8 | integer typeOfDihedral, classOfPotential
|
---|
9 | integer indexOfFirstDihedral, indexOfLastDihedral
|
---|
10 | real*8 totalEnergy, coulombEnergy, rij2, rij4, rij6, rij12, E0
|
---|
11 | real*8 vdWEnergy, vdWDiff, torsionEnergy,xatI, yatI, zatI, cgI, tmpE
|
---|
12 | real*8 hydrogenBondEnergy
|
---|
13 |
|
---|
14 | indexOfMolecule = 1
|
---|
15 | indexOfFirstResidue = irsml1(indexOfMolecule)
|
---|
16 | indexOfLastResidue = irsml2(indexOfMolecule)
|
---|
17 | indexOfFirstAtom = iatrs1(indexOfFirstResidue)
|
---|
18 | indexOfLastAtom = iatrs2(indexOfLastResidue)
|
---|
19 | numberOfAtoms = indexOfLastAtom - indexOfFirstAtom
|
---|
20 | !$acc region
|
---|
21 | do i = indexOfFirstAtom, indexOfLastAtom
|
---|
22 | xatI = xat(i)
|
---|
23 | yatI = yat(i)
|
---|
24 | zatI = zat(i)
|
---|
25 | cgI = cgat(i)
|
---|
26 | atomTypeI = ityat(i)
|
---|
27 |
|
---|
28 | do j = i + 1, indexOfLastAtom
|
---|
29 | rij2 = (xatI - xat(j)) * (xatI - xat(j)) + (yatI - yat(j)) * (yatI - yat(j)) &
|
---|
30 | + (zatI - zat(j)) * (zatI - zat(j))
|
---|
31 | rij4 = rij2 * rij2
|
---|
32 | rij6 = rij4 * rij2
|
---|
33 |
|
---|
34 | ! Calculate Coulomb energy
|
---|
35 | coulombEnergy = coulombEnergy + cgI * cgat(j) / sqrt(rij2)
|
---|
36 | ! Calculate van-der-Waals and hydrogen bond energy
|
---|
37 | vdWEnergy = vdWEnergy + aij(atomTypeI, ityat(j)) / (rij6 * rij6) - cij(atomTypeI, ityat(j)) / rij6
|
---|
38 | ! Calculate hydrogen bond energy
|
---|
39 | hydrogenBondEnergy = hydrogenBondEnergy + ahb(ityat(i),ityat(j))/(rij6*rij6) &
|
---|
40 | - chb(ityat(i),ityat(j))/(rij6*rij4)
|
---|
41 | end do
|
---|
42 |
|
---|
43 | do i14=i14at1(i),i14at2(i)
|
---|
44 | j = l14at(i14)
|
---|
45 | rij2 = (xat(i) - xat(j)) * (xat(i) - xat(j)) + (yat(i) - yat(j)) * (yat(i) - yat(j)) &
|
---|
46 | + (zat(i) - zat(j)) * (zat(i) - zat(j))
|
---|
47 | rij4 = rij2 * rij2
|
---|
48 | rij12 = rij4 * rij4 * rij4
|
---|
49 |
|
---|
50 | vdWDiff = vdWDiff - 0.5 * aij(ityat(i), ityat(j)) / (rij12)
|
---|
51 | end do
|
---|
52 |
|
---|
53 | end do
|
---|
54 |
|
---|
55 | indexOfFirstDihedral = ivrml1(indexOfMolecule)
|
---|
56 | indexOfLastDihedral = indexOfFirstDihedral + nvrml(indexOfMolecule)
|
---|
57 |
|
---|
58 | do i = indexOfFirstDihedral, indexOfLastDihedral
|
---|
59 | indexOfAtom = iatvr(i)
|
---|
60 | typeOfDihedral = ityvr(i)
|
---|
61 | classOfPotential = iclvr(i)
|
---|
62 | if (typeOfDihedral == 3) then
|
---|
63 | if (e0to(classOfPotential) == 0.0) then
|
---|
64 | tmpE = torsionEnergy + e0to(classOfPotential) * (1.0 + sgto(classOfPotential) &
|
---|
65 | * cos(toat(indexOfAtom) * rnto(classOfPotential)))
|
---|
66 | endif
|
---|
67 | endif
|
---|
68 |
|
---|
69 | torsionEnergy = tmpE
|
---|
70 | end do
|
---|
71 |
|
---|
72 | totalEnergy = coulombEnergy + vdWEnergy + vdWDiff + torsionEnergy
|
---|
73 |
|
---|
74 | enyshe_simple = totalEnergy
|
---|
75 | !$acc end region
|
---|
76 | return
|
---|
77 |
|
---|
78 | end function enyshe_simple
|
---|