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)
351 sumQ(k,iprot)=sumQ(k,iprot)+
352 & sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight
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)
365 c-------------------------------------------------------------------------------
366 subroutine sum_thermal(iprot,elowest_all,entbase_all)
369 include "DIMENSIONS.ZSCOPT"
372 integer IERROR,ErrCode
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)
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
398 double precision elowest_all,entbase_all,RgasT,Temper
401 if (nparmset.eq.0) then
402 nazwa=prefix(:ilen(prefix))
403 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
405 write (liczba,'(bz,i4.4)') myparm
406 nazwa=prefix(:ilen(prefix))//'_par'//liczba
407 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
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)
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,
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,
443 if (me1.eq.master) then
444 open (istat,file=nazwa)
445 write (istat,'(a)') '# Thermal characteristics of folding'
447 Temper = GridT0+i*deltaT
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
466 write (iout,*) "Temper",temper," RgasT",RgasT
467 write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
468 & " sumEbis",sumEbis(i,iprot)
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,
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),
482 write (iout,*) "zvtot",
483 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
490 open (istat,file=nazwa)
491 write (istat,'(/a)') '# Thermal characteristics of folding'
493 Temper = GridT0+i*deltaT
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,
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),
518 & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)