Index: INCL.H
===================================================================
--- INCL.H	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ INCL.H	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -1,5 +1,4 @@
-!      implicit none
-      implicit integer*4 (i-n)
-      implicit real*8 (a-h,o-z)
+      implicit none
+      
       character*255 version
 
Index: bgs.f
===================================================================
--- bgs.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ bgs.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -2,5 +2,5 @@
 ! Trial version implementing the semi-local conformational update
 ! BGS (Biased Gaussian Steps). This file presently contains the
-! functions initlund, bgsposs and bgs. 
+! functions initlund, bgsposs and bgs.
 !
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
@@ -9,9 +9,11 @@
 
 
-! Checks if it is possible to perform a BGS update starting at the 
+! Checks if it is possible to perform a BGS update starting at the
 ! variable indexed ipos. Calls: none.
       logical function bgsposs(ips)
       include 'INCL.H'
       include 'incl_lund.h'
+      integer jv, ips, iaa, nursvr, nnonfx, i
+
       logical ians
 
@@ -20,20 +22,20 @@
       ians=.true.
 !      print *,'evaluating bgs possibility for ',ips,nmvr(jv)
-      if (nmvr(jv).ne.'phi') then 
+      if (nmvr(jv).ne.'phi') then
 !         print *,'bgs not possible because variable name is ',nmvr(jv)
          ians=.false.
-      else if (iaa.gt.(irsml2(mlvr(jv))-3)) then 
+      else if (iaa.gt.(irsml2(mlvr(jv))-3)) then
 !         print *,'bgs impossible, residue too close to end'
 !         print *,'iaa = ',iaa,' end = ',irsml2(mlvr(jv))
          ians=.false.
-      else 
-         nnonfx=0 
+      else
+         nnonfx=0
          do i=iaa,iaa+3
-            if (iphi(i).gt.0) then 
+            if (iphi(i).gt.0) then
                if (.not.fxvr(iphi(i))) then
                   nnonfx=nnonfx+1
                endif
             endif
-            if (.not.fxvr(ipsi(i))) then 
+            if (.not.fxvr(ipsi(i))) then
                nnonfx=nnonfx+1
             endif
@@ -52,15 +54,15 @@
 
 ! Biased Gaussian Steps. Implements a semi-local conformational update
-! which modifies the protein backbone locally in a certain range of 
-! amino acids. The 'down-stream' parts of the molecule outside the 
+! which modifies the protein backbone locally in a certain range of
+! amino acids. The 'down-stream' parts of the molecule outside the
 ! region of update get small rigid body translations and rotations.
 !
-! Use the update sparingly. It is rather local, and is not of great 
+! Use the update sparingly. It is rather local, and is not of great
 ! value if we need big changes in the conformation. It is recommended
-! that this update be used to refine a structure around a low energy 
+! that this update be used to refine a structure around a low energy
 ! candidate structure. Even at low energies, if you always
-! perform BGS, the chances of coming out of that minimum are small. 
-! So, there is a probability bgsprob, which decides whether BGS or the 
-! normal single angle update is used. 
+! perform BGS, the chances of coming out of that minimum are small.
+! So, there is a probability bgsprob, which decides whether BGS or the
+! normal single angle update is used.
 !
 ! Calls: energy, dummy (function provided as argument), addang, (rand)
@@ -70,4 +72,10 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision grnd, xiv, bv, ab, rv, dv, a, sum, p, r1, r2
+      double precision ppsi, wfw, addang, enw, energy, wbw, rd, delta
+      double precision dummy, eol1
+
+      integer ivar, i, nph, jv, ia, nursvr, icurraa, j, k, l
+
       external dummy
       dimension xiv(8,3),bv(8,3),rv(3,3),dv(3,8,3)
@@ -76,5 +84,5 @@
 ! Initialize
 !      print *,'using BGS on angle ',nmvr(idvr(ivar))
-      if (bgsnvar.eq.0) then 
+      if (bgsnvar.eq.0) then
          bgs=0
          goto 171
@@ -91,5 +99,5 @@
       do i=1,4
          icurraa=ia+i-1
-         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 
+         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
             nph=nph+1
             xiv(nph,1)=xat(iCa(icurraa))
@@ -102,5 +110,5 @@
      &           +bv(nph,3)*bv(nph,3)
             iph(nph)=iphi(icurraa)
-         endif 
+         endif
          if (.not.fxvr(ipsi(icurraa))) then
             nph=nph+1
@@ -136,5 +144,5 @@
          enddo
       enddo
-      do i=1,nph 
+      do i=1,nph
          do j=i,nph
             A(i,j)=0
@@ -157,7 +165,7 @@
                sum=sum-A(i,k)*A(j,k)
             enddo
-            if (i.eq.j) then 
-               p(i)=sqrt(sum) 
-            else 
+            if (i.eq.j) then
+               p(i)=sqrt(sum)
+            else
                A(j,i)=sum/p(i)
             endif
@@ -177,5 +185,5 @@
          dph(i)=0
       enddo
-! Solve lower triangular matrix to get dphi proposals 
+! Solve lower triangular matrix to get dphi proposals
       do i=nph,1,-1
          sum=ppsi(i)
@@ -190,5 +198,5 @@
       do i=1,nph
          sum=sum+ppsi(i)*ppsi(i)
-      enddo 
+      enddo
       wfw=exp(-sum)
       do i=1,nph
@@ -208,5 +216,5 @@
       do i=1,4
          icurraa=ia+i-1
-         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 
+         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
             nph=nph+1
             xiv(nph,1)=xat(iCa(icurraa))
@@ -219,6 +227,6 @@
      &           +bv(nph,3)*bv(nph,3)
             iph(nph)=iphi(icurraa)
-         endif 
-         if (.not.fxvr(ipsi(icurraa))) then 
+         endif
+         if (.not.fxvr(ipsi(icurraa))) then
             nph=nph+1
             xiv(nph,1)=xat(iC(icurraa))
@@ -253,5 +261,5 @@
          enddo
       enddo
-      do i=1,nph 
+      do i=1,nph
          do j=i,nph
             A(i,j)=0
@@ -262,5 +270,5 @@
             enddo
             A(i,j)=bbgs*A(i,j)
-            if (i.eq.j) then 
+            if (i.eq.j) then
                A(i,j)=A(i,j)+1
             endif
@@ -274,7 +282,7 @@
                sum=sum-A(i,k)*A(j,k)
             enddo
-            if (i.eq.j) then 
-               p(i)=sqrt(sum) 
-            else 
+            if (i.eq.j) then
+               p(i)=sqrt(sum)
+            else
                A(j,i)=sum/p(i)
             endif
@@ -285,5 +293,5 @@
          do j=i+1,nph
             ppsi(i)=ppsi(i)+A(j,i)*dph(j)
-         enddo  
+         enddo
       enddo
       sum=0
@@ -291,5 +299,5 @@
          sum=sum+ppsi(i)*ppsi(i)
       enddo
-      wbw=exp(-sum) 
+      wbw=exp(-sum)
       do i=1,nph
          wbw=wbw*p(i)
@@ -306,13 +314,13 @@
 !       call outpdb(0,'after.pdb')
 !      print *,'after outpdb for after.pdb'
-!      do i=1,nph 
+!      do i=1,nph
 !         print *,'BGS>',i,iph(i),vlvr(iph(i)),dph(i)
 !      enddo
-      if (rd.ge.delta) then 
+      if (rd.ge.delta) then
 !     accept
          eol1=enw
          bgs=1
 !         print *,'BGS move accepted'
-      else 
+      else
 !     reject
          vlvr = ovr
Index: cnteny.f
===================================================================
--- cnteny.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ cnteny.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2005       Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -20,5 +20,12 @@
       include 'INCL.H'
 
-      parameter (coeycn=2.d0)  ! min. cont. energy to display 
+      double precision coeycn, coey, eyatcn, cqi, xi, yi, zi, xij, yij
+      double precision zij, rij2, rij4, rij6, rij, sr, ep, ey
+
+      integer ieyel, ntlvr, nml, iat1, nat, i, ifivr, i1s, io, iv, ia
+      integer it, ic, i2s, ims, i1, i2, ii, ity, ivw, j, jj, jty, i14
+      integer nbc, ir, nursat, idxat
+
+      parameter (coeycn=2.d0)  ! min. cont. energy to display
 
       dimension eyatcn(mxat),idxat(mxat)
Index: difang.f
===================================================================
--- difang.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ difang.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -21,6 +21,6 @@
 ! ......................................................
 
-      implicit real*8 (a-h,o-z)
-
+      implicit none
+      double precision pi, pi2, d, a2, a1
       parameter (pi=3.141592653589793d0,
      &           pi2=2.d0*pi)
@@ -47,5 +47,6 @@
 ! ......................................................
 
