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
8 double precision quot,quotl,f,T0 /300.0d0/
9 character*8 ename(maxene)
12 read(1,*) nene,sigma2,temper(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)
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)
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)
36 c Transfer weights to variables
38 call funclik(n,x,nf,f,uiparm,rmsave,fdum)
40 write (2,'(i5,2e15.5)') i,ener0(i),ener(i)
42 write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",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)
48 call funclik(n,x,nf,f,uiparm,rmsave,fdum)
49 write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
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)
63 c-------------------------------------------------------------------------------
64 subroutine funclik(n,x,nf,f,uiparm,rmsave,ufparm)
67 include "COMMON.CALC-single"
68 character*8 ename(maxene)
74 double precision ufparm
76 double precision ww(maxene),sumlik(maxT),rmsave(maxT)
78 double precision beta,over,boltz,sumQ,emin,ee,
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),
93 c write (2,*) "beta",beta
98 ener(i)=ener(i)+ww(j)*enetb(j,i)
100 ee = ener(i)-entfac(i)/beta
101 if (ee.lt.emin) emin=ee
108 over=dexp(-0.5d0*rmstab(i)**2/sigma2)
109 c if (temper(iT).gt.340.0d0) over=1.0d0-over
111 boltz=-beta*(ener(i)-emin)+entfac(i)
112 c write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
114 sumQ=sumQ+dexp(boltz)
115 rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
116 sumlik(iT)=sumlik(iT)+over*boltz
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
126 c-------------------------------------------------------------------------------
127 subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
129 integer n,nf,uiparm(1)
130 double precision x(n),g(n),urparm(1),ufparm
133 double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
137 call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
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)
145 c-------------------------------------------------------------------------------
146 double precision function fdum()
150 c-------------------------------------------------------------------------------
151 logical function stopx(idum)
156 c-------------------------------------------------------------------------------
159 include "COMMON.CALC-single"
161 double precision x(n)
162 double precision fabs
165 if (mask(i).gt.0) then
167 weight(i)=fabs(x(ii))
172 c-------------------------------------------------------------------------------
175 include "COMMON.CALC-single"
177 double precision x(maxene)
180 if (mask(i).gt.0) then
188 c-------------------------------------------------------------------------------
189 double precision function fabs(x)
191 double precision a /1.0d4/
192 if (dabs(x).gt.1.0d-2) then
195 fabs = dlog(dexp(a*x)+dexp(-a*x))/a