source/maxlik/src_MD_T_maxlik-NEWCORR/make_distrib.f

   1       subroutine make_distrib(icl,iprot)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'COMMON.ENERGIES'
   6       include 'COMMON.ALLPROT'
   7       include 'COMMON.CLASSES'
   8       include 'COMMON.WEIGHTDER'
   9       include 'COMMON.VMCPAR'
  10       include 'COMMON.PROTNAME'
  11       character*6 liczba
  12       character*2 lbatch
  13       character*256 nazwa
  14       double precision emin,emax,erange,binsize,aux,sumpart,fac,
  15      &  distr(0:1000,maxT,maxbatch)
  16       integer i,ind,j,k,nbin,ib,iprot,icl,ibatch
  17       integer ilen
  18       external ilen
  19       logical not_done
  20       print *,"Make_Distrib: iprot",iprot," icl",icl
  21       emin=1.0e20
  22       emax=-1.0e20
  23       do i=1,ntot(iprot)
  24         if (iscore(i,0,iprot).eq.icl) then
  25           if (emin.gt.e_total(i,iprot)) then
  26             emin=e_total(i,iprot)
  27           endif
  28           if (emax.lt.e_total(i,iprot)) then
  29             emax=e_total(i,iprot)
  30           endif
  31         endif
  32       enddo
  33       do ibatch=1,nbatch(iprot)
  34       do ib=1,nbeta(ibatch,iprot)
  35       fac=betaT(2*ib,ibatch,iprot)
  36       erange=emax-emin
  37       binsize=1.0d0
  38       nbin=erange/binsize
  39       print *,"emin",emin," emax",emax," binsize",binsize," nbin",nbin
  40       do i=0,nbin
  41         distr(i,ib,ibatch)=0.0D0
  42       enddo
  43       sumpart=0.0d0
  44       do i=1,ntot(iprot)
  45         if (iscore(i,0,iprot).eq.icl) then
  46           ind=dint((e_total(i,iprot)-emin)/binsize)
  47           aux=fac*(e_total(i,iprot)-emin)
  48           if (aux.lt.50.0d0) then
  49             aux=dexp(-aux)
  50             distr(ind,ib,ibatch)=distr(ind,ib,ibatch)+aux
  51             sumpart=sumpart+aux
  52           endif
  53         endif
  54       enddo
  55       do i=0,nbin
  56         distr(i,ib,ibatch)=distr(i,ib,ibatch)/sumpart
  57       enddo
  58       enddo
  59       enddo
  60 c
  61       not_done=.true.
  62       do while (not_done)
  63         not_done=.false.
  64         do ibatch=1,nbatch(iprot)
  65         do ib=1,nbeta(ibatch,iprot)
  66           if (distr(nbin,ib,ibatch).ge.1.0d-5) goto 10
  67         enddo
  68         enddo
  69         nbin=nbin-1
  70         not_done=.true.
  71    10   continue
  72       enddo
  73       write (liczba,'(bz,i6.6)') icl
  74       do ibatch=1,nbatch(iprot)
  75       write(lbatch,'(bz,i2.2)') ibatch
  76       nazwa=protname(iprot)(:ilen(protname(iprot)))
  77      &   //'.'//lbatch//'.'//liczba//
  78      &   '.distr'
  79       open(88,file=nazwa,status="unknown")
  80       write (88,'(a10,20f10.1)') '# Energy',(betaT(ib,ibatch,iprot),
  81      &  ib=1,nbeta(ibatch,iprot))
  82       do i=0,nbin
  83         write (88,'(21f10.3)') binsize*floor(emin/binsize)
  84      &   +i*binsize,(distr(i,ib,ibatch),ib=1,nbeta(ibatch,iprot))
  85       enddo
  86       close(88)
  87       enddo ! ibatch
  88       return
  89       end