Index: enylun.f
===================================================================
--- enylun.f	(revision 87ddcf7a314d90070d607f0ec7ac0a700b843abc)
+++ enylun.f	(revision e40e335d8909b50cd784019735fcdb0c28ee19bb)
@@ -1,78 +1,53 @@
-! *******************************************************************
-! 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
-! file.
-!
-! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
-!                      Jan H. Meinke, Sandipan Mohanty
-!
-      subroutine set_local_constr(ires,frst,lst,root)
-      include 'INCL.H'
-      include 'incl_lund.h'
-      integer iat1, iat2
-
-          integer ires,frst,lst,root
-          do iat1=iCa(ires)+frst,iCa(ires)+lst
-             do iat2=iat1+1,iCa(ires)+lst
-                matcon(iat2-iat1,iat1)=0
-                matcon(iat1-iat2,iat2)=0
-             enddo
-          enddo
-          iat1=iCa(ires)+root
-          do iat2=iCa(ires)+root,iCa(ires)+lst
-                matcon(iat2-iat1,iat1)=0
-                matcon(iat1-iat2,iat2)=0
-          enddo
-      end subroutine set_local_constr
-
+c *******************************************************************
+c SMMP version of Anders Irback's force field, to be called the Lund
+c force field. This file contains the function enylun, which in turn 
+c calls all the terms in the energy function. The terms Bias (ebias), 
+c Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the 
+c Excluded volume (eexvol and eloexv) are also implemented in this 
+c file.
+c
+c Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
+c                      Jan H. Meinke, Sandipan Mohanty
+c
       subroutine init_lundff
       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
 
       print *,'initializing Lund forcefield'
-!     Some parameters in the Lund force field.
-!     The correspondence between internal energy scale and kcal/mol
-      eunit=1.3315d0
-!      eunit=1.0
-!     Bias
-      kbias=100.0d0*eunit
-!      print *,'Bias'
-!     Hydrogen bonds
-      epshb1=3.1d0*eunit
-      epshb2=2.0d0*eunit
-      sighb=2.0d0
-      cthb=4.5d0
+c     Some parameters in the Lund force field.
+c     The correspondence between internal energy scale and kcal/mol
+      eunit=1.3315
+c     Bias
+      kbias=100.0*eunit
+c      print *,'Bias'
+c     Hydrogen bonds
+      epshb1=3.1*eunit
+      epshb2=2.0*eunit
+      sighb=2.0
+      cthb=4.5
       cthb2=cthb*cthb
-      powa=0.5d0
-      powb=0.5d0
-      blhb=-30.0d0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
+      powa=0.5
+      powb=0.5
+      blhb=-30.0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
       alhb=-(5*((sighb/cthb)**12)-6*((sighb/cthb)**10))-blhb*cthb2
       sighb2=sighb*sighb
-      cdon=1.0d0
-      cacc=(1.0d0/1.23d0)**powb
-      csacc=(1.0d0/1.25d0)**powb
-!      print *,'Hydrogen bonds'
-!     Hydrophobicity
-!      print *,'Hydrophobicity with nhptyp = ',nhptyp
-
-      hpstrg(1)=0.0d0*eunit
-      hpstrg(2)=0.1d0*eunit
-      hpstrg(3)=0.1d0*eunit
-      hpstrg(4)=0.1d0*eunit
-      hpstrg(5)=0.9d0*eunit
-      hpstrg(6)=2.8d0*eunit
-      hpstrg(7)=0.1d0*eunit
-      hpstrg(8)=2.8d0*eunit
-      hpstrg(9)=3.2d0*eunit
+      cdon=1.0
+      cacc=(1.0/1.23)**powb
+      csacc=(1.0/1.25)**powb
+c      print *,'Hydrogen bonds'
+c     Hydrophobicity
+c      print *,'Hydrophobicity with nhptyp = ',nhptyp
+
+      hpstrg(1)=0.0*eunit
+      hpstrg(2)=0.1*eunit
+      hpstrg(3)=0.1*eunit
+      hpstrg(4)=0.1*eunit
+      hpstrg(5)=0.9*eunit
+      hpstrg(6)=2.8*eunit
+      hpstrg(7)=0.1*eunit
+      hpstrg(8)=2.8*eunit
+      hpstrg(9)=3.2*eunit
 
       do i=1,mxrs
@@ -86,13 +61,12 @@
          call tolost(mynm)
          if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
-     &           .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
-     &           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
+     #           .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
+     #           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
             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
@@ -112,5 +86,5 @@
             ihpat(i,2)=iCa(i)+5
             ihpat(i,3)=iCa(i)+7
-            ihpat(i,4)=iCa(i)+11
+            ihpat(i,3)=iCa(i)+11
          else if (mynm.eq.'ile') then
             nhpat(i)=4
@@ -118,5 +92,5 @@
             ihpat(i,2)=iCa(i)+4
             ihpat(i,3)=iCa(i)+8
-            ihpat(i,4)=iCa(i)+11
+            ihpat(i,3)=iCa(i)+11
          else if (mynm.eq.'met') then
             nhpat(i)=4
