1 subroutine averages(iprot)
4 include "DIMENSIONS.ZSCOPT"
11 include "COMMON.WEIGHTS"
12 include "COMMON.WEIGHTDER"
13 include "COMMON.ENERGIES"
14 include "COMMON.IOUNITS"
15 include "COMMON.VMCPAR"
16 include "COMMON.NAMES"
17 include "COMMON.INTERACT"
18 include "COMMON.TIME1"
19 include "COMMON.CHAIN"
20 include "COMMON.PROTFILES"
21 include "COMMON.COMPAR"
22 include "COMMON.CLASSES"
23 C Define local variables
24 integer i,ii,iii,inat,jj,j,k,kk,l,ik,jk,iprot,ib
25 integer ipass_conf,istart_conf,iend_conf
26 double precision energia(0:max_ene)
27 double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk
28 double precision aux,fac,ef1,ef2,em1,em2,var1,var2
29 double precision gnmr,tcpu_ini,tcpu
38 do ib=1,nbeta(1,iprot)
40 sumlikder(k,ib,iprot)=0.0d0
41 sumqder(k,ib,iprot)=0.0d0
44 sumlikeps(k,ib,iprot)=0.0d0
45 sumqeps(k,ib,iprot)=0.0d0
48 sumliktor(k,ib,iprot)=0.0d0
49 sumqtor(k,ib,iprot)=0.0d0
51 sumlik(ib,iprot)=0.0d0
52 efree(ib,1,iprot)=0.0d0
55 do ib=1,nbeta(2,iprot)
57 eave_pftot(k,ib,iprot)=0.0d0
58 eave_pfprimtot(k,ib,iprot)=0.0d0
59 eave_pfbistot(k,ib,iprot)=0.0d0
60 emix_pftot(k,ib,iprot)=0.0d0
61 emix_pfprimtot(k,ib,iprot)=0.0d0
62 emix_pfbistot(k,ib,iprot)=0.0d0
63 emixsq_pftot(k,ib,iprot)=0.0d0
66 enepsave_ftot(k,ib,iprot)=0.0d0
67 eneps_mix_ftot(k,ib,iprot)=0.0d0
68 eneps_mix_fbistot(k,ib,iprot)=0.0d0
69 eneps_mixsq_ftot(k,ib,iprot)=0.0d0
72 entorave_ftot(k,ib,iprot)=0.0d0
73 entorave_fprimtot(k,ib,iprot)=0.0d0
74 entorave_fbistot(k,ib,iprot)=0.0d0
75 entor_mix_ftot(k,ib,iprot)=0.0d0
76 entor_mix_fprimtot(k,ib,iprot)=0.0d0
77 entor_mix_fbistot(k,ib,iprot)=0.0d0
78 entor_mixsq_ftot(k,ib,iprot)=0.0d0
80 emean_ftot(ib,iprot)=0.0d0
81 ebis_ftot(ib,iprot)=0.0d0
82 esquare_ftot(ib,iprot)=0.0d0
84 c Conformation-dependent averages
85 do inat=1,natlike(iprot)
86 do ib=1,nbeta(inat+2,iprot)
87 do l=1,natdim(inat,iprot)
89 nu_pf(k,l,ib,inat,iprot)=0.0d0
92 nuepsave_f(k,l,ib,inat,iprot)=0.0d0
95 nutorave_f(k,l,ib,inat,iprot)=0.0d0
97 nuave(l,ib,inat,iprot)=0.0d0
102 c Calculate the contributions to averages from each processor or all
103 c contributions if calculations are run in uniprocessor mode.
104 c The derivatives of energy in epsilons are dumped to disk, if needed.
107 write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot)
109 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
112 do i=1,ntot_work(iprot)
116 do i=indstart(me1,iprot),indend(me1,iprot)
121 do i=1,ntot_work(iprot)
126 do i=1,ntot_work(iprot)
127 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
135 open (ientin,file=benefiles(iprot),status="old",
136 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
139 do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
140 iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
142 do istart_conf=1,ntot_work(iprot),maxstr_proc
143 iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
146 c Read the chunk of energies and derivatives off a DA scratchfile.
148 ipass_conf=ipass_conf+1
149 do i=1,ntot_work(iprot)
153 do i=istart_conf,iend_conf
158 write (iout,*) "ipass_conf",ipass_conf,
159 & " istart_conf",istart_conf," iend_conf",iend_conf
160 do i=1,ntot_work(iprot)
161 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
165 call daread_ene(iprot,istart_conf,iend_conf)
172 c Put tohether all contributions.
177 c-------------------------------------------------------------------------------
178 subroutine ave_eval(iprot)
181 include "DIMENSIONS.ZSCOPT"
184 integer IERROR,ErrCode
188 include "COMMON.WEIGHTS"
189 include "COMMON.WEIGHTDER"
190 include "COMMON.ENERGIES"
191 include "COMMON.IOUNITS"
192 include "COMMON.VMCPAR"
193 include "COMMON.NAMES"
194 include "COMMON.INTERACT"
195 include "COMMON.TIME1"
196 include "COMMON.CHAIN"
197 include "COMMON.PROTFILES"
198 include "COMMON.THERMAL"
200 include "COMMON.CLASSES"
201 include "COMMON.COMPAR"
202 include "COMMON.TORSION"
203 include "COMMON.SCCOR"
204 include "COMMON.PROTNAME"
205 C Define local variables
206 integer i,ii,iii,jj,j,k,kk,l,ll,m,ik,jk,iprot,ib,inat,
210 integer ipass_conf,istart_conf,iend_conf
211 double precision energia(0:max_ene)
212 double precision etoti,sumpart,esquarei,emeani,
213 & elowesti,enepsjk,eprim,ebis,eprimi,ebisi,etoti_orig,
214 & entorjk,eave_pft(max_ene),emix_pft(max_ene),
215 & esquare_ft,efree_t,emixsq_pft(max_ene),
216 & eneps_mixt(nntyp),eneps_meant(nntyp),
217 & enepsave_ft(nntyp),eneps_mix_ft(nntyp),
218 & eneps_mixsq_ft(nntyp),emean_ft
219 double precision aux,auxf,fac,facf,ef1,ef2,em1,em2,var1,var2,
220 & deltei,deltei_orig,temper,elowest_all
221 double precision gnmr,tcpu_ini,tcpu
222 double precision ftune_epsprim
223 external ftune_epsprim
229 C Compute the likelihood sum
230 DO IB=1,NBETA(1,IPROT)
231 elowest_all=elowest(ib,1,iprot)
233 write(csa_bank,'(a,f4.0,4hMlik,i3.3)')
234 & protname(iprot)(:ilen(protname(iprot))),
235 & 1.0d0/(0.001987*betaT(ib,1,iprot)),me
236 open (icsa_bank,file=csa_bank,status="unknown")
239 write(iout,*) "iprot",iprot," ib",ib,
240 & " elowest",elowest_l(ib,iprot)
242 fac=betaT(ib,1,iprot)
243 temper=1.0d0/(fac*Rgas)
245 sumlik(ib,iprot)=0.0d0
247 sumlikder(k,ib,iprot)=0.0d0
248 sumqder(k,ib,iprot)=0.0d0
250 do j=1,ntot_work(iprot)
253 if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
254 & " indstart",indstart(me1,iprot),
255 & " indend",indend(me1,iprot)
260 c Temperature-dependent energy
263 etoti=etoti+ww(k)*escal(k,ib,1,iprot)
266 deltei=etoti-elowest(ib,1,iprot)
268 write (iout,*) "etoti",etoti," deltei",deltei
269 write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
270 write (iout,'(20f8.3)') (enetb(jj,k,iprot),
273 aux=entfac(i,iprot)+fac*deltei
274 sumlik(ib,iprot)=sumlik(ib,iprot)+aux*Ptab(jj,ib,iprot)
276 write (iout,'(2i5,f5.2,7(a,e15.5))')
277 & i,ib,fac," e_total",etoti,
278 & " eini",eini(i,iprot)," entfac",entfac(i,iprot),
279 & " e_lowest",elowest(ib,1,iprot)," aux",aux
280 write (iout,'(i5,5f12.5)')
281 & i,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot),
282 & aux*Ptab(jj,ib,iprot)
285 write (icsa_bank,'(2i5,4f12.5,e15.5)')
286 & i,jj,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot),
290 sumlikder(k,ib,iprot)=sumlikder(k,ib,iprot)+
291 & enetb(jj,k,iprot)*escal(k,ib,1,iprot)*Ptab(jj,ib,iprot)
297 sumlikeps(k,ib,iprot)=sumlikeps(k,ib,iprot)+
298 & (ftune_epsprim(eps(ik,jk))*
299 & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))
300 & *Ptab(jj,ib,iprot)*ww(1)
304 sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot)
305 & +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(13,ib,1,iprot)
307 do k=nbacktor_var+1,ntor_var
308 sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot)
309 & +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(19,ib,1,iprot)
311 if (aux.le.150.0) then
315 sumqder(k,ib,iprot)=sumqder(k,ib,iprot)
316 & +aux*enetb(jj,k,iprot)*escal(k,ib,1,iprot)
322 sumqeps(k,ib,iprot)=sumqeps(k,ib,iprot)+
323 & ww(1)*(ftune_epsprim(eps(ik,jk))*
324 & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))*aux
328 sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot)
329 & +aux*entor(k,jj,iprot)*escal(13,ib,1,iprot)
331 do k=nbacktor_var+1,ntor_var
332 sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot)
333 & +aux*entor(k,jj,iprot)*escal(19,ib,1,iprot)
339 efree(ib,1,iprot)=sumpart
341 write (iout,*) " ib",ib," iprot",iprot,
342 & " sumlik",sumlik(ib,iprot),
343 & " sumq",efree(ib,1,iprot)
349 C Calculation of heat capacity
350 DO IB=1,NBETA(2,IPROT)
351 elowest_all=elowest(ib,2,iprot)
353 write(iout,*) "iprot",iprot," ib",ib,
354 & " elowest",elowest(ib,2,iprot)
355 write(iout,*) "escalbis",(escalbis(k,ib,2,iprot),k=1,n_ene)
357 fac=betaT(ib,2,iprot)
358 temper=1.0d0/(fac*Rgas)
364 write (iout,*) "Before sum, ib=",ib
365 write (iout,*) "efree",efree(ib,2,iprot),sumpart,
366 & " emean",emean_ftot(ib,iprot),
367 & " esquare",esquare_ftot(ib,iprot),
368 & " ebis",ebis_ftot(ib,iprot),
369 & " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
371 do j=1,ntot_work(iprot)
374 if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
375 & " indstart",indstart(me1,iprot),
376 & " indend",indend(me1,iprot)
382 write (iout,*) "ib",ib," j",j,
383 write (iout,*) "nu",(nu(k,jj,iprot),k=1,
388 c Temperature-dependent energy
394 etoti=etoti+ww(k)*escal(k,ib,2,iprot)
396 eprim=eprim+ww(k)*escalprim(k,ib,2,iprot)
398 ebis=ebis+ww(k)*escalbis(k,ib,2,iprot)
401 deltei=etoti-elowest(ib,2,iprot)
403 write (iout,*) "etoti",etoti," deltei",deltei
404 write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
405 write (iout,'(20f8.3)') (enetb(jj,k,iprot),
408 aux=entfac(i,iprot)+fac*deltei
410 write (iout,'(f5.2,7(a,e15.5))')
411 & fac," e_total",etoti,
412 & " eini",eini(i,iprot)," entfac",entfac(i,iprot),
413 & " eprim",eprim," ebis",ebis,
414 & " e_lowest",elowest(ib,2,iprot)," aux",aux
416 if (aux.le.150.0) then
419 etoti=etoti-temper*eprim
420 emeani=emeani+etoti*aux
421 esquarei=esquarei+aux*etoti**2
424 eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
425 & +aux*enetb(jj,k,iprot)*escal(k,ib,2,iprot)
426 c write (iout,*) "eave_pf",eave_pf(k,ii,ib,iprot)
427 eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
428 & +aux*enetb(jj,k,iprot)*(escal(k,ib,2,iprot)-
429 & temper*escalprim(k,ib,2,iprot))
430 eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
431 & +aux*enetb(jj,k,iprot)*escalbis(k,ib,2,iprot)
432 emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)
433 & +aux*enetb(jj,k,iprot)*etoti*escal(k,ib,2,iprot)
434 emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
435 & +aux*enetb(jj,k,iprot)*etoti*(escal(k,ib,2,iprot)
436 & -temper*escalprim(k,ib,2,iprot))
437 emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
438 & +aux*enetb(jj,k,iprot)*ebis*escal(k,ib,2,iprot)
439 emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
440 & +aux*enetb(jj,k,iprot)*etoti**2*escal(k,ib,2,iprot)
446 enepsjk=ftune_epsprim(eps(ik,jk))*
447 & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
448 enepsjk=enepsjk*ww(1)
449 enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
451 eneps_mix_ftot(k,ib,iprot)=
452 & eneps_mix_ftot(k,ib,iprot)+aux*enepsjk*etoti
453 eneps_mix_fbistot(k,ib,iprot)=
454 & eneps_mix_fbistot(k,ib,iprot)+aux*enepsjk*ebis
455 eneps_mixsq_ftot(k,ib,iprot)=
456 & eneps_mixsq_ftot(k,ib,iprot)+aux*enepsjk*etoti**2
460 entorjk=entor(k,jj,iprot)
461 c write (iout,*) " k"," ik",ik," jk",jk,
462 c & " entor",entorjk," contirb",
463 c & aux*entorjk*escal(13,ib,2,iprot)
464 entorave_ftot(k,ib,iprot)=
465 & entorave_ftot(k,ib,iprot)
466 & +aux*entorjk*escal(13,ib,2,iprot)
467 c write (iout,*) "entorave_f",
468 c & entorave_ftot(k,ib,iprot)
469 entorave_fprimtot(k,ib,iprot)=
470 & entorave_fprimtot(k,ib,iprot)
471 & +aux*entorjk*(escal(13,ib,2,iprot)-
472 & temper*escalprim(13,ib,2,iprot))
473 entorave_fbistot(k,ib,iprot)=
474 & entorave_fbistot(k,ib,iprot)
475 & +aux*entorjk*escalbis(13,ib,2,iprot)
476 entor_mix_ftot(k,ib,iprot)=
477 & entor_mix_ftot(k,ib,iprot)
478 & +aux*entorjk*etoti*escal(13,ib,2,iprot)
479 entor_mix_fprimtot(k,ib,iprot)=
480 & entor_mix_fprimtot(k,ib,iprot)
481 & +aux*entorjk*etoti*(escal(13,ib,2,iprot)-
482 & temper*escalprim(13,ib,2,iprot))
483 entor_mix_fbistot(k,ib,iprot)=
484 & entor_mix_fbistot(k,ib,iprot)
485 & +aux*entorjk*ebis*escal(13,ib,2,iprot)
486 entor_mixsq_ftot(k,ib,iprot)=
487 & entor_mixsq_ftot(k,ib,iprot)
488 & +aux*entorjk*etoti**2*escal(13,ib,2,iprot)
490 do k=nbacktor_var+1,ntor_var
491 entorjk=entor(k,jj,iprot)
492 c write (iout,*) " k"," ik",ik," jk",jk,
493 c & " entor",entorjk," contirb",
494 c & aux*entorjk*escal(19,ib,2,iprot)
495 entorave_ftot(k,ib,iprot)=
496 & entorave_ftot(k,ib,iprot)
497 & +aux*entorjk*escal(19,ib,2,iprot)
498 c write (iout,*) "entorave_f",
499 c & entorave_ftot(k,ib,iprot)
500 entorave_fprimtot(k,ib,iprot)=
501 & entorave_fprimtot(k,ib,iprot)
502 & +aux*entorjk*(escal(19,ib,2,iprot)-
503 & temper*escalprim(19,ib,2,iprot))
504 entorave_fbistot(k,ib,iprot)=
505 & entorave_fbistot(k,ib,iprot)
506 & +aux*entorjk*escalbis(19,ib,2,iprot)
507 entor_mix_ftot(k,ib,iprot)=
508 & entor_mix_ftot(k,ib,iprot)
509 & +aux*entorjk*etoti*escal(19,ib,2,iprot)
510 entor_mix_fprimtot(k,ib,iprot)=
511 & entor_mix_fprimtot(k,ib,iprot)
512 & +aux*entorjk*etoti*(escal(19,ib,2,iprot)-
513 & temper*escalprim(13,ib,2,iprot))
514 entor_mix_fbistot(k,ib,iprot)=
515 & entor_mix_fbistot(k,ib,iprot)
516 & +aux*entorjk*ebis*escal(19,ib,2,iprot)
517 entor_mixsq_ftot(k,ib,iprot)=
518 & entor_mixsq_ftot(k,ib,iprot)
519 & +aux*entorjk*etoti**2*escal(19,ib,2,iprot)
524 efree(ib,2,iprot)=sumpart
525 emean_ftot(ib,iprot)=emeani
526 ebis_ftot(ib,iprot)=ebisi
527 esquare_ftot(ib,iprot)=esquarei
529 write (iout,*) "After sum, ib=",ib
530 write (iout,*) "efree",efree(ib,2,iprot),sumpart,
531 & " emean",emean_ftot(ib,iprot),
532 & " esquare",esquare_ftot(ib,iprot),
533 & " ebis",ebis_ftot(ib,iprot),
534 & " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene)
536 ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)*temper
538 eave_pfbistot(k,ib,iprot)=
539 & eave_pfbistot(k,ib,iprot)*temper
540 emix_pfbistot(k,ib,iprot)=
541 & emix_pfbistot(k,ib,iprot)*temper
544 eneps_mix_fbistot(k,ib,iprot)=
545 & eneps_mix_fbistot(k,ib,iprot)*temper
548 entorave_fbistot(k,ib,iprot)=
549 & entorave_fbistot(k,ib,iprot)*temper
550 entor_mix_fbistot(k,ib,iprot)=
551 & entor_mix_fbistot(k,ib,iprot)*temper
554 write (iout,*) "eave_pf",(eave_pftot(k,ib,iprot),
556 write (iout,*) "entorave_f",(entorave_ftot(k,ib,iprot),
560 write (iout,*) "ib",ib," temper",temper,
561 & " ebis",ebis_ftot(ib,iprot)
564 C Calculation of conformation-dependent averages
565 DO INAT=1,NATLIKE(IPROT)
566 DO IB=1,NBETA(I+2,IPROT)
567 elowest_all=elowest(ib,inat+2,iprot)
569 write(iout,*) "iprot",iprot," ib",ib,
570 & " elowest",elowest(ib,iprot)
572 fac=betaT(ib,2,iprot)
573 temper=1.0d0/(fac*Rgas)
575 do j=1,ntot_work(iprot)
578 if (i.gt.0) write (iout,*) "i",i," iprot",iprot,
579 & " indstart",indstart(me1,iprot),
580 & " indend",indend(me1,iprot)
586 write (iout,*) "ib",ib," j",j,
587 write (iout,*) "nu",(nu(k,jj,iprot),k=1,
592 c Temperature-dependent energy
595 etoti=etoti+ww(k)*escal(k,ib,inat+2,iprot)
598 deltei=etoti-elowest(ib,inat+2,iprot)
600 write (iout,*) "etoti",etoti," deltei",deltei
601 write (iout,'(20f8.3)') (ww(k),k=1,n_ene)
602 write (iout,'(20f8.3)') (enetb(jj,k,iprot),
605 aux=entfac(i,iprot)+fac*deltei
607 write (iout,'(f5.2,7(a,e15.5))')
608 & fac," e_total",etoti,
609 & " eini",eini(i,iprot)," entfac",entfac(i,iprot),
610 & " eprim",eprim," ebis",ebis,
611 & " e_lowest",elowest(ib,inat+2,iprot)," aux",aux
613 if (aux.le.150.0) then
616 c 4/13/04 AL Components of the conformation-dependent averages
618 eave_nat_pftot(k,ib,inat,iprot)=
619 & eave_nat_pftot(k,ib,inat,iprot)
620 & +aux*enetb(jj,k,iprot)*escal(k,ib,inat+2,iprot)
626 enepsjk=ftune_epsprim(eps(ik,jk))*
627 & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)
628 enepsjk=enepsjk*ww(1)
629 enepsave_nat_ftot(k,ib,inat,iprot)=
630 & enepsave_nat_ftot(k,ib,inat,iprot)
636 entorjk=entor(k,jj,iprot)
637 c write (iout,*) " k",
638 c & " entor",entorjk," contirb",
639 c & aux*entorjk*escal(13,ib,inat+2,iprot)
640 entorave_nat_ftot(k,ib,inat,iprot)=
641 & entorave_nat_ftot(k,ib,inat,iprot)
642 & +aux*entorjk*escal(13,ib,inat+2,iprot)
644 do k=nbacktor_var+1,ntor_var
646 entorjk=entor(k,jj,iprot)
647 c write (iout,*) " k",
648 c & " entor",entorjk," contirb",
649 c & aux*entorjk*escal(19,ib,inat+2,iprot)
650 entorave_nat_ftot(k,ib,inat,iprot)=
651 & entorave_nat_ftot(k,ib,inat,iprot)
652 & +aux*entorjk*escal(19,ib,inat+2,iprot)
655 do l=1,natdim(inat,iprot)
656 nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
657 & +aux*nu(l,inat,jj,iprot)
659 nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
660 & +aux*enetb(jj,k,iprot)*nu(k,inat,jj,iprot)*
661 & escal(k,ib,inat+2,iprot)
667 nuepsave_f(k,l,ib,inat,iprot)=
668 & nuepsave_f(k,l,ib,inat,iprot)+aux*enepsjk*
669 & nu(l,inat,jj,iprot)
674 do jk=-ntortyp,ntortyp
676 if (mask_tor(0,jk,ik,iblock).gt.0) then
678 nutorave_f(k,l,ib,inat,iprot)=
679 & nutorave_f(k,l,ib,inat,iprot)
680 & +aux*entorjk*escal(13,ib,inat+2,iprot)*
681 & nu(l,inat,jj,iprot)
689 if (mask_tor(ll,jk,ik,1).gt.0) then
691 nutorave_f(k,l,ib,inat,iprot)=
692 & nutorave_f(k,l,ib,inat,iprot)
693 & +aux*entorjk*escal(13,ib,inat+2,iprot)*
694 & nu(l,inat,jj,iprot)
703 write (iout,*) "iprot",iprot," ib",ib,
704 write (iout,*) "nuave"
705 write (iout,'(20f10.5)') (nuave(k,ib,iprot),
706 & k=1,natconstr(iprot))
709 write (iout,*) "inat",inat,
710 & " efree_nat",efree_nat(ib,inat,iprot)
714 write (iout,*) "iprot",iprot," ib",ib
715 write (iout,*) "nuave0"
716 write (iout,'(20f10.5)') (nuave(k,ib,inat,iprot),
717 & k=1,natdim(inat,iprot))
723 c-------------------------------------------------------------------------------
724 subroutine ave_sum(iprot)
727 include "DIMENSIONS.ZSCOPT"
730 integer IERROR,ErrCode
734 include "COMMON.WEIGHTS"
735 include "COMMON.WEIGHTDER"
736 include "COMMON.ENERGIES"
737 include "COMMON.IOUNITS"
738 include "COMMON.VMCPAR"
739 include "COMMON.NAMES"
740 include "COMMON.INTERACT"
741 include "COMMON.TIME1"
742 include "COMMON.CHAIN"
743 include "COMMON.PROTFILES"
745 include "COMMON.COMPAR"
746 C Define local variables
747 integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,inat
748 integer ipass_conf,istart_conf,iend_conf
749 double precision energia(0:max_ene)
750 double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk,
751 & eave_pft(max_ene),emix_pft(max_ene),eave_pfprimt(max_ene),
752 & eave_pfbist(max_ene),emix_pfprimt(max_ene),emix_pfbist(max_ene),
753 & esquare_ft,efree_t,
754 & emixsq_pft(max_ene),eneps_mixt(nntyp),eneps_meant(nntyp),
755 & enepsave_ft(nntyp),eneps_mix_ft(nntyp),entorave_ft(maxtor_var),
756 & entor_mix_ft(maxtor_var),
757 & eneps_mix_fbist(nntyp),entorave_fbist(maxtor_var),
758 & entorave_fprimt(maxtor_var),entor_mix_fprimt(maxtor_var),
759 & entor_mix_fbist(maxtor_var),eneps_mixsq_ft(nntyp),
760 & entor_mixsq_ft(maxtor_var),emean_ft,ebis_ft,nuave_t(maxdimnat),
761 & nu_pft(max_ene,maxdimnat),
762 & nuepsave_ft(nntyp,maxdimnat),
763 & nutorave_ft(maxtor_var,maxdimnat),
764 & sumlik_t,sumlikder_t(max_ene),sumlikeps_t(nntyp),
765 & sumliktor_t(maxtor_var),
766 & sumq_t,sumqder_t(max_ene),sumqeps_t(nntyp),sumqtor_t(maxtor_var)
768 double precision aux,fac,ef1,ef2,em1,em2,var1,var2,efree_tot,facF,
770 double precision gnmr,tcpu_ini,tcpu
776 c Maximum likelihood contribution
777 DO IB=1,NBETA(1,IPROT)
778 fac=betaT(ib,1,iprot)
781 write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
782 write (iout,*) "iprot",iprot," ib",ib
783 write (iout,*) "sumlik",sumlik(ib,iprot),
784 & " sumq",efree(ib,1,iprot)
787 call MPI_Reduce( sumlik(ib,iprot), sumlik_t,1,
788 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
789 call MPI_Reduce( efree(ib,1,iprot), sumq_t,1,
790 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
792 write (iout,*) "sumlik",sumlik(ib,iprot)," sumlik_t",sumlik_t
793 write (iout,*) "sumq",efree(ib,1,iprot)," sumq_t",sumq_t
796 call MPI_Reduce(sumlikder(1,ib,iprot),sumlikder_t(1),
798 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
799 call MPI_Reduce(sumlikeps(1,ib,iprot),
800 & sumlikeps_t(1), nntyp,
801 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
802 call MPI_Reduce(sumliktor(1,ib,iprot),
803 & sumliktor_t(1), ntor_var,
804 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
805 call MPI_Reduce(sumqder(1,ib,iprot),sumqder_t(1),n_ene,
806 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
807 call MPI_Reduce(sumqeps(1,ib,iprot),sumqeps_t(1), nntyp,
808 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
809 call MPI_Reduce(sumqtor(1,ib,iprot),sumqtor_t(1), ntor_var,
810 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
812 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
815 if (me1.eq.master) then
817 write (iout,*) "ib",ib,
818 & " elowest",elowest(ib,1,iprot),
819 & " sumlik",sumlik_t," qsum",sumq_t," fac",fac
822 sumlikder(k,ib,iprot)=fac*(
824 & -sumqder_t(k)/sumq_t)
827 sumlikeps(k,ib,iprot)=fac*(
828 & sumlikeps_t(k)-sumqeps_t(k)/sumq_t)
830 c write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
831 c & eave_pftot(19,ib,iprot)
833 sumliktor(k,ib,iprot)=fac*(
834 & sumliktor_t(k)-sumqtor_t(k)/sumq_t)
836 sumlik(ib,iprot)=sumlik_t+dlog(sumq_t)
837 efree(ib,1,iprot)=sumq_t
839 write (iout,*) "ib",ib," iprot",iprot," final sumlik",
840 & sumlik(ib,iprot)," sumq",efree(ib,1,iprot)
845 sumlikder(k,ib,iprot)=fac*(
846 & sumlikder(k,ib,iprot)
847 & -sumqder(k,ib,iprot)/efree(ib,1,iprot))
850 sumlikeps(k,ib,iprot)=fac*(
851 & sumlikeps(k,ib,iprot)
852 & -sumqeps(k,ib,iprot)/efree(ib,1,iprot))
855 sumliktor(k,ib,iprot)=fac*(
856 & sumliktor(k,ib,iprot)
857 & -sumqtor(k,ib,iprot)/efree(ib,1,iprot))
859 sumlik(ib,iprot)=sumlik(ib,iprot)+dlog(efree(ib,1,iprot)
860 c & -elowest(nbeta(iprot)+ib,iprot)*fac
864 if (imask(k).gt.0) then
866 sumlikder(ii,ib,iprot)=sumlikder(k,ib,iprot)
870 c Heat capacity and averages
871 DO IB=1,NBETA(2,IPROT)
872 fac=betaT(ib,2,iprot)
875 write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3"
876 write (iout,*) "iprot",iprot," ib",ib
879 call MPI_Reduce( efree(ib,2,iprot), efree_t,1,
880 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
882 write (iout,*) "efree",efree(iprot
883 write (iout,*) "efree_t",efree_t
886 call MPI_Reduce(emean_ftot(ib,iprot),emean_ft,1,
887 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
888 call MPI_Reduce(ebis_ftot(ib,iprot),ebis_ft,1,
889 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
890 call MPI_Reduce(esquare_ftot(ib,iprot),esquare_ft,1,
891 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
892 call MPI_Reduce(eave_pftot(1,ib,iprot),eave_pft,
894 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
895 call MPI_Reduce(eave_pfprimtot(1,ib,iprot),
896 & eave_pfprimt,n_ene,
897 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
898 call MPI_Reduce(eave_pfbistot(1,ib,iprot),
900 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
901 c write (iout,*) "eave_pf",(eave_pf(k,iprot),
903 c write (iout,*) "eave_pft",(eave_pft(k),k=1,n_ene)
904 call MPI_Reduce(emix_pftot(1,ib,iprot),emix_pft,
906 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
907 call MPI_Reduce(emix_pfprimtot(1,ib,iprot),
908 & emix_pfprimt,n_ene,
909 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
910 call MPI_Reduce(emix_pfbistot(1,ib,iprot),
912 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
913 call MPI_Reduce(emixsq_pftot(1,ib,iprot),
915 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
916 call MPI_Reduce(enepsave_ftot(1,ib,iprot),
917 & enepsave_ft, nntyp,
918 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
919 call MPI_Reduce(eneps_mix_ftot(1,ib,iprot),
920 & eneps_mix_ft,nntyp,
921 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
922 call MPI_Reduce(eneps_mix_fbistot(1,ib,iprot),
923 & eneps_mix_fbist,nntyp,
924 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
925 call MPI_Reduce(eneps_mixsq_ftot(1,ib,iprot),
926 & eneps_mixsq_ft,nntyp,
927 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
928 c write (iout,*) "enepsave_f",(enepsave_f(k,iprot),
930 c write (iout,*) "enepsave_ft",(enepsave_ft(k),
932 call MPI_Reduce(entorave_ftot(1,ib,iprot),
933 & entorave_ft, ntor_var,
934 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
935 call MPI_Reduce(entorave_fprimtot(1,ib,iprot),
936 & entorave_fprimt, ntor_var,
937 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
938 call MPI_Reduce(entorave_fbistot(1,ib,iprot),
939 & entorave_fbist, ntor_var,
940 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
941 call MPI_Reduce(entor_mix_ftot(1,ib,iprot),
942 & entor_mix_ft,ntor_var,
943 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
944 call MPI_Reduce(entor_mix_fprimtot(1,ib,iprot),
945 & entor_mix_fprimt,ntor_var,
946 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
947 call MPI_Reduce(entor_mix_fbistot(1,ib,iprot),
948 & entor_mix_fbist,ntor_var,
949 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
950 call MPI_Reduce(entor_mixsq_ftot(1,ib,iprot),
951 & entor_mixsq_ft,ntor_var,
952 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
953 c write (iout,*) "entorave_f",(entorave_f(k,iprot),
955 c write (iout,*) "entorave_ft",(entorave_ft(k),
958 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
961 if (me1.eq.master) then
962 elowest_all=elowest(ib,2,iprot)
964 write (iout,*) "ib",ib,
965 & " elowest",elowest(ib,iprot),
966 & "efree",efree_t," fac",fac," facF",facF,
967 & " efree_tot",efree_tot
970 eave_pftot(k,ib,iprot)=eave_pft(k)/efree_t
971 eave_pfprimtot(k,ib,iprot)=eave_pfprimt(k)/efree_t
972 eave_pfbistot(k,ib,iprot)=eave_pfbist(k)/efree_t
973 emix_pftot(k,ib,iprot)=emix_pft(k)/efree_t
974 emix_pfprimtot(k,ib,iprot)=emix_pfprimt(k)/efree_t
975 emix_pfbistot(k,ib,iprot)=emix_pfbist(k)/efree_t
976 emixsq_pftot(k,ib,iprot)=emixsq_pft(k)/efree_t
979 enepsave_ftot(k,ib,iprot)=enepsave_ft(k)/efree_t
980 eneps_mix_ftot(k,ib,iprot)=eneps_mix_ft(k)/efree_t
981 eneps_mix_fbistot(k,ib,iprot)=eneps_mix_fbist(k)/efree_t
982 eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ft(k)/efree_t
984 c write (iout,*) "eavetor",eave_pftot(13,ib,iprot),
985 c & eave_pftot(19,ib,iprot)
987 entorave_ftot(k,ib,iprot)=entorave_ft(k)/
989 c write (iout,*) "iprot",iprot," ib",ib," k",k,
990 c & " entorave_ftot", entorave_ftot(k,ib,iprot)
991 entorave_fprimtot(k,ib,iprot)=entorave_fprimt(k)/efree_t
992 entorave_fbistot(k,ib,iprot)=entorave_fbist(k)/efree_t
993 entor_mix_ftot(k,ib,iprot)=entor_mix_ft(k)/efree_t
994 entor_mix_fprimtot(k,ib,iprot)=entor_mix_fprimt(k)/efree_t
995 entor_mix_fbistot(k,ib,iprot)=entor_mix_fbist(k)/efree_t
996 entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ft(k)/efree_t
998 emean_ftot(ib,iprot)=emean_ft/efree_t
999 ebis_ftot(ib,iprot)=ebis_ft/efree_t
1000 esquare_ftot(ib,iprot)=esquare_ft/efree_t
1001 efree(ib,2,iprot)=-dlog(efree_t)/fac+elowest(ib,2,iprot)
1005 eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot)
1006 & /efree(ib,2,iprot)
1007 eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot)
1008 & /efree(ib,2,iprot)
1009 eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot)
1010 & /efree(ib,2,iprot)
1011 emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)/efree(ib,iprot)
1012 emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot)
1013 & /efree(ib,2,iprot)
1014 emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot)
1015 & /efree(ib,2,iprot)
1016 emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot)
1017 & /efree(ib,2,iprot)
1020 enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot)
1021 & /efree(ib,2,iprot)
1022 eneps_mix_ftot(k,ib,iprot)=eneps_mix_ftot(k,ib,iprot)
1023 & /efree(ib,2,iprot)
1024 eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ftot(k,ib,iprot)
1025 & /efree(ib,2,iprot)
1028 entorave_ftot(k,ib,iprot)=entorave_f(k,ib,iprot)
1029 & /efree(ib,2,iprot)
1030 entor_mix_ftot(k,ib,iprot)=entor_mix_ftot(k,ib,iprot)
1031 & /efree(ib,2,iprot)
1032 entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ftot(k,ib,iprot)
1033 & /efree(ib,2,iprot)
1035 emean_ftot(ib,iprot)=emean_ftot(ib,iprot)/efree(ib,2,iprot)
1036 ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)/efree(ib,2,iprot)
1037 esquare_ftot(ib,iprot)=esquare_ftot(ib,iprot)/efree(ib,2,iprot)
1040 c 4/13/04 AL Components of the correlation coefficients and their derivatives
1041 DO INAT=1,NATLIKE(IPROT)
1042 DO IB=1,NBETA(INAT+2,IPROT)
1044 call MPI_Reduce( efree(ib,inat+2,iprot), efree_t,1,
1045 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1046 call MPI_Reduce( eave_nat_pftot(1, ib,inat,iprot), eave_pft,
1048 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1049 call MPI_Reduce( enepsave_nat_ftot(1, ib,inat,iprot),
1050 & enepsave_ft, nntyp,
1051 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1052 call MPI_Reduce(entorave_nat_ftot(1,ib,inat,iprot),
1053 & entorave_ft, ntor_var,
1054 & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR)
1055 call MPI_Reduce(nuave(1,ib,inat,iprot),nuave_t(1),
1056 & natdim(inat,iprot), MPI_DOUBLE_PRECISION, MPI_SUM,
1057 & Master, Comm1, IERROR)
1059 write (iout,*) "After REDUCE nuave",iprot,ib
1060 write (iout,'(20f10.5)')
1061 & (nuave(l,ib,iprot),l=1,natconstr(iprot))
1063 call MPI_Reduce(nu_pf(1,1,ib,inat,iprot),nu_pft(1,1),
1064 & max_ene*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1065 & MPI_SUM, Master, Comm1, IERROR)
1066 call MPI_Reduce(nuepsave_f(1,1,ib,inat,iprot),
1068 & nntyp*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1069 & MPI_SUM, Master, Comm1, IERROR)
1070 call MPI_Reduce(nutorave_f(1,1,ib,inat,iprot),
1072 & maxtor_var*natdim(inat,iprot), MPI_DOUBLE_PRECISION,
1073 & MPI_SUM, Master, Comm1, IERROR)
1075 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3"
1078 if (me1.eq.master) then
1080 write (iout,*) "ib",ib,
1081 & " elowest",elowest(ib,iprot),
1082 & "efree",efree_t," fac",fac,
1083 & " efree_tot",efree_tot
1085 do l=1,natdim(inat,iprot)
1087 nu_pf(k,l,ib,inat,iprot)=nu_pft(k,l)/efree_t
1090 nuepsave_f(k,l,ib,inat,iprot)=nuepsave_ft(k,l)/efree_t
1093 nutorave_f(k,l,ib,inat,iprot)=nutorave_ft(k,l)/efree_t
1095 nuave(l,ib,inat,iprot)=nuave_t(l)/efree_t
1098 eave_nat_pftot(k,ib,inat,iprot)=eave_pft(k)/efree_t
1101 enepsave_nat_ftot(k,ib,inat,iprot)=enepsave_ft(k)/efree_t
1104 entorave_nat_ftot(k,ib,inat,iprot)=entorave_ft(k)/efree_t
1108 do l=1,natdim(inat,iprot)
1110 nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot)
1111 & /efree(ib,inat+2,iprot)
1114 nuepsave_f(k,l,ib,inat,iprot)=nuepsave_f(k,l,ib,inat,iprot)
1115 & /efree(ib,inat+2,iprot)
1118 nutorave_ftot(k,l,ib,inat,iprot)=
1119 & nutorave_ftot(k,l,ib,inat,iprot)
1120 & /efree(ib,inat+2,iprot)
1122 nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot)
1123 & /efree(ib,inat+2,iprot)
1126 eave_nat_pftot(k,ib,inat,iprot)=
1127 & eave_nat_pftot(k,ib,inat,iprot)
1128 & /efree(ib,inat+2,iprot)
1131 enepsave_nat_ftot(k,ib,inat,iprot)=
1132 & enepsave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1135 enetorave_nat_ftot(k,ib,inat,iprot)=
1136 & enetorave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot)
1141 write (iout,*) "ib",ib," efree_tot",efree_tot