--- /dev/null
+ 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,rjunk
+ external fdum,funclik
+ double precision quot,quotl,f,T0 /300.0d0/
+ character*8 ename(maxene)
+ common /names/ ename
+c print *,"start"
+ read(1,*) nene,sigma2,temper(1)
+ nT = 1
+ write (2,*) "nene",nene," nT",nT," sigma",sigma2
+ read(1,*) (ename(i),i=1,nene)
+ read(1,*) (weight(i),i=1,nene)
+c read(1,*) (iweight(i),i=1,nene)
+ read(1,*) (mask(i),i=1,nene)
+c read(1,*) (temper(i),i=1,nT)
+ i=0
+ do
+c print *,"i=",i
+ read(1,*,end=10,err=10) ii,(enetb(j,i+1),j=1,nene),ener0(i+1),
+ & rmstab(i+1),rjunk,rjunk,rjunk,qtab(i+1)
+ i=i+1
+ enddo
+ 10 continue
+ nconf=i
+ do i=1,nconf
+ entfac(i)=0.0d0
+ enddo
+ write (2,*) "nconf",nconf
+ write (2,'(i5,2x,a4,f10.5,2i5)')
+ & (i,ename(i),weight(i),iweight(i),mask(i),i=1,nene)
+ sigma2=sigma2*sigma2
+c Transfer weights to variables
+ call w2x(n,x)
+ call funclik(n,x,nf,f,uiparm,rmsave,fdum)
+ do i=1,nconf
+ write (2,'(i5,2e15.5)') i,ener0(i),ener(i)
+ enddo
+ write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
+ call minsumsl(n,x,f)
+ write (2,*) "n",n," x",(x(i),i=1,n)
+ write (2,'(i5,2x,a4,f10.5,i5)')
+ & (j,ename(j),weight(j),iweight(j),j=1,nene)
+ write (2,*) "f",f
+ call funclik(n,x,nf,f,uiparm,rmsave,fdum)
+ write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
+
+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
+ double precision x(n)
+ integer uiparm
+ double precision ufparm
+ external ufparm
+ double precision ww(maxene),sumlik(maxT),rmsave(maxT)
+ integer it,i,j
+ double precision beta,over,boltz,sumQ,emin,ee,
+ & sumover
+ call x2w(n,x)
+ f=0.0d0
+ do iT=1,nT
+ do i=1,nene
+ ww(i)=weight(i)
+ enddo
+ write (2,*) "iT",iT," temper",temper(iT)," beta",beta
+ write (2,'(i5,2x,a4,2f10.5,i5)')
+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),
+ & (i,ename(i),weight(i),ww(i),iweight(i),
+ & i=21,20+n)
+ beta=1.0d0/(temper(iT)*1.987D-3)
+c write (2,*) "beta",beta
+ emin=1.0d10
+ do i=1,nconf
+ 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
+ sumQ=0.0d0
+ sumover=0.0d0
+ do i=1,nconf
+crms over=dexp(-0.5d0*rmstab(i)**2/sigma2)
+ over=dexp(-0.5d0*(1-qtab(i))**2/sigma2)
+c if (temper(iT).gt.340.0d0) over=1.0d0-over
+ sumover=sumover+over
+ boltz=-beta*(ener(i)-emin)+entfac(i)
+c write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
+c & dexp(boltz)
+ sumQ=sumQ+dexp(boltz)
+crms rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
+ rmsave(iT)=rmsave(iT)+(1-qtab(i))*dexp(boltz)
+ sumlik(iT)=sumlik(iT)+over*boltz
+ enddo
+ sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
+ rmsave(iT)=rmsave(iT)/sumQ
+ write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
+c write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover
+ f=f-sumlik(iT)
+ enddo
+ 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"
+ integer n,i,ii
+ double precision x(n)
+ double precision fabs
+ ii=0
+ do i=1,nene
+ if (mask(i).gt.0) then
+ ii=ii+1
+ weight(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)
+ ii=0
+ do i=1,nene
+ if (mask(i).gt.0) then
+ ii=ii+1
+ x(ii)=weight(i)
+ endif
+ enddo
+ n=ii
+ return
+ end
+c-------------------------------------------------------------------------------
+ double precision function fabs(x)
+ double precision x
+ double precision a /100.0d0/
+ if (dabs(x).gt.1.0d-4) then
+ fabs = dabs(x)
+ else
+ fabs = dlog(dexp(a*x)+dexp(-a*x))/a
+ endif
+ return
+ end
projects / unres.git / blobdiff