@@ -124,5 +98,5 @@
             ihpat(i,2)=iCa(i)+5
             ihpat(i,3)=iCa(i)+8
-            ihpat(i,4)=iCa(i)+9
+            ihpat(i,3)=iCa(i)+9
          else if (mynm.eq.'phe') then
             nhpat(i)=6
@@ -130,7 +104,7 @@
             ihpat(i,2)=iCa(i)+6
             ihpat(i,3)=iCa(i)+8
-            ihpat(i,4)=iCa(i)+10
-            ihpat(i,5)=iCa(i)+12
-            ihpat(i,6)=iCa(i)+14
+            ihpat(i,3)=iCa(i)+10
+            ihpat(i,3)=iCa(i)+12
+            ihpat(i,3)=iCa(i)+14
          else if (mynm.eq.'tyr') then
             nhpat(i)=6
@@ -138,7 +112,7 @@
             ihpat(i,2)=iCa(i)+6
             ihpat(i,3)=iCa(i)+8
-            ihpat(i,4)=iCa(i)+10
-            ihpat(i,5)=iCa(i)+13
-            ihpat(i,6)=iCa(i)+15
+            ihpat(i,3)=iCa(i)+10
+            ihpat(i,3)=iCa(i)+13
+            ihpat(i,3)=iCa(i)+15
          else if (mynm.eq.'trp') then
             nhpat(i)=6
@@ -146,45 +120,45 @@
             ihpat(i,2)=iCa(i)+11
             ihpat(i,3)=iCa(i)+13
-            ihpat(i,4)=iCa(i)+15
-            ihpat(i,5)=iCa(i)+17
-            ihpat(i,6)=iCa(i)+19
-         endif
-      enddo
-!      print *,'Hydrophobicity'
-
-!     Excluded volume and local pair excluded volume terms
+            ihpat(i,3)=iCa(i)+15
+            ihpat(i,3)=iCa(i)+17
+            ihpat(i,3)=iCa(i)+19
+         endif
+      enddo
+c      print *,'Hydrophobicity'
+
+c     Excluded volume and local pair excluded volume terms
       exvk=0.1*eunit
       exvcut=4.3
       exvcut2=exvcut*exvcut
 
-      sigsa(1)=1.0d0  ! hydrogen
-      sigsa(2)=1.0d0  ! hydrogen
-      sigsa(3)=1.0d0  ! hydrogen
-      sigsa(4)=1.0d0  ! hydrogen
-      sigsa(5)=1.0d0  ! hydrogen
-      sigsa(6)=1.0d0  ! hydrogen
-      sigsa(7)=1.75d0 ! carbon
-      sigsa(8)=1.75d0 ! carbon
-      sigsa(9)=1.75d0 ! carbon
-      sigsa(10)=1.42d0 ! oxygen
-      sigsa(11)=1.42d0 ! oxygen
-      sigsa(12)=1.42d0 ! oxygen
-      sigsa(13)=1.55d0 ! nitrogen
-      sigsa(14)=1.55d0 ! nitrogen
-      sigsa(15)=1.55d0 ! nitrogen
-      sigsa(16)=1.77d0 ! sulfur
-      sigsa(17)=1.0d0  ! hydrogen
-      sigsa(18)=1.75d0 ! carbon
+      sigsa(1)=1.0  ! hydrogen
+      sigsa(2)=1.0  ! hydrogen
+      sigsa(3)=1.0  ! hydrogen
+      sigsa(4)=1.0  ! hydrogen
+      sigsa(5)=1.0  ! hydrogen
+      sigsa(6)=1.0  ! hydrogen
+      sigsa(7)=1.75 ! carbon
+      sigsa(8)=1.75 ! carbon
+      sigsa(9)=1.75 ! carbon
+      sigsa(10)=1.42 ! oxygen
+      sigsa(11)=1.42 ! oxygen
+      sigsa(12)=1.42 ! oxygen
+      sigsa(13)=1.55 ! nitrogen
+      sigsa(14)=1.55 ! nitrogen
+      sigsa(15)=1.55 ! nitrogen
+      sigsa(16)=1.77 ! sulfur
+      sigsa(17)=1.0  ! hydrogen
+      sigsa(18)=1.75 ! carbon
 
       do i=1,mxtyat
          do j=1,mxtyat
             sig2lcp(i,j)=(sigsa(i)+sigsa(j))**2
-            asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0d0)
-            bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0d0)/exvcut2
-         enddo
-      enddo
-!      print *,'Local pair excluded volume constants'
-
-      exvlam=0.75d0
+            asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0)
+            bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0)/exvcut2
+         enddo
+      enddo
+c      print *,'Local pair excluded volume constants'
+
+      exvlam=0.75
       exvcutg=exvcut*exvlam
       exvcutg2=exvcutg*exvcutg
@@ -197,16 +171,16 @@
          enddo
       enddo
