source/maxlik/src_FPy/averages_sc.F.safe

   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 C Define local variables
 203       integer i,ii,iii,jj,j,k,kk,l,m,ik,jk,iprot,ib,inat,
 204      &  ic
 205       integer ipass_conf,istart_conf,iend_conf
 206       double precision energia(0:max_ene)
 207       double precision etoti,sumpart,esquarei,emeani,
 208      &  elowesti,enepsjk,eprim,ebis,eprimi,ebisi,etoti_orig,
 209      &  entorjk,eave_pft(max_ene),emix_pft(max_ene),
 210      &  esquare_ft,efree_t,emixsq_pft(max_ene),
 211      &  eneps_mixt(nntyp),eneps_meant(nntyp),
 212      &  enepsave_ft(nntyp),eneps_mix_ft(nntyp),
 213      &  eneps_mixsq_ft(nntyp),emean_ft
 214       double precision aux,auxf,fac,facf,ef1,ef2,em1,em2,var1,var2,
 215      &  deltei,deltei_orig,temper,elowest_all
 216       double precision gnmr,tcpu_ini,tcpu
 217       double precision ftune_epsprim
 218       external ftune_epsprim
 219       logical lprn
 220       integer icant
 221       external icant
 222
 223       lprn=.false.
 224 C Compute the likelihood sum
 225       DO IB=1,NBETA(1,IPROT)
 226         elowest_all=elowest(ib,1,iprot)
 227 #ifdef DEBUG
 228           write(iout,*) "iprot",iprot," ib",ib,
 229      &      " elowest",elowest_l(ib,iprot)
 230 #endif
 231         fac=betaT(ib,1,iprot)
 232         temper=1.0d0/(fac*Rgas)
 233         sumpart=0.0d0
 234         sumlik(ib,iprot)=0.0d0
 235         do k=1,n_ene
 236           sumlikder(k,ib,iprot)=0.0d0
 237           sumqder(k,ib,iprot)=0.0d0
 238         enddo
 239         do j=1,ntot_work(iprot)
 240           i = tsil(j,iprot)
 241 #ifdef DEBUG
 242           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 243      &      " indstart",indstart(me1,iprot),
 244      &      " indend",indend(me1,iprot)
 245 #endif
 246           if (i.eq.0) goto 221
 247           jj = i2ii(i,iprot)
 248           if (jj.gt.0) then
 249 c Temperature-dependent energy
 250             etoti=0.0d0
 251             do k=1,n_ene
 252               etoti=etoti+ww(k)*escal(k,ib,1,iprot)
 253      &            *enetb(jj,k,iprot)
 254             enddo
 255             deltei=etoti-elowest(ib,1,iprot)
 256 #ifdef DEBUG
 257             write (iout,*) "etoti",etoti," deltei",deltei
 258             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 259             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 260      &        k=1,n_ene)
 261 #endif
 262             aux=entfac(i,iprot)+fac*deltei
 263 #ifdef DEBUG
 264             write (iout,'(2i5,f5.2,7(a,e15.5))')
 265      &        i,ib,fac," e_total",etoti,
 266      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 267      &        " e_lowest",elowest(ii,ib,iprot)," aux",aux
 268 #endif
 269             sumlik(ib,iprot)=sumlik(ib,iprot)+aux*Ptab(jj,ib,iprot)
 270             do k=1,n_ene
 271               sumlikder(k,ib,iprot)=sumlikder(k,ib,iprot)+
 272      &        +enetb(jj,k,iprot)*escal(k,ib,1,iprot)*Ptab(jj,ib,iprot)
 273             enddo
 274             k=0
 275             do ik=1,ntyp
 276               do jk=1,ik
 277                 k=k+1
 278                 sumlikeps(k,ib,iprot)=sumlikeps(k,ib,iprot)+
 279      &              (ftune_epsprim(eps(ik,jk))*
 280      &              eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))
 281      &              *Ptab(jj,ib,iprot)
 282               enddo
 283             enddo
 284             k=0
 285             do ik=1,ntyp
 286               do jk=1,ntyp
 287                 if (mask_tor(jk,ik).gt.0) then
 288                 k=k+1
 289                 entorjk=entor(k,jj,iprot)
 290                 sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot)
 291      &                +Ptab(jj,ib,iprot)*entorjk*escal(13,ib,1,iprot)
 292                 endif
 293               enddo
 294             enddo
 295             if (aux.le.150.0) then
 296               aux=dexp(-aux)
 297               sumpart=sumpart+aux
 298               do k=1,n_ene
 299                 sumqder(k,ib,iprot)=sumqder(k,ib,iprot)
 300      &            +aux*enetb(jj,k,iprot)*escal(k,ib,1,iprot)
 301               enddo
 302               k=0
 303               do ik=1,ntyp
 304                 do jk=1,ik
 305                   k=k+1
 306                   sumqeps(k,ib,iprot)=sumqeps(k,ib,iprot)+
 307      &                (ftune_epsprim(eps(ik,jk))*
 308      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))*aux
 309                 enddo
 310               enddo
 311               k=0
 312               do ik=1,ntyp
 313                 do jk=1,ntyp
 314                   if (mask_tor(jk,ik).gt.0) then
 315                   k=k+1
 316                   entorjk=entor(k,jj,iprot)
 317                   sumqtor(k,ib,iprot)=
 318      &             sumqtor(k,ib,iprot)+aux*entorjk*escal(13,ib,1,iprot)
 319                   endif
 320                 enddo
 321               enddo
 322             endif
 323           endif
 324   221     continue
 325         enddo ! j
 326         efree(ib,1,iprot)=sumpart
 327 #ifdef DEBUG
 328           write (iout,*) " ib",ib," iprot",iprot,
 329      &     " sumlik",sumlik(ib,iprot),
 330      &     " sumq",efree(ib,1,iprot)
 331 #endif
 332       ENDDO ! ib
 333 C Calculation of heat capacity
 334       DO IB=1,NBETA(2,IPROT)
 335         elowest_all=elowest(ib,2,iprot)
 336 #ifdef DEBUG
 337         write(iout,*) "iprot",iprot," ib",ib,
 338      &      " elowest",elowest(ib,2,iprot)
 339         write(iout,*) "escalbis",(escalbis(k,ib,2,iprot),k=1,n_ene)
 340 #endif
 341         fac=betaT(ib,2,iprot)
 342         temper=1.0d0/(fac*Rgas)
 343         sumpart=0.0d0
 344         emeani=0.0d0
 345         esquarei=0.0d0
 346         ebisi=0.0d0
 347 #ifdef DEBUG
 348         write (iout,*) "Before sum, ib=",ib
 349         write (iout,*) "efree",efree(ib,2,iprot),sumpart,
 350      &     " emean",emean_ftot(ib,iprot),
 351      &     " esquare",esquare_ftot(ib,iprot),
 352      &     " ebis",ebis_ftot(ib,iprot),
 353      &     " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
 354 #endif
 355         do j=1,ntot_work(iprot)
 356           i = tsil(j,iprot)
 357 #ifdef DEBUG
 358           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 359      &      " indstart",indstart(me1,iprot),
 360      &      " indend",indend(me1,iprot)
 361 #endif
 362           if (i.eq.0) cycle
 363           jj = i2ii(i,iprot)
 364 #ifdef DEBUG
 365           if (jj.gt.0) then
 366             write (iout,*) "ib",ib," j",j,
 367             write (iout,*) "nu",(nu(k,jj,iprot),k=1,
 368      &        natconstr(iprot))
 369           endif
 370 #endif
 371           if (jj.gt.0) then
 372 c Temperature-dependent energy
 373             etoti=0.0d0
 374             etoti_orig=0.0d0
 375             eprim=0.0d0
 376             ebis=0.0d0
 377             do k=1,n_ene
 378               etoti=etoti+ww(k)*escal(k,ib,2,iprot)
 379      &            *enetb(jj,k,iprot)
 380               eprim=eprim+ww(k)*escalprim(k,ib,2,iprot)
 381      &            *enetb(jj,k,iprot)
 382               ebis=ebis+ww(k)*escalbis(k,ib,2,iprot)
 383      &            *enetb(jj,k,iprot)
 384             enddo
 385             deltei=etoti-elowest(ib,2,iprot)
 386 #ifdef DEBUG
 387             write (iout,*) "etoti",etoti," deltei",deltei
 388             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 389             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 390      &        k=1,n_ene)
 391 #endif
 392             aux=entfac(i,iprot)+fac*deltei
 393 #ifdef DEBUG
 394             write (iout,'(f5.2,7(a,e15.5))')
 395      &        fac," e_total",etoti,
 396      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 397      &        " eprim",eprim," ebis",ebis,
 398      &        " e_lowest",elowest(ib,2,iprot)," aux",aux
 399 #endif
 400             if (aux.le.150.0) then
 401               aux=dexp(-aux)
 402               sumpart=sumpart+aux
 403               etoti=etoti-temper*eprim
 404               emeani=emeani+etoti*aux
 405               esquarei=esquarei+aux*etoti**2
 406               ebisi=ebisi+aux*ebis
 407               do k=1,n_ene
 408                 eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
 409      &           +aux*enetb(jj,k,iprot)*escal(k,ib,2,iprot)
 410 c                write (iout,*) "eave_pf",eave_pf(k,ii,ib,iprot)
 411                 eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
 412      &           +aux*enetb(jj,k,iprot)*(escal(k,ib,2,iprot)-
 413      &            temper*escalprim(k,ib,2,iprot))
 414                 eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
 415      &           +aux*enetb(jj,k,iprot)*escalbis(k,ib,2,iprot)
 416                 emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)
 417      &           +aux*enetb(jj,k,iprot)*etoti*escal(k,ib,2,iprot)
 418                 emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
 419      &           +aux*enetb(jj,k,iprot)*etoti*(escal(k,ib,2,iprot)
 420      &            -temper*escalprim(k,ib,2,iprot))
 421                 emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
 422      &           +aux*enetb(jj,k,iprot)*ebis*escal(k,ib,2,iprot)
 423                 emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
 424      &           +aux*enetb(jj,k,iprot)*etoti**2*escal(k,ib,2,iprot)
 425               enddo
 426               k=0
 427               do ik=1,ntyp
 428                 do jk=1,ik
 429                   k=k+1
 430                   enepsjk=ftune_epsprim(eps(ik,jk))*
 431      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
 432                   enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
 433      &             +aux*enepsjk
 434                   eneps_mix_ftot(k,ib,iprot)=
 435      &             eneps_mix_ftot(k,ib,iprot)+aux*enepsjk*etoti
 436                   eneps_mix_fbistot(k,ib,iprot)=
 437      &             eneps_mix_fbistot(k,ib,iprot)+aux*enepsjk*ebis
 438                   eneps_mixsq_ftot(k,ib,iprot)=
 439      &             eneps_mixsq_ftot(k,ib,iprot)+aux*enepsjk*etoti**2
 440                 enddo
 441               enddo
 442               k=0
 443               do ik=1,ntyp
 444                 do jk=1,ntyp
 445                   if (mask_tor(jk,ik).gt.0) then
 446                   k=k+1
 447                   entorjk=entor(k,jj,iprot)
 448 c                  write (iout,*) " k"," ik",ik," jk",jk,
 449 c     &             " entor",entorjk," contirb",
 450 c     &             aux*entorjk*escal(13,ib,2,iprot)
 451                   entorave_ftot(k,ib,iprot)=
 452      &             entorave_ftot(k,ib,iprot)
 453      &             +aux*entorjk*escal(13,ib,2,iprot)
 454 c                  write (iout,*) "entorave_f",
 455 c     &              entorave_ftot(k,ib,iprot)
 456                   entorave_fprimtot(k,ib,iprot)=
 457      &             entorave_fprimtot(k,ib,iprot)
 458      &             +aux*entorjk*(escal(13,ib,2,iprot)-
 459      &              temper*escalprim(13,ib,2,iprot))
 460                   entorave_fbistot(k,ib,iprot)=
 461      &             entorave_fbistot(k,ib,iprot)
 462      &             +aux*entorjk*escalbis(13,ib,2,iprot)
 463                   entor_mix_ftot(k,ib,iprot)=
 464      &             entor_mix_ftot(k,ib,iprot)
 465      &             +aux*entorjk*etoti*escal(13,ib,2,iprot)
 466                   entor_mix_fprimtot(k,ib,iprot)=
 467      &             entor_mix_fprimtot(k,ib,iprot)
 468      &             +aux*entorjk*etoti*(escal(13,ib,2,iprot)-
 469      &            temper*escalprim(13,ib,2,iprot))
 470                   entor_mix_fbistot(k,ib,iprot)=
 471      &             entor_mix_fbistot(k,ib,iprot)
 472      &             +aux*entorjk*ebis*escal(13,ib,2,iprot)
 473                   entor_mixsq_ftot(k,ib,iprot)=
 474      &             entor_mixsq_ftot(k,ib,iprot)
 475      &             +aux*entorjk*etoti**2*escal(13,ib,2,iprot)
 476                   endif
 477                 enddo
 478               enddo
 479             endif
 480           endif
 481         enddo
 482         efree(ib,2,iprot)=sumpart
 483         emean_ftot(ib,iprot)=emeani
 484         ebis_ftot(ib,iprot)=ebisi
 485         esquare_ftot(ib,iprot)=esquarei
 486 #ifdef DEBUG
 487         write (iout,*) "After sum, ib=",ib
 488         write (iout,*) "efree",efree(ib,2,iprot),sumpart,
 489      &     " emean",emean_ftot(ib,iprot),
 490      &     " esquare",esquare_ftot(ib,iprot),
 491      &     " ebis",ebis_ftot(ib,iprot),
 492      &     " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
 493 #endif
 494         ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)*temper
 495         do k=1,n_ene
 496           eave_pfbistot(k,ib,iprot)=
 497      &      eave_pfbistot(k,ib,iprot)*temper
 498           emix_pfbistot(k,ib,iprot)=
 499      &      emix_pfbistot(k,ib,iprot)*temper
 500         enddo
 501         do k=1,nntyp
 502           eneps_mix_fbistot(k,ib,iprot)=
 503      &     eneps_mix_fbistot(k,ib,iprot)*temper
 504         enddo
 505         do k=1,ntor_var
 506           entorave_fbistot(k,ib,iprot)=
 507      &     entorave_fbistot(k,ib,iprot)*temper
 508           entor_mix_fbistot(k,ib,iprot)=
 509      &     entor_mix_fbistot(k,ib,iprot)*temper
 510         enddo
 511 #ifdef DEBUG
 512         write (iout,*) "eave_pf",(eave_pftot(k,ib,iprot),
 513      &      k=1,n_ene)
 514         write (iout,*) "entorave_f",(entorave_ftot(k,ib,iprot),
 515      &      k=1,ntor_var)
 516 #endif
 517 #ifdef DEBUG
 518         write (iout,*) "ib",ib," temper",temper,
 519      &       " ebis",ebis_ftot(ib,iprot)
 520 #endif
 521       ENDDO ! ib
 522 C Calculation of conformation-dependent averages
 523       DO INAT=1,NATLIKE(IPROT)
 524       DO IB=1,NBETA(I+2,IPROT)
 525         elowest_all=elowest(ib,inat+2,iprot)
 526 #ifdef DEBUG
 527           write(iout,*) "iprot",iprot," ib",ib,
 528      &      " elowest",elowest(ib,iprot)
 529 #endif
 530         fac=betaT(ib,2,iprot)
 531         temper=1.0d0/(fac*Rgas)
 532         emeani=0.0d0
 533         do j=1,ntot_work(iprot)
 534           i = tsil(j,iprot)
 535 #ifdef DEBUG
 536           if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
 537      &      " indstart",indstart(me1,iprot),
 538      &      " indend",indend(me1,iprot)
 539 #endif
 540           if (i.eq.0) cycle
 541           jj = i2ii(i,iprot)
 542 #ifdef DEBUG
 543           if (jj.gt.0) then
 544             write (iout,*) "ib",ib," j",j,
 545             write (iout,*) "nu",(nu(k,jj,iprot),k=1,
 546      &        natconstr(iprot))
 547           endif
 548 #endif
 549           if (jj.gt.0) then
 550 c Temperature-dependent energy
 551             etoti=0.0d0
 552             do k=1,n_ene
 553               etoti=etoti+ww(k)*escal(k,ib,inat+2,iprot)
 554      &            *enetb(jj,k,iprot)
 555             enddo
 556             deltei=etoti-elowest(ib,inat+2,iprot)
 557 #ifdef DEBUG
 558             write (iout,*) "etoti",etoti," deltei",deltei
 559             write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
 560             write (iout,'(20f8.3)') (enetb(jj,k,iprot),
 561      &        k=1,n_ene)
 562 #endif
 563             aux=entfac(i,iprot)+fac*deltei
 564 #ifdef DEBUG
 565             write (iout,'(f5.2,7(a,e15.5))')
 566      &        fac," e_total",etoti,
 567      &        " eini",eini(i,iprot)," entfac",entfac(i,iprot),
 568      &        " eprim",eprim," ebis",ebis,
 569      &        " e_lowest",elowest(ib,inat+2,iprot)," aux",aux
 570 #endif
 571             if (aux.le.150.0) then
 572               aux=dexp(-aux)
 573               sumpart=sumpart+aux
 574 c 4/13/04 AL Components of the conformation-dependent averages
 575               do k=1,n_ene
 576                 eave_nat_pftot(k,ib,inat,iprot)=
 577      &           eave_nat_pftot(k,ib,inat,iprot)
 578      &           +aux*enetb(jj,k,iprot)*escal(k,ib,inat+2,iprot)
 579               enddo
 580               k=0
 581               do ik=1,ntyp
 582                 do jk=1,ik
 583                   k=k+1
 584                   enepsjk=ftune_epsprim(eps(ik,jk))*
 585      &                eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
 586                   enepsave_nat_ftot(k,ib,inat,iprot)=
 587      &             enepsave_nat_ftot(k,ib,inat,iprot)
 588      &             +aux*enepsjk
 589                 enddo
 590               enddo
 591               k=0
 592               do ik=1,ntyp
 593                 do jk=1,ntyp
 594                   if (mask_tor(jk,ik).gt.0) then
 595                   k=k+1
 596                   entorjk=entor(k,jj,iprot)
 597 c                  write (iout,*) " k"," ik",ik," jk",jk,
 598 c     &             " entor",entorjk," contirb",
 599 c     &             aux*entorjk*escal(13,ib,inat+2,iprot)
 600                   entorave_nat_ftot(k,ib,inat,iprot)=
 601      &             entorave_nat_ftot(k,ib,inat,iprot)
 602      &             +aux*entorjk*escal(13,ib,inat+2,iprot)
 603                   endif
 604                 enddo
 605               enddo
 606               do l=1,natdim(inat,iprot)
 607                 nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
 608      &            +aux*nu(l,inat,jj,iprot)
 609                 do k=1,n_ene
 610                   nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
 611      &           +aux*enetb(jj,k,iprot)*nu(k,inat,jj,iprot)*
 612      &            escal(k,ib,inat+2,iprot)
 613                 enddo
 614                 k=0
 615                 do ik=1,ntyp
 616                   do jk=1,ik
 617                     k=k+1
 618                     nuepsave_f(k,l,ib,inat,iprot)=
 619      &               nuepsave_f(k,l,ib,inat,iprot)+aux*enepsjk*
 620      &               nu(l,inat,jj,iprot)
 621                   enddo
 622                 enddo
 623                 k=0
 624                 do ik=1,ntyp
 625                   do jk=1,ntyp
 626                     if (mask_tor(jk,ik).gt.0) then
 627                       k=k+1
 628                       nutorave_f(k,l,ib,inat,iprot)=
 629      &                 nutorave_f(k,l,ib,inat,iprot)
 630      &                 +aux*entorjk*escal(13,ib,inat+2,iprot)*
 631      &                  nu(l,inat,jj,iprot)
 632                     endif
 633                   enddo
 634                 enddo
 635               enddo
 636             endif
 637           endif
 638 #ifdef DEBUG
 639           write (iout,*) "iprot",iprot," ib",ib,
 640           write (iout,*) "nuave"
 641           write (iout,'(20f10.5)') (nuave(k,ib,iprot),
 642      &            k=1,natconstr(iprot))
 643 #endif
 644 #ifdef DEBUG
 645           write (iout,*) "inat",inat,
 646      &      " efree_nat",efree_nat(ib,inat,iprot)
 647 #endif
 648         enddo
 649 #ifdef DEBUG
 650           write (iout,*) "iprot",iprot," ib",ib
 651           write (iout,*) "nuave0"
 652           write (iout,'(20f10.5)') (nuave(k,ib,inat,iprot),
 653      &            k=1,natdim(inat,iprot))
 654 #endif
 655       ENDDO ! ib
 656       ENDDO ! INAT
 657       return
 658       end
 659 c-------------------------------------------------------------------------------
 660       subroutine ave_sum(iprot)
 661       implicit none
 662       include "DIMENSIONS"
 663       include "DIMENSIONS.ZSCOPT"
 664 #ifdef MPI
 665       include "mpif.h"
 666       integer IERROR,ErrCode
 667       include "COMMON.MPI"
 668 #endif
 669       integer n,nf,What
 670       include "COMMON.WEIGHTS"
 671       include "COMMON.WEIGHTDER"
 672       include "COMMON.ENERGIES"
 673       include "COMMON.IOUNITS"
 674       include "COMMON.VMCPAR"
 675       include "COMMON.NAMES"
 676       include "COMMON.INTERACT"
 677       include "COMMON.TIME1"
 678       include "COMMON.CHAIN"
 679       include "COMMON.PROTFILES"
 680       include "COMMON.VAR"
 681       include "COMMON.COMPAR"
 682 C Define local variables
 683       integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,inat
 684       integer ipass_conf,istart_conf,iend_conf
 685       double precision energia(0:max_ene)
 686       double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk,
 687      &  eave_pft(max_ene),emix_pft(max_ene),eave_pfprimt(max_ene),
 688      &  eave_pfbist(max_ene),emix_pfprimt(max_ene),emix_pfbist(max_ene),
 689      &  esquare_ft,efree_t,
 690      &  emixsq_pft(max_ene),eneps_mixt(nntyp),eneps_meant(nntyp),
 691      &  enepsave_ft(nntyp),eneps_mix_ft(nntyp),entorave_ft(maxtor_var),
 692      &  entor_mix_ft(maxtor_var),
 693      &  eneps_mix_fbist(nntyp),entorave_fbist(maxtor_var),
 694      &  entorave_fprimt(maxtor_var),entor_mix_fprimt(maxtor_var),
 695      &  entor_mix_fbist(maxtor_var),eneps_mixsq_ft(nntyp),
 696      &  entor_mixsq_ft(maxtor_var),emean_ft,ebis_ft,nuave_t(maxdimnat),
 697      &  nu_pft(max_ene,maxdimnat),
 698      &  nuepsave_ft(nntyp,maxdimnat),
 699      &  nutorave_ft(maxtor_var,maxdimnat),
 700      &  sumlik_t,sumlikder_t(max_ene),sumlikeps_t(nntyp),
 701      &  sumliktor_t(maxtor_var),
 702      &  sumq_t,sumqder_t(max_ene),sumqeps_t(nntyp),sumqtor_t(maxtor_var)
 703
 704       double precision aux,fac,ef1,ef2,em1,em2,var1,var2,efree_tot,facF,
 705      &  elowest_all
 706       double precision gnmr,tcpu_ini,tcpu
 707       logical lprn
 708       integer icant
 709       external icant
 710
 711       lprn=.true.
 712 c Maximum likelihood contribution
 713       DO IB=1,NBETA(1,IPROT)
 714         fac=betaT(ib,1,iprot)
 715 #ifdef MPI
 716 #ifdef DEBUG
 717         write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
 718         write (iout,*) "iprot",iprot," ib",ib
 719         write (iout,*) "sumlik",sumlik(ib,iprot),
 720      &     " sumq",efree(ib,1,iprot)
 721         call flush(iout)
 722 #endif
 723         call MPI_Reduce( sumlik(ib,iprot), sumlik_t,1,
 724      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 725         call MPI_Reduce( efree(ib,1,iprot), sumq_t,1,
 726      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 727 #ifdef DEBUG
 728         write (iout,*) "sumlik",sumlik(ib,iprot)," sumlik_t",sumlik_t
 729         write (iout,*) "sumq",efree(ib,1,iprot)," sumq_t",sumq_t
 730         call flush(iout)
 731 #endif
 732         call MPI_Reduce(sumlikder(1,ib,iprot),sumlikder_t(1),
 733      &      n_ene,
 734      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 735         call MPI_Reduce(sumlikeps(1,ib,iprot),
 736      &      sumlikeps_t(1), nntyp,
 737      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 738         call MPI_Reduce(sumliktor(1,ib,iprot),
 739      &      sumliktor_t(1), ntor_var,
 740      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 741         call MPI_Reduce(sumqder(1,ib,iprot),sumqder_t(1),n_ene,
 742      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 743         call MPI_Reduce(sumqeps(1,ib,iprot),sumqeps_t(1), nntyp,
 744      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 745         call MPI_Reduce(sumqtor(1,ib,iprot),sumqtor_t(1), ntor_var,
 746      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 747 #ifdef DEBUG
 748         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
 749         call flush(iout)
 750 #endif
 751         if (me1.eq.master) then
 752 #ifdef DEBUG
 753           write (iout,*) "ib",ib,
 754      &      " elowest",elowest(ib,1,iprot),
 755      &      " sumlik",sumlik_t," qsum",sumq_t," fac",fac
 756 #endif
 757           do k=1,n_ene
 758             sumlikder(k,ib,iprot)=fac*(
 759      &          sumlikder_t(k)
 760      &          -sumqder_t(k)/sumq_t)
 761           enddo
 762           do k=1,nntyp
 763             sumlikeps(k,ib,iprot)=fac*(
 764      &          sumlikeps_t(k)-sumqeps_t(k)/sumq_t)
 765           enddo
 766 c          write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
 767 c     &          eave_pftot(19,ib,iprot)
 768           do k=1,ntor_var
 769             sumliktor(k,ib,iprot)=fac*(
 770      &          sumliktor_t(k)-sumqtor_t(k)/sumq_t)
 771           enddo
 772           sumlik(ib,iprot)=sumlik_t+dlog(sumq_t)
 773           efree(ib,1,iprot)=sumq_t
 774 #ifdef DEBUG
 775           write (iout,*) "ib",ib," iprot",iprot," final sumlik",
 776      &       sumlik(ib,iprot)," sumq",efree(ib,1,iprot)
 777 #endif
 778         endif
 779 #else
 780         do k=1,n_ene
 781           sumlikder(k,ib,iprot)=fac*(
 782      &          sumlikder(k,ib,iprot)
 783      &         -sumqder(k,ib,iprot)/efree(ib,1,iprot))
 784         enddo
 785         do k=1,nntyp
 786           sumlikeps(k,ib,iprot)=fac*(
 787      &        sumlikeps(k,ib,iprot)
 788      &       -sumqeps(k,ib,iprot)/efree(ib,1,iprot))
 789         enddo
 790         do k=1,ntor_var
 791           sumliktor(k,ib,iprot)=fac*(
 792      &      sumliktor(k,ib,iprot)
 793      &     -sumqtor(k,ib,iprot)/efree(ib,1,iprot))
 794         enddo
 795         sumlik(ib,iprot)=sumlik(ib,iprot)+dlog(efree(ib,1,iprot)
 796 c     &      -elowest(nbeta(iprot)+ib,iprot)*fac
 797 #endif
 798         ii=0
 799         do k=1,n_ene
 800           if (imask(k).gt.0) then
 801             ii=ii+1
 802             sumlikder(ii,ib,iprot)=sumlikder(k,ib,iprot)
 803           endif
 804         enddo
 805       ENDDO ! ib
 806 c Heat capacity and averages
 807       DO IB=1,NBETA(2,IPROT)
 808         fac=betaT(ib,2,iprot)
 809 #ifdef MPI
 810 #ifdef DEBUG
 811         write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
 812         write (iout,*) "iprot",iprot," ib",ib
 813         call flush(iout)
 814 #endif
 815         call MPI_Reduce( efree(ib,2,iprot), efree_t,1,
 816      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 817 #ifdef DEBUG
 818         write (iout,*) "efree",efree(iprot
 819         write (iout,*) "efree_t",efree_t
 820         call flush(iout)
 821 #endif
 822         call MPI_Reduce(emean_ftot(ib,iprot),emean_ft,1,
 823      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 824         call MPI_Reduce(ebis_ftot(ib,iprot),ebis_ft,1,
 825      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 826         call MPI_Reduce(esquare_ftot(ib,iprot),esquare_ft,1,
 827      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 828         call MPI_Reduce(eave_pftot(1,ib,iprot),eave_pft,
 829      &    n_ene,
 830      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 831         call MPI_Reduce(eave_pfprimtot(1,ib,iprot),
 832      &    eave_pfprimt,n_ene,
 833      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 834         call MPI_Reduce(eave_pfbistot(1,ib,iprot),
 835      &    eave_pfbist,n_ene,
 836      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 837 c          write (iout,*) "eave_pf",(eave_pf(k,iprot),
 838 c     &      k=1,n_ene)
 839 c          write (iout,*) "eave_pft",(eave_pft(k),k=1,n_ene)
 840         call MPI_Reduce(emix_pftot(1,ib,iprot),emix_pft,
 841      &    n_ene,
 842      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 843         call MPI_Reduce(emix_pfprimtot(1,ib,iprot),
 844      &    emix_pfprimt,n_ene,
 845      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 846         call MPI_Reduce(emix_pfbistot(1,ib,iprot),
 847      &    emix_pfbist,n_ene,
 848      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 849         call MPI_Reduce(emixsq_pftot(1,ib,iprot),
 850      &    emixsq_pft, n_ene,
 851      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 852         call MPI_Reduce(enepsave_ftot(1,ib,iprot),
 853      &    enepsave_ft, nntyp,
 854      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 855         call MPI_Reduce(eneps_mix_ftot(1,ib,iprot),
 856      &    eneps_mix_ft,nntyp,
 857      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 858         call MPI_Reduce(eneps_mix_fbistot(1,ib,iprot),
 859      &    eneps_mix_fbist,nntyp,
 860      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 861         call MPI_Reduce(eneps_mixsq_ftot(1,ib,iprot),
 862      &    eneps_mixsq_ft,nntyp,
 863      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 864 c          write (iout,*) "enepsave_f",(enepsave_f(k,iprot),
 865 c     &      k=1,nntyp)
 866 c          write (iout,*) "enepsave_ft",(enepsave_ft(k),
 867 c     &      k=1,nntyp)
 868         call MPI_Reduce(entorave_ftot(1,ib,iprot),
 869      &    entorave_ft, ntor_var,
 870      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 871         call MPI_Reduce(entorave_fprimtot(1,ib,iprot),
 872      &    entorave_fprimt, ntor_var,
 873      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 874         call MPI_Reduce(entorave_fbistot(1,ib,iprot),
 875      &    entorave_fbist, ntor_var,
 876      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 877         call MPI_Reduce(entor_mix_ftot(1,ib,iprot),
 878      &    entor_mix_ft,ntor_var,
 879      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 880         call MPI_Reduce(entor_mix_fprimtot(1,ib,iprot),
 881      &    entor_mix_fprimt,ntor_var,
 882      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 883         call MPI_Reduce(entor_mix_fbistot(1,ib,iprot),
 884      &   entor_mix_fbist,ntor_var,
 885      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 886         call MPI_Reduce(entor_mixsq_ftot(1,ib,iprot),
 887      &    entor_mixsq_ft,ntor_var,
 888      &    MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 889 c          write (iout,*) "entorave_f",(entorave_f(k,iprot),
 890 c     &      k=1,ntor_var)
 891 c          write (iout,*) "entorave_ft",(entorave_ft(k),
 892 c     &      k=1,ntor_var)
 893 #ifdef DEBUG
 894         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
 895         call flush(iout)
 896 #endif
 897         if (me1.eq.master) then
 898           elowest_all=elowest(ib,2,iprot)
 899 #ifdef DEBUG
 900           write (iout,*) "ib",ib,
 901      &      " elowest",elowest(ib,iprot),
 902      &      "efree",efree_t," fac",fac," facF",facF,
 903      &      " efree_tot",efree_tot
 904 #endif
 905           do k=1,n_ene
 906             eave_pftot(k,ib,iprot)=eave_pft(k)/efree_t
 907             eave_pfprimtot(k,ib,iprot)=eave_pfprimt(k)/efree_t
 908             eave_pfbistot(k,ib,iprot)=eave_pfbist(k)/efree_t
 909             emix_pftot(k,ib,iprot)=emix_pft(k)/efree_t
 910             emix_pfprimtot(k,ib,iprot)=emix_pfprimt(k)/efree_t
 911             emix_pfbistot(k,ib,iprot)=emix_pfbist(k)/efree_t
 912             emixsq_pftot(k,ib,iprot)=emixsq_pft(k)/efree_t
 913           enddo
 914           do k=1,nntyp
 915             enepsave_ftot(k,ib,iprot)=enepsave_ft(k)/efree_t
 916             eneps_mix_ftot(k,ib,iprot)=eneps_mix_ft(k)/efree_t
 917             eneps_mix_fbistot(k,ib,iprot)=eneps_mix_fbist(k)/efree_t
 918             eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ft(k)/efree_t
 919           enddo
 920 c          write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
 921 c     &          eave_pftot(19,ib,iprot)
 922           do k=1,ntor_var
 923             entorave_ftot(k,ib,iprot)=entorave_ft(k)/
 924      &        efree_t
 925 c            write (iout,*) "iprot",iprot," ib",ib," k",k,
 926 c     &         " entorave_ftot", entorave_ftot(k,ib,iprot)
 927             entorave_fprimtot(k,ib,iprot)=entorave_fprimt(k)/efree_t
 928             entorave_fbistot(k,ib,iprot)=entorave_fbist(k)/efree_t
 929             entor_mix_ftot(k,ib,iprot)=entor_mix_ft(k)/efree_t
 930             entor_mix_fprimtot(k,ib,iprot)=entor_mix_fprimt(k)/efree_t
 931             entor_mix_fbistot(k,ib,iprot)=entor_mix_fbist(k)/efree_t
 932             entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ft(k)/efree_t
 933           enddo
 934           emean_ftot(ib,iprot)=emean_ft/efree_t
 935           ebis_ftot(ib,iprot)=ebis_ft/efree_t
 936           esquare_ftot(ib,iprot)=esquare_ft/efree_t
 937           efree(ib,2,iprot)=-dlog(efree_t)/fac+elowest(ib,2,iprot)
 938         endif
 939 #else
 940         do k=1,n_ene
 941           eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
 942      &      /efree(ib,2,iprot)
 943           eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
 944      &      /efree(ib,2,iprot)
 945           eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
 946      &      /efree(ib,2,iprot)
 947           emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)/efree(ib,iprot)
 948           emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
 949      &      /efree(ib,2,iprot)
 950           emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
 951      &      /efree(ib,2,iprot)
 952           emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
 953      &      /efree(ib,2,iprot)
 954         enddo
 955         do k=1,nntyp
 956           enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
 957      &      /efree(ib,2,iprot)
 958           eneps_mix_ftot(k,ib,iprot)=eneps_mix_ftot(k,ib,iprot)
 959      &      /efree(ib,2,iprot)
 960           eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ftot(k,ib,iprot)
 961      &      /efree(ib,2,iprot)
 962         enddo
 963         do k=1,ntor_var
 964           entorave_ftot(k,ib,iprot)=entorave_f(k,ib,iprot)
 965      &      /efree(ib,2,iprot)
 966           entor_mix_ftot(k,ib,iprot)=entor_mix_ftot(k,ib,iprot)
 967      &      /efree(ib,2,iprot)
 968           entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ftot(k,ib,iprot)
 969      &      /efree(ib,2,iprot)
 970         enddo
 971         emean_ftot(ib,iprot)=emean_ftot(ib,iprot)/efree(ib,2,iprot)
 972         ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)/efree(ib,2,iprot)
 973         esquare_ftot(ib,iprot)=esquare_ftot(ib,iprot)/efree(ib,2,iprot)
 974 #endif
 975       ENDDO ! IB
 976 c 4/13/04 AL Components of the correlation coefficients and their derivatives
 977       DO INAT=1,NATLIKE(IPROT)
 978         DO IB=1,NBETA(INAT+2,IPROT)
 979 #ifdef MPI
 980           call MPI_Reduce( efree(ib,inat+2,iprot), efree_t,1,
 981      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 982           call MPI_Reduce( eave_nat_pftot(1, ib,inat,iprot), eave_pft,
 983      &      n_ene,
 984      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 985           call MPI_Reduce( enepsave_nat_ftot(1, ib,inat,iprot),
 986      &      enepsave_ft, nntyp,
 987      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 988           call MPI_Reduce(entorave_nat_ftot(1,ib,inat,iprot),
 989      &      entorave_ft, ntor_var,
 990      &      MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
 991           call MPI_Reduce(nuave(1,ib,inat,iprot),nuave_t(1),
 992      &      natdim(inat,iprot), MPI_DOUBLE_PRECISION, MPI_SUM,
 993      &      Master, Comm1, IERROR)
 994 #ifdef DEBUG
 995           write (iout,*) "After REDUCE nuave",iprot,ib
 996           write (iout,'(20f10.5)')
 997      &      (nuave(l,ib,iprot),l=1,natconstr(iprot))
 998 #endif
 999           call MPI_Reduce(nu_pf(1,1,ib,inat,iprot),nu_pft(1,1),
1000      &      max_ene*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1001      &      MPI_SUM, Master, Comm1, IERROR)
1002           call MPI_Reduce(nuepsave_f(1,1,ib,inat,iprot),
1003      &      nuepsave_ft(1,1),
1004      &      nntyp*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1005      &      MPI_SUM, Master, Comm1, IERROR)
1006           call MPI_Reduce(nutorave_f(1,1,ib,inat,iprot),
1007      &      nutorave_ft(1,1),
1008      &      maxtor_var*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1009      &      MPI_SUM, Master, Comm1, IERROR)
1010 #ifdef DEBUG
1011         write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
1012         call flush(iout)
1013 #endif
1014           if (me1.eq.master) then
1015 #ifdef DEBUG
1016             write (iout,*) "ib",ib,
1017      &      " elowest",elowest(ib,iprot),
1018      &      "efree",efree_t," fac",fac,
1019      &      " efree_tot",efree_tot
1020 #endif
1021             do l=1,natdim(inat,iprot)
1022               do k=1,n_ene
1023                 nu_pf(k,l,ib,inat,iprot)=nu_pft(k,l)/efree_t
1024               enddo
1025               do k=1,nntyp
1026                 nuepsave_f(k,l,ib,inat,iprot)=nuepsave_ft(k,l)/efree_t
1027               enddo
1028               do k=1,ntor_var
1029                 nutorave_f(k,l,ib,inat,iprot)=nutorave_ft(k,l)/efree_t
1030               enddo
1031               nuave(l,ib,inat,iprot)=nuave_t(l)/efree_t
1032             enddo
1033             do k=1,n_ene
1034               eave_nat_pftot(k,ib,inat,iprot)=eave_pft(k)/efree_t
1035             enddo
1036             do k=1,nntyp
1037               enepsave_nat_ftot(k,ib,inat,iprot)=enepsave_ft(k)/efree_t
1038             enddo
1039             do k=1,ntor_var
1040               entorave_nat_ftot(k,ib,inat,iprot)=entorave_ft(k)/efree_t
1041             enddo
1042           endif
1043 #else
1044           do l=1,natdim(inat,iprot)
1045             do k=1,n_ene
1046               nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
1047      &         /efree(ib,inat+2,iprot)
1048             enddo
1049             do k=1,nntyp
1050              nuepsave_f(k,l,ib,inat,iprot)=nuepsave_f(k,l,ib,inat,iprot)
1051      &            /efree(ib,inat+2,iprot)
1052             enddo
1053             do k=1,ntor_var
1054               nutorave_ftot(k,l,ib,inat,iprot)=
1055      &          nutorave_ftot(k,l,ib,inat,iprot)
1056      &          /efree(ib,inat+2,iprot)
1057             enddo
1058             nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
1059      &        /efree(ib,inat+2,iprot)
1060           enddo ! l
1061           do k=1,n_ene
1062             eave_nat_pftot(k,ib,inat,iprot)=
1063      &       eave_nat_pftot(k,ib,inat,iprot)
1064      &       /efree(ib,inat+2,iprot)
1065           enddo
1066           do k=1,nntyp
1067             enepsave_nat_ftot(k,ib,inat,iprot)=
1068      &         enepsave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1069           enddo
1070           do k=1,ntor_var
1071             enetorave_nat_ftot(k,ib,inat,iprot)=
1072      &        enetorave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1073           enddo
1074 #endif
1075         ENDDO ! IB
1076 #ifdef DEBUG
1077         write (iout,*) "ib",ib," efree_tot",efree_tot
1078 #endif
1079       ENDDO ! INAT
1080
1081       return
1082       end