source/maxlik/src_CSA/maxlik-opt-multprot.f

   1       implicit none
   2       include "DIMENSIONS"
   3       include "COMMON.CALC"
   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
   7       external fdum,funclik
   8       double precision quot,quotl,f,T0 /300.0d0/
   9       character*8 ename(maxene)
  10       character*480 karta
  11       common /names/ ename
  12       character*32 protfile(maxprot)
  13 c      print *,"start"
  14       read(1,*) nprot,nene,sigma2,temper(1),ienecheck
  15       nT = 1
  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)
  23
  24       do iprot=1,nprot
  25
  26       read (1,'(2a)') protname(iprot)
  27       read (1,'(2a)') protfile(iprot)
  28
  29       write (2,*) "Reading data for protein ",protname(iprot)
  30       write (2,*) "File: ",protfile(iprot)
  31
  32       open (7,file=protfile(iprot),status='old')
  33
  34       i=0
  35       do
  36 c        print *,"i=",i
  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)
  41         i=i+1
  42       enddo
  43    10 continue
  44       nconf(iprot)=i
  45       do i=1,nconf(iprot)
  46         entfac(i,iprot)=0.0d0
  47       enddo
  48       write (2,*) "Protein:",iprot, " nconf",nconf(iprot)
  49
  50       close(7)
  51
  52       enddo ! iprot
  53
  54        write (2,'(i5,2x,a4,f10.5,2i5)')
  55      &   (i,ename(i),weight(i),iweight(i),mask(i),i=1,nene)
  56       sigma2=sigma2*sigma2
  57 c Transfer weights to variables
  58       call w2x(n,x)
  59 c      write (2,*) "BEFORE funclik: x",(x(i),i=1,n)
  60       call funclik(n,x,nf,f,uiparm,urparm,fdum)
  61
  62       if (ienecheck.gt.0) then
  63       write (2,*) "Checking energies"
  64       do iprot=1,nprot
  65       write (2,*) "Protein",iprot," name",protname(iprot)
  66       write (2,'(a5,2a15)') "Conf","UNRES-calc E","Initial E"
  67       do i=1,nconf(iprot)
  68         write (2,'(i5,2e15.5)') i,ener0(i,iprot),ener(i,iprot)
  69       enddo
  70       enddo
  71       endif
  72
  73       do iprot=1,nprot
  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)
  77       enddo
  78       write (2,*) "Initial target function f",f
  79
  80       call minsumsl(n,x,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)
  84       write (2,*) "f",f
  85       call funclik(n,x,nf,f,uiparm,urparm,fdum)
  86
  87       do iprot=1,nprot
  88       write (2,*) "Protein",iprot," name",protname(iprot)
  89       write (2,*) "rmsave",(rmsave(i,iprot),i=1,nT)
  90       enddo
  91       write (2,*) "Final target function f",f
  92
  93 c      do i=10,30
  94 c      do k=10,30
  95 c      weight(6)=0.1d0*i
  96 c      weight(1)=0.1d0*k
  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)
 100 c      write (2,*) "f",f
 101 c      enddo
 102 c      enddo
 103       stop
 104       end
 105 c-------------------------------------------------------------------------------
 106       subroutine funclik(n,x,nf,f,uiparm,urparm,ufparm)
 107       implicit none
 108       include "DIMENSIONS"
 109       include "COMMON.CALC"
 110       character*8 ename(maxene)
 111       common /names/ ename
 112       integer n,nf,iprot
 113       double precision f
 114       double precision x(n)
 115       integer uiparm
 116       double precision ufparm
 117       external ufparm
 118       double precision ww(maxene)
 119       double precision urparm(1)
 120       integer it,i,j
 121       double precision beta,over,boltz,sumQ,emin,ee,
 122      &   sumover
 123 c      write (2,*) "funclik: x",(x(i),i=1,n)
 124       call x2w(n,x)
 125       f=0.0d0
 126
 127
 128 c      do iT=1,nT
 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),
 132 c     &    i=21,20+n)
 133 c      enddo
 134
 135       do iprot=1,nprot
 136
 137       do iT=1,nT
 138         do i=1,nene
 139           ww(i)=weight(i)
 140         enddo
 141         beta=1.0d0/(temper(iT)*1.987D-3)
 142 c        write (2,*) "iprot",iprot," iT",iT," temper",temper(iT),
 143 c     &       " beta",beta
 144 c        write (2,*) "beta",beta
 145         emin=1.0d10
 146         do i=1,nconf(iprot)
 147           ener(i,iprot)=0.0d0
 148           do j=1,nene
 149             ener(i,iprot)=ener(i,iprot)+ww(j)*enetb(j,i,iprot)
 150           enddo
 151           ee = ener(i,iprot)-entfac(i,iprot)/beta
 152           if (ee.lt.emin) emin=ee
 153         enddo
 154         rmsave(it,iprot)=0.0d0
 155         sumlik(it,iprot)=0.0d0
 156         sumQ=0.0d0
 157         sumover=0.0d0
 158         do i=1,nconf(iprot)
 159           over=dexp(-0.5d0*rmstab(i,iprot)**2/sigma2)
 160 c          if (temper(iT).gt.340.0d0) over=1.0d0-over
 161           sumover=sumover+over
 162           boltz=-beta*(ener(i,iprot)-emin)+entfac(i,iprot)
 163 c          write (2,*) i,ener(i),entfac(i),rmstab(i),over,boltz,
 164 c     &      dexp(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
 168         enddo
 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),
 174 c     &     sumQ,sumover
 175         f=f-sumlik(iT,iprot)
 176       enddo
 177
 178       enddo ! iprot
 179
 180       return
 181       end
 182 c-------------------------------------------------------------------------------
 183       subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
 184       implicit none
 185       integer n,nf,uiparm(1)
 186       double precision x(n),g(n),urparm(1),ufparm
 187       external ufparm
 188       integer i
 189       double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
 190       do i=1,n
 191         xi=x(i)
 192         x(i)=xi+delta
 193         call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
 194         x(i)=xi-delta
 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)
 198       enddo
 199       return
 200       end
 201 c-------------------------------------------------------------------------------
 202       double precision function fdum()
 203       fdum=0.0d0
 204       return
 205       end
 206 c-------------------------------------------------------------------------------
 207       logical function stopx(idum)
 208       integer idum
 209       stopx=.false.
 210       return
 211       end
 212 c-------------------------------------------------------------------------------
 213       subroutine x2w(n,x)
 214       include "DIMENSIONS"
 215       include "COMMON.CALC"
 216       integer n,i,ii
 217       double precision x(n)
 218       double precision fabs
 219       ii=0
 220       do i=1,nene
 221         if (mask(i).gt.0) then
 222           ii=ii+1
 223           weight(i)=fabs(x(ii))
 224         endif
 225       enddo
 226       return
 227       end
 228 c-------------------------------------------------------------------------------
 229       subroutine w2x(n,x)
 230       include "DIMENSIONS"
 231       include "COMMON.CALC"
 232       integer n,i,ii
 233       double precision x(maxene)
 234       ii=0
 235       do i=1,nene
 236         if (mask(i).gt.0) then
 237           ii=ii+1
 238           x(ii)=weight(i)
 239         endif
 240       enddo
 241       n=ii
 242       return
 243       end
 244 c-------------------------------------------------------------------------------
 245       double precision function fabs(x)
 246       double precision x
 247       double precision a /1.0d4/
 248       if (dabs(x).gt.1.0d-2) then
 249         fabs = dabs(x)
 250       else
 251         fabs = dlog(dexp(a*x)+dexp(-a*x))/a
 252       endif
 253       return
 254       end