+++ /dev/null
- implicit none
- include "DIMENSIONS"
- include "COMMON.CALC"
- integer i,j,k,l,ii,iprot,nf,n,uiparm(1),ienecheck
- double precision x(maxene)
- double precision urparm(1),fdum,rjunk
- external fdum,funclik
- double precision quot,quotl,f,T0 /300.0d0/
- character*8 ename(maxene)
- character*480 karta
- common /names/ ename
- character*32 protfile(maxprot)
-c print *,"start"
- read(1,*) nprot,nene,sigma2,temper(1),ienecheck
- nT = 1
- write (2,*) "nene",nene," nT",nT," sigma",sigma2,
- & " ienechek",ienecheck
- 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)
-
- do iprot=1,nprot
-
- read (1,'(2a)') protname(iprot)
- read (1,'(2a)') protfile(iprot)
-
- write (2,*) "Reading data for protein ",protname(iprot)
- write (2,*) "File: ",protfile(iprot)
-
- open (7,file=protfile(iprot),status='old')
-
- i=0
- do
-c print *,"i=",i
- read(7,'(a)',end=10) karta
- if (index(karta,'#').gt.0) cycle
- read(karta,*,end=10,err=10) ii,(enetb(j,i+1,iprot),j=1,nene),
- & ener0(i+1,iprot),rmstab(i+1,iprot),rjunk,rjunk,qtab(i+1,iprot)
- i=i+1
- enddo
- 10 continue
- nconf(iprot)=i
- do i=1,nconf(iprot)
- entfac(i,iprot)=0.0d0
- enddo
- write (2,*) "Protein:",iprot, " nconf",nconf(iprot)
-
- close(7)
-
- enddo ! iprot
-
- 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)
-c write (2,*) "BEFORE funclik: x",(x(i),i=1,n)
- call funclik(n,x,nf,f,uiparm,urparm,fdum)
-
- if (ienecheck.gt.0) then
- write (2,*) "Checking energies"
- do iprot=1,nprot
- write (2,*) "Protein",iprot," name",protname(iprot)
- write (2,'(a5,2a15)') "Conf","UNRES-calc E","Initial E"
- do i=1,nconf(iprot)
- write (2,'(i5,2e15.5)') i,ener0(i,iprot),ener(i,iprot)
- enddo
- enddo
- endif
-
- do iprot=1,nprot
- write (2,*) "Protein",iprot," name",protname(iprot)
- write (2,*) "Initial average rmsd(s)"
- write (2,*) "rmsave",(rmsave(i,iprot),i=1,nT)
- enddo
- write (2,*) "Initial target function 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,urparm,fdum)
-
- do iprot=1,nprot
- write (2,*) "Protein",iprot," name",protname(iprot)
- write (2,*) "rmsave",(rmsave(i,iprot),i=1,nT)
- enddo
- write (2,*) "Final target function 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,urparm,ufparm)
- implicit none
- include "DIMENSIONS"
- include "COMMON.CALC"
- character*8 ename(maxene)
- common /names/ ename
- integer n,nf,iprot
- double precision f
- double precision x(n)
- integer uiparm
- double precision ufparm
- external ufparm
- double precision ww(maxene)
- double precision urparm(1)
- integer it,i,j
- double precision beta,over,boltz,sumQ,emin,ee,
- & sumover
-c write (2,*) "funclik: x",(x(i),i=1,n)
- call x2w(n,x)
- f=0.0d0
-
-
-c do iT=1,nT
-c write (2,'(i5,2x,a4,f10.5,i5)')
-c & (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
-c & (i,ename(i),weight(i),iweight(i),
-c & i=21,20+n)
-c enddo
-
- do iprot=1,nprot
-
- do iT=1,nT
- do i=1,nene
- ww(i)=weight(i)
- enddo
- beta=1.0d0/(temper(iT)*1.987D-3)
-c write (2,*) "iprot",iprot," iT",iT," temper",temper(iT),
-c & " beta",beta
-c write (2,*) "beta",beta
- emin=1.0d10
- do i=1,nconf(iprot)
- ener(i,iprot)=0.0d0
- do j=1,nene
- ener(i,iprot)=ener(i,iprot)+ww(j)*enetb(j,i,iprot)
- enddo
- ee = ener(i,iprot)-entfac(i,iprot)/beta
- if (ee.lt.emin) emin=ee
- enddo
- rmsave(it,iprot)=0.0d0
- sumlik(it,iprot)=0.0d0
- sumQ=0.0d0
- sumover=0.0d0
- do i=1,nconf(iprot)
- over=dexp(-0.5d0*rmstab(i,iprot)**2/sigma2)
-c if (temper(iT).gt.340.0d0) over=1.0d0-over
- sumover=sumover+over
- boltz=-beta*(ener(i,iprot)-emin)+entfac(i,iprot)
-c write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
-c & dexp(boltz)
- sumQ=sumQ+dexp(boltz)
- rmsave(iT,iprot)=rmsave(iT,iprot)+rmstab(i,iprot)*dexp(boltz)
- sumlik(iT,iprot)=sumlik(iT,iprot)+over*boltz
- enddo
- sumlik(it,iprot)=sumlik(iT,iprot)-dlog(sumQ)*sumover
- rmsave(iT,iprot)=rmsave(iT,iprot)/sumQ
-c write (2,*) iprot,iT,temper(iT),rmsave(iT,iprot),
-c & sumlik(iT,iprot),sumQ,sumover
-c write (2,*) iT,temper(iT),rmsave(iT,iprot),sumlik(iT,iprot),
-c & sumQ,sumover
- f=f-sumlik(iT,iprot)
- enddo
-
- enddo ! iprot
-
- 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 /1.0d4/
- if (dabs(x).gt.1.0d-2) then
- fabs = dabs(x)
- else
- fabs = dlog(dexp(a*x)+dexp(-a*x))/a
- endif
- return
- end
projects / unres.git / blobdiff