subroutine parmread
C
C Read the parameters of the probability distributions of the virtual-bond
C valence angles and the side chains and energy parameters.
C
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'sizesclu.dat'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.TORSION'
      include 'COMMON.FFIELD'
      include 'COMMON.NAMES'
      include 'COMMON.SBRIDGE'
      include 'COMMON.SCCOR'
      include 'COMMON.SCROT'
      character*1 t1,t2,t3
      character*1 onelett(4) /"G","A","P","D"/
      logical lprint
      dimension blower(3,3,maxlob)
      double precision ip,mp
C
C Body
C
C Set LPRINT=.TRUE. for debugging
      dwa16=2.0d0**(1.0d0/6.0d0)
      lprint=.false.
      itypro=20
C Assign virtual-bond length
      vbl=3.8D0
      vblinv=1.0D0/vbl
      vblinv2=vblinv*vblinv
#ifdef CRYST_BOND
      read (ibond,*) vbldp0,akp
      do i=1,ntyp
        nbondterm(i)=1
        read (ibond,*) vbldsc0(1,i),aksc(1,i)
        dsc(i) = vbldsc0(1,i)
        if (i.eq.10) then
          dsc_inv(i)=0.0D0
        else
          dsc_inv(i)=1.0D0/dsc(i)
        endif
      enddo
#else
      read (ibond,*) ijunk,vbldp0,akp,rjunk
      do i=1,ntyp
        read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
     &   j=1,nbondterm(i))
        dsc(i) = vbldsc0(1,i)
        if (i.eq.10) then
          dsc_inv(i)=0.0D0
        else
          dsc_inv(i)=1.0D0/dsc(i)
        endif
      enddo
#endif
      if (lprint) then
        write(iout,'(/a/)')"Force constants virtual bonds:"
        write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
     &   'inertia','Pstok'
        write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
        do i=1,ntyp
          write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
     &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
          do j=2,nbondterm(i)
            write (iout,'(13x,3f10.5)')
     &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
          enddo
        enddo
      endif
#ifdef CRYST_THETA
C
C Read the parameters of the probability distribution/energy expression 
C of the virtual-bond valence angles theta
C
      do i=1,ntyp
        read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
        read (ithep,*) (polthet(j,i),j=0,3)
	read (ithep,*) (gthet(j,i),j=1,3)
	read (ithep,*) theta0(i),sig0(i),sigc0(i)
	sigc0(i)=sigc0(i)**2
      enddo
      close (ithep)
      if (lprint) then
c       write (iout,'(a)') 
c    &    'Parameters of the virtual-bond valence angles:'
c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
c    & '    ATHETA0   ','         A1   ','        A2    ',
c    & '        B1    ','         B2   '        
c       do i=1,ntyp
c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
c       enddo
c       write (iout,'(/a/9x,5a/79(1h-))') 
c    & 'Parameters of the expression for sigma(theta_c):',
c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
c    & '     ALPH3    ','    SIGMA0C   '        
c       do i=1,ntyp
c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
c    &      (polthet(j,i),j=0,3),sigc0(i) 
c       enddo
c       write (iout,'(/a/9x,5a/79(1h-))') 
c    & 'Parameters of the second gaussian:',
c    & '    THETA0    ','     SIGMA0   ','        G1    ',
c    & '        G2    ','         G3   '        
c       do i=1,ntyp
c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
c    &       sig0(i),(gthet(j,i),j=1,3)
c       enddo
	write (iout,'(a)') 
     &    'Parameters of the virtual-bond valence angles:'
        write (iout,'(/a/9x,5a/79(1h-))') 
     & 'Coefficients of expansion',
     & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
     & '   b1*10^1    ','    b2*10^1   '        
        do i=1,ntyp
          write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
     &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
        enddo
	write (iout,'(/a/9x,5a/79(1h-))') 
     & 'Parameters of the expression for sigma(theta_c):',
     & ' alpha0       ','  alph1       ',' alph2        ',
     & ' alhp3        ','   sigma0c    '        
	do i=1,ntyp
          write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
     &      (polthet(j,i),j=0,3),sigc0(i) 
	enddo
	write (iout,'(/a/9x,5a/79(1h-))') 
     & 'Parameters of the second gaussian:',
     & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
     & '        G2    ','   G3*10^1    '        
	do i=1,ntyp
          write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
     &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
	enddo
      endif