-!      print *,'General excluded volume constants'
-
-!     Initialization of the connections matrix matcon(i,j). The index
-!     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
-!     molecule made of natural amino acids, atoms with indices separated
-!     by more than 35 can not be connected by three covalent bonds.
+c      print *,'General excluded volume constants'
+
+c     Initialization of the connections matrix matcon(i,j). The index
+c     i runs from -mxconr to +mxconr, and j from 1 to mxat. 
+c     matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed
+c                      = 2, if atoms i1 and i2 are separated by 3 covalent
+c                           bonds and their distance can change
+c                      = 1, for all other pairs
+c     if abs(i2-i1) > mxconr, the atoms are assumed to be separated by 
+c     many bonds, and with no restriction on their distances. On a protein 
+c     molecule made of natural amino acids, atoms with indices separated
+c     by more than 35 can not be connected by three covalent bonds.
 
       do i=1,mxat
@@ -216,5 +190,5 @@
          matcon(0,i)=0
       enddo
-!     continued...
+c     continued...
       do iml=1,ntlml
          do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
@@ -235,5 +209,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 !
@@ -250,17 +224,14 @@
          enddo
 
-!         print *,'going to initialize connections for first residue'
-!         print *,'iN,iCa,iC =',iN(irsml1(iml)),
-!     #        iCa(irsml1(iml)),iC(irsml1(iml))
+c         print *,'going to initialize connections for first residue'
+c         print *,'iN,iCa,iC =',iN(irsml1(iml)),
+c     #        iCa(irsml1(iml)),iC(irsml1(iml))
          do iat1=iN(irsml1(iml))+1,iCa(irsml1(iml))-1
-!            print *,'connections for iat1 = ',iat1
+c            print *,'connections for iat1 = ',iat1
             matcon(iat1-iN(irsml1(iml)),iN(irsml1(iml)))=0
             matcon(iN(irsml1(iml))-iat1,iat1)=0
             matcon(iat1-iCa(irsml1(iml)),iCa(irsml1(iml)))=0
             matcon(iCa(irsml1(iml))-iat1,iat1)=0
-            do iat2=iat1,iCa(irsml1(iml))-1
-                matcon(iat2-iat1,iat1)=0
-                matcon(iat1-iat2,iat2)=0
-            enddo
+
             matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
             matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
@@ -271,12 +242,12 @@
          enddo
 
-!     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.
+c     Below: for certain residues, some atoms separated by 3 or more bonds
+c     do not change distance. So, the connection matrix term for such pairs
+c     should be zero. 
 
          do irs=irsml1(iml),irsml2(iml)
             if (irs.eq.irsml1(iml)) then
                iatoff=1
-            else
+            else 
                iatoff=0
             endif
@@ -288,58 +259,74 @@
             mynm=seq(irs)
             call tolost(mynm)
-            if ((mynm.eq.'pro')) then
+            if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
+     #              .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
+     #              .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
                prlvr=.true.     ! residue i is a proline variant
-            else
+            else 
                prlvr=.false.
