4 integer i,j,k,l,ii,iprot,nf,n,uiparm(1),ienecheck
5 double precision x(maxene)
6 double precision urparm(1),fdum,rjunk
8 double precision quot,quotl,f,T0 /300.0d0/
9 character*8 ename(maxene)
12 character*32 protfile(maxprot)
14 read(1,*) nprot,nene,sigma2,temper(1),ienecheck
16 write (2,*) "nene",nene," nT",nT," sigma",sigma2,
17 & " ienechek",ienecheck
18 read(1,*) (ename(i),i=1,nene)
19 read(1,*) (weight(i),i=1,nene)
20 c read(1,*) (iweight(i),i=1,nene)
21 read(1,*) (mask(i),i=1,nene)
22 c read(1,*) (temper(i),i=1,nT)
26 read (1,'(2a)') protname(iprot)
27 read (1,'(2a)') protfile(iprot)
29 write (2,*) "Reading data for protein ",protname(iprot)
30 write (2,*) "File: ",protfile(iprot)
32 open (7,file=protfile(iprot),status='old')
37 read(7,'(a)',end=10) karta
38 if (index(karta,'#').gt.0) cycle
39 read(karta,*,end=10,err=10) ii,(enetb(j,i+1,iprot),j=1,nene),
40 & ener0(i+1,iprot),rmstab(i+1,iprot),rjunk,rjunk,qtab(i+1,iprot)
48 write (2,*) "Protein:",iprot, " nconf",nconf(iprot)
54 write (2,'(i5,2x,a4,f10.5,2i5)')
55 & (i,ename(i),weight(i),iweight(i),mask(i),i=1,nene)
57 c Transfer weights to variables
59 c write (2,*) "BEFORE funclik: x",(x(i),i=1,n)
60 call funclik(n,x,nf,f,uiparm,urparm,fdum)
62 if (ienecheck.gt.0) then
63 write (2,*) "Checking energies"
65 write (2,*) "Protein",iprot," name",protname(iprot)
66 write (2,'(a5,2a15)') "Conf","UNRES-calc E","Initial E"
68 write (2,'(i5,2e15.5)') i,ener0(i,iprot),ener(i,iprot)
74 write (2,*) "Protein",iprot," name",protname(iprot)
75 write (2,*) "Initial average rmsd(s)"
76 write (2,*) "rmsave",(rmsave(i,iprot),i=1,nT)
78 write (2,*) "Initial target function f",f
81 write (2,*) "n",n," x",(x(i),i=1,n)
82 write (2,'(i5,2x,a4,f10.5,i5)')
83 & (j,ename(j),weight(j),iweight(j),j=1,nene)
85 call funclik(n,x,nf,f,uiparm,urparm,fdum)
88 write (2,*) "Protein",iprot," name",protname(iprot)
89 write (2,*) "rmsave",(rmsave(i,iprot),i=1,nT)
91 write (2,*) "Final target function f",f
97 c write (2,'(i5,2x,a4,f10.5,i5)')
98 c & (j,ename(j),weight(j),iweight(j),j=1,nene)
99 c call funclik(nene,weight,nf,f,uiparm,rmsave,fdum)
105 c-------------------------------------------------------------------------------
106 subroutine funclik(n,x,nf,f,uiparm,urparm,ufparm)
109 include "COMMON.CALC"
110 character*8 ename(maxene)
114 double precision x(n)
116 double precision ufparm
118 double precision ww(maxene)
119 double precision urparm(1)
121 double precision beta,over,boltz,sumQ,emin,ee,
123 c write (2,*) "funclik: x",(x(i),i=1,n)
129 c write (2,'(i5,2x,a4,f10.5,i5)')
130 c & (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
131 c & (i,ename(i),weight(i),iweight(i),
141 beta=1.0d0/(temper(iT)*1.987D-3)
142 c write (2,*) "iprot",iprot," iT",iT," temper",temper(iT),
144 c write (2,*) "beta",beta
149 ener(i,iprot)=ener(i,iprot)+ww(j)*enetb(j,i,iprot)
151 ee = ener(i,iprot)-entfac(i,iprot)/beta
152 if (ee.lt.emin) emin=ee
154 rmsave(it,iprot)=0.0d0
155 sumlik(it,iprot)=0.0d0
159 over=dexp(-0.5d0*rmstab(i,iprot)**2/sigma2)
160 c if (temper(iT).gt.340.0d0) over=1.0d0-over
162 boltz=-beta*(ener(i,iprot)-emin)+entfac(i,iprot)
163 c write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
165 sumQ=sumQ+dexp(boltz)
166 rmsave(iT,iprot)=rmsave(iT,iprot)+rmstab(i,iprot)*dexp(boltz)
167 sumlik(iT,iprot)=sumlik(iT,iprot)+over*boltz
169 sumlik(it,iprot)=sumlik(iT,iprot)-dlog(sumQ)*sumover
170 rmsave(iT,iprot)=rmsave(iT,iprot)/sumQ
171 c write (2,*) iprot,iT,temper(iT),rmsave(iT,iprot),
172 c & sumlik(iT,iprot),sumQ,sumover
173 c write (2,*) iT,temper(iT),rmsave(iT,iprot),sumlik(iT,iprot),
182 c-------------------------------------------------------------------------------
183 subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
185 integer n,nf,uiparm(1)
186 double precision x(n),g(n),urparm(1),ufparm
189 double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
193 call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
195 call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
196 g(i)=(fplus-fminus)/delta2
197 c write(2,*) i,fplus,fminus,g(i)
201 c-------------------------------------------------------------------------------
202 double precision function fdum()
206 c-------------------------------------------------------------------------------
207 logical function stopx(idum)
212 c-------------------------------------------------------------------------------
215 include "COMMON.CALC"
217 double precision x(n)
218 double precision fabs
221 if (mask(i).gt.0) then
223 weight(i)=fabs(x(ii))
228 c-------------------------------------------------------------------------------
231 include "COMMON.CALC"
233 double precision x(maxene)
236 if (mask(i).gt.0) then
244 c-------------------------------------------------------------------------------
245 double precision function fabs(x)
247 double precision a /1.0d4/
248 if (dabs(x).gt.1.0d-2) then
251 fabs = dlog(dexp(a*x)+dexp(-a*x))/a