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.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/
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)
188 c write (iout,*) "The list array"
189 c do i=1,ntot_work(iprot)
190 c write (iout,*) i,list_conf(i)
192 c write (iout,*) "iprot",iprot," natlike",natlike(iprot)
193 do i=1,ntot_work(iprot)
196 c etoti=e_total(i,iprot)-elowest_all
198 c write (iout,*) "etoti",etoti+elowest_all," auxf",auxf," nu",
200 c write (iout,'(2i5,20f8.2)') i,jj,(enetb(jj,kk,iprot),
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)
207 temper=GridT0+k*deltaT
208 fac=1.0d0/(1.987D-3*temper)
210 if (rescale_mode.eq.0) then
216 else if (rescale_mode.eq.1) then
223 denom=kfacl-1.0d0+quotl
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)
229 else if (rescale_mode.eq.2) then
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)
246 write (iout,*) "Unknown RESCALE_MODE",rescale_mode
255 if (shield_mode.gt.0) then
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)
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)
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)
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)
311 etoti=etoti-elowest_all
314 write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj,
316 & etoti+elowest_all,e_total(i,iprot),
317 & elowest_all,fac*etoti,auxf,aux
319 if (aux.le.150.0) then
320 etoti=etoti-temper*eprim
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
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)
341 sumRMS(k,iprot)=sumRMS(k,iprot)+rmstb(i,iprot)*weight
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)
347 tt=dabs(1.0d0/(1.987D-3*betaT(l,1,iprot))-temper)
352 write (iout,*) "jj",jj," sigma2",sigma2(iprot),"Ptab",
353 & Ptab(jj,ib,iprot)," weight",weight
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
359 sumQ(k,iprot)=sumQ(k,iprot)+1.0d2*weight
366 write (iout,*) "calc_thermal"
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),
375 c-------------------------------------------------------------------------------
376 subroutine sum_thermal(iprot,elowest_all,entbase_all)
379 include "DIMENSIONS.ZSCOPT"
382 integer IERROR,ErrCode
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)
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
408 double precision elowest_all,entbase_all,RgasT,Temper
411 if (nparmset.eq.0) then
412 nazwa=prefix(:ilen(prefix))
413 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
415 write (liczba,'(bz,i4.4)') myparm
416 nazwa=prefix(:ilen(prefix))//'_par'//liczba
417 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
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)
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,
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,
453 if (me1.eq.master) then
454 open (istat,file=nazwa)
455 write (istat,'(a)') '# Thermal characteristics of folding'
457 Temper = GridT0+i*deltaT
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
476 write (iout,*) "Temper",temper," RgasT",RgasT
477 write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
478 & " sumEbis",sumEbis(i,iprot)
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,
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),
492 write (iout,*) "zvtot",
493 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
500 open (istat,file=nazwa)
501 write (istat,'(/a)') '# Thermal characteristics of folding'
503 Temper = GridT0+i*deltaT
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,
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),
528 & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)