#else
C
C Read the parameters of Utheta determined from ab initio surfaces
C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
C
      read (ithep,*) nthetyp,ntheterm,ntheterm2,
     &  ntheterm3,nsingle,ndouble
      nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
      read (ithep,*) (ithetyp(i),i=1,ntyp1)
      do i=1,maxthetyp
        do j=1,maxthetyp
          do k=1,maxthetyp
            aa0thet(i,j,k)=0.0d0
            do l=1,ntheterm
              aathet(l,i,j,k)=0.0d0
            enddo
            do l=1,ntheterm2
              do m=1,nsingle
                bbthet(m,l,i,j,k)=0.0d0
                ccthet(m,l,i,j,k)=0.0d0
                ddthet(m,l,i,j,k)=0.0d0
                eethet(m,l,i,j,k)=0.0d0
              enddo
            enddo
            do l=1,ntheterm3
              do m=1,ndouble
                do mm=1,ndouble
                 ffthet(mm,m,l,i,j,k)=0.0d0
                 ggthet(mm,m,l,i,j,k)=0.0d0
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo
      do i=1,nthetyp
        do j=1,nthetyp
          do k=1,nthetyp
            read (ithep,'(3a)') res1,res2,res3
            read (ithep,*) aa0thet(i,j,k)
            read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
            read (ithep,*)
     &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
     &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
     &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
     &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
            read (ithep,*)
     &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
     &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
     &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
          enddo
        enddo
      enddo
C
C For dummy ends assign glycine-type coefficients of theta-only terms; the
C coefficients of theta-and-gamma-dependent terms are zero.
C
      do i=1,nthetyp
        do j=1,nthetyp
          do l=1,ntheterm
            aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
            aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
          enddo
          aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
          aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
        enddo
        do l=1,ntheterm
          aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
        enddo
        aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
      enddo
C
C Control printout of the coefficients of virtual-bond-angle potentials
C
      if (lprint) then
        write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
        do i=1,nthetyp+1
          do j=1,nthetyp+1
            do k=1,nthetyp+1
              write (iout,'(//4a)')
     &         'Type ',onelett(i),onelett(j),onelett(k)
              write (iout,'(//a,10x,a)') " l","a[l]"
              write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
              write (iout,'(i2,1pe15.5)')
     &           (l,aathet(l,i,j,k),l=1,ntheterm)
            do l=1,ntheterm2
              write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
     &          "b",l,"c",l,"d",l,"e",l
              do m=1,nsingle
                write (iout,'(i2,4(1pe15.5))') m,
     &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
     &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
              enddo
            enddo
            do l=1,ntheterm3
              write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
     &          "f+",l,"f-",l,"g+",l,"g-",l
              do m=2,ndouble
                do n=1,m-1
                  write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
     &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
     &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo
#ifdef AIX
      call flush_(iout)
#else
      call flush(iout)
#endif
      endif
#endif

#ifdef CRYST_SC
C
C Read the parameters of the probability distribution/energy expression
C of the side chains.
C
      do i=1,ntyp
	read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
        if (i.eq.10) then
          dsc_inv(i)=0.0D0
        else
          dsc_inv(i)=1.0D0/dsc(i)
        endif
	if (i.ne.10) then
        do j=1,nlob(i)
          do k=1,3
            do l=1,3
              blower(l,k,j)=0.0D0
            enddo
          enddo
        enddo  
	bsc(1,i)=0.0D0
        read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
	do j=2,nlob(i)
	  read (irotam,*) bsc(j,i)
	  read (irotam,*) (censc(k,j,i),k=1,3),
     &                                 ((blower(k,l,j),l=1,k),k=1,3)
        enddo
	do j=1,nlob(i)
	  do k=1,3
	    do l=1,k
	      akl=0.0D0
	      do m=1,3
		akl=akl+blower(k,m,j)*blower(l,m,j)
              enddo
	      gaussc(k,l,j,i)=akl
	      gaussc(l,k,j,i)=akl
            enddo
          enddo 
	enddo
	endif
      enddo
      close (irotam)
      if (lprint) then
	write (iout,'(/a)') 'Parameters of side-chain local geometry'
	do i=1,ntyp
	  nlobi=nlob(i)
          if (nlobi.gt.0) then
          write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
     &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
