Index: Makefile
===================================================================
--- Makefile	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
+++ Makefile	(revision 7137e5d1b97bd3e5b69aa77ddb1a56f2cbf42453)
@@ -17,9 +17,12 @@
 	$(F90) $(F_FLAGS) $<
 
+.SUFFIXES:	.o .F90
+.F90.o:
+	$(F90) $(F_FLAGS) $<
 # ==================================== Variables for compiling and linking
 L_FLAGS=-g
 # Linker flags
 # Default flags, e.g., for gfortran, g77
-#L_FLAGS=-O2 -g
+L_FLAGS=-O2 -g
 
 # Flags for Intel's ifort
@@ -37,10 +40,10 @@
 # Compiler flags
 # Default flags, e.g., for gfortran, g77
-#F_FLAGS=-c -O2 -g
+F_FLAGS=-c -O2 -g
 # Flags for Intel's ifort
 #F_FLAGS=-c -O3 -axT -g
 #F_FLAGS=-c -fast -parallel
 # Flags for debugging build
- F_FLAGS=-c -O0 -g
+# F_FLAGS=-c -O0 -g
 # FLAGS for Cray
 # F_FLAGS=-c -fastsse -Mipa
@@ -52,8 +55,8 @@
 
 # Select your compiler
-F90=ifort
+#F90=ifort
 #F90=pgf90
 #F90=xlf
-#F90=gfortran
+F90=gfortran
 MPIF90=mpif90
 
@@ -69,5 +72,5 @@
 init_energy.o init_molecule.o gradient.o energy.o \
 regul.o nursvr.o twister.o eninteract.o bgs.o eyabgn.o enylun.o utilities.o\
-partem_s.o esolan.o opesol.o init_lund.o
+partem_s.o esolan.o opesol.o init_lund.o enyshe_simple.o timer.o
 #_________Serial implementation of ECEPP/3 and solvent
 SOBJ=enyshe.o enysol.o
@@ -92,4 +95,5 @@
 # metropolis
 # ============================================== Linking
+all: $(PROG)
 
 $(PROG):	$(OBJ) $(SOBJ) main.o
