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
174 double precision aux,auxf,fac,weight,etoti,eprim,ebis,denom,
175 & 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 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)
205 temper=GridT0+k*deltaT
206 fac=1.0d0/(1.987D-3*temper)
208 if (rescale_mode.eq.0) then
214 else if (rescale_mode.eq.1) then
221 denom=kfacl-1.0d0+quotl
223 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
224 ftbis(l)=l*kfacl*quotl1*
225 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
227 else if (rescale_mode.eq.2) then
234 logfac=1.0d0/dlog(eplus+eminus)
235 tanhT=(eplus-eminus)/(eplus+eminus)
236 fT(l)=1.12692801104297249644d0*logfac
237 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
238 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
239 & 2*l*quotl1/T0*logfac*
240 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
244 write (iout,*) "Unknown RESCALE_MODE",rescale_mode
264 escal1prim(3)=ftprim(1)
265 escal1prim(4)=ftprim(3)
266 escal1prim(5)=ftprim(4)
267 escal1prim(6)=ftprim(5)
268 escal1prim(7)=ftprim(2)
269 escal1prim(8)=ftprim(2)
270 escal1prim(9)=ftprim(3)
271 escal1prim(10)=ftprim(5)
272 escal1prim(13)=ftprim(1)
273 escal1prim(14)=ftprim(2)
274 escal1prim(19)=ftprim(1)
275 escal1bis(3)=ftbis(1)
276 escal1bis(4)=ftbis(3)
277 escal1bis(5)=ftbis(4)
278 escal1bis(6)=ftbis(5)
279 escal1bis(7)=ftbis(2)
280 escal1bis(8)=ftbis(2)
281 escal1bis(9)=ftbis(3)
282 escal1bis(10)=ftbis(5)
283 escal1bis(13)=ftbis(1)
284 escal1bis(14)=ftbis(2)
285 escal1bis(19)=ftbis(1)
290 etoti=etoti+ww(l)*escal1(l)*enetb(jj,l,iprot)
291 eprim=eprim+ww(l)*escal1prim(l)*enetb(jj,l,iprot)
292 ebis=ebis+ww(l)*escal1bis(l)*enetb(jj,l,iprot)
295 write (iout,'(4i5,3f10.5)') i,ii,jj,iprot,GridT0+k*deltaT,
296 & etoti,nu(ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)
298 etoti=etoti-elowest_all
301 write (iout,'(a,2i5,f7.1,6f10.3)') "etoti",i,jj,
303 & etoti+elowest_all,e_total(i,iprot),
304 & elowest_all,fac*etoti,auxf,aux
306 if (aux.le.150.0) then
307 etoti=etoti-temper*eprim
309 sumW(k,iprot)=sumW(k,iprot)+weight
310 sumE(k,iprot)=sumE(k,iprot)+etoti*weight
311 sumEbis(k,iprot)=sumEbis(k,iprot)+ebis*weight
312 sumEsq(k,iprot)=sumEsq(k,iprot)+etoti**2*weight
313 sumEEE(k,iprot)=sumEEE(k,iprot)+etoti**3*weight
315 sumQ(k,iprot)=sumQ(k,iprot)+nu(1,ii,jj,iprot)*weight
316 sumRMS(k,iprot)=sumRMS(k,iprot)
317 & +nu(1,ii+1,jj,iprot)*(nct_zs(iprot)-nnt_zs(iprot)+1)*weight
318 sumRGY(k,iprot)=sumRGY(k,iprot)
319 & +nu(1,ii+2,jj,iprot)*weight
320 sumQsq(k,iprot)=sumQsq(k,iprot)
321 & +nu(1,ii,jj,iprot)**2*weight
322 sumEQ(k,iprot)=sumEQ(k,iprot)
323 & +etoti*nu(1,ii,jj,iprot)*weight
324 sumEEQ(k,iprot)=sumEEQ(k,iprot)+etoti**2*nu(1,ii,jj,iprot)
332 c write (iout,'(i5,7e15.5)') k,sumW(k,iprot),sumE(k,iprot),
333 c & sumEsq(k,iprot),sumQ(k,iprot),sumQsq(k,iprot),sumEQ(k,iprot)
337 c-------------------------------------------------------------------------------
338 subroutine sum_thermal(iprot,elowest_all,entbase_all)
341 include "DIMENSIONS.ZSCOPT"
344 integer IERROR,ErrCode
346 double precision sumW_t(0:MaxGridT),sumE_t(0:MaxGridT),
347 & sumEbis_t(0:MaxGridT),sumRMS_t(0:MaxGridT),sumRGY_t(0:MaxGridT),
348 & sumEsq_t(0:MaxGridT),sumQ_t(0:MaxGridT),sumQsq_t(0:MaxGridT),
349 & sumEQ_t(0:MaxGridT),sumEEQ_t(0:MaxGridT),sumEEE_t(0:MaxGridT)
351 include "COMMON.THERMAL"
352 include "COMMON.PROTNAME"
353 include "COMMON.WEIGHTS"
354 include "COMMON.WEIGHTDER"
355 include "COMMON.CLASSES"
356 include "COMMON.ENERGIES"
357 include "COMMON.IOUNITS"
358 include "COMMON.VMCPAR"
359 include "COMMON.NAMES"
360 include "COMMON.INTERACT"
361 include "COMMON.TIME1"
362 include "COMMON.CHAIN"
363 include "COMMON.PROTFILES"
364 include "COMMON.OPTIM"
365 C Define local variables
369 double precision elowest_all,entbase_all,RgasT,Temper
372 if (nparmset.eq.0) then
373 nazwa=prefix(:ilen(prefix))
374 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
376 write (liczba,'(bz,i3.3)') myparm
377 nazwa=prefix(:ilen(prefix))//'_par'//liczba
378 & //'.'//protname(iprot)(:ilen(protname(iprot)))//'.'//'thermal'
381 call MPI_Reduce(sumW(0,iprot),sumW_t(0),nGridT+1,
382 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
383 call MPI_Reduce(sumE(0,iprot),sumE_t(0),nGridT+1,
384 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
385 call MPI_Reduce(sumEbis(0,iprot),sumEbis_t(0),nGridT+1,
386 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
387 call MPI_Reduce(sumEsq(0,iprot),sumEsq_t(0),nGridT+1,
388 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
389 call MPI_Reduce(sumEEE(0,iprot),sumEEE_t(0),nGridT+1,
390 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
391 call MPI_Reduce(sumQ(0,iprot),sumQ_t(0),nGridT+1,
392 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
393 call MPI_Reduce(sumRMS(0,iprot),sumRMS_t(0),nGridT+1,
394 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
395 call MPI_Reduce(sumRGY(0,iprot),sumRGY_t(0),nGridT+1,
396 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
397 call MPI_Reduce(sumQsq(0,iprot),sumQsq_t(0),nGridT+1,
398 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
399 call MPI_Reduce(sumEQ(0,iprot),sumEQ_t(0),nGridT+1,
400 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
401 call MPI_Reduce(sumEEQ(0,iprot),sumEEQ_t(0),nGridT+1,
402 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,Comm1,IERROR)
403 call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
404 & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
405 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
407 if (me1.eq.master) then
408 open (istat,file=nazwa)
409 write (istat,'(a)') '# Thermal characteristics of folding'
411 Temper = GridT0+i*deltaT
413 c write (iout,*) "i",i," Temper",Temper," RgasT",RgasT,
414 c & " entbase_all",entbase_all
415 sumE(i,iprot)=sumE_t(i)/sumW_t(i)
416 sumEbis(i,iprot)=Temper*sumEbis_t(i)/sumW_t(i)
417 sumQ(i,iprot)=sumQ_t(i)/sumW_t(i)
418 sumRMS(i,iprot)=sumRMS_t(i)/sumW_t(i)
419 sumRGY(i,iprot)=sumRGY_t(i)/sumW_t(i)
420 sumEsq(i,iprot)=sumEsq_t(i)/sumW_t(i)
421 sumEEE(i,iprot)=sumEEE_t(i)/sumW_t(i)
422 & -3*sumE(i,iprot)*sumEsq(i,iprot)+2*sumE(i,iprot)**3
423 sumQsq(i,iprot)=sumQsq_t(i)/sumW_t(i)-sumQ(i,iprot)**2
424 sumEQ(i,iprot)=sumEQ_t(i)/sumW_t(i)
425 sumEEQ(i,iprot)=sumEEQ_t(i)/sumW_t(i)-
426 & sumEsq(i,iprot)*sumQ(i,iprot)-
427 & 2*sumEQ(i,iprot)*sumE(i,iprot)+
428 & 2*sumQ(i,iprot)*sumE(i,iprot)**2
430 write (iout,*) "Temper",temper," RgasT",RgasT
431 write (iout,*) "sumE",sumE(i,iprot)," sumEsq",sumEsq(i,iprot),
432 & " sumEbis",sumEbis(i,iprot)
434 sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
435 sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
436 sumW(i,iprot)=-dlog(sumW_t(i))*RgasT+elowest_all
437 sumE(i,iprot)=sumE(i,iprot)+elowest_all
438 write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
440 & sumE(i,iprot),(sumE(i,iprot)-sumW(i,iprot))/RgasT
441 & +entbase_all,sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
442 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
443 & sumQsq(i,iprot),sumEQ(i,iprot),sumEEE(i,iprot),
446 write (iout,*) "zvtot",
447 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot)
454 open (istat,file=nazwa)
455 write (istat,'(/a)') '# Thermal characteristics of folding'
457 Temper = GridT0+i*deltaT
459 sumE(i,iprot)=sumE(i,iprot)/sumW(i,iprot)
460 sumEbis(i,iprot)=Temper*sumE(i,iprot)/sumW(i,iprot)
461 sumQ(i,iprot)=sumQ(i,iprot)/sumW(i,iprot)
462 sumRMS(i,iprot)=sumRMS(i,iprot)/sumW(i,iprot)
463 sumRGY(i,iprot)=sumRGY(i,iprot)/sumW(i,iprot)
464 sumEsq(i,iprot)=sumEsq(i,iprot)/sumW(i,iprot)
465 sumQsq(i,iprot)=sumQsq(i,iprot)/sumW(i,iprot)-sumQ(i,iprot)**2
466 sumEQ(i,iprot)=sumEQ(i,iprot)/sumW(i,iprot)
467 sumEEQ(i,iprot)=sumEEQ(i,iprot)/sumW(i,iprot)-
468 & sumEsq(i,iprot)*sumQ(i,iprot)-
469 & 2*sumEQ(i,iprot)*sumE(i,iprot)+
470 & 2*sumQ(i,iprot)*sumE(i,iprot)**2
471 sumEQ(i,iprot)=sumEQ(i,iprot)-sumQ(i,iprot)*sumE(i,iprot)
472 sumEsq(i,iprot)=sumEsq(i,iprot)-sumE(i,iprot)**2
473 sumW(i,iprot)=-dlog(sumW(i,iprot))*RgasT+elowest_all
474 sumE(i,iprot)=sumE(i,iprot)+elowest_all
475 write (istat,'(f7.1,3f15.5,3f10.5,5e15.5)') Temper,
478 & (sumE(i,iprot)-sumW(i,iprot))/RgasT+entbase_all,
479 & sumQ(i,iprot),sumRMS(i,iprot),sumRGY(i,iprot),
480 & sumEsq(i,iprot)/(RgasT*Temper)-sumEbis(i,iprot),
482 & sumEQ(i,iprot),sumEEE(i,iprot),sumEEQ(i,iprot)