source/maxlik/src_FPy.org/averages_sc.F

   1       subroutine averages(iprot)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5 #ifdef MPI
   6       include "mpif.h"
   7       integer IERROR,ErrCode
   8       include "COMMON.MPI"
   9 #endif
  10       integer n,nf,What
  11       include "COMMON.WEIGHTS"
  12       include "COMMON.WEIGHTDER"
  13       include "COMMON.ENERGIES"
  14       include "COMMON.IOUNITS"
  15       include "COMMON.VMCPAR"
  16       include "COMMON.NAMES"
  17       include "COMMON.INTERACT"
  18       include "COMMON.TIME1"
  19       include "COMMON.CHAIN"
  20       include "COMMON.PROTFILES"
  21       include "COMMON.COMPAR"
  22       include "COMMON.CLASSES"
  23 C Define local variables
  24       integer i,ii,iii,inat,jj,j,k,kk,l,ik,jk,iprot,ib
  25       integer ipass_conf,istart_conf,iend_conf
  26       double precision energia(0:max_ene)
  27       double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk
  28       double precision aux,fac,ef1,ef2,em1,em2,var1,var2
  29       double precision gnmr,tcpu_ini,tcpu
  30       logical lprn
  31       integer icant
  32       external icant
  33       lprn=.true.
  34 c
  35 c Zero out tables.
  36 c
  37 c Mazimum likeliood
  38       do ib=1,nbeta(1,iprot)
  39         do k=1,n_ene
  40           sumlikder(k,ib,iprot)=0.0d0
  41           sumqder(k,ib,iprot)=0.0d0
  42         enddo
  43         do k=1,nntyp
  44           sumlikeps(k,ib,iprot)=0.0d0
  45           sumqeps(k,ib,iprot)=0.0d0
  46         enddo
  47         do k=1,maxtor_var
  48           sumliktor(k,ib,iprot)=0.0d0
  49           sumqtor(k,ib,iprot)=0.0d0
  50         enddo
  51         sumlik(ib,iprot)=0.0d0
  52         efree(ib,1,iprot)=0.0d0
  53       enddo
  54 c Heat capacity
  55       do ib=1,nbeta(2,iprot)
  56         do k=1,n_ene
  57           eave_pftot(k,ib,iprot)=0.0d0
  58           eave_pfprimtot(k,ib,iprot)=0.0d0
  59           eave_pfbistot(k,ib,iprot)=0.0d0
  60           emix_pftot(k,ib,iprot)=0.0d0
  61           emix_pfprimtot(k,ib,iprot)=0.0d0
  62           emix_pfbistot(k,ib,iprot)=0.0d0
  63           emixsq_pftot(k,ib,iprot)=0.0d0
  64         enddo
  65         do k=1,nntyp
  66           enepsave_ftot(k,ib,iprot)=0.0d0
  67           eneps_mix_ftot(k,ib,iprot)=0.0d0
  68           eneps_mix_fbistot(k,ib,iprot)=0.0d0
  69           eneps_mixsq_ftot(k,ib,iprot)=0.0d0
  70         enddo
  71         do k=1,maxtor_var
  72           entorave_ftot(k,ib,iprot)=0.0d0
  73           entorave_fprimtot(k,ib,iprot)=0.0d0
  74           entorave_fbistot(k,ib,iprot)=0.0d0
  75           entor_mix_ftot(k,ib,iprot)=0.0d0
  76           entor_mix_fprimtot(k,ib,iprot)=0.0d0
  77           entor_mix_fbistot(k,ib,iprot)=0.0d0
  78           entor_mixsq_ftot(k,ib,iprot)=0.0d0
  79         enddo
  80         emean_ftot(ib,iprot)=0.0d0
  81         ebis_ftot(ib,iprot)=0.0d0
  82         esquare_ftot(ib,iprot)=0.0d0
  83       enddo
  84 c Conformation-dependent averages
  85       do inat=1,natlike(iprot)
  86         do ib=1,nbeta(inat+2,iprot)
  87           do l=1,natdim(inat,iprot)
  88             do k=1,n_ene
  89               nu_pf(k,l,ib,inat,iprot)=0.0d0
  90             enddo
  91             do k=1,nntyp
  92               nuepsave_f(k,l,ib,inat,iprot)=0.0d0
  93             enddo
  94             do k=1,maxtor_var
  95               nutorave_f(k,l,ib,inat,iprot)=0.0d0
  96             enddo
  97             nuave(l,ib,inat,iprot)=0.0d0
  98           enddo
  99         enddo
 100       enddo
 101 c
 102 c Calculate the contributions to averages from each processor or all
 103 c contributions if calculations are run in uniprocessor mode.
 104 c The derivatives of energy in epsilons are dumped to disk, if needed.
 105 c
 106 #ifdef DEBUG
 107       write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot)
 108 #endif
 109       IF (NCHUNK_CONF(IPROT).EQ.1) THEN
 110
 111 #ifdef MPI
 112       do i=1,ntot_work(iprot)
 113         i2ii(i,iprot)=0
 114       enddo
 115       ii=0
 116       do i=indstart(me1,iprot),indend(me1,iprot)
 117         ii=ii+1
 118         i2ii(i,iprot)=ii
 119       enddo
 120 #else
 121       do i=1,ntot_work(iprot)
 122         i2ii(i,iprot)=i
 123       endif
 124 #endif
 125 #ifdef DEBUG
 126       do i=1,ntot_work(iprot)
 127         write (iout,*) "i",i," i2ii",i2ii(i,iprot)
 128       enddo
 129       call flush(iout)
 130 #endif
 131       call ave_eval(iprot)
 132
 133       ELSE
 134
 135       open (ientin,file=benefiles(iprot),status="old",
 136      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 137       ipass_conf=0
 138 #ifdef MPI
 139       do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 140         iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
 141 #else
 142       do istart_conf=1,ntot_work(iprot),maxstr_proc
 143         iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
 144 #endif
 145 c
 146 c Read the chunk of energies and derivatives off a DA scratchfile.
 147 c
 148         ipass_conf=ipass_conf+1
 149         do i=1,ntot_work(iprot)
 150           i2ii(i,iprot)=0
 151         enddo
 152         ii=0
 153         do i=istart_conf,iend_conf
 154           ii=ii+1
 155           i2ii(i,iprot)=ii
 156         enddo
 157 #ifdef DEBUG
 158         write (iout,*) "ipass_conf",ipass_conf,
 159      &    " istart_conf",istart_conf," iend_conf",iend_conf
 160         do i=1,ntot_work(iprot)
 161           write (iout,*) "i",i," i2ii",i2ii(i,iprot)
 162         enddo
 163         call flush(iout)
 164 #endif
 165         call daread_ene(iprot,istart_conf,iend_conf)
 166         call ave_eval(iprot)
 167       enddo
 168
 169       close(ientin)
 170       ENDIF
 171 c
 172 c Put tohether all contributions.
 173 c
 174       call ave_sum(iprot)
 175       return
 176       end
 177 c-------------------------------------------------------------------------------
 178       subroutine ave_eval(iprot)
 179       implicit none
 180       include "DIMENSIONS"
 181       include "DIMENSIONS.ZSCOPT"
 182 #ifdef MPI
 183       include "mpif.h"
 184       integer IERROR,ErrCode
 185       include "COMMON.MPI"
 186 #endif
 187       integer n,nf,What
 188       include "COMMON.WEIGHTS"
 189       include "COMMON.WEIGHTDER"
 190       include "COMMON.ENERGIES"
 191       include "COMMON.IOUNITS"
 192       include "COMMON.VMCPAR"
 193       include "COMMON.NAMES"
 194       include "COMMON.INTERACT"
 195       include "COMMON.TIME1"
 196       include "COMMON.CHAIN"
 197       include "COMMON.PROTFILES"
 198       include "COMMON.THERMAL"
 199       include "COMMON.VAR"
 200       include "COMMON.CLASSES"
 201       include "COMMON.COMPAR"
 202       include "COMMON.TORSION"
 203       include "COMMON.SCCOR"
 204       include "COMMON.PROTNAME"
 205 C Define local variables
 206       integer i,ii,iii,jj,j,k,kk,l,ll,m,ik,jk,iprot,ib,inat,
 207      &  ic,iblock
 208       integer ilen
 209       external ilen
 210       integer ipass_conf,istart_conf,iend_conf
 211       double precision energia(0:max_ene)
 212       double precision etoti,sumpart,esquarei,emeani,
 213      &  elowesti,enepsjk,eprim,ebis,eprimi,ebisi,etoti_orig,
 214      &  entorjk,eave_pft(max_ene),emix_pft(max_ene),
 215      &  esquare_ft,efree_t,emixsq_pft(max_ene),
 216      &  eneps_mixt(nntyp),eneps_meant(nntyp),
 217      &  enepsave_ft(nntyp),eneps_mix_ft(nntyp),
 218      &  eneps_mixsq_ft(nntyp),emean_ft
 219       double precision aux,auxf,fac,facf,ef1,ef2,em1,em2,var1,var2,
 220      &  deltei,deltei_orig,temper,elowest_all
 221       double precision gnmr,tcpu_ini,tcpu
 222       double precision ftune_epsprim
 223       external ftune_epsprim
 224       logical lprn
 225       integer icant
 226       external icant
 227
 228       lprn=.false.
 229 C Compute the likelihood sum
 230       DO IB=1,NBETA(1,IPROT)
 231         elowest_all=elowest(ib,1,iprot)
 232 #ifdef OUT_LIK
 233       write(csa_bank,'(a,f4.0,4hMlik,i3.3)')
 234      &  protname(iprot)(:ilen(protname(iprot))),
 235      &  1.0d0/(0.001987*betaT(ib,1,iprot)),me
 236       open (icsa_bank,file=csa_bank,status="unknown")
 237 #endif
 238 #ifdef DEBUG
 239           write(iout,*) "iprot",iprot," ib",ib,
 240      &      " elowest",elowest_l(ib,iprot)
 241 #endif
 242         fac=betaT(ib,1,iprot)
 243         temper=1.0d0/(fac*Rgas)
 244         sumpart=0.0d0
 245         sumlik(ib,iprot)=0.0d0
 246         do k=1,n_ene
 247           sumlikder(k,ib,iprot)=0.0d0
 248           sumqder(k,ib,iprot)=0.0d0
 249         enddo
 250         do j=1,ntot_work(iprot)
 251           i = tsil(j,iprot)
 252 #ifdef DEBUG
 253           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 254      &      " indstart",indstart(me1,iprot),
 255      &      " indend",indend(me1,iprot)
 256 #endif
 257           if (i.eq.0) goto 221
 258           jj = i2ii(i,iprot)
 259           if (jj.gt.0) then
 260 c Temperature-dependent energy
 261             etoti=0.0d0
 262             do k=1,n_ene
 263               etoti=etoti+ww(k)*escal(k,ib,1,iprot)
 264      &            *enetb(jj,k,iprot)
 265             enddo
 266             deltei=etoti-elowest(ib,1,iprot)
 267 #ifdef DEBUG
 268             write (iout,*) "etoti",etoti," deltei",deltei
 269             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 270             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 271      &        k=1,n_ene)
 272 #endif
 273             aux=entfac(i,iprot)+fac*deltei
 274             sumlik(ib,iprot)=sumlik(ib,iprot)+aux*Ptab(jj,ib,iprot)
 275 #ifdef DEBUG
 276             write (iout,'(2i5,f5.2,7(a,e15.5))')
 277      &        i,ib,fac," e_total",etoti,
 278      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 279      &        " e_lowest",elowest(ib,1,iprot)," aux",aux
 280             write (iout,'(i5,5f12.5)')
 281      &             i,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot),
 282      &             aux*Ptab(jj,ib,iprot)
 283 #endif
 284 #ifdef OUT_LIK
 285             write (icsa_bank,'(2i5,4f12.5,e15.5)')
 286      &            i,jj,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot),
 287      &            Ptab(jj,ib,iprot)
 288 #endif
 289             do k=1,n_ene
 290               sumlikder(k,ib,iprot)=sumlikder(k,ib,iprot)+
 291      &        enetb(jj,k,iprot)*escal(k,ib,1,iprot)*Ptab(jj,ib,iprot)
 292             enddo
 293             k=0
 294             do ik=1,ntyp
 295               do jk=1,ik
 296                 k=k+1
 297                 sumlikeps(k,ib,iprot)=sumlikeps(k,ib,iprot)+
 298      &              (ftune_epsprim(eps(ik,jk))*
 299      &              eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))
 300      &              *Ptab(jj,ib,iprot)*ww(1)
 301               enddo
 302             enddo
 303             do k=1,nbacktor_var
 304               sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot)
 305      &         +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(13,ib,1,iprot)
 306             enddo
 307             do k=nbacktor_var+1,ntor_var
 308               sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot)
 309      &         +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(19,ib,1,iprot)
 310             enddo
 311             if (aux.le.150.0) then
 312               aux=dexp(-aux)
 313               sumpart=sumpart+aux
 314               do k=1,n_ene
 315                 sumqder(k,ib,iprot)=sumqder(k,ib,iprot)
 316      &            +aux*enetb(jj,k,iprot)*escal(k,ib,1,iprot)
 317               enddo
 318               k=0
 319               do ik=1,ntyp
 320                 do jk=1,ik
 321                   k=k+1
 322                   sumqeps(k,ib,iprot)=sumqeps(k,ib,iprot)+
 323      &                ww(1)*(ftune_epsprim(eps(ik,jk))*
 324      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))*aux
 325                 enddo
 326               enddo
 327               do k=1,nbacktor_var
 328                 sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot)
 329      &            +aux*entor(k,jj,iprot)*escal(13,ib,1,iprot)
 330               enddo
 331               do k=nbacktor_var+1,ntor_var
 332                 sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot)
 333      &            +aux*entor(k,jj,iprot)*escal(19,ib,1,iprot)
 334               enddo
 335             endif
 336           endif
 337   221     continue
 338         enddo ! j
 339         efree(ib,1,iprot)=sumpart
 340 #ifdef DEBUG
 341           write (iout,*) " ib",ib," iprot",iprot,
 342      &     " sumlik",sumlik(ib,iprot),
 343      &     " sumq",efree(ib,1,iprot)
 344 #endif
 345 #ifdef OUT_LIK
 346        close(icsa_bank)
 347 #endif
 348       ENDDO ! ib
 349 C Calculation of heat capacity
 350       DO IB=1,NBETA(2,IPROT)
 351         elowest_all=elowest(ib,2,iprot)
 352 #ifdef DEBUG
 353         write(iout,*) "iprot",iprot," ib",ib,
 354      &      " elowest",elowest(ib,2,iprot)
 355         write(iout,*) "escalbis",(escalbis(k,ib,2,iprot),k=1,n_ene)
 356 #endif
 357         fac=betaT(ib,2,iprot)
 358         temper=1.0d0/(fac*Rgas)
 359         sumpart=0.0d0
 360         emeani=0.0d0
 361         esquarei=0.0d0
 362         ebisi=0.0d0
 363 #ifdef DEBUG
 364         write (iout,*) "Before sum, ib=",ib
 365         write (iout,*) "efree",efree(ib,2,iprot),sumpart,
 366      &     " emean",emean_ftot(ib,iprot),
 367      &     " esquare",esquare_ftot(ib,iprot),
 368      &     " ebis",ebis_ftot(ib,iprot),
 369      &     " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
 370 #endif
 371         do j=1,ntot_work(iprot)
 372           i = tsil(j,iprot)
 373 #ifdef DEBUG
 374           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 375      &      " indstart",indstart(me1,iprot),
 376      &      " indend",indend(me1,iprot)
 377 #endif
 378           if (i.eq.0) cycle
 379           jj = i2ii(i,iprot)
 380 #ifdef DEBUG
 381           if (jj.gt.0) then
 382             write (iout,*) "ib",ib," j",j,
 383             write (iout,*) "nu",(nu(k,jj,iprot),k=1,
 384      &        natconstr(iprot))
 385           endif
 386 #endif
 387           if (jj.gt.0) then
 388 c Temperature-dependent energy
 389             etoti=0.0d0
 390             etoti_orig=0.0d0
 391             eprim=0.0d0
 392             ebis=0.0d0
 393             do k=1,n_ene
 394               etoti=etoti+ww(k)*escal(k,ib,2,iprot)
 395      &            *enetb(jj,k,iprot)
 396               eprim=eprim+ww(k)*escalprim(k,ib,2,iprot)
 397      &            *enetb(jj,k,iprot)
 398               ebis=ebis+ww(k)*escalbis(k,ib,2,iprot)
 399      &            *enetb(jj,k,iprot)
 400             enddo
 401             deltei=etoti-elowest(ib,2,iprot)
 402 #ifdef DEBUG
 403             write (iout,*) "etoti",etoti," deltei",deltei
 404             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 405             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 406      &        k=1,n_ene)
 407 #endif
 408             aux=entfac(i,iprot)+fac*deltei
 409 #ifdef DEBUG
 410             write (iout,'(f5.2,7(a,e15.5))')
 411      &        fac," e_total",etoti,
 412      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 413      &        " eprim",eprim," ebis",ebis,
 414      &        " e_lowest",elowest(ib,2,iprot)," aux",aux
 415 #endif
 416             if (aux.le.150.0) then
 417               aux=dexp(-aux)
 418               sumpart=sumpart+aux
 419               etoti=etoti-temper*eprim
 420               emeani=emeani+etoti*aux
 421               esquarei=esquarei+aux*etoti**2
 422               ebisi=ebisi+aux*ebis
 423               do k=1,n_ene
 424                 eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
 425      &           +aux*enetb(jj,k,iprot)*escal(k,ib,2,iprot)
 426 c                write (iout,*) "eave_pf",eave_pf(k,ii,ib,iprot)
 427                 eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
 428      &           +aux*enetb(jj,k,iprot)*(escal(k,ib,2,iprot)-
 429      &            temper*escalprim(k,ib,2,iprot))
 430                 eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
 431      &           +aux*enetb(jj,k,iprot)*escalbis(k,ib,2,iprot)
 432                 emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)
 433      &           +aux*enetb(jj,k,iprot)*etoti*escal(k,ib,2,iprot)
 434                 emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
 435      &           +aux*enetb(jj,k,iprot)*etoti*(escal(k,ib,2,iprot)
 436      &            -temper*escalprim(k,ib,2,iprot))
 437                 emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
 438      &           +aux*enetb(jj,k,iprot)*ebis*escal(k,ib,2,iprot)
 439                 emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
 440      &           +aux*enetb(jj,k,iprot)*etoti**2*escal(k,ib,2,iprot)
 441               enddo
 442               k=0
 443               do ik=1,ntyp
 444                 do jk=1,ik
 445                   k=k+1
 446                   enepsjk=ftune_epsprim(eps(ik,jk))*
 447      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
 448                   enepsjk=enepsjk*ww(1)
 449                   enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
 450      &             +aux*enepsjk
 451                   eneps_mix_ftot(k,ib,iprot)=
 452      &             eneps_mix_ftot(k,ib,iprot)+aux*enepsjk*etoti
 453                   eneps_mix_fbistot(k,ib,iprot)=
 454      &             eneps_mix_fbistot(k,ib,iprot)+aux*enepsjk*ebis
 455                   eneps_mixsq_ftot(k,ib,iprot)=
 456      &             eneps_mixsq_ftot(k,ib,iprot)+aux*enepsjk*etoti**2
 457                 enddo
 458               enddo
 459               do k=1,nbacktor_var
 460                   entorjk=entor(k,jj,iprot)
 461 c                  write (iout,*) " k"," ik",ik," jk",jk,
 462 c     &             " entor",entorjk," contirb",
 463 c     &             aux*entorjk*escal(13,ib,2,iprot)
 464                   entorave_ftot(k,ib,iprot)=
 465      &             entorave_ftot(k,ib,iprot)
 466      &             +aux*entorjk*escal(13,ib,2,iprot)
 467 c                  write (iout,*) "entorave_f",
 468 c     &              entorave_ftot(k,ib,iprot)
 469                   entorave_fprimtot(k,ib,iprot)=
 470      &             entorave_fprimtot(k,ib,iprot)
 471      &             +aux*entorjk*(escal(13,ib,2,iprot)-
 472      &              temper*escalprim(13,ib,2,iprot))
 473                   entorave_fbistot(k,ib,iprot)=
 474      &             entorave_fbistot(k,ib,iprot)
 475      &             +aux*entorjk*escalbis(13,ib,2,iprot)
 476                   entor_mix_ftot(k,ib,iprot)=
 477      &             entor_mix_ftot(k,ib,iprot)
 478      &             +aux*entorjk*etoti*escal(13,ib,2,iprot)
 479                   entor_mix_fprimtot(k,ib,iprot)=
 480      &             entor_mix_fprimtot(k,ib,iprot)
 481      &             +aux*entorjk*etoti*(escal(13,ib,2,iprot)-
 482      &            temper*escalprim(13,ib,2,iprot))
 483                   entor_mix_fbistot(k,ib,iprot)=
 484      &             entor_mix_fbistot(k,ib,iprot)
 485      &             +aux*entorjk*ebis*escal(13,ib,2,iprot)
 486                   entor_mixsq_ftot(k,ib,iprot)=
 487      &             entor_mixsq_ftot(k,ib,iprot)
 488      &             +aux*entorjk*etoti**2*escal(13,ib,2,iprot)
 489               enddo
 490               do k=nbacktor_var+1,ntor_var
 491                   entorjk=entor(k,jj,iprot)
 492 c                  write (iout,*) " k"," ik",ik," jk",jk,
 493 c     &             " entor",entorjk," contirb",
 494 c     &             aux*entorjk*escal(19,ib,2,iprot)
 495                   entorave_ftot(k,ib,iprot)=
 496      &             entorave_ftot(k,ib,iprot)
 497      &             +aux*entorjk*escal(19,ib,2,iprot)
 498 c                  write (iout,*) "entorave_f",
 499 c     &              entorave_ftot(k,ib,iprot)
 500                   entorave_fprimtot(k,ib,iprot)=
 501      &             entorave_fprimtot(k,ib,iprot)
 502      &             +aux*entorjk*(escal(19,ib,2,iprot)-
 503      &              temper*escalprim(19,ib,2,iprot))
 504                   entorave_fbistot(k,ib,iprot)=
 505      &             entorave_fbistot(k,ib,iprot)
 506      &             +aux*entorjk*escalbis(19,ib,2,iprot)
 507                   entor_mix_ftot(k,ib,iprot)=
 508      &             entor_mix_ftot(k,ib,iprot)
 509      &             +aux*entorjk*etoti*escal(19,ib,2,iprot)
 510                   entor_mix_fprimtot(k,ib,iprot)=
 511      &             entor_mix_fprimtot(k,ib,iprot)
 512      &             +aux*entorjk*etoti*(escal(19,ib,2,iprot)-
 513      &            temper*escalprim(13,ib,2,iprot))
 514                   entor_mix_fbistot(k,ib,iprot)=
 515      &             entor_mix_fbistot(k,ib,iprot)
 516      &             +aux*entorjk*ebis*escal(19,ib,2,iprot)
 517                   entor_mixsq_ftot(k,ib,iprot)=
 518      &             entor_mixsq_ftot(k,ib,iprot)
 519      &             +aux*entorjk*etoti**2*escal(19,ib,2,iprot)
 520               enddo
 521             endif
 522           endif
 523         enddo
 524         efree(ib,2,iprot)=sumpart
 525         emean_ftot(ib,iprot)=emeani
 526         ebis_ftot(ib,iprot)=ebisi
 527         esquare_ftot(ib,iprot)=esquarei
 528 #ifdef DEBUG
 529         write (iout,*) "After sum, ib=",ib
 530         write (iout,*) "efree",efree(ib,2,iprot),sumpart,
 531      &     " emean",emean_ftot(ib,iprot),
 532      &     " esquare",esquare_ftot(ib,iprot),
 533      &     " ebis",ebis_ftot(ib,iprot),
 534      &     " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
 535 #endif
 536         ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)*temper
 537         do k=1,n_ene
 538           eave_pfbistot(k,ib,iprot)=
 539      &      eave_pfbistot(k,ib,iprot)*temper
 540           emix_pfbistot(k,ib,iprot)=
 541      &      emix_pfbistot(k,ib,iprot)*temper
 542         enddo
 543         do k=1,nntyp
 544           eneps_mix_fbistot(k,ib,iprot)=
 545      &     eneps_mix_fbistot(k,ib,iprot)*temper
 546         enddo
 547         do k=1,ntor_var
 548           entorave_fbistot(k,ib,iprot)=
 549      &     entorave_fbistot(k,ib,iprot)*temper
 550           entor_mix_fbistot(k,ib,iprot)=
 551      &     entor_mix_fbistot(k,ib,iprot)*temper
 552         enddo
 553 #ifdef DEBUG
 554         write (iout,*) "eave_pf",(eave_pftot(k,ib,iprot),
 555      &      k=1,n_ene)
 556         write (iout,*) "entorave_f",(entorave_ftot(k,ib,iprot),
 557      &      k=1,ntor_var)
 558 #endif
 559 #ifdef DEBUG
 560         write (iout,*) "ib",ib," temper",temper,
 561      &       " ebis",ebis_ftot(ib,iprot)
 562 #endif
 563       ENDDO ! ib
 564 C Calculation of conformation-dependent averages
 565       DO INAT=1,NATLIKE(IPROT)
 566       DO IB=1,NBETA(I+2,IPROT)
 567         elowest_all=elowest(ib,inat+2,iprot)
 568 #ifdef DEBUG
 569           write(iout,*) "iprot",iprot," ib",ib,
 570      &      " elowest",elowest(ib,iprot)
 571 #endif
 572         fac=betaT(ib,2,iprot)
 573         temper=1.0d0/(fac*Rgas)
 574         emeani=0.0d0
 575         do j=1,ntot_work(iprot)
 576           i = tsil(j,iprot)
 577 #ifdef DEBUG
 578           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 579      &      " indstart",indstart(me1,iprot),
 580      &      " indend",indend(me1,iprot)
 581 #endif
 582           if (i.eq.0) cycle
 583           jj = i2ii(i,iprot)
 584 #ifdef DEBUG
 585           if (jj.gt.0) then
 586             write (iout,*) "ib",ib," j",j,
 587             write (iout,*) "nu",(nu(k,jj,iprot),k=1,
 588      &        natconstr(iprot))
 589           endif
 590 #endif
 591           if (jj.gt.0) then
 592 c Temperature-dependent energy
 593             etoti=0.0d0
 594             do k=1,n_ene
 595               etoti=etoti+ww(k)*escal(k,ib,inat+2,iprot)
 596      &            *enetb(jj,k,iprot)
 597             enddo
 598             deltei=etoti-elowest(ib,inat+2,iprot)
 599 #ifdef DEBUG
 600             write (iout,*) "etoti",etoti," deltei",deltei
 601             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 602             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 603      &        k=1,n_ene)
 604 #endif
 605             aux=entfac(i,iprot)+fac*deltei
 606 #ifdef DEBUG
 607             write (iout,'(f5.2,7(a,e15.5))')
 608      &        fac," e_total",etoti,
 609      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 610      &        " eprim",eprim," ebis",ebis,
 611      &        " e_lowest",elowest(ib,inat+2,iprot)," aux",aux
 612 #endif
 613             if (aux.le.150.0) then
 614               aux=dexp(-aux)
 615               sumpart=sumpart+aux
 616 c 4/13/04 AL Components of the conformation-dependent averages
 617               do k=1,n_ene
 618                 eave_nat_pftot(k,ib,inat,iprot)=
 619      &           eave_nat_pftot(k,ib,inat,iprot)
 620      &           +aux*enetb(jj,k,iprot)*escal(k,ib,inat+2,iprot)
 621               enddo
 622               k=0
 623               do ik=1,ntyp
 624                 do jk=1,ik
 625                   k=k+1
 626                   enepsjk=ftune_epsprim(eps(ik,jk))*
 627      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
 628                   enepsjk=enepsjk*ww(1)
 629                   enepsave_nat_ftot(k,ib,inat,iprot)=
 630      &             enepsave_nat_ftot(k,ib,inat,iprot)
 631      &             +aux*enepsjk
 632                 enddo
 633               enddo
 634               k=0
 635               do k=1,nbacktor_var
 636                   entorjk=entor(k,jj,iprot)
 637 c                  write (iout,*) " k",
 638 c     &             " entor",entorjk," contirb",
 639 c     &             aux*entorjk*escal(13,ib,inat+2,iprot)
 640                   entorave_nat_ftot(k,ib,inat,iprot)=
 641      &             entorave_nat_ftot(k,ib,inat,iprot)
 642      &             +aux*entorjk*escal(13,ib,inat+2,iprot)
 643               enddo
 644               do k=nbacktor_var+1,ntor_var
 645                 do l=1,3
 646                   entorjk=entor(k,jj,iprot)
 647 c                  write (iout,*) " k",
 648 c     &             " entor",entorjk," contirb",
 649 c     &             aux*entorjk*escal(19,ib,inat+2,iprot)
 650                   entorave_nat_ftot(k,ib,inat,iprot)=
 651      &             entorave_nat_ftot(k,ib,inat,iprot)
 652      &             +aux*entorjk*escal(19,ib,inat+2,iprot)
 653                 enddo
 654               enddo
 655               do l=1,natdim(inat,iprot)
 656                 nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
 657      &            +aux*nu(l,inat,jj,iprot)
 658                 do k=1,n_ene
 659                   nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
 660      &           +aux*enetb(jj,k,iprot)*nu(k,inat,jj,iprot)*
 661      &            escal(k,ib,inat+2,iprot)
 662                 enddo
 663                 k=0
 664                 do ik=1,ntyp
 665                   do jk=1,ik
 666                     k=k+1
 667                     nuepsave_f(k,l,ib,inat,iprot)=
 668      &               nuepsave_f(k,l,ib,inat,iprot)+aux*enepsjk*
 669      &               nu(l,inat,jj,iprot)
 670                   enddo
 671                 enddo
 672                 k=0
 673                 do ik=0,ntortyp-1
 674                   do jk=-ntortyp,ntortyp
 675                     do iblock=1,2
 676                     if (mask_tor(0,jk,ik,iblock).gt.0) then
 677                       k=k+1
 678                       nutorave_f(k,l,ib,inat,iprot)=
 679      &                 nutorave_f(k,l,ib,inat,iprot)
 680      &                 +aux*entorjk*escal(13,ib,inat+2,iprot)*
 681      &                  nu(l,inat,jj,iprot)
 682                     endif
 683                     enddo
 684                   enddo
 685                 enddo
 686                 do ik=1,nsccortyp
 687                   do jk=1,nsccortyp
 688                     do ll=1,3
 689                     if (mask_tor(ll,jk,ik,1).gt.0) then
 690                       k=k+1
 691                       nutorave_f(k,l,ib,inat,iprot)=
 692      &                 nutorave_f(k,l,ib,inat,iprot)
 693      &                 +aux*entorjk*escal(13,ib,inat+2,iprot)*
 694      &                  nu(l,inat,jj,iprot)
 695                     endif
 696                     enddo
 697                   enddo
 698                 enddo
 699               enddo
 700             endif
 701           endif
 702 #ifdef DEBUG
 703           write (iout,*) "iprot",iprot," ib",ib,
 704           write (iout,*) "nuave"
 705           write (iout,'(20f10.5)') (nuave(k,ib,iprot),
 706      &            k=1,natconstr(iprot))
 707 #endif
 708 #ifdef DEBUG
 709           write (iout,*) "inat",inat,
 710      &      " efree_nat",efree_nat(ib,inat,iprot)
 711 #endif
 712         enddo
 713 #ifdef DEBUG
 714           write (iout,*) "iprot",iprot," ib",ib
 715           write (iout,*) "nuave0"
 716           write (iout,'(20f10.5)') (nuave(k,ib,inat,iprot),
 717      &            k=1,natdim(inat,iprot))
 718 #endif
 719       ENDDO ! ib
 720       ENDDO ! INAT
 721       return
 722       end
 723 c-------------------------------------------------------------------------------
 724       subroutine ave_sum(iprot)
 725       implicit none
 726       include "DIMENSIONS"
 727       include "DIMENSIONS.ZSCOPT"
 728 #ifdef MPI
 729       include "mpif.h"
 730       integer IERROR,ErrCode
 731       include "COMMON.MPI"
 732 #endif
 733       integer n,nf,What
 734       include "COMMON.WEIGHTS"
 735       include "COMMON.WEIGHTDER"
 736       include "COMMON.ENERGIES"
 737       include "COMMON.IOUNITS"
 738       include "COMMON.VMCPAR"
 739       include "COMMON.NAMES"
 740       include "COMMON.INTERACT"
 741       include "COMMON.TIME1"
 742       include "COMMON.CHAIN"
 743       include "COMMON.PROTFILES"
 744       include "COMMON.VAR"
 745       include "COMMON.COMPAR"
 746 C Define local variables
 747       integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,inat
 748       integer ipass_conf,istart_conf,iend_conf
 749       double precision energia(0:max_ene)
 750       double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk,
 751      &  eave_pft(max_ene),emix_pft(max_ene),eave_pfprimt(max_ene),
 752      &  eave_pfbist(max_ene),emix_pfprimt(max_ene),emix_pfbist(max_ene),
 753      &  esquare_ft,efree_t,
 754      &  emixsq_pft(max_ene),eneps_mixt(nntyp),eneps_meant(nntyp),
 755      &  enepsave_ft(nntyp),eneps_mix_ft(nntyp),entorave_ft(maxtor_var),
 756      &  entor_mix_ft(maxtor_var),
 757      &  eneps_mix_fbist(nntyp),entorave_fbist(maxtor_var),
 758      &  entorave_fprimt(maxtor_var),entor_mix_fprimt(maxtor_var),
 759      &  entor_mix_fbist(maxtor_var),eneps_mixsq_ft(nntyp),
 760      &  entor_mixsq_ft(maxtor_var),emean_ft,ebis_ft,nuave_t(maxdimnat),
 761      &  nu_pft(max_ene,maxdimnat),
 762      &  nuepsave_ft(nntyp,maxdimnat),
 763      &  nutorave_ft(maxtor_var,maxdimnat),
 764      &  sumlik_t,sumlikder_t(max_ene),sumlikeps_t(nntyp),
 765      &  sumliktor_t(maxtor_var),
 766      &  sumq_t,sumqder_t(max_ene),sumqeps_t(nntyp),sumqtor_t(maxtor_var)
 767
 768       double precision aux,fac,ef1,ef2,em1,em2,var1,var2,efree_tot,facF,
 769      &  elowest_all
 770       double precision gnmr,tcpu_ini,tcpu
 771       logical lprn
 772       integer icant
 773       external icant
 774
 775       lprn=.true.
 776 c Maximum likelihood contribution
 777       DO IB=1,NBETA(1,IPROT)
 778         fac=betaT(ib,1,iprot)
 779 #ifdef MPI
 780 #ifdef DEBUG
 781         write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
 782         write (iout,*) "iprot",iprot," ib",ib
 783         write (iout,*) "sumlik",sumlik(ib,iprot),
 784      &     " sumq",efree(ib,1,iprot)
 785         call flush(iout)
 786 #endif
 787         call MPI_Reduce( sumlik(ib,iprot), sumlik_t,1,
 788      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 789         call MPI_Reduce( efree(ib,1,iprot), sumq_t,1,
 790      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 791 #ifdef DEBUG
 792         write (iout,*) "sumlik",sumlik(ib,iprot)," sumlik_t",sumlik_t
 793         write (iout,*) "sumq",efree(ib,1,iprot)," sumq_t",sumq_t
 794         call flush(iout)
 795 #endif
 796         call MPI_Reduce(sumlikder(1,ib,iprot),sumlikder_t(1),
 797      &      n_ene,
 798      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 799         call MPI_Reduce(sumlikeps(1,ib,iprot),
 800      &      sumlikeps_t(1), nntyp,
 801      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 802         call MPI_Reduce(sumliktor(1,ib,iprot),
 803      &      sumliktor_t(1), ntor_var,
 804      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 805         call MPI_Reduce(sumqder(1,ib,iprot),sumqder_t(1),n_ene,
 806      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 807         call MPI_Reduce(sumqeps(1,ib,iprot),sumqeps_t(1), nntyp,
 808      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 809         call MPI_Reduce(sumqtor(1,ib,iprot),sumqtor_t(1), ntor_var,
 810      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 811 #ifdef DEBUG
 812         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
 813         call flush(iout)
 814 #endif
 815         if (me1.eq.master) then
 816 #ifdef DEBUG
 817           write (iout,*) "ib",ib,
 818      &      " elowest",elowest(ib,1,iprot),
 819      &      " sumlik",sumlik_t," qsum",sumq_t," fac",fac
 820 #endif
 821           do k=1,n_ene
 822             sumlikder(k,ib,iprot)=fac*(
 823      &          sumlikder_t(k)
 824      &          -sumqder_t(k)/sumq_t)
 825           enddo
 826           do k=1,nntyp
 827             sumlikeps(k,ib,iprot)=fac*(
 828      &          sumlikeps_t(k)-sumqeps_t(k)/sumq_t)
 829           enddo
 830 c          write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
 831 c     &          eave_pftot(19,ib,iprot)
 832           do k=1,ntor_var
 833             sumliktor(k,ib,iprot)=fac*(
 834      &          sumliktor_t(k)-sumqtor_t(k)/sumq_t)
 835           enddo
 836           sumlik(ib,iprot)=sumlik_t+dlog(sumq_t)
 837           efree(ib,1,iprot)=sumq_t
 838 #ifdef DEBUG
 839           write (iout,*) "ib",ib," iprot",iprot," final sumlik",
 840      &       sumlik(ib,iprot)," sumq",efree(ib,1,iprot)
 841 #endif
 842         endif
 843 #else
 844         do k=1,n_ene
 845           sumlikder(k,ib,iprot)=fac*(
 846      &          sumlikder(k,ib,iprot)
 847      &         -sumqder(k,ib,iprot)/efree(ib,1,iprot))
 848         enddo
 849         do k=1,nntyp
 850           sumlikeps(k,ib,iprot)=fac*(
 851      &        sumlikeps(k,ib,iprot)
 852      &       -sumqeps(k,ib,iprot)/efree(ib,1,iprot))
 853         enddo
 854         do k=1,ntor_var
 855           sumliktor(k,ib,iprot)=fac*(
 856      &      sumliktor(k,ib,iprot)
 857      &     -sumqtor(k,ib,iprot)/efree(ib,1,iprot))
 858         enddo
 859         sumlik(ib,iprot)=sumlik(ib,iprot)+dlog(efree(ib,1,iprot)
 860 c     &      -elowest(nbeta(iprot)+ib,iprot)*fac
 861 #endif
 862         ii=0
 863         do k=1,n_ene
 864           if (imask(k).gt.0) then
 865             ii=ii+1
 866             sumlikder(ii,ib,iprot)=sumlikder(k,ib,iprot)
 867           endif
 868         enddo
 869       ENDDO ! ib
 870 c Heat capacity and averages
 871       DO IB=1,NBETA(2,IPROT)
 872         fac=betaT(ib,2,iprot)
 873 #ifdef MPI
 874 #ifdef DEBUG
 875         write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
 876         write (iout,*) "iprot",iprot," ib",ib
 877         call flush(iout)
 878 #endif
 879         call MPI_Reduce( efree(ib,2,iprot), efree_t,1,
 880      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 881 #ifdef DEBUG
 882         write (iout,*) "efree",efree(iprot
 883         write (iout,*) "efree_t",efree_t
 884         call flush(iout)
 885 #endif
 886         call MPI_Reduce(emean_ftot(ib,iprot),emean_ft,1,
 887      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 888         call MPI_Reduce(ebis_ftot(ib,iprot),ebis_ft,1,
 889      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 890         call MPI_Reduce(esquare_ftot(ib,iprot),esquare_ft,1,
 891      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 892         call MPI_Reduce(eave_pftot(1,ib,iprot),eave_pft,
 893      &    n_ene,
 894      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 895         call MPI_Reduce(eave_pfprimtot(1,ib,iprot),
 896      &    eave_pfprimt,n_ene,
 897      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 898         call MPI_Reduce(eave_pfbistot(1,ib,iprot),
 899      &    eave_pfbist,n_ene,
 900      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 901 c          write (iout,*) "eave_pf",(eave_pf(k,iprot),
 902 c     &      k=1,n_ene)
 903 c          write (iout,*) "eave_pft",(eave_pft(k),k=1,n_ene)
 904         call MPI_Reduce(emix_pftot(1,ib,iprot),emix_pft,
 905      &    n_ene,
 906      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 907         call MPI_Reduce(emix_pfprimtot(1,ib,iprot),
 908      &    emix_pfprimt,n_ene,
 909      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 910         call MPI_Reduce(emix_pfbistot(1,ib,iprot),
 911      &    emix_pfbist,n_ene,
 912      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 913         call MPI_Reduce(emixsq_pftot(1,ib,iprot),
 914      &    emixsq_pft, n_ene,
 915      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 916         call MPI_Reduce(enepsave_ftot(1,ib,iprot),
 917      &    enepsave_ft, nntyp,
 918      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 919         call MPI_Reduce(eneps_mix_ftot(1,ib,iprot),
 920      &    eneps_mix_ft,nntyp,
 921      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 922         call MPI_Reduce(eneps_mix_fbistot(1,ib,iprot),
 923      &    eneps_mix_fbist,nntyp,
 924      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 925         call MPI_Reduce(eneps_mixsq_ftot(1,ib,iprot),
 926      &    eneps_mixsq_ft,nntyp,
 927      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 928 c          write (iout,*) "enepsave_f",(enepsave_f(k,iprot),
 929 c     &      k=1,nntyp)
 930 c          write (iout,*) "enepsave_ft",(enepsave_ft(k),
 931 c     &      k=1,nntyp)
 932         call MPI_Reduce(entorave_ftot(1,ib,iprot),
 933      &    entorave_ft, ntor_var,
 934      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 935         call MPI_Reduce(entorave_fprimtot(1,ib,iprot),
 936      &    entorave_fprimt, ntor_var,
 937      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 938         call MPI_Reduce(entorave_fbistot(1,ib,iprot),
 939      &    entorave_fbist, ntor_var,
 940      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 941         call MPI_Reduce(entor_mix_ftot(1,ib,iprot),
 942      &    entor_mix_ft,ntor_var,
 943      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 944         call MPI_Reduce(entor_mix_fprimtot(1,ib,iprot),
 945      &    entor_mix_fprimt,ntor_var,
 946      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 947         call MPI_Reduce(entor_mix_fbistot(1,ib,iprot),
 948      &   entor_mix_fbist,ntor_var,
 949      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 950         call MPI_Reduce(entor_mixsq_ftot(1,ib,iprot),
 951      &    entor_mixsq_ft,ntor_var,
 952      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 953 c          write (iout,*) "entorave_f",(entorave_f(k,iprot),
 954 c     &      k=1,ntor_var)
 955 c          write (iout,*) "entorave_ft",(entorave_ft(k),
 956 c     &      k=1,ntor_var)
 957 #ifdef DEBUG
 958         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
 959         call flush(iout)
 960 #endif
 961         if (me1.eq.master) then
 962           elowest_all=elowest(ib,2,iprot)
 963 #ifdef DEBUG
 964           write (iout,*) "ib",ib,
 965      &      " elowest",elowest(ib,iprot),
 966      &      "efree",efree_t," fac",fac," facF",facF,
 967      &      " efree_tot",efree_tot
 968 #endif
 969           do k=1,n_ene
 970             eave_pftot(k,ib,iprot)=eave_pft(k)/efree_t
 971             eave_pfprimtot(k,ib,iprot)=eave_pfprimt(k)/efree_t
 972             eave_pfbistot(k,ib,iprot)=eave_pfbist(k)/efree_t
 973             emix_pftot(k,ib,iprot)=emix_pft(k)/efree_t
 974             emix_pfprimtot(k,ib,iprot)=emix_pfprimt(k)/efree_t
 975             emix_pfbistot(k,ib,iprot)=emix_pfbist(k)/efree_t
 976             emixsq_pftot(k,ib,iprot)=emixsq_pft(k)/efree_t
 977           enddo
 978           do k=1,nntyp
 979             enepsave_ftot(k,ib,iprot)=enepsave_ft(k)/efree_t
 980             eneps_mix_ftot(k,ib,iprot)=eneps_mix_ft(k)/efree_t
 981             eneps_mix_fbistot(k,ib,iprot)=eneps_mix_fbist(k)/efree_t
 982             eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ft(k)/efree_t
 983           enddo
 984 c          write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
 985 c     &          eave_pftot(19,ib,iprot)
 986           do k=1,ntor_var
 987             entorave_ftot(k,ib,iprot)=entorave_ft(k)/
 988      &        efree_t
 989 c            write (iout,*) "iprot",iprot," ib",ib," k",k,
 990 c     &         " entorave_ftot", entorave_ftot(k,ib,iprot)
 991             entorave_fprimtot(k,ib,iprot)=entorave_fprimt(k)/efree_t
 992             entorave_fbistot(k,ib,iprot)=entorave_fbist(k)/efree_t
 993             entor_mix_ftot(k,ib,iprot)=entor_mix_ft(k)/efree_t
 994             entor_mix_fprimtot(k,ib,iprot)=entor_mix_fprimt(k)/efree_t
 995             entor_mix_fbistot(k,ib,iprot)=entor_mix_fbist(k)/efree_t
 996             entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ft(k)/efree_t
 997           enddo
 998           emean_ftot(ib,iprot)=emean_ft/efree_t
 999           ebis_ftot(ib,iprot)=ebis_ft/efree_t
1000           esquare_ftot(ib,iprot)=esquare_ft/efree_t
1001           efree(ib,2,iprot)=-dlog(efree_t)/fac+elowest(ib,2,iprot)
1002         endif
1003 #else
1004         do k=1,n_ene
1005           eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
1006      &      /efree(ib,2,iprot)
1007           eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
1008      &      /efree(ib,2,iprot)
1009           eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
1010      &      /efree(ib,2,iprot)
1011           emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)/efree(ib,iprot)
1012           emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
1013      &      /efree(ib,2,iprot)
1014           emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
1015      &      /efree(ib,2,iprot)
1016           emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
1017      &      /efree(ib,2,iprot)
1018         enddo
1019         do k=1,nntyp
1020           enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
1021      &      /efree(ib,2,iprot)
1022           eneps_mix_ftot(k,ib,iprot)=eneps_mix_ftot(k,ib,iprot)
1023      &      /efree(ib,2,iprot)
1024           eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ftot(k,ib,iprot)
1025      &      /efree(ib,2,iprot)
1026         enddo
1027         do k=1,ntor_var
1028           entorave_ftot(k,ib,iprot)=entorave_f(k,ib,iprot)
1029      &      /efree(ib,2,iprot)
1030           entor_mix_ftot(k,ib,iprot)=entor_mix_ftot(k,ib,iprot)
1031      &      /efree(ib,2,iprot)
1032           entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ftot(k,ib,iprot)
1033      &      /efree(ib,2,iprot)
1034         enddo
1035         emean_ftot(ib,iprot)=emean_ftot(ib,iprot)/efree(ib,2,iprot)
1036         ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)/efree(ib,2,iprot)
1037         esquare_ftot(ib,iprot)=esquare_ftot(ib,iprot)/efree(ib,2,iprot)
1038 #endif
1039       ENDDO ! IB
1040 c 4/13/04 AL Components of the correlation coefficients and their derivatives
1041       DO INAT=1,NATLIKE(IPROT)
1042         DO IB=1,NBETA(INAT+2,IPROT)
1043 #ifdef MPI
1044           call MPI_Reduce( efree(ib,inat+2,iprot), efree_t,1,
1045      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1046           call MPI_Reduce( eave_nat_pftot(1, ib,inat,iprot), eave_pft,
1047      &      n_ene,
1048      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1049           call MPI_Reduce( enepsave_nat_ftot(1, ib,inat,iprot),
1050      &      enepsave_ft, nntyp,
1051      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1052           call MPI_Reduce(entorave_nat_ftot(1,ib,inat,iprot),
1053      &      entorave_ft, ntor_var,
1054      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1055           call MPI_Reduce(nuave(1,ib,inat,iprot),nuave_t(1),
1056      &      natdim(inat,iprot), MPI_DOUBLE_PRECISION, MPI_SUM,
1057      &      Master, Comm1, IERROR)
1058 #ifdef DEBUG
1059           write (iout,*) "After REDUCE nuave",iprot,ib
1060           write (iout,'(20f10.5)')
1061      &      (nuave(l,ib,iprot),l=1,natconstr(iprot))
1062 #endif
1063           call MPI_Reduce(nu_pf(1,1,ib,inat,iprot),nu_pft(1,1),
1064      &      max_ene*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1065      &      MPI_SUM, Master, Comm1, IERROR)
1066           call MPI_Reduce(nuepsave_f(1,1,ib,inat,iprot),
1067      &      nuepsave_ft(1,1),
1068      &      nntyp*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1069      &      MPI_SUM, Master, Comm1, IERROR)
1070           call MPI_Reduce(nutorave_f(1,1,ib,inat,iprot),
1071      &      nutorave_ft(1,1),
1072      &      maxtor_var*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1073      &      MPI_SUM, Master, Comm1, IERROR)
1074 #ifdef DEBUG
1075         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
1076         call flush(iout)
1077 #endif
1078           if (me1.eq.master) then
1079 #ifdef DEBUG
1080             write (iout,*) "ib",ib,
1081      &      " elowest",elowest(ib,iprot),
1082      &      "efree",efree_t," fac",fac,
1083      &      " efree_tot",efree_tot
1084 #endif
1085             do l=1,natdim(inat,iprot)
1086               do k=1,n_ene
1087                 nu_pf(k,l,ib,inat,iprot)=nu_pft(k,l)/efree_t
1088               enddo
1089               do k=1,nntyp
1090                 nuepsave_f(k,l,ib,inat,iprot)=nuepsave_ft(k,l)/efree_t
1091               enddo
1092               do k=1,ntor_var
1093                 nutorave_f(k,l,ib,inat,iprot)=nutorave_ft(k,l)/efree_t
1094               enddo
1095               nuave(l,ib,inat,iprot)=nuave_t(l)/efree_t
1096             enddo
1097             do k=1,n_ene
1098               eave_nat_pftot(k,ib,inat,iprot)=eave_pft(k)/efree_t
1099             enddo
1100             do k=1,nntyp
1101               enepsave_nat_ftot(k,ib,inat,iprot)=enepsave_ft(k)/efree_t
1102             enddo
1103             do k=1,ntor_var
1104               entorave_nat_ftot(k,ib,inat,iprot)=entorave_ft(k)/efree_t
1105             enddo
1106           endif
1107 #else
1108           do l=1,natdim(inat,iprot)
1109             do k=1,n_ene
1110               nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
1111      &         /efree(ib,inat+2,iprot)
1112             enddo
1113             do k=1,nntyp
1114              nuepsave_f(k,l,ib,inat,iprot)=nuepsave_f(k,l,ib,inat,iprot)
1115      &            /efree(ib,inat+2,iprot)
1116             enddo
1117             do k=1,ntor_var
1118               nutorave_ftot(k,l,ib,inat,iprot)=
1119      &          nutorave_ftot(k,l,ib,inat,iprot)
1120      &          /efree(ib,inat+2,iprot)
1121             enddo
1122             nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
1123      &        /efree(ib,inat+2,iprot)
1124           enddo ! l
1125           do k=1,n_ene
1126             eave_nat_pftot(k,ib,inat,iprot)=
1127      &       eave_nat_pftot(k,ib,inat,iprot)
1128      &       /efree(ib,inat+2,iprot)
1129           enddo
1130           do k=1,nntyp
1131             enepsave_nat_ftot(k,ib,inat,iprot)=
1132      &         enepsave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1133           enddo
1134           do k=1,ntor_var
1135             enetorave_nat_ftot(k,ib,inat,iprot)=
1136      &        enetorave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1137           enddo
1138 #endif
1139         ENDDO ! IB
1140 #ifdef DEBUG
1141         write (iout,*) "ib",ib," efree_tot",efree_tot
1142 #endif
1143       ENDDO ! INAT
1144
1145       return
1146       end