source/maxlik/src_FPy.org/proc.F.org

   1       subroutine proc_data(nvarr,x,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
   8       include "COMMON.MPI"
   9 #endif
  10       include "COMMON.CONTROL"
  11       include "COMMON.ENERGIES"
  12       include "COMMON.VMCPAR"
  13       include "COMMON.OPTIM"
  14       include "COMMON.WEIGHTS"
  15       include "COMMON.PROTNAME"
  16       include "COMMON.IOUNITS"
  17       include "COMMON.FFIELD"
  18       include "COMMON.PROTFILES"
  19       include "COMMON.COMPAR"
  20       include "COMMON.CHAIN"
  21       include "COMMON.CLASSES"
  22       include "COMMON.THERMAL"
  23       include "COMMON.WEIGHTDER"
  24       include "COMMON.EREF"
  25       include "COMMON.SBRIDGE"
  26       include "COMMON.INTERACT"
  27       double precision fac0,t0
  28       double precision fT(5),ftprim(5),ftbis(5),quot,quotl,quotl1,
  29      &  denom,kfac,kfacl,logfac,eplus,eminus,tanhT,entfacmax,
  30      &  sument,sumentsq,aux
  31       real csingle(3,maxres2)
  32       double precision creff(3,maxres2),cc(3,maxres2)
  33       double precision przes(3),obrot(3,3),etoti
  34       logical non_conv
  35       integer wykl
  36       integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng,
  37      &  inn,in,iref,ib,ibatch,num
  38       integer l1,l2,ncon_all
  39       double precision rcyfr
  40       integer jj,istart_conf,iend_conf,ipass_conf,iii
  41       integer ilen,iroof
  42       external ilen,iroof
  43       double precision ej,emaxj,sumP,sumP_t
  44       double precision fac
  45       double precision x(max_paropt)
  46       double precision wconf(maxstr_proc)
  47       logical LPRN
  48       integer maxl1,maxl2
  49       double precision rms,rmsnat,qwolynes,qwolyness
  50       T0=3.0d2
  51       fac0=1.0d0/(T0*Rgas)
  52       kfac=2.4d0
  53       write (iout,*) "fac0",fac0
  54 c
  55 c Determine the number of variables and the initial vector of optimizable
  56 c parameters
  57 c
  58       lprn=.false.
  59       if (restart_flag) then
  60         call w2x(nvarr,x_orig,*11)
  61         do i=1,n_ene
  62           ww_oorig(i)=ww(i)
  63         enddo
  64         open (88,file=prefix(:ilen(prefix))//'.restin',status='unknown')
  65         read(88,*,end=99,err=99) (x(i),i=1,nvarr)
  66         close(88)
  67         write (iout,*)
  68         write (iout,*)
  69      &    "Variables replaced with values from the restart file."
  70         do i=1,nvarr
  71           write (iout,'(i5,f15.5)') i,x(i)
  72         enddo
  73         write (iout,*)
  74         call x2w(nvarr,x)
  75         goto 88
  76   99    write (iout,*)
  77      &    "Error - restart file nonexistent, incompatible or damaged."
  78 #ifdef MPI
  79         call MPI_Finalize(Ierror)
  80 #endif
  81         stop
  82      &    "Error - restart file nonexistent, incompatible or damaged."
  83       else
  84         call w2x(nvarr,x,*11)
  85         do i=1,n_ene
  86           ww_oorig(i)=ww(i)
  87         enddo
  88         do i=1,nvarr
  89           x_orig(i)=x(i)
  90         enddo
  91       endif
  92   88  continue
  93 c
  94 c Setup weights for UNRES
  95 c
  96       wsc=ww(1)
  97       wscp=ww(2)
  98       welec=ww(3)
  99       wcorr=ww(4)
 100       wcorr5=ww(5)
 101       wcorr6=ww(6)
 102       wel_loc=ww(7)
 103       wturn3=ww(8)
 104       wturn4=ww(9)
 105       wturn6=ww(10)
 106       wang=ww(11)
 107       wscloc=ww(12)
 108       wtor=ww(13)
 109       wtor_d=ww(14)
 110       wstrain=ww(15)
 111       wvdwpp=ww(16)
 112       wbond=ww(17)
 113       wsccor=ww(19)
 114 #ifdef SCP14
 115       scal14=ww(17)/ww(2)
 116 #endif
 117       do i=1,n_ene
 118         weights(i)=ww(i)
 119       enddo
 120
 121       DO IPROT=1,NPROT
 122
 123       call restore_molinfo(iprot)
 124
 125       DO IBATCH=1,natlike(iprot)+2
 126
 127 c Calculate temperature-dependent weight scaling factors
 128       do ib=1,nbeta(ibatch,iprot)
 129         fac=betaT(ib,ibatch,iprot)
 130         quot=fac0/fac
 131         if (rescale_mode.eq.0) then
 132           do l=1,5
 133             fT(l)=1.0d0
 134             fTprim(l)=0.0d0
 135             fTbis(l)=0.0d0
 136           enddo
 137         else if (rescale_mode.eq.1) then
 138           quotl=1.0d0
 139           kfacl=1.0d0
 140           do l=1,5
 141             quotl1=quotl
 142             quotl=quotl*quot
 143             kfacl=kfacl*kfac
 144             denom=kfacl-1.0d0+quotl
 145             fT(l)=kfacl/denom
 146             ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 147             ftbis(l)=l*kfacl*quotl1*
 148      &       (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 149             write (iout,*) "ib",ib," beta",betaT(ib,ibatch,iprot),
 150      &      " quot",quot," l",l,fT(l),fTprim(l),fTbis(l)
 151           enddo
 152         else if (rescale_mode.eq.2) then
 153           quotl=1.0d0
 154           do l=1,5
 155             quotl1=quotl
 156             quotl=quotl*quot
 157             eplus=dexp(quotl)
 158             eminus=dexp(-quotl)
 159             logfac=1.0d0/dlog(eplus+eminus)
 160             tanhT=(eplus-eminus)/(eplus+eminus)
 161             fT(l)=1.12692801104297249644d0*logfac
 162             ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 163             ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 164      &      2*l*quotl1/T0*logfac*
 165      &      (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 166      &      +ftprim(l)*tanhT)
 167           enddo
 168         else
 169           write (iout,*) "Unknown RESCALE_MODE",rescale_mode
 170           call flush(iout)
 171           stop
 172         endif
 173         do k=1,n_ene
 174           escal(k,ib,ibatch,iprot)=1.0d0
 175           escalprim(k,ib,ibatch,iprot)=0.0d0
 176           escalbis(k,ib,ibatch,iprot)=0.0d0
 177         enddo
 178         escal(3,ib,ibatch,iprot)=ft(1)
 179         escal(4,ib,ibatch,iprot)=ft(3)
 180         escal(5,ib,ibatch,iprot)=ft(4)
 181         escal(6,ib,ibatch,iprot)=ft(5)
 182         escal(7,ib,ibatch,iprot)=ft(2)
 183         escal(8,ib,ibatch,iprot)=ft(2)
 184         escal(9,ib,ibatch,iprot)=ft(3)
 185         escal(10,ib,ibatch,iprot)=ft(5)
 186         escal(13,ib,ibatch,iprot)=ft(1)
 187         escal(14,ib,ibatch,iprot)=ft(2)
 188         escal(19,ib,ibatch,iprot)=ft(1)
 189         escalprim(3,ib,ibatch,iprot)=ftprim(1)
 190         escalprim(4,ib,ibatch,iprot)=ftprim(3)
 191         escalprim(5,ib,ibatch,iprot)=ftprim(4)
 192         escalprim(6,ib,ibatch,iprot)=ftprim(5)
 193         escalprim(7,ib,ibatch,iprot)=ftprim(2)
 194         escalprim(8,ib,ibatch,iprot)=ftprim(2)
 195         escalprim(9,ib,ibatch,iprot)=ftprim(3)
 196         escalprim(10,ib,ibatch,iprot)=ftprim(5)
 197         escalprim(13,ib,ibatch,iprot)=ftprim(1)
 198         escalprim(14,ib,ibatch,iprot)=ftprim(2)
 199         escalprim(19,ib,ibatch,iprot)=ftprim(1)
 200         escalbis(3,ib,ibatch,iprot)=ftbis(1)
 201         escalbis(4,ib,ibatch,iprot)=ftbis(3)
 202         escalbis(5,ib,ibatch,iprot)=ftbis(4)
 203         escalbis(6,ib,ibatch,iprot)=ftbis(5)
 204         escalbis(7,ib,ibatch,iprot)=ftbis(2)
 205         escalbis(8,ib,ibatch,iprot)=ftbis(2)
 206         escalbis(9,ib,ibatch,iprot)=ftbis(3)
 207         escalbis(10,ib,ibatch,iprot)=ftbis(5)
 208         escalbis(13,ib,ibatch,iprot)=ftbis(1)
 209         escalbis(14,ib,ibatch,iprot)=ftbis(2)
 210         escalbis(19,ib,ibatch,iprot)=ftbis(1)
 211       enddo
 212
 213       betmin(ibatch,iprot)=betaT(1,ibatch,iprot)
 214       betmax(ibatch,iprot)=betaT(1,ibatch,iprot)
 215       do ib=2,nbeta(ibatch,iprot)
 216         if (betaT(ib,ibatch,iprot).lt.betmin(ibatch,iprot))
 217      &     betmin(ibatch,iprot)=betaT(ib,ibatch,iprot)
 218         if (betaT(ib,ibatch,iprot).gt.betmax(ibatch,iprot))
 219      &     betmax(ibatch,iprot)=betaT(ib,ibatch,iprot)
 220         enddo
 221        write (iout,'(2(a,i5,2x),2(a,f7.2,2x))')
 222      &    "Protein",iprot," batch",ibatch,
 223      &    " betmin",betmin(ibatch,iprot)," betmax",betmax(ibatch,iprot)
 224
 225       ENDDO ! IBATCH
 226
 227       ENDDO ! IPROT
 228 c
 229 c Compute energy and entropy histograms to filter out conformations
 230 c with potntially too low or too high weights.
 231 c
 232       call determine_histograms
 233 c
 234 c Make the working list of conformations
 235 c
 236       call make_list(.true.,.true.,nvarr,x)
 237 c 3/8/2013 Calculate RMSDs/Qs and the Ptab array
 238       do iprot=1,nprot
 239
 240       call restore_molinfo(iprot)
 241
 242       open (icbase,file=bprotfiles(iprot),status="old",
 243      &    form="unformatted",access="direct",recl=lenrec(iprot))
 244       entfacmax=entfac(1,iprot)
 245       sument=entfac(1,iprot)
 246       sumentsq=entfac(1,iprot)**2
 247       do i=2,ntot_work(iprot)
 248         if (entfac(i,iprot).gt.entfacmax) entfacmax=entfac(i,iprot)
 249         sument=sument+entfac(i,iprot)
 250         sumentsq=sumentsq+entfac(i,iprot)**2
 251       enddo
 252       sument=sument/ntot_work(iprot)
 253       sumentsq=dsqrt(sumentsq/ntot_work(iprot)-sument**2)
 254       write (2,*) "entfacmax",entfacmax," entave",sument,
 255      &   " stdent",sumentsq
 256 #ifdef WEIDIST
 257 #ifdef MPI
 258       nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
 259       ipass_conf=0
 260       jj=0
 261       sumP=0.0d0
 262       do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 263       ipass_conf=ipass_conf+1
 264       write (iout,*) "Pass",ipass_conf
 265       istart_conf=i
 266       iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
 267 #else
 268       nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
 269       ipass_conf=0
 270       sumP=0.0d0
 271       do i=1,ntot(iprot),maxstr_proc
 272       ipass_conf=ipass_conf+1
 273       write (iout,*) "Pass",ipass_conf
 274       istart_conf=i
 275       iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
 276 #endif
 277 c
 278 c Read the chunk of conformations off a DA scratchfile.
 279 c
 280       call daread_ccoords(iprot,istart_conf,iend_conf)
 281
 282       IF (GEOM_AND_WHAM_WEIGHTS) THEN
 283
 284 c Compute energies
 285       do ib=1,nbeta(1,iprot)
 286         write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
 287         ii=0
 288         do iii=istart_conf,iend_conf
 289           ii=ii+1
 290           kk=tsil(iii,iprot)
 291           if (kk.eq.0) cycle
 292           etoti=0.0d0
 293           do k=1,n_ene
 294             etoti=etoti+ww(k)*escal(k,ib,1,iprot)
 295      &            *enetb(ii,k,iprot)
 296           enddo
 297           e_total_T(iii,ib)=entfac(kk,iprot)
 298      &             +betaT(ib,1,iprot)*(etoti-elowest(ib,1,iprot))
 299         enddo
 300       enddo ! ib
 301 #ifdef DEBUG
 302       write (iout,*) "Energies before Gather"
 303       do iii=1,ntot_work(iprot)
 304         write (iout,'(i5,100E12.5)') iii,
 305      &         (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
 306       enddo
 307 #endif
 308 #ifdef MPI
 309       do ib=1,nbeta(1,iprot)
 310         call MPI_AllGatherv(e_total_T(indstart(me1,iprot),ib),
 311      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
 312      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
 313      &     Comm1, IERROR)
 314       enddo ! ib
 315 #ifdef DEBUG
 316       write (iout,*) "Energies after Gather"
 317       do iii=1,ntot_work(iprot)
 318         write (iout,'(i5,100E12.5)') iii,
 319      &         (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
 320       enddo
 321 #endif
 322 #endif
 323
 324       ENDIF ! GEOM_AND_WHAM_WEIGHTS
 325
 326       ii=0
 327 c Compute distances
 328       do iii=istart_conf,iend_conf
 329         ii=ii+1
 330         k=tsil(iii,iprot)
 331         if (k.eq.0) cycle
 332         call restore_ccoords(iprot,ii)
 333 c Calculate the rmsds of the current conformations and convert them into weights
 334 c to approximate equal weighting of all points of the conformational space
 335         do k=1,2*nres
 336           do l=1,3
 337             creff(l,k)=c(l,k)
 338           enddo
 339         enddo
 340         do ib=1,nbeta(1,iprot)
 341          weirms(iii,ib)=0.0d0
 342         enddo
 343         do iref=1,ntot_work(iprot)
 344           if (iref.ne.iii) then
 345             read(icbase,rec=iref) ((csingle(l,k),l=1,3),k=1,nres),
 346      &      ((csingle(l,k),l=1,3),k=nnt+nres,nct+nres),
 347      &      nss,(ihpb(k),jhpb(k),k=1,nss),eini(iref,iprot),
 348      &      entfac(iref,iprot),
 349      &      ((nu(k,l,ii,iprot),k=1,maxdimnat),l=1,natlike(iprot))
 350             do k=1,2*nres
 351               do l=1,3
 352                 c(l,k)=csingle(l,k)
 353               enddo
 354             enddo
 355 #ifdef DEBUG
 356             do k=1,nres
 357               write (iout,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
 358      &                      (creff(l,k),l=1,3)
 359             enddo
 360 #endif
 361             if (rmscomp(iprot)) then
 362             call fitsq(rms,c(1,nstart_sup(iprot)),
 363      &                 creff(1,nstart_sup(iprot)),
 364      &                 nsup(iprot),przes,obrot,non_conv)
 365             if (non_conv) then
 366               print *,'Error: FITSQ non-convergent, iref',iref
 367               rms=1.0d2
 368             else if (rms.lt.-1.0d-6) then
 369               print *,'Error: rms^2 = ',rms,iref
 370               rms = 1.0d2
 371             else if (rms.lt.0) then
 372               rms=0.0d0
 373             else
 374               rms = dsqrt(rms)
 375             endif
 376 #ifdef DEBUG
 377             write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
 378      &          e_total_T(iref,ib),
 379      &          e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
 380      &          " rms",rms," fac",dexp(-0.5d0*(rms/sigma2(iprot))**2),
 381      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
 382      &          dexp(-0.5d0*(rms/sigma2(iprot))**2)*
 383      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 384 #endif
 385             aux=dexp(-0.5d0*(rms/sigma2(iprot))**2)
 386             else
 387             rms=qwolyness(creff,iprot)
 388 #ifdef DEBUG
 389             write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
 390      &          e_total_T(iref,ib),
 391      &          e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
 392      &          " rms",-dlog(rms)," fac",rms,
 393      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
 394      &          rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 395 #endif
 396             aux=rms
 397             endif ! rmscomp
 398           else
 399             aux=1.0d0
 400           endif
 401           do ib=1,nbeta(1,iprot)
 402             if (GEOM_AND_WHAM_WEIGHTS) then
 403               weirms(iii,ib) = weirms(iii,ib)
 404      &          +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 405             else
 406               weirms(iii,ib) = weirms(iii,ib)+aux
 407             endif
 408           enddo ! ib
 409         enddo ! iref
 410       enddo ! iii
 411 #ifdef DEBUG
 412       write (iout,*) "weirms"
 413       do iii=istart_conf,iend_conf
 414         write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1,
 415      &       nbeta(1,iprot))
 416       enddo
 417 #endif
 418       enddo ! i
 419 #endif
 420
 421 #ifdef MPI
 422       nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
 423       DO IB=1,NBETA(1,IPROT)
 424       ipass_conf=0
 425       jj=0
 426       sumP=0.0d0
 427       do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 428       ipass_conf=ipass_conf+1
 429       write (iout,*) "Pass",ipass_conf
 430       istart_conf=i
 431       iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
 432 #else
 433       nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
 434       ipass_conf=0
 435       sumP=0.0d0
 436       do i=1,ntot(iprot),maxstr_proc
 437       ipass_conf=ipass_conf+1
 438       write (iout,*) "Pass",ipass_conf
 439       istart_conf=i
 440       iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
 441 #endif
 442       ii=0
 443 c
 444 c Read the chunk of conformations off a DA scratchfile.
 445 c
 446       call daread_ccoords(iprot,istart_conf,iend_conf)
 447       write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
 448       do iii=istart_conf,iend_conf
 449         ii=ii+1
 450         jj=jj+1
 451         k=tsil(iii,iprot)
 452         call restore_ccoords(iprot,ii)
 453         Ptab(jj,ib,iprot)=0.0d0
 454         if (rmscomp(iprot)) then
 455           do iref=1,nref(ib,iprot)
 456 #ifdef DEBUG
 457             write (2,*) "nstart_sup",nstart_sup(iprot),
 458      &          " nend_sup",nend_sup(iprot)
 459             do k=nstart_sup(iprot),nend_sup(iprot)
 460               write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
 461      &          (cref(l,k,iref,ib,iprot),l=1,3)
 462             enddo
 463             do k=nstart_sup(iprot),nend_sup(iprot)
 464               write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3),
 465      &          (cref(l,k+nres,iref,ib,iprot),l=1,3)
 466             enddo
 467 #endif
 468             rms=rmsnat(ib,jj,iref,iprot)
 469 #ifdef DEBUG
 470             write (iout,*) "iii",iii," ii",ii," jj",jj," rms",rms,
 471      &            "efree",efreeref(iref,ib,iprot)
 472 #endif
 473             Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)
 474      &          +dexp(-efreeref(iref,ib,iprot)
 475      &                -0.5d0*(rms/sigma2(iprot))**2)
 476           enddo
 477 #if defined(WEIENT)
 478 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 479 c              search of conformational space.
 480           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 481           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 482 #ifdef DEBUG
 483           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 484 c     &     " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot)
 485      &     "entfac",entfac(k,iprot)," weirms",weirms(k,iprot),
 486      &     " Ptab",Ptab(jj,ib,iprot)
 487 #endif
 488 #elif defined(WEIDIST)
 489           write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 490      &     " weirms",weirms(iii,ib)," Ptab/weirms",
 491      &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 492           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 493 #elif defined(WEICLUST)
 494           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot)
 495 #endif
 496           sumP=sumP+Ptab(jj,ib,iprot)
 497         else
 498           do iref=1,nref(ib,iprot)
 499             rms = qwolynes(ib,iref,iprot)
 500 #ifdef DEBUG
 501             write (iout,*) "iii",iii," ii",ii," jj",jj,
 502      &        "iref",iref," rms",rms,-dlog(rms),
 503      &           "efree",efreeref(iref,ib,iprot)
 504 #endif
 505              Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms*
 506      &                        *dexp(-efreeref(iref,ib,iprot))
 507           enddo
 508 #if defined(WEIENT)
 509 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 510 c              search of conformational space.
 511           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 512           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 513 #ifdef DEBUG
 514           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 515      &    " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot)
 516 #endif
 517 #elif defined(WEICLUST)
 518           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot)
 519 #elif defined(WEIDIST)
 520           write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 521      &     " weirms",weirms(iii,ib)," Ptab/weirms",
 522      &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 523           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 524 #endif
 525           sumP=sumP+Ptab(jj,ib,iprot)
 526         endif
 527 #ifdef DEBUG
 528         write (iout,*) "ii",ii," jj",jj," ib",ib,
 529      &    "Ptab",Ptab(jj,ib,iprot)," sumP",sumP
 530 #endif
 531       enddo ! iii
 532       enddo ! i
 533 #ifdef DEBUG
 534       write (iout,*) "ib",ib," sumP before REDUCE",sumP
 535 #endif
 536 #ifdef MPI
 537       call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION,
 538      &        MPI_SUM,Comm1,IERROR)
 539       sumP=sumP_t
 540 #endif
 541 #ifdef DEBUG
 542       write (iout,*) "ib",ib," sumP after REDUCE",sumP
 543 #endif
 544       jj=0
 545       do iii=indstart(me1,iprot),indend(me1,iprot)
 546       jj=jj+1
 547 #ifdef DEBUG
 548       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 549 #endif
 550       Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP
 551 #ifdef DEBUG
 552       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 553 #endif
 554 #undef DEBUG
 555       enddo
 556       ENDDO ! IB
 557       close(icbase)
 558       enddo ! iprot
 559 #ifdef DEBUG
 560       write (iout,*) "ELOWEST at the end of PROC_DATA"
 561       do iprot=1,nprot
 562         do ibatch=1,natlike(iprot)+2
 563           do ib=1,nbeta(ibatch,iprot)
 564             write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
 565      &         " elowest",elowest(ib,ibatch,iprot)
 566           enddo
 567         enddo
 568       enddo
 569 #endif
 570       write (iout,*)
 571       return
 572    11 return1
 573       end
 574 c-------------------------------------------------------------------------------
 575       subroutine determine_histograms
 576       implicit none
 577       include "DIMENSIONS"
 578       include "DIMENSIONS.ZSCOPT"
 579       integer maxind
 580       parameter (maxind=1000)
 581 #ifdef MPI
 582       include "mpif.h"
 583       integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
 584       include "COMMON.MPI"
 585 #endif
 586       include "COMMON.COMPAR"
 587       include "COMMON.ENERGIES"
 588       include "COMMON.VMCPAR"
 589       include "COMMON.WEIGHTS"
 590       include "COMMON.PROTNAME"
 591       include "COMMON.IOUNITS"
 592       include "COMMON.FFIELD"
 593       include "COMMON.PROTFILES"
 594       include "COMMON.CHAIN"
 595       include "COMMON.CONTROL"
 596       include "COMMON.CLASSES"
 597       integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt
 598       double precision eminh,emaxh,tminh,tmaxh
 599       integer he(0:maxind),
 600      &  ht(0:maxind),hemax,htmax
 601 c
 602 c Construct energy and entropy histograms
 603 c
 604       write (iout,*) "Calling DETERMINE HISTOGRAMS"
 605       call flush(iout)
 606       do iprot=1,nprot
 607           eminh=1.0d20
 608           emaxh=-1.0d20
 609           tminh=1.0d20
 610           tmaxh=-1.0d20
 611           do i=0,maxind
 612             he(i)=0
 613             ht(i)=0
 614           enddo
 615           do i=1,ntot(iprot)
 616             iii=i
 617             if (eini(iii,iprot).lt.eminh)
 618      &        eminh=eini(iii,iprot)
 619             if (eini(iii,iprot).gt.emaxh)
 620      &        emaxh=eini(iii,iprot)
 621              if (entfac(iii,iprot).lt.tminh)
 622      &        tminh=entfac(iii,iprot)
 623             if (entfac(iii,iprot).gt.tmaxh)
 624      &        tmaxh=entfac(iii,iprot)
 625           enddo
 626 c          write (2,*)  "DETERMINE HISTOGRAMS 1"
 627           call flush(iout)
 628           do i=1,ntot(iprot)
 629             iii=i
 630             inde=eini(iii,iprot)-eminh
 631             indt=entfac(iii,iprot)-tminh
 632             if (inde.le.maxind) he(inde)=he(inde)+1
 633             if (indt.le.maxind) ht(indt)=ht(indt)+1
 634           enddo
 635 c          write (2,*)  "DETERMINE HISTOGRAMS 2"
 636             idelte=emaxh-eminh
 637             ideltt=tmaxh-tminh
 638             hemax=0
 639 c            write (iout,*) "idelte",idelte," ideltt",ideltt
 640             call flush(iout)
 641             do i=0,idelte
 642 c              write (iout,*) "i",i," he",he(i)," hemax",hemax
 643               call flush(iout)
 644               if (he(i).gt.hemax) hemax=he(i)
 645             enddo
 646             htmax=0
 647             do i=0,ideltt
 648 c              write (iout,*) "i",i," ht",ht(i)," htmax",htmax
 649               call flush(iout)
 650               if (ht(i).gt.htmax) htmax=ht(i)
 651             enddo
 652 c          write (2,*)  "DETERMINE HISTOGRAMS 3"
 653           call flush(iout)
 654             hemax=hemax/hefac(iprot)
 655             if (hemax.lt.hemax_low(iprot))
 656      &        hemax=hemax_low(iprot)
 657             htmax=htmax/htfac(iprot)
 658             if (htmax.lt.htmax_low(iprot))
 659      &        htmax=htmax_low(iprot)
 660             i=0
 661             do while (i.lt.idelte .and. he(i).lt.hemax)
 662               i=i+1
 663             enddo
 664 c          write (2,*)  "DETERMINE HISTOGRAMS 4"
 665           call flush(iout)
 666             e_lowb(iprot)=eminh+i
 667             i=0
 668             do while (i.lt.ideltt .and. ht(i).lt.htmax)
 669               i=i+1
 670             enddo
 671             t_lowb(iprot)=tminh+i
 672 #ifdef DEBUG
 673             write (iout,*) "protein",iprot
 674             write (iout,*) "energy histogram"
 675             do i=0,idelte
 676               write (iout,'(f15.5,i10)') eminh+i,he(i)
 677             enddo
 678             write (iout,*) "entropy histogram"
 679             do i=0,ideltt
 680               write (iout,'(f15.5,i10)') tminh+i,ht(i)
 681             enddo
 682             write (iout,*) "emin",eminh," tmin",tminh,
 683      &       " hemax",hemax," htmax",htmax,
 684      &       " e_lowb",e_lowb(iprot),
 685      &       " t_lowb",t_lowb(iprot)
 686             call flush(iout)
 687 #endif
 688       enddo
 689       return
 690       end
 691 c------------------------------------------------------------------------------
 692       subroutine imysort(n, x, ipermut)
 693       implicit none
 694       integer n
 695       integer x(n),xtemp
 696       integer ipermut(n)
 697       integer i,j,imax,ipm
 698       do i=1,n
 699         xtemp=x(i)
 700         imax=i
 701         do j=i+1,n
 702           if (x(j).lt.xtemp) then
 703             imax=j
 704             xtemp=x(j)
 705           endif
 706         enddo
 707         x(imax)=x(i)
 708         x(i)=xtemp
 709         ipm=ipermut(imax)
 710         ipermut(imax)=ipermut(i)
 711         ipermut(i)=ipm
 712       enddo
 713       return
 714       end
 715 c--------------------------------------------------------------------------
 716       subroutine write_conf_count
 717       implicit none
 718       include "DIMENSIONS"
 719       include "DIMENSIONS.ZSCOPT"
 720       include "COMMON.CLASSES"
 721       include "COMMON.COMPAR"
 722       include "COMMON.VMCPAR"
 723       include "COMMON.PROTNAME"
 724       include "COMMON.IOUNITS"
 725       include "COMMON.PROTFILES"
 726       integer iprot,ibatch,i,ii,j,kk
 727       integer ilen
 728       external ilen
 729       logical LPRN
 730       write (iout,*)
 731       write (iout,*)"Numbers of conformations after applying the cutoff"
 732       write (iout,*)
 733       do iprot=1,nprot
 734         write (iout,'(a,2x,a,i10)') "Protein",
 735      &    protname(iprot)(:ilen(protname(iprot))),
 736      &    ntot_work(iprot)
 737       enddo
 738       return
 739       end
 740 c-------------------------------------------------------------------------
 741       double precision function fdum()
 742       fdum=0.0D0
 743       return
 744       end