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