+            endif 
+
+            if ((mynm.eq.'arg').or.(mynm.eq.'arg+')) then
+               do iat1=iatoff+iatrs1(irs)+13,iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               enddo
+            else if ((mynm.eq.'his').or.(mynm.eq.'hise')
+     #              .or.(mynm.eq.'hisd').or.(mynm.eq.'his+')) then
+               do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               enddo
+            else if (mynm.eq.'phe') then 
+               do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               enddo
+            else if (mynm.eq.'tyr') then
+               do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     if (iat2.ne.iatrs1(irs)+15) then
+                        matcon(iat2-iat1,iat1)=0
+                        matcon(iat1-iat2,iat2)=0
+                     endif
+                  enddo
+               enddo
+            else if (mynm.eq.'trp') then
+               do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               enddo
+            else if (prlvr) then
+c           Proline. Many more distances are fixed because of the fixed 
+c           phi angle
+               do iat1=iatoff+iatrs1(irs),iatrs2(irs)-2-iatmrg
+                  do iat2=iat1+1,iatrs2(irs)-2-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               enddo
+c           distances to the C' atom of the previous residue are also fixed
+               if (irs.ne.irsml1(iml)) then
+                  iat1=iowat(iatrs1(irs))
+                  do iat2=iatrs1(irs),iatrs2(irs)-1-iatmrg
+                     matcon(iat2-iat1,iat1)=0
+                     matcon(iat1-iat2,iat2)=0
+                  enddo
+               endif
             endif
-            if ((mynm.eq.'asn')) then
-                call set_local_constr(irs,5,9,2)
-            else if ((mynm.eq.'gln')) then
-                call set_local_constr(irs,8,12,5)
-            else if ((mynm.eq.'arg')) then
-                call set_local_constr(irs,11,19,8)
-            else if ((mynm.eq.'his')) then
-                call set_local_constr(irs,5,12,2)
-            else if (mynm.eq.'phe') then
-                call set_local_constr(irs,5,15,2)
-            else if (mynm.eq.'tyr') then
-                call set_local_constr(irs,5,16,2)
-                iat1=iCa(irs)+12 ! H_h
-                iat2=iCa(irs)+13 ! C_e2
-                iat3=iCa(irs)+8  ! C_e1
-                matcon(iat2-iat1,iat1)=2
-                matcon(iat1-iat2,iat2)=2
-                matcon(iat3-iat1,iat1)=2
-                matcon(iat1-iat3,iat3)=2
-            else if (mynm.eq.'trp') then
-                call set_local_constr(irs,5,19,2)
-            else if (prlvr) then
-!           Proline. Many more distances are fixed because of the fixed
-!           phi angle
-                call set_local_constr(irs,-3,11,-3)
-                matcon(iCa(irs-1)-iCa(irs)-8,iCa(irs)+8)=0
-                matcon(iCa(irs)+8-iCa(irs-1),iCa(irs-1))=0
-            endif
-         enddo
-
-!        Distances fixed because of the constant omega angle
-         do irs=irsml1(iml),irsml2(iml)-1
-            matcon(iCa(irs+1)-iC(irs)-1,iC(irs)+1)=0 ! O_i -- Ca_{i+1}
-            matcon(-iCa(irs+1)+iC(irs)+1,iCa(irs+1))=0 ! O_i -- Ca_{i+1}
-
-            matcon(iN(irs+1)-iC(irs),iC(irs)+1)=0    ! O_i -- H_{i+1}
-            matcon(-iN(irs+1)+iC(irs),iN(irs+1)+1)=0    ! O_i -- H_{i+1}
-
-            matcon(iCa(irs+1)-iCa(irs),iCa(irs))=0   ! Ca_i -- Ca_{i+1}
-            matcon(-iCa(irs+1)+iCa(irs),iCa(irs+1))=0   ! Ca_i -- Ca_{i+1}
-
-            matcon(iN(irs+1)+1-iCa(irs),iCa(irs))=0  ! Ca_i -- H_{i+1}
-            matcon(-iN(irs+1)-1+iCa(irs),iN(irs+1)+1)=0  ! Ca_i -- H_{i+1}
-         enddo
-
-      enddo
-!     finished initializing matrix conmat
-!      print *,'Connections matrix'
-!     Local pair excluded volume
+         enddo
+      enddo
+c     finished initializing matrix conmat
+c      print *,'Connections matrix'
+
+c     Local pair excluded volume
       do i=1,mxml
          ilpst(i)=1
@@ -355,5 +342,5 @@
             do iat2=iat1+1,iatrs2(irsml2(iml))
                if ((iat2-iat1.le.mxconr).and.
-     &                 matcon(iat2-iat1,iat1).eq.2) then
+     #                 matcon(iat2-iat1,iat1).eq.2) then
                   ilp=ilp+1
                   lcp1(ilp)=iat1
@@ -364,13 +351,13 @@
 
          ilpnd(iml)=ilp
-         if (iml.lt.ntlml) then
+         if (iml.lt.ntlml) then 
             ilpst(iml+1)=ilp+1
          endif
-         print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
-!        print *,'local pair list'
+c         print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
+c         print *,'local pair list'
          do lci=ilpst(iml),ilpnd(iml)
             iat1=lcp1(lci)
             iat2=lcp2(lci)
-!            print *,lci,iat1,iat2,ityat(iat1),ityat(iat2)
+c            print *,lci,iat1,iat2,matcon(iat2-iat1,iat1)
          enddo
       enddo
@@ -381,6 +368,4 @@
       integer function ihptype(iaa)
       include 'INCL.H'
-      integer iaa, ityp
-
       character mynm*4
       mynm=seq(iaa)
@@ -390,12 +375,14 @@
          ityp=1
       else if ((mynm.eq.'val').or.(mynm.eq.'leu').or.(mynm.eq.'ile')
-     &        .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then
+     #        .or.(mynm.eq.'met').or.(mynm.eq.'pro').or.(mynm.eq.'cpro')
+     #        .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
+     #        .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
          ityp=2
-      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
-     &        then
+      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp')) 
+     #        then
          ityp=3
       endif
       ihptype=ityp
-      return
+      return 
       end
 
@@ -403,11 +390,7 @@
       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
@@ -425,13 +408,9 @@
       return
       end
-!     Evaluates backbone backbone hydrogen bond strength for residues
-!     i and j, taking the donor from residue i and acceptor from residue j
+c     Evaluates backbone backbone hydrogen bond strength for residues 
+c     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
@@ -444,7 +423,7 @@
       dz=zat(a1)-zat(d2)
       r2=dx*dx+dy*dy+dz*dz
-      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
+      if (r2.gt.cthb2) then 
+c         print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2
+c         print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb
          ehbmmrs=0
          return
@@ -453,10 +432,10 @@
       cb=(xat(a2)-xat(a1))*dx+(yat(a2)-yat(a1))*dy+(zat(a2)-zat(a1))*dz
       if (powa.gt.0.and.ca.le.0) then
-!         print *,'hbmm, returning 0 because of angle a'
+c         print *,'hbmm, returning 0 because of angle a'
          ehbmmrs=0
          return
       endif
       if (powb.gt.0.and.cb.le.0) then
-!         print *,'hbmm, returning 0 because of angle b'
+c         print *,'hbmm, returning 0 because of angle b'
          ehbmmrs=0
          return
@@ -465,26 +444,15 @@
       r4=r6*r6
       r6=r6*r4
-      rdon2=(xat(d2)-xat(d1))**2+(yat(d2)-yat(d1))**2+
-     & (zat(d2)-zat(d1))**2
-      racc2=(xat(a2)-xat(a1))**2+(yat(a2)-yat(a1))**2+
-     & (zat(a2)-zat(a1))**2
-      evlu=((ca*ca/(r2*rdon2))**(0.5*powa))*
-     & ((cb*cb/(r2*racc2))**(0.5*powb))
+      evlu=((ca*ca/r2)**(0.5*powa))*((cb*cb/r2)**(0.5*powb))
       evlu=evlu*(r6*(5*r6-6*r4)+alhb+blhb*r2)
-!      print *,'found hbmm contribution ',i,j,epshb1,evlu
+c      print *,'found hbmm contribution ',evlu
       ehbmmrs=epshb1*evlu
       return
       end
       real*8 function enylun(nml)
-!     nml = 1 .. ntlml. No provision exists to handle out of range values
-!     for nml inside this function.
+c     nml = 1 .. ntlml. No provision exists to handle out of range values
+c     for nml inside this function.
       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
@@ -494,50 +462,47 @@
       eyvr=0.0   ! Local pair excluded volume, in a sense a variable potential
       eyvw=0.0   ! atom-atom repulsion, excluded volume
-!     atom-atom repulsion is calculated on a system wide basis, instead of
-!     molecule by molecule for efficiency. Look into function exvlun.
+c     atom-atom repulsion is calculated on a system wide basis, instead of
+c     molecule by molecule for efficiency. Look into function exvlun.
 
       istres=irsml1(nml)
       indres=irsml2(nml)
 
-!     First, all terms that can be calculated on a residue by residue basis
+c     First, all terms that can be calculated on a residue by residue basis
       do i=istres,indres
          mynm=seq(i)
          call tolost(mynm)
          if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
-     &           .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
-     &           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
+     #           .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
+     #           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
             prlvr=.true.        ! residue i is a proline variant
-         else
+         else 
             prlvr=.false.
-         endif
-
-!     Bias, or local electrostatic term. Excluded from the list are
-!     residues at the ends of the chain, glycine and all proline variants
+         endif 
+
+c     Bias, or local electrostatic term. Excluded from the list are
+c     residues at the ends of the chain, glycine and all proline variants
          if ((i.ne.istres).and.(i.ne.indres).and.
-     &           .not.prlvr.and.mynm.ne.'gly') then
+     #           .not.prlvr.and.mynm.ne.'gly') then
             eyel=eyel+ebiasrs(i)
          endif
-!     Backbone--backbone hydrogen bonds
+c     Backbone--backbone hydrogen bonds
          shbm1=1.0
          shbm2=1.0
          if ((i.eq.istres).or.(i.eq.indres)) shbm1=0.5
-!     Residue i contributes the donor, and j, the acceptor, so both i and
-!     j run over the whole set of amino acids.
-!     No terms for residue i, if it is a proline variant.
-         if (.not.prlvr) then
+c     Residue i contributes the donor, and j, the acceptor, so both i and
+c     j run over the whole set of amino acids.
+c     No terms for residue i, if it is a proline variant.
+         if (.not.prlvr) then  
             do j=istres,indres
-               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
+               if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
+               etmp=ehbmmrs(i,j)
+               eyhb=eyhb+shbm1*shbm2*etmp
             enddo
          endif
-!     Hydrophobicity, only if residue i is hydrophobic to start with
+c     Hydrophobicity, only if residue i is hydrophobic to start with
          ihpi=ihptype(i)
-         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.
+         if (ihpi.ge.0) then 
+c        Unlike hydrogen bonds, the hydrophobicity potential is symmetric
+c        in i and j. So, the loop for j runs from i+1 to the end.
 
             do j=i+1,indres
@@ -547,16 +512,15 @@
                   if (j.eq.(i+1)) etmp=0
                   if (j.eq.(i+2)) etmp=0.5*etmp
-!                  print *, 'hydrophobicity contribution ',etmp,i,j
                   eysl=eysl+etmp
-               endif
+               endif 
             enddo
          endif
       enddo
 
-!     Terms that are not calculated residue by residue ...
-
-!     Local pair or third-neighbour excluded volume
-!     Numerically this is normally the term with largest positive
-!     contribution to the energy in an equilibrated stystem.
+c     Terms that are not calculated residue by residue ...
+
+c     Local pair or third-neighbour excluded volume
+c     Numerically this is normally the term with largest positive
+c     contribution to the energy in an equilibrated stystem. 
 
       i1=ilpst(nml)
@@ -574,15 +538,14 @@
          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)
-            if (etmp1.ge.0) then
-!               print *,'local pair contribution ',iat1,iat2,iatt1,
-!     & iatt2,sqrt(sig2lcp(iatt1,iatt2)),etmp1
+            if (etmp1.ge.500) then
+               print *,'local pair contribution ',iat1,iat2,etmp1
             endif
             etmp=etmp+etmp1
          endif
-!         print *,'pair : ',iat1,iat2,' contribution ',etmp1
-!         print *,exvcut2,r2
+c         print *,'pair : ',iat1,iat2,' contribution ',etmp1 
+c         print *,exvcut2,r2
       enddo
       eyvr=exvk*etmp
@@ -602,14 +565,10 @@
       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)
 
       b2=20.25
       a2=12.25