Index: enyshe_simple.f90
===================================================================
--- enyshe_simple.f90	(revision 7137e5d1b97bd3e5b69aa77ddb1a56f2cbf42453)
+++ enyshe_simple.f90	(revision 7137e5d1b97bd3e5b69aa77ddb1a56f2cbf42453)
@@ -0,0 +1,78 @@
+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
Index: timer.F90
===================================================================
--- timer.F90	(revision 7137e5d1b97bd3e5b69aa77ddb1a56f2cbf42453)
+++ timer.F90	(revision 7137e5d1b97bd3e5b69aa77ddb1a56f2cbf42453)
@@ -0,0 +1,141 @@
+! timer.f90
+!
+! (c)2007 Jan H. Meinke
+!
+#ifdef BGL
+module rts
+! interface definition of rts_get_timebase, Fortran2003, available in XLF Version 9
+  interface
+    function rts_get_timebase() bind(c)
+      use, intrinsic       :: iso_c_binding
+      integer(c_long_long) :: rts_get_timebase
+    end function rts_get_timebase
+  end interface
+end module rts
+#endif
+
+!! This module contains a minimal set of routines to time sections of code.
+!  It can keep track of several timers at once. Each timer is associated with
+!  an integer ID. Individual timing events are not stored.
+!  
+!  Currently only timerIDs between 1 and 10 are valid.
+module timer
+#ifdef BGL
+    use rts
+    use, intrinsic :: iso_c_binding
+#endif
+    implicit none
+    
+    !! Contains the minimum time ever measured, the maximum time ever measured,
+    !  the average time, the standard deviation, and the number of measurements.
+    type TimingStructure
+        integer :: counter
+#ifdef BGL
+        integer(c_long_long) :: t0, t1
+#else
+        real(kind=8) :: t0, t1
+#endif
+        real(kind=8) :: minTime, maxTime
+        real(kind=8) :: sumOfTimes, sumOfSquares
+        real(kind=8) :: average, sigma
+    end type TimingStructure
+            
+    !! Contains the timing data for all timers.
+    type(TimingStructure), allocatable :: timingData(:)
+            
+    integer :: initialTimerCount = 10
+#ifdef BGL
+    real(kind=8) :: clockspeed = 700.0d6
+#endif
+            
+    contains    
+        !! Initializes timer
+        !  This routine must be called before starting any timer.
+        !  
+        !  @author Jan H. Meinke
+        subroutine init_timer()
+            integer :: i
+            
+            ! Initialize timers
+            allocate (timingData(initialTimerCount))
+            forall(i = 1:initialTimerCount)
+                timingData(i) = TimingStructure(0, -1.0, -1.0, -1.0, -1.0, 0.0, &
+                        0.0, -1.0, -1.0)
+            end forall
+        end subroutine init_timer
+    
+        !! Start a new measurement. 
+        !
+        !  @param id ID of the timer
+        !
+        !  @author Jan H. Meinke
+        subroutine start_timer(id)
+            integer, intent(in) :: id
+            timingData(id)%t0 = getCurrentTime()
+        end subroutine start_timer
+                
+        !! Take a measurement and stop the timer with ID id
+        !
+        !  @param id ID of the timer
+        !
+        !  @author Jan H. Meinke
+        subroutine stop_timer(id)
+            integer, intent(in) :: id
+            real(kind=8) :: dt
+                    
+            timingData(id)%t1 = getCurrentTime()
+#ifdef BGL
+            dt = (timingData(id)%t1 - timingData(id)%t0) / clockspeed
+#else 
+            dt = (timingData(id)%t1 - timingData(id)%t0)
+#endif
+            timingData(id)%counter = timingData(id)%counter + 1
+            timingData(id)%maxTime = max(dt, timingData(id)%maxTime)
+            if (timingData(id)%minTime > 0) then
+                timingData(id)%minTime = min(dt, timingData(id)%minTime)
+            else 
+                timingData(id)%minTime = dt
+            end if
+            timingData(id)%sumOfTimes = timingData(id)%sumOfTimes + dt
+            timingData(id)%sumOfSquares = timingData(id)%sumOfSquares + dt ** 2
+                    
+        end subroutine stop_timer
+              
+        subroutine evaluate(id)
+            integer, intent(in) :: id
+            if (timingData(id)%counter > 0) then
+                timingData(id)%average = timingData(id)%sumOfTimes / timingData(id)%counter
+            end if
+            if (timingData(id)%counter > 1) then
+!                write (*,*) timingData(id)%average, timingData(id)%sumOfSquares,&
+!                        timingData(id)%counter
+                timingData(id)%sigma = sqrt((abs(timingData(id)%average ** 2 &
+                        - timingData(id)%sumOfSquares / timingData(id)%counter)) &
+                        / (timingData(id)%counter - 1))
+            end if
+        end subroutine evaluate
+        
+        !! Wrapper around the timing function used.
+        !  If we are using MPI, we can use MPI_WTime, otherwise we use cpu_time(time)
+#ifdef BGL
+        integer(c_long_long) function getCurrentTime()
+            integer(c_long_long) :: time
+            time = rts_get_timebase()
+            getCurrentTime = time
+        end function getCurrentTime
+#elif MPI
+        real(kind=8) function getCurrentTime()
+            real(kind=8) :: time
+            include 'mpif.h'
+            time = MPI_Wtime()
+            getCurrentTime = time
+        end function getCurrentTime
+#else
+        real(kind=8) function getCurrentTime()
+            real(kind=8) :: time
+            call cpu_time(time)
+            getCurrentTime = time
+        end function getCurrentTime
+#endif 
+        
+end module timer
