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