implicit none
      include "DIMENSIONS"
      include "COMMON.CALC"
      integer i,j,k,l,ii,nf,n,uiparm(1)
      double precision x(maxene)
      double precision rmsave(maxT),fdum,rr
      external fdum,funclik
      double precision quot,quotl,f,T0 /300.0d0/
      character*8 ename(maxene)
      common /names/ ename
      character*1 restyp(nbase)
      data restyp /'A','G','C','T','U'/
      character maskchar(0:1) /' ','*'/
      character*80 Naresfile,Fracfile
      double precision g(100)
c      print *,"start"
      read(1,*) nene,sigma2,wsq
      write (2,*) "nene",nene," nT",nT," sigma",sigma2
      read(1,*) (ename(i),i=1,nene)
      read(1,*) (weight(i),i=1,nene)
      read(1,*) (iweight(i),i=1,nene)
      read(1,*) (mask(i),i=1,nene)
      do i=1,nnbase
        read(1,*) sig0(i),(weightel(3*(i-1)+j),maskel(3*(i-1)+j),j=1,3)
      enddo
      read(1,'(a)') Naresfile
      read(1,'(a)') Fracfile
      close(1)
      open(7,file=Naresfile,status='old')
      i=0
      do 
c        print *,"i=",i
        read(7,*,end=10,err=10) ii,(enetb(j,i+1),j=1,nene),
     &             ((eletb(3*(j-1)+k,i+1),j=1,nnbase),k=1,3),
     &                          entfac(i+1),
     &                          qtab(i+1),rmstab(i+1),rgytab(i+1),
     &                          fpair(i+1),rr,fdup(i+1)
        i=i+1
      enddo
   10 continue
      nconf=i
      close(7)
      write (2,*) "nconf",nconf
      write (2,'(/"Initial energy-term weights (* optimized)")')
      write (2,'(i5,2x,a4,f10.5,1x,a1,i5)') 
     &   (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
      ii=0
      write (2,*) 
     & "Initial base-dipole-interaction parameters (* optimizable)"
      do i=1,nbase
        do j=1,i
          ii=ii+1
          write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
     &     sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
     &     k=1,3)
        enddo
      enddo
      write (2,*) "sigma",sigma2," wsq",wsq
      sigma2=sigma2*sigma2
      close(7)
      open(7,file=Fracfile,status='old')
      i=0
      do 
        read(7,*,end=11,err=11) temper(i+1),frac(i+1)
        i=i+1
      enddo
   11 continue
      nT=i
      close(7)
      write (2,*) "Fractions of base pairs"
      do i=1,nT
        write (2,'(i5,f8.1,f10.5)') i,temper(i),frac(i)
      enddo
c Transfer weights to variables
      call w2x(n,x)
      write (2,*) "n",n
c Compute the temperature scale factors
      do i=1,nT
        ft(1,i)=1.0d0
        quot=temper(i)/T0
        quotl=1.0d0
        do l=2,2
          quotl=quotl*quot
          fT(l,i)=1.12692801104297249644d0/
     &             dlog(dexp(quotl)+dexp(-quotl))
        enddo
      enddo
c      do i=1,nT
c        write (2,'(i5,f8.3,f10.5)') (i,(ft(j,i),j=1,2),i=1,nT)
c      enddo
      
      call funclik(n,x,nf,f,uiparm,rmsave,fdum)
      write (2,*) "f",f
      write(2,'(/a3,a8,3a10)')" Nr","    Temp","    rmsave",
     &   "      fave"," fave(exp)"
      do i=1,nT
        write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
     &   fave(i),frac(i)
      enddo
c      call grad(n,x,nf,g,uiparm,rmsave,fdum)
c      write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
c      stop
      call minsumsl(n,x,f)
      call funclik(n,x,nf,f,uiparm,rmsave,fdum)
      write (2,'(/"Final parameters (",i5,")")') n
      do i=1,n
        write (2,'(i5,f10.5)') i,x(i)
      enddo
      write (2,'(/"Energy-term weights (* optimized)")')
      write (2,'(i5,2x,a4,f10.5,1x,a1,i5)') 
     &   (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
      write (2,*) "Base-dipole-interaction parameters (* optimized)"
      ii=0
      do i=1,nbase
        do j=1,i
          ii=ii+1
          write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
     &     sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
     &     k=1,3)
        enddo
      enddo
      write (2,*) "f",f
      write(2,'(/a3,a8,3a10)')" Nr","    Temp","    rmsave",
     &   "      fave"," fave(exp)"
      do i=1,nT
        write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
     &   fave(i),frac(i)
      enddo

c      do i=10,30
c      do k=10,30
c      weight(6)=0.1d0*i
c      weight(1)=0.1d0*k
c       write (2,'(i5,2x,a4,f10.5,i5)') 
c     &   (j,ename(j),weight(j),iweight(j),j=1,nene)
c      call funclik(nene,weight,nf,f,uiparm,rmsave,fdum)
c      write (2,*) "f",f
c      enddo
c      enddo
      stop
      end
c-------------------------------------------------------------------------------
      subroutine funclik(n,x,nf,f,uiparm,rmsave,ufparm) 
      implicit none
      include "DIMENSIONS"
      include "COMMON.CALC"
      character*8 ename(maxene)
      common /names/ ename
      integer n,nf
      double precision f,loglik,chisq
      double precision x(n)
      integer uiparm
      double precision ufparm
      external ufparm
      double precision ww(maxene),sumlik(maxT),rmsave(maxT)
      integer it,i,j,k
      double precision beta,over,boltz,sumQ,ener(maxconf),emin,ee,enel,
     &   sumover,eboltz
      call x2w(n,x)
      loglik=0.0d0
      chisq=0.0d0
      do iT=1,nT
        do i=1,nene
          ww(i)=weight(i)*ft(iweight(i),iT)
        enddo
