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 if (Ptab(jj,ib,iprot).gt.0.0d0) then
352 sumQ(k,iprot)=sumQ(k,iprot)+
353 & sqrt(-2*sigma2(iprot)*dlog(Ptab(jj,ib,iprot)))*weight
355 sumQ(k,iprot)=sumQ(k,iprot)+1.0d2*weight
363 write (iout,'(i5,7e15.5)') k,sumW(k,iprot),sumE(k,iprot),
364 & sumEsq(k,iprot),sumRMS(k,iprot),sumQsq(k,iprot),sumEQ(k,iprot)
369 c-------------------------------------------------------------------------------
370 subroutine sum_thermal(iprot,elowest_all,entbase_all)
373 include "DIMENSIONS.ZSCOPT"
376 integer IERROR,ErrCode
378 double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT),
379 & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT),
380 & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT),
381 & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT)
382 double precision rmstb_(maxstr_proc)
384 include "COMMON.THERMAL"
385 include "COMMON.PROTNAME"
386 include "COMMON.WEIGHTS"
387 include "COMMON.WEIGHTDER"
388 include "COMMON.CLASSES"
389 include "COMMON.ENERGIES"
390 include "COMMON.IOUNITS"
391 include "COMMON.VMCPAR"
392 include "COMMON.NAMES"
393 include "COMMON.INTERACT"
394 include "COMMON.TIME1"
395 include "COMMON.CHAIN"
396 include "COMMON.PROTFILES"
397 include "COMMON.OPTIM"
398 C Define local variables
402 double precision elowest_all,entbase_all,RgasT,Temper
405 if (nparmset.eq.0) then
406 nazwa=prefix(:ilen(prefix))
407 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
409 write (liczba,'(bz,i4.4)') myparm
410 nazwa=prefix(:ilen(prefix))//'_par'//liczba
411 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
414 call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1,
415 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
416 call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1,
417 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
418 call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1,
419 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
420 call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1,
421 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
422 call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1,
423 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
424 call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1,
425 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
426 call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1,
427 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
428 call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1,
429 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
430 call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1,
431 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
432 call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1,
433 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
434 call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1,
435 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
436 do i=1,scount(me,iprot)
437 rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
439 c call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
440 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
441 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
443 call MPI_Gatherv(rmstb_(1),
444 & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
445 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
447 if (me1.eq.master) then
448 open (istat,file=nazwa)
449 write (istat,'(a)') '# Thermal characteristics of folding'
451 Temper = GridT0+i*deltaT
453 c write (iout,*) "i",i," Temper",Temper," RgasT",RgasT,
454 c & " entbase_all",entbase_all
455 sumE(i,iprot)=sumE_t(i)/sumW_t(i)
456 sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i)
457 sumQ(i,iprot)=sumQ_t(i)/sumW_t(i)
458 sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i)
459 sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i)
460 sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i)
461 sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i)
462 & -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3
463 sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2
464 sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i)
465 sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)-
466 & sumEsq(i,iprot)*sumQ(i,iprot)-
467 & 2*sumEQ(i,iprot)*sumE(i,iprot)+
468 & 2*sumQ(i,iprot)*sumE(i,iprot)**2
470 write (iout,*) "Temper",temper," RgasT",RgasT
471 write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
472 & " sumEbis",sumEbis(i,iprot)
474 sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
475 sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
476 sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all
477 sumE(i,iprot)=sumE(i,iprot)+elowest_all
478 write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
480 & sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT
481 & +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
482 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
483 & sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot),
486 write (iout,*) "zvtot",
487 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
494 open (istat,file=nazwa)
495 write (istat,'(/a)') '# Thermal characteristics of folding'
497 Temper = GridT0+i*deltaT
499 sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot)
500 sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot)
501 sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot)
502 sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot)
503 sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot)
504 sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot)
505 sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2
506 sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot)
507 sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)-
508 & sumEsq(i,iprot)*sumQ(i,iprot)-
509 & 2*sumEQ(i,iprot)*sumE(i,iprot)+
510 & 2*sumQ(i,iprot)*sumE(i,iprot)**2
511 sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
512 sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
513 sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all
514 sumE(i,iprot)=sumE(i,iprot)+elowest_all
515 write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
518 & (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all,
519 & sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
520 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
522 & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)