source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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 #ifdef DEBUG
 431             write (iout,*) "iref",iref," iii",iii," rms",rms," aux",aux
 432 #endif
 433             else
 434             rms=qwolyness(creff,iprot)
 435 #ifdef DEBUG
 436             write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
 437      &          e_total_T(iref,ib),
 438      &          e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
 439      &          " rms",-dlog(rms)," fac",rms,
 440      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
 441      &          rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 442 #endif
 443             aux=rms
 444             endif ! rmscomp
 445 c 7/6/16 AL Add angle comparison
 446             if (anglecomp(iprot).or.sidecomp(iprot))
 447      &        call int_from_cart(.true.,.false.)
 448             if (anglecomp(iprot)) then
 449               rmsthet = 0.0d0
 450               do k=nnt+2,nct
 451                 dtheta = theta(k)-thetareff(k)
 452                 rmsthet = rmsthet+dtheta*dtheta
 453               enddo
 454               rmsthet=rmsthet/(nct-nnt-2)
 455               rmsphi = 0.0d0
 456               do k=nnt+3,nct
 457                 dphi = pinorm(phi(k)-phireff(k))
 458                 rmsphi = rmsphi + dphi*dphi
 459               enddo
 460               rmsphi=rmsphi/(nct-nnt-3)
 461               aux=aux*dexp(-0.5d0*(rmsthet+rmsphi)/sigmaang2(iprot)**2)
 462 #ifdef DEBUG
 463               write (iout,*) "rmsthet",dsqrt(rmsthet)," rmsphi",
 464      &           dsqrt(rmsphi)," auxang",aux
 465 #endif
 466             endif
 467             if (sidecomp(iprot)) then
 468               call sc_loc_geom(.false.)
 469               rmsside = 0.0d0
 470               do k=nnt+1,nct-1
 471                 dxref = xxtab(k)-xxreff(k)
 472                 dyref = yytab(k)-yyreff(k)
 473                 dzref = zztab(k)-zzreff(k)
 474                 rmsside = rmsside + dx*dx+dy*dy+dz*dz
 475               enddo
 476               rmsside=rmsside/(nct-nnt+1)
 477               aux=aux*dexp(-0.5d0*rmsside/sigmaside2(iprot)**2)
 478             endif
 479           else
 480             aux=1.0d0
 481           endif
 482           do ib=1,nbeta(1,iprot)
 483             if (GEOM_AND_WHAM_WEIGHTS) then
 484               weirms(iii,ib) = weirms(iii,ib)
 485      &          +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 486             else
 487               weirms(iii,ib) = weirms(iii,ib)+aux
 488             endif
 489           enddo ! ib
 490         enddo ! iref
 491       enddo ! iii
 492 #ifdef DEBUG
 493       write (iout,*) "weirms"
 494       do iii=istart_conf,iend_conf
 495         write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1,
 496      &       nbeta(1,iprot))
 497       enddo
 498 #endif
 499       enddo ! i
 500 #endif
 501
 502 #ifdef MPI
 503       nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
 504       DO IB=1,NBETA(1,IPROT)
 505 #ifdef OUT_PTAB
 506       write(csa_bank,'(a,f4.0,4hPtab,i4.4)')
 507      &  protname(iprot)(:ilen(protname(iprot))),
 508      &  1.0d0/(0.001987*betaT(ib,1,iprot)),me
 509       open (icsa_bank,file=csa_bank,status="unknown")
 510 #endif
 511       ipass_conf=0
 512       jj=0
 513       sumP=0.0d0
 514       do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 515       ipass_conf=ipass_conf+1
 516       write (iout,*) "Pass",ipass_conf
 517       istart_conf=i
 518       iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
 519 #else
 520       nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
 521       ipass_conf=0
 522       sumP=0.0d0
 523       do i=1,ntot(iprot),maxstr_proc
 524       ipass_conf=ipass_conf+1
 525       write (iout,*) "Pass",ipass_conf
 526       istart_conf=i
 527       iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
 528 #endif
 529       ii=0
 530 c
 531 c Read the chunk of conformations off a DA scratchfile.
 532 c
 533       call daread_ccoords(iprot,istart_conf,iend_conf)
 534       write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
 535       do iii=istart_conf,iend_conf
 536 c        if (ib.eq.1) then
 537           rmsave=0.0d0
 538           rmsmin=1.0d10
 539 c        endif
 540         ii=ii+1
 541         jj=jj+1
 542         k=tsil(iii,iprot)
 543         call restore_ccoords(iprot,ii)
 544 c        write (iout,*) "zeroing Ptab",iii,ii,jj
 545         Ptab(jj,ib,iprot)=0.0d0
 546         if (rmscomp(iprot)) then
 547           do iref=1,nref(ib,iprot)
 548 #ifdef DEBUG
 549             itmp=ipdb
 550             ipdb=iout
 551             call pdbout(0.0d0,"conf")
 552 c            write (2,*) "nstart_sup",nstart_sup(iprot),
 553 c     &          " nend_sup",nend_sup(iprot)
 554 c            do k=nstart_sup(iprot),nend_sup(iprot)
 555 c              write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
 556 c     &          (cref(l,k,iref,ib,iprot),l=1,3)
 557 c            enddo
 558 c            do k=nstart_sup(iprot),nend_sup(iprot)
 559 c              write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3),
 560 c     &          (cref(l,k+nres,iref,ib,iprot),l=1,3)
 561 c            enddo
 562             do k=1,2*nres
 563               do l=1,3
 564                 cc(l,k)=c(l,k)
 565                 c(l,k)=cref(l,k,iref,ib,iprot)
 566               enddo
 567             enddo
 568             call pdbout(0.0d0,"ref")
 569             do k=1,2*nres
 570               do l=1,3
 571                 c(l,k)=cc(l,k)
 572               enddo
 573             enddo
 574             ipdb=itmp
 575 #endif
 576             rms=rmsnat(ib,jj,iref,iprot)
 577 c 7/6/16 AL Add angle comparison
 578             if (anglecomp(iprot).or.sidecomp(iprot))
 579      &        call int_from_cart(.true.,.false.)
 580             if (anglecomp(iprot)) then
 581 c              write (iout,*) "jj",jj," iref",iref
 582               rmsthet = 0.0d0
 583               do k=nnt+2,nct
 584                 dtheta = theta(k)-theta_ref(k,iref,ib,iprot)
 585 c                write (iout,*) k,theta(k),theta_ref(k,iref,ib,iprot),
 586 c     &            dtheta
 587                 rmsthet = rmsthet+dtheta*dtheta
 588               enddo
 589               rmsthet=rmsthet/(nct-nnt-2)
 590               rmsphi = 0.0d0
 591               do k=nnt+3,nct
 592                 dphi = pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
 593 c                write (iout,*) k,phi(k),phi_ref(k,iref,ib,iprot),
 594 c     &           pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
 595                 rmsphi = rmsphi + dphi*dphi
 596               enddo
 597               rmsphi=rmsphi/(nct-nnt-3)
 598             endif
 599             if (sidecomp(iprot)) then
 600               call sc_loc_geom(.false.)
 601               rmsside = 0.0d0
 602               do k=nnt+1,nct-1
 603                 dxref = xxtab(k)-xxref(k,iref,ib,iprot)
 604                 dyref = yytab(k)-yyref(k,iref,ib,iprot)
 605                 dzref = zztab(k)-zzref(k,iref,ib,iprot)
 606                 rmsside = rmsside + dx*dx+dy*dy+dz*dz
 607               enddo
 608               rmsside=rmsside/(nct-nnt+1)
 609             endif
 610 #ifdef DEBUG
 611             write (iout,*) "nnt",nnt," nct",nct," nthet",nct-nnt-2,
 612      &       " nphi",nct-nnt-3
 613             write (iout,*) "ii",ii," jj",jj," iref",iref,
 614      &     " rms",rms," rmsthet",dsqrt(rmsthet)," rmsphi",dsqrt(rmsphi),
 615      &     "efree",efreeref(iref,ib,iprot)
 616 #endif
 617 c            if (ib.eq.1) then
 618               rmsave=rmsave+rms
 619               if (rms.lt.rmsmin) then
 620                 imin=iref
 621                 rmsmin=rms
 622 c              endif
 623             endif
 624             Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)
 625      &          +dexp(-efreeref(iref,ib,iprot)
 626      &                -0.5d0*(rms/sigma2(iprot))**2)
 627             if (anglecomp(iprot))
 628      &      Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
 629      &                dexp(-0.5d0*(rmsthet+rmsphi)/sigmaang2(iprot)**2)
 630             if (sidecomp(iprot))
 631      &      Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
 632      &                dexp(-0.5d0*rmsside/sigmaside2(iprot)**2)
 633           enddo
 634 c          if (ib.eq.1) then
 635             rmsave=rmsave/nref(ib,iprot)
 636 c          endif
 637 #if defined(WEIENT)
 638 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 639 c              search of conformational space.
 640           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 641           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 642 #ifdef DEBUG
 643           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 644 c     &     " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot)
 645      &     "entfac",entfac(k,iprot)," weirms",weirms(iii,ib),
 646      &     " Ptab",Ptab(jj,ib,iprot)
 647 #endif
 648 #elif defined(WEIDIST)
 649 c          write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 650 c     &     " weirms",weirms(iii,ib)," Ptab/weirms",
 651 c     &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 652 #ifdef OUT_PTAB
 653           write (icsa_bank,'(2i6,3f10.5,2f8.4)')
 654      &      iii,imin,-dlog(Ptab(jj,ib,iprot)),
 655      &      -dlog(weirms(iii,ib)),
 656      &      -dlog(Ptab(jj,ib,iprot)/weirms(iii,ib)),
 657      &      rmsave,rmsmin
 658 #endif
 659           if (ib.eq.1) rmstb(iii,iprot)=rmsmin
 660           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 661 #elif defined(WEICLUST)
 662           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot)
 663 #endif
 664           sumP=sumP+Ptab(jj,ib,iprot)
 665         else
 666           do iref=1,nref(ib,iprot)
 667             rms = qwolynes(ib,iref,iprot)
 668 #ifdef DEBUG
 669             write (iout,*) "iii",iii," ii",ii," jj",jj,
 670      &        "iref",iref," rms",rms,-dlog(rms),
 671      &           "efree",efreeref(iref,ib,iprot)
 672 #endif
 673              Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms*
 674      &                        *dexp(-efreeref(iref,ib,iprot))
 675           enddo
 676 #if defined(WEIENT)
 677 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 678 c              search of conformational space.
 679           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 680           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 681 #ifdef DEBUG
 682           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 683      &    " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot)
 684 #endif
 685 #elif defined(WEICLUST)
 686           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot)
 687 #elif defined(WEIDIST)
 688           write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 689      &     " weirms",weirms(iii,ib)," Ptab/weirms",
 690      &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 691           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 692 #endif
 693           sumP=sumP+Ptab(jj,ib,iprot)
 694         endif
 695 #ifdef DEBUG
 696         write (iout,*) "ii",ii," jj",jj," ib",ib," iprot",iprot,
 697      &    " Ptab",Ptab(jj,ib,iprot)," sumP",sumP
 698 #endif
 699       enddo ! iii
 700       enddo ! i
 701 #ifdef DEBUG
 702       write (iout,*) "iprot",iprot," ib",ib," sumP before REDUCE",sumP
 703 #endif
 704 #ifdef MPI
 705       call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION,
 706      &        MPI_SUM,Comm1,IERROR)
 707       sumP=sumP_t
 708 #endif
 709 #ifdef DEBUG
 710       write (iout,*) "iprot",iprot," ib",ib," sumP after REDUCE",sumP
 711 #endif
 712       jj=0
 713       do iii=indstart(me1,iprot),indend(me1,iprot)
 714       jj=jj+1
 715 #ifdef DEBUG
 716       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 717 #endif
 718       Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP
 719 #ifdef DEBUG
 720       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 721 #endif
 722       enddo
 723 #ifdef OUT_PTAB
 724       close(icsa_bank)
 725 #endif
 726       ENDDO ! IB
 727       close(icbase)
 728       enddo ! iprot
 729 #ifdef DEBUG
 730       write (iout,*) "ELOWEST at the end of PROC_DATA"
 731       do iprot=1,nprot
 732         do ibatch=1,natlike(iprot)+2
 733           do ib=1,nbeta(ibatch,iprot)
 734             write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
 735      &         " elowest",elowest(ib,ibatch,iprot)
 736           enddo
 737         enddo
 738       enddo
 739       write (iout,*) "Ptab at the end of PROC"
 740       do iprot=1,nprot
 741         do ib=1,nbeta(1,iprot)
 742         write (iout,*) "Protein",iprot," beta",ib
 743         jj=0
 744         do i=indstart(me1,iprot),indend(me1,iprot)
 745           jj=jj+1
 746           write (iout,*) iprot,ib,ii,jj,Ptab(jj,ib,iprot)
 747         enddo
 748         enddo
 749       enddo
 750 #endif
 751       write (iout,*)
 752       return
 753    11 return1
 754       end
 755 c-------------------------------------------------------------------------------
 756       subroutine determine_histograms
 757       implicit none
 758       include "DIMENSIONS"
 759       include "DIMENSIONS.ZSCOPT"
 760       integer maxind
 761       parameter (maxind=1000)
 762 #ifdef MPI
 763       include "mpif.h"
 764       integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
 765       include "COMMON.MPI"
 766 #endif
 767       include "COMMON.COMPAR"
 768       include "COMMON.ENERGIES"
 769       include "COMMON.VMCPAR"
 770       include "COMMON.WEIGHTS"
 771       include "COMMON.PROTNAME"
 772       include "COMMON.IOUNITS"
 773       include "COMMON.FFIELD"
 774       include "COMMON.PROTFILES"
 775       include "COMMON.CHAIN"
 776       include "COMMON.CONTROL"
 777       include "COMMON.CLASSES"
 778       integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt
 779       double precision eminh,emaxh,tminh,tmaxh
 780       integer he(0:maxind),
 781      &  ht(0:maxind),hemax,htmax
 782 c
 783 c Construct energy and entropy histograms
 784 c
 785       write (iout,*) "Calling DETERMINE HISTOGRAMS"
 786       call flush(iout)
 787       do iprot=1,nprot
 788           eminh=1.0d20
 789           emaxh=-1.0d20
 790           tminh=1.0d20
 791           tmaxh=-1.0d20
 792           do i=0,maxind
 793             he(i)=0
 794             ht(i)=0
 795           enddo
 796           do i=1,ntot(iprot)
 797             iii=i
 798             if (eini(iii,iprot).lt.eminh)
 799      &        eminh=eini(iii,iprot)
 800             if (eini(iii,iprot).gt.emaxh)
 801      &        emaxh=eini(iii,iprot)
 802              if (entfac(iii,iprot).lt.tminh)
 803      &        tminh=entfac(iii,iprot)
 804             if (entfac(iii,iprot).gt.tmaxh)
 805      &        tmaxh=entfac(iii,iprot)
 806           enddo
 807 c          write (2,*)  "DETERMINE HISTOGRAMS 1"
 808           call flush(iout)
 809           do i=1,ntot(iprot)
 810             iii=i
 811             inde=eini(iii,iprot)-eminh
 812             indt=entfac(iii,iprot)-tminh
 813             if (inde.le.maxind) he(inde)=he(inde)+1
 814             if (indt.le.maxind) ht(indt)=ht(indt)+1
 815           enddo
 816 c          write (2,*)  "DETERMINE HISTOGRAMS 2"
 817             idelte=emaxh-eminh
 818             ideltt=tmaxh-tminh
 819             hemax=0
 820 c            write (iout,*) "idelte",idelte," ideltt",ideltt
 821             call flush(iout)
 822             do i=0,idelte
 823 c              write (iout,*) "i",i," he",he(i)," hemax",hemax
 824               call flush(iout)
 825               if (he(i).gt.hemax) hemax=he(i)
 826             enddo
 827             htmax=0
 828             do i=0,ideltt
 829 c              write (iout,*) "i",i," ht",ht(i)," htmax",htmax
 830               call flush(iout)
 831               if (ht(i).gt.htmax) htmax=ht(i)
 832             enddo
 833 c          write (2,*)  "DETERMINE HISTOGRAMS 3"
 834           call flush(iout)
 835             hemax=hemax/hefac(iprot)
 836             if (hemax.lt.hemax_low(iprot))
 837      &        hemax=hemax_low(iprot)
 838             htmax=htmax/htfac(iprot)
 839             if (htmax.lt.htmax_low(iprot))
 840      &        htmax=htmax_low(iprot)
 841             i=0
 842             do while (i.lt.idelte .and. he(i).lt.hemax)
 843               i=i+1
 844             enddo
 845 c          write (2,*)  "DETERMINE HISTOGRAMS 4"
 846           call flush(iout)
 847             e_lowb(iprot)=eminh+i
 848             i=0
 849             do while (i.lt.ideltt .and. ht(i).lt.htmax)
 850               i=i+1
 851             enddo
 852             t_lowb(iprot)=tminh+i
 853 #ifdef DEBUG
 854             write (iout,*) "protein",iprot
 855             write (iout,*) "energy histogram"
 856             do i=0,idelte
 857               write (iout,'(f15.5,i10)') eminh+i,he(i)
 858             enddo
 859             write (iout,*) "entropy histogram"
 860             do i=0,ideltt
 861               write (iout,'(f15.5,i10)') tminh+i,ht(i)
 862             enddo
 863             write (iout,*) "emin",eminh," tmin",tminh,
 864      &       " hemax",hemax," htmax",htmax,
 865      &       " e_lowb",e_lowb(iprot),
 866      &       " t_lowb",t_lowb(iprot)
 867             call flush(iout)
 868 #endif
 869       enddo
 870       return
 871       end
 872 c------------------------------------------------------------------------------
 873       subroutine imysort(n, x, ipermut)
 874       implicit none
 875       integer n
 876       integer x(n),xtemp
 877       integer ipermut(n)
 878       integer i,j,imax,ipm
 879       do i=1,n
 880         xtemp=x(i)
 881         imax=i
 882         do j=i+1,n
 883           if (x(j).lt.xtemp) then
 884             imax=j
 885             xtemp=x(j)
 886           endif
 887         enddo
 888         x(imax)=x(i)
 889         x(i)=xtemp
 890         ipm=ipermut(imax)
 891         ipermut(imax)=ipermut(i)
 892         ipermut(i)=ipm
 893       enddo
 894       return
 895       end
 896 c--------------------------------------------------------------------------
 897       subroutine write_conf_count
 898       implicit none
 899       include "DIMENSIONS"
 900       include "DIMENSIONS.ZSCOPT"
 901       include "COMMON.CLASSES"
 902       include "COMMON.COMPAR"
 903       include "COMMON.VMCPAR"
 904       include "COMMON.PROTNAME"
 905       include "COMMON.IOUNITS"
 906       include "COMMON.PROTFILES"
 907       integer iprot,ibatch,i,ii,j,kk
 908       integer ilen
 909       external ilen
 910       logical LPRN
 911       write (iout,*)
 912       write (iout,*)"Numbers of conformations after applying the cutoff"
 913       write (iout,*)
 914       do iprot=1,nprot
 915         write (iout,'(a,2x,a,i10)') "Protein",
 916      &    protname(iprot)(:ilen(protname(iprot))),
 917      &    ntot_work(iprot)
 918       enddo
 919       return
 920       end
 921 c-------------------------------------------------------------------------
 922       double precision function fdum()
 923       fdum=0.0D0
 924       return
 925       end