c          write (iout,'(a,f10.4,4(16x,f10.4))')
c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
           write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
     &                             'log h',(bsc(j,i),j=1,nlobi)
           write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
     &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
c          write (iout,'(a)')
c         do j=1,nlobi
c           ind=0
c           do k=1,3
c             do l=1,k
c              ind=ind+1
c              blower(k,l,j)=gaussc(ind,j,i)
c             enddo
c           enddo
c         enddo
	  do k=1,3
            write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
     &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
          enddo
	  endif
        enddo
      endif
#else
C
C Read scrot parameters for potentials determined from all-atom AM1 calculations
C added by Urszula Kozlowska 07/11/2007
C
      do i=1,ntyp
        read (irotam,*)
       if (i.eq.10) then
         read (irotam,*)
       else
         do j=1,65
           read(irotam,*) sc_parmin(j,i)
         enddo
       endif
      enddo
#endif
      close(irotam)
#ifdef CRYST_TOR
C
C Read torsional parameters in old format
C
      read (itorp,*) ntortyp,nterm_old
      write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
      read (itorp,*) (itortyp(i),i=1,ntyp)
      do i=1,ntortyp
	do j=1,ntortyp
	  read (itorp,'(a)')
	  do k=1,nterm_old
	    read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
          enddo
        enddo
      enddo
      close (itorp)
      if (lprint) then
	write (iout,'(/a/)') 'Torsional constants:'
	do i=1,ntortyp
	  do j=1,ntortyp
	    write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
	    write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
          enddo
        enddo
      endif
#else
C
C Read torsional parameters
C
      read (itorp,*) ntortyp
      read (itorp,*) (itortyp(i),i=1,ntyp)
      write (iout,*) 'ntortyp',ntortyp
      do i=1,ntortyp
	do j=1,ntortyp
	  read (itorp,*) nterm(i,j),nlor(i,j)
          v0ij=0.0d0
          si=-1.0d0
	  do k=1,nterm(i,j)
	    read (itorp,*) kk,v1(k,i,j),v2(k,i,j) 
            v0ij=v0ij+si*v1(k,i,j)
            si=-si
          enddo
	  do k=1,nlor(i,j)
	    read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) 
            v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
          enddo
          v0(i,j)=v0ij
        enddo
      enddo
      close (itorp)
      if (lprint) then
	write (iout,'(/a/)') 'Torsional constants:'
	do i=1,ntortyp
	  do j=1,ntortyp
            write (iout,*) 'ityp',i,' jtyp',j
            write (iout,*) 'Fourier constants'
            do k=1,nterm(i,j)
	      write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
            enddo
            write (iout,*) 'Lorenz constants'
            do k=1,nlor(i,j)
	      write (iout,'(3(1pe15.5))') 
     &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
            enddo
          enddo
        enddo
      endif
