subroutine ave_property(ncon,*) ! construct the conformational ensembles at REMD temperatures implicit none include "DIMENSIONS" include "sizesclu.dat" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,errcode,status(MPI_STATUS_SIZE) #endif include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.CLUSTER" include "COMMON.PROPERTY" double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/ double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc, & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, & evdw_t,beta,w integer i,ii,ik,iproc,iscor,j,k,l,ib,ncon double precision sumprob,eini,efree,rmsdev character*80 bxname character*2 licz1 character*5 ctemper integer ilen external ilen real*8 Fmin,Fmin_all double precision energia(0:max_ene) double precision qfree(600),ave_rms_chain(600),ave_assoc(600), & ave_hbond_intra(600),ave_hbond_inter(600),ave_assoc_hbond(600) #ifdef MPI double precision qfree_part(600),ave_rms_chain_part(600), & ave_assoc_part(600),ave_hbond_intra_part(600), & ave_hbond_inter_part(600),ave_assoc_hbond_part(600) c write (iout,*) me," indstart",indstart(me)," indend",indend(me) call daread_ccoords(indstart(me),indend(me)) #endif #ifdef MPI do i=1,scount(me) #else do i=1,ncon #endif do j=1,nres do k=1,3 c(k,j)=allcart(k,j,i) c(k,j+nres)=allcart(k,j+nres,i) enddo enddo do k=1,3 c(k,nres+1)=c(k,1) c(k,nres+nres)=c(k,nres) enddo nss=nss_all(i) do j=1,nss ihpb(j)=ihpb_all(j,i) jhpb(j)=jhpb_all(j,i) enddo call int_from_cart1(.false.) do k=1,6 ft(k)=1.0d0 enddo call etotal(energia(0),fT) #ifdef DEBUG call enerprint(energia,ft) #endif do k=1,max_ene enetb(k,i)=energia(k) enddo call rmsd_intra(i) call monomer_contact(i) call hbonds_intra_inter(i) #ifdef DEBUG write (iout,*) "jcon",i write (iout,*) (rms_intrachain(j,i),j=1,nchain) write (iout,*) (npept_cont_interchain(j,i),j=1,nchain) write (iout,*) (npept_cont_intrachain(j,i),j=1,nchain) write (iout,*) (iassoc(j,i),j=1,nchain) call flush(iout) #endif enddo #ifdef DEBUG do i=1,ncon write (iout,*) "jcon",i write (iout,*) (rms_intrachain(j,i),j=1,nchain) write (iout,*) (npept_cont_interchain(j,i),j=1,nchain) write (iout,*) (npept_cont_intrachain(j,i),j=1,nchain) write (iout,*) (iassoc(j,i),j=1,nchain) enddo #endif do ib=200,500 beta=1.0d0/(1.987d-3*ib) Fmin=1.0d10 #ifdef MPI do i=1,scount(me) ii=i+indstart(me)-1 #else do i=1,ncon ii=i #endif evdw=enetb(1,i) #ifdef SCP14 evdw2_14=enetb(17,i) evdw2=enetb(2,i)+evdw2_14 #else evdw2=enetb(2,i) evdw2_14=0.0d0 #endif #ifdef SPLITELE ees=enetb(3,i) evdw1=enetb(16,i) #else ees=enetb(3,i) evdw1=0.0d0 #endif ecorr=enetb(4,i) ecorr5=enetb(5,i) ecorr6=enetb(6,i) eel_loc=enetb(7,i) eello_turn3=enetb(8,i) eello_turn4=enetb(9,i) eturn6=enetb(10,i) ebe=enetb(11,i) escloc=enetb(12,i) etors=enetb(13,i) etors_d=enetb(14,i) ehpb=enetb(15,i) estr=enetb(18,i) esccor=enetb(19,i) edihcnstr=enetb(20,i) evdw_t=enetb(21,i) #ifdef DEBUG write (iout,*) "evdw", wsc, evdw,evdw_t write (iout,*) "evdw2", wscp, evdw2 write (iout,*) "welec", ft(1),welec,ees write (iout,*) "evdw1", wvdwpp,evdw1 write (iout,*) "ebe" ebe,wang #endif #ifdef DEBUG write (iout,*) "fdim calc", i,ii, & totfree(i),entfac(ii),Fdimless(i) #endif temper=ib+0.0d0 if (rescale_mode.eq.1) then quot=1.0d0/(T0*beta*1.987D-3) quotl=1.0d0 kfacl=1.0d0 do l=1,5 quotl1=quotl quotl=quotl*quot kfacl=kfacl*kfac fT(l)=kfacl/(kfacl-1.0d0+quotl) enddo #if defined(FUNCTH) ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ & 320.0d0 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2) ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3) #elif defined(FUNCT) fT(6)=betaT/T0 ftprim(6)=1.0d0/T0 ftbis(6)=0.0d0 #else fT(6)=1.0d0 ftprim(6)=0.0d0 ftbis(6)=0.0d0 #endif else if (rescale_mode.eq.2) then quot=1.0d0/(T0*beta*1.987D-3) quotl=1.0d0 do l=1,5 quotl=quotl*quot fT(l)=1.12692801104297249644d0/ & dlog(dexp(quotl)+dexp(-quotl)) enddo #if defined(FUNCTH) ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ & 320.0d0 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2) ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3) #elif defined(FUNCT) fT(6)=betaT/T0 ftprim(6)=1.0d0/T0 ftbis(6)=0.0d0 #else fT(6)=1.0d0 ftprim(6)=0.0d0 ftbis(6)=0.0d0 #endif endif c write (iout,*) "ft",(ft(k),k=1,6) #ifdef SPLITELE etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+welec*ft(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc & +wstrain*ehpb+nss*ebr+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d & +wbond*estr+wsccor*ft(1)*esccor #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +welec*ft(1)*(ees+evdw1) & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc & +wstrain*ehpb+nss*ebr+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d & +wbond*estr+wsccor*ft(1)*esccor #endif #ifdef DEBUG write (iout,10) (evdw+ft(6)*evdw_t),wsc,evdw2,wscp,ees, & welec*ft(1),evdw1, & wvdwpp, & estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*ft(1), & etors_d,wtor_d*ft(2),ehpb,wstrain, & ecorr,wcorr*ft(3),ecorr5,wcorr5*ft(4),ecorr6,wcorr6*ft(5), & eel_loc,wel_loc*ft(2),eello_turn3,wturn3*ft(2), & eello_turn4,wturn4*ft(3),eturn6,wturn6*ft(5), & esccor,wsccor*ft(1),edihcnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p elec)'/ & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/ & 'ESTR= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/ & 'EBE= ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/ & 'ESC= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/ & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/ & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/ & 'EHBP= ',1pE16.6,' WEIGHT=',1pD16.6, & ' (SS bridges & dist. cnstr.)'/ & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/ & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/ & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/ & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif totfree(i)=beta*etot+entfac(ii) #ifdef DEBUG write (iout,*) "ib",ib," beta",beta," i",i," ii",ii, & " etot",etot,entfac(ii),totfree(i) #endif if (totfree(i).lt.Fmin) Fmin=totfree(i) enddo #ifdef MPI call MPI_Allreduce(Fmin,Fmin_all,1, & MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,IERROR) #else Fmin_all=Fmin #endif qfree_part(ib)=0.0d0 ave_rms_chain_part(ib)=0.0d0 ave_assoc_part(ib)=0.0d0 ave_hbond_intra_part(ib)=0.0d0 ave_hbond_inter_part(ib)=0.0d0 ave_assoc_hbond_part(ib)=0.0d0 #ifdef MPI do i=1,scount(me) w=dexp(-totfree(i)+Fmin_all) c write (iout,*) "i",i," w",w qfree_part(ib)=qfree_part(ib)+w w=w/nchain do j=1,nchain ave_rms_chain_part(ib)=ave_rms_chain_part(ib) & +w*rms_intrachain(j,i) ave_assoc_part(ib)=ave_assoc_part(ib)+w*iassoc(j,i) ave_hbond_intra_part(ib)=ave_hbond_intra_part(ib) & +w*npept_cont_intrachain(j,i) ave_hbond_inter_part(ib)=ave_hbond_inter_part(ib) & +w*npept_cont_interchain(j,i) if (npept_cont_interchain(j,i).gt.0) & ave_assoc_hbond_part(ib)=ave_assoc_hbond_part(ib)+w enddo enddo enddo call MPI_Reduce(qfree_part(1),qfree(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) call MPI_Reduce(ave_rms_chain_part(1),ave_rms_chain(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) call MPI_Reduce(ave_assoc_part(1),ave_assoc(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) call MPI_Reduce(ave_hbond_intra_part(1),ave_hbond_intra(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) call MPI_Reduce(ave_hbond_inter_part(1),ave_hbond_inter(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) call MPI_Reduce(ave_assoc_hbond_part(1),ave_assoc_hbond(1),600, & MPI_DOUBLE_PRECISION,MPI_SUM,Master,MPI_COMM_WORLD,IERROR) #else qfree(ib)=0.0d0 ave_rms_chain(ib)=0.0d0 ave_assoc(ib)=0.0d0 ave_hbond_intra(ib)=0.0d0 ave_hbond_inter(ib)=0.0d0 ave_assoc_hbond(ib)=0.0d0 do i=1,ncon w=dexp(-totfree(i)+Fmin_all) qfree(ib)=qfree(ib)+w w=w/nchain do j=1,nchain ave_rms_chain(ib)=ave_rms_chain(ib)+w*rms_intrachain(j,i) ave_assoc(ib)=ave_assoc(ib)+w*iassoc(j,i) ave_hbond_intra(ib)=ave_hbond_intra(ib) & +w*npept_cont_intrachain(j,i) ave_hbond_inter(ib)=ave_hbond_inter(ib) & +w*npept_cont_interchain(j,i) if (npept_cont_interchain(j,i).gt.0) & ave_assoc_hbond(ib)=ave_assoc_hbond(ib)+w enddo enddo enddo #endif #ifdef MPI if (me.eq.Master) then #endif do ib=200,500 write (istat,'(i5,5f10.5)') ib,ave_rms_chain(ib)/qfree(ib), & ave_assoc(ib)/qfree(ib),ave_assoc_hbond(ib)/qfree(ib), & ave_hbond_intra(ib)/qfree(ib),ave_hbond_inter(ib)/qfree(ib) enddo call flush(istat) #ifdef MPI endif #endif return end !-------------------------------------------------- subroutine mysort1(n, x, ipermut) implicit none integer i,j,imax,ipm,n real x(n) integer ipermut(n) real xtemp 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 !-------------------------------------------------- subroutine rmsd_intra(jcon) implicit none include "DIMENSIONS" include "sizesclu.dat" include "COMMON.CONTROL" include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.CLUSTER" include "COMMON.PROPERTY" integer i,ii,j,k double precision xx(3,maxres),yy(3,maxres) double precision rms integer jcon logical non_conv double precision przes(3),obrot(3,3) if (lside) then do i=1,nchain ii=0 do k=ichain(1,i),ichain(2,i) ii=ii+1 do j=1,3 xx(j,ii)=allcart(j,i,jcon) yy(j,ii)=cref(j,i,1) enddo enddo do k=ichain(1,i),ichain(2,i) c if (itype(i).ne.10) then ii=ii+1 do j=1,3 xx(j,ii)=allcart(j,i+nres,jcon) yy(j,ii)=cref(j,i+nres,1) enddo c endif enddo call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv) rms_intrachain(i,jcon)=dsqrt(rms) enddo else do i=nnt,nct do j=1,3 xx(j,i)=allcart(j,i,jcon) enddo enddo #ifdef DEBUG write (iout,*) "Conformation",jcon #endif do i=1,nchain #ifdef DEBUG write (iout,*) "chain",i #endif do j=1,ichain(2,i)-ichain(1,i)+1 #ifdef DEBUG write (iout,'(i5,3f10.5,5x,3f10.5)') j, & (xx(k,ichain(1,i)+j-1),k=1,3), & (cref(k,ichain(1,i)+j-1,1),k=1,3) #endif enddo call fitsq(rms,xx(1,ichain(1,i)),cref(1,ichain(1,i),1), & ichain(2,i)-ichain(1,i)+1,przes,obrot,non_conv) rms_intrachain(i,jcon)=dsqrt(rms) #ifdef DEBUG write (iout,*) "rms",rms #endif enddo endif return end !-------------------------------------------------- subroutine hbonds_intra_inter(jcon) implicit none include "DIMENSIONS" include "sizesclu.dat" include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.CLUSTER" include "COMMON.PROPERTY" integer ncont,icont(2,maxcont) integer i,j,jcon,nr1,nr2 do i=1,nchain npept_cont_interchain(i,jcon)=0 npept_cont_intrachain(i,jcon)=0 enddo do i=1,nres do j=1,3 c(j,i)=allcart(j,i,jcon) enddo enddo call elecont(.false.,ncont,icont,nnt,nct) do i=1,ncont nr1=nres_chain(icont(1,i)) nr2=nres_chain(icont(2,i)) if (nr1.eq.nr2) then npept_cont_intrachain(nr1,jcon)= & npept_cont_intrachain(nr1,jcon)+1 else npept_cont_interchain(nr1,jcon)= & npept_cont_interchain(nr1,jcon)+1 npept_cont_interchain(nr2,jcon)= & npept_cont_interchain(nr2,jcon)+1 endif enddo return end !-------------------------------------------------- subroutine monomer_contact(jcon) implicit none include "DIMENSIONS" include "sizesclu.dat" include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.CLUSTER" include "COMMON.PROPERTY" integer ncont,icont(2,maxcont) integer i,j,ic1,ic2,jcon do i=1,nchain iassoc(i,jcon)=0 enddo do i=1,nres do j=1,3 c(j,i)=allcart(j,i,jcon) enddo enddo call contact(.false.,ncont,icont,nnt,nct) do i=1,ncont ic1=nres_chain(icont(1,i)) ic2=nres_chain(icont(2,i)) if (ic1.ne.ic2) then iassoc(ic1,jcon)=1 iassoc(ic2,jcon)=1 endif enddo return end