-      implicit real*8 (a-h,o-z)
+      implicit none
+      double precision pi, pi2, d, a1, a2
 
       parameter (pi=3.141592653589793d0,
Index: energy.f
===================================================================
--- energy.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ energy.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -22,6 +22,11 @@
 
       include 'INCL.H'
+      double precision esm, teysl, enyshe, enyflx, enylun, enyreg
+      double precision enysol, esolan, exvlun, eyabgn, eninteract
+
+      integer i
+
       double precision teysm, teyel, teyvw, teyhb, teyvr
-!      print *,'energy function with ientyp  = ',ientyp   
+!      print *,'energy function with ientyp  = ',ientyp
       esm = 0.d0
       teysm = 0.d0
@@ -32,5 +37,5 @@
       teysl = 0.d0
 
-      do i = 1,ntlml  
+      do i = 1,ntlml
          eysm=0
          eyel=0
@@ -44,5 +49,5 @@
         if (ientyp.eq.0.or.ientyp.eq.3) then
           esm=esm+enyshe(i)
-        else if (ientyp.eq.1) then 
+        else if (ientyp.eq.1) then
           esm=esm+enyflx(i)
         else if (ientyp.eq.2) then
@@ -57,15 +62,15 @@
 !     The Lund term stores the hydrophobicity energy in eysl
            teysl = teysl + eysl
-        else 
+        else
 !     .. and the excluded volume term in eyvw, which is calculated once.
            teyvw = teyvw + eyvw
         endif
 
-        if (ireg.eq.1)  eyrg=enyreg(i)        
+        if (ireg.eq.1)  eyrg=enyreg(i)
 
       enddo
 
-      if (ientyp.ne.2) then 
-!     Don't touch eysl if using Lund potential, as enylun stores 
+      if (ientyp.ne.2) then
+!     Don't touch eysl if using Lund potential, as enylun stores
 !     its hydrophobicity term there.
          if (itysol.gt.0) then
@@ -79,5 +84,5 @@
             eysl=0.d0
          endif
-      else 
+      else
 !     Add excluded volume term and save it in eyvw
          esm=esm+exvlun(0)
@@ -85,10 +90,10 @@
       endif
 
-! The Abagyan entropic corrections depend on the area exposed to the 
+! The Abagyan entropic corrections depend on the area exposed to the
 ! solvent for each residue. So, this term has to be evaluated after the
 ! solvent term.
       eyab=0.0
       if (ientyp.eq.3) then
-         do i = 1,ntlml  
+         do i = 1,ntlml
             eyab=eyab+eyabgn(i)
          enddo
@@ -103,5 +108,5 @@
       eyvr = teyvr
       eysl = teysl
-      
+
       if (ientyp.ne.2) then
 !     This is temporary. eninteract() does not yet know how to calculate
@@ -116,5 +121,5 @@
 !c Calculates the internal energy for a single molecule.
 !  All the partial energies are thus set to their values for molecule
-!  nml. 
+!  nml.
 !
 !  @param nml the ID of the molecule
@@ -125,9 +130,14 @@
 
 !f2py intent(in) nml
-      
+
       include 'INCL.H'
+      double precision esm, enyshe, enyflx, enylun, enyreg, enysol
+      double precision esolan, exvlun, eyabgn
+
+      integer nml, i
+
       esm = 0.d0
 
-      call setvar(nml,vlvr)  
+      call setvar(nml,vlvr)
 
       if (ientyp.eq.0.or.ientyp.eq.3) then
@@ -149,8 +159,8 @@
             eysl=0.d0
          endif
-      else 
+      else
          esm=esm+exvlun(nml)
       endif
-! The Abagyan entropic corrections depend on the area exposed to the 
+! The Abagyan entropic corrections depend on the area exposed to the
 ! solvent for each residue. So, this term has to be evaluated after the
 ! solvent term.
Index: eninteract.f
===================================================================
--- eninteract.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ eninteract.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -15,5 +15,10 @@
 
         include 'INCL.H'
-                
+
+      double precision cqi, xi, yi, zi, xj, yj, zj, xij, yij, zij, rij2
+      double precision rij4, rij6, rij, sr, ep
+
+      integer iml, jml, ires, jres, iat, ity, jat, jty
+
         eysmi=0.0
         eyeli=0.0
@@ -42,5 +47,5 @@
                     yj=yat(jat)
                     zj=zat(jat)
-                    
+
                     xij=xat(jat)-xi
                     yij=yat(jat)-yi
@@ -69,5 +74,5 @@
      &                          -chb(ity,jty)/(rij6*rij4)
                     endif
-                  enddo 
+                  enddo
                 enddo
               enddo
@@ -78,4 +83,4 @@
         eysmi = eyeli + eyvwi + eyhbi
         eninteract = eysmi
-        return 
-      end   
+        return
+      end
Index: enylun.f
===================================================================
--- enylun.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enylun.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -1,8 +1,8 @@
 ! *******************************************************************
 ! SMMP version of Anders Irback's force field, to be called the Lund
-! force field. This file contains the function enylun, which in turn 
-! calls all the terms in the energy function. The terms Bias (ebias), 
-! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the 
-! Excluded volume (eexvol and eloexv) are also implemented in this 
+! force field. This file contains the function enylun, which in turn
+! calls all the terms in the energy function. The terms Bias (ebias),
+! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the
+! Excluded volume (eexvol and eloexv) are also implemented in this
 ! file.
 !
@@ -13,4 +13,6 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      integer iat1, iat2
+
           integer ires,frst,lst,root
           do iat1=iCa(ires)+frst,iCa(ires)+lst
@@ -30,4 +32,9 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision eunit, csacc
+
+      integer i, j, iml, iat1, iat2, k, iat3, l, iat4, m, iat3p, irs
+      integer iatoff, iatmrg, ilp, lci
+
       character mynm*4
       logical prlvr
@@ -83,9 +90,9 @@
             prlvr=.true.        ! residue i is a proline variant
             print *, 'proline variant ',mynm,i
-         else 
+         else
             prlvr=.false.
-         endif 
-
-         if (mynm.eq.'ala') then 
+         endif
+
+         if (mynm.eq.'ala') then
             nhpat(i)=1
             ihpat(i,1)=iCa(i)+2
@@ -193,11 +200,11 @@
 
 !     Initialization of the connections matrix matcon(i,j). The index
-!     i runs from -mxconr to +mxconr, and j from 1 to mxat. 
+!     i runs from -mxconr to +mxconr, and j from 1 to mxat.
 !     matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed
 !                      = 2, if atoms i1 and i2 are separated by 3 covalent
 !                           bonds and their distance can change
 !                      = 1, for all other pairs
-!     if abs(i2-i1) > mxconr, the atoms are assumed to be separated by 
-!     many bonds, and with no restriction on their distances. On a protein 
+!     if abs(i2-i1) > mxconr, the atoms are assumed to be separated by
+!     many bonds, and with no restriction on their distances. On a protein
 !     molecule made of natural amino acids, atoms with indices separated
 !     by more than 35 can not be connected by three covalent bonds.
@@ -228,5 +235,5 @@
                      do m=1,nbdat(iat2)
                         iat3p=ibdat(m,iat2)
-                        if (iat3p.ne.iat3) then 
+                        if (iat3p.ne.iat3) then
                            matcon(iat4-iat3p,iat3p)=2 ! 3 covalent bonds
                            matcon(iat3p-iat4,iat4)=2 !
@@ -255,5 +262,5 @@
                 matcon(iat2-iat1,iat1)=0
                 matcon(iat1-iat2,iat2)=0
-            enddo 
+            enddo
             matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
             matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
@@ -266,10 +273,10 @@
 !     Below: for certain residues, some atoms separated by 3 or more bonds
 !     do not change distance. So, the connection matrix term for such pairs
-!     should be zero. 
+!     should be zero.
 
          do irs=irsml1(iml),irsml2(iml)
             if (irs.eq.irsml1(iml)) then
                iatoff=1
-            else 
+            else
                iatoff=0
             endif
@@ -283,10 +290,10 @@
             if ((mynm.eq.'pro')) then
                prlvr=.true.     ! residue i is a proline variant
-            else 
+            else
                prlvr=.false.
-            endif 
-            if ((mynm.eq.'asn')) then 
+            endif
+            if ((mynm.eq.'asn')) then
                 call set_local_constr(irs,5,9,2)
-            else if ((mynm.eq.'gln')) then 
+            else if ((mynm.eq.'gln')) then
                 call set_local_constr(irs,8,12,5)
             else if ((mynm.eq.'arg')) then
@@ -294,5 +301,5 @@
             else if ((mynm.eq.'his')) then
                 call set_local_constr(irs,5,12,2)
-            else if (mynm.eq.'phe') then 
+            else if (mynm.eq.'phe') then
                 call set_local_constr(irs,5,15,2)
             else if (mynm.eq.'tyr') then
@@ -308,5 +315,5 @@
                 call set_local_constr(irs,5,19,2)
             else if (prlvr) then
-!           Proline. Many more distances are fixed because of the fixed 
+!           Proline. Many more distances are fixed because of the fixed
 !           phi angle
                 call set_local_constr(irs,-3,11,-3)
@@ -347,10 +354,9 @@
          do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
             do iat2=iat1+1,iatrs2(irsml2(iml))
-               if (iat2-iat1.le.mxconr) then 
-                  if (matcon(iat2-iat1,iat1).eq.2) then
-                    ilp=ilp+1
-                    lcp1(ilp)=iat1
-                    lcp2(ilp)=iat2
-                  endif
+               if ((iat2-iat1.le.mxconr).and.
+     &                 matcon(iat2-iat1,iat1).eq.2) then
+                  ilp=ilp+1
+                  lcp1(ilp)=iat1
+                  lcp2(ilp)=iat2
                endif
             enddo
@@ -358,5 +364,5 @@
 
          ilpnd(iml)=ilp
-         if (iml.lt.ntlml) then 
+         if (iml.lt.ntlml) then
             ilpst(iml+1)=ilp+1
          endif
@@ -369,4 +375,5 @@
          enddo
       enddo
+
       print *,'finished initializing Lund force field'
       end
@@ -374,4 +381,6 @@
       integer function ihptype(iaa)
       include 'INCL.H'
+      integer iaa, ityp
+
       character mynm*4
       mynm=seq(iaa)
@@ -383,10 +392,10 @@
      &        .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then
          ityp=2
-      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp')) 
+      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
      &        then
          ityp=3
       endif
       ihptype=ityp
-      return 
+      return
       end
 
@@ -394,7 +403,11 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision q1, q2, et, xij, yij, zij
+
+      integer i, iat1, irsd, j, iat2
+
       dimension q1(2),q2(2)
       data q1/-0.2,0.2/
-      data q2/0.42,-0.42/ 
+      data q2/0.42,-0.42/
 
       et=0.0
@@ -412,9 +425,13 @@
       return
       end
-!     Evaluates backbone backbone hydrogen bond strength for residues 
+!     Evaluates backbone backbone hydrogen bond strength for residues
 !     i and j, taking the donor from residue i and acceptor from residue j
       real*8 function ehbmmrs(i,j)
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision rdon2, racc2, evlu
+
+      integer i, j
+
       double precision r2,r4,r6,dx,dy,dz,ca,cb
       integer d1,d2,a1,a2
@@ -427,5 +444,5 @@
       dz=zat(a1)-zat(d2)
       r2=dx*dx+dy*dy+dz*dz
-      if (r2.gt.cthb2) then 
+      if (r2.gt.cthb2) then
 !         print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2
 !         print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb
@@ -464,4 +481,10 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision ebiasrs, shbm1, shbm2, etmp, ehbmmrs, ehp, etmp1
+      double precision xij, yij, zij, r2, r6
+
+      integer istres, nml, indres, i, j, ihpi, ihptype, ihpj, i1, i2
+      integer iat1, iat2, iatt1, iatt2
+
       character mynm*4
       logical prlvr
@@ -485,7 +508,7 @@
      &           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
             prlvr=.true.        ! residue i is a proline variant
-         else 
+         else
             prlvr=.false.
-         endif 
+         endif
 
 !     Bias, or local electrostatic term. Excluded from the list are
@@ -502,17 +525,17 @@
 !     j run over the whole set of amino acids.
 !     No terms for residue i, if it is a proline variant.
-         if (.not.prlvr) then  
+         if (.not.prlvr) then
             do j=istres,indres
-               if ((j.lt.(i-2)).or.(j.gt.(i+1))) then 
+               if ((j.lt.(i-2)).or.(j.gt.(i+1))) then
                     shbm2=1.0
                     if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
                     etmp=ehbmmrs(i,j)
                     eyhb=eyhb+shbm1*shbm2*etmp
-                endif 
+                endif
             enddo
          endif
 !     Hydrophobicity, only if residue i is hydrophobic to start with
          ihpi=ihptype(i)
-         if (ihpi.ge.0) then 
+         if (ihpi.ge.0) then
 !        Unlike hydrogen bonds, the hydrophobicity potential is symmetric
 !        in i and j. So, the loop for j runs from i+1 to the end.
@@ -526,5 +549,5 @@
 !                  print *, 'hydrophobicity contribution ',etmp,i,j
                   eysl=eysl+etmp
-               endif 
+               endif
             enddo
          endif
@@ -535,5 +558,5 @@
 !     Local pair or third-neighbour excluded volume
 !     Numerically this is normally the term with largest positive
-!     contribution to the energy in an equilibrated stystem. 
+!     contribution to the energy in an equilibrated stystem.
 
       i1=ilpst(nml)
@@ -551,5 +574,5 @@
          r2=(xij*xij + yij*yij + zij*zij)
          if (r2.le.exvcut2) then
-            r6=sig2lcp(iatt1,iatt2)/r2 
+            r6=sig2lcp(iatt1,iatt2)/r2
             r6=r6*r6*r6
             etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2)
@@ -560,5 +583,5 @@
             etmp=etmp+etmp1
          endif
-!         print *,'pair : ',iat1,iat2,' contribution ',etmp1 
+!         print *,'pair : ',iat1,iat2,' contribution ',etmp1
 !         print *,exvcut2,r2
       enddo
@@ -579,4 +602,8 @@
       include 'INCL.H'
       include 'incl_lund.h'
+      double precision b2, a2, r2min, xij, yij, zij, dtmp, ssum
+
+      integer ihp1, ihp2, ni, i1, nj, i2, i, k, j, l
+
       dimension r2min(12)
 
@@ -608,8 +635,8 @@
       ssum=0
       do i=1,ni+nj
-         if (r2min(i).le.b2) then 
-            if (r2min(i).lt.a2) then 
+         if (r2min(i).le.b2) then
+            if (r2min(i).lt.a2) then
                ssum=ssum+1
-            else 
+            else
                ssum=ssum+(b2-r2min(i))/(b2-a2)
             endif
@@ -626,18 +653,26 @@
       include 'incl_lund.h'
 !     For multi-chain systems it makes little sense to split the calculation
-!     of this term into an 'interaction part' and a contribution from 
+!     of this term into an 'interaction part' and a contribution from
 !     individual molecules. So, normally this should always be called with
 !     argument nml=0. Only for diagnostic reasons, you might want to find
-!     the contribution from one molecule in a multi-chain system assuming 
+!     the contribution from one molecule in a multi-chain system assuming
 !     there was no other molecule.
+      double precision xmin, ymin, zmin, xmax, ymax, zmax, sizex, sizey
+      double precision sizez, shiftx, shifty, shiftz, etmp, xij, yij
+      double precision zij, r2, r6, etmp1
+
+      integer nml, istat, indat, i, incell, ndx, ndy, ndz, nxy, ncell
+      integer nocccl, locccl, j, mx, my, mz, icellj, icell, jj, isort
+      integer icl, lcell, ix, iz, iy, nex, ney, nez, nsame, nngbr, jx
+      integer jy, jz, jcl, ii, ngbr, i1, iat1, i2, iat2, iatt1, iatt2
+
       dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
       dimension icell(mxat)
       integer :: jx1, jy1, jz1
-      logical doit
-
-      if (nml.eq.0) then 
+
+      if (nml.eq.0) then
          istat=iatrs1(irsml1(1))
          indat=iatrs2(irsml2(ntlml))
-      else 
+      else
          istat=iatrs1(irsml1(nml))
          indat=iatrs2(irsml2(nml))
@@ -645,7 +680,7 @@
 
       eyvw=0.0
-!     The beginning part of this implementation is very similar to the 
-!     assignment of cells to the atoms during calculation of solvent 
-!     accessible surface area. So, much of that part is similar. But 
+!     The beginning part of this implementation is very similar to the
+!     assignment of cells to the atoms during calculation of solvent
+!     accessible surface area. So, much of that part is similar. But
 !     unlike the accessible surface calculations, this term is symmetric
 !     in any two participating atoms. So, the part after the assignment
@@ -724,9 +759,9 @@
          icellj=mx+my*ndx+mz*nxy+1
          icell(j)=icellj
-         if (icellj.gt.mxcell) then 
+         if (icellj.gt.mxcell) then
             print *,'exvlun> bad cell index ',icellj,' for atom ',j
             stop
-         else 
-            if (incell(icellj).eq.0) then 
+         else
+            if (incell(icellj).eq.0) then
 !           previously unoccupied cell
                nocccl=nocccl+1
@@ -756,5 +791,5 @@
          ix=mod(lcell-1,ndx)
          iz = (lcell - 1) / nxy
-         iy = ((lcell - 1) - iz * nxy) / ndx 
+         iy = ((lcell - 1) - iz * nxy) / ndx
 
 c         print *,'icl=',icl,'absolute index of cell = ',lcell
@@ -767,23 +802,23 @@
          nngbr=0
          do jx=ix,nex
-            if (jx.ge.ndx) then 
+            if (jx.ge.ndx) then
                jx1=0
-            else 
+            else
                jx1=jx
             endif
             do jy=iy,ney
-               if (jy.ge.ndy) then 
+               if (jy.ge.ndy) then
                   jy1=0
-               else 
+               else
                   jy1=jy
