source/maxlik/src_CSA/maxlik-opt-tmscore.f

   1       implicit none
   2       include "DIMENSIONS"
   3       include "COMMON.CALC"
   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,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"
  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         write (2,*) "iT",iT," temper",temper(iT)," beta",beta
  87         write (2,'(i5,2x,a4,2f10.5,i5)')
  88 c        write (2,'(i5,2x,a4,2f10.5,i5,f10.5)')
  89 c     &   (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
  90      &   (i,ename(i),weight(i),ww(i),iweight(i),
  91      &    i=21,20+n)
  92         beta=1.0d0/(temper(iT)*1.987D-3)
  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 crms          over=dexp(-0.5d0*rmstab(i)**2/sigma2)
 109           over=dexp(-0.5d0*(1-qtab(i))**2/sigma2)
 110 c          if (temper(iT).gt.340.0d0) over=1.0d0-over
 111           sumover=sumover+over
 112           boltz=-beta*(ener(i)-emin)+entfac(i)
 113 c          write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
 114 c     &      dexp(boltz)
 115           sumQ=sumQ+dexp(boltz)
 116 crms          rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
 117           rmsave(iT)=rmsave(iT)+(1-qtab(i))*dexp(boltz)
 118           sumlik(iT)=sumlik(iT)+over*boltz
 119         enddo
 120         sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
 121         rmsave(iT)=rmsave(iT)/sumQ
 122         write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
 123 c        write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
 124         f=f-sumlik(iT)
 125       enddo
 126       return
 127       end
 128 c-------------------------------------------------------------------------------
 129       subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
 130       implicit none
 131       integer n,nf,uiparm(1)
 132       double precision x(n),g(n),urparm(1),ufparm
 133       external ufparm
 134       integer i
 135       double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
 136       do i=1,n
 137         xi=x(i)
 138         x(i)=xi+delta
 139         call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
 140         x(i)=xi-delta
 141         call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
 142         g(i)=(fplus-fminus)/delta2
 143 c        write(2,*) i,fplus,fminus,g(i)
 144       enddo
 145       return
 146       end
 147 c-------------------------------------------------------------------------------
 148       double precision function fdum()
 149       fdum=0.0d0
 150       return
 151       end
 152 c-------------------------------------------------------------------------------
 153       logical function stopx(idum)
 154       integer idum
 155       stopx=.false.
 156       return
 157       end
 158 c-------------------------------------------------------------------------------
 159       subroutine x2w(n,x)
 160       include "DIMENSIONS"
 161       include "COMMON.CALC"
 162       integer n,i,ii
 163       double precision x(n)
 164       double precision fabs
 165       ii=0
 166       do i=1,nene
 167         if (mask(i).gt.0) then
 168           ii=ii+1
 169           weight(i)=fabs(x(ii))
 170         endif
 171       enddo
 172       return
 173       end
 174 c-------------------------------------------------------------------------------
 175       subroutine w2x(n,x)
 176       include "DIMENSIONS"
 177       include "COMMON.CALC"
 178       integer n,i,ii
 179       double precision x(maxene)
 180       ii=0
 181       do i=1,nene
 182         if (mask(i).gt.0) then
 183           ii=ii+1
 184           x(ii)=weight(i)
 185         endif
 186       enddo
 187       n=ii
 188       return
 189       end
 190 c-------------------------------------------------------------------------------
 191       double precision function fabs(x)
 192       double precision x
 193       double precision a /100.0d0/
 194       if (dabs(x).gt.1.0d-4) then
 195         fabs = dabs(x)
 196       else
 197         fabs = dlog(dexp(a*x)+dexp(-a*x))/a
 198       endif
 199       return
 200       end