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

   1       implicit none
   2       include "DIMENSIONS"
   3       include "COMMON.CALC"
   4       integer i,j,k,l,ii,nf,n,uiparm(1)
   5       double precision x(maxene)
   6       double precision rmsave(maxT),fdum,rr
   7       external fdum,funclik
   8       double precision quot,quotl,f,T0 /300.0d0/
   9       character*8 ename(maxene)
  10       common /names/ ename
  11       character*1 restyp(nbase)
  12       data restyp /'A','G','C','T','U'/
  13       character maskchar(0:1) /' ','*'/
  14       character*80 Naresfile,Fracfile
  15       double precision g(100)
  16 c      print *,"start"
  17       read(1,*) nene,sigma2,wsq
  18       write (2,*) "nene",nene," nT",nT," sigma",sigma2
  19       read(1,*) (ename(i),i=1,nene)
  20       read(1,*) (weight(i),i=1,nene)
  21       read(1,*) (iweight(i),i=1,nene)
  22       read(1,*) (mask(i),i=1,nene)
  23       do i=1,nnbase
  24         read(1,*) sig0(i),(weightel(3*(i-1)+j),maskel(3*(i-1)+j),j=1,3)
  25       enddo
  26       read(1,'(a)') Naresfile
  27       read(1,'(a)') Fracfile
  28       close(1)
  29       open(7,file=Naresfile,status='old')
  30       i=0
  31       do
  32 c        print *,"i=",i
  33         read(7,*,end=10,err=10) ii,(enetb(j,i+1),j=1,nene),
  34      &             ((eletb(3*(j-1)+k,i+1),j=1,nnbase),k=1,3),
  35      &                          entfac(i+1),
  36      &                          qtab(i+1),rmstab(i+1),rgytab(i+1),
  37      &                          fpair(i+1),rr,fdup(i+1)
  38         i=i+1
  39       enddo
  40    10 continue
  41       nconf=i
  42       close(7)
  43       write (2,*) "nconf",nconf
  44       write (2,'(/"Initial energy-term weights (* optimized)")')
  45       write (2,'(i5,2x,a4,f10.5,1x,a1,i5)')
  46      &   (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
  47       ii=0
  48       write (2,*)
  49      & "Initial base-dipole-interaction parameters (* optimizable)"
  50       do i=1,nbase
  51         do j=1,i
  52           ii=ii+1
  53           write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
  54      &     sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
  55      &     k=1,3)
  56         enddo
  57       enddo
  58       write (2,*) "sigma",sigma2," wsq",wsq
  59       sigma2=sigma2*sigma2
  60       close(7)
  61       open(7,file=Fracfile,status='old')
  62       i=0
  63       do
  64         read(7,*,end=11,err=11) temper(i+1),frac(i+1)
  65         i=i+1
  66       enddo
  67    11 continue
  68       nT=i
  69       close(7)
  70       write (2,*) "Fractions of base pairs"
  71       do i=1,nT
  72         write (2,'(i5,f8.1,f10.5)') i,temper(i),frac(i)
  73       enddo
  74 c Transfer weights to variables
  75       call w2x(n,x)
  76       write (2,*) "n",n
  77 c Compute the temperature scale factors
  78       do i=1,nT
  79         ft(1,i)=1.0d0
  80         quot=temper(i)/T0
  81         quotl=1.0d0
  82         do l=2,2
  83           quotl=quotl*quot
  84           fT(l,i)=1.12692801104297249644d0/
  85      &             dlog(dexp(quotl)+dexp(-quotl))
  86         enddo
  87       enddo
  88 c      do i=1,nT
  89 c        write (2,'(i5,f8.3,f10.5)') (i,(ft(j,i),j=1,2),i=1,nT)
  90 c      enddo
  91
  92       call funclik(n,x,nf,f,uiparm,rmsave,fdum)
  93       write (2,*) "f",f
  94       write(2,'(/a3,a8,3a10)')" Nr","    Temp","    rmsave",
  95      &   "      fave"," fave(exp)"
  96       do i=1,nT
  97         write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
  98      &   fave(i),frac(i)
  99       enddo
 100 c      call grad(n,x,nf,g,uiparm,rmsave,fdum)
 101 c      write (2,*) "rmsave",(rmsave(i),i=1,nT),"f",f
 102 c      stop
 103       call minsumsl(n,x,f)
 104       call funclik(n,x,nf,f,uiparm,rmsave,fdum)
 105       write (2,'(/"Final parameters (",i5,")")') n
 106       do i=1,n
 107         write (2,'(i5,f10.5)') i,x(i)
 108       enddo
 109       write (2,'(/"Energy-term weights (* optimized)")')
 110       write (2,'(i5,2x,a4,f10.5,1x,a1,i5)')
 111      &   (j,ename(j),weight(j),maskchar(mask(j)),iweight(j),j=1,nene)
 112       write (2,*) "Base-dipole-interaction parameters (* optimized)"
 113       ii=0
 114       do i=1,nbase
 115         do j=1,i
 116           ii=ii+1
 117           write (2,'(2a2,f10.5,3(f10.5,1x,a1))') restyp(i),restyp(j),
 118      &     sig0(ii),(weightel(3*(ii-1)+k),maskchar(maskel(3*(ii-1)+k)),
 119      &     k=1,3)
 120         enddo
 121       enddo
 122       write (2,*) "f",f
 123       write(2,'(/a3,a8,3a10)')" Nr","    Temp","    rmsave",
 124      &   "      fave"," fave(exp)"
 125       do i=1,nT
 126         write (2,'(i3,f8.1,3f10.5)') i,temper(i),rmsave(i),
 127      &   fave(i),frac(i)
 128       enddo
 129
 130 c      do i=10,30
 131 c      do k=10,30
 132 c      weight(6)=0.1d0*i
 133 c      weight(1)=0.1d0*k
 134 c       write (2,'(i5,2x,a4,f10.5,i5)')
 135 c     &   (j,ename(j),weight(j),iweight(j),j=1,nene)
 136 c      call funclik(nene,weight,nf,f,uiparm,rmsave,fdum)
 137 c      write (2,*) "f",f
 138 c      enddo
 139 c      enddo
 140       stop
 141       end
 142 c-------------------------------------------------------------------------------
 143       subroutine funclik(n,x,nf,f,uiparm,rmsave,ufparm)
 144       implicit none
 145       include "DIMENSIONS"
 146       include "COMMON.CALC"
 147       character*8 ename(maxene)
 148       common /names/ ename
 149       integer n,nf
 150       double precision f,loglik,chisq
 151       double precision x(n)
 152       integer uiparm
 153       double precision ufparm
 154       external ufparm
 155       double precision ww(maxene),sumlik(maxT),rmsave(maxT)
 156       integer it,i,j,k
 157       double precision beta,over,boltz,sumQ,ener(maxconf),emin,ee,enel,
 158      &   sumover,eboltz
 159       call x2w(n,x)
 160       loglik=0.0d0
 161       chisq=0.0d0
 162       do iT=1,nT
 163         do i=1,nene
 164           ww(i)=weight(i)*ft(iweight(i),iT)
 165         enddo
 166 c        write (2,*) "iT",iT," temper",temper(iT)," beta",beta
 167 c        write (2,'(i5,2x,a4,2f10.5,i5,f10.5)')
 168 c     &   (i,ename(i),weight(i),ww(i),iweight(i),ft(iweight(i),iT),
 169 c     &    i=1,n)
 170         beta=1.0d0/(temper(iT)*1.987D-3)
 171 c        write (2,*) "beta",beta
 172         emin=1.0d10
 173         do i=1,nconf
 174           enel=0.0d0
 175           do j=1,nnbase
 176 c            write (2,*)
 177 c     &       i,j,sig0(j),(weightel(3*(j-1)+k),eletb(3*(j-1)+k,i),k=1,3)
 178             enel=enel+weightel(3*(j-1)+1)*sig0(j)**6*eletb(3*(j-1)+1,i)
 179      &               +weightel(3*(j-1)+2)*sig0(j)**3*eletb(3*(j-1)+2,i)
 180      &               +weightel(3*j)*sig0(j)**6*eletb(3*j,i)
 181           enddo
 182 c          write (2,*) i,enel,enetb(6,i)
 183           enetb(6,i)=enel
 184           ener(i)=0.0d0
 185           do j=1,nene
 186             ener(i)=ener(i)+ww(j)*enetb(j,i)
 187           enddo
 188           ee = ener(i)-entfac(i)/beta
 189           if (ee.lt.emin) emin=ee
 190         enddo
 191         rmsave(it)=0.0d0
 192         sumlik(it)=0.0d0
 193         fave(it)=0.0d0
 194         sumQ=0.0d0
 195         sumover=0.0d0
 196         do i=1,nconf
 197           boltz=-beta*(ener(i)-emin)+entfac(i)
 198           eboltz=dexp(boltz)
 199           rmsave(iT)=rmsave(iT)+rmstab(i)*dexp(boltz)
 200           over=dexp(-0.5d0*rmstab(i)**2/sigma2)
 201           if (frac(iT).lt.0.5d0) over=1.0d0-over
 202           sumover=sumover+over
 203 c          write (2,*) i,ener(i),entfac(i),rmstab(i),fpair(i),fdup(i),
 204 c     &      over,boltz,eboltz
 205           fave(iT)=fave(iT)+fdup(i)*eboltz
 206           sumQ=sumQ+eboltz
 207           sumlik(iT)=sumlik(iT)+over*boltz
 208         enddo
 209         fave(iT)=fave(iT)/sumq
 210         sumlik(it)=sumlik(iT)-dlog(sumQ)*sumover
 211         rmsave(iT)=rmsave(iT)/sumQ
 212         if (frac(iT).gt.0.95d0 .or. frac(iT).lt.0.05d0) then
 213           loglik=loglik-sumlik(iT)
 214         endif
 215         chisq=chisq+(frac(iT)-fave(iT))**2
 216 c        write (2,*) iT,temper(iT),rmsave(iT),sumlik(iT),sumQ,sumover,
 217 c     &   fave(iT),frac(iT)
 218 c        call flush(2)
 219       enddo
 220       f = loglik/nconf+wsq*chisq
 221 c      write (2,*) "loglik",loglik/nconf," chisq",chisq," f",f
 222       return
 223       end
 224 c-------------------------------------------------------------------------------
 225       subroutine grad(n,x,nf,g,uiparm,urparm,ufparm)
 226       implicit none
 227       integer n,nf,uiparm(1)
 228       double precision x(n),g(n),urparm(1),ufparm
 229       external ufparm
 230       integer i
 231       double precision xi,fplus,fminus,delta/1.0d-9/,delta2/2.0d-9/
 232       do i=1,n
 233         xi=x(i)
 234         x(i)=xi+delta
 235         call funclik(n,x,nf,fplus,uiparm,urparm,ufparm)
 236         x(i)=xi-delta
 237         call funclik(n,x,nf,fminus,uiparm,urparm,ufparm)
 238         g(i)=(fplus-fminus)/delta2
 239 c        write(2,*) i,fplus,fminus,g(i)
 240       enddo
 241       return
 242       end
 243 c-------------------------------------------------------------------------------
 244       double precision function fdum()
 245       fdum=0.0d0
 246       return
 247       end
 248 c-------------------------------------------------------------------------------
 249       logical function stopx(idum)
 250       integer idum
 251       stopx=.false.
 252       return
 253       end
 254 c-------------------------------------------------------------------------------
 255       subroutine x2w(n,x)
 256       include "DIMENSIONS"
 257       include "COMMON.CALC"
 258       double precision fabs
 259       integer n,i,ii
 260       double precision x(n)
 261       ii=0
 262       do i=1,nene
 263         if (mask(i).gt.0) then
 264           ii=ii+1
 265           weight(i)=fabs(x(ii))
 266         endif
 267       enddo
 268       do i=1,nnbase
 269         if (maskel(3*(i-1)+1).gt.0) then
 270           ii=ii+1
 271           weightel(3*(i-1)+1)=-fabs(x(ii))
 272         endif
 273         if (maskel(3*(i-1)+2).gt.0) then
 274           ii=ii+1
 275           weightel(3*(i-1)+2)=x(ii)
 276         endif
 277         if (maskel(3*i).gt.0) then
 278           ii=ii+1
 279           weightel(3*i)=-fabs(x(ii))
 280         endif
 281       enddo
 282       return
 283       end
 284 c-------------------------------------------------------------------------------
 285       subroutine w2x(n,x)
 286       include "DIMENSIONS"
 287       include "COMMON.CALC"
 288       integer n,i,ii
 289       double precision x(maxene+3*nnbase)
 290       ii=0
 291       do i=1,nene
 292 c        write (2,*) "i",i," mask",mask(i)," ii",ii
 293         if (mask(i).gt.0) then
 294           ii=ii+1
 295           x(ii)=weight(i)
 296         endif
 297       enddo
 298       do i=1,3*nnbase
 299 c        write (2,*) "i",i," maskel",maskel(i)," ii",ii
 300         if (maskel(i).gt.0) then
 301           ii=ii+1
 302           x(ii)=weightel(i)
 303         endif
 304       enddo
 305       n=ii
 306 c      write (2,*) "W2X: n=",n
 307       return
 308       end
 309 c------------------------------------------------------------------------------------
 310       double precision function fabs(x)
 311       double precision x
 312       double precision a /100.0d0/
 313       if (dabs(x).gt.1.0d-2) then
 314         fabs = dabs(x)
 315       else
 316         fabs = dlog(dexp(a*x)+dexp(-a*x))/a
 317       endif
 318       return
 319       end