subroutine thermal(iprot) c Calculate the average energy, average Q, specific heat, and Q dispersion c curves in temperature for protein IPROT. implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What include "COMMON.THERMAL" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.COMPAR" C Define local variables integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,l1,l2,ll1,ll2,ibatch integer ipass_conf,istart_conf,iend_conf logical lprn integer icant external icant double precision elowest_all,elowest_ent_all,entbase_all c write (iout,*) "Thermal" c call flush(iout) lprn=.true. elowest_all=1.0d20 elowest_ent_all=1.0d20 do ibatch=1,natlike(iprot)+2 do ib=1,nbeta(ibatch,iprot) c write (iout,*) "ibatch",ibatch," ib",ib call flush(iout) if (elowest(ib,ibatch,iprot).lt.elowest_all) & elowest_all=elowest(ib,ibatch,iprot) c write (iout,*) "ibatch",ibatch," ib",ib," after if" call flush(iout) if (elowest_ent(ib,ibatch,iprot).lt.elowest_ent_all) then elowest_ent_all = elowest_ent(ib,ibatch,iprot) endif enddo enddo c write (iout,*) "THERMAL: elowest_all",elowest_all call flush(iout) c c Zero out tables. c do i=0,nGridT sumW(i,iprot)=0.0d0 sumE(i,iprot)=0.0d0 sumEbis(i,iprot)=0.0d0 sumEsq(i,iprot)=0.0d0 sumEEE(i,iprot)=0.0d0 sumQ(i,iprot)=0.0d0 sumRMS(i,iprot)=0.0d0 sumRGY(i,iprot)=0.0d0 sumQsq(i,iprot)=0.0d0 sumEQ(i,iprot)=0.0d0 sumEEQ(i,iprot)=0.0d0 enddo c c Calculate the contributions to thermodynamic quantities from all processors. c #ifdef DEBUG write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot) #endif IF (NCHUNK_CONF(IPROT).EQ.1) THEN #ifdef MPI do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=indstart(me1,iprot),indend(me1,iprot) ii=ii+1 i2ii(i,iprot)=ii enddo #else do i=1,ntot_work(iprot) i2ii(i,iprot)=i endif #endif #ifdef DEBUG do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif call calc_thermal(iprot,elowest_all,*10) ELSE open (ientin,file=benefiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) ipass_conf=0 #ifdef MPI do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot)) #else do istart_conf=1,ntot_work(iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot)) #endif c c Read the chunk of energies and derivatives off a DA scratchfile. c ipass_conf=ipass_conf+1 do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=istart_conf,iend_conf ii=ii+1 i2ii(i,iprot)=ii enddo #ifdef DEBUG write (iout,*) "ipass_conf",ipass_conf, & " istart_conf",istart_conf," iend_conf",iend_conf do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif call daread_ene(iprot,istart_conf,iend_conf) call calc_thermal(iprot,elowest_all,*10) enddo close(ientin) ENDIF 10 continue c c Put tohether all contributions. c call sum_thermal(iprot,elowest_all,entbase_all) return end c------------------------------------------------------------------------------- subroutine calc_thermal(iprot,elowest_all,*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What include "COMMON.CONTROL" include "COMMON.THERMAL" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.CLASSES" include "COMMON.COMPAR" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.FFIELD" include "COMMON.PROTFILES" include "COMMON.SBRIDGE" include "COMMON.ALLPROT" C Define local variables integer i,ii,jj,j,k,kk,l,m,iprot,ibatch,istart_conf,iend_conf,ib double precision aux,auxf,fac,weight,etoti,eprim,ebis,denom, & tt,temper,quot,quotl,quotl1,kfacl,logfac,eplus,eminus,tanhT double precision elowest_all,entbase_all double precision ft(5),ftprim(5),ftbis(5),escal1(max_ene), & escal1prim(max_ene),escal1bis(max_ene) double precision T0 /300.0d0/, kfac /2.4d0/ ii=natlike(iprot)-2 c write (iout,*) "ii",ii c write (iout,*) "Thermal: the i2ii array" c do i=1,ntot_work(iprot) c write (iout,*) i,i2ii(i) c enddo c write (iout,*) "The list array" c do i=1,ntot_work(iprot) c write (iout,*) i,list_conf(i) c enddo c write (iout,*) "iprot",iprot," natlike",natlike(iprot) do i=1,ntot_work(iprot) jj = i2ii(i,iprot) if (jj.gt.0) then c etoti=e_total(i,iprot)-elowest_all auxf=entfac(i,iprot) c write (iout,*) "etoti",etoti+elowest_all," auxf",auxf," nu", c & nu(ii,jj,iprot) c write (iout,'(2i5,20f8.2)') i,jj,(enetb(jj,kk,iprot), c & kk=1,n_ene) c if (ii.gt.0) rmstb(i,iprot)=nu(1,ii,jj,iprot) c write (iout,*) "i",i," jj",jj," nu",(nu(k,jj,iprot), c & k=1,natlike(iprot)+2) c write (iout,*) "i",i," rmstb",rmstb(i,iprot) do k=0,nGridT temper=GridT0+k*deltaT fac=1.0d0/(1.987D-3*temper) quot=temper/T0 if (rescale_mode.eq.0) then do l=1,5 fT(l)=1.0d0 fTprim(l)=0.0d0 fTbis(l)=0.0d0 enddo else if (rescale_mode.eq.1) then quotl=1.0d0 kfacl=1.0d0 do l=1,5 quotl1=quotl quotl=quotl*quot kfacl=kfacl*kfac denom=kfacl-1.0d0+quotl fT(l)=kfacl/denom ftprim(l)=-l*ft(l)*quotl1/(T0*denom) ftbis(l)=l*kfacl*quotl1* & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3) enddo else if (rescale_mode.eq.2) then quotl=1.0d0 do l=1,5 quotl1=quotl quotl=quotl*quot eplus=dexp(quotl) eminus=dexp(-quotl) logfac=1.0d0/dlog(eplus+eminus) tanhT=(eplus-eminus)/(eplus+eminus) fT(l)=1.12692801104297249644d0*logfac ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- & 2*l*quotl1/T0*logfac* & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) & +ftprim(l)*tanhT) enddo else write (iout,*) "Unknown RESCALE_MODE",rescale_mode call flush(iout) return1 endif do l=1,n_ene escal1(l)=1.0d0 escal1prim(l)=0.0d0 escal1bis(l)=0.0d0 enddo if (shield_mode.gt.0) then escal1(1)=ft(1) escal1(2)=ft(1) escal1(16)=ft(1) escal1prim(1)=ftprim(1) escal1prim(2)=ftprim(1) escal1prim(16)=ftprim(1) escal1bis(1)=ftbis(1) escal1bis(2)=ftbis(1) escal1bis(16)=ftbis(1) endif escal1(3)=ft(1) escal1(4)=ft(3) escal1(5)=ft(4) escal1(6)=ft(5) escal1(7)=ft(2) escal1(8)=ft(2) escal1(9)=ft(3) escal1(10)=ft(5) escal1(13)=ft(1) escal1(14)=ft(2) escal1(19)=ft(1) escal1prim(3)=ftprim(1) escal1prim(4)=ftprim(3) escal1prim(5)=ftprim(4) escal1prim(6)=ftprim(5) escal1prim(7)=ftprim(2) escal1prim(8)=ftprim(2) escal1prim(9)=ftprim(3) escal1prim(10)=ftprim(5) escal1prim(13)=ftprim(1) escal1prim(14)=ftprim(2) escal1prim(19)=ftprim(1) escal1bis(3)=ftbis(1) escal1bis(4)=ftbis(3) escal1bis(5)=ftbis(4) escal1bis(6)=ftbis(5) escal1bis(7)=ftbis(2) escal1bis(8)=ftbis(2) escal1bis(9)=ftbis(3) escal1bis(10)=ftbis(5) escal1bis(13)=ftbis(1) escal1bis(14)=ftbis(2) escal1bis(19)=ftbis(1) etoti=0.0d0 eprim=0.0d0 ebis=0.0d0 do l=1,n_ene etoti=etoti+ww(l)*escal1(l)*enetb(jj,l,iprot) eprim=eprim+ww(l)*escal1prim(l)*enetb(jj,l,iprot) ebis=ebis+ww(l)*escal1bis(l)*enetb(jj,l,iprot) enddo #ifdef DEBUG write (iout,'(4i5,3f10.5)') i,ii,jj,iprot,GridT0+k*deltaT, & etoti,nu(ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1) #endif etoti=etoti-elowest_all aux=fac*etoti+auxf #ifdef DEBUG write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj, & GridT0+k*deltaT, & etoti+elowest_all,e_total(i,iprot), & elowest_all,fac*etoti,auxf,aux #endif if (aux.le.150.0) then etoti=etoti-temper*eprim weight=dexp(-aux) sumW(k,iprot)=sumW(k,iprot)+weight sumE(k,iprot)=sumE(k,iprot)+etoti*weight sumEbis(k,iprot)=sumEbis(k,iprot)+ebis*weight sumEsq(k,iprot)=sumEsq(k,iprot)+etoti**2*weight sumEEE(k,iprot)=sumEEE(k,iprot)+etoti**3*weight if (ii.gt.0) then sumQ(k,iprot)=sumQ(k,iprot)+nu(1,ii,jj,iprot)*weight sumRMS(k,iprot)=sumRMS(k,iprot) & +rmstb(i,iprot)*weight c & +nu(1,ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)*weight sumRGY(k,iprot)=sumRGY(k,iprot) & +nu(1,ii+2,jj,iprot)*weight sumQsq(k,iprot)=sumQsq(k,iprot) & +nu(1,ii,jj,iprot)**2*weight sumEQ(k,iprot)=sumEQ(k,iprot) & +etoti*nu(1,ii,jj,iprot)*weight sumEEQ(k,iprot)=sumEEQ(k,iprot)+etoti**2*nu(1,ii,jj,iprot) & *weight endif sumRMS(k,iprot)=sumRMS(k,iprot)+rmstb(i,iprot)*weight ib=1 tt=dabs(1.0d0/(1.987D-3*betaT(1,1,iprot))-temper) do l=2,nbeta(1,iprot) if(dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper).lt.tt) & then tt=dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper) ib=l endif enddo sumQ(k,iprot)=sumQ(k,iprot)+ & sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight endif enddo endif enddo #ifdef DEBUG do k=0,nGridT write (iout,'(i5,7e15.5)') k,sumW(k,iprot),sumE(k,iprot), & sumEsq(k,iprot),sumRMS(k,iprot),sumQsq(k,iprot),sumEQ(k,iprot) enddo #endif return end c------------------------------------------------------------------------------- subroutine sum_thermal(iprot,elowest_all,entbase_all) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT), & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT), & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT), & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT) double precision rmstb_(maxstr_proc) #endif include "COMMON.THERMAL" include "COMMON.PROTNAME" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.OPTIM" C Define local variables integer i,iprot integer ilen external ilen double precision elowest_all,entbase_all,RgasT,Temper character*64 nazwa character*4 liczba if (nparmset.eq.0) then nazwa=prefix(:ilen(prefix)) & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal' else write (liczba,'(bz,i4.4)') myparm nazwa=prefix(:ilen(prefix))//'_par'//liczba & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal' endif #ifdef MPI call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR) do i=1,scount(me,iprot) rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot) enddo c call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, c & Comm1, IERROR) call MPI_Gatherv(rmstb_(1), & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) if (me1.eq.master) then open (istat,file=nazwa) write (istat,'(a)') '# Thermal characteristics of folding' do i=0,NGridT Temper = GridT0+i*deltaT RgasT = Rgas*Temper c write (iout,*) "i",i," Temper",Temper," RgasT",RgasT, c & " entbase_all",entbase_all sumE(i,iprot)=sumE_t(i)/sumW_t(i) sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i) sumQ(i,iprot)=sumQ_t(i)/sumW_t(i) sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i) sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i) sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i) sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i) & -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3 sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2 sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i) sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)- & sumEsq(i,iprot)*sumQ(i,iprot)- & 2*sumEQ(i,iprot)*sumE(i,iprot)+ & 2*sumQ(i,iprot)*sumE(i,iprot)**2 #ifdef DEBUG write (iout,*) "Temper",temper," RgasT",RgasT write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot), & " sumEbis",sumEbis(i,iprot) #endif sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot) sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2 sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all sumE(i,iprot)=sumE(i,iprot)+elowest_all write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper, & sumW(i,iprot), & sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT & +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot), & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot), & sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot), & sumEEQ(i,iprot) #ifdef DEBUG write (iout,*) "zvtot", & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot) #endif call flush(istat) enddo close (istat) endif #else open (istat,file=nazwa) write (istat,'(/a)') '# Thermal characteristics of folding' do i=0,NGridT Temper = GridT0+i*deltaT RgasT = Rgas*Temper sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot) sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot) sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot) sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot) sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot) sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot) sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2 sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot) sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)- & sumEsq(i,iprot)*sumQ(i,iprot)- & 2*sumEQ(i,iprot)*sumE(i,iprot)+ & 2*sumQ(i,iprot)*sumE(i,iprot)**2 sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot) sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2 sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all sumE(i,iprot)=sumE(i,iprot)+elowest_all write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper, & sumW(i,iprot), & sumE(i,iprot), & (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all, & sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot), & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot), & sumQsq(i,iprot), & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot) enddo close (istat) #endif return end