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