source/maxlik/src_FPy/proc.F

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