source/maxlik/src_MD_T_maxlik-NEWCORR.safe/thermal.F.safe

   1       subroutine thermal(iprot)
   2 c Calculate the average energy, average Q, specific heat, and Q dispersion
   3 c curves in temperature for protein IPROT.
   4       implicit none
   5       include "DIMENSIONS"
   6       include "DIMENSIONS.ZSCOPT"
   7 #ifdef MPI
   8       include "mpif.h"
   9       integer IERROR,ErrCode
  10       include "COMMON.MPI"
  11 #endif
  12       integer n,nf,What
  13       include "COMMON.THERMAL"
  14       include "COMMON.WEIGHTS"
  15       include "COMMON.WEIGHTDER"
  16       include "COMMON.CLASSES"
  17       include "COMMON.ENERGIES"
  18       include "COMMON.IOUNITS"
  19       include "COMMON.VMCPAR"
  20       include "COMMON.NAMES"
  21       include "COMMON.INTERACT"
  22       include "COMMON.TIME1"
  23       include "COMMON.CHAIN"
  24       include "COMMON.PROTFILES"
  25       include "COMMON.COMPAR"
  26 C Define local variables
  27       integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,l1,l2,ll1,ll2,ibatch
  28       integer ipass_conf,istart_conf,iend_conf
  29       logical lprn
  30       integer icant
  31       external icant
  32       double precision elowest_all,elowest_ent_all,entbase_all
  33
  34 c      write (iout,*) "Thermal"
  35 c      call flush(iout)
  36       lprn=.true.
  37       elowest_all=1.0d20
  38       elowest_ent_all=1.0d20
  39       do ibatch=1,natlike(iprot)+2
  40         do ib=1,nbeta(ibatch,iprot)
  41 c          write (iout,*) "ibatch",ibatch," ib",ib
  42           call flush(iout)
  43           if (elowest(ib,ibatch,iprot).lt.elowest_all)
  44      &      elowest_all=elowest(ib,ibatch,iprot)
  45 c          write (iout,*) "ibatch",ibatch," ib",ib," after if"
  46           call flush(iout)
  47           if (elowest_ent(ib,ibatch,iprot).lt.elowest_ent_all) then
  48             elowest_ent_all = elowest_ent(ib,ibatch,iprot)
  49           endif
  50         enddo
  51       enddo
  52 c      write (iout,*) "THERMAL: elowest_all",elowest_all
  53       call flush(iout)
  54 c
  55 c Zero out tables.
  56 c
  57       do i=0,nGridT
  58         sumW(i,iprot)=0.0d0
  59         sumE(i,iprot)=0.0d0
  60         sumEbis(i,iprot)=0.0d0
  61         sumEsq(i,iprot)=0.0d0
  62         sumEEE(i,iprot)=0.0d0
  63         sumQ(i,iprot)=0.0d0
  64         sumRMS(i,iprot)=0.0d0
  65         sumRGY(i,iprot)=0.0d0
  66         sumQsq(i,iprot)=0.0d0
  67         sumEQ(i,iprot)=0.0d0
  68         sumEEQ(i,iprot)=0.0d0
  69       enddo
  70 c
  71 c Calculate the contributions to thermodynamic quantities from all processors.
  72 c
  73 #ifdef DEBUG
  74       write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot)
  75 #endif
  76       IF (NCHUNK_CONF(IPROT).EQ.1) THEN
  77
  78 #ifdef MPI
  79       do i=1,ntot_work(iprot)
  80         i2ii(i,iprot)=0
  81       enddo
  82       ii=0
  83       do i=indstart(me1,iprot),indend(me1,iprot)
  84         ii=ii+1
  85         i2ii(i,iprot)=ii
  86       enddo
  87 #else
  88       do i=1,ntot_work(iprot)
  89         i2ii(i,iprot)=i
  90       endif
  91 #endif
  92 #ifdef DEBUG
  93       do i=1,ntot_work(iprot)
  94         write (iout,*) "i",i," i2ii",i2ii(i,iprot)
  95       enddo
  96       call flush(iout)
  97 #endif
  98       call calc_thermal(iprot,elowest_all,*10)
  99
 100       ELSE
 101
 102       open (ientin,file=benefiles(iprot),status="old",
 103      &    form="unformatted",access="direct",recl=lenrec_ene(iprot))
 104       ipass_conf=0
 105 #ifdef MPI
 106       do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
 107         iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
 108 #else
 109       do istart_conf=1,ntot_work(iprot),maxstr_proc
 110         iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
 111 #endif
 112 c
 113 c Read the chunk of energies and derivatives off a DA scratchfile.
 114 c
 115         ipass_conf=ipass_conf+1
 116         do i=1,ntot_work(iprot)
 117           i2ii(i,iprot)=0
 118         enddo
 119         ii=0
 120         do i=istart_conf,iend_conf
 121           ii=ii+1
 122           i2ii(i,iprot)=ii
 123         enddo
 124 #ifdef DEBUG
 125         write (iout,*) "ipass_conf",ipass_conf,
 126      &    " istart_conf",istart_conf," iend_conf",iend_conf
 127         do i=1,ntot_work(iprot)
 128           write (iout,*) "i",i," i2ii",i2ii(i,iprot)
 129         enddo
 130         call flush(iout)
 131 #endif
 132         call daread_ene(iprot,istart_conf,iend_conf)
 133         call calc_thermal(iprot,elowest_all,*10)
 134       enddo
 135
 136       close(ientin)
 137       ENDIF
 138  10   continue
 139 c
 140 c Put tohether all contributions.
 141 c
 142       call sum_thermal(iprot,elowest_all,entbase_all)
 143       return
 144       end
 145 c-------------------------------------------------------------------------------
 146       subroutine calc_thermal(iprot,elowest_all,*)
 147       implicit none
 148       include "DIMENSIONS"
 149       include "DIMENSIONS.ZSCOPT"
 150 #ifdef MPI
 151       include "mpif.h"
 152       integer IERROR,ErrCode
 153       include "COMMON.MPI"
 154 #endif
 155       integer n,nf,What
 156       include "COMMON.THERMAL"
 157       include "COMMON.WEIGHTS"
 158       include "COMMON.WEIGHTDER"
 159       include "COMMON.CLASSES"
 160       include "COMMON.COMPAR"
 161       include "COMMON.ENERGIES"
 162       include "COMMON.IOUNITS"
 163       include "COMMON.VMCPAR"
 164       include "COMMON.NAMES"
 165       include "COMMON.INTERACT"
 166       include "COMMON.TIME1"
 167       include "COMMON.CHAIN"
 168       include "COMMON.FFIELD"
 169       include "COMMON.PROTFILES"
 170       include "COMMON.SBRIDGE"
 171       include "COMMON.ALLPROT"
 172 C Define local variables
 173       integer i,ii,jj,j,k,kk,l,m,iprot,ibatch,istart_conf,iend_conf
 174       double precision aux,auxf,fac,weight,etoti,eprim,ebis,denom,
 175      &  temper,quot,quotl,quotl1,kfacl,logfac,eplus,eminus,tanhT
 176       double precision elowest_all,entbase_all
 177       double precision ft(5),ftprim(5),ftbis(5),escal1(max_ene),
 178      &  escal1prim(max_ene),escal1bis(max_ene)
 179       double precision T0 /300.0d0/, kfac /2.4d0/
 180
 181       ii=natlike(iprot)-2
 182 c      write (iout,*) "ii",ii
 183 c      write (iout,*) "Thermal: the i2ii array"
 184 c      do i=1,ntot_work(iprot)
 185 c        write (iout,*) i,i2ii(i)
 186 c      enddo
 187 c      write (iout,*) "The list array"
 188 c      do i=1,ntot_work(iprot)
 189 c        write (iout,*) i,list_conf(i)
 190 c      enddo
 191 c      write (iout,*) "iprot",iprot," natlike",natlike(iprot)
 192       do i=1,ntot_work(iprot)
 193         jj = i2ii(i,iprot)
 194         if (jj.gt.0) then
 195 c          etoti=e_total(i,iprot)-elowest_all
 196           auxf=entfac(i,iprot)
 197 c          write (iout,*) "etoti",etoti+elowest_all," auxf",auxf," nu",
 198 c     &     nu(ii,jj,iprot)
 199 c          write (iout,'(2i5,20f8.2)') i,jj,(enetb(jj,kk,iprot),
 200 c     &      kk=1,n_ene)
 201           if (ii.gt.0) rmstb(i,iprot)=nu(1,ii,jj,iprot)
 202 c          write (iout,*) "i",i," jj",jj," nu",(nu(k,jj,iprot),
 203 c     &     k=1,natlike(iprot)+2)
 204           do k=0,nGridT
 205             temper=GridT0+k*deltaT
 206             fac=1.0d0/(1.987D-3*temper)
 207             quot=temper/T0
 208             if (rescale_mode.eq.0) then
 209               do l=1,5
 210                 fT(l)=1.0d0
 211                 fTprim(l)=0.0d0
 212                 fTbis(l)=0.0d0
 213               enddo
 214             else if (rescale_mode.eq.1) then
 215               quotl=1.0d0
 216               kfacl=1.0d0
 217               do l=1,5
 218                 quotl1=quotl
 219                 quotl=quotl*quot
 220                 kfacl=kfacl*kfac
 221                 denom=kfacl-1.0d0+quotl
 222                 fT(l)=kfacl/denom
 223                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 224                 ftbis(l)=l*kfacl*quotl1*
 225      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 226               enddo
 227             else if (rescale_mode.eq.2) then
 228               quotl=1.0d0
 229               do l=1,5
 230                 quotl1=quotl
 231                 quotl=quotl*quot
 232                 eplus=dexp(quotl)
 233                 eminus=dexp(-quotl)
 234                 logfac=1.0d0/dlog(eplus+eminus)
 235                 tanhT=(eplus-eminus)/(eplus+eminus)
 236                 fT(l)=1.12692801104297249644d0*logfac
 237                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 238                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 239      &          2*l*quotl1/T0*logfac*
 240      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 241      &          +ftprim(l)*tanhT)
 242               enddo
 243             else
 244               write (iout,*) "Unknown RESCALE_MODE",rescale_mode
 245               call flush(iout)
 246               return1
 247             endif
 248             do l=1,n_ene
 249               escal1(l)=1.0d0
 250               escal1prim(l)=0.0d0
 251               escal1bis(l)=0.0d0
 252             enddo
 253             escal1(3)=ft(1)
 254             escal1(4)=ft(3)
 255             escal1(5)=ft(4)
 256             escal1(6)=ft(5)
 257             escal1(7)=ft(2)
 258             escal1(8)=ft(2)
 259             escal1(9)=ft(3)
 260             escal1(10)=ft(5)
 261             escal1(13)=ft(1)
 262             escal1(14)=ft(2)
 263             escal1(19)=ft(1)
 264             escal1prim(3)=ftprim(1)
 265             escal1prim(4)=ftprim(3)
 266             escal1prim(5)=ftprim(4)
 267             escal1prim(6)=ftprim(5)
 268             escal1prim(7)=ftprim(2)
 269             escal1prim(8)=ftprim(2)
 270             escal1prim(9)=ftprim(3)
 271             escal1prim(10)=ftprim(5)
 272             escal1prim(13)=ftprim(1)
 273             escal1prim(14)=ftprim(2)
 274             escal1prim(19)=ftprim(1)
 275             escal1bis(3)=ftbis(1)
 276             escal1bis(4)=ftbis(3)
 277             escal1bis(5)=ftbis(4)
 278             escal1bis(6)=ftbis(5)
 279             escal1bis(7)=ftbis(2)
 280             escal1bis(8)=ftbis(2)
 281             escal1bis(9)=ftbis(3)
 282             escal1bis(10)=ftbis(5)
 283             escal1bis(13)=ftbis(1)
 284             escal1bis(14)=ftbis(2)
 285             escal1bis(19)=ftbis(1)
 286             etoti=0.0d0
 287             eprim=0.0d0
 288             ebis=0.0d0
 289             do l=1,n_ene
 290               etoti=etoti+ww(l)*escal1(l)*enetb(jj,l,iprot)
 291               eprim=eprim+ww(l)*escal1prim(l)*enetb(jj,l,iprot)
 292               ebis=ebis+ww(l)*escal1bis(l)*enetb(jj,l,iprot)
 293             enddo
 294 #ifdef DEBUG
 295             write (iout,'(4i5,3f10.5)') i,ii,jj,iprot,GridT0+k*deltaT,
 296      &        etoti,nu(ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)
 297 #endif
 298             etoti=etoti-elowest_all
 299             aux=fac*etoti+auxf
 300 #ifdef DEBUG
 301             write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj,
 302      &        GridT0+k*deltaT,
 303      &        etoti+elowest_all,e_total(i,iprot),
 304      &        elowest_all,fac*etoti,auxf,aux
 305 #endif
 306             if (aux.le.150.0) then
 307               etoti=etoti-temper*eprim
 308               weight=dexp(-aux)
 309               sumW(k,iprot)=sumW(k,iprot)+weight
 310               sumE(k,iprot)=sumE(k,iprot)+etoti*weight
 311               sumEbis(k,iprot)=sumEbis(k,iprot)+ebis*weight
 312               sumEsq(k,iprot)=sumEsq(k,iprot)+etoti**2*weight
 313               sumEEE(k,iprot)=sumEEE(k,iprot)+etoti**3*weight
 314               if (ii.gt.0) then
 315               sumQ(k,iprot)=sumQ(k,iprot)+nu(1,ii,jj,iprot)*weight
 316               sumRMS(k,iprot)=sumRMS(k,iprot)
 317      &       +nu(1,ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)*weight
 318               sumRGY(k,iprot)=sumRGY(k,iprot)
 319      &           +nu(1,ii+2,jj,iprot)*weight
 320               sumQsq(k,iprot)=sumQsq(k,iprot)
 321      &           +nu(1,ii,jj,iprot)**2*weight
 322               sumEQ(k,iprot)=sumEQ(k,iprot)
 323      &           +etoti*nu(1,ii,jj,iprot)*weight
 324               sumEEQ(k,iprot)=sumEEQ(k,iprot)+etoti**2*nu(1,ii,jj,iprot)
 325      &          *weight
 326               endif
 327             endif
 328           enddo
 329         endif
 330       enddo
 331 c      do k=0,nGridT
 332 c        write (iout,'(i5,7e15.5)') k,sumW(k,iprot),sumE(k,iprot),
 333 c     &   sumEsq(k,iprot),sumQ(k,iprot),sumQsq(k,iprot),sumEQ(k,iprot)
 334 c      enddo
 335       return
 336       end
 337 c-------------------------------------------------------------------------------
 338       subroutine sum_thermal(iprot,elowest_all,entbase_all)
 339       implicit none
 340       include "DIMENSIONS"
 341       include "DIMENSIONS.ZSCOPT"
 342 #ifdef MPI
 343       include "mpif.h"
 344       integer IERROR,ErrCode
 345       include "COMMON.MPI"
 346       double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT),
 347      & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT),
 348      & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT),
 349      & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT)
 350 #endif
 351       include "COMMON.THERMAL"
 352       include "COMMON.PROTNAME"
 353       include "COMMON.WEIGHTS"
 354       include "COMMON.WEIGHTDER"
 355       include "COMMON.CLASSES"
 356       include "COMMON.ENERGIES"
 357       include "COMMON.IOUNITS"
 358       include "COMMON.VMCPAR"
 359       include "COMMON.NAMES"
 360       include "COMMON.INTERACT"
 361       include "COMMON.TIME1"
 362       include "COMMON.CHAIN"
 363       include "COMMON.PROTFILES"
 364       include "COMMON.OPTIM"
 365 C Define local variables
 366       integer i,iprot
 367       integer ilen
 368       external ilen
 369       double precision elowest_all,entbase_all,RgasT,Temper
 370       character*64 nazwa
 371       character*3 liczba
 372       if (nparmset.eq.0) then
 373         nazwa=prefix(:ilen(prefix))
 374      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 375       else
 376         write (liczba,'(bz,i3.3)') myparm
 377         nazwa=prefix(:ilen(prefix))//'_par'//liczba
 378      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 379       endif
 380 #ifdef MPI
 381       call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1,
 382      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 383       call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1,
 384      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 385       call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1,
 386      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 387       call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1,
 388      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 389       call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1,
 390      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 391       call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1,
 392      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 393       call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1,
 394      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 395       call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1,
 396      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 397       call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1,
 398      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 399       call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1,
 400      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 401       call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1,
 402      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 403       call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
 404      &   scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 405      &   scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 406      &   Comm1, IERROR)
 407       if (me1.eq.master) then
 408         open (istat,file=nazwa)
 409         write (istat,'(a)') '# Thermal characteristics of folding'
 410         do i=0,NGridT
 411           Temper = GridT0+i*deltaT
 412           RgasT = Rgas*Temper
 413 c          write (iout,*) "i",i," Temper",Temper," RgasT",RgasT,
 414 c     &     " entbase_all",entbase_all
 415           sumE(i,iprot)=sumE_t(i)/sumW_t(i)
 416           sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i)
 417           sumQ(i,iprot)=sumQ_t(i)/sumW_t(i)
 418           sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i)
 419           sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i)
 420           sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i)
 421           sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i)
 422      &     -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3
 423           sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2
 424           sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i)
 425           sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)-
 426      &      sumEsq(i,iprot)*sumQ(i,iprot)-
 427      &      2*sumEQ(i,iprot)*sumE(i,iprot)+
 428      &      2*sumQ(i,iprot)*sumE(i,iprot)**2
 429 #ifdef DEBUG
 430           write (iout,*) "Temper",temper," RgasT",RgasT
 431           write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
 432      &      " sumEbis",sumEbis(i,iprot)
 433 #endif
 434           sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 435           sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 436           sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all
 437           sumE(i,iprot)=sumE(i,iprot)+elowest_all
 438           write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 439      &      sumW(i,iprot),
 440      &      sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT
 441      &      +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 442      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 443      &      sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot),
 444      &      sumEEQ(i,iprot)
 445 #ifdef DEBUG
 446           write (iout,*) "zvtot",
 447      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
 448 #endif
 449           call flush(istat)
 450         enddo
 451         close (istat)
 452       endif
 453 #else
 454       open (istat,file=nazwa)
 455       write (istat,'(/a)') '# Thermal characteristics of folding'
 456       do i=0,NGridT
 457         Temper = GridT0+i*deltaT
 458         RgasT = Rgas*Temper
 459         sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot)
 460         sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot)
 461         sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot)
 462         sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot)
 463         sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot)
 464         sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot)
 465         sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2
 466         sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot)
 467         sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)-
 468      &    sumEsq(i,iprot)*sumQ(i,iprot)-
 469      &    2*sumEQ(i,iprot)*sumE(i,iprot)+
 470      &    2*sumQ(i,iprot)*sumE(i,iprot)**2
 471         sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 472         sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 473         sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all
 474         sumE(i,iprot)=sumE(i,iprot)+elowest_all
 475         write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 476      &   sumW(i,iprot),
 477      &   sumE(i,iprot),
 478      &   (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all,
 479      &   sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 480      &   sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 481      &   sumQsq(i,iprot),
 482      &   sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)
 483       enddo
 484       close (istat)
 485 #endif
 486       return
 487       end