source/maxlik/src-Fmatch/ptab_calc.F

   1       subroutine ptab_calc
   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       real csingle(3,maxres2)
  32       double precision entfacmax,sument,sumentsq,aux
  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       integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng,
  38      &  inn,in,iref,ib,ibatch,num,itmp,imin,iperm
  39       integer jj,istart_conf,iend_conf,ipass_conf,iii
  40       integer ilen,iroof
  41       external ilen,iroof
  42       double precision ej,emaxj,sumP,sumP_t
  43       double precision fac
  44       double precision dx,dy,dz,dxref,dyref,dzref,
  45      &  dtheta,dphi,rmsthet,rmsphi,rmsside
  46 #ifdef MPI
  47       double precision e_total_T_(maxstr_proc)
  48 #endif
  49       double precision rms,rmsave,rmsmin,rmsnat,rmscalc,rmscalc_thet,
  50      &   rmscalc_phi,rmscalc_side,qwolynes,qwolyness
  51
  52       do iprot=1,nprot
  53
  54       call restore_molinfo(iprot)
  55
  56       open (icbase,file=bprotfiles(iprot),status="old",
  57      &    form="unformatted",access="direct",recl=lenrec(iprot))
  58       entfacmax=entfac(1,iprot)
  59       sument=entfac(1,iprot)
  60       sumentsq=entfac(1,iprot)**2
  61       do i=2,ntot_work(iprot)
  62         if (entfac(i,iprot).gt.entfacmax) entfacmax=entfac(i,iprot)
  63         sument=sument+entfac(i,iprot)
  64         sumentsq=sumentsq+entfac(i,iprot)**2
  65       enddo
  66       sument=sument/ntot_work(iprot)
  67       sumentsq=dsqrt(sumentsq/ntot_work(iprot)-sument**2)
  68       write (2,*) "entfacmax",entfacmax," entave",sument,
  69      &   " stdent",sumentsq
  70 #ifdef WEIDIST
  71 #ifdef MPI
  72       nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
  73       ipass_conf=0
  74       jj=0
  75       sumP=0.0d0
  76       do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
  77       ipass_conf=ipass_conf+1
  78       write (iout,*) "Pass",ipass_conf
  79       istart_conf=i
  80       iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
  81 #else
  82       nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
  83       ipass_conf=0
  84       sumP=0.0d0
  85       do i=1,ntot(iprot),maxstr_proc
  86       ipass_conf=ipass_conf+1
  87       write (iout,*) "Pass",ipass_conf
  88       istart_conf=i
  89       iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
  90 #endif
  91 c
  92 c Read the chunk of conformations off a DA scratchfile.
  93 c
  94       call daread_ccoords(iprot,istart_conf,iend_conf)
  95
  96       IF (GEOM_AND_WHAM_WEIGHTS) THEN
  97
  98 c Compute energies
  99       do ib=1,nbeta(1,iprot)
 100         write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
 101         ii=0
 102         do iii=istart_conf,iend_conf
 103           ii=ii+1
 104           kk=tsil(iii,iprot)
 105           if (kk.eq.0) cycle
 106           etoti=0.0d0
 107           do k=1,n_ene
 108             etoti=etoti+ww(k)*escal(k,ib,1,iprot)
 109      &            *enetb(ii,k,iprot)
 110           enddo
 111           e_total_T(iii,ib)=entfac(kk,iprot)
 112      &             +betaT(ib,1,iprot)*(etoti-elowest(ib,1,iprot))
 113         enddo
 114       enddo ! ib
 115 #ifdef DEBUG
 116       write (iout,*) "Energies before Gather"
 117       do iii=1,ntot_work(iprot)
 118         write (iout,'(i5,100E12.5)') iii,
 119      &         (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
 120       enddo
 121 #endif
 122 #ifdef MPI
 123       do ib=1,nbeta(1,iprot)
 124         do k=1,scount(me1,iprot)
 125           e_total_T_(k)=e_total_T(indstart(me1,iprot)+k-1,ib)
 126         enddo
 127 c        call MPI_AllGatherv(e_total_T(indstart(me1,iprot),ib),
 128 c     &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
 129 c     &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
 130 c     &     Comm1, IERROR)
 131         call MPI_AllGatherv(e_total_T_(1),
 132      &     scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
 133      &     scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
 134      &     Comm1, IERROR)
 135       enddo ! ib
 136 #ifdef DEBUG
 137       write (iout,*) "Energies after Gather"
 138       do iii=1,ntot_work(iprot)
 139         write (iout,'(i5,100E12.5)') iii,
 140      &         (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
 141       enddo
 142 #endif
 143 #endif
 144
 145       ENDIF ! GEOM_AND_WHAM_WEIGHTS
 146       ii=0
 147 c Compute distances
 148 c      write (iout,*) "Start calculating weirms",istart_conf,iend_conf
 149       do iii=istart_conf,iend_conf
 150         ii=ii+1
 151         k=tsil(iii,iprot)
 152 c        write (iout,*) "Calculating weirms iii",iii," ii",ii," k",k
 153         if (k.eq.0) cycle
 154         call restore_ccoords(iprot,ii)
 155 c Calculate the rmsds of the current conformations and convert them into weights
 156 c to approximate equal weighting of all points of the conformational space
 157         do k=1,2*nres
 158           do l=1,3
 159             creff(l,k)=c(l,k)
 160           enddo
 161         enddo
 162 c 7/6/16 AL Add angle comparison
 163         if (anglecomp(iprot) .or. sidecomp(iprot)) then
 164           call int_from_cart(.true.,.false.)
 165           call sc_loc_geom(.false.)
 166           do k=1,nres
 167             phireff(k)=phi(k)
 168             thetareff(k)=theta(k)
 169             xxreff(k)=xxtab(k)
 170             yyreff(k)=yytab(k)
 171             zzreff(k)=zztab(k)
 172           enddo
 173         endif
 174         do ib=1,nbeta(1,iprot)
 175          weirms(iii,ib)=0.0d0
 176         enddo
 177 C-----------TEST
 178 c            caonly(iprot)=.true.
 179 C-----------TEST
 180         do iref=1,ntot_work(iprot)
 181           if (iref.ne.iii) then
 182             read(icbase,rec=iref) ((csingle(l,k),l=1,3),k=1,nres),
 183      &      ((csingle(l,k),l=1,3),k=nnt+nres,nct+nres),
 184      &      nss,(ihpb(k),jhpb(k),k=1,nss),eini(iref,iprot),
 185      &      entfac(iref,iprot),
 186      &      ((nu(k,l,ii,iprot),k=1,maxdimnat),l=1,natlike(iprot))
 187             do k=1,2*nres
 188               do l=1,3
 189                 c(l,k)=csingle(l,k)
 190               enddo
 191             enddo
 192 #ifdef DEBUG
 193             do k=1,nres
 194               write (iout,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
 195      &                      (creff(l,k),l=1,3)
 196             enddo
 197 #endif
 198             rms = 0.0d0
 199             iperm=1
 200             if (rmscomp(iprot)) then
 201               rms=rmscalc(c(1,1),creff(1,1),iref,iii,iprot,iperm)
 202 #ifdef DEBUG
 203               write (iout,*) "iref",iref," iii",iii," rms",rms
 204 #endif
 205 #ifdef DEBUG
 206               write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
 207      &          e_total_T(iref,ib),
 208      &          e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
 209      &          " rms",rms," fac",dexp(-0.5d0*(rms/sigma2(iprot))**2),
 210      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
 211      &          dexp(-0.5d0*(rms/sigma2(iprot))**2)*
 212      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 213 #endif
 214               aux=dexp(-0.5d0*(rms/sigma2(iprot))**2)
 215 #ifdef DEBUG
 216               write (iout,*)"iref",iref," iii",iii," rms",rms," aux",aux
 217 #endif
 218             else
 219               rms=qwolyness(creff,iprot)
 220 #ifdef DEBUG
 221               write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
 222      &          e_total_T(iref,ib),
 223      &          e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
 224      &          " rms",-dlog(rms)," fac",rms,
 225      &          dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
 226      &          rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 227 #endif
 228               aux=rms
 229             endif ! rmscomp
 230 c 7/6/16 AL Add angle comparison
 231             if (anglecomp(iprot).or.sidecomp(iprot))
 232      &        call int_from_cart(.true.,.false.)
 233             if (anglecomp(iprot)) then
 234 c              write (iout,*) "jj",jj," iref",iref
 235               rmsthet=rmscalc_thet(theta(1),thetareff(1),iperm)
 236               rmsphi=rmscalc_phi(phi(1),phireff(1),iperm)
 237 #ifdef DEBUG
 238               write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi
 239 #endif
 240              aux=aux*dexp(-0.5d0*(rmsthet**2+rmsphi**2)/
 241      &         sigmaang2(iprot)**2)
 242 #ifdef DEBUG
 243              write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi,
 244      &        " auxang",aux
 245 #endif
 246             endif
 247             if (sidecomp(iprot)) then
 248               call sc_loc_geom(.false.)
 249               rmsside=rmscalc_side(xxtab(1),yytab(1),zztab(1),
 250      &          xxreff(1),yyreff(1),zzreff(1),iperm)
 251 #ifdef DEBUG
 252               write (iout,*) "rmsside",rmsside
 253 #endif
 254               aux=aux*dexp(-0.5d0*(rmsside/sigmaside2(iprot))**2)
 255             endif
 256           else
 257             aux=1.0d0
 258           endif
 259           do ib=1,nbeta(1,iprot)
 260             if (GEOM_AND_WHAM_WEIGHTS) then
 261               weirms(iii,ib) = weirms(iii,ib)
 262      &          +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
 263             else
 264               weirms(iii,ib) = weirms(iii,ib)+aux
 265             endif
 266           enddo ! ib
 267         enddo ! iref
 268       enddo ! iii
 269 C-----------TEST
 270 c            caonly(iprot)=.false.
 271 C-----------TEST
 272 #ifdef DEBUG
 273       write (iout,*) "weirms"
 274       do iii=istart_conf,iend_conf
 275         write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1,
 276      &       nbeta(1,iprot))
 277       enddo
 278 #endif
 279       enddo ! i
 280 #endif
 281
 282 #ifdef MPI
 283       nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
 284       DO IB=1,NBETA(1,IPROT)
 285 #ifdef OUT_PTAB
 286       write(csa_bank,'(a,f4.0,4hPtab,i4.4)')
 287      &  protname(iprot)(:ilen(protname(iprot))),
 288      &  1.0d0/(0.001987*betaT(ib,1,iprot)),me
 289       open (icsa_bank,file=csa_bank,status="unknown")
 290 #endif
 291       ipass_conf=0
 292       jj=0
 293       sumP=0.0d0
 294       do i=indstart(me1,iprot),indend(me1,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,indend(me1,iprot))
 299 #else
 300       nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
 301       ipass_conf=0
 302       sumP=0.0d0
 303       do i=1,ntot(iprot),maxstr_proc
 304       ipass_conf=ipass_conf+1
 305       write (iout,*) "Pass",ipass_conf
 306       istart_conf=i
 307       iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
 308 #endif
 309       ii=0
 310 c
 311 c Read the chunk of conformations off a DA scratchfile.
 312 c
 313       call daread_ccoords(iprot,istart_conf,iend_conf)
 314       write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
 315       do iii=istart_conf,iend_conf
 316 c        if (ib.eq.1) then
 317           rmsave=0.0d0
 318           rmsmin=1.0d10
 319 c        endif
 320         ii=ii+1
 321         jj=jj+1
 322         k=tsil(iii,iprot)
 323         call restore_ccoords(iprot,ii)
 324 c        write (iout,*) "zeroing Ptab",iii,ii,jj
 325         Ptab(jj,ib,iprot)=0.0d0
 326         if (rmscomp(iprot)) then
 327           do iref=1,nref(ib,iprot)
 328 #ifdef DEBUG
 329             itmp=ipdb
 330             ipdb=iout
 331             call pdbout(0.0d0,"conf")
 332 c            write (2,*) "nstart_sup",nstart_sup(iprot),
 333 c     &          " nend_sup",nend_sup(iprot)
 334 c            do k=nstart_sup(iprot),nend_sup(iprot)
 335 c              write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
 336 c     &          (cref(l,k,iref,ib,iprot),l=1,3)
 337 c            enddo
 338 c            do k=nstart_sup(iprot),nend_sup(iprot)
 339 c              write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3),
 340 c     &          (cref(l,k+nres,iref,ib,iprot),l=1,3)
 341 c            enddo
 342             do k=1,2*nres
 343               do l=1,3
 344                 cc(l,k)=c(l,k)
 345                 c(l,k)=cref(l,k,iref,ib,iprot)
 346               enddo
 347             enddo
 348             call pdbout(0.0d0,"ref")
 349             do k=1,2*nres
 350               do l=1,3
 351                 c(l,k)=cc(l,k)
 352               enddo
 353             enddo
 354             ipdb=itmp
 355 #endif
 356             rms=rmsnat(ib,jj,iref,iprot,iperm)
 357 #ifdef DEBUG
 358             write (iout,*) iii,iref," rmsnat",rms
 359 #endif
 360 c 7/6/16 AL Add angle comparison
 361 c            write(iout,*) "anglecomp",anglecomp
 362             if (anglecomp(iprot).or.sidecomp(iprot))
 363      &        call int_from_cart(.true.,.false.)
 364             if (anglecomp(iprot)) then
 365 c              write (iout,*) "jj",jj," iref",iref
 366               rmsthet=rmscalc_thet(theta(1),
 367      &          theta_ref(1,iref,ib,iprot),iperm)
 368               rmsphi=rmscalc_phi(phi(1),phi_ref(1,iref,ib,iprot),
 369      &          iperm)
 370 #ifdef DEBUG
 371               write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi
 372 #endif
 373             endif
 374             if (sidecomp(iprot)) then
 375               call sc_loc_geom(.false.)
 376               rmsside=rmscalc_side(xxtab(1),yytab(1),zztab(1),
 377      &          xxref(1,iref,ib,iprot),yyref(1,iref,ib,iprot),
 378      &          zzref(1,iref,ib,iprot),iperm)
 379 #ifdef DEBUG
 380               write (iout,*) "rmsside",rmsside
 381 #endif
 382             endif
 383 #ifdef DEBUG
 384             write (iout,*) "ii",ii," jj",jj," iref",iref,
 385      &            " rms",rms," rmsthet",rmsthet," rmsphi",rmsphi,
 386      &            " efree",efreeref(iref,ib,iprot)
 387 #endif
 388 c            if (ib.eq.1) then
 389               rmsave=rmsave+rms
 390               if (rms.lt.rmsmin) then
 391                 imin=iref
 392                 rmsmin=rms
 393 c              endif
 394             endif
 395             Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)
 396      &          +dexp(-efreeref(iref,ib,iprot)
 397      &                -0.5d0*(rms/sigma2(iprot))**2)
 398             if (anglecomp(iprot))
 399      &      Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
 400      &         dexp(-0.5d0*(rmsthet**2+rmsphi**2)/sigmaang2(iprot)**2)
 401             if (sidecomp(iprot))
 402      &      Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
 403      &              dexp(-0.5d0*(rmsside/sigmaside2(iprot))**2)
 404           enddo
 405 c          if (ib.eq.1) then
 406             rmsave=rmsave/nref(ib,iprot)
 407 c          endif
 408 #if defined(WEIENT)
 409 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 410 c              search of conformational space.
 411           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 412           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 413 #ifdef DEBUG
 414           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 415 c     &     " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot)
 416      &     "entfac",entfac(k,iprot)," weirms",weirms(iii,ib),
 417      &     " Ptab",Ptab(jj,ib,iprot)
 418 #endif
 419 #elif defined(WEIDIST)
 420 c          write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 421 c     &     " weirms",weirms(iii,ib)," Ptab/weirms",
 422 c     &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 423 #ifdef OUT_PTAB
 424           write (icsa_bank,'(2i6,3f10.5,2f8.4)')
 425      &      iii,imin,-dlog(Ptab(jj,ib,iprot)),
 426      &      -dlog(weirms(iii,ib)),
 427      &      -dlog(Ptab(jj,ib,iprot)/weirms(iii,ib)),
 428      &      rmsave,rmsmin
 429 #endif
 430           if (ib.eq.1) rmstb(iii,iprot)=rmsmin
 431           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 432 #elif defined(WEICLUST)
 433           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot)
 434 #endif
 435           sumP=sumP+Ptab(jj,ib,iprot)
 436         else
 437           do iref=1,nref(ib,iprot)
 438             rms = qwolynes(ib,iref,iprot)
 439 #ifdef DEBUG
 440             write (iout,*) "iii",iii," ii",ii," jj",jj,
 441      &        "iref",iref," rms",rms,-dlog(rms),
 442      &           "efree",efreeref(iref,ib,iprot)
 443 #endif
 444              Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms*
 445      &                        *dexp(-efreeref(iref,ib,iprot))
 446           enddo
 447 #if defined(WEIENT)
 448 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
 449 c              search of conformational space.
 450           aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
 451           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
 452 #ifdef DEBUG
 453           write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
 454      &    " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot)
 455 #endif
 456 #elif defined(WEICLUST)
 457           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot)
 458 #elif defined(WEIDIST)
 459           write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
 460      &     " weirms",weirms(iii,ib)," Ptab/weirms",
 461      &     Ptab(jj,ib,iprot)/weirms(iii,ib)
 462           Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
 463 #endif
 464           sumP=sumP+Ptab(jj,ib,iprot)
 465         endif
 466 #ifdef DEBUG
 467         write (iout,*) "ii",ii," jj",jj," ib",ib," iprot",iprot,
 468      &    " Ptab",Ptab(jj,ib,iprot)," sumP",sumP
 469 #endif
 470       enddo ! iii
 471       enddo ! i
 472 #ifdef DEBUG
 473       write (iout,*) "iprot",iprot," ib",ib," sumP before REDUCE",sumP
 474 #endif
 475 #ifdef MPI
 476       call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION,
 477      &        MPI_SUM,Comm1,IERROR)
 478       sumP=sumP_t
 479 #endif
 480 #ifdef DEBUG
 481       write (iout,*) "iprot",iprot," ib",ib," sumP after REDUCE",sumP
 482 #endif
 483       jj=0
 484       do iii=indstart(me1,iprot),indend(me1,iprot)
 485       jj=jj+1
 486 #ifdef DEBUG
 487       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 488 #endif
 489       Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP
 490 #ifdef DEBUG
 491       write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
 492 #endif
 493       enddo
 494 #ifdef OUT_PTAB
 495       close(icsa_bank)
 496 #endif
 497       ENDDO ! IB
 498       close(icbase)
 499       enddo ! iprot
 500 #ifdef DEBUG
 501       write (iout,*) "ELOWEST at the end of PROC_DATA"
 502       do iprot=1,nprot
 503         do ibatch=1,natlike(iprot)+2
 504           do ib=1,nbeta(ibatch,iprot)
 505             write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
 506      &         " elowest",elowest(ib,ibatch,iprot)
 507           enddo
 508         enddo
 509       enddo
 510       write (iout,*) "Ptab at the end of PROC"
 511       do iprot=1,nprot
 512         do ib=1,nbeta(1,iprot)
 513         write (iout,*) "Protein",iprot," beta",ib
 514         jj=0
 515         do i=indstart(me1,iprot),indend(me1,iprot)
 516           jj=jj+1
 517           write (iout,*) iprot,ib,ii,jj,Ptab(jj,ib,iprot)
 518         enddo
 519         enddo
 520       enddo
 521 #endif
 522       return
 523       end