source/maxlik/src-Fmatch_safe/averages_sc.F.org

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