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