-!      ihp1=ihptype(i1)
-!      ihp2=ihptype(i2)
+c      ihp1=ihptype(i1)
+c      ihp2=ihptype(i2)
       if ((ihp1.le.0).or.(ihp2.le.0)) then
          ehp=0.0
@@ -633,17 +592,15 @@
          enddo
       enddo
-      ssum=0
+      sum=0
       do i=1,ni+nj
-         if (r2min(i).le.b2) then
-            if (r2min(i).lt.a2) then
-               ssum=ssum+1
-            else
-               ssum=ssum+(b2-r2min(i))/(b2-a2)
+         if (r2min(i).le.b2) then 
+            if (r2min(i).lt.a2) then 
+               sum=sum+1
+            else 
+               sum=sum+(b2-r2min(i))/(b2-a2)
             endif
          endif
-!         hpstrgth=hpstrg((ihp1-1)*nhptyp+ihp2)
-!         print *, 'hp diagnosis ',ni,nj,ssum/(ni+nj),hpstrgth
-      enddo
-      ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*ssum/(ni+nj)
+      enddo
+      ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*sum/(ni+nj)
       return
       end
@@ -652,27 +609,16 @@
       include 'INCL.H'
       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
-!     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
-!     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
-
+c     For multi-chain systems it makes little sense to split the calculation
+c     of this term into an 'interaction part' and a contribution from 
+c     individual molecules. So, normally this should always be called with
+c     argument nml=0. Only for diagnostic reasons, you might want to find
+c     the contribution from one molecule in a multi-chain system assuming 
+c     there was no other molecule.
       dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
       dimension icell(mxat)
