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