subroutine ptab_calc implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag include "COMMON.MPI" #endif include "COMMON.CONTROL" include "COMMON.ENERGIES" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.WEIGHTS" include "COMMON.PROTNAME" include "COMMON.IOUNITS" include "COMMON.FFIELD" include "COMMON.PROTFILES" include "COMMON.COMPAR" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.THERMAL" include "COMMON.WEIGHTDER" include "COMMON.EREF" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.VAR" include "COMMON.TIME1" include "COMMON.FMATCH" double precision fac0,t0 real csingle(3,maxres2) double precision entfacmax,sument,sumentsq,aux double precision creff(3,maxres2),cc(3,maxres2) double precision thetareff(maxres2),phireff(maxres2), & xxreff(maxres2),yyreff(maxres2),zzreff(maxres2) double precision przes(3),obrot(3,3),etoti integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng, & inn,in,iref,ib,ibatch,num,itmp,imin,iperm integer jj,istart_conf,iend_conf,ipass_conf,iii integer ilen,iroof external ilen,iroof double precision ej,emaxj,sumP,sumP_t double precision fac double precision dx,dy,dz,dxref,dyref,dzref, & dtheta,dphi,rmsthet,rmsphi,rmsside #ifdef MPI double precision e_total_T_(maxstr_proc) #endif double precision rms,rmsave,rmsmin,rmsnat,rmscalc,rmscalc_thet, & rmscalc_phi,rmscalc_side,qwolynes,qwolyness do iprot=1,nprot call restore_molinfo(iprot) open (icbase,file=bprotfiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec(iprot)) entfacmax=entfac(1,iprot) sument=entfac(1,iprot) sumentsq=entfac(1,iprot)**2 do i=2,ntot_work(iprot) if (entfac(i,iprot).gt.entfacmax) entfacmax=entfac(i,iprot) sument=sument+entfac(i,iprot) sumentsq=sumentsq+entfac(i,iprot)**2 enddo sument=sument/ntot_work(iprot) sumentsq=dsqrt(sumentsq/ntot_work(iprot)-sument**2) write (2,*) "entfacmax",entfacmax," entave",sument, & " stdent",sumentsq #ifdef WEIDIST #ifdef MPI nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc) ipass_conf=0 jj=0 sumP=0.0d0 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc ipass_conf=ipass_conf+1 write (iout,*) "Pass",ipass_conf istart_conf=i iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot)) #else nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc) ipass_conf=0 sumP=0.0d0 do i=1,ntot(iprot),maxstr_proc ipass_conf=ipass_conf+1 write (iout,*) "Pass",ipass_conf istart_conf=i iend_conf=min0(i+maxstr_proc-1,ntot(iprot)) #endif c c Read the chunk of conformations off a DA scratchfile. c call daread_ccoords(iprot,istart_conf,iend_conf) IF (GEOM_AND_WHAM_WEIGHTS) THEN c Compute energies do ib=1,nbeta(1,iprot) write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot) ii=0 do iii=istart_conf,iend_conf ii=ii+1 kk=tsil(iii,iprot) if (kk.eq.0) cycle etoti=0.0d0 do k=1,n_ene etoti=etoti+ww(k)*escal(k,ib,1,iprot) & *enetb(ii,k,iprot) enddo e_total_T(iii,ib)=entfac(kk,iprot) & +betaT(ib,1,iprot)*(etoti-elowest(ib,1,iprot)) enddo enddo ! ib #ifdef DEBUG write (iout,*) "Energies before Gather" do iii=1,ntot_work(iprot) write (iout,'(i5,100E12.5)') iii, & (e_total_T(iii,ib),ib=1,nbeta(1,iprot)) enddo #endif #ifdef MPI do ib=1,nbeta(1,iprot) do k=1,scount(me1,iprot) e_total_T_(k)=e_total_T(indstart(me1,iprot)+k-1,ib) enddo c call MPI_AllGatherv(e_total_T(indstart(me1,iprot),ib), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION, c & Comm1, IERROR) call MPI_AllGatherv(e_total_T_(1), & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION, & Comm1, IERROR) enddo ! ib #ifdef DEBUG write (iout,*) "Energies after Gather" do iii=1,ntot_work(iprot) write (iout,'(i5,100E12.5)') iii, & (e_total_T(iii,ib),ib=1,nbeta(1,iprot)) enddo #endif #endif ENDIF ! GEOM_AND_WHAM_WEIGHTS ii=0 c Compute distances c write (iout,*) "Start calculating weirms",istart_conf,iend_conf do iii=istart_conf,iend_conf ii=ii+1 k=tsil(iii,iprot) c write (iout,*) "Calculating weirms iii",iii," ii",ii," k",k if (k.eq.0) cycle call restore_ccoords(iprot,ii) c Calculate the rmsds of the current conformations and convert them into weights c to approximate equal weighting of all points of the conformational space do k=1,2*nres do l=1,3 creff(l,k)=c(l,k) enddo enddo c 7/6/16 AL Add angle comparison if (anglecomp(iprot) .or. sidecomp(iprot)) then call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) do k=1,nres phireff(k)=phi(k) thetareff(k)=theta(k) xxreff(k)=xxtab(k) yyreff(k)=yytab(k) zzreff(k)=zztab(k) enddo endif do ib=1,nbeta(1,iprot) weirms(iii,ib)=0.0d0 enddo C-----------TEST c caonly(iprot)=.true. C-----------TEST do iref=1,ntot_work(iprot) if (iref.ne.iii) then read(icbase,rec=iref) ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k),l=1,3),k=nnt+nres,nct+nres), & nss,(ihpb(k),jhpb(k),k=1,nss),eini(iref,iprot), & entfac(iref,iprot), & ((nu(k,l,ii,iprot),k=1,maxdimnat),l=1,natlike(iprot)) do k=1,2*nres do l=1,3 c(l,k)=csingle(l,k) enddo enddo #ifdef DEBUG do k=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3), & (creff(l,k),l=1,3) enddo #endif rms = 0.0d0 iperm=1 if (rmscomp(iprot)) then rms=rmscalc(c(1,1),creff(1,1),iref,iii,iprot,iperm) #ifdef DEBUG write (iout,*) "iref",iref," iii",iii," rms",rms #endif #ifdef DEBUG write (2,*) "iii",iii," iref=",iref," ib",ib," etotal", & e_total_T(iref,ib), & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib), & " rms",rms," fac",dexp(-0.5d0*(rms/sigma2(iprot))**2), & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)), & dexp(-0.5d0*(rms/sigma2(iprot))**2)* & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)) #endif aux=dexp(-0.5d0*(rms/sigma2(iprot))**2) #ifdef DEBUG write (iout,*)"iref",iref," iii",iii," rms",rms," aux",aux #endif else rms=qwolyness(creff,iprot) #ifdef DEBUG write (2,*) "iii",iii," iref=",iref," ib",ib," etotal", & e_total_T(iref,ib), & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib), & " rms",-dlog(rms)," fac",rms, & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)), & rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)) #endif aux=rms endif ! rmscomp c 7/6/16 AL Add angle comparison if (anglecomp(iprot).or.sidecomp(iprot)) & call int_from_cart(.true.,.false.) if (anglecomp(iprot)) then c write (iout,*) "jj",jj," iref",iref rmsthet=rmscalc_thet(theta(1),thetareff(1),iperm) rmsphi=rmscalc_phi(phi(1),phireff(1),iperm) #ifdef DEBUG write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi #endif aux=aux*dexp(-0.5d0*(rmsthet**2+rmsphi**2)/ & sigmaang2(iprot)**2) #ifdef DEBUG write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi, & " auxang",aux #endif endif if (sidecomp(iprot)) then call sc_loc_geom(.false.) rmsside=rmscalc_side(xxtab(1),yytab(1),zztab(1), & xxreff(1),yyreff(1),zzreff(1),iperm) #ifdef DEBUG write (iout,*) "rmsside",rmsside #endif aux=aux*dexp(-0.5d0*(rmsside/sigmaside2(iprot))**2) endif else aux=1.0d0 endif do ib=1,nbeta(1,iprot) if (GEOM_AND_WHAM_WEIGHTS) then weirms(iii,ib) = weirms(iii,ib) & +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)) else weirms(iii,ib) = weirms(iii,ib)+aux endif enddo ! ib enddo ! iref enddo ! iii C-----------TEST c caonly(iprot)=.false. C-----------TEST #ifdef DEBUG write (iout,*) "weirms" do iii=istart_conf,iend_conf write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1, & nbeta(1,iprot)) enddo #endif enddo ! i #endif #ifdef MPI nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc) DO IB=1,NBETA(1,IPROT) #ifdef OUT_PTAB write(csa_bank,'(a,f4.0,4hPtab,i4.4)') & protname(iprot)(:ilen(protname(iprot))), & 1.0d0/(0.001987*betaT(ib,1,iprot)),me open (icsa_bank,file=csa_bank,status="unknown") #endif ipass_conf=0 jj=0 sumP=0.0d0 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc ipass_conf=ipass_conf+1 write (iout,*) "Pass",ipass_conf istart_conf=i iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot)) #else nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc) ipass_conf=0 sumP=0.0d0 do i=1,ntot(iprot),maxstr_proc ipass_conf=ipass_conf+1 write (iout,*) "Pass",ipass_conf istart_conf=i iend_conf=min0(i+maxstr_proc-1,ntot(iprot)) #endif ii=0 c c Read the chunk of conformations off a DA scratchfile. c call daread_ccoords(iprot,istart_conf,iend_conf) write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot) do iii=istart_conf,iend_conf c if (ib.eq.1) then rmsave=0.0d0 rmsmin=1.0d10 c endif ii=ii+1 jj=jj+1 k=tsil(iii,iprot) call restore_ccoords(iprot,ii) c write (iout,*) "zeroing Ptab",iii,ii,jj Ptab(jj,ib,iprot)=0.0d0 if (rmscomp(iprot)) then do iref=1,nref(ib,iprot) #ifdef DEBUG itmp=ipdb ipdb=iout call pdbout(0.0d0,"conf") c write (2,*) "nstart_sup",nstart_sup(iprot), c & " nend_sup",nend_sup(iprot) c do k=nstart_sup(iprot),nend_sup(iprot) c write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3), c & (cref(l,k,iref,ib,iprot),l=1,3) c enddo c do k=nstart_sup(iprot),nend_sup(iprot) c write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3), c & (cref(l,k+nres,iref,ib,iprot),l=1,3) c enddo do k=1,2*nres do l=1,3 cc(l,k)=c(l,k) c(l,k)=cref(l,k,iref,ib,iprot) enddo enddo call pdbout(0.0d0,"ref") do k=1,2*nres do l=1,3 c(l,k)=cc(l,k) enddo enddo ipdb=itmp #endif rms=rmsnat(ib,jj,iref,iprot,iperm) #ifdef DEBUG write (iout,*) iii,iref," rmsnat",rms #endif c 7/6/16 AL Add angle comparison c write(iout,*) "anglecomp",anglecomp if (anglecomp(iprot).or.sidecomp(iprot)) & call int_from_cart(.true.,.false.) if (anglecomp(iprot)) then c write (iout,*) "jj",jj," iref",iref rmsthet=rmscalc_thet(theta(1), & theta_ref(1,iref,ib,iprot),iperm) rmsphi=rmscalc_phi(phi(1),phi_ref(1,iref,ib,iprot), & iperm) #ifdef DEBUG write (iout,*) "rmsthet",rmsthet," rmsphi",rmsphi #endif endif if (sidecomp(iprot)) then call sc_loc_geom(.false.) rmsside=rmscalc_side(xxtab(1),yytab(1),zztab(1), & xxref(1,iref,ib,iprot),yyref(1,iref,ib,iprot), & zzref(1,iref,ib,iprot),iperm) #ifdef DEBUG write (iout,*) "rmsside",rmsside #endif endif #ifdef DEBUG write (iout,*) "ii",ii," jj",jj," iref",iref, & " rms",rms," rmsthet",rmsthet," rmsphi",rmsphi, & " efree",efreeref(iref,ib,iprot) #endif c if (ib.eq.1) then rmsave=rmsave+rms if (rms.lt.rmsmin) then imin=iref rmsmin=rms c endif endif Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot) & +dexp(-efreeref(iref,ib,iprot) & -0.5d0*(rms/sigma2(iprot))**2) if (anglecomp(iprot)) & Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)* & dexp(-0.5d0*(rmsthet**2+rmsphi**2)/sigmaang2(iprot)**2) if (sidecomp(iprot)) & Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)* & dexp(-0.5d0*(rmsside/sigmaside2(iprot))**2) enddo c if (ib.eq.1) then rmsave=rmsave/nref(ib,iprot) c endif #if defined(WEIENT) c AL 1/24/2014 Correct the weight of a conformation to represent energy-free c search of conformational space. aux=dmin1(entfacmax-entfac(k,iprot),60.0d0) Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux) #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k, c & " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot) & "entfac",entfac(k,iprot)," weirms",weirms(iii,ib), & " Ptab",Ptab(jj,ib,iprot) #endif #elif defined(WEIDIST) c write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot), c & " weirms",weirms(iii,ib)," Ptab/weirms", c & Ptab(jj,ib,iprot)/weirms(iii,ib) #ifdef OUT_PTAB write (icsa_bank,'(2i6,3f10.5,2f8.4)') & iii,imin,-dlog(Ptab(jj,ib,iprot)), & -dlog(weirms(iii,ib)), & -dlog(Ptab(jj,ib,iprot)/weirms(iii,ib)), & rmsave,rmsmin #endif if (ib.eq.1) rmstb(iii,iprot)=rmsmin Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib) #elif defined(WEICLUST) Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot) #endif sumP=sumP+Ptab(jj,ib,iprot) else do iref=1,nref(ib,iprot) rms = qwolynes(ib,iref,iprot) #ifdef DEBUG write (iout,*) "iii",iii," ii",ii," jj",jj, & "iref",iref," rms",rms,-dlog(rms), & "efree",efreeref(iref,ib,iprot) #endif Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms* & *dexp(-efreeref(iref,ib,iprot)) enddo #if defined(WEIENT) c AL 1/24/2014 Correct the weight of a conformation to represent energy-free c search of conformational space. aux=dmin1(entfacmax-entfac(k,iprot),60.0d0) Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux) #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k, & " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot) #endif #elif defined(WEICLUST) Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot) #elif defined(WEIDIST) write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot), & " weirms",weirms(iii,ib)," Ptab/weirms", & Ptab(jj,ib,iprot)/weirms(iii,ib) Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib) #endif sumP=sumP+Ptab(jj,ib,iprot) endif #ifdef DEBUG write (iout,*) "ii",ii," jj",jj," ib",ib," iprot",iprot, & " Ptab",Ptab(jj,ib,iprot)," sumP",sumP #endif enddo ! iii enddo ! i #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib," sumP before REDUCE",sumP #endif #ifdef MPI call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION, & MPI_SUM,Comm1,IERROR) sumP=sumP_t #endif #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib," sumP after REDUCE",sumP #endif jj=0 do iii=indstart(me1,iprot),indend(me1,iprot) jj=jj+1 #ifdef DEBUG write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot) #endif Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP #ifdef DEBUG write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot) #endif enddo #ifdef OUT_PTAB close(icsa_bank) #endif ENDDO ! IB close(icbase) enddo ! iprot #ifdef DEBUG write (iout,*) "ELOWEST at the end of PROC_DATA" do iprot=1,nprot do ibatch=1,natlike(iprot)+2 do ib=1,nbeta(ibatch,iprot) write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib, & " elowest",elowest(ib,ibatch,iprot) enddo enddo enddo write (iout,*) "Ptab at the end of PROC" do iprot=1,nprot do ib=1,nbeta(1,iprot) write (iout,*) "Protein",iprot," beta",ib jj=0 do i=indstart(me1,iprot),indend(me1,iprot) jj=jj+1 write (iout,*) iprot,ib,ii,jj,Ptab(jj,ib,iprot) enddo enddo enddo #endif return end