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