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" include "COMMON.TIME1" include "COMMON.FMATCH" 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 Fmean,SStot 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 cutoffeval=.false. 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 if (.not.maxlik(iprot)) cycle 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 if (shield_mode.gt.0) then escal(1,ib,ibatch,iprot)=ft(1) escal(2,ib,ibatch,iprot)=ft(1) escal(16,ib,ibatch,iprot)=ft(1) escalprim(1,ib,ibatch,iprot)=ft(1) escalprim(2,ib,ibatch,iprot)=ft(1) escalprim(16,ib,ibatch,iprot)=ft(1) escalbis(1,ib,ibatch,iprot)=ft(1) escalbis(2,ib,ibatch,iprot)=ft(1) escalbis(16,ib,ibatch,iprot)=ft(1) endif 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) call make_list_MD(.true.) c 3/8/2013 Calculate RMSDs/Qs and the Ptab array call PTab_calc c AL 5/13/2019 Compute variances of MD forces do iprot=1,nprot if (.not.fmatch(iprot)) cycle open (icbase,file=bprotfiles_MD(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_MD(iprot)) Fmean=0.0d0 do ii=1,ntot_work_MD(iprot),maxstr_proc istart_conf=ii iend_conf=min0(ii+maxstr_proc-1,ntot_MD(iprot)) call daread_MDcoords(iprot,istart_conf,iend_conf) do j=istart_conf,iend_conf c jj=list_conf_MD(j,iprot) c write (iout,*) "iprot",iprot," j",j," jj",jj do i=1,nsite(iprot) do k=1,3 Fmean=Fmean+CGforce_MD(k,i,j,iprot) enddo enddo enddo enddo Fmean=Fmean/(3*nsite(iprot)*ntot_work_MD(iprot)) write(iout,*) "PROC_DATA Protein",iprot," Fmean",Fmean SStot=0.0d0 do j=1,ntot_work_MD(iprot) c jj=list_conf_MD(j,iprot) do i=1,nsite(iprot) do k=1,3 SStot=SStot+(CGforce_MD(k,i,j,iprot)-Fmean)**2 enddo enddo enddo varforce(iprot)=SStot/(3*nsite(iprot)*ntot_work_MD(iprot)) meanforce(iprot)=Fmean write (iout,*) "PROC_DATA Protein",iprot," varforce", & varforce(iprot)," mean force",meanforce(iprot) close(icbase) enddo ! iprot 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 if (.not.maxlik(iprot)) cycle 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