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.
6 include "DIMENSIONS.ZSCOPT"
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
32 double precision elowest_all,elowest_ent_all,entbase_all
34 c write (iout,*) "Thermal"
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
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"
47 if (elowest_ent(ib,ibatch,iprot).lt.elowest_ent_all) then
48 elowest_ent_all = elowest_ent(ib,ibatch,iprot)
52 c write (iout,*) "THERMAL: elowest_all",elowest_all
60 sumEbis(i,iprot)=0.0d0
71 c Calculate the contributions to thermodynamic quantities from all processors.
74 write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot)
76 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
79 do i=1,ntot_work(iprot)
83 do i=indstart(me1,iprot),indend(me1,iprot)
88 do i=1,ntot_work(iprot)
93 do i=1,ntot_work(iprot)
94 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
98 call calc_thermal(iprot,elowest_all,*10)
102 open (ientin,file=benefiles(iprot),status="old",
103 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
106 do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
107 iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
109 do istart_conf=1,ntot_work(iprot),maxstr_proc
110 iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
113 c Read the chunk of energies and derivatives off a DA scratchfile.
115 ipass_conf=ipass_conf+1
116 do i=1,ntot_work(iprot)
120 do i=istart_conf,iend_conf
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)
132 call daread_ene(iprot,istart_conf,iend_conf)
133 call calc_thermal(iprot,elowest_all,*10)
140 c Put tohether all contributions.
142 call sum_thermal(iprot,elowest_all,entbase_all)
145 c-------------------------------------------------------------------------------
146 subroutine calc_thermal(iprot,elowest_all,*)
149 include "DIMENSIONS.ZSCOPT"
152 integer IERROR,ErrCode
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/
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)
187 c write (iout,*) "The list array"
188 c do i=1,ntot_work(iprot)
189 c write (iout,*) i,list_conf(i)
191 c write (iout,*) "iprot",iprot," natlike",natlike(iprot)
192 do i=1,ntot_work(iprot)
195 c etoti=e_total(i,iprot)-elowest_all
197 c write (iout,*) "etoti",etoti+elowest_all," auxf",auxf," nu",
199 c write (iout,'(2i5,20f8.2)') i,jj,(enetb(jj,kk,iprot),
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)
206 temper=GridT0+k*deltaT
207 fac=1.0d0/(1.987D-3*temper)
209 if (rescale_mode.eq.0) then
215 else if (rescale_mode.eq.1) then
222 denom=kfacl-1.0d0+quotl
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)
228 else if (rescale_mode.eq.2) then
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)
245 write (iout,*) "Unknown RESCALE_MODE",rescale_mode
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)
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)
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)
299 etoti=etoti-elowest_all
302 write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj,
304 & etoti+elowest_all,e_total(i,iprot),
305 & elowest_all,fac*etoti,auxf,aux
307 if (aux.le.150.0) then
308 etoti=etoti-temper*eprim
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
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)
329 sumRMS(k,iprot)=sumRMS(k,iprot)+rmstb(i,iprot)*weight
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)
335 tt=dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper)
339 sumQ(k,iprot)=sumQ(k,iprot)+
340 & sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight
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)
353 c-------------------------------------------------------------------------------
354 subroutine sum_thermal(iprot,elowest_all,entbase_all)
357 include "DIMENSIONS.ZSCOPT"
360 integer IERROR,ErrCode
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)
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
386 double precision elowest_all,entbase_all,RgasT,Temper
389 if (nparmset.eq.0) then
390 nazwa=prefix(:ilen(prefix))
391 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
393 write (liczba,'(bz,i3.3)') myparm
394 nazwa=prefix(:ilen(prefix))//'_par'//liczba
395 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
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)
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,
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,
431 if (me1.eq.master) then
432 open (istat,file=nazwa)
433 write (istat,'(a)') '# Thermal characteristics of folding'
435 Temper = GridT0+i*deltaT
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
454 write (iout,*) "Temper",temper," RgasT",RgasT
455 write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
456 & " sumEbis",sumEbis(i,iprot)
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,
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),
470 write (iout,*) "zvtot",
471 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
478 open (istat,file=nazwa)
479 write (istat,'(/a)') '# Thermal characteristics of folding'
481 Temper = GridT0+i*deltaT
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,
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),
506 & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)