X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Fmaxlik%2Fsrc_CSA%2Fmaxlik-opt-el.f;fp=source%2Fmaxlik%2Fsrc_CSA%2Fmaxlik-opt-el.f;h=0000000000000000000000000000000000000000;hp=7ff3c244727e975084887df723138b438d6e1ae4;hb=af72f8e89a5d33f0d86ba898d6c5bbbda4b25b84;hpb=f9e536df1fd1627429123fe8990edfcdc2cc1a9a

diff --git a/source/maxlik/src_CSA/maxlik-opt-el.f b/source/maxlik/src_CSA/maxlik-opt-el.f
deleted file mode 100644
index 7ff3c24..0000000
--- a/source/maxlik/src_CSA/maxlik-opt-el.f
+++ /dev/null
@@ -1,319 +0,0 @@
-      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