C
C 6/23/01 Read parameters for double torsionals
C
      do i=1,ntortyp
        do j=1,ntortyp
          do k=1,ntortyp
            read (itordp,'(3a1)') t1,t2,t3
            if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
     &        .or. t3.ne.onelett(k)) then
              write (iout,*) "Error in double torsional parameter file",
     &         i,j,k,t1,t2,t3
               stop "Error in double torsional parameter file"
            endif
            read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
            read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
            read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
            read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
            read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
            read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
     &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
          enddo
        enddo
      enddo
      if (lprint) then
      write (iout,*) 
      write (iout,*) 'Constants for double torsionals'
      do i=1,ntortyp
        do j=1,ntortyp 
          do k=1,ntortyp
            write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
     &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
            write (iout,*)
            write (iout,*) 'Single angles:'
            do l=1,ntermd_1(i,j,k)
              write (iout,'(i5,2f10.5,5x,2f10.5)') l,
     &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
     &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
            enddo
            write (iout,*)
            write (iout,*) 'Pairs of angles:'
            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
            do l=1,ntermd_2(i,j,k)
              write (iout,'(i5,20f10.5)') 
     &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
            enddo
            write (iout,*)
            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
            do l=1,ntermd_2(i,j,k)
              write (iout,'(i5,20f10.5)') 
     &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
            enddo
            write (iout,*)
          enddo
        enddo
      enddo
      endif
#endif
C Read of Side-chain backbone correlation parameters
C Modified 11 May 2012 by Adasko
CCC
C
      read (isccor,*) nsccortyp
      read (isccor,*) (isccortyp(i),i=1,ntyp)
c      write (iout,*) 'ntortyp',ntortyp
      maxinter=3
cc maxinter is maximum interaction sites
      do l=1,maxinter
      do i=1,nsccortyp
        do j=1,nsccortyp
          read (isccor,*)
     &nterm_sccor(i,j),nlor_sccor(i,j)
          v0ijsccor=0.0d0
          si=-1.0d0

          do k=1,nterm_sccor(i,j)
            read (isccor,*) kk,v1sccor(k,l,i,j)
     &    ,v2sccor(k,l,i,j)
            v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
            si=-si
          enddo
          do k=1,nlor_sccor(i,j)
            read (isccor,*) kk,vlor1sccor(k,i,j),
     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
            v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
     &(1+vlor3sccor(k,i,j)**2)
          enddo
          v0sccor(i,j)=v0ijsccor
        enddo
      enddo
      enddo
      close (isccor)

      if (lprint) then
        write (iout,'(/a/)') 'Torsional constants:'
        do l=1,maxinter
        do i=1,nsccortyp
          do j=1,nsccortyp
            write (iout,*) 'ityp',i,' jtyp',j
            write (iout,*) 'Fourier constants'
            do k=1,nterm_sccor(i,j)
              write (iout,'(2(1pe15.5))')
     & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
            enddo
            write (iout,*) 'Lorenz constants'
            do k=1,nlor_sccor(i,j)
              write (iout,'(3(1pe15.5))')
     &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
            enddo
          enddo
        enddo
        enddo
      endif

C
C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
C         interaction energy of the Gly, Ala, and Pro prototypes.
C
      read (ifourier,*) nloctyp
      do i=1,nloctyp
        read (ifourier,*)
        read (ifourier,*) (b(ii,i),ii=1,13)
        if (lprint) then
        write (iout,*) 'Type',i
        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
        endif
        B1(1,i)  = b(3,i)
        B1(2,i)  = b(5,i)
        B1tilde(1,i) = b(3,i)
        B1tilde(2,i) =-b(5,i) 
        B2(1,i)  = b(2,i)
        B2(2,i)  = b(4,i)
        CC(1,1,i)= b(7,i)
        CC(2,2,i)=-b(7,i)
        CC(2,1,i)= b(9,i)
        CC(1,2,i)= b(9,i)
        Ctilde(1,1,i)=b(7,i)
        Ctilde(1,2,i)=b(9,i)
        Ctilde(2,1,i)=-b(9,i)
        Ctilde(2,2,i)=b(7,i)
        DD(1,1,i)= b(6,i)
        DD(2,2,i)=-b(6,i)
        DD(2,1,i)= b(8,i)
        DD(1,2,i)= b(8,i)
        Dtilde(1,1,i)=b(6,i)
        Dtilde(1,2,i)=b(8,i)
        Dtilde(2,1,i)=-b(8,i)
        Dtilde(2,2,i)=b(6,i)
        EE(1,1,i)= b(10,i)+b(11,i)
        EE(2,2,i)=-b(10,i)+b(11,i)
        EE(2,1,i)= b(12,i)-b(13,i)
        EE(1,2,i)= b(12,i)+b(13,i)
      enddo
      if (lprint) then
      do i=1,nloctyp
        write (iout,*) 'Type',i
        write (iout,*) 'B1'