-               endif 
+               endif
                do jz=iz,nez
-                  if (jz.ge.ndz) then 
+                  if (jz.ge.ndz) then
                      jz1=0
-                  else 
+                  else
                      jz1=jz
                   endif
                   jcl=jx1 + ndx*jy1 + nxy*jz1 + 1
-!                  write(*,*)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
+!                  write (logString, *)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
                   do ii=incell(jcl)+1,incell(jcl+1)
 !     count the total number of neighbours
@@ -799,8 +834,8 @@
          enddo
 !     A few more cells need to be searched, so that we cover 13 of the 26
-!     neighbouring cells. 
+!     neighbouring cells.
 !        1
          jx=ix+1
-         if (jx.ge.ndx) jx=0 
+         if (jx.ge.ndx) jx=0
          jy=iy
          jz=iz-1
@@ -883,12 +918,7 @@
                r2=(xij*xij+yij*yij+zij*zij)
 
-               if (r2.le.exvcutg2) then 
-                  doit=.false.
-                  if (abs(iat2-iat1).gt.mxconr ) then
-                    doit=.true.
-                  else if (matcon(iat2-iat1,iat1).eq.1) then
-                    doit=.true.
-                  endif
-                  if (doit) then
+               if (r2.le.exvcutg2) then
+                  if (abs(iat2-iat1).gt.mxconr.or.
+     &                 matcon(iat2-iat1,iat1).eq.1) then
                      iatt1=ityat(iat1)
                      iatt2=ityat(iat2)
@@ -931,4 +961,8 @@
       include 'incl_lund.h'
 
+      double precision etmp, etmp1, xij, yij, zij, r2, r6
+
+      integer iat1, iat2, iatt1, iatt2
+
       etmp=0.0
       etmp1=0.0
@@ -941,5 +975,5 @@
             r2=(xij*xij+yij*yij+zij*zij)
 
-            if (r2.le.exvcutg2) then 
+            if (r2.le.exvcutg2) then
                if (abs(iat2-iat1).gt.mxconr.or.
      &                 matcon(iat2-iat1,iat1).eq.1) then
@@ -951,5 +985,5 @@
      &                 +bsaexv(iatt1,iatt2)*r2
                   etmp=etmp+etmp1
-               else 
+               else
 !                  print *,'atoms ', iat1,' and ',iat2,' were close',
 !     #                 'but matcon is ',matcon(iat2-iat1,iat1)
Index: enyreg.f
===================================================================
--- enyreg.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enyreg.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -21,4 +21,8 @@
 
 
+      double precision eny
+
+      integer i, nml, j
+
       eny = 0.d0
 
Index: enyshe.f
===================================================================
--- enyshe.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enyshe.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -1,8 +1,8 @@
 ! **************************************************************
 !
-! This file contains the subroutines: enyshe 
+! This file contains the subroutines: enyshe
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -20,5 +20,5 @@
 !
 ! The function loops over all moving sets within the molecule. Within
-! this loop it loops over the van-der-Waals domains of each atom in the 
+! this loop it loops over the van-der-Waals domains of each atom in the
 ! moving set and finally over the atoms that belong to the 1-4 interaction
 ! set.
@@ -28,4 +28,10 @@
 
 ! If nml == 0 calculate the interaction between all pairs.
+      double precision e0, vr, cqi, xi, yi, zi, xij, yij, zij, rij2
+      double precision rij4, rij6, rij, sr, ep
+
+      integer nml, ntlvr, ifivr, i1s, io, iv, ia, it, ic, i2s, ims, i1
+      integer i2, i, ity, ivw, j, jty, i14
+
       if (nml.eq.0) then
           ntlvr = nvr
@@ -33,5 +39,5 @@
         ntlvr=nvrml(nml)
       endif
-      
+
       if (ntlvr.eq.0) then
         write (*,'(a,i4)')
@@ -48,5 +54,5 @@
         ifivr = ivrml1(1)
         i1s = imsml1(ntlml) + nmsml(ntlml)
-      else 
+      else
 ! Index of first variable in molecule.
         ifivr=ivrml1(nml)
@@ -54,28 +60,28 @@
         i1s=imsml1(nml)+nmsml(nml)
       endif
-! Loop over moving sets/variables in reverse order     
-      do io=ifivr+ntlvr-1,ifivr,-1  
+! Loop over moving sets/variables in reverse order
+      do io=ifivr+ntlvr-1,ifivr,-1
 ! The array iorvr contains the variables in an "apropriate" order.
-        iv=iorvr(io)       
+        iv=iorvr(io)
 ! Index of the primary moving atom for the variable with index iv
-        ia=iatvr(iv)       
+        ia=iatvr(iv)
 ! Get the type of variable iv (valence length, valence angle, dihedral angle)
-        it=ityvr(iv)       
+        it=ityvr(iv)
 ! Class of variable iv's potential  (Q: What are they)
-        ic=iclvr(iv)       
+        ic=iclvr(iv)
 ! If iv is a dihedral angle ...
-        if (it.eq.3) then      
+        if (it.eq.3) then
 ! Barrier height * 1/2 of the potential of iv.
           e0=e0to(ic)
 ! Calculate the periodic potential term. sgto is the sign of the barrier, rnto is
 ! the periodicity and toat is torsion angle(?) associate with atom ia.
-          if (e0.ne.0.) 
+          if (e0.ne.0.)
      &         eyvr=eyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic)))
 ! else if iv is a valence angle ...
-        elseif (it.eq.2) then  
+        elseif (it.eq.2) then
 ! vr is the valence angle of ia
           vr=baat(ia)
 ! else if iv is a valence length...
-        elseif (it.eq.1) then  
+        elseif (it.eq.1) then
 ! vr is the length of the valence bond
           vr=blat(ia)
@@ -86,7 +92,7 @@
         i2s=i1s-1
 ! index of first moving set associated with iv
-        i1s=imsvr1(iv) 
+        i1s=imsvr1(iv)
 ! Loop over all moving sets starting from the one associated with vr to the end.
-        do ims=i1s,i2s  
+        do ims=i1s,i2s
 ! First atom of the current moving set
           i1=latms1(ims)
@@ -94,5 +100,5 @@
           i2=latms2(ims)
 ! Loop over all atoms of the current moving set.
-          do i=i1,i2  
+          do i=i1,i2
 ! Atom class of current atom
             ity=ityat(i)
@@ -104,9 +110,9 @@
             zi=zat(i)
 ! Loop over the atoms of the van der Waals domain belonging to atom i
-            do ivw=ivwat1(i),ivwat2(i)  
-! Loop over the atoms of the van der Waals domain of the atoms of the 
+            do ivw=ivwat1(i),ivwat2(i)
+! Loop over the atoms of the van der Waals domain of the atoms of the
 ! van der Waals domain of atom i
 ! Q: Which atoms are in these domains?
-              do j=lvwat1(ivw),lvwat2(ivw)  
+              do j=lvwat1(ivw),lvwat2(ivw)
 ! Atom type of partner
                 jty=ityat(j)
@@ -139,11 +145,11 @@
                 endif
 
-              enddo  
-            enddo  
-            
+              enddo
+            enddo
+
 ! Loop over 1-4 interaction partners
 ! The interactions between atoms that are three bonds apart in the protein are
 ! dominated by quantum mechanical effects. They are treated separately.
-            do i14=i14at1(i),i14at2(i)   
+            do i14=i14at1(i),i14at2(i)
               j=l14at(i14)
 
Index: enyshe_p.f
===================================================================
--- enyshe_p.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enyshe_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -1,8 +1,8 @@
 ! **************************************************************
 !
-! This file contains the subroutines: enyshe 
+! This file contains the subroutines: enyshe
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -14,11 +14,11 @@
 
 !     ............................................................................
-!     
+!
 !     PURPOSE: Calculate internal energy of molecule 'nml' with ECEPP parameters
-!     
+!
 !     CALLS: none
-!     
+!
 !     The function loops over all moving sets within the molecule. Within
-!     this loop it loops over the van-der-Waals domains of each atom in the 
+!     this loop it loops over the van-der-Waals domains of each atom in the
 !     moving set and finally over the atoms that belong to the 1-4 interaction
 !     set.
@@ -31,4 +31,12 @@
 
 ! If nml == 0 calculate the interaction between all pairs.
+      double precision eysmsum, teysm, teyel, teyvw, teyhb, teyvr
+      double precision startwtime, e0, vr, cqi, xi, yi, zi, xij, yij
+      double precision zij, rij2, rij4, rij6, rij, sr, ep, endwtime
+
+      integer nml, ntlvr, ifivr, i1s, iend, istart, loopcounter, io, iv
+      integer ia, it, ic, i2s, ims, i1, i2, i, ity, ivw, j, jty, i14
+      integer ierror
+
       if (nml.eq.0) then
          ntlvr = nvr
@@ -36,5 +44,5 @@
          ntlvr=nvrml(nml)
       endif
-      
+
       if (ntlvr.eq.0) then
          write (*,'(a,i4)')
@@ -56,5 +64,5 @@
          ifivr = ivrml1(1)
          i1s = imsml1(ntlml) + nmsml(ntlml)
-      else 
+      else
 !     Index of first variable in molecule.
          ifivr=ivrml1(nml)
@@ -62,5 +70,5 @@
          i1s=imsml1(nml)+nmsml(nml)
       endif
-!     Loop over variables in reverse order     
+!     Loop over variables in reverse order
 !     This is the first loop to parallize. We'll just split the moving sets
 !     over the number of available processors and sum the energy up in the end.
@@ -72,6 +80,6 @@
       startwtime = MPI_Wtime()
       loopcounter = 0
-!      do io=ifivr+ntlvr-1,ifivr,-1  
-      do io = workPerProcessor(nml, myrank) - 1, 
+!      do io=ifivr+ntlvr-1,ifivr,-1
+      do io = workPerProcessor(nml, myrank) - 1,
      &        workPerProcessor(nml, myrank+1), -1
          if (io.lt.istart) then
@@ -79,25 +87,25 @@
          endif
 !     The array iorvr contains the variables in an "apropriate" order.
-         iv=iorvr(io)       
+         iv=iorvr(io)
 !     Index of the primary moving atom for the variable with index iv
-         ia=iatvr(iv)       
+         ia=iatvr(iv)
 !     Get the type of variable iv (valence length, valence angle, dihedral angle)
-         it=ityvr(iv)       
+         it=ityvr(iv)
 !     Class of variable iv's potential  (Q: What are they)
-         ic=iclvr(iv)       
+         ic=iclvr(iv)
 !     If iv is a dihedral angle ...
-         if (it.eq.3) then      
+         if (it.eq.3) then
 !     Barrier height * 1/2 of the potential of iv.
             e0=e0to(ic)
 !     Calculate the periodic potential term. sgto is the sign of the barrier, rnto is
 !     the periodicity and toat is torsion angle(?) associate with atom ia.
-            if (e0.ne.0.) 
+            if (e0.ne.0.)
      &           teyvr=teyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic)))
 !     else if iv is a valence angle ...
-         elseif (it.eq.2) then  
+         elseif (it.eq.2) then
 !     vr is the valence angle of ia
             vr=baat(ia)
 !     else if iv is a valence length...
-         elseif (it.eq.1) then  
+         elseif (it.eq.1) then
 !     vr is the length of the valence bond
             vr=blat(ia)
@@ -108,7 +116,7 @@
          i2s=i1s-1
 !     index of first moving set associated with iv
-         i1s=imsvr1(iv) 
+         i1s=imsvr1(iv)
 !     Loop over all moving sets starting from the one associated with vr to the end.
-         do ims=i1s,i2s  
+         do ims=i1s,i2s
 !     First atom of the current moving set
             i1=latms1(ims)
@@ -116,5 +124,5 @@
             i2=latms2(ims)
 !     Loop over all atoms of the current moving set.
-            do i=i1,i2  
+            do i=i1,i2
 !     Atom class of current atom
                ity=ityat(i)
@@ -126,9 +134,9 @@
                zi=zat(i)
 !     Loop over the atoms of the van der Waals domain belonging to atom i
-               do ivw=ivwat1(i),ivwat2(i)  
-!     Loop over the atoms of the van der Waals domain of the atoms of the 
+               do ivw=ivwat1(i),ivwat2(i)
+!     Loop over the atoms of the van der Waals domain of the atoms of the
 !     van der Waals domain of atom i
 !     Q: Which atoms are in these domains?
-                  do j=lvwat1(ivw),lvwat2(ivw)  
+                  do j=lvwat1(ivw),lvwat2(ivw)
 
                      loopcounter = loopcounter + 1
@@ -163,11 +171,11 @@
                      endif
 
-                  enddo  
-               enddo  
-               
+                  enddo
+               enddo
+
 !     Loop over 1-4 interaction partners
 !     The interactions between atoms that are three bonds apart in the protein are
 !     dominated by quantum mechanical effects. They are treated separately.
-               do i14=i14at1(i),i14at2(i)   
+               do i14=i14at1(i),i14at2(i)
                   loopcounter = loopcounter + 1
                   j=l14at(i14)
@@ -210,13 +218,13 @@
 
 !     Collect energies from all nodes and sum them up
-      call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION,  
+      call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION,
      &     MPI_SUM, my_mpi_comm, ierror)
       call MPI_ALLREDUCE(teyel, eyel, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
      &     my_mpi_comm, ierror)
