subroutine proc_data(nvarr,x,*) 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" double precision fac0,t0 double precision fT(5),ftprim(5),ftbis(5),quot,quotl,quotl1, & denom,kfac,kfacl,logfac,eplus,eminus,tanhT,entfacmax, & sument,sumentsq,aux real csingle(3,maxres2) 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 logical non_conv integer wykl integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng, & inn,in,iref,ib,ibatch,num,itmp,imin integer l1,l2,ncon_all double precision rcyfr 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 x(max_paropt) double precision wconf(maxstr_proc) double precision dx,dy,dz,dxref,dyref,dzref, & dtheta,dphi,rmsthet,rmsphi,rmsside double precision pinorm external pinorm #ifdef MPI double precision e_total_T_(maxstr_proc) #endif logical LPRN integer maxl1,maxl2 double precision rms,rmsave,rmsmin,rmsnat,qwolynes,qwolyness T0=3.0d2 fac0=1.0d0/(T0*Rgas) kfac=2.4d0 write (iout,*) "fac0",fac0 write(iout,*) "torsional and valence angle mode",tor_mode c c Determine the number of variables and the initial vector of optimizable c parameters c lprn=.false. if (restart_flag) then call w2x(nvarr,x_orig,*11) do i=1,n_ene ww_oorig(i)=ww(i) enddo open (88,file=prefix(:ilen(prefix))//'.restin',status='unknown') read(88,*,end=99,err=99) (x(i),i=1,nvarr) close(88) write (iout,*) write (iout,*) & "Variables replaced with values from the restart file." do i=1,nvarr write (iout,'(i5,f15.5)') i,x(i) enddo write (iout,*) call x2w(nvarr,x) goto 88 99 write (iout,*) & "Error - restart file nonexistent, incompatible or damaged." #ifdef MPI call MPI_Finalize(Ierror) #endif stop & "Error - restart file nonexistent, incompatible or damaged." else call w2x(nvarr,x,*11) write (iout,*) "PROC: after W2X nvarr",nvarr do i=1,n_ene ww_oorig(i)=ww(i) enddo do i=1,nvarr x_orig(i)=x(i) enddo endif 88 continue c c Setup weights for UNRES c wsc=ww(1) wscp=ww(2) welec=ww(3) wcorr=ww(4) wcorr5=ww(5) wcorr6=ww(6) wel_loc=ww(7) wturn3=ww(8) wturn4=ww(9) wturn6=ww(10) wang=ww(11) wscloc=ww(12) wtor=ww(13) wtor_d=ww(14) wstrain=ww(15) wvdwpp=ww(16) wbond=ww(17) wsccor=ww(19) #ifdef SCP14 scal14=ww(17)/ww(2) #endif do i=1,n_ene weights(i)=ww(i) enddo DO IPROT=1,NPROT call restore_molinfo(iprot) DO IBATCH=1,natlike(iprot)+2 c Calculate temperature-dependent weight scaling factors do ib=1,nbeta(ibatch,iprot) fac=betaT(ib,ibatch,iprot) quot=fac0/fac if (rescale_mode.eq.0) then do l=1,5 fT(l)=1.0d0 fTprim(l)=0.0d0 fTbis(l)=0.0d0 enddo else if (rescale_mode.eq.1) then quotl=1.0d0 kfacl=1.0d0 do l=1,5 quotl1=quotl quotl=quotl*quot kfacl=kfacl*kfac denom=kfacl-1.0d0+quotl fT(l)=kfacl/denom ftprim(l)=-l*ft(l)*quotl1/(T0*denom) ftbis(l)=l*kfacl*quotl1* & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3) write (iout,*) "ib",ib," beta",betaT(ib,ibatch,iprot), & " quot",quot," l",l,fT(l),fTprim(l),fTbis(l) enddo else if (rescale_mode.eq.2) then quotl=1.0d0 do l=1,5 quotl1=quotl quotl=quotl*quot eplus=dexp(quotl) eminus=dexp(-quotl) logfac=1.0d0/dlog(eplus+eminus) tanhT=(eplus-eminus)/(eplus+eminus) fT(l)=1.12692801104297249644d0*logfac ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- & 2*l*quotl1/T0*logfac* & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) & +ftprim(l)*tanhT) enddo else write (iout,*) "Unknown RESCALE_MODE",rescale_mode call flush(iout) stop endif do k=1,n_ene escal(k,ib,ibatch,iprot)=1.0d0 escalprim(k,ib,ibatch,iprot)=0.0d0 escalbis(k,ib,ibatch,iprot)=0.0d0 enddo escal(3,ib,ibatch,iprot)=ft(1) escal(4,ib,ibatch,iprot)=ft(3) escal(5,ib,ibatch,iprot)=ft(4) escal(6,ib,ibatch,iprot)=ft(5) escal(7,ib,ibatch,iprot)=ft(2) escal(8,ib,ibatch,iprot)=ft(2) escal(9,ib,ibatch,iprot)=ft(3) escal(10,ib,ibatch,iprot)=ft(5) escal(13,ib,ibatch,iprot)=ft(1) escal(14,ib,ibatch,iprot)=ft(2) escal(19,ib,ibatch,iprot)=ft(1) escalprim(3,ib,ibatch,iprot)=ftprim(1) escalprim(4,ib,ibatch,iprot)=ftprim(3) escalprim(5,ib,ibatch,iprot)=ftprim(4) escalprim(6,ib,ibatch,iprot)=ftprim(5) escalprim(7,ib,ibatch,iprot)=ftprim(2) escalprim(8,ib,ibatch,iprot)=ftprim(2) escalprim(9,ib,ibatch,iprot)=ftprim(3) escalprim(10,ib,ibatch,iprot)=ftprim(5) escalprim(13,ib,ibatch,iprot)=ftprim(1) escalprim(14,ib,ibatch,iprot)=ftprim(2) escalprim(19,ib,ibatch,iprot)=ftprim(1) escalbis(3,ib,ibatch,iprot)=ftbis(1) escalbis(4,ib,ibatch,iprot)=ftbis(3) escalbis(5,ib,ibatch,iprot)=ftbis(4) escalbis(6,ib,ibatch,iprot)=ftbis(5) escalbis(7,ib,ibatch,iprot)=ftbis(2) escalbis(8,ib,ibatch,iprot)=ftbis(2) escalbis(9,ib,ibatch,iprot)=ftbis(3) escalbis(10,ib,ibatch,iprot)=ftbis(5) escalbis(13,ib,ibatch,iprot)=ftbis(1) escalbis(14,ib,ibatch,iprot)=ftbis(2) escalbis(19,ib,ibatch,iprot)=ftbis(1) enddo betmin(ibatch,iprot)=betaT(1,ibatch,iprot) betmax(ibatch,iprot)=betaT(1,ibatch,iprot) do ib=2,nbeta(ibatch,iprot) if (betaT(ib,ibatch,iprot).lt.betmin(ibatch,iprot)) & betmin(ibatch,iprot)=betaT(ib,ibatch,iprot) if (betaT(ib,ibatch,iprot).gt.betmax(ibatch,iprot)) & betmax(ibatch,iprot)=betaT(ib,ibatch,iprot) enddo write (iout,'(2(a,i5,2x),2(a,f7.2,2x))') & "Protein",iprot," batch",ibatch, & " betmin",betmin(ibatch,iprot)," betmax",betmax(ibatch,iprot) ENDDO ! IBATCH ENDDO ! IPROT c c Compute energy and entropy histograms to filter out conformations c with potntially too low or too high weights. c call determine_histograms c c Make the working list of conformations c call make_list(.true.,.true.,nvarr,x) c 3/8/2013 Calculate RMSDs/Qs and the Ptab array 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 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 if (rmscomp(iprot)) then call fitsq(rms,c(1,nstart_sup(iprot)), & creff(1,nstart_sup(iprot)), & nsup(iprot),przes,obrot,non_conv) if (non_conv) then print *,'Error: FITSQ non-convergent, iref',iref rms=1.0d2 else if (rms.lt.-1.0d-6) then print *,'Error: rms^2 = ',rms,iref rms = 1.0d2 else if (rms.lt.0) then rms=0.0d0 else rms = dsqrt(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) 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 rmsthet = 0.0d0 do k=nnt+2,nct dtheta = theta(k)-thetareff(k) rmsthet = rmsthet+dtheta*dtheta enddo rmsthet=rmsthet/(nct-nnt-2) rmsphi = 0.0d0 do k=nnt+3,nct dphi = pinorm(phi(k)-phireff(k)) rmsphi = rmsphi + dphi*dphi enddo rmsphi=rmsphi/(nct-nnt-3) aux=aux*dexp(-0.5d0*(rmsthet+rmsphi)/sigmaang2(iprot)**2) endif if (sidecomp(iprot)) then call sc_loc_geom(.false.) rmsside = 0.0d0 do k=nnt+1,nct-1 dxref = xxtab(k)-xxreff(k) dyref = yytab(k)-yyreff(k) dzref = zztab(k)-zzreff(k) rmsside = rmsside + dx*dx+dy*dy+dz*dz enddo rmsside=rmsside/(nct-nnt+1) 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 #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) 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 = 0.0d0 do k=nnt+2,nct dtheta = theta(k)-theta_ref(k,iref,ib,iprot) c write (iout,*) k,theta(k),theta_ref(k,iref,ib,iprot), c & dtheta rmsthet = rmsthet+dtheta*dtheta enddo rmsthet=rmsthet/(nct-nnt-2) rmsphi = 0.0d0 do k=nnt+3,nct dphi = pinorm(phi(k)-phi_ref(k,iref,ib,iprot)) c write (iout,*) k,phi(k),phi_ref(k,iref,ib,iprot), c & pinorm(phi(k)-phi_ref(k,iref,ib,iprot)) rmsphi = rmsphi + dphi*dphi enddo rmsphi=rmsphi/(nct-nnt-3) endif if (sidecomp(iprot)) then call sc_loc_geom(.false.) rmsside = 0.0d0 do k=nnt+1,nct-1 dxref = xxtab(k)-xxref(k,iref,ib,iprot) dyref = yytab(k)-yyref(k,iref,ib,iprot) dzref = zztab(k)-zzref(k,iref,ib,iprot) rmsside = rmsside + dx*dx+dy*dy+dz*dz enddo rmsside=rmsside/(nct-nnt+1) endif #ifdef DEBUG write (iout,*) "ii",ii," jj",jj," iref",iref, & " rms",rms, & "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+rmsphi)/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 write (iout,*) return 11 return1 end c------------------------------------------------------------------------------- subroutine determine_histograms implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" integer maxind parameter (maxind=1000) #ifdef MPI include "mpif.h" integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag include "COMMON.MPI" #endif include "COMMON.COMPAR" include "COMMON.ENERGIES" include "COMMON.VMCPAR" include "COMMON.WEIGHTS" include "COMMON.PROTNAME" include "COMMON.IOUNITS" include "COMMON.FFIELD" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.CONTROL" include "COMMON.CLASSES" integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt double precision eminh,emaxh,tminh,tmaxh integer he(0:maxind), & ht(0:maxind),hemax,htmax c c Construct energy and entropy histograms c write (iout,*) "Calling DETERMINE HISTOGRAMS" call flush(iout) do iprot=1,nprot eminh=1.0d20 emaxh=-1.0d20 tminh=1.0d20 tmaxh=-1.0d20 do i=0,maxind he(i)=0 ht(i)=0 enddo do i=1,ntot(iprot) iii=i if (eini(iii,iprot).lt.eminh) & eminh=eini(iii,iprot) if (eini(iii,iprot).gt.emaxh) & emaxh=eini(iii,iprot) if (entfac(iii,iprot).lt.tminh) & tminh=entfac(iii,iprot) if (entfac(iii,iprot).gt.tmaxh) & tmaxh=entfac(iii,iprot) enddo c write (2,*) "DETERMINE HISTOGRAMS 1" call flush(iout) do i=1,ntot(iprot) iii=i inde=eini(iii,iprot)-eminh indt=entfac(iii,iprot)-tminh if (inde.le.maxind) he(inde)=he(inde)+1 if (indt.le.maxind) ht(indt)=ht(indt)+1 enddo c write (2,*) "DETERMINE HISTOGRAMS 2" idelte=emaxh-eminh ideltt=tmaxh-tminh hemax=0 c write (iout,*) "idelte",idelte," ideltt",ideltt call flush(iout) do i=0,idelte c write (iout,*) "i",i," he",he(i)," hemax",hemax call flush(iout) if (he(i).gt.hemax) hemax=he(i) enddo htmax=0 do i=0,ideltt c write (iout,*) "i",i," ht",ht(i)," htmax",htmax call flush(iout) if (ht(i).gt.htmax) htmax=ht(i) enddo c write (2,*) "DETERMINE HISTOGRAMS 3" call flush(iout) hemax=hemax/hefac(iprot) if (hemax.lt.hemax_low(iprot)) & hemax=hemax_low(iprot) htmax=htmax/htfac(iprot) if (htmax.lt.htmax_low(iprot)) & htmax=htmax_low(iprot) i=0 do while (i.lt.idelte .and. he(i).lt.hemax) i=i+1 enddo c write (2,*) "DETERMINE HISTOGRAMS 4" call flush(iout) e_lowb(iprot)=eminh+i i=0 do while (i.lt.ideltt .and. ht(i).lt.htmax) i=i+1 enddo t_lowb(iprot)=tminh+i #ifdef DEBUG write (iout,*) "protein",iprot write (iout,*) "energy histogram" do i=0,idelte write (iout,'(f15.5,i10)') eminh+i,he(i) enddo write (iout,*) "entropy histogram" do i=0,ideltt write (iout,'(f15.5,i10)') tminh+i,ht(i) enddo write (iout,*) "emin",eminh," tmin",tminh, & " hemax",hemax," htmax",htmax, & " e_lowb",e_lowb(iprot), & " t_lowb",t_lowb(iprot) call flush(iout) #endif enddo return end c------------------------------------------------------------------------------ subroutine imysort(n, x, ipermut) implicit none integer n integer x(n),xtemp integer ipermut(n) integer i,j,imax,ipm do i=1,n xtemp=x(i) imax=i do j=i+1,n if (x(j).lt.xtemp) then imax=j xtemp=x(j) endif enddo x(imax)=x(i) x(i)=xtemp ipm=ipermut(imax) ipermut(imax)=ipermut(i) ipermut(i)=ipm enddo return end c-------------------------------------------------------------------------- subroutine write_conf_count implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CLASSES" include "COMMON.COMPAR" include "COMMON.VMCPAR" include "COMMON.PROTNAME" include "COMMON.IOUNITS" include "COMMON.PROTFILES" integer iprot,ibatch,i,ii,j,kk integer ilen external ilen logical LPRN write (iout,*) write (iout,*)"Numbers of conformations after applying the cutoff" write (iout,*) do iprot=1,nprot write (iout,'(a,2x,a,i10)') "Protein", & protname(iprot)(:ilen(protname(iprot))), & ntot_work(iprot) enddo return end c------------------------------------------------------------------------- double precision function fdum() fdum=0.0D0 return end