c        write (iout,'(f10.5)') B1(:,i)
        write(iout,*) B1(1,i),B1(2,i)
        write (iout,*) 'B2'
c        write (iout,'(f10.5)') B2(:,i)
        write(iout,*) B2(1,i),B2(2,i)
        write (iout,*) 'CC'
        do j=1,2
          write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
        enddo
        write(iout,*) 'DD'
        do j=1,2
          write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
        enddo
        write(iout,*) 'EE'
        do j=1,2
          write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
        enddo
      enddo
      endif
C 
C Read electrostatic-interaction parameters
C
      if (lprint) then
	write (iout,'(/a)') 'Electrostatic interaction constants:'
	write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
     &            'IT','JT','APP','BPP','AEL6','AEL3'
      endif
      read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
      read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
      read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
      read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
      close (ielep)
      do i=1,2
        do j=1,2
        rri=rpp(i,j)**6
        app (i,j)=epp(i,j)*rri*rri 
        bpp (i,j)=-2.0D0*epp(i,j)*rri
        ael6(i,j)=elpp6(i,j)*4.2D0**6
        ael3(i,j)=elpp3(i,j)*4.2D0**3
        if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
     &                    ael6(i,j),ael3(i,j)
        enddo
      enddo
C
C Read side-chain interaction parameters.
C
      read (isidep,*) ipot,expon
      if (ipot.lt.1 .or. ipot.gt.5) then
        write (iout,'(2a)') 'Error while reading SC interaction',
     &               'potential file - unknown potential type.'
        stop
      endif
      expon2=expon/2
      write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
     & ', exponents are ',expon,2*expon 
      goto (10,20,30,30,40) ipot
C----------------------- LJ potential ---------------------------------
   10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
      if (lprint) then
	write (iout,'(/a/)') 'Parameters of the LJ potential:'
	write (iout,'(a/)') 'The epsilon array:'
	call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
	write (iout,'(/a)') 'One-body parameters:'
	write (iout,'(a,4x,a)') 'residue','sigma'
	write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
      endif
      goto 50
C----------------------- LJK potential --------------------------------
   20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
     &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
      if (lprint) then
	write (iout,'(/a/)') 'Parameters of the LJK potential:'
	write (iout,'(a/)') 'The epsilon array:'
	call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
	write (iout,'(/a)') 'One-body parameters:'
	write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
        write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
     &        i=1,ntyp)
      endif
      goto 50
C---------------------- GB or BP potential -----------------------------
   30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
     &  (alp(i),i=1,ntyp)
C For the GB potential convert sigma'**2 into chi'
      if (ipot.eq.4) then
	do i=1,ntyp
	  chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
        enddo
      endif
      if (lprint) then
	write (iout,'(/a/)') 'Parameters of the BP potential:'
	write (iout,'(a/)') 'The epsilon array:'
	call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
	write (iout,'(/a)') 'One-body parameters:'
	write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
     &       '    chip  ','    alph  '
	write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
     &                     chip(i),alp(i),i=1,ntyp)
      endif
      goto 50
C--------------------- GBV potential -----------------------------------
   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
     &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
     &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
      if (lprint) then
	write (iout,'(/a/)') 'Parameters of the GBV potential:'
	write (iout,'(a/)') 'The epsilon array:'
	call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
	write (iout,'(/a)') 'One-body parameters:'
	write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
     &      's||/s_|_^2','    chip  ','    alph  '
	write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
     &           sigii(i),chip(i),alp(i),i=1,ntyp)
      endif
   50 continue
      close (isidep)