-      call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
-     &     my_mpi_comm, ierror)
-      call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
-     &     my_mpi_comm, ierror)
-      call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
+      call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
+     &     my_mpi_comm, ierror)
+      call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
+     &     my_mpi_comm, ierror)
+      call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
      &     my_mpi_comm, ierror)
 
Index: enysol.f
===================================================================
--- enysol.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enysol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -10,5 +10,5 @@
 ! **************************************************************
 
-      
+
       real*8 function enysol(nmol)
 
@@ -17,5 +17,5 @@
 !
 !     Double Cubic Lattice algorithm for calculating the
-!     solvation energy of proteins using 
+!     solvation energy of proteins using
 !     solvent accessible area method.
 !
@@ -24,8 +24,8 @@
 !
 ! -------------------------------------------------------------
-! TODO: Check the solvent energy for multiple molecules     
+! TODO: Check the solvent energy for multiple molecules
 !     arguments
       integer nmol
-      
+
 !     functions
       integer nursat
@@ -41,5 +41,5 @@
       double precision sizes, trad, zmin, xmax, xmin, ymax, zmax
       double precision sdr, sdd, volume
-      
+
       dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat)
       dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox)
@@ -47,5 +47,5 @@
 
 !      common/ressurf/surfres(mxrs)
-      
+
       eyslh = 0.0
       eyslp = 0.0
@@ -61,5 +61,5 @@
       do i=nrslow,nrshi
        surfres(i) = 0.0d0
-      end do 
+      end do
 
       numat= nup - nlow + 1
@@ -68,5 +68,5 @@
          inbox(i)=0
       end do
-     
+
       asa=0.0d0
       vdvol=0.0d0
@@ -135,5 +135,5 @@
        stop
       end if
-        
+
 ! Let us shift the borders to home all boxes
 
@@ -175,6 +175,6 @@
         inbox(i+1)=inbox(i+1)+inbox(i)
       end do
-         
-        
+
+
 !   Sorting the atoms by the their box numbers
 
@@ -184,6 +184,6 @@
          indsort(jj)=i
          inbox(j)=jj-1
-      end do   
-         
+      end do
+
 ! Getting started
 
@@ -205,5 +205,5 @@
              ney=min(iy+1,ndy-1)
              nez=min(iz+1,ndz-1)
-                     
+
 !  Atoms in the boxes around
 
@@ -221,6 +221,6 @@
                end do
               end do
-             end do     
-                             
+             end do
+
              do  ia=inbox(ibox)+1,inbox(ibox+1)
                jbi=indsort(ia)
@@ -274,5 +274,5 @@
                      end do
  99                  continue
-    
+
                      if(ik.gt.nnei)then
                        surfc(il)=.true.
@@ -316,5 +316,5 @@
                  surfres(jres) = surfres(jres) + area
                end if
-             end do 
+             end do
            end if
           end do
@@ -341,5 +341,5 @@
       character lin*80
       integer i
-!    Skipping comment lines, which begin with '!'  
+!    Skipping comment lines, which begin with '!'
       read(20,'(a)') lin
       do while(lin(1:1).eq.'!')
@@ -352,5 +352,5 @@
 !       write(*,'(a,i5)') 'the number of points---->',npnt
 
-!    Read the surface points   
+!    Read the surface points
 
       do i=1,npnt
@@ -358,7 +358,7 @@
 !          write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
       end do
- 
+
       return
- 
+
       end
 
Index: enysol_p.f
===================================================================
--- enysol_p.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ enysol_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -10,5 +10,5 @@
 ! **************************************************************
 
-      
+
       real*8 function enysol(nmol)
 
@@ -19,5 +19,5 @@
 !
 !     Double Cubic Lattice algorithm for calculating the
-!     solvation energy of proteins using 
+!     solvation energy of proteins using
 !     solvent accessible area method.
 !
@@ -26,5 +26,18 @@
 !
 ! -------------------------------------------------------------
-! TODO: Check the solvent energy for multiple molecules     
+! TODO: Check the solvent energy for multiple molecules
+      double precision startwtime, avr_x, avr_y, avr_z, xmin, ymin, zmin
+      double precision xmax, ymax, zmax, rmax, diamax, sizex, sizey
+      double precision sizez, shiftx, shifty, shiftz, boxpp, trad, dx
+      double precision dy, dz, dd, akrad, dr, xyz, radb, radb2, sdd, sdr
+      double precision area, volume, eyslsum, endwtime
+
+      integer nmol, nrslow, nrshi, nlow, nup, i, numat, inbox, j, ndx
+      integer ndy, ndz, nqxy, ncbox, mx, my, mz, nboxj, numbox, jj
+      integer indsort, iboxmin, iboxmax, ibox, iz, iy, ix, lbn, nsx, nsy
+      integer nsz, nex, ney, nez, jcnt, jz, jy, jx, jbox, ii, look, ia
+      integer jbi, nnei, ib, jtk, il, lst, ilk, ik, icount, jres, nursat
+      integer ierror, nhx, mhx, nbt, mbt
+
       dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat)
       dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox)
@@ -33,5 +46,5 @@
 
 !       common/ressurf/surfres(mxrs)
-      real*8 tsurfres(mxrs) 
+      real*8 tsurfres(mxrs)
       startwtime = MPI_Wtime()
       root = 0
@@ -50,5 +63,5 @@
       do i=nrslow,nrshi
        surfres(i) = 0.0d0
-      end do 
+      end do
 
       numat= nup - nlow + 1
@@ -57,5 +70,5 @@
          inbox(i)=0
       end do
-     
+
       asa=0.0d0
       vdvol=0.0d0
@@ -124,5 +137,5 @@
        stop
       end if
-        
+
 ! Let us shift the borders to home all boxes
 
@@ -164,6 +177,6 @@
         inbox(i+1)=inbox(i+1)+inbox(i)
       end do
-         
-        
+
+
 !   Sorting the atoms by the their box numbers
 
@@ -173,6 +186,6 @@
          indsort(jj)=i
          inbox(j)=jj-1
-      end do   
-         
+      end do
+
 !    Getting started
 !    We have to loop over ncbox boxes and have no processors available
@@ -202,5 +215,5 @@
              ney=min(iy+1,ndy-1)
              nez=min(iz+1,ndz-1)
-                     
+
 !  Atoms in the boxes around
 
@@ -218,6 +231,6 @@
                end do
               end do
-             end do     
-                             
+             end do
+
              do  ia=inbox(ibox)+1,inbox(ibox+1)
                jbi=indsort(ia)
@@ -271,5 +284,5 @@
                      end do
  99                  continue
-    
+
                      if(ik.gt.nnei)then
                        surfc(il)=.true.
@@ -313,17 +326,17 @@
                  surfres(jres) = surfres(jres) + area
                end if
-             end do 
+             end do
            end if
 !           end do
 !          end do
        end do
-      call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION, 
+      call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION,
      &      MPI_SUM,my_mpi_comm, ierror)
 !       write(*,*) 'enysol>', myrank, eysl, eyslsum
       tsurfres = surfres
-      call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION, 
+      call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION,
      &      MPI_SUM,my_mpi_comm, ierror)
       eysl = eyslsum
-      
+
       endwtime = MPI_Wtime()
       if (myrank.le.-1) then
@@ -349,7 +362,9 @@
       subroutine tessel
       include 'INCL.H'
+      integer i
+
       character lin*80
 
-!    Skipping comment lines, which begin with '!'  
+!    Skipping comment lines, which begin with '!'
 
       read(20,'(a)') lin
@@ -363,14 +378,14 @@
 !        write(*,'(a,i5)') 'the number of points---->',npnt
 
-!    Read the surface points   
+!    Read the surface points
 
       do i=1,npnt
          read(20,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
-         
+
 !        write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
       end do
- 
+
       return
- 
+
       end
 
Index: esolan.f
===================================================================
--- esolan.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ esolan.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -14,8 +14,8 @@
 ! -----------------------------------------------------------------
 !
-! Calculates the solvation energy of the protein with 
-! solavtion parameters model E=\sum\sigma_iA_i. 
+! Calculates the solvation energy of the protein with
+! solavtion parameters model E=\sum\sigma_iA_i.
 ! The solvent accessible surface area per atom
-! and its gradients with respect to the Cartesian coordinates of 
+! and its gradients with respect to the Cartesian coordinates of
 ! the atoms are calculated using the projection method described in
 !
@@ -25,8 +25,8 @@
 !
 ! INPUT: nmol - the order number of the protein chain.
-!      The atomic coordinates, radii and solvation parameters  
+!      The atomic coordinates, radii and solvation parameters
 !      are taken from the global arrays of SMMP
 !
-! OUTPUT: global array gradan(mxat,3) contains the gradients of 
+! OUTPUT: global array gradan(mxat,3) contains the gradients of
 !         solvation energy per atom
 !         local array as(mxat) contains accessible surface area per atom
@@ -42,9 +42,9 @@
 
       include 'INCL.H'
-      
-      
+
+
       integer nmol, ks0, ks2
 Cf2py intent(in) nmol
-      
+
       parameter (ks0=mxat*(mxat+1))
       parameter (ks2=mxat+mxat)
@@ -53,9 +53,9 @@
 
       integer neib, neibp, ivrx, iv, ial, i, ij, j, iii, idi,ivr, ifu,
-     &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll, 
+     &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll,
      &          nb, numat, nyx
       real*8 vertex, ax, pol, as, ayx, ayx1, probe, dd, ddat, dadx, gp,
      &          dalp,dbet, daalp, dabet, vrx, dv, dx, dy, dz, dt, di,
-     &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold, 
+     &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold,
      &          energy, sss, a1, a, ab1, a2, aa, ac2, ac1, ab2, am2,
      &          am, aia, ak, alp, am1, ap1, amibe, amabe, amaal, amial,
@@ -63,8 +63,8 @@
      &          ba, bcd, bc, bet, c1, c, c2, cc, sfcg, caba, cg, cfsg,
      &          cfcg, daba, ct, cs, d1, d2, dbe, dal, dcotbma, dcbet,
-     &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp, 
-     &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22, 
+     &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp,
+     &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22,
      &          r22p, r2, pom, prob, r, r42p, ratom, rr, s, sfsg, sf,
-     &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv, 
+     &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv,
      &          zz, zv
       dimension neib(0:mxat),neibp(0:mxat),ivrx(ks0)
@@ -78,10 +78,10 @@
       integer ta2(0:mxat),ta3(0:mxat),fullarc(0:mxat),al(0:ks2)
       real*8 neibor(mxat,4)
-      
+
       real*8, dimension(:, :, :), allocatable :: grad
-      
+
       allocate(grad(mxat, mxat, 3))
-      
-      
+
+
 !     open(3,file='solvation.dat')! Output file with solvation data
 
@@ -90,5 +90,5 @@
       if (nmol.eq.0) then
          numat=iatrs2(irsml2(ntlml))-iatrs1(irsml1(1))+1 ! Number of atoms
-      else 
+      else
          numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms
       endif
@@ -96,5 +96,5 @@
       do i=1,numat
         do j=1,3
-         gradan(i,j)=0.0d0 
+         gradan(i,j)=0.0d0
       end do
       end do
@@ -102,10 +102,10 @@
       do i=1,numat
          xold(i)=xat(i)! Redefine the coordinates just for safety
-         yold(i)=yat(i)! 
-         zold(i)=zat(i)! 
-      pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main' 
+         yold(i)=yat(i)!
+         zold(i)=zat(i)!
+      pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main'
       As(i)=pi4*pol(i)! Initially the whole surface of the atom is accessible.
       end do
-      
+
              do iii=1,numat
               do jjj=1,numat
@@ -120,8 +120,8 @@
 !***************************
 
-      do 1 i=1,numat ! Lop over atoms "i"           
+      do 1 i=1,numat ! Lop over atoms "i"
 ! ----------------------------------------------------
-         
-       R=rvdw(i) 
+
+       R=rvdw(i)
        jj=0
         do j=1,numat !Find the neighbours of "i"
@@ -132,5 +132,5 @@
          write(*,*)'ERROR in data: centres of two atoms coincide!'
          write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
-     &                   xold(j),yold(j),zold(j),rvdw(j) 
+     &                   xold(j),yold(j),zold(j),rvdw(j)
          stop 'Centres of atoms coincide!'
        endif
@@ -146,5 +146,5 @@
          end if
          end do!Neighbours
-       
+
        neib(0)=jj   ! The number of neighbors of "i"
 
@@ -197,5 +197,5 @@
 !
         uga=grnd() ! Random \gamma angle
-        ufi=grnd() ! Random \pi angle 
+        ufi=grnd() ! Random \pi angle
         uga=pi*uga/2
         ufi=pi*ufi/2
@@ -219,6 +219,6 @@
         end do
       end do
-!       
-      
+!
+
        if(jjj.ne.0) then ! Rotation is necessary
          sfsg=sf*sg
@@ -243,5 +243,5 @@
       do jj=1,neib(0)
       j=neib(jj)
-       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+      
+       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+
      &            (ddat(j,4)-R)**2-pol(j)           ! a
         neibor(j,2)=-pom*ddat(j,2)                  ! b
@@ -259,5 +259,5 @@
        do while(k.gt.1)               ! B
         k=k-1
-! Analyse mutual disposition of every pair of neighbours                                                      
+! Analyse mutual disposition of every pair of neighbours
         do 13 L=neib(0),k+1,-1        ! A
 
@@ -283,5 +283,5 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.      ! a01 02
      &    -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! The circle neib(k) encloses circle neib(L) and the later is discarded 
+! The circle neib(k) encloses circle neib(L) and the later is discarded
            neib(0)=neib(0)-1
            do j=L,neib(0)
@@ -326,5 +326,5 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.  ! a03 02
      &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! Don't exclude neib(k) 
+! Don't exclude neib(k)
            neib(0)=neib(0)-1
            do j=k,neib(0)
@@ -370,5 +370,5 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.   ! a07 02
      &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! discard the circle neib(L) 
+! discard the circle neib(L)
            neib(0)=neib(0)-1
            do j=L,neib(0)
@@ -385,5 +385,5 @@
            p2=(b1-b2)*pom
 !-- Assign t and s coordinates and the order number of circles
-           vertex(iv,1)=(t+p1)/am 
+           vertex(iv,1)=(t+p1)/am
            vertex(iv,2)=(s-p2)/am
            vertex(iv,3)=neib(k)
@@ -394,5 +394,5 @@
            vertex(iv,3)=neib(k)
            vertex(iv,4)=neib(L)
-!--  
+!--
           endif                                                 ! 10
 
@@ -491,5 +491,5 @@
             dbe=amabe-amibe
             if(dal.lt.pi) then
-             As(i)=R22*(pi-dal)  
+             As(i)=R22*(pi-dal)
             elseif(dbe.lt.pi) then
              As(i)=R22*(pi-dbe)
@@ -510,5 +510,5 @@
        enddo
        do j=1,neib(0)
-! "iv" is the number of intersection point. 
+! "iv" is the number of intersection point.
         do k=1,iv
        if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30
@@ -532,5 +532,5 @@
        do k=1,ifu
          a=neibor(fullarc(k),1)
-         if(a.lt.0.d0) then 
+         if(a.lt.0.d0) then
           aia=-1.d0
          else
@@ -579,5 +579,5 @@
          ak=a*a
          daba=dabs(a)
-         if(a.lt.0d0) then 
+         if(a.lt.0d0) then
           aia=-1d0
           else
@@ -679,5 +679,5 @@
           do j=1,nyx
 
-!  
+!
 !  Escape 'bad' intersections.
            if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40
@@ -686,5 +686,5 @@
            ap1=ct+rr*dcos(prob) ! The middle point of the arc.
            ap2=cs+rr*dsin(prob)
-! Verify if the middle point belongs to a covered part. 
+! Verify if the middle point belongs to a covered part.
 ! If yes then omit.
            do ij=1,ial
@@ -754,12 +754,12 @@
              dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2
              dabet(4)=(2d0*ak*a2)/am2
-               dv(2)=R42*b*vv1      
-               dv(3)=R42*c*vv1      
-               dv(4)=(d-R42*a)*vv1      
-               dv(1)=-dv(4)*R42      
-               dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2      
-             dx(2)=-(d+R42*a)*dv(2)*vv2      
-             dx(3)=-(d+R42*a)*dv(3)*vv2      
-             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2      
+               dv(2)=R42*b*vv1
+               dv(3)=R42*c*vv1
+               dv(4)=(d-R42*a)*vv1
+               dv(1)=-dv(4)*R42
+               dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2
+             dx(2)=-(d+R42*a)*dv(2)*vv2
+             dx(3)=-(d+R42*a)*dv(3)*vv2
+             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2
              dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
      &             R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
@@ -769,5 +769,5 @@
      &             R42*ak)*aia*dv(3)/a*vv2
                dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0*
