source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/thermal.F

   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.CONTROL"
 157       include "COMMON.THERMAL"
 158       include "COMMON.WEIGHTS"
 159       include "COMMON.WEIGHTDER"
 160       include "COMMON.CLASSES"
 161       include "COMMON.COMPAR"
 162       include "COMMON.ENERGIES"
 163       include "COMMON.IOUNITS"
 164       include "COMMON.VMCPAR"
 165       include "COMMON.NAMES"
 166       include "COMMON.INTERACT"
 167       include "COMMON.TIME1"
 168       include "COMMON.CHAIN"
 169       include "COMMON.FFIELD"
 170       include "COMMON.PROTFILES"
 171       include "COMMON.SBRIDGE"
 172       include "COMMON.ALLPROT"
 173 C Define local variables
 174       integer i,ii,jj,j,k,kk,l,m,iprot,ibatch,istart_conf,iend_conf,ib
 175       double precision aux,auxf,fac,weight,etoti,eprim,ebis,denom,
 176      &  tt,temper,quot,quotl,quotl1,kfacl,logfac,eplus,eminus,tanhT
 177       double precision elowest_all,entbase_all
 178       double precision ft(5),ftprim(5),ftbis(5),escal1(max_ene),
 179      &  escal1prim(max_ene),escal1bis(max_ene)
 180       double precision T0 /300.0d0/, kfac /2.4d0/
 181
 182       ii=natlike(iprot)-2
 183 c      write (iout,*) "ii",ii
 184 c      write (iout,*) "Thermal: the i2ii array"
 185 c      do i=1,ntot_work(iprot)
 186 c        write (iout,*) i,i2ii(i)
 187 c      enddo
 188 c      write (iout,*) "The list array"
 189 c      do i=1,ntot_work(iprot)
 190 c        write (iout,*) i,list_conf(i)
 191 c      enddo
 192 c      write (iout,*) "iprot",iprot," natlike",natlike(iprot)
 193       do i=1,ntot_work(iprot)
 194         jj = i2ii(i,iprot)
 195         if (jj.gt.0) then
 196 c          etoti=e_total(i,iprot)-elowest_all
 197           auxf=entfac(i,iprot)
 198 c          write (iout,*) "etoti",etoti+elowest_all," auxf",auxf," nu",
 199 c     &     nu(ii,jj,iprot)
 200 c          write (iout,'(2i5,20f8.2)') i,jj,(enetb(jj,kk,iprot),
 201 c     &      kk=1,n_ene)
 202 c          if (ii.gt.0) rmstb(i,iprot)=nu(1,ii,jj,iprot)
 203 c          write (iout,*) "i",i," jj",jj," nu",(nu(k,jj,iprot),
 204 c     &     k=1,natlike(iprot)+2)
 205 c          write (iout,*) "i",i," rmstb",rmstb(i,iprot)
 206           do k=0,nGridT
 207             temper=GridT0+k*deltaT
 208             fac=1.0d0/(1.987D-3*temper)
 209             quot=temper/T0
 210             if (rescale_mode.eq.0) then
 211               do l=1,5
 212                 fT(l)=1.0d0
 213                 fTprim(l)=0.0d0
 214                 fTbis(l)=0.0d0
 215               enddo
 216             else if (rescale_mode.eq.1) then
 217               quotl=1.0d0
 218               kfacl=1.0d0
 219               do l=1,5
 220                 quotl1=quotl
 221                 quotl=quotl*quot
 222                 kfacl=kfacl*kfac
 223                 denom=kfacl-1.0d0+quotl
 224                 fT(l)=kfacl/denom
 225                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 226                 ftbis(l)=l*kfacl*quotl1*
 227      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 228               enddo
 229             else if (rescale_mode.eq.2) then
 230               quotl=1.0d0
 231               do l=1,5
 232                 quotl1=quotl
 233                 quotl=quotl*quot
 234                 eplus=dexp(quotl)
 235                 eminus=dexp(-quotl)
 236                 logfac=1.0d0/dlog(eplus+eminus)
 237                 tanhT=(eplus-eminus)/(eplus+eminus)
 238                 fT(l)=1.12692801104297249644d0*logfac
 239                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 240                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 241      &          2*l*quotl1/T0*logfac*
 242      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 243      &          +ftprim(l)*tanhT)
 244               enddo
 245             else
 246               write (iout,*) "Unknown RESCALE_MODE",rescale_mode
 247               call flush(iout)
 248               return1
 249             endif
 250             do l=1,n_ene
 251               escal1(l)=1.0d0
 252               escal1prim(l)=0.0d0
 253               escal1bis(l)=0.0d0
 254             enddo
 255             if (shield_mode.gt.0) then
 256               escal1(1)=ft(1)
 257               escal1(2)=ft(1)
 258               escal1(16)=ft(1)
 259               escal1prim(1)=ftprim(1)
 260               escal1prim(2)=ftprim(1)
 261               escal1prim(16)=ftprim(1)
 262               escal1bis(1)=ftbis(1)
 263               escal1bis(2)=ftbis(1)
 264               escal1bis(16)=ftbis(1)
 265             endif
 266             escal1(3)=ft(1)
 267             escal1(4)=ft(3)
 268             escal1(5)=ft(4)
 269             escal1(6)=ft(5)
 270             escal1(7)=ft(2)
 271             escal1(8)=ft(2)
 272             escal1(9)=ft(3)
 273             escal1(10)=ft(5)
 274             escal1(13)=ft(1)
 275             escal1(14)=ft(2)
 276             escal1(19)=ft(1)
 277             escal1prim(3)=ftprim(1)
 278             escal1prim(4)=ftprim(3)
 279             escal1prim(5)=ftprim(4)
 280             escal1prim(6)=ftprim(5)
 281             escal1prim(7)=ftprim(2)
 282             escal1prim(8)=ftprim(2)
 283             escal1prim(9)=ftprim(3)
 284             escal1prim(10)=ftprim(5)
 285             escal1prim(13)=ftprim(1)
 286             escal1prim(14)=ftprim(2)
 287             escal1prim(19)=ftprim(1)
 288             escal1bis(3)=ftbis(1)
 289             escal1bis(4)=ftbis(3)
 290             escal1bis(5)=ftbis(4)
 291             escal1bis(6)=ftbis(5)
 292             escal1bis(7)=ftbis(2)
 293             escal1bis(8)=ftbis(2)
 294             escal1bis(9)=ftbis(3)
 295             escal1bis(10)=ftbis(5)
 296             escal1bis(13)=ftbis(1)
 297             escal1bis(14)=ftbis(2)
 298             escal1bis(19)=ftbis(1)
 299             etoti=0.0d0
 300             eprim=0.0d0
 301             ebis=0.0d0
 302             do l=1,n_ene
 303               etoti=etoti+ww(l)*escal1(l)*enetb(jj,l,iprot)
 304               eprim=eprim+ww(l)*escal1prim(l)*enetb(jj,l,iprot)
 305               ebis=ebis+ww(l)*escal1bis(l)*enetb(jj,l,iprot)
 306             enddo
 307 #ifdef DEBUG
 308             write (iout,'(4i5,3f10.5)') i,ii,jj,iprot,GridT0+k*deltaT,
 309      &        etoti,nu(ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)
 310 #endif
 311             etoti=etoti-elowest_all
 312             aux=fac*etoti+auxf
 313 #ifdef DEBUG
 314             write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj,
 315      &        GridT0+k*deltaT,
 316      &        etoti+elowest_all,e_total(i,iprot),
 317      &        elowest_all,fac*etoti,auxf,aux
 318 #endif
 319             if (aux.le.150.0) then
 320               etoti=etoti-temper*eprim
 321               weight=dexp(-aux)
 322               sumW(k,iprot)=sumW(k,iprot)+weight
 323               sumE(k,iprot)=sumE(k,iprot)+etoti*weight
 324               sumEbis(k,iprot)=sumEbis(k,iprot)+ebis*weight
 325               sumEsq(k,iprot)=sumEsq(k,iprot)+etoti**2*weight
 326               sumEEE(k,iprot)=sumEEE(k,iprot)+etoti**3*weight
 327               if (ii.gt.0) then
 328               sumQ(k,iprot)=sumQ(k,iprot)+nu(1,ii,jj,iprot)*weight
 329               sumRMS(k,iprot)=sumRMS(k,iprot)
 330      &          +rmstb(i,iprot)*weight
 331 c     &       +nu(1,ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)*weight
 332               sumRGY(k,iprot)=sumRGY(k,iprot)
 333      &           +nu(1,ii+2,jj,iprot)*weight
 334               sumQsq(k,iprot)=sumQsq(k,iprot)
 335      &           +nu(1,ii,jj,iprot)**2*weight
 336               sumEQ(k,iprot)=sumEQ(k,iprot)
 337      &           +etoti*nu(1,ii,jj,iprot)*weight
 338               sumEEQ(k,iprot)=sumEEQ(k,iprot)+etoti**2*nu(1,ii,jj,iprot)
 339      &          *weight
 340               endif
 341               sumRMS(k,iprot)=sumRMS(k,iprot)+rmstb(i,iprot)*weight
 342               ib=1
 343               tt=dabs(1.0d0/(1.987D-3*betaT(1,1,iprot))-temper)
 344               do l=2,nbeta(1,iprot)
 345                 if(dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper).lt.tt)
 346      &          then
 347                   tt=dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper)
 348                   ib=l
 349                 endif
 350               enddo
 351 #ifdef DEBUG
 352               write (iout,*) "jj",jj," sigma2",sigma2(iprot),"Ptab",
 353      &           Ptab(jj,ib,iprot)," weight",weight
 354 #endif
 355               if (Ptab(jj,ib,iprot).gt.0.0d0) then
 356                 sumQ(k,iprot)=sumQ(k,iprot)+
 357      &           sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight
 358               else
 359                 sumQ(k,iprot)=sumQ(k,iprot)+1.0d2*weight
 360               endif
 361             endif
 362           enddo
 363         endif
 364       enddo
 365 #ifdef DEBUG
 366       write (iout,*) "calc_thermal"
 367       do k=0,nGridT
 368         write (iout,'(i5,8e15.5)') k,sumW(k,iprot),sumE(k,iprot),
 369      &   sumQ(k,iprot),sumEsq(k,iprot),sumRMS(k,iprot),sumQsq(k,iprot),
 370      &   sumEQ(k,iprot)
 371       enddo
 372 #endif
 373       return
 374       end
 375 c-------------------------------------------------------------------------------
 376       subroutine sum_thermal(iprot,elowest_all,entbase_all)
 377       implicit none
 378       include "DIMENSIONS"
 379       include "DIMENSIONS.ZSCOPT"
 380 #ifdef MPI
 381       include "mpif.h"
 382       integer IERROR,ErrCode
 383       include "COMMON.MPI"
 384       double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT),
 385      & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT),
 386      & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT),
 387      & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT)
 388       double precision rmstb_(maxstr_proc)
 389 #endif
 390       include "COMMON.THERMAL"
 391       include "COMMON.PROTNAME"
 392       include "COMMON.WEIGHTS"
 393       include "COMMON.WEIGHTDER"
 394       include "COMMON.CLASSES"
 395       include "COMMON.ENERGIES"
 396       include "COMMON.IOUNITS"
 397       include "COMMON.VMCPAR"
 398       include "COMMON.NAMES"
 399       include "COMMON.INTERACT"
 400       include "COMMON.TIME1"
 401       include "COMMON.CHAIN"
 402       include "COMMON.PROTFILES"
 403       include "COMMON.OPTIM"
 404 C Define local variables
 405       integer i,iprot
 406       integer ilen
 407       external ilen
 408       double precision elowest_all,entbase_all,RgasT,Temper
 409       character*64 nazwa
 410       character*4 liczba
 411       if (nparmset.eq.0) then
 412         nazwa=prefix(:ilen(prefix))
 413      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 414       else
 415         write (liczba,'(bz,i4.4)') myparm
 416         nazwa=prefix(:ilen(prefix))//'_par'//liczba
 417      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 418       endif
 419 #ifdef MPI
 420       call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1,
 421      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 422       call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1,
 423      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 424       call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1,
 425      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 426       call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1,
 427      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 428       call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1,
 429      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 430       call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1,
 431      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 432       call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1,
 433      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 434       call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1,
 435      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 436       call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1,
 437      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 438       call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1,
 439      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 440       call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1,
 441      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 442       do i=1,scount(me,iprot)
 443         rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
 444       enddo
 445 c      call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
 446 c     &   scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 447 c     &   scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 448 c     &   Comm1, IERROR)
 449       call MPI_Gatherv(rmstb_(1),
 450      &   scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 451      &   scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 452      &   Comm1, IERROR)
 453       if (me1.eq.master) then
 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 c          write (iout,*) "i",i," Temper",Temper," RgasT",RgasT,
 460 c     &     " entbase_all",entbase_all
 461           sumE(i,iprot)=sumE_t(i)/sumW_t(i)
 462           sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i)
 463           sumQ(i,iprot)=sumQ_t(i)/sumW_t(i)
 464           sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i)
 465           sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i)
 466           sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i)
 467           sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i)
 468      &     -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3
 469           sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2
 470           sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i)
 471           sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)-
 472      &      sumEsq(i,iprot)*sumQ(i,iprot)-
 473      &      2*sumEQ(i,iprot)*sumE(i,iprot)+
 474      &      2*sumQ(i,iprot)*sumE(i,iprot)**2
 475 #ifdef DEBUG
 476           write (iout,*) "Temper",temper," RgasT",RgasT
 477           write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
 478      &      " sumEbis",sumEbis(i,iprot)
 479 #endif
 480           sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 481           sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 482           sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all
 483           sumE(i,iprot)=sumE(i,iprot)+elowest_all
 484           write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 485      &      sumW(i,iprot),
 486      &      sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT
 487      &      +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 488      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 489      &      sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot),
 490      &      sumEEQ(i,iprot)
 491 #ifdef DEBUG
 492           write (iout,*) "zvtot",
 493      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
 494 #endif
 495           call flush(istat)
 496         enddo
 497         close (istat)
 498       endif
 499 #else
 500       open (istat,file=nazwa)
 501       write (istat,'(/a)') '# Thermal characteristics of folding'
 502       do i=0,NGridT
 503         Temper = GridT0+i*deltaT
 504         RgasT = Rgas*Temper
 505         sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot)
 506         sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot)
 507         sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot)
 508         sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot)
 509         sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot)
 510         sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot)
 511         sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2
 512         sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot)
 513         sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)-
 514      &    sumEsq(i,iprot)*sumQ(i,iprot)-
 515      &    2*sumEQ(i,iprot)*sumE(i,iprot)+
 516      &    2*sumQ(i,iprot)*sumE(i,iprot)**2
 517         sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 518         sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 519         sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all
 520         sumE(i,iprot)=sumE(i,iprot)+elowest_all
 521         write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 522      &   sumW(i,iprot),
 523      &   sumE(i,iprot),
 524      &   (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all,
 525      &   sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 526      &   sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 527      &   sumQsq(i,iprot),
 528      &   sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)
 529       enddo
 530       close (istat)
 531 #endif
 532       return
 533       end