real*8 function enyshe_simple()

    include "INCL.H"

    integer i,j,numberOfAtoms,indexOfMolecule, i14
    integer indexOfFirstResidue, indexOfLastResidue
    integer indexOfFirstAtom, indexOfLastAtom, atomTypeI, indexOfAtom
    integer typeOfDihedral, classOfPotential
    integer indexOfFirstDihedral, indexOfLastDihedral
    real*8 totalEnergy, coulombEnergy, rij2, rij4, rij6, rij12, E0
    real*8 vdWEnergy, vdWDiff, torsionEnergy,xatI, yatI, zatI, cgI, tmpE
    real*8 hydrogenBondEnergy

    indexOfMolecule = 1
    indexOfFirstResidue = irsml1(indexOfMolecule)
    indexOfLastResidue = irsml2(indexOfMolecule)
    indexOfFirstAtom = iatrs1(indexOfFirstResidue)
    indexOfLastAtom = iatrs2(indexOfLastResidue)
    numberOfAtoms = indexOfLastAtom - indexOfFirstAtom
!$acc region
    do i = indexOfFirstAtom, indexOfLastAtom
        xatI = xat(i)
        yatI = yat(i)
        zatI = zat(i)
        cgI = cgat(i)
        atomTypeI = ityat(i)
        
        do j = i + 1, indexOfLastAtom
            rij2 = (xatI - xat(j)) * (xatI - xat(j)) + (yatI - yat(j)) * (yatI - yat(j)) &
                 + (zatI - zat(j)) * (zatI - zat(j))
            rij4 = rij2 * rij2
            rij6 = rij4 * rij2

            ! Calculate Coulomb energy
            coulombEnergy = coulombEnergy + cgI * cgat(j) / sqrt(rij2)
            ! Calculate van-der-Waals and hydrogen bond energy
            vdWEnergy = vdWEnergy + aij(atomTypeI, ityat(j)) / (rij6 * rij6) - cij(atomTypeI, ityat(j)) / rij6
            ! Calculate hydrogen bond energy
            hydrogenBondEnergy = hydrogenBondEnergy + ahb(ityat(i),ityat(j))/(rij6*rij6) &
                                              - chb(ityat(i),ityat(j))/(rij6*rij4)
        end do

        do i14=i14at1(i),i14at2(i)
            j = l14at(i14)
            rij2 = (xat(i) - xat(j)) * (xat(i) - xat(j)) + (yat(i) - yat(j)) * (yat(i) - yat(j)) &
                + (zat(i) - zat(j)) * (zat(i) - zat(j))
            rij4 = rij2 * rij2
            rij12 = rij4 * rij4 * rij4

            vdWDiff = vdWDiff - 0.5 * aij(ityat(i), ityat(j)) / (rij12)
        end do

    end do

    indexOfFirstDihedral = ivrml1(indexOfMolecule)
    indexOfLastDihedral = indexOfFirstDihedral + nvrml(indexOfMolecule)

    do i = indexOfFirstDihedral, indexOfLastDihedral
        indexOfAtom = iatvr(i)
        typeOfDihedral = ityvr(i)
        classOfPotential = iclvr(i)
        if (typeOfDihedral == 3) then
            if (e0to(classOfPotential) == 0.0) then
                tmpE = torsionEnergy + e0to(classOfPotential) * (1.0 + sgto(classOfPotential) &
                         * cos(toat(indexOfAtom) * rnto(classOfPotential)))
            endif
        endif
        
        torsionEnergy = tmpE
        end do

        totalEnergy = coulombEnergy + vdWEnergy + vdWDiff + torsionEnergy

        enyshe_simple = totalEnergy
!$acc end region
        return

end function enyshe_simple