-     &             R42*ak)*aia*dv(4)/a*vv2      
+     &             R42*ak)*aia*dv(4)/a*vv2
              dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2
              dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2
@@ -820,5 +820,5 @@
              dii(4,2)=dii(1,2)*R42
              dii(4,3)=2d0*(ddat(jk,4)+R)*R42
-             
+
              do idi=1,3
               grad(i,jk,idi)=grad(i,jk,idi)+
@@ -826,5 +826,5 @@
      &                   di(3)*dii(3,idi)+di(4)*dii(4,idi)
 
-             enddo             
+             enddo
              do idi=1,4
                dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp
@@ -838,5 +838,5 @@
              di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt))
             di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt))
-             enddo             
+             enddo
 
              dii(1,1)=2d0*ddat(kk,2)
@@ -866,5 +866,5 @@
 40           continue
           enddo
-         
+
          else ! analyse the case with lines (not necessary after rotation).
 
@@ -920,5 +920,5 @@
            enddo
            probe(nyx+1)=ayx1(nyx)+1d0
-           
+
 
            bc=b*b+c*c
@@ -960,5 +960,5 @@
 22           continue
             endif
-         enddo       
+         enddo
 11      continue
 
@@ -1017,9 +1017,9 @@
 111   continue
 
-!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',     
+!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',
 !    & '   gradz   Rad    Atom'
 !     write(3,*)
 
-      
+
       do i=1,numat
         enr=As(i)*sigma(i)
@@ -1038,5 +1038,5 @@
 
       esolan=sss
-      
+
       deallocate(grad)
 !     close(3)
@@ -1046,4 +1046,4 @@
 200   format(2i4,3f16.6)
 700   format(i5,7f8.3,5x,a4)
-         
+
       end
Index: eyabgn.f
===================================================================
--- eyabgn.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ eyabgn.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -5,10 +5,10 @@
 !                      Jan H. Meinke, Sandipan Mohanty
 !
-! Corrections to ECEPP energy terms due to R. A. Abagyan et al. 
-! 
+! Corrections to ECEPP energy terms due to R. A. Abagyan et al.
+!
 ! Two terms are calculated: eyrccr and eyentr, representing respectively
-! c a term to slightly shift the backbone dihedral angle preferences in 
+! c a term to slightly shift the backbone dihedral angle preferences in
 ! the ECEPP potential slightly away from the helix region, and another
-! term to estimate the side-chain entropy from a given configuration. 
+! term to estimate the side-chain entropy from a given configuration.
 !
 !
@@ -16,4 +16,8 @@
       real*8 function eyrccr(nml)
       include 'INCL.H'
+      double precision et, rsscl
+
+      integer nml, istres, indres, i, ipsi, in, ica, ic, mlvr, iphi
+
       dimension iN(mxrs),iCa(mxrs),iC(mxrs),mlvr(mxvr)
       dimension iphi(mxrs),ipsi(mxrs)
@@ -34,7 +38,7 @@
          call tolost(mynm)
          if ((mynm.eq.'val').or.(mynm.eq.'ile').or.
-     &           (mynm.eq.'thr')) then 
+     &           (mynm.eq.'thr')) then
             rsscl=1.0
-         else 
+         else
             rsscl=0.5
          endif
@@ -48,5 +52,5 @@
       eyrccr=et
       return
-      end 
+      end
 
       subroutine init_abgn
@@ -54,4 +58,8 @@
 !       dimension rsstrg(mxrs)
 !       common /abgncor/rsstrg
+      double precision xarea, estrg
+
+      integer i, istres, indres, imytyp, j
+
       dimension xarea(nrsty),estrg(nrsty)
       character mynm*4
@@ -102,17 +110,17 @@
 !            print *,'comparing ',mynm,' with ',rsnmcd(j),imytyp
          enddo
-         if (imytyp.eq.0) then 
+         if (imytyp.eq.0) then
             print *, 'Unknown residue type ',seq(i)
             print *, 'Abagyan term strength set to 0'
             rsstrg(i)=0.0
-         else 
+         else
             rsstrg(i)=estrg(imytyp)/xarea(imytyp)
          endif
 !         print *,'residue ',i,seq(i),' type ',imytyp
-!         print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i) 
+!         print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i)
       enddo
       print *, 'initialized Abagyan corrections to ECEPP force field'
       end
-      
+
       real*8 function eyentr(nml)
       include 'INCL.H'
@@ -121,4 +129,8 @@
 !       common/ressurf/surfres(mxrs)
 !      common/bet/beta
+      double precision eentr, aars, strh
+
+      integer nml, istres, indres, i
+
       eentr=0
       if (nml.eq.0) then
@@ -136,5 +148,5 @@
 !        The maximal burial entropies were estimated at temperature 300k
 !        The values in the array estrg are k_B * T (=300k) * Entropy
-!        Presently we need it at temperature 1/beta, so we need to 
+!        Presently we need it at temperature 1/beta, so we need to
 !        multiply the strengths in estrg with (1/beta)/(300 kelvin)
 !        300 kelvin is approximately 0.59576607 kcal/mol.
@@ -149,7 +161,11 @@
       return
       end
-      
-      real*8 function eyabgn(nml) 
+
+      real*8 function eyabgn(nml)
       include 'INCL.H'
+      double precision eyrccr, eyentr
+
+      integer nml
+
       eyabgn=eyrccr(nml)+eyentr(nml)
 !      print *,'Abagyan term = ',eyabgn
Index: gradient.f
===================================================================
--- gradient.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ gradient.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -21,4 +21,8 @@
       include 'INCL.H'
 
+
+      double precision esm
+
+      integer i, ivr1, ivr2, j
 
       esm = 0.d0
Index: init_energy.f
===================================================================
--- init_energy.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ init_energy.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -5,5 +5,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -17,5 +17,5 @@
 ! PURPOSE: initialize energy parameters
 !        0  => ECEPP2 or ECEPP3 depending on the value of sh2
-!        1  => FLEX 
+!        1  => FLEX
 !        2  => Lund force field
 !        3  => ECEPP with Abagyan corrections
@@ -29,12 +29,14 @@
       include 'INCL.H'
 
+      integer ll, iendst, its
+
       character libdir*(*),tesfil*80
-      
+
       if (ientyp.eq.1) then
           flex = .true.
-      else 
+      else
           flex = .false.
       end if
-      
+
       lunlib=10
       ll=iendst(libdir)
@@ -56,5 +58,5 @@
       endif
 
-      if (ientyp .eq. 2) then 
+      if (ientyp .eq. 2) then
         reslib=libdir(1:ll)//'lib.lun'
       endif
@@ -78,5 +80,5 @@
 
       else
- 
+
         if (itysol.ne.0) then
           write(*,'(a)') ' init_energy>  undefined solvent type !'
@@ -110,4 +112,8 @@
 
       include 'INCL.H'
+
+      double precision hbc, hba, ri, ai, aei, aic, a, aj, c, rij
+
+      integer i, j, iac, ido, jac, jdo
 
       dimension hbc(mxhbdo,mxhbac),hba(mxhbdo,mxhbac)
@@ -308,4 +314,6 @@
 
       include 'INCL.H'
+      
+      integer i
 
 !  Atom types ------------------------------------------------------------
@@ -517,5 +525,5 @@
      & 'ile ','leu ','lys ','lys+','met ','phe ','cpro','pro ','cpru',
      & 'prou','pron','pro+','ser ','thr ','trp ','tyr ','val ' /
-  
+
       data (onltcd(i),i=1,nrsty)/  ! One-letter codes for amino acid types
      & 'A',   'R',   'R',   'N',   'D',   'D',   'C',   'C',   'Q',
@@ -527,11 +535,11 @@
 !      coefficients for their solvation free energy (kcal/molxA**2)
 
-!  Method: 
-
-!  itysol=1 : OONS --> T.Ooi, et al, 
+!  Method:
+
+!  itysol=1 : OONS --> T.Ooi, et al,
 !                      Proc. Natl. Acad. Sci. USA 8 (1987) 3086-3090.
-!  itysol=2 : JRF  --> J.Vila, et al, 
+!  itysol=2 : JRF  --> J.Vila, et al,
 !                      PROTEINS: Struct Funct Genet 10(1991) 199-218.
-!  itysol=3 : WE92 --> L.Wesson, D.Eisenberg, 
+!  itysol=3 : WE92 --> L.Wesson, D.Eisenberg,
 !                      Protein Science 1 (1992) 227-235.
 !  itysol=4 : SCH1 --> D.Eisenberg, et al,
@@ -539,9 +547,9 @@
 !  itysol=5 : SCH2 --> A.H.Juffer, et al,
 !                      Proteine Science 4 (1995) 2499-2509.
-!  itysol=6 : SCH3 --> L.Wesson, D.Eisenberg, 
+!  itysol=6 : SCH3 --> L.Wesson, D.Eisenberg,
 !                      Protein Science 1 (1992) 227-235.
 !  itysol=7 : SCH4 --> C.A. Schiffer, et al,
 !                      Mol. Simul. 10(1993) 121-149.
-!  itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan, 
+!  itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan,
 !                      Nature 319 (1986) 199-203.
 !  itysol=9 : BM   --> B. Freyberg, et al,
Index: init_molecule.f
===================================================================
--- init_molecule.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ init_molecule.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -16,5 +16,5 @@
 ! PURPOSE: construct starting structure of molecule(s)
 !
-!          iabin = 1  : ab Initio using sequence & 
+!          iabin = 1  : ab Initio using sequence &
 !                       variables given in input files
 !          iabin != 1 : sequence, variable information
@@ -32,9 +32,12 @@
       include 'INCP.H'
 
+      integer iabin, iendst, ntl, i, j, l, it, ier, ir, nursvr, i1, i2
+      integer its
+
 Cf2py character*80 optional, intent(in) :: seqfile = ' '
-Cf2py character*80 optional, intent(in) :: varfile = ' ' 
-      
+Cf2py character*80 optional, intent(in) :: varfile = ' '
+
       character grpn*4,grpc*4