-      integer :: jx1, jy1, jz1
-
-      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))
@@ -680,16 +626,16 @@
 
       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
-!     unlike the accessible surface calculations, this term is symmetric
-!     in any two participating atoms. So, the part after the assignment
-!     of cells differs even in the structure of the code.
+c     The beginning part of this implementation is very similar to the 
+c     assignment of cells to the atoms during calculation of solvent 
+c     accessible surface area. So, much of that part is similar. But 
+c     unlike the accessible surface calculations, this term is symmetric
+c     in any two participating atoms. So, the part after the assignment
+c     of cells differs even in the structure of the code.
 
       do i=1,mxcell
          incell(i)=0
       enddo
-!      print *,'evaluating general excluded volume :',istat,',',indat
-!     Find minimal containing box
+c      print *,'evaluating general excluded volume :',istat,',',indat
+c     Find minimal containing box
       xmin=xat(istat)
       ymin=yat(istat)
@@ -720,5 +666,5 @@
       sizey=ymax-ymin
       sizez=zmax-zmin
-!     Number of cells along each directions that fit into the box.
+c     Number of cells along each directions that fit into the box.
       ndx=int(sizex/exvcutg)+1
       ndy=int(sizey/exvcutg)+1
@@ -727,12 +673,12 @@
       nxy=ndx*ndy
       ncell=nxy*ndz
-!      print *,'Number of cells along x,y,z = ',ndx,',',ndy,',',ndz
+c      print *,'Number of cells along x,y,z = ',ndx,',',ndy,',',ndz
       if (ncell.ge.mxcell) then
          print *,'exvlun> required number of cells',ncell,
-     &        ' exceeded the limit ',mxcell
+     #        ' exceeded the limit ',mxcell
          print *,'recompile with a higher mxcell.'
          stop
       endif
-!     Expand box to contain an integral number of cells along each direction
+c     Expand box to contain an integral number of cells along each direction
       shiftx=(dble(ndx)*exvcutg-sizex)/2.0
       shifty=(dble(ndy)*exvcutg-sizey)/2.0
@@ -745,6 +691,6 @@
       zmax=zmax+shiftz
 
-!     Set occupied cells to zero. Note that the maximum number of occupied
-!     cells is the number of atoms in the system.
+c     Set occupied cells to zero. Note that the maximum number of occupied
+c     cells is the number of atoms in the system.
       nocccl=0
       do i=1,mxat
@@ -752,5 +698,5 @@
       enddo
 
-!     Put atoms in cells
+c     Put atoms in cells
       do j=istat,indat
          mx=min(int(max((xat(j)-xmin)/exvcutg,0.0d0)),ndx-1)
@@ -759,10 +705,10 @@
          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
