source/maxlik/src_CSA/maxlik-opt.f

   1       implicit none
   2       include "DIMENSIONS"
   3       include "COMMON.CALC-single"
   4       integer i,j,k,l,ii,nf,n,uiparm(1)
   5       double precision x(maxene)
   6       double precision rmsave(maxT),fdum,rjunk
   7       external fdum,funclik
   8       double precision quot,quotl,f,T0 /300.0d0/
   9       character*8 ename(maxene)
  10       common /names/ ename
  11 c      print *,"start"
  12       read(1,*) nene,sigma2,temper(1)
  13       nT = 1
  14       write (2,*) "nene",nene," nT",nT," sigma",sigma2
  15       read(1,*) (ename(i),i=1,nene)
  16       read(1,*) (weight(i),i=1,nene)
  17 c      read(1,*) (iweight(i),i=1,nene)
  18       read(1,*) (mask(i),i=1,nene)
  19 c      read(1,*) (temper(i),i=1,nT)
  20       i=0
  21       do
  22 c        print *,"i=",i
  23         read(1,*,end=10,err=10) ii,(enetb(j,i+1),j=1,nene),ener0(i+1),
  24      &                          rmstab(i+1),rjunk,rjunk,qtab(i+1)
  25         i=i+1
  26       enddo
  27    10 continue
  28        nconf=i
  29        do i=1,nconf
  30          entfac(i)=0.0d0
  31        enddo
  32        write (2,*) "nconf",nconf
  33        write (2,'(i5,2x,a4,f10.5,2i5)')
  34      &   (i,ename(i),weight(i),iweight(i),mask(i),i=1,nene)
  35       sigma2=sigma2*sigma2
  36 c Transfer weights to variables
  37       call w2x(n,x)
  38       call funclik(n,x,nf,f,uiparm,rmsave,fdum)
  39       do i=1,nconf
  40         write (2,'(i5,2e15.5)') i,ener0(i),ener(i)
  41       enddo
  42       write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
  43       call minsumsl(n,x,f)
  44       write (2,*) "n",n," x",(x(i),i=1,n)
  45        write (2,'(i5,2x,a4,f10.5,i5)')
  46      &   (j,ename(j),weight(j),iweight(j),j=1,nene)
  47       write (2,*) "f",f
  48       call funclik(n,x,nf,f,uiparm,rmsave,fdum)
  49       write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
  50
  51 c      do i=10,30
  52 c      do k=10,30
  53 c      weight(6)=0.1d0*i
  54 c      weight(1)=0.1d0*k
  55 c       write (2,'(i5,2x,a4,f10.5,i5)')
  56 c     &   (j,ename(j),weight(j),iweight(j),j=1,nene)
  57 c      call funclik(nene,weight,nf,f,uiparm,rmsave,fdum)
  58 c      write (2,*) "f",f
  59 c      enddo
  60 c      enddo
  61       stop
  62       end
  63 c-------------------------------------------------------------------------------
  64       subroutine funclik(n,x,nf,f,uiparm,rmsave,ufparm)
  65       implicit none
  66       include "DIMENSIONS"
  67       include "COMMON.CALC-single"
  68       character*8 ename(maxene)
  69       common /names/ ename
  70       integer n,nf
  71       double precision f
  72       double precision x(n)
  73       integer uiparm
  74       double precision ufparm
  75       external ufparm
  76       double precision ww(maxene),sumlik(maxT),rmsave(maxT)
  77       integer it,i,j
  78       double precision beta,over,boltz,sumQ,emin,ee,
  79      &   sumover
  80       call x2w(n,x)
  81       f=0.0d0
  82       do iT=1,nT
  83         do i=1,nene
  84           ww(i)=weight(i)
  85         enddo
  86         beta=1.0d0/(temper(iT)*1.987D-3)
  87         write (2,*) "iT",iT," temper",temper(iT)," beta",beta
  88         write (2,'(i5,2x,a4,2f10.5,i5)')
  89 c        write (2,'(i5,2x,a4,2f10.5,i5,f10.5)')
  90 c     &   (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
  91      &   (i,ename(i),weight(i),ww(i),iweight(i),
  92      &    i=21,20+n)
  93 c        write (2,*) "beta",beta
  94         emin=1.0d10
  95         do i=1,nconf
  96           ener(i)=0.0d0
  97           do j=1,nene
  98             ener(i)=ener(i)+ww(j)*enetb(j,i)
  99           enddo
 100           ee = ener(i)-entfac(i)/beta
 101           if (ee.lt.emin) emin=ee
 102         enddo
 103         rmsave(it)=0.0d0
 104         sumlik(it)=0.0d0
 105         sumQ=0.0d0
 106         sumover=0.0d0
 107         do i=1,nconf
 108           over=dexp(-0.5d0*rmstab(i)**2/sigma2)
 109 c          if (temper(iT).gt.340.0d0) over=1.0d0-over
 110           sumover=sumover+over
 111           boltz=-beta*(ener(i)-emin)+entfac(i)
 112 c          write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
 113 c     &      dexp(boltz)
 114           sumQ=sumQ+dexp(boltz)
 115           rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
 116           sumlik(iT)=sumlik(iT)+over*boltz
 117         enddo
 118         sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
 119         rmsave(iT)=rmsave(iT)/sumQ
 120         write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
 121 c        write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
 122         f=f-sumlik(iT)
 123       enddo
 124       return
 125       end
 126 c-------------------------------------------------------------------------------
 127       subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
 128       implicit none
 129       integer n,nf,uiparm(1)
 130       double precision x(n),g(n),urparm(1),ufparm
 131       external ufparm
 132       integer i
 133       double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
 134       do i=1,n
 135         xi=x(i)
 136         x(i)=xi+delta
 137         call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
 138         x(i)=xi-delta
 139         call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
 140         g(i)=(fplus-fminus)/delta2
 141 c        write(2,*) i,fplus,fminus,g(i)
 142       enddo
 143       return
 144       end
 145 c-------------------------------------------------------------------------------
 146       double precision function fdum()
 147       fdum=0.0d0
 148       return
 149       end
 150 c-------------------------------------------------------------------------------
 151       logical function stopx(idum)
 152       integer idum
 153       stopx=.false.
 154       return
 155       end
 156 c-------------------------------------------------------------------------------
 157       subroutine x2w(n,x)
 158       include "DIMENSIONS"
 159       include "COMMON.CALC-single"
 160       integer n,i,ii
 161       double precision x(n)
 162       double precision fabs
 163       ii=0
 164       do i=1,nene
 165         if (mask(i).gt.0) then
 166           ii=ii+1
 167           weight(i)=fabs(x(ii))
 168         endif
 169       enddo
 170       return
 171       end
 172 c-------------------------------------------------------------------------------
 173       subroutine w2x(n,x)
 174       include "DIMENSIONS"
 175       include "COMMON.CALC-single"
 176       integer n,i,ii
 177       double precision x(maxene)
 178       ii=0
 179       do i=1,nene
 180         if (mask(i).gt.0) then
 181           ii=ii+1
 182           x(ii)=weight(i)
 183         endif
 184       enddo
 185       n=ii
 186       return
 187       end
 188 c-------------------------------------------------------------------------------
 189       double precision function fabs(x)
 190       double precision x
 191       double precision a /1.0d4/
 192       if (dabs(x).gt.1.0d-2) then
 193         fabs = dabs(x)
 194       else
 195         fabs = dlog(dexp(a*x)+dexp(-a*x))/a
 196       endif
 197       return
 198       end