C-----------------------------------------------------------------------
C Calculate the "working" parameters of SC interactions.
      do i=2,ntyp
        do j=1,i-1
	  eps(i,j)=eps(j,i)
        enddo
      enddo
      do i=1,ntyp
        do j=i,ntyp
          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
          sigma(j,i)=sigma(i,j)
          rs0(i,j)=dwa16*sigma(i,j)
          rs0(j,i)=rs0(i,j)
        enddo
      enddo
      if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
     & 'Working parameters of the SC interactions:',
     & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
     & '  chi1   ','   chi2   ' 
      do i=1,ntyp
	do j=i,ntyp
	  epsij=eps(i,j)
	  if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
	    rrij=sigma(i,j)
          else
	    rrij=rr0(i)+rr0(j)
          endif
	  r0(i,j)=rrij
	  r0(j,i)=rrij
	  rrij=rrij**expon
	  epsij=eps(i,j)
	  sigeps=dsign(1.0D0,epsij)
	  epsij=dabs(epsij)
	  aa(i,j)=epsij*rrij*rrij
	  bb(i,j)=-sigeps*epsij*rrij
	  aa(j,i)=aa(i,j)
	  bb(j,i)=bb(i,j)
	  if (ipot.gt.2) then
	    sigt1sq=sigma0(i)**2
	    sigt2sq=sigma0(j)**2
	    sigii1=sigii(i)
	    sigii2=sigii(j)
            ratsig1=sigt2sq/sigt1sq
	    ratsig2=1.0D0/ratsig1
	    chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
	    if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
            rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
          else
	    rsum_max=sigma(i,j)
          endif
c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
            sigmaii(i,j)=rsum_max
            sigmaii(j,i)=rsum_max 
c         else
c           sigmaii(i,j)=r0(i,j)
c           sigmaii(j,i)=r0(i,j)
c         endif
cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
          if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
            r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
            augm(i,j)=epsij*r_augm**(2*expon)
c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
	    augm(j,i)=augm(i,j)
          else
	    augm(i,j)=0.0D0
	    augm(j,i)=0.0D0
          endif
	  if (lprint) then
            write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
     &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
     &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
	  endif
        enddo
      enddo
C
C Define the SC-p interaction constants
C
      do i=1,20
        do j=1,2
          eps_scp(i,j)=-1.5d0
          rscp(i,j)=4.0d0
        enddo
      enddo
#ifdef OLDSCP
      do i=1,20
C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
C helix formation)
c       aad(i,1)=0.3D0*4.0D0**12
C Following line for constants currently implemented
C "Hard" SC-p repulsion (gives correct turn spacing in helices)
        aad(i,1)=1.5D0*4.0D0**12
c       aad(i,1)=0.17D0*5.6D0**12
        aad(i,2)=aad(i,1)
C "Soft" SC-p repulsion
        bad(i,1)=0.0D0
C Following line for constants currently implemented
c       aad(i,1)=0.3D0*4.0D0**6
C "Hard" SC-p repulsion
        bad(i,1)=3.0D0*4.0D0**6
c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
        bad(i,2)=bad(i,1)
c       aad(i,1)=0.0D0
c       aad(i,2)=0.0D0
c       bad(i,1)=1228.8D0
c       bad(i,2)=1228.8D0
      enddo
#else
C
C 8/9/01 Read the SC-p interaction constants from file
C
      do i=1,ntyp
        read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
      enddo
      do i=1,ntyp
        aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
        aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
        bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
        bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
      enddo

      if (lprint) then
        write (iout,*) "Parameters of SC-p interactions:"
        do i=1,20
          write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
     &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
        enddo
      endif
#endif
C
C Define the constants of the disulfide bridge
C
C      ebr=-5.50D0
c
c Old arbitrary potential - commented out.
c
c      dbr= 4.20D0
c      fbr= 3.30D0
c
c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
c energy surface of diethyl disulfide.
c A. Liwo and U. Kozlowska, 11/24/03
c
C      D0CM = 3.78d0
C      AKCM = 15.1d0
C      AKTH = 11.0d0
C      AKCT = 12.0d0
C      V1SS =-1.08d0
C      V2SS = 7.61d0
C      V3SS = 13.7d0
      return
      end