-!           previously unoccupied cell
+         else 
+            if (incell(icellj).eq.0) then 
+c           previously unoccupied cell
                nocccl=nocccl+1
                locccl(nocccl)=icellj
@@ -771,11 +717,11 @@
          endif
       enddo
-!      print *,'finished assigning cells. nocccl = ',nocccl
-!     Cummulative occupancy of i'th cell
+c      print *,'finished assigning cells. nocccl = ',nocccl
+c     Cummulative occupancy of i'th cell
       do i=1,ncell
          incell(i+1)=incell(i+1)+incell(i)
       enddo
-!      print *,'finished making cumulative cell sums'
-!     Sorting atoms by their cell index
+c      print *,'finished making cumulative cell sums'
+c     Sorting atoms by their cell index
       do i=istat,indat
          j=icell(i)
@@ -784,46 +730,29 @@
          incell(j)=jj-1
       enddo
-!      print *,'sorted atoms by cell index'
+c      print *,'sorted atoms by cell index'
       etmp=0.0
       do icl=1,nocccl
-!     loop through occupied cells
+c     loop through occupied cells
          lcell=locccl(icl)
          ix=mod(lcell-1,ndx)
-         iz = (lcell - 1) / nxy
-         iy = ((lcell - 1) - iz * nxy) / ndx
-
+         iy=(mod(lcell-1,nxy)-ix)/ndx
+         iz=(lcell-1-ix-ndx*iy)/nxy
 c         print *,'icl=',icl,'absolute index of cell = ',lcell
 c         print *,'iz,iy,ix = ',iz,iy,ix
 c     find all atoms in current cell and all its forward-going neighbours
-         nex=ix+1
-         ney=iy+1
-         nez=iz+1
+         nex=min(ix+1,ndx-1)
+         ney=min(iy+1,ndy-1)
+         nez=min(iz+1,ndz-1)
          nsame=0
          nngbr=0
          do jx=ix,nex
-            if (jx.ge.ndx) then
-               jx1=0
-            else
-               jx1=jx
-            endif
             do jy=iy,ney
-               if (jy.ge.ndy) then
-                  jy1=0
-               else
-                  jy1=jy
-               endif
                do jz=iz,nez
-                  if (jz.ge.ndz) then
-                     jz1=0
-                  else
-                     jz1=jz
-                  endif
-                  jcl=jx1 + ndx*jy1 + nxy*jz1 + 1
-!                  write (logString, *)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
+                  jcl=jx+ndx*jy+nxy*jz+1
                   do ii=incell(jcl)+1,incell(jcl+1)
-!     count the total number of neighbours
+c                    count the total number of neighbours
                      nngbr=nngbr+1
                      if (jx.eq.ix.and.jy.eq.iy.and.jz.eq.iz) then
-!     count how many neighbours are from the same cell
+c                    count how many neighbours are from the same cell
                         nsame=nsame+1
                      endif
@@ -833,12 +762,10 @@
             enddo
          enddo
-!     A few more cells need to be searched, so that we cover 13 of the 26
-!     neighbouring cells.
-!        1
+c     A few more cells need to be searched, so that we cover 13 of the 26
+c     neighbouring cells. 
+c        1
          jx=ix+1
-         if (jx.ge.ndx) jx=0
          jy=iy
          jz=iz-1
-         if (jz.lt.0) jz=ndz-1
          jcl=jx+ndx*jy+nxy*jz+1
          do ii=incell(jcl)+1,incell(jcl+1)
@@ -846,10 +773,8 @@
             ngbr(nngbr)=isort(ii)
          enddo
-!        2
+c        2
          jx=ix
          jy=iy-1
-         if (jy.lt.0) jy =ndy-1
          jz=iz+1
-         if (jz.ge.ndz) jz=0
          jcl=jx+ndx*jy+nxy*jz+1
          do ii=incell(jcl)+1,incell(jcl+1)
@@ -857,9 +782,7 @@
             ngbr(nngbr)=isort(ii)
          enddo
-!        3
+c        3
          jx=ix-1
-         if (jx.lt.0)jx =ndx-1
          jy=iy+1
-         if (jy.ge.ndy) jy=0
          jz=iz
          jcl=jx+ndx*jy+nxy*jz+1
@@ -868,11 +791,8 @@
             ngbr(nngbr)=isort(ii)
          enddo
-!        4
+c        4
          jx=ix+1
-         if (jx.ge.ndx) jx=0
          jy=iy+1
-         if (jy.ge.ndy) jy=0
          jz=iz-1
-         if (jz.lt.0) jz=ndz-1
          jcl=jx+ndx*jy+nxy*jz+1
          do ii=incell(jcl)+1,incell(jcl+1)
@@ -880,11 +800,8 @@
             ngbr(nngbr)=isort(ii)
          enddo
-!        5
+c        5
          jx=ix+1
-         if (jx.ge.ndx) jx = 0
          jy=iy-1
-         if (jy.lt.0) jy=ndy-1
          jz=iz+1
