subroutine averages(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.COMPAR" include "COMMON.CLASSES" C Define local variables integer i,ii,iii,inat,jj,j,k,kk,l,ik,jk,iprot,ib integer ipass_conf,istart_conf,iend_conf double precision energia(0:max_ene) double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk double precision aux,fac,ef1,ef2,em1,em2,var1,var2 double precision gnmr,tcpu_ini,tcpu logical lprn integer icant external icant lprn=.true. c c Zero out tables. c c Mazimum likeliood do ib=1,nbeta(1,iprot) do k=1,n_ene sumlikder(k,ib,iprot)=0.0d0 sumqder(k,ib,iprot)=0.0d0 enddo do k=1,nntyp sumlikeps(k,ib,iprot)=0.0d0 sumqeps(k,ib,iprot)=0.0d0 enddo do k=1,maxtor_var sumliktor(k,ib,iprot)=0.0d0 sumqtor(k,ib,iprot)=0.0d0 enddo sumlik(ib,iprot)=0.0d0 efree(ib,1,iprot)=0.0d0 enddo c Heat capacity do ib=1,nbeta(2,iprot) do k=1,n_ene eave_pftot(k,ib,iprot)=0.0d0 eave_pfprimtot(k,ib,iprot)=0.0d0 eave_pfbistot(k,ib,iprot)=0.0d0 emix_pftot(k,ib,iprot)=0.0d0 emix_pfprimtot(k,ib,iprot)=0.0d0 emix_pfbistot(k,ib,iprot)=0.0d0 emixsq_pftot(k,ib,iprot)=0.0d0 enddo do k=1,nntyp enepsave_ftot(k,ib,iprot)=0.0d0 eneps_mix_ftot(k,ib,iprot)=0.0d0 eneps_mix_fbistot(k,ib,iprot)=0.0d0 eneps_mixsq_ftot(k,ib,iprot)=0.0d0 enddo do k=1,maxtor_var entorave_ftot(k,ib,iprot)=0.0d0 entorave_fprimtot(k,ib,iprot)=0.0d0 entorave_fbistot(k,ib,iprot)=0.0d0 entor_mix_ftot(k,ib,iprot)=0.0d0 entor_mix_fprimtot(k,ib,iprot)=0.0d0 entor_mix_fbistot(k,ib,iprot)=0.0d0 entor_mixsq_ftot(k,ib,iprot)=0.0d0 enddo emean_ftot(ib,iprot)=0.0d0 ebis_ftot(ib,iprot)=0.0d0 esquare_ftot(ib,iprot)=0.0d0 enddo c Conformation-dependent averages do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) do l=1,natdim(inat,iprot) do k=1,n_ene nu_pf(k,l,ib,inat,iprot)=0.0d0 enddo do k=1,nntyp nuepsave_f(k,l,ib,inat,iprot)=0.0d0 enddo do k=1,maxtor_var nutorave_f(k,l,ib,inat,iprot)=0.0d0 enddo nuave(l,ib,inat,iprot)=0.0d0 enddo enddo enddo c c Calculate the contributions to averages from each processor or all c contributions if calculations are run in uniprocessor mode. c The derivatives of energy in epsilons are dumped to disk, if needed. c #ifdef DEBUG write (iout,*) "Protein",iprot," nchunk_conf",nchunk_conf(iprot) #endif IF (NCHUNK_CONF(IPROT).EQ.1) THEN #ifdef MPI do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=indstart(me1,iprot),indend(me1,iprot) ii=ii+1 i2ii(i,iprot)=ii enddo #else do i=1,ntot_work(iprot) i2ii(i,iprot)=i endif #endif #ifdef DEBUG do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif call ave_eval(iprot) ELSE open (ientin,file=benefiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) ipass_conf=0 #ifdef MPI do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot)) #else do istart_conf=1,ntot_work(iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot)) #endif c c Read the chunk of energies and derivatives off a DA scratchfile. c ipass_conf=ipass_conf+1 do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=istart_conf,iend_conf ii=ii+1 i2ii(i,iprot)=ii enddo #ifdef DEBUG write (iout,*) "ipass_conf",ipass_conf, & " istart_conf",istart_conf," iend_conf",iend_conf do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif call daread_ene(iprot,istart_conf,iend_conf) call ave_eval(iprot) enddo close(ientin) ENDIF c c Put tohether all contributions. c call ave_sum(iprot) return end c------------------------------------------------------------------------------- subroutine ave_eval(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.THERMAL" include "COMMON.VAR" include "COMMON.CLASSES" include "COMMON.COMPAR" include "COMMON.TORSION" include "COMMON.SCCOR" include "COMMON.PROTNAME" C Define local variables integer i,ii,iii,jj,j,k,kk,l,ll,m,ik,jk,iprot,ib,inat, & ic,iblock integer ilen external ilen integer ipass_conf,istart_conf,iend_conf double precision energia(0:max_ene) double precision etoti,sumpart,esquarei,emeani, & elowesti,enepsjk,eprim,ebis,eprimi,ebisi,etoti_orig, & entorjk,eave_pft(max_ene),emix_pft(max_ene), & esquare_ft,efree_t,emixsq_pft(max_ene), & eneps_mixt(nntyp),eneps_meant(nntyp), & enepsave_ft(nntyp),eneps_mix_ft(nntyp), & eneps_mixsq_ft(nntyp),emean_ft double precision aux,auxf,fac,facf,ef1,ef2,em1,em2,var1,var2, & deltei,deltei_orig,temper,elowest_all double precision gnmr,tcpu_ini,tcpu double precision ftune_epsprim external ftune_epsprim logical lprn integer icant external icant lprn=.false. C Compute the likelihood sum DO IB=1,NBETA(1,IPROT) elowest_all=elowest(ib,1,iprot) #ifdef OUT_LIK write(csa_bank,'(a,f4.0,4hMlik,i3.3)') & protname(iprot)(:ilen(protname(iprot))), & 1.0d0/(0.001987*betaT(ib,1,iprot)),me open (icsa_bank,file=csa_bank,status="unknown") #endif #ifdef DEBUG write(iout,*) "iprot",iprot," ib",ib, & " elowest",elowest_l(ib,iprot) #endif fac=betaT(ib,1,iprot) temper=1.0d0/(fac*Rgas) sumpart=0.0d0 sumlik(ib,iprot)=0.0d0 do k=1,n_ene sumlikder(k,ib,iprot)=0.0d0 sumqder(k,ib,iprot)=0.0d0 enddo do j=1,ntot_work(iprot) i = tsil(j,iprot) #ifdef DEBUG if (i.gt.0) write (iout,*) "i",i," iprot",iprot, & " indstart",indstart(me1,iprot), & " indend",indend(me1,iprot) #endif if (i.eq.0) goto 221 jj = i2ii(i,iprot) if (jj.gt.0) then c Temperature-dependent energy etoti=0.0d0 do k=1,n_ene etoti=etoti+ww(k)*escal(k,ib,1,iprot) & *enetb(jj,k,iprot) enddo deltei=etoti-elowest(ib,1,iprot) #ifdef DEBUG write (iout,*) "etoti",etoti," deltei",deltei write (iout,'(20f8.3)') (ww(k),k=1,n_ene) write (iout,'(20f8.3)') (enetb(jj,k,iprot), & k=1,n_ene) #endif aux=entfac(i,iprot)+fac*deltei sumlik(ib,iprot)=sumlik(ib,iprot)+aux*Ptab(jj,ib,iprot) #ifdef DEBUG write (iout,'(2i5,f5.2,7(a,e15.5))') & i,ib,fac," e_total",etoti, & " eini",eini(i,iprot)," entfac",entfac(i,iprot), & " e_lowest",elowest(ib,1,iprot)," aux",aux write (iout,'(i5,5f12.5)') & i,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot), & aux*Ptab(jj,ib,iprot) #endif #ifdef OUT_LIK write (icsa_bank,'(2i5,4f12.5,e15.5)') & i,jj,etoti,aux,dlog(Ptab(jj,ib,iprot)),rmstb(i,iprot), & Ptab(jj,ib,iprot) #endif do k=1,n_ene sumlikder(k,ib,iprot)=sumlikder(k,ib,iprot)+ & enetb(jj,k,iprot)*escal(k,ib,1,iprot)*Ptab(jj,ib,iprot) enddo k=0 do ik=1,ntyp do jk=1,ik k=k+1 sumlikeps(k,ib,iprot)=sumlikeps(k,ib,iprot)+ & (ftune_epsprim(eps(ik,jk))* & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot)) & *Ptab(jj,ib,iprot)*ww(1) enddo enddo do k=1,nbacktor_var sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot) & +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(13,ib,1,iprot) enddo do k=nbacktor_var+1,ntor_var sumliktor(k,ib,iprot)=sumliktor(k,ib,iprot) & +Ptab(jj,ib,iprot)*entor(k,jj,iprot)*escal(19,ib,1,iprot) enddo if (aux.le.150.0) then aux=dexp(-aux) sumpart=sumpart+aux do k=1,n_ene sumqder(k,ib,iprot)=sumqder(k,ib,iprot) & +aux*enetb(jj,k,iprot)*escal(k,ib,1,iprot) enddo k=0 do ik=1,ntyp do jk=1,ik k=k+1 sumqeps(k,ib,iprot)=sumqeps(k,ib,iprot)+ & ww(1)*(ftune_epsprim(eps(ik,jk))* & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot))*aux enddo enddo do k=1,nbacktor_var sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot) & +aux*entor(k,jj,iprot)*escal(13,ib,1,iprot) enddo do k=nbacktor_var+1,ntor_var sumqtor(k,ib,iprot)=sumqtor(k,ib,iprot) & +aux*entor(k,jj,iprot)*escal(19,ib,1,iprot) enddo endif endif 221 continue enddo ! j efree(ib,1,iprot)=sumpart #ifdef DEBUG write (iout,*) " ib",ib," iprot",iprot, & " sumlik",sumlik(ib,iprot), & " sumq",efree(ib,1,iprot) #endif #ifdef OUT_LIK close(icsa_bank) #endif ENDDO ! ib C Calculation of heat capacity DO IB=1,NBETA(2,IPROT) elowest_all=elowest(ib,2,iprot) #ifdef DEBUG write(iout,*) "iprot",iprot," ib",ib, & " elowest",elowest(ib,2,iprot) write(iout,*) "escalbis",(escalbis(k,ib,2,iprot),k=1,n_ene) #endif fac=betaT(ib,2,iprot) temper=1.0d0/(fac*Rgas) sumpart=0.0d0 emeani=0.0d0 esquarei=0.0d0 ebisi=0.0d0 #ifdef DEBUG write (iout,*) "Before sum, ib=",ib write (iout,*) "efree",efree(ib,2,iprot),sumpart, & " emean",emean_ftot(ib,iprot), & " esquare",esquare_ftot(ib,iprot), & " ebis",ebis_ftot(ib,iprot), & " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene) #endif do j=1,ntot_work(iprot) i = tsil(j,iprot) #ifdef DEBUG if (i.gt.0) write (iout,*) "i",i," iprot",iprot, & " indstart",indstart(me1,iprot), & " indend",indend(me1,iprot) #endif if (i.eq.0) cycle jj = i2ii(i,iprot) #ifdef DEBUG if (jj.gt.0) then write (iout,*) "ib",ib," j",j, write (iout,*) "nu",(nu(k,jj,iprot),k=1, & natconstr(iprot)) endif #endif if (jj.gt.0) then c Temperature-dependent energy etoti=0.0d0 etoti_orig=0.0d0 eprim=0.0d0 ebis=0.0d0 do k=1,n_ene etoti=etoti+ww(k)*escal(k,ib,2,iprot) & *enetb(jj,k,iprot) eprim=eprim+ww(k)*escalprim(k,ib,2,iprot) & *enetb(jj,k,iprot) ebis=ebis+ww(k)*escalbis(k,ib,2,iprot) & *enetb(jj,k,iprot) enddo deltei=etoti-elowest(ib,2,iprot) #ifdef DEBUG write (iout,*) "etoti",etoti," deltei",deltei write (iout,'(20f8.3)') (ww(k),k=1,n_ene) write (iout,'(20f8.3)') (enetb(jj,k,iprot), & k=1,n_ene) #endif aux=entfac(i,iprot)+fac*deltei #ifdef DEBUG write (iout,'(f5.2,7(a,e15.5))') & fac," e_total",etoti, & " eini",eini(i,iprot)," entfac",entfac(i,iprot), & " eprim",eprim," ebis",ebis, & " e_lowest",elowest(ib,2,iprot)," aux",aux #endif if (aux.le.150.0) then aux=dexp(-aux) sumpart=sumpart+aux etoti=etoti-temper*eprim emeani=emeani+etoti*aux esquarei=esquarei+aux*etoti**2 ebisi=ebisi+aux*ebis do k=1,n_ene eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*escal(k,ib,2,iprot) c write (iout,*) "eave_pf",eave_pf(k,ii,ib,iprot) eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*(escal(k,ib,2,iprot)- & temper*escalprim(k,ib,2,iprot)) eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*escalbis(k,ib,2,iprot) emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*etoti*escal(k,ib,2,iprot) emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*etoti*(escal(k,ib,2,iprot) & -temper*escalprim(k,ib,2,iprot)) emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*ebis*escal(k,ib,2,iprot) emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot) & +aux*enetb(jj,k,iprot)*etoti**2*escal(k,ib,2,iprot) enddo k=0 do ik=1,ntyp do jk=1,ik k=k+1 enepsjk=ftune_epsprim(eps(ik,jk))* & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot) enepsjk=enepsjk*ww(1) enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot) & +aux*enepsjk eneps_mix_ftot(k,ib,iprot)= & eneps_mix_ftot(k,ib,iprot)+aux*enepsjk*etoti eneps_mix_fbistot(k,ib,iprot)= & eneps_mix_fbistot(k,ib,iprot)+aux*enepsjk*ebis eneps_mixsq_ftot(k,ib,iprot)= & eneps_mixsq_ftot(k,ib,iprot)+aux*enepsjk*etoti**2 enddo enddo do k=1,nbacktor_var entorjk=entor(k,jj,iprot) c write (iout,*) " k"," ik",ik," jk",jk, c & " entor",entorjk," contirb", c & aux*entorjk*escal(13,ib,2,iprot) entorave_ftot(k,ib,iprot)= & entorave_ftot(k,ib,iprot) & +aux*entorjk*escal(13,ib,2,iprot) c write (iout,*) "entorave_f", c & entorave_ftot(k,ib,iprot) entorave_fprimtot(k,ib,iprot)= & entorave_fprimtot(k,ib,iprot) & +aux*entorjk*(escal(13,ib,2,iprot)- & temper*escalprim(13,ib,2,iprot)) entorave_fbistot(k,ib,iprot)= & entorave_fbistot(k,ib,iprot) & +aux*entorjk*escalbis(13,ib,2,iprot) entor_mix_ftot(k,ib,iprot)= & entor_mix_ftot(k,ib,iprot) & +aux*entorjk*etoti*escal(13,ib,2,iprot) entor_mix_fprimtot(k,ib,iprot)= & entor_mix_fprimtot(k,ib,iprot) & +aux*entorjk*etoti*(escal(13,ib,2,iprot)- & temper*escalprim(13,ib,2,iprot)) entor_mix_fbistot(k,ib,iprot)= & entor_mix_fbistot(k,ib,iprot) & +aux*entorjk*ebis*escal(13,ib,2,iprot) entor_mixsq_ftot(k,ib,iprot)= & entor_mixsq_ftot(k,ib,iprot) & +aux*entorjk*etoti**2*escal(13,ib,2,iprot) enddo do k=nbacktor_var+1,ntor_var entorjk=entor(k,jj,iprot) c write (iout,*) " k"," ik",ik," jk",jk, c & " entor",entorjk," contirb", c & aux*entorjk*escal(19,ib,2,iprot) entorave_ftot(k,ib,iprot)= & entorave_ftot(k,ib,iprot) & +aux*entorjk*escal(19,ib,2,iprot) c write (iout,*) "entorave_f", c & entorave_ftot(k,ib,iprot) entorave_fprimtot(k,ib,iprot)= & entorave_fprimtot(k,ib,iprot) & +aux*entorjk*(escal(19,ib,2,iprot)- & temper*escalprim(19,ib,2,iprot)) entorave_fbistot(k,ib,iprot)= & entorave_fbistot(k,ib,iprot) & +aux*entorjk*escalbis(19,ib,2,iprot) entor_mix_ftot(k,ib,iprot)= & entor_mix_ftot(k,ib,iprot) & +aux*entorjk*etoti*escal(19,ib,2,iprot) entor_mix_fprimtot(k,ib,iprot)= & entor_mix_fprimtot(k,ib,iprot) & +aux*entorjk*etoti*(escal(19,ib,2,iprot)- & temper*escalprim(13,ib,2,iprot)) entor_mix_fbistot(k,ib,iprot)= & entor_mix_fbistot(k,ib,iprot) & +aux*entorjk*ebis*escal(19,ib,2,iprot) entor_mixsq_ftot(k,ib,iprot)= & entor_mixsq_ftot(k,ib,iprot) & +aux*entorjk*etoti**2*escal(19,ib,2,iprot) enddo endif endif enddo efree(ib,2,iprot)=sumpart emean_ftot(ib,iprot)=emeani ebis_ftot(ib,iprot)=ebisi esquare_ftot(ib,iprot)=esquarei #ifdef DEBUG write (iout,*) "After sum, ib=",ib write (iout,*) "efree",efree(ib,2,iprot),sumpart, & " emean",emean_ftot(ib,iprot), & " esquare",esquare_ftot(ib,iprot), & " ebis",ebis_ftot(ib,iprot), & " emix_pfprim",(emix_pfprimtot(k,ib,iprot),k=1,n_ene) #endif ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)*temper do k=1,n_ene eave_pfbistot(k,ib,iprot)= & eave_pfbistot(k,ib,iprot)*temper emix_pfbistot(k,ib,iprot)= & emix_pfbistot(k,ib,iprot)*temper enddo do k=1,nntyp eneps_mix_fbistot(k,ib,iprot)= & eneps_mix_fbistot(k,ib,iprot)*temper enddo do k=1,ntor_var entorave_fbistot(k,ib,iprot)= & entorave_fbistot(k,ib,iprot)*temper entor_mix_fbistot(k,ib,iprot)= & entor_mix_fbistot(k,ib,iprot)*temper enddo #ifdef DEBUG write (iout,*) "eave_pf",(eave_pftot(k,ib,iprot), & k=1,n_ene) write (iout,*) "entorave_f",(entorave_ftot(k,ib,iprot), & k=1,ntor_var) #endif #ifdef DEBUG write (iout,*) "ib",ib," temper",temper, & " ebis",ebis_ftot(ib,iprot) #endif ENDDO ! ib C Calculation of conformation-dependent averages DO INAT=1,NATLIKE(IPROT) DO IB=1,NBETA(I+2,IPROT) elowest_all=elowest(ib,inat+2,iprot) #ifdef DEBUG write(iout,*) "iprot",iprot," ib",ib, & " elowest",elowest(ib,iprot) #endif fac=betaT(ib,2,iprot) temper=1.0d0/(fac*Rgas) emeani=0.0d0 do j=1,ntot_work(iprot) i = tsil(j,iprot) #ifdef DEBUG if (i.gt.0) write (iout,*) "i",i," iprot",iprot, & " indstart",indstart(me1,iprot), & " indend",indend(me1,iprot) #endif if (i.eq.0) cycle jj = i2ii(i,iprot) #ifdef DEBUG if (jj.gt.0) then write (iout,*) "ib",ib," j",j, write (iout,*) "nu",(nu(k,jj,iprot),k=1, & natconstr(iprot)) endif #endif if (jj.gt.0) then c Temperature-dependent energy etoti=0.0d0 do k=1,n_ene etoti=etoti+ww(k)*escal(k,ib,inat+2,iprot) & *enetb(jj,k,iprot) enddo deltei=etoti-elowest(ib,inat+2,iprot) #ifdef DEBUG write (iout,*) "etoti",etoti," deltei",deltei write (iout,'(20f8.3)') (ww(k),k=1,n_ene) write (iout,'(20f8.3)') (enetb(jj,k,iprot), & k=1,n_ene) #endif aux=entfac(i,iprot)+fac*deltei #ifdef DEBUG write (iout,'(f5.2,7(a,e15.5))') & fac," e_total",etoti, & " eini",eini(i,iprot)," entfac",entfac(i,iprot), & " eprim",eprim," ebis",ebis, & " e_lowest",elowest(ib,inat+2,iprot)," aux",aux #endif if (aux.le.150.0) then aux=dexp(-aux) sumpart=sumpart+aux c 4/13/04 AL Components of the conformation-dependent averages do k=1,n_ene eave_nat_pftot(k,ib,inat,iprot)= & eave_nat_pftot(k,ib,inat,iprot) & +aux*enetb(jj,k,iprot)*escal(k,ib,inat+2,iprot) enddo k=0 do ik=1,ntyp do jk=1,ik k=k+1 enepsjk=ftune_epsprim(eps(ik,jk))* & eneps(1,k,jj,iprot)+eneps(2,k,jj,iprot) enepsjk=enepsjk*ww(1) enepsave_nat_ftot(k,ib,inat,iprot)= & enepsave_nat_ftot(k,ib,inat,iprot) & +aux*enepsjk enddo enddo k=0 do k=1,nbacktor_var entorjk=entor(k,jj,iprot) c write (iout,*) " k", c & " entor",entorjk," contirb", c & aux*entorjk*escal(13,ib,inat+2,iprot) entorave_nat_ftot(k,ib,inat,iprot)= & entorave_nat_ftot(k,ib,inat,iprot) & +aux*entorjk*escal(13,ib,inat+2,iprot) enddo do k=nbacktor_var+1,ntor_var do l=1,3 entorjk=entor(k,jj,iprot) c write (iout,*) " k", c & " entor",entorjk," contirb", c & aux*entorjk*escal(19,ib,inat+2,iprot) entorave_nat_ftot(k,ib,inat,iprot)= & entorave_nat_ftot(k,ib,inat,iprot) & +aux*entorjk*escal(19,ib,inat+2,iprot) enddo enddo do l=1,natdim(inat,iprot) nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot) & +aux*nu(l,inat,jj,iprot) do k=1,n_ene nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot) & +aux*enetb(jj,k,iprot)*nu(k,inat,jj,iprot)* & escal(k,ib,inat+2,iprot) enddo k=0 do ik=1,ntyp do jk=1,ik k=k+1 nuepsave_f(k,l,ib,inat,iprot)= & nuepsave_f(k,l,ib,inat,iprot)+aux*enepsjk* & nu(l,inat,jj,iprot) enddo enddo k=0 do ik=0,ntortyp-1 do jk=-ntortyp,ntortyp do iblock=1,2 if (mask_tor(0,jk,ik,iblock).gt.0) then k=k+1 nutorave_f(k,l,ib,inat,iprot)= & nutorave_f(k,l,ib,inat,iprot) & +aux*entorjk*escal(13,ib,inat+2,iprot)* & nu(l,inat,jj,iprot) endif enddo enddo enddo do ik=1,nsccortyp do jk=1,nsccortyp do ll=1,3 if (mask_tor(ll,jk,ik,1).gt.0) then k=k+1 nutorave_f(k,l,ib,inat,iprot)= & nutorave_f(k,l,ib,inat,iprot) & +aux*entorjk*escal(13,ib,inat+2,iprot)* & nu(l,inat,jj,iprot) endif enddo enddo enddo enddo endif endif #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib, write (iout,*) "nuave" write (iout,'(20f10.5)') (nuave(k,ib,iprot), & k=1,natconstr(iprot)) #endif #ifdef DEBUG write (iout,*) "inat",inat, & " efree_nat",efree_nat(ib,inat,iprot) #endif enddo #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib write (iout,*) "nuave0" write (iout,'(20f10.5)') (nuave(k,ib,inat,iprot), & k=1,natdim(inat,iprot)) #endif ENDDO ! ib ENDDO ! INAT return end c------------------------------------------------------------------------------- subroutine ave_sum(iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.VAR" include "COMMON.COMPAR" C Define local variables integer i,ii,iii,jj,j,k,l,ik,jk,iprot,ib,inat integer ipass_conf,istart_conf,iend_conf double precision energia(0:max_ene) double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk, & eave_pft(max_ene),emix_pft(max_ene),eave_pfprimt(max_ene), & eave_pfbist(max_ene),emix_pfprimt(max_ene),emix_pfbist(max_ene), & esquare_ft,efree_t, & emixsq_pft(max_ene),eneps_mixt(nntyp),eneps_meant(nntyp), & enepsave_ft(nntyp),eneps_mix_ft(nntyp),entorave_ft(maxtor_var), & entor_mix_ft(maxtor_var), & eneps_mix_fbist(nntyp),entorave_fbist(maxtor_var), & entorave_fprimt(maxtor_var),entor_mix_fprimt(maxtor_var), & entor_mix_fbist(maxtor_var),eneps_mixsq_ft(nntyp), & entor_mixsq_ft(maxtor_var),emean_ft,ebis_ft,nuave_t(maxdimnat), & nu_pft(max_ene,maxdimnat), & nuepsave_ft(nntyp,maxdimnat), & nutorave_ft(maxtor_var,maxdimnat), & sumlik_t,sumlikder_t(max_ene),sumlikeps_t(nntyp), & sumliktor_t(maxtor_var), & sumq_t,sumqder_t(max_ene),sumqeps_t(nntyp),sumqtor_t(maxtor_var) double precision aux,fac,ef1,ef2,em1,em2,var1,var2,efree_tot,facF, & elowest_all double precision gnmr,tcpu_ini,tcpu logical lprn integer icant external icant lprn=.true. c Maximum likelihood contribution DO IB=1,NBETA(1,IPROT) fac=betaT(ib,1,iprot) #ifdef MPI #ifdef DEBUG write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3" write (iout,*) "iprot",iprot," ib",ib write (iout,*) "sumlik",sumlik(ib,iprot), & " sumq",efree(ib,1,iprot) call flush(iout) #endif call MPI_Reduce( sumlik(ib,iprot), sumlik_t,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce( efree(ib,1,iprot), sumq_t,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) #ifdef DEBUG write (iout,*) "sumlik",sumlik(ib,iprot)," sumlik_t",sumlik_t write (iout,*) "sumq",efree(ib,1,iprot)," sumq_t",sumq_t call flush(iout) #endif call MPI_Reduce(sumlikder(1,ib,iprot),sumlikder_t(1), & n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(sumlikeps(1,ib,iprot), & sumlikeps_t(1), nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(sumliktor(1,ib,iprot), & sumliktor_t(1), ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(sumqder(1,ib,iprot),sumqder_t(1),n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(sumqeps(1,ib,iprot),sumqeps_t(1), nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(sumqtor(1,ib,iprot),sumqtor_t(1), ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) #ifdef DEBUG write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3" call flush(iout) #endif if (me1.eq.master) then #ifdef DEBUG write (iout,*) "ib",ib, & " elowest",elowest(ib,1,iprot), & " sumlik",sumlik_t," qsum",sumq_t," fac",fac #endif do k=1,n_ene sumlikder(k,ib,iprot)=fac*( & sumlikder_t(k) & -sumqder_t(k)/sumq_t) enddo do k=1,nntyp sumlikeps(k,ib,iprot)=fac*( & sumlikeps_t(k)-sumqeps_t(k)/sumq_t) enddo c write (iout,*) "eavetor",eave_pftot(13,ib,iprot), c & eave_pftot(19,ib,iprot) do k=1,ntor_var sumliktor(k,ib,iprot)=fac*( & sumliktor_t(k)-sumqtor_t(k)/sumq_t) enddo sumlik(ib,iprot)=sumlik_t+dlog(sumq_t) efree(ib,1,iprot)=sumq_t #ifdef DEBUG write (iout,*) "ib",ib," iprot",iprot," final sumlik", & sumlik(ib,iprot)," sumq",efree(ib,1,iprot) #endif endif #else do k=1,n_ene sumlikder(k,ib,iprot)=fac*( & sumlikder(k,ib,iprot) & -sumqder(k,ib,iprot)/efree(ib,1,iprot)) enddo do k=1,nntyp sumlikeps(k,ib,iprot)=fac*( & sumlikeps(k,ib,iprot) & -sumqeps(k,ib,iprot)/efree(ib,1,iprot)) enddo do k=1,ntor_var sumliktor(k,ib,iprot)=fac*( & sumliktor(k,ib,iprot) & -sumqtor(k,ib,iprot)/efree(ib,1,iprot)) enddo sumlik(ib,iprot)=sumlik(ib,iprot)+dlog(efree(ib,1,iprot) c & -elowest(nbeta(iprot)+ib,iprot)*fac #endif ii=0 do k=1,n_ene if (imask(k).gt.0) then ii=ii+1 sumlikder(ii,ib,iprot)=sumlikder(k,ib,iprot) endif enddo ENDDO ! ib c Heat capacity and averages DO IB=1,NBETA(2,IPROT) fac=betaT(ib,2,iprot) #ifdef MPI #ifdef DEBUG write (iout,*) "Processor",me,me1," calling MPI_Reduce: 3" write (iout,*) "iprot",iprot," ib",ib call flush(iout) #endif call MPI_Reduce( efree(ib,2,iprot), efree_t,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) #ifdef DEBUG write (iout,*) "efree",efree(iprot write (iout,*) "efree_t",efree_t call flush(iout) #endif call MPI_Reduce(emean_ftot(ib,iprot),emean_ft,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(ebis_ftot(ib,iprot),ebis_ft,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(esquare_ftot(ib,iprot),esquare_ft,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eave_pftot(1,ib,iprot),eave_pft, & n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eave_pfprimtot(1,ib,iprot), & eave_pfprimt,n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eave_pfbistot(1,ib,iprot), & eave_pfbist,n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) c write (iout,*) "eave_pf",(eave_pf(k,iprot), c & k=1,n_ene) c write (iout,*) "eave_pft",(eave_pft(k),k=1,n_ene) call MPI_Reduce(emix_pftot(1,ib,iprot),emix_pft, & n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(emix_pfprimtot(1,ib,iprot), & emix_pfprimt,n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(emix_pfbistot(1,ib,iprot), & emix_pfbist,n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(emixsq_pftot(1,ib,iprot), & emixsq_pft, n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(enepsave_ftot(1,ib,iprot), & enepsave_ft, nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eneps_mix_ftot(1,ib,iprot), & eneps_mix_ft,nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eneps_mix_fbistot(1,ib,iprot), & eneps_mix_fbist,nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(eneps_mixsq_ftot(1,ib,iprot), & eneps_mixsq_ft,nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) c write (iout,*) "enepsave_f",(enepsave_f(k,iprot), c & k=1,nntyp) c write (iout,*) "enepsave_ft",(enepsave_ft(k), c & k=1,nntyp) call MPI_Reduce(entorave_ftot(1,ib,iprot), & entorave_ft, ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entorave_fprimtot(1,ib,iprot), & entorave_fprimt, ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entorave_fbistot(1,ib,iprot), & entorave_fbist, ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entor_mix_ftot(1,ib,iprot), & entor_mix_ft,ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entor_mix_fprimtot(1,ib,iprot), & entor_mix_fprimt,ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entor_mix_fbistot(1,ib,iprot), & entor_mix_fbist,ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entor_mixsq_ftot(1,ib,iprot), & entor_mixsq_ft,ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) c write (iout,*) "entorave_f",(entorave_f(k,iprot), c & k=1,ntor_var) c write (iout,*) "entorave_ft",(entorave_ft(k), c & k=1,ntor_var) #ifdef DEBUG write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3" call flush(iout) #endif if (me1.eq.master) then elowest_all=elowest(ib,2,iprot) #ifdef DEBUG write (iout,*) "ib",ib, & " elowest",elowest(ib,iprot), & "efree",efree_t," fac",fac," facF",facF, & " efree_tot",efree_tot #endif do k=1,n_ene eave_pftot(k,ib,iprot)=eave_pft(k)/efree_t eave_pfprimtot(k,ib,iprot)=eave_pfprimt(k)/efree_t eave_pfbistot(k,ib,iprot)=eave_pfbist(k)/efree_t emix_pftot(k,ib,iprot)=emix_pft(k)/efree_t emix_pfprimtot(k,ib,iprot)=emix_pfprimt(k)/efree_t emix_pfbistot(k,ib,iprot)=emix_pfbist(k)/efree_t emixsq_pftot(k,ib,iprot)=emixsq_pft(k)/efree_t enddo do k=1,nntyp enepsave_ftot(k,ib,iprot)=enepsave_ft(k)/efree_t eneps_mix_ftot(k,ib,iprot)=eneps_mix_ft(k)/efree_t eneps_mix_fbistot(k,ib,iprot)=eneps_mix_fbist(k)/efree_t eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ft(k)/efree_t enddo c write (iout,*) "eavetor",eave_pftot(13,ib,iprot), c & eave_pftot(19,ib,iprot) do k=1,ntor_var entorave_ftot(k,ib,iprot)=entorave_ft(k)/ & efree_t c write (iout,*) "iprot",iprot," ib",ib," k",k, c & " entorave_ftot", entorave_ftot(k,ib,iprot) entorave_fprimtot(k,ib,iprot)=entorave_fprimt(k)/efree_t entorave_fbistot(k,ib,iprot)=entorave_fbist(k)/efree_t entor_mix_ftot(k,ib,iprot)=entor_mix_ft(k)/efree_t entor_mix_fprimtot(k,ib,iprot)=entor_mix_fprimt(k)/efree_t entor_mix_fbistot(k,ib,iprot)=entor_mix_fbist(k)/efree_t entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ft(k)/efree_t enddo emean_ftot(ib,iprot)=emean_ft/efree_t ebis_ftot(ib,iprot)=ebis_ft/efree_t esquare_ftot(ib,iprot)=esquare_ft/efree_t efree(ib,2,iprot)=-dlog(efree_t)/fac+elowest(ib,2,iprot) endif #else do k=1,n_ene eave_pftot(k,ib,iprot)=eave_pftot(k,ib,iprot) & /efree(ib,2,iprot) eave_pfprimtot(k,ib,iprot)=eave_pfprimtot(k,ib,iprot) & /efree(ib,2,iprot) eave_pfbistot(k,ib,iprot)=eave_pfbistot(k,ib,iprot) & /efree(ib,2,iprot) emix_pftot(k,ib,iprot)=emix_pftot(k,ib,iprot)/efree(ib,iprot) emix_pfprimtot(k,ib,iprot)=emix_pfprimtot(k,ib,iprot) & /efree(ib,2,iprot) emix_pfbistot(k,ib,iprot)=emix_pfbistot(k,ib,iprot) & /efree(ib,2,iprot) emixsq_pftot(k,ib,iprot)=emixsq_pftot(k,ib,iprot) & /efree(ib,2,iprot) enddo do k=1,nntyp enepsave_ftot(k,ib,iprot)=enepsave_ftot(k,ib,iprot) & /efree(ib,2,iprot) eneps_mix_ftot(k,ib,iprot)=eneps_mix_ftot(k,ib,iprot) & /efree(ib,2,iprot) eneps_mixsq_ftot(k,ib,iprot)=eneps_mixsq_ftot(k,ib,iprot) & /efree(ib,2,iprot) enddo do k=1,ntor_var entorave_ftot(k,ib,iprot)=entorave_f(k,ib,iprot) & /efree(ib,2,iprot) entor_mix_ftot(k,ib,iprot)=entor_mix_ftot(k,ib,iprot) & /efree(ib,2,iprot) entor_mixsq_ftot(k,ib,iprot)=entor_mixsq_ftot(k,ib,iprot) & /efree(ib,2,iprot) enddo emean_ftot(ib,iprot)=emean_ftot(ib,iprot)/efree(ib,2,iprot) ebis_ftot(ib,iprot)=ebis_ftot(ib,iprot)/efree(ib,2,iprot) esquare_ftot(ib,iprot)=esquare_ftot(ib,iprot)/efree(ib,2,iprot) #endif ENDDO ! IB c 4/13/04 AL Components of the correlation coefficients and their derivatives DO INAT=1,NATLIKE(IPROT) DO IB=1,NBETA(INAT+2,IPROT) #ifdef MPI call MPI_Reduce( efree(ib,inat+2,iprot), efree_t,1, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce( eave_nat_pftot(1, ib,inat,iprot), eave_pft, & n_ene, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce( enepsave_nat_ftot(1, ib,inat,iprot), & enepsave_ft, nntyp, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(entorave_nat_ftot(1,ib,inat,iprot), & entorave_ft, ntor_var, & MPI_DOUBLE_PRECISION, MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(nuave(1,ib,inat,iprot),nuave_t(1), & natdim(inat,iprot), MPI_DOUBLE_PRECISION, MPI_SUM, & Master, Comm1, IERROR) #ifdef DEBUG write (iout,*) "After REDUCE nuave",iprot,ib write (iout,'(20f10.5)') & (nuave(l,ib,iprot),l=1,natconstr(iprot)) #endif call MPI_Reduce(nu_pf(1,1,ib,inat,iprot),nu_pft(1,1), & max_ene*natdim(inat,iprot), MPI_DOUBLE_PRECISION, & MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(nuepsave_f(1,1,ib,inat,iprot), & nuepsave_ft(1,1), & nntyp*natdim(inat,iprot), MPI_DOUBLE_PRECISION, & MPI_SUM, Master, Comm1, IERROR) call MPI_Reduce(nutorave_f(1,1,ib,inat,iprot), & nutorave_ft(1,1), & maxtor_var*natdim(inat,iprot), MPI_DOUBLE_PRECISION, & MPI_SUM, Master, Comm1, IERROR) #ifdef DEBUG write (iout,*) "Processor",me,me1," finished MPI_Reduce: 3" call flush(iout) #endif if (me1.eq.master) then #ifdef DEBUG write (iout,*) "ib",ib, & " elowest",elowest(ib,iprot), & "efree",efree_t," fac",fac, & " efree_tot",efree_tot #endif do l=1,natdim(inat,iprot) do k=1,n_ene nu_pf(k,l,ib,inat,iprot)=nu_pft(k,l)/efree_t enddo do k=1,nntyp nuepsave_f(k,l,ib,inat,iprot)=nuepsave_ft(k,l)/efree_t enddo do k=1,ntor_var nutorave_f(k,l,ib,inat,iprot)=nutorave_ft(k,l)/efree_t enddo nuave(l,ib,inat,iprot)=nuave_t(l)/efree_t enddo do k=1,n_ene eave_nat_pftot(k,ib,inat,iprot)=eave_pft(k)/efree_t enddo do k=1,nntyp enepsave_nat_ftot(k,ib,inat,iprot)=enepsave_ft(k)/efree_t enddo do k=1,ntor_var entorave_nat_ftot(k,ib,inat,iprot)=entorave_ft(k)/efree_t enddo endif #else do l=1,natdim(inat,iprot) do k=1,n_ene nu_pf(k,l,ib,inat,iprot)=nu_pf(k,l,ib,inat,iprot) & /efree(ib,inat+2,iprot) enddo do k=1,nntyp nuepsave_f(k,l,ib,inat,iprot)=nuepsave_f(k,l,ib,inat,iprot) & /efree(ib,inat+2,iprot) enddo do k=1,ntor_var nutorave_ftot(k,l,ib,inat,iprot)= & nutorave_ftot(k,l,ib,inat,iprot) & /efree(ib,inat+2,iprot) enddo nuave(l,ib,inat,iprot)=nuave(l,ib,inat,iprot) & /efree(ib,inat+2,iprot) enddo ! l do k=1,n_ene eave_nat_pftot(k,ib,inat,iprot)= & eave_nat_pftot(k,ib,inat,iprot) & /efree(ib,inat+2,iprot) enddo do k=1,nntyp enepsave_nat_ftot(k,ib,inat,iprot)= & enepsave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot) enddo do k=1,ntor_var enetorave_nat_ftot(k,ib,inat,iprot)= & enetorave_nat_ftot(k,ib,inat,iprot)/efree(ib,inat+2,iprot) enddo #endif ENDDO ! IB #ifdef DEBUG write (iout,*) "ib",ib," efree_tot",efree_tot #endif ENDDO ! INAT return end