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