X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Fmaxlik%2Fsrc_CSA%2Fmaxlik-opt-multprot.f;fp=source%2Fmaxlik%2Fsrc_CSA%2Fmaxlik-opt-multprot.f;h=0000000000000000000000000000000000000000;hp=e51674f9d1fa1e5bdf16f54256623ac1ac7bad89;hb=94c9f54ff76d22351c0cd2f1b01e809515163940;hpb=2a226bfc86eabc6e4eae0c3ad1cbc3cb5417a05a diff --git a/source/maxlik/src_CSA/maxlik-opt-multprot.f b/source/maxlik/src_CSA/maxlik-opt-multprot.f deleted file mode 100644 index e51674f..0000000 --- a/source/maxlik/src_CSA/maxlik-opt-multprot.f +++ /dev/null @@ -1,254 +0,0 @@ - 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