c        write (2,*) "iT",iT," temper",temper(iT)," beta",beta
c        write (2,'(i5,2x,a4,2f10.5,i5,f10.5)') 
c     &   (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
c     &    i=1,n)
        beta=1.0d0/(temper(iT)*1.987D-3)
c        write (2,*) "beta",beta
        emin=1.0d10
        do i=1,nconf
          enel=0.0d0
          do j=1,nnbase
c            write (2,*) 
c     &       i,j,sig0(j),(weightel(3*(j-1)+k),eletb(3*(j-1)+k,i),k=1,3)
            enel=enel+weightel(3*(j-1)+1)*sig0(j)**6*eletb(3*(j-1)+1,i)
     &               +weightel(3*(j-1)+2)*sig0(j)**3*eletb(3*(j-1)+2,i)
     &               +weightel(3*j)*sig0(j)**6*eletb(3*j,i)
          enddo
c          write (2,*) i,enel,enetb(6,i)
          enetb(6,i)=enel
          ener(i)=0.0d0
          do j=1,nene
            ener(i)=ener(i)+ww(j)*enetb(j,i)
          enddo
          ee = ener(i)-entfac(i)/beta
          if (ee.lt.emin) emin=ee
        enddo
        rmsave(it)=0.0d0
        sumlik(it)=0.0d0
        fave(it)=0.0d0
        sumQ=0.0d0
        sumover=0.0d0
        do i=1,nconf
          boltz=-beta*(ener(i)-emin)+entfac(i)
          eboltz=dexp(boltz)
          rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
          over=dexp(-0.5d0*rmstab(i)**2/sigma2)
          if (frac(iT).lt.0.5d0) over=1.0d0-over
          sumover=sumover+over
c          write (2,*) i,ener(i),entfac(i),rmstab(i),fpair(i),fdup(i),
c     &      over,boltz,eboltz
          fave(iT)=fave(iT)+fdup(i)*eboltz
          sumQ=sumQ+eboltz
          sumlik(iT)=sumlik(iT)+over*boltz 
        enddo
        fave(iT)=fave(iT)/sumq
        sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
        rmsave(iT)=rmsave(iT)/sumQ
        if (frac(iT).gt.0.95d0 .or. frac(iT).lt.0.05d0) then
          loglik=loglik-sumlik(iT)
        endif
        chisq=chisq+(frac(iT)-fave(iT))**2
c        write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover,
c     &   fave(iT),frac(iT)
c        call flush(2)
      enddo
      f = loglik/nconf+wsq*chisq 
c      write (2,*) "loglik",loglik/nconf," chisq",chisq," f",f
      return
      end
c-------------------------------------------------------------------------------
      subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
      implicit none
      integer n,nf,uiparm(1)
      double precision x(n),g(n),urparm(1),ufparm
      external ufparm
      integer i
      double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
      do i=1,n
        xi=x(i) 
        x(i)=xi+delta
        call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
        x(i)=xi-delta
        call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
        g(i)=(fplus-fminus)/delta2 
c        write(2,*) i,fplus,fminus,g(i)
      enddo
      return
      end
c-------------------------------------------------------------------------------
      double precision function fdum()
      fdum=0.0d0
      return
      end
c-------------------------------------------------------------------------------
      logical function stopx(idum)
      integer idum
      stopx=.false.
      return
      end
c-------------------------------------------------------------------------------
      subroutine x2w(n,x)
      include "DIMENSIONS"
      include "COMMON.CALC"
      double precision fabs
      integer n,i,ii
      double precision x(n)
      ii=0
      do i=1,nene
        if (mask(i).gt.0) then
          ii=ii+1
          weight(i)=fabs(x(ii))
        endif
      enddo
      do i=1,nnbase
        if (maskel(3*(i-1)+1).gt.0) then
          ii=ii+1
          weightel(3*(i-1)+1)=-fabs(x(ii))
        endif
        if (maskel(3*(i-1)+2).gt.0) then
          ii=ii+1
          weightel(3*(i-1)+2)=x(ii)
        endif
        if (maskel(3*i).gt.0) then
          ii=ii+1
          weightel(3*i)=-fabs(x(ii))
        endif
      enddo
      return
      end
c-------------------------------------------------------------------------------
      subroutine w2x(n,x)
      include "DIMENSIONS"
      include "COMMON.CALC"
      integer n,i,ii
      double precision x(maxene+3*nnbase)
      ii=0
      do i=1,nene
c        write (2,*) "i",i," mask",mask(i)," ii",ii
        if (mask(i).gt.0) then
          ii=ii+1
          x(ii)=weight(i)
        endif
      enddo
      do i=1,3*nnbase
c        write (2,*) "i",i," maskel",maskel(i)," ii",ii
        if (maskel(i).gt.0) then
          ii=ii+1
          x(ii)=weightel(i)
        endif
      enddo
      n=ii
c      write (2,*) "W2X: n=",n
      return
      end
c------------------------------------------------------------------------------------
      double precision function fabs(x)
      double precision x
      double precision a /100.0d0/
      if (dabs(x).gt.1.0d-2) then
        fabs = dabs(x)
      else
        fabs = dlog(dexp(a*x)+dexp(-a*x))/a
      endif
      return
      end