source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/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               sumQ(k,iprot)=sumQ(k,iprot)+
 352      &           sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight
 353             endif
 354           enddo
 355         endif
 356       enddo
 357 #ifdef DEBUG
 358       do k=0,nGridT
 359         write (iout,'(i5,7e15.5)') k,sumW(k,iprot),sumE(k,iprot),
 360      &   sumEsq(k,iprot),sumRMS(k,iprot),sumQsq(k,iprot),sumEQ(k,iprot)
 361       enddo
 362 #endif
 363       return
 364       end
 365 c-------------------------------------------------------------------------------
 366       subroutine sum_thermal(iprot,elowest_all,entbase_all)
 367       implicit none
 368       include "DIMENSIONS"
 369       include "DIMENSIONS.ZSCOPT"
 370 #ifdef MPI
 371       include "mpif.h"
 372       integer IERROR,ErrCode
 373       include "COMMON.MPI"
 374       double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT),
 375      & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT),
 376      & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT),
 377      & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT)
 378       double precision rmstb_(maxstr_proc)
 379 #endif
 380       include "COMMON.THERMAL"
 381       include "COMMON.PROTNAME"
 382       include "COMMON.WEIGHTS"
 383       include "COMMON.WEIGHTDER"
 384       include "COMMON.CLASSES"
 385       include "COMMON.ENERGIES"
 386       include "COMMON.IOUNITS"
 387       include "COMMON.VMCPAR"
 388       include "COMMON.NAMES"
 389       include "COMMON.INTERACT"
 390       include "COMMON.TIME1"
 391       include "COMMON.CHAIN"
 392       include "COMMON.PROTFILES"
 393       include "COMMON.OPTIM"
 394 C Define local variables
 395       integer i,iprot
 396       integer ilen
 397       external ilen
 398       double precision elowest_all,entbase_all,RgasT,Temper
 399       character*64 nazwa
 400       character*4 liczba
 401       if (nparmset.eq.0) then
 402         nazwa=prefix(:ilen(prefix))
 403      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 404       else
 405         write (liczba,'(bz,i4.4)') myparm
 406         nazwa=prefix(:ilen(prefix))//'_par'//liczba
 407      &   //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
 408       endif
 409 #ifdef MPI
 410       call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1,
 411      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 412       call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1,
 413      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 414       call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1,
 415      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 416       call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1,
 417      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 418       call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1,
 419      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 420       call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1,
 421      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 422       call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1,
 423      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 424       call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1,
 425      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 426       call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1,
 427      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 428       call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1,
 429      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 430       call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1,
 431      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
 432       do i=1,scount(me,iprot)
 433         rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
 434       enddo
 435 c      call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
 436 c     &   scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 437 c     &   scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 438 c     &   Comm1, IERROR)
 439       call MPI_Gatherv(rmstb_(1),
 440      &   scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
 441      &   scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
 442      &   Comm1, IERROR)
 443       if (me1.eq.master) then
 444         open (istat,file=nazwa)
 445         write (istat,'(a)') '# Thermal characteristics of folding'
 446         do i=0,NGridT
 447           Temper = GridT0+i*deltaT
 448           RgasT = Rgas*Temper
 449 c          write (iout,*) "i",i," Temper",Temper," RgasT",RgasT,
 450 c     &     " entbase_all",entbase_all
 451           sumE(i,iprot)=sumE_t(i)/sumW_t(i)
 452           sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i)
 453           sumQ(i,iprot)=sumQ_t(i)/sumW_t(i)
 454           sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i)
 455           sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i)
 456           sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i)
 457           sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i)
 458      &     -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3
 459           sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2
 460           sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i)
 461           sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)-
 462      &      sumEsq(i,iprot)*sumQ(i,iprot)-
 463      &      2*sumEQ(i,iprot)*sumE(i,iprot)+
 464      &      2*sumQ(i,iprot)*sumE(i,iprot)**2
 465 #ifdef DEBUG
 466           write (iout,*) "Temper",temper," RgasT",RgasT
 467           write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
 468      &      " sumEbis",sumEbis(i,iprot)
 469 #endif
 470           sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 471           sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 472           sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all
 473           sumE(i,iprot)=sumE(i,iprot)+elowest_all
 474           write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 475      &      sumW(i,iprot),
 476      &      sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT
 477      &      +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 478      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 479      &      sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot),
 480      &      sumEEQ(i,iprot)
 481 #ifdef DEBUG
 482           write (iout,*) "zvtot",
 483      &      sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
 484 #endif
 485           call flush(istat)
 486         enddo
 487         close (istat)
 488       endif
 489 #else
 490       open (istat,file=nazwa)
 491       write (istat,'(/a)') '# Thermal characteristics of folding'
 492       do i=0,NGridT
 493         Temper = GridT0+i*deltaT
 494         RgasT = Rgas*Temper
 495         sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot)
 496         sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot)
 497         sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot)
 498         sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot)
 499         sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot)
 500         sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot)
 501         sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2
 502         sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot)
 503         sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)-
 504      &    sumEsq(i,iprot)*sumQ(i,iprot)-
 505      &    2*sumEQ(i,iprot)*sumE(i,iprot)+
 506      &    2*sumQ(i,iprot)*sumE(i,iprot)**2
 507         sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
 508         sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
 509         sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all
 510         sumE(i,iprot)=sumE(i,iprot)+elowest_all
 511         write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
 512      &   sumW(i,iprot),
 513      &   sumE(i,iprot),
 514      &   (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all,
 515      &   sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
 516      &   sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
 517      &   sumQsq(i,iprot),
 518      &   sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)
 519       enddo
 520       close (istat)
 521 #endif
 522       return
 523       end