-      character navr*3, nars*4  
+      character navr*3, nars*4
       character seqfile*80, varfile*80
       integer ontlml
@@ -44,7 +47,7 @@
       readFromStdin = .false.
 
-      write (*,*) 'init_molecule: Solvent: ', itysol
-      if (iabin.eq.1) then  
-         
+      write (logString, *) 'init_molecule: Solvent: ', itysol
+      if (iabin.eq.1) then
+
 !     ----------------------------------------- get sequence for molecule(s)
          lunseq=11
@@ -53,40 +56,41 @@
          endif
          if (iendst(seqfile).le.1.or.seqfile.eq.' ') then
- 1          write (*,'(/,a,$)') ' file with SEQUENCE:'
+ 1          write (logString, '(/,a,$)') ' file with SEQUENCE:'
             seqfil=' '
             read (*,'(a)',err=1) seqfil
             readFromStdin = .true.
-         else 
+         else
             seqfil = seqfile
-         endif 
+         endif
          call redseq
-         
-         write (*,*) 'File with sequence is ', seqfil(1:iendst(seqfil))
-         
+
+         write (logString, *) 'File with sequence is ', 
+     &      seqfil(1:iendst(seqfil))
+
 !     --------------------------------- read & assemble data from libraries
 !     initial coordinates, interaction lists
-         
+
          ntl = ntlml
          do i=ontlml, ntl
-            
+
             call getmol(i)      ! assemble data from libraries
-            
+
             do j=1,6            ! initialize global parameters
                gbpr(j,i)=0.d0
             enddo
-            
+
             call bldmol(i)      ! co-ordinates
-            
+
             ntlml = i
             call addend(i,grpn,grpc) ! modify ends
             call setmvs(i) ! determine sets of moving atoms for given variables
             call mklist(i)      ! compile lists of interaction partners
-            
+
          enddo
-         
+
 !     --------------------------- Read the initial conformation if necessary
-         if(readFromStdin) then 
-            write (*,'(a,$)') ' file with VARIABLES:'
-!     
+         if(readFromStdin) then
+            write (logString, '(a,$)') ' file with VARIABLES:'
+!
             varfil=' '
             read(*,'(a)',end=2,err=2) varfil
@@ -96,46 +100,12 @@
          l=iendst(varfil)
          if (l.gt.0.and.varfil.ne.' ') then
-            write (*,'(1x,a,/)') varfil(1:l)
+            write (logString, '(1x,a,/)') varfil(1:l)
             lunvar=13
-            
+
             call redvar         ! get vars. and rebuild
-            
-         endif
-         
- 2       write(*,*) ' '
-         
-         ireg = 0
-         
-      else                      ! =========================== from PDB
-         if (iendst(seqfile).le.1) then
- 3          write (*,'(/,a,$)') ' PDB-file:'
-            seqfil=' '
-            read (*,'(a)',err=3) seqfil
-         else
-            seqfil = seqfile
-         endif
-         write (*,*) 'PDB structure ',seqfil(1:iendst(seqfil))
-         print *, 'calling readpdb with ',seqfile
-         call pdbread(seqfil,ier)
-         
-         if (ier.ne.0) stop
-         
-         call pdbvars()
-         
-         ireg = 1
-         
-      endif
-      
-! If Lund force field is in use, keep omega angles fixed
-      if (ientyp.eq.2) then
-         do iv=1,nvrml(ntlml)
-            if ((nmvr(iv)(1:2).eq.'om')) then
-                vlvr(iv)=pi
-                toat(iatvr(iv))=pi
-                fxvr(iv)=.true.
-                print *, 'Fixed variable ',iv,nmvr(iv),vlvr(iv)
-            endif
-         enddo
-      endif
+
+         endif
+
+ 2       write (logString, *) ' '
 
 !     -------------------- get: nvr,idvr, vlvr, olvlvr
@@ -161,69 +131,91 @@
          enddo
 
+         ireg = 0
+
+      else                      ! =========================== from PDB
+         if (iendst(seqfile).le.1) then
+ 3          write (logString, '(/,a,$)') ' PDB-file:'
+            seqfil=' '
+            read (*,'(a)',err=3) seqfil
+         else
+            seqfil = seqfile
+         endif
+         write (logString, *) 'PDB structure ',seqfil(1:iendst(seqfil))
+         print *, 'calling readpdb with ',seqfile
+         call pdbread(seqfil,ier)
+
+         if (ier.ne.0) stop
+
+         call pdbvars()
+
+         ireg = 1
+
+      endif
+
 !     -------------------------- set var. amplitudes for simulations
-      
+
       do i=1,ivrml1(ntlml)+nvrml(ntlml)-1
-         
+
          if (ityvr(i).eq.3.and..not.fxvr(i)) then ! torsion
-            
+
             navr = nmvr(i)
-            
+
             ir = nursvr(i)
             nars = seq(ir)
-            
+
             if (                         navr(1:2).eq.'om'
-            
+
      &     .or.nars(1:3).eq.'arg'.and.(navr(1:2).eq.'x5'
      &           .or.navr(1:2).eq.'x6')
-            
+
      &           .or.(nars(1:3).eq.'asn'.or.nars(1:3).eq.'asp')
      &           .and.navr(1:2).eq.'x3'
-            
+
      &           .or.(nars(1:3).eq.'gln'.or.nars(1:3).eq.'glu')
      &           .and.navr(1:2).eq.'x4'
-            
+
      &           ) then
-            
+
 !     axvr(i) = pi/9.d0  ! 20 deg.
             axvr(i) = pi2       ! Trying out 360 deg. for these as well
-            
+
          else
             axvr(i) = pi2       ! 360 deg.
          endif
-         
+
       else
          axvr(i) = 0.d0
       endif
-      
+
       enddo                     ! vars.
-      
+
 !     --------------------- initialize solvation pars. if necessary
 
       if (itysol.ne.0) then
-         
+
          i1=iatrs1(irsml1(1))   ! 1st atom of 1st molecule
          i2=iatrs2(irsml2(ntlml)) ! last atom of last molecule
-         
+
          its = iabs(itysol)
-         
+
          do i=i1,i2             ! all atoms
             it=ityat(i)
             sigma(i)=coef_sl(its,it)
             rvdw(i) =rad_vdw(its,it)
-            
+
             if (nmat(i)(1:1).ne.'h') rvdw(i)=rvdw(i)+rwater
-            
+
          enddo
-         
+
       endif
 ! Initialize calpha array
       do i=ontlml, ntlml
-         call c_alfa(i,1)    
+         call c_alfa(i,1)
       enddo
 
 !     Initialize arrays used in the BGS update
-      call init_lund()      
+      call init_lund()
       return
       end
-      
-      
+
+
Index: main.f
===================================================================
--- main.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ main.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -16,9 +16,14 @@
 ! $Id: main.f 334 2007-08-07 09:23:59Z meinke $
 ! **************************************************************
-      
+
       program main
 
       include 'INCL.H'
       include 'INCP.H'
+      double precision eps, temp, e, energy
+
+      integer maxloglevel, logfileunit, iabin, imin, maxit, nequi
+      integer nsweep, nmes, ncalls, nacalls
+
       common/updstats/ncalls(5),nacalls(5)
       character*80 libdir, seqfile, varfile
@@ -32,7 +37,7 @@
 !            Directory for SMMP libraries
 !     Change the following directory path to where you want to put SMMP
-!     libraries of residues. 
+!     libraries of residues.
       libdir='./SMMP/'
-      
+
 !     Set the maximum log level. The larger the number the more detailed
 !     the log.
@@ -48,5 +53,5 @@
       ientyp = 0
 !        0  => ECEPP2 or ECEPP3 depending on the value of sh2
-!        1  => FLEX 
+!        1  => FLEX
 !        2  => Lund force field
 !        3  => ECEPP with Abagyan corrections
@@ -70,13 +75,14 @@
       iabin = 1  ! =0: read from PDB-file
                  ! =1: ab Initio from sequence (& variables)
-      seqfile='EXAMPLES/enkefa.seq'
-      varfile='EXAMPLES/enkefa.ann'
-!       varfile = ' '
-      
+      seqfile='polyq.seq'
+!      seqfile='polyA.pdb'
+      varfile='polyq.var'
+       varfile = ' '
+
       ntlml = 0
       write (*,*) 'Solvent: ', itysol
 !     Initialize random number generator.
       call sgrnd(31433)
-      
+
       if (itysol.eq.0.and.ientyp.eq.3) then
          print *,'Can not use Abagyan entropic corrections without '
@@ -87,8 +93,8 @@
       call init_molecule(iabin,grpn,grpc,seqfile,varfile)
 
-! Decide if and when to use BGS, and initialize Lund data structures 
+! Decide if and when to use BGS, and initialize Lund data structures
       bgsprob=0.75   ! Prob for BGS, given that it is possible
-! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 
-! 2 => temperature dependent choice 
+! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
+! 2 => temperature dependent choice
       upchswitch=1
       rndord=.true.
@@ -96,5 +102,5 @@
       if (ientyp.eq.2) call init_lundff
       if (ientyp.eq.3) call init_abgn
-      
+
 
 ! ========================================  Add your task down here
@@ -103,14 +109,17 @@
       maxit = 15000 ! maximum number of iterations in minimization
       eps = 1.0d-7 ! requested precision
-      call minim(imin, maxit, eps)
-      call outvar(0, ' ')
+!      call minim(imin, maxit, eps)
+!      call outvar(0, ' ')
 !     To do a canonical Monte Carlo simulation uncomment the lines below
-!       nequi = 100
-!       nsweep = 50000
-!       nmes = 10
-!       temp = 300.0
-!       lrand = .true.
+       nequi = 100
+       nsweep = 50000
+       nmes = 10
+       temp = 300.0
+       lrand = .true.
+       E = energy()
+       write(*,*) E, eyel,eyvw,eyhb,eyvr
+       call outpdb(1, "polyA.pdb")
 !      Canonical Monte Carlo
-!       call canon(nequi, nsweep, nmes, temp, lrand)
+!        call canon(nequi, nsweep, nmes, temp, lrand)
 
 !      For simulated annealing uncomment the lines below
@@ -118,4 +127,4 @@
 !      tmax = 500.0
 !      call anneal(nequi, nsweep, nmes, tmax, tmin, lrand);
-! ========================================  End of main      
+! ========================================  End of main
        end
Index: main_p.f
===================================================================
--- main_p.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ main_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -1,16 +1,16 @@
 !     **************************************************************
-!     
+!
 !     This file contains the   main (PARALLEL TEMPERING  JOBS ONLY,
 !     FOR SINGULAR PROCESSOR JOBS USE main)
-!     
+!
 !     This file contains also the subroutine: p_init_molecule
-!     
+!
 !     Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!     Shura Hayryan, Chin-Ku 
+!     Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
-!     
+!
 !     CALLS init_energy,p_init_molecule,partem_p
-!     
+!
 !     **************************************************************
       program pmain
@@ -21,4 +21,9 @@
       include 'mpif.h'
 
+      double precision startwtime, group_world, error, endwtime
+
+      integer ierr, num_proc, iabin, nequi, nswp, nmes, nsave, ifrm, j
+      integer i, nml, nresi, my_pt_rank, ncalls, nacalls
+
       character*80 libdir
       character*80 in_fil,ou_fil,filebase, varfile
@@ -28,5 +33,5 @@
       logical newsta
 
-!c    Number of replicas 
+!c    Number of replicas
       integer num_replica
 !c    Number of processors per replica
@@ -62,5 +67,5 @@
 
 !     =================================================== Energy setup
-      libdir='SMMP/'      
+      libdir='SMMP/'
 !     Directory for SMMP libraries
 
@@ -71,8 +76,8 @@
       ientyp = 0
 !     0  => ECEPP2 or ECEPP3 depending on the value of sh2
-!     1  => FLEX 
+!     1  => FLEX
 !     2  => Lund force field
 !     3  => ECEPP with Abagyan corrections
-!     
+!
 
       sh2=.false.               ! .true. for ECEPP/2; .false. for ECEPP3
@@ -106,28 +111,28 @@
                             !   temperatures must have at least as many
                             !   entries
-      nequi=10              ! Number of MC sweeps before measurements 
-                            !   and replica exchanges are started 
+      nequi=10              ! Number of MC sweeps before measurements
+                            !   and replica exchanges are started
       nswp=500000           ! Number of sweeps
       nmes=10               ! Interval for measurements and replica exchange
       nsave=1000            ! Not used at the moment
-                    
-      switch = -1           ! How should the configuration be   
+
+      switch = -1           ! How should the configuration be
                             !   initialized?
-                            ! -1 stretched chain 
+                            ! -1 stretched chain
                             !  0 don't do anything
                             !  1 initialize each angle to a random value
-      
+
       ifrm=0
       ntlml = 0
 
-! Decide if and when to use BGS, and initialize Lund data structures 
+! Decide if and when to use BGS, and initialize Lund data structures
       bgsprob=0.6    ! Prob for BGS, given that it is possible
-! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 
-! 2 => temperature dependent choice 
+! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
+! 2 => temperature dependent choice
       upchswitch=1
       rndord=.true.
 !     =================================================================
 !     Distribute nodes to parallel tempering tasks
-!     I assume that the number of nodes available is an integer 
+!     I assume that the number of nodes available is an integer
 !     multiple n of the number of replicas. Each replica then gets n
 !     processors to do its energy calculation.