-         if (jz.ge.ndz) jz=0
          jcl=jx+ndx*jy+nxy*jz+1
          do ii=incell(jcl)+1,incell(jcl+1)
@@ -892,11 +809,8 @@
             ngbr(nngbr)=isort(ii)
          enddo
-!        6
+c        6
          jx=ix+1
-         if (jx.ge.ndx) jx=0
          jy=iy-1
-         if (jy.lt.0) jy=ndy-1
          jz=iz-1
-         if (jz.lt.0) jz=ndy-1
          jcl=jx+ndx*jy+nxy*jz+1
          do ii=incell(jcl)+1,incell(jcl+1)
@@ -905,11 +819,11 @@
          enddo
 
-!         print *,'atoms in same cell ',nsame
-!         print *,'atoms in neighbouring cells ',nngbr
+c         print *,'atoms in same cell ',nsame
+c         print *,'atoms in neighbouring cells ',nngbr
          do i1=1,nsame
-!        Over all atoms from the original cell
+c        Over all atoms from the original cell
             iat1=ngbr(i1)
             do i2=i1,nngbr
-!           Over all atoms in the original+neighbouring cells
+c           Over all atoms in the original+neighbouring cells
                iat2=ngbr(i2)
                xij=xat(iat1)-xat(iat2)
@@ -918,7 +832,7 @@
                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
+     #                 matcon(iat2-iat1,iat1).eq.1) then
                      iatt1=ityat(iat1)
                      iatt2=ityat(iat2)
@@ -926,11 +840,11 @@
                      r6=r6*r6*r6
                      etmp1=r6*r6+asaexv(iatt1,iatt2)
-     &                    +bsaexv(iatt1,iatt2)*r2
+     #                    +bsaexv(iatt1,iatt2)*r2
                      etmp=etmp+etmp1
-!                     if (etmp1.ge.2000) then
-!                        print *,'contribution ',iat1,iat2,etmp1
-!                        call outpdb(1,'EXAMPLES/clash.pdb')
-!                        stop
-!                     endif
+c                     if (etmp1.ge.2000) then
+c                        print *,'contribution ',iat1,iat2,etmp1
+c                        call outpdb(1,'EXAMPLES/clash.pdb')
+c                        stop
+c                     endif
                   endif
                endif
@@ -938,16 +852,16 @@
          enddo
       enddo
-!      irs=1
-!      do iat=iatrs1(irs),iatrs2(irs)
-!         do j=-mxconr,mxconr
-!            print *,iat,j,':',matcon(j,iat)
-!         enddo
-!      enddo
-!      irs=irsml2(1)
-!      do iat=iatrs1(irs),iatrs2(irs)
-!         do j=-mxconr,mxconr
-!            print *,iat,j,':',matcon(j,iat)
-!         enddo
-!      enddo
+c      irs=1
+c      do iat=iatrs1(irs),iatrs2(irs)
+c         do j=-mxconr,mxconr
+c            print *,iat,j,':',matcon(j,iat)
+c         enddo
+c      enddo
+c      irs=irsml2(1)
+c      do iat=iatrs1(irs),iatrs2(irs)
+c         do j=-mxconr,mxconr
+c            print *,iat,j,':',matcon(j,iat)
+c         enddo
+c      enddo
 
       eyvw=exvk*etmp
@@ -957,11 +871,7 @@
 
       real*8 function exvbrfc()
-!     Brute force excluded volume evaluation
+c     Brute force excluded volume evaluation
       include 'INCL.H'
       include 'incl_lund.h'
-
-      double precision etmp, etmp1, xij, yij, zij, r2, r6
-
-      integer iat1, iat2, iatt1, iatt2
 
       etmp=0.0
@@ -975,7 +885,7 @@
             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
+     #                 matcon(iat2-iat1,iat1).eq.1) then
                   iatt1=ityat(iat1)
                   iatt2=ityat(iat2)
@@ -983,9 +893,20 @@
                   r6=r6*r6*r6
                   etmp1=r6*r6+asaexv(iatt1,iatt2)
-     &                 +bsaexv(iatt1,iatt2)*r2
+     #                 +bsaexv(iatt1,iatt2)*r2
                   etmp=etmp+etmp1
-               else
-!                  print *,'atoms ', iat1,' and ',iat2,' were close',
-!     #                 'but matcon is ',matcon(iat2-iat1,iat1)
+                  if (iat1.eq.43.and.iat2.eq.785) then
+                     print *,'contribution ',iat1,iat2,etmp1
+                     print *,'r2 = ',r2
+                     print *,'atom types : ',iatt1,iatt2
+                     print *,'sig2exv = ',sig2exv(iatt1,iatt2)
+                     print *,'asa = ',asaexv(iatt1,iatt2)
+                     print *,'bsa = ',bsaexv(iatt1,iatt2)
+
+c     call outpdb(1,'EXAMPLES/clash.pdb')
+c     stop
+                  endif
+               else 
+c                  print *,'atoms ', iat1,' and ',iat2,' were close',
+c     #                 'but matcon is ',matcon(iat2-iat1,iat1)
                endif
             endif
