4 integer i,j,k,l,ii,nf,n,uiparm(1)
5 double precision x(maxene)
6 double precision rmsave(maxT),fdum,rr
8 double precision quot,quotl,f,T0 /300.0d0/
9 character*8 ename(maxene)
11 character*1 restyp(nbase)
12 data restyp /'A','G','C','T','U'/
13 character maskchar(0:1) /' ','*'/
14 character*80 Naresfile,Fracfile
15 double precision g(100)
17 read(1,*) nene,sigma2,wsq
18 write (2,*) "nene",nene," nT",nT," sigma",sigma2
19 read(1,*) (ename(i),i=1,nene)
20 read(1,*) (weight(i),i=1,nene)
21 read(1,*) (iweight(i),i=1,nene)
22 read(1,*) (mask(i),i=1,nene)
24 read(1,*) sig0(i),(weightel(3*(i-1)+j),maskel(3*(i-1)+j),j=1,3)
26 read(1,'(a)') Naresfile
27 read(1,'(a)') Fracfile
29 open(7,file=Naresfile,status='old')
33 read(7,*,end=10,err=10) ii,(enetb(j,i+1),j=1,nene),
34 & ((eletb(3*(j-1)+k,i+1),j=1,nnbase),k=1,3),
36 & qtab(i+1),rmstab(i+1),rgytab(i+1),
37 & fpair(i+1),rr,fdup(i+1)
43 write (2,*) "nconf",nconf
44 write (2,'(/"Initial energy-term weights (* optimized)")')
45 write (2,'(i5,2x,a4,f10.5,1x,a1,i5)')
46 & (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
49 & "Initial base-dipole-interaction parameters (* optimizable)"
53 write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
54 & sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
58 write (2,*) "sigma",sigma2," wsq",wsq
61 open(7,file=Fracfile,status='old')
64 read(7,*,end=11,err=11) temper(i+1),frac(i+1)
70 write (2,*) "Fractions of base pairs"
72 write (2,'(i5,f8.1,f10.5)') i,temper(i),frac(i)
74 c Transfer weights to variables
77 c Compute the temperature scale factors
84 fT(l,i)=1.12692801104297249644d0/
85 & dlog(dexp(quotl)+dexp(-quotl))
89 c write (2,'(i5,f8.3,f10.5)') (i,(ft(j,i),j=1,2),i=1,nT)
92 call funclik(n,x,nf,f,uiparm,rmsave,fdum)
94 write(2,'(/a3,a8,3a10)')" Nr"," Temp"," rmsave",
95 & " fave"," fave(exp)"
97 write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
100 c call grad(n,x,nf,g,uiparm,rmsave,fdum)
101 c write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
104 call funclik(n,x,nf,f,uiparm,rmsave,fdum)
105 write (2,'(/"Final parameters (",i5,")")') n
107 write (2,'(i5,f10.5)') i,x(i)
109 write (2,'(/"Energy-term weights (* optimized)")')
110 write (2,'(i5,2x,a4,f10.5,1x,a1,i5)')
111 & (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
112 write (2,*) "Base-dipole-interaction parameters (* optimized)"
117 write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
118 & sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
123 write(2,'(/a3,a8,3a10)')" Nr"," Temp"," rmsave",
124 & " fave"," fave(exp)"
126 write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
134 c write (2,'(i5,2x,a4,f10.5,i5)')
135 c & (j,ename(j),weight(j),iweight(j),j=1,nene)
136 c call funclik(nene,weight,nf,f,uiparm,rmsave,fdum)
142 c-------------------------------------------------------------------------------
143 subroutine funclik(n,x,nf,f,uiparm,rmsave,ufparm)
146 include "COMMON.CALC"
147 character*8 ename(maxene)
150 double precision f,loglik,chisq
151 double precision x(n)
153 double precision ufparm
155 double precision ww(maxene),sumlik(maxT),rmsave(maxT)
157 double precision beta,over,boltz,sumQ,ener(maxconf),emin,ee,enel,
164 ww(i)=weight(i)*ft(iweight(i),iT)
166 c write (2,*) "iT",iT," temper",temper(iT)," beta",beta
167 c write (2,'(i5,2x,a4,2f10.5,i5,f10.5)')
168 c & (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
170 beta=1.0d0/(temper(iT)*1.987D-3)
171 c write (2,*) "beta",beta
177 c & i,j,sig0(j),(weightel(3*(j-1)+k),eletb(3*(j-1)+k,i),k=1,3)
178 enel=enel+weightel(3*(j-1)+1)*sig0(j)**6*eletb(3*(j-1)+1,i)
179 & +weightel(3*(j-1)+2)*sig0(j)**3*eletb(3*(j-1)+2,i)
180 & +weightel(3*j)*sig0(j)**6*eletb(3*j,i)
182 c write (2,*) i,enel,enetb(6,i)
186 ener(i)=ener(i)+ww(j)*enetb(j,i)
188 ee = ener(i)-entfac(i)/beta
189 if (ee.lt.emin) emin=ee
197 boltz=-beta*(ener(i)-emin)+entfac(i)
199 rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
200 over=dexp(-0.5d0*rmstab(i)**2/sigma2)
201 if (frac(iT).lt.0.5d0) over=1.0d0-over
203 c write (2,*) i,ener(i),entfac(i),rmstab(i),fpair(i),fdup(i),
204 c & over,boltz,eboltz
205 fave(iT)=fave(iT)+fdup(i)*eboltz
207 sumlik(iT)=sumlik(iT)+over*boltz
209 fave(iT)=fave(iT)/sumq
210 sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
211 rmsave(iT)=rmsave(iT)/sumQ
212 if (frac(iT).gt.0.95d0 .or. frac(iT).lt.0.05d0) then
213 loglik=loglik-sumlik(iT)
215 chisq=chisq+(frac(iT)-fave(iT))**2
216 c write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover,
217 c & fave(iT),frac(iT)
220 f = loglik/nconf+wsq*chisq
221 c write (2,*) "loglik",loglik/nconf," chisq",chisq," f",f
224 c-------------------------------------------------------------------------------
225 subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
227 integer n,nf,uiparm(1)
228 double precision x(n),g(n),urparm(1),ufparm
231 double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
235 call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
237 call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
238 g(i)=(fplus-fminus)/delta2
239 c write(2,*) i,fplus,fminus,g(i)
243 c-------------------------------------------------------------------------------
244 double precision function fdum()
248 c-------------------------------------------------------------------------------
249 logical function stopx(idum)
254 c-------------------------------------------------------------------------------
257 include "COMMON.CALC"
258 double precision fabs
260 double precision x(n)
263 if (mask(i).gt.0) then
265 weight(i)=fabs(x(ii))
269 if (maskel(3*(i-1)+1).gt.0) then
271 weightel(3*(i-1)+1)=-fabs(x(ii))
273 if (maskel(3*(i-1)+2).gt.0) then
275 weightel(3*(i-1)+2)=x(ii)
277 if (maskel(3*i).gt.0) then
279 weightel(3*i)=-fabs(x(ii))
284 c-------------------------------------------------------------------------------
287 include "COMMON.CALC"
289 double precision x(maxene+3*nnbase)
292 c write (2,*) "i",i," mask",mask(i)," ii",ii
293 if (mask(i).gt.0) then
299 c write (2,*) "i",i," maskel",maskel(i)," ii",ii
300 if (maskel(i).gt.0) then
306 c write (2,*) "W2X: n=",n
309 c------------------------------------------------------------------------------------
310 double precision function fabs(x)
312 double precision a /100.0d0/
313 if (dabs(x).gt.1.0d-2) then
316 fabs = dlog(dexp(a*x)+dexp(-a*x))/a