source/maxlik/src_FPy.org/proc.F

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