@@ -139,6 +144,6 @@
 !     could just use i * num_ppr but this way it's more flexible.
       j = 0
-      do i = 1, num_replica 
-         ranks(i) = j 
+      do i = 1, num_replica
+         ranks(i) = j
          proc_range(1) = j
          proc_range(2) = j + num_ppr - 1
@@ -146,5 +151,5 @@
          call mpi_group_range_incl(group_world, 1, proc_range, group(i)
      &                              ,error)
-         write (*,*) "Assigning rank ", j, proc_range, 
+         write (*,*) "Assigning rank ", j, proc_range,
      &               "to group", group(i)
          call flush(6)
@@ -167,5 +172,5 @@
       call mpi_group_incl(group_world, num_replica, ranks, group_partem,
      &                    error)
-      call mpi_comm_create(mpi_comm_world, group_partem, partem_comm, 
+      call mpi_comm_create(mpi_comm_world, group_partem, partem_comm,
      &                     error)
 
@@ -184,5 +189,5 @@
          varfile = 'EXAMPLES/1bdd.var'
          call init_molecule(iabin, grpn, grpc,in_fil,varfile)
-      else 
+      else
          filebase = "conf_0000.var"
          call init_molecule(iabin, grpn, grpc,in_fil,
@@ -220,5 +225,5 @@
       nml = 1
       call distributeWorkLoad(no, nml)
-      
+
       call partem_p(num_replica, nequi, nswp, nmes, nsave, newsta,
      &              switch, rep_id, partem_comm)
Index: metropolis.f
===================================================================
--- metropolis.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ metropolis.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -20,4 +20,10 @@
       include 'INCP.H'
       include 'incl_lund.h'
+      
+      double precision eol, energy, grnd, dv, addang, enw, delta, dummy
+      double precision rd, ex, acz, eol1
+
+      integer nsw, iupstate, iupt, ivar, jv, iml, i, ncalls, nacalls
+
       external dummy
 !      common/bet/beta
@@ -178,4 +184,5 @@
       
       subroutine accanalyze(iuptype,iupdstate)
+      integer ncalls, nacalls, iupdstate, iuptype
       common/updstats/ncalls(5),nacalls(5)
       ncalls(5)=ncalls(5)+1
@@ -202,4 +209,5 @@
       integer iiii
       double precision bbbb,curprob, grnd,up2bmax,up2bmin
+
       common/updtparam/up2bmax,up2bmin
       curprob=(bbbb-up2bmin)/(up2bmax-up2bmin)
@@ -211,4 +219,5 @@
       end
       block data updtchs
+      integer ncalls, nacalls
       double precision up2bmax, up2bmin
       common/updtparam/up2bmax,up2bmin
Index: mincjg.f
===================================================================
--- mincjg.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ mincjg.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -40,6 +40,11 @@
 ! ....................................................................
 
-      implicit real*8 (a-h,o-z)
-      implicit integer*4 (i-n)
+      double precision amf, tol, eps, gsqrd, fmin, f, t, g, ga, d, xa, x
+      double precision gsq2, gnew, dfpr, stmin, fuit, gdit, gt, gdmi
+      double precision sbound, stepch, step, wo, sum, fch, acur, ddspln
+      double precision gspln, beta, gama, yt, dt, gamden, di
+
+      integer mxfcon, maxlin, mxn, ier, iter, iterfm, iterrs, nfun
+      integer nfopt, i, n, nfbeg, iretry, maxfun
 
       parameter (AMF = 10.d0,
@@ -102,5 +107,5 @@
 
     2 step = stmin + stepch
-      wo = 0.d0 
+      wo = 0.d0
 
       do i=1,n
@@ -182,5 +187,5 @@
 
       if ( fch .ne. 0.d0 )  ddspln = ddspln + (wo + wo) / stepch
- 
+
       if ( gdmi .eq. 0.d0 )  goto 6
 
Index: nursvr.f
===================================================================
--- nursvr.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ nursvr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -18,4 +18,6 @@
 ! ...........................................................
       include 'INCL.H'
+
+      integer i, ifirs, ivr, j
 
       do i=ntlml,1,-1
@@ -45,4 +47,6 @@
       include 'INCL.H'
 
+      integer i, ifirs, ilars, iat, j
+
       do i=1,ntlml
 
Index: opereg.f
===================================================================
--- opereg.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ opereg.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -33,4 +33,13 @@
       include 'INCL.H'
       include 'INCP.H'
+
+      double precision xfvr, yfvr, zfvr, xfrvr, yfrvr, zfrvr, x1, y1, z1
+      double precision a, sa, ca, xk, yk, zk, xi, yi, zi, xji, yji, zji
+      double precision dx, dy, dz, x, y, z, xfat, yfat, zfat, xfrat
+      double precision yfrat, zfrat, xb, yb, zb, ex, ey, ez, xfiv, yfiv
+      double precision zfiv, xfriv, yfriv, zfriv
+
+      integer ntlvr, nml, ix2, ifivr, ilavr, i, ii, j, i1s, i1a, io, iv
+      integer it, ia, ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, i14
 
       dimension xfat(mxat),yfat(mxat),zfat(mxat),
@@ -100,5 +109,5 @@
 
           eyrg = eyrg + xji**2 + yji**2 + zji**2  ! The regularization energy is just
-                                                  ! the sum over the atom distances 
+                                                  ! the sum over the atom distances
                                                   ! squared.
 
@@ -303,4 +312,8 @@
       include 'INCL.H'
 
+      double precision del, pro, gdn, enyreg
+
+      integer ii, nml, i
+
       parameter (del=1.d-7)
 
@@ -339,4 +352,8 @@
 
       include 'INCL.H'
+
+      double precision del, vlvrx, ovr, eynw, enyreg, gdn, gda
+
+      integer i, it, iv, nml
 
       parameter (del=1.d-6)
Index: opesol.f
===================================================================
--- opesol.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ opesol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -39,5 +39,5 @@
       integer ntlvr, nml, ix2, ifivr, ilavr, i, i1s, i1a, io, iv, it, ia
       integer ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, j, i14
-      
+
       dimension xfat(mxat),yfat(mxat),zfat(mxat),
      &          xfrat(mxat),yfrat(mxat),zfrat(mxat),
@@ -251,5 +251,5 @@
 
       include 'INCL.H'
-      
+
       double precision del, vlvrx, ovr, eynw, esolan, gda, gdn
 
Index: outpdb.f
===================================================================
--- outpdb.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ outpdb.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -24,4 +24,9 @@
 
       include 'INCL.H'
+
+      double precision occ, bva
+
+      integer i0, i9, nml, im1, im2, ibegst, iout, iat, iml, irs, ifirs
+      integer ifiat, nrs, i, j, iendst, nfi, ibd, jj, nbd
 
       dimension ibd(4)
@@ -99,5 +104,5 @@
             call toupst(atnm)
             call toupst(res)
- 
+
             linout = ' '
             write (linout,1) linty,iat,atnm,res(1:3),chid,irs,cdin,
@@ -170,5 +175,5 @@
       write (linout,'(a3)') 'END'
       write(iout,'(a80)') linout
-      
+
       close(iout)
 
Index: outvar.f
===================================================================
--- outvar.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ outvar.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,15 +4,15 @@
 !
 ! Copyright 2005       Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
 !
 ! **************************************************************
- 
+
       subroutine outvar(nml,fileName)
-  
+
 !--------------------------------------------------------------------
-!       Output of variables of the current protein conformation 
-!- 
+!       Output of variables of the current protein conformation
+!-
 !       nml != 0       :     molecule index
 !       nml  = 0       :     all molecules
@@ -24,4 +24,6 @@
 
       include 'INCL.H'
+
+      integer nml, ibegst, iout
 
       character*(*) fileName
@@ -49,4 +51,6 @@
 
       include 'INCL.H'
+
+      integer nml, im1, im2, iml, i, iout, ibegst, ifivr, is, nursvr, it
 
       character mlfd*7,strg(6)*17
Index: partem_p.f
===================================================================
--- partem_p.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ partem_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -43,8 +43,8 @@
 
       integer ifrrm, nmes, nswp, num_rep, i, j, nresi, iold, inode
-      integer intem, iv, jold, idum1, idum2, idum3, mpi_integer
-      integer mpi_comm_world, ierr, mpi_double_precision, nsw, nequi
+      integer intem, iv, jold, idum1, idum2, idum3
+      integer ierr, nsw, nequi
       integer nml, nhel, mhel, nbet, mbet, mhb, imhb, nctot, ncnat
-      integer mpi_comm_null, k1, k, nu, no1, in, jn
+      integer k1, k, nu, no1, in, jn
       
       dimension  eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC),
@@ -225,5 +225,5 @@
             call flush(6)
          endif
-         if(mod(iold,nmes).eq.0) then
+         if(mod(iold,nmes).eq.0.and.myrank.eq.0) then
             if ((rep_id + 1).eq.trackID.and.myrank.eq.0) then
                frame = iold /nmes
Index: partem_s.f
===================================================================
--- partem_s.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ partem_s.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -10,5 +10,5 @@
 !
 ! **************************************************************
-      
+
       subroutine partem_s(num_rep, nequi, nswp, nmes, nsave, newsta,
      &                     switch)
@@ -20,4 +20,11 @@
 !
       include 'INCL.H'
+      
+      double precision gasc, dv, grnd, vr, addang, energ, energy, acz
+      double precision ee, rgy, temp0, delta, ra
+
+      integer i, j, nml, k, nstart, iv, nsw, iswitch, k1, num_rep1, nu
+      integer in, jn
+
       external can_weight
 ! TODO Store global coordinates in pgbpr
@@ -44,5 +51,5 @@
 !     ipoi:   Points to replica ipoi(k) which is currently
 !             at inverse temperature pbe(k)
-!     eol:    energy of each replica 
+!     eol:    energy of each replica
 !     acc:    accepatance rate of each replica
 !
@@ -55,5 +62,5 @@
       allocate(temp(num_rep), pbe(num_rep), eol(num_rep))
       allocate(acc(num_rep), ipoi(num_rep))
-      
+
       if (.not.(allocated(coor_G).and.allocated(temp)
      &          .and. allocated(pbe).and. allocated(eol)
@@ -77,5 +84,5 @@
          ipoi(i) = 0
       end do
-        
+
 
 !     READ IN TEMPERATURES
@@ -90,5 +97,5 @@
 
       open(13,file='time.d',status='unknown')
-      
+
       if(.not.newsta) then
 ! READ START Values
@@ -122,7 +129,7 @@
                   pgbpr(j, nml, i) = gbpr(j, nml)
                end do
-            end do            
-         end do
-      
+            end do
+         end do
+
 !      Equilibrization for each replica (No replica exchange move)
          do k=1,num_rep
@@ -141,5 +148,5 @@
                CALL METROPOLIS(energ,acz,can_weight)
             end do
-            write(*,*) 'Start energy after equilibration for replica:', 
+            write(*,*) 'Start energy after equilibration for replica:',
      &                 k, energ
             do i=1,nvr
@@ -155,5 +162,5 @@
          end do
       end if
-            
+
 ! Now begins the simulation with Multiple Markov Chains
       iswitch = 1
@@ -176,5 +183,5 @@
             end do
             CALL METROPOLIS(energ,acz,can_weight)
-! 
+!
             if(mod(nsw,nmes).eq.0) then
 ! Measure and store here all quantities you want to analyse later
@@ -195,5 +202,5 @@
                end do
             end do
-        
+
             eol(k) = energ
             acc(k) = acz
@@ -214,5 +221,5 @@
             j=i+1
             if(i.eq.num_rep) j=1
-            in=ipoi(i)      
+            in=ipoi(i)
             jn=ipoi(j)
             delta=-pbe(i)*eol(jn)-pbe(j)*eol(in)
Index: pdbread.f
===================================================================
--- pdbread.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ pdbread.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -21,8 +21,9 @@
 ! ......................................................
 
-      implicit real*8 (a-h,o-z)
-      implicit integer*4 (i-n)
-
       include 'INCP.H'
+      
+      double precision cor
+
+      integer ier, l, iendst, lunpdb, io, iopfil, iat, i
 
 ! -------------------------- input
@@ -47,5 +48,5 @@
         lunpdb = 99
       else
-        write (*,'(a)') 
+        write (*,'(a)')
      &    ' pdbread> empty file name to read pdb-structure'
 
@@ -56,5 +57,5 @@
 
       if (io.le.0) then
-        write (*,'(a,/,a)') 
+        write (*,'(a,/,a)')
      &  ' pdbread> ERROR opening file to read pdb-structure: ',
      &  pdbfil(1:iendst(pdbfil))
@@ -69,5 +70,5 @@
 
       if ( line(17:17).ne.' ' )  then
-        write (*,'(a,/,a,/,a,/,2a)') 
+        write (*,'(a,/,a,/,a,/,2a)')
      &  ' pdbread> found alternate atom location: ',
      &  '                !',
@@ -85,5 +86,5 @@
 
       if ((natp+1).gt.MXATP) then
-        write (*,'(a,i5,a,/,a)') 
+        write (*,'(a,i5,a,/,a)')
      &  ' pdbread>  >MXATP (',MXATP,') ATOM lines in PDB file ',
      &  pdbfil(1:iendst(pdbfil))
@@ -96,5 +97,5 @@
 
         if ((nchp+1).gt.MXCHP) then
-          write (*,'(a,i3,a,/,a)') 
+          write (*,'(a,i3,a,/,a)')
      &    ' pdbread>  >MXCHP (',MXCHP,') chains in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -105,5 +106,5 @@
 
         if ((nrsp+1).gt.MXRSP) then
-          write (*,'(a,i3,a,/,a)') 
+          write (*,'(a,i3,a,/,a)')
      &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -140,5 +141,5 @@
 
         if ((nrsp+1).gt.MXRSP) then
-          write (*,'(a,i3,a,/,a)') 
+          write (*,'(a,i3,a,/,a)')
      &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -171,5 +172,5 @@
       goto 1
 
-    2 write (*,'(a,/,a,/,2a)') 
+    2 write (*,'(a,/,a,/,2a)')
      &  ' pdbread> ERROR reading ATOM line ',
      &  line(:l),
@@ -194,5 +195,5 @@
       else
 
-        write (*,'(a,/,a)') 
+        write (*,'(a,/,a)')
      &  ' pdbread> NO atom coordinates selected from file ',
      &  pdbfil(1:iendst(pdbfil))
@@ -226,4 +227,9 @@
       include 'INCL.H'
       include 'INCP.H'
+
+      double precision dihedr, rm, av1, av2, rmsd, h, x, y, z
+
+      integer nml, nrs, nc, irb, ire, irs, i, ii, it, i1, i2, i3, i4, j1
+      integer j2, j3, j4, inew, j, ix
 
       character res*3
@@ -451,4 +457,6 @@
       include 'INCP.H'
 
+      integer irs, nml, i1, i2, iat, ix, i
+
       character*4 atm
 
@@ -476,5 +484,5 @@
               endif
             enddo
-               
+
 !            write(*,'(8a)') ' pdbvars> ',atm,' not found in '
 !     #       ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs)
@@ -484,5 +492,5 @@
     1     ixatp(iat)=ix
 
-        enddo   ! SMMP atoms of 'irs' 
+        enddo   ! SMMP atoms of 'irs'
       enddo   ! residues
 
@@ -494,4 +502,9 @@
       include 'INCL.H'
 
+      double precision tol, h1, h2, h3, x1, x2, x3, d, z1, z2, z3, y1
+      double precision y2, y3, yp1, yp2
+
+      integer i1, nml, i2, i3, i
+
       parameter (TOL = 1.d-12)
 
@@ -499,5 +512,5 @@
 ! -> determine global parameters: shifts dX,dY,dZ
 ! & angles alpha,beta,gamma [rad], put into 'gbpr'
-! 
+!
 ! CALLS: none
 !
@@ -547,5 +560,5 @@
       y2=z3*x1-z1*x3
       y3=z1*x2-z2*x1
-      
+
       if ( ( 1.d0 - abs(y3) ) .gt. TOL )  then       ! ============ |beta| < PI/2
 
Index: regul.f
===================================================================
--- regul.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ regul.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -4,5 +4,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -10,5 +10,5 @@
 ! **************************************************************
 
-      
+
       subroutine regul(nml, iter, nsteps, acc)
 
@@ -26,9 +26,13 @@
       include 'INCL.H'
       include 'INCP.H'
-      
+
 !f2py intent(in) nml
 !f2py intent(in) iter
 !f2py intent(in) nsteps
 !f2py intent(in) acc
+
+      double precision acc, rm, av1, av2, rmsd, dn
+
+      integer nsteps, nml, nrs, i, n, iter, it
 
       dimension rm(3,3),av1(3),av2(3)
@@ -60,5 +64,5 @@
 
 
-      do i = ivrml1(nml),nvrml(nml) 
+      do i = ivrml1(nml),nvrml(nml)
         fxvro(i) = fxvr(i)  ! save
         if (isrfvr(i)) fxvr(i) = .true.  ! fix vars. defined in ref.str.
@@ -77,5 +81,5 @@
       call cnteny(nml)
 
-      do i = ivrml1(nml),nvrml(nml) 
+      do i = ivrml1(nml),nvrml(nml)
         fxvr(i) = fxvro(i)  ! restore
       enddo  ! vars.
@@ -86,5 +90,5 @@
       wtey = 0.d0
 
-      n=iter            
+      n=iter
       dn=1.d0/dble(n)
 
Index: rgyr.f
===================================================================
--- rgyr.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ rgyr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -25,4 +25,11 @@
 !
       include 'INCL.H'
+      
+      double precision dn, dnp, dnh, dx, dxp, dy, dyp, dyh, dz, dzp, dzh
+      double precision d2, d2p, d2h, xi, yi, zi, dxh, rg2, rg2p, rg2h
+      double precision rgy, ee
+
+      integer nml, nml1, nml2, nat, i, i1, i2
+
 !f2py intent(in) nml
 !f2py intent(out) rgy
Index: rmsdfun.f
===================================================================
--- rmsdfun.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ rmsdfun.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -5,5 +5,5 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -26,4 +26,8 @@
 !
 ! Input
+      double precision xrf, yrf, zrf, rm, av1, av2, rssd
+
+      integer nml, ir1, ir2, ixat, isl
+
       dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
 ! Local
@@ -57,9 +61,14 @@
 !
 !      NB  uncomment last lines in 'fitmol' to return coordinates
-!          in 'x2' after fitting the ref. str. onto SMMP structure 
+!          in 'x2' after fitting the ref. str. onto SMMP structure
 ! ----------------------------------------------------------------
 
       include 'INCL.H'
       include 'INCP.H'
+      
+      double precision x1, x2, xrf, yrf, zrf, rm, av1, av2, rmsd
+
+      integer nml, nr, na, n, im, ir, ia, ir1, ir2, ix, ixat, isl
+
 
 !-------------------------------------------------------- input
@@ -147,10 +156,15 @@
 ! .......................................................
 !f2py intent(out) rmsd
-                    
+
       include 'INCL.H'
 !      implicit real*8 (a-h,o-z)
 !      implicit integer*4 (i-n)
- 
+
 ! ------------------------------------------- input/output
+      double precision dn, a1, a2, x1, x2, q, dm, dp, dxm, dym, dzm, dxp
+      double precision dyp, dzp, e, v, em, rmsd, rm
+
+      integer n, i, j, ndim4, im
+
       dimension x1(3,mxat),x2(3,mxat)
 ! -------------------------------------------------- local
@@ -270,6 +284,8 @@
 !f2py intent(out) d
 !f2py intent(out) v
+      integer nmax
+
       parameter (NMAX=500)
-      
+
       integer n,nrot,i,ip,iq,j
 
@@ -398,5 +414,5 @@
 !
 !------------------------------------------------------------------------------
-! Reads in pdb-file 'string' into INCP.H and initalizes 
+! Reads in pdb-file 'string' into INCP.H and initalizes
 ! the files that 'rmdsopt' needs to calculate the rmsd
 ! of a configuration with the pdb-configuration
@@ -409,6 +425,8 @@
       include 'INCP.H'
 
+      integer i, nml, ier
+
       character string*(*)
- 
+
       if(string.eq.'smmp') then
 !
@@ -428,5 +446,5 @@
 !
          call pdbread(string,ier)
-         if(ier.ne.0) stop 
+         if(ier.ne.0) stop
          call atixpdb(nml)
 !
Index: twister.f
===================================================================
--- twister.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ twister.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -47,5 +47,5 @@
 * May 03, 2000
 *
-* compile with: f90 -i32 -dp -o mt19937 mt19937.f 
+* compile with: f90 -i32 -dp -o mt19937 mt19937.f
 * Changes:
 *    Do not use UMASK as parameter due to restrictions of CF90
@@ -55,5 +55,5 @@
 *
 *    Now the T3E produces the same random number as a RS6000 running AIX
-*      
+*
 * from other programs sgrnd() must be called with integer*4
 *
@@ -64,5 +64,9 @@
       subroutine sgrnd(seed)
 *
-      implicit integer(a-z)
+      integer seed, tmaskb, tmaskc, tshftu, tshfts, tshftt
+      integer tshftl, y
+
+      integer N, N1, mti, mt, m, mata, lmask, mag01
+      integer kk
 *
 * Period parameters
@@ -88,4 +92,5 @@
 
       block data twbloks
+      integer n,n1,mt,mti
       parameter(N     =  624)
       parameter(N1    =  N+1)
@@ -99,5 +104,8 @@
       double precision function grnd()
 *
-      implicit integer(a-z)
+      integer tmaskb, tmaskc, tshftu, tshfts, tshftt, tshftl, y
+
+      integer n, n1, mti, mt, m, mata, lmask, mag01, kk
+
 *
 * Period parameters
Index: utilities.f
===================================================================
--- utilities.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ utilities.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -8,5 +8,5 @@
 
 !! Calculate the best way to distribute the work load across processors.
-!  It calculates the average number of interactions and then tries to 
+!  It calculates the average number of interactions and then tries to
 !  assign a number of interactions to each processor that is as close
 !  as possible to the average. The routine should be called once for
@@ -21,4 +21,6 @@
       include 'INCL.H'
 
+      integer i1ms, io, iv, i2ms, ms
+
       integer num_ppr, nml
       integer idxOfFirstVariable, idxOfLastVariable
@@ -26,5 +28,5 @@
       integer totalct, irank, itarget
       double precision ipps
-      
+
       if (nml.eq.0) then
          idxOfFirstVariable = ivrml1(1)
@@ -36,5 +38,5 @@
             end do
          end do
-      else 
+      else
          idxOfFirstVariable = ivrml1(nml)
          idxOfLastVariable = ivrml1(nml) + nvrml(nml) - 1
@@ -44,5 +46,5 @@
          end do
       end if
-      
+
       isum = 0
       do io = idxOfLastVariable, idxOfFirstVariable, - 1
@@ -52,10 +54,10 @@
          do ms = i1ms, i2ms
             do at = latms1(ms), latms2(ms)
-               do ivw=ivwat1(at),ivwat2(at)  
+               do ivw=ivwat1(at),ivwat2(at)
                   do j=lvwat1(ivw),lvwat2(ivw)
                      isum = isum + 1
                   end do
                end do
-               do i14=i14at1(at),i14at2(at)   
+               do i14=i14at1(at),i14at2(at)
                   isum = isum + 1
                end do
@@ -66,5 +68,5 @@
       write (*,*) "Total number of interactions:", isum
       write (*,*) "Average # of interactions per processor", ipps
-      
+
       totalct = 0
       irank = 1
@@ -72,9 +74,9 @@
       if (nml.eq.0) then
          i1ms = imsml1(ntlml)+ nmsml(ntlml)
-      else 
+      else
          i1ms = imsml1(nml)+ nmsml(nml)
       end if
       do io = idxOfLastVariable, idxOfFirstVariable, - 1
-         isum = 0 
+         isum = 0
          iv = iorvr(io)
          i2ms = i1ms - 1
@@ -83,10 +85,10 @@
             do at = latms1(ms), latms2(ms)
                atct = atct + 1
-               do ivw=ivwat1(at),ivwat2(at)  
+               do ivw=ivwat1(at),ivwat2(at)
                   do j=lvwat1(ivw),lvwat2(ivw)
                      isum = isum + 1
                   end do
                end do
-               do i14=i14at1(at),i14at2(at)   
+               do i14=i14at1(at),i14at2(at)
                   isum = isum + 1
                end do
@@ -113,5 +115,5 @@
 
       end subroutine distributeWorkLoad
-      
+
 !-----------------------------------------------------------------------
 !     The function fileNameMP takes a template of a file name in the
@@ -128,10 +130,10 @@
 !     \endcode
 !     will output base_0011.dat.
-!     
+!
 !     @param base the base file name, e.g., base_0000.dat.
 !     @param i1 index of the first character that may be replaced
 !     @param i2 index of the last character that may be replaced
 !     @param rank the number that should be inserted into the file name.
-!     
+!
 !     @return file name for rank
 !-----------------------------------------------------------------------
@@ -165,5 +167,5 @@
          write(fileNameMP(i2-5:i2), '(i6)') rank
       endif
-      end function fileNameMP 
+      end function fileNameMP
 !     End fileNameMP
 
@@ -172,23 +174,25 @@
 !     Add messages to log. This routine takes the log (debugging) mes-
 !     sages and writes them to the log file if the log level is less or
-!     equal to the maximum log level given by the global variable 
+!     equal to the maximum log level given by the global variable
 !     MAXLOGLEVEL.
 !
 !     @author Jan H. Meinke
 !
-!     @param loglevel level at which this message should be added to 
+!     @param loglevel level at which this message should be added to
 !            the log.
 !     @param message message to be written to the log.
-!     @param rank global rank of this node if running an MPI job zero 
+!     @param rank global rank of this node if running an MPI job zero
 !            otherwise.
 !----------------------------------------------------------------------
       subroutine addLogMessage(loglevel, message, rank)
-      
+
+      integer maxloglevel, logfileunit
+
          integer :: loglevel, rank
          character(LEN=*) :: message
-         
+
          if (loglevel <= MAXLOGLEVEL) then
             write(LOGFILEUNIT, *) message
          end if
-         
+
       end subroutine addLogMessage
Index: zimmer.f
===================================================================
--- zimmer.f	(revision 6650a56b8f862383b669bf41c64d328acf404e72)
+++ zimmer.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
@@ -40,4 +40,7 @@
 !f2py intent(in) nresi
       character*1 zim
+      double precision xphi, xpsi
+
+      integer j, nresi, i, iv, nres, nursvr
 
      
