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