module conform_compar !----------------------------------------------------------------------------- use names use io_units use geometry_data, only:nres use math, only:pinorm use geometry, only:dist use regularize_, only:fitsq ! use wham_data #ifndef CLUSTER use w_compar_data #endif #ifdef MPI use MPI_data ! include "COMMON.MPI" #endif implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains #ifndef CLUSTER !----------------------------------------------------------------------------- ! conf_compar.F !----------------------------------------------------------------------------- subroutine conf_compar(jcon,lprn,print_class) ! implicit real*8 (a-h,o-z) use energy_data, only:icont,ncont,nnt,nct,maxcont!,& ! nsccont_frag_ref,isccont_frag_ref #ifdef MPI include "mpif.h" #endif ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'DIMENSIONS.FREE' ! include 'COMMON.CONTROL' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.PEPTCONT' ! include 'COMMON.CONTACTS1' ! include 'COMMON.HEADER' ! include 'COMMON.FREE' ! include 'COMMON.ENERGIES' !#ifdef MPI ! include 'COMMON.MPI' !#endif ! integer ilen ! external ilen logical :: lprn,print_class integer :: ncont_frag(mmaxfrag),& icont_frag(2,maxcont,mmaxfrag),ncontsc,& icontsc(1,maxcont),nsccont_frag(mmaxfrag),& isccont_frag(2,maxcont,mmaxfrag) integer :: isecstr(nres) integer :: itemp(maxfrag) character(len=4) :: liczba real(kind=8) :: Epot,rms integer :: jcon,i,j,ind,ncnat,nsec_match,ishift,ishif1,ishif2,& nc_match,ncon_match,iclass_rms,ishifft_rms,ishiff,ishif integer :: k,kk,iclass_con,iscor,ik,ishifft_con,idig,iex,im ! print *,"Enter conf_compar",jcon call angnorm12(rmsang) ! Level 1: check secondary and supersecondary structure call elecont(lprn,ncont,icont,nnt,nct) if (lprn) then write (iout,*) "elecont finished" call flush(iout) endif call secondary2(lprn,.false.,ncont,icont,isecstr) if (lprn) then write (iout,*) "secondary2 finished" call flush(iout) endif call contact(lprn,ncontsc,icontsc,nnt,nct) if (lprn) then write(iout,*) "Assigning electrostatic contacts" call flush(iout) endif call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,& icont_frag) if (lprn) then write(iout,*) "Assigning sidechain contacts" call flush(iout) endif call contacts_between_fragments(lprn,3,ncontsc,icontsc,& nsccont_frag,isccont_frag) if (lprn) then write(iout,*) "--> After contacts_between_fragments" call flush(iout) endif do i=1,nlevel do j=1,isnfrag(nlevel+1) iclass(j,i)=0 enddo enddo do j=1,nfrag(1) ind = icant(j,j) if (lprn) then write (iout,'(80(1h=))') write (iout,*) "Level",1," fragment",j write (iout,'(80(1h=))') endif call flush(iout) rmsfrag(j,1)=rmscalc(0,1,j,jcon,lprn) ! Compare electrostatic contacts in the current conf with that in the native ! structure. if (lprn) write (iout,*) & "Comparing electrostatic contact map and local structure" call flush(iout) ncnat=ncont_frag_ref(ind) ! write (iout,*) "before match_contact:",nc_fragm(j,1), ! & nc_req_setf(j,1) ! call flush(iout) call match_secondary(j,isecstr,nsec_match,lprn) if (lprn) write (iout,*) "Fragment",j," nsec_match",& nsec_match," length",len_frag(j,1)," min_len",& frac_sec*len_frag(j,1) if (nsec_match.lt.frac_sec*len_frag(j,1)) then iclass(j,1)=0 if (lprn) write (iout,*) "Fragment",j,& " has incorrect secondary structure" else iclass(j,1)=1 if (lprn) write (iout,*) "Fragment",j,& " has correct secondary structure" endif if (ielecont(j,1).gt.0) then call match_contact(ishif1,ishif2,nc_match,ncon_match,& ncont_frag_ref(ind),icont_frag_ref(1,1,ind),& ncont_frag(ind),icont_frag(1,1,ind),& j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),& nc_req_setf(j,1),istruct(j),.true.,lprn) else if (isccont(j,1).gt.0) then call match_contact(ishif1,ishif2,nc_match,ncon_match,& nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),& nsccont_frag(ind),isccont_frag(1,1,ind),& j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),& nc_req_setf(j,1),istruct(j),.true.,lprn) else if (iloc(j).gt.0) then ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1) call match_contact(ishif1,ishif2,nc_match,ncon_match,& 0,icont_frag_ref(1,1,ind),& ncont_frag(ind),icont_frag(1,1,ind),& j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),& 0,istruct(j),.true.,lprn) ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1) else ishif=0 nc_match=1 endif if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2 ishif=ishif1 qfrag(j,1)=qwolynes(1,j) if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2 if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match ! write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1) if (irms(j,1).gt.0) then if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then iclass_rms=2 ishifft_rms=0 else ishiff=0 rms=1.0d2 iclass_rms=0 do while (rms.gt.rmscutfrag(1,j,1) .and. & ishiff.lt.n_shift(1,j,1)) ishiff=ishiff+1 rms=rmscalc(-ishiff,1,j,jcon,lprn) ! write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff, ! & " rms",rms," rmscut",rmscutfrag(1,j,1) if (lprn) write (iout,*) "rms",rmsfrag(j,1) if (rms.gt.rmscutfrag(1,j,1)) then rms=rmscalc(ishiff,1,j,jcon,lprn) ! write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff, ! & " rms",rms endif if (lprn) write (iout,*) "rms",rmsfrag(j,1) enddo ! write (iout,*) "After loop: rms",rms, ! & " rmscut",rmscutfrag(1,j,1) ! write (iout,*) "iclass_rms",iclass_rms if (rms.le.rmscutfrag(1,j,1)) then ishifft_rms=ishiff rmsfrag(j,1)=rms iclass_rms=1 endif ! write (iout,*) "iclass_rms",iclass_rms endif ! write (iout,*) "ishif",ishif if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms else iclass_rms=1 endif ! write (iout,*) "ishif",ishif," iclass",iclass(j,1), ! & " iclass_rms",iclass_rms if (nc_match.gt.0 .and. iclass_rms.gt.0) then if (ishif.eq.0) then iclass(j,1)=iclass(j,1)+6 else iclass(j,1)=iclass(j,1)+2 endif endif ncont_nat(1,j,1)=nc_match ncont_nat(2,j,1)=ncon_match ishifft(j,1)=ishif ! write (iout,*) "iclass",iclass(j,1) enddo ! Next levels: Check arrangements of elementary fragments. do i=2,nlevel do j=1,nfrag(i) if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i)) if (lprn) then write (iout,'(80(1h=))') write (iout,*) "Level",i," fragment",j write (iout,'(80(1h=))') endif ! If an elementary fragment doesn't exist, don't check higher hierarchy levels. do k=1,npiece(j,i) ik=ipiece(k,j,i) if (iclass(ik,1).eq.0) then iclass(j,i)=0 goto 12 endif enddo if (i.eq.2 .and. ielecont(j,i).gt.0) then iclass_con=0 ishifft_con=0 if (lprn) write (iout,*) & "Comparing electrostatic contact map: fragments",& ipiece(1,j,i),ipiece(2,j,i)," ind",ind call match_contact(ishif1,ishif2,nc_match,ncon_match,& ncont_frag_ref(ind),icont_frag_ref(1,1,ind),& ncont_frag(ind),icont_frag(1,1,ind),& j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),& nc_req_setf(j,i),2,.false.,lprn) ishif=ishif1 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2 if (nc_match.gt.0) then if (ishif.eq.0) then iclass_con=2 else iclass_con=1 endif endif ncont_nat(1,j,i)=nc_match ncont_nat(2,j,i)=ncon_match ishifft_con=ishif else if (i.eq.2 .and. isccont(j,i).gt.0) then iclass_con=0 ishifft_con=0 if (lprn) write (iout,*) & "Comparing sidechain contact map: fragments",& ipiece(1,j,i),ipiece(2,j,i)," ind",ind call match_contact(ishif1,ishif2,nc_match,ncon_match,& nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),& nsccont_frag(ind),isccont_frag(1,1,ind),& j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),& nc_req_setf(j,i),2,.false.,lprn) ishif=ishif1 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2 if (nc_match.gt.0) then if (ishif.eq.0) then iclass_con=2 else iclass_con=1 endif endif ncont_nat(1,j,i)=nc_match ncont_nat(2,j,i)=ncon_match ishifft_con=ishif else if (i.eq.2) then iclass_con=2 ishifft_con=0 endif if (i.eq.2) qfrag(j,2)=qwolynes(2,j) if (lprn) write (iout,*) & "Comparing rms: fragments",& (ipiece(k,j,i),k=1,npiece(j,i)) rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn) if (irms(j,i).gt.0) then iclass_rms=0 ishifft_rms=0 if (lprn) write (iout,*) "rms",rmsfrag(j,i) ! write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i), ! & " rmscutfrag",rmscutfrag(1,j,i) if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then iclass_rms=2 ishifft_rms=0 else ishif=0 rms=1.0d2 do while (rms.gt.rmscutfrag(1,j,i) .and. & ishif.lt.n_shift(1,j,i)) ishif=ishif+1 rms=rmscalc(-ishif,i,j,jcon,lprn) ! print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms if (lprn) write (iout,*) "rms",rmsfrag(j,i) if (rms.gt.rmscutfrag(1,j,i)) then rms=rmscalc(ishif,i,j,jcon,lprn) ! print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms endif if (lprn) write (iout,*) "rms",rms enddo if (rms.le.rmscutfrag(1,j,i)) then ishifft_rms=ishif rmsfrag(j,i)=rms iclass_rms=1 endif endif endif if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and. & isccont(j,i).eq.0 ) then write (iout,*) "Error: no measure of comparison specified:",& " level",i," part",j stop endif if (lprn) & write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms if (i.eq.2) then iclass(j,i) = min0(iclass_con,iclass_rms) if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then ishifft(j,i)=ishifft_rms else ishifft(j,i)=ishifft_con endif else if (i.gt.2) then iclass(j,i) = iclass_rms ishifft(j,i)= ishifft_rms endif 12 continue enddo enddo rms_nat=rmsnat(jcon) qnat=qwolynes(0,0) ! Compute the structural class iscor=0 IF (.NOT. BINARY) THEN do i=1,nlevel IF (I.EQ.1) THEN do j=1,nfrag(i) itemp(j)=iclass(j,i) enddo do kk=-1,1 do j=1,nfrag(i) idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j iex = 2**idig im=mod(itemp(j),2) itemp(j)=itemp(j)/2 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex, ! & " iclass",iclass(j,i)," im",im iscor=iscor+im*iex enddo enddo ELSE do j=1,nfrag(i) idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j iex = 2**idig if (iclass(j,i).gt.0) then im=1 else im=0 endif ! write (iout,*) "i",i," j",j," idig",idig," iex",iex, ! & " iclass",iclass(j,i)," im",im iscor=iscor+im*iex enddo do j=1,nfrag(i) idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j iex = 2**idig if (iclass(j,i).gt.1) then im=1 else im=0 endif ! write (iout,*) "i",i," j",j," idig",idig," iex",iex, ! & " iclass",iclass(j,i)," im",im iscor=iscor+im*iex enddo ENDIF enddo iscore=iscor ENDIF if (print_class) then #ifdef MPI write(istat,'(i6,$)') jcon+indstart(me)-1 write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),& -entfac(jcon) #else write(istat,'(i6,$)') jcon write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),& -entfac(jcon) #endif write (istat,'(f8.3,2f6.3,$)') & rms_nat,qnat,rmsang/(nres-3) do j=1,nlevel write(istat,'(1x,$,20(i3,$))') & (ncont_nat(1,k,j),k=1,nfrag(j)) if (j.lt.3) then write(istat,'(1x,$,20(f5.1,f5.2$))') & (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j)) else write(istat,'(1x,$,20(f5.1$))') & (rmsfrag(k,j),k=1,nfrag(j)) endif write(istat,'(1x,$,20(i1,$))') & (iclass(k,j),k=1,nfrag(j)) enddo if (binary) then write (istat,'(" ",$)') do j=1,nlevel write (istat,'(100(i1,$))')(iclass(k,j),& k=1,nfrag(j)) if (j.lt.nlevel) write(iout,'(".",$)') enddo write (istat,*) else write (istat,'(i10)') iscore endif endif RETURN END subroutine conf_compar !----------------------------------------------------------------------------- ! angnorm.f !----------------------------------------------------------------------------- subroutine add_angpair(ici,icj,nang_pair,iang_pair) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' integer :: ici,icj,nang_pair,iang_pair(2,nres) integer :: i,ian1,ian2 ! write (iout,*) "add_angpair: ici",ici," icj",icj, ! & " nang_pair",nang_pair ian1=ici+2 if (ian1.lt.4 .or. ian1.gt.nres) return ian2=icj+2 ! write (iout,*) "ian1",ian1," ian2",ian2 if (ian2.lt.4 .or. ian2.gt.nres) return do i=1,nang_pair if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return enddo nang_pair=nang_pair+1 iang_pair(1,nang_pair)=ian1 iang_pair(2,nang_pair)=ian2 return end subroutine add_angpair !------------------------------------------------------------------------- subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,lprn) use geometry_data, only:nstart_sup,nend_sup,phi,theta,& rad2deg,dwapi ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' real(kind=8) :: pinorm,deltang logical :: lprn integer :: jfrag,ishif1,ishif2,nn,npart,nn4,nne real(kind=8) :: diffang_max,angn,fract,ff integer :: i,j,nbeg,nend,ll,longest if (lprn) write (iout,'(80(1h*))') angn=0.0d0 nn = 0 fract = 1.0d0 npart = npiece(jfrag,1) nn4 = nstart_sup+3 nne = min0(nend_sup,nres) if (lprn) write (iout,*) "nn4",nn4," nne",nne do i=1,npart nbeg = ifrag(1,i,jfrag) + 3 - ishif1 if (nbeg.lt.nn4) nbeg=nn4 nend = ifrag(2,i,jfrag) + 1 - ishif2 if (nend.gt.nne) nend=nne if (nend.ge.nbeg) then nn = nn + nend - nbeg + 1 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,& " nn",nn," ishift1",ishif1," ishift2",ishif2 if (lprn) write (iout,*) "angles" longest=0 ll = 0 do j=nbeg,nend ! deltang = pinorm(phi(j)-phi_ref(j+ishif1)) deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),& theta_ref(j+ishif1),phi(j),theta(j-1),theta(j)) if (dabs(deltang).gt.diffang_max) then if (ll.gt.longest) longest = ll ll = 0 else ll=ll+1 endif if (ll.gt.longest) longest = ll if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),& rad2deg*phi_ref(j+ishif1),rad2deg*deltang angn=angn+dabs(deltang) enddo longest=longest+3 ff = dfloat(longest)/dfloat(nend - nbeg + 4) if (lprn) write (iout,*)"segment",i," longest fragment within",& diffang_max*rad2deg,":",longest," fraction",ff if (ff.lt.fract) fract = ff endif enddo if (nn.gt.0) then angn = angn/nn else angn = dwapi endif if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,& " fract",fract return end subroutine angnorm !------------------------------------------------------------------------- subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,& diffang_max,anorm,fract) use geometry_data, only:nstart_sup,nend_sup,phi,theta,& rad2deg ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' integer :: ncont,icont(2,ncont),longest real(kind=8) :: anorm,diffang_max,fract integer :: npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece) real(kind=8) :: pinorm logical :: lprn integer :: jfrag,ishif1,ishif2 integer :: nn,nn4,nne,npart,i,j,jstart,jend,ic1,ic2,idi,iic integer :: nbeg,nend,ll real(kind=8) :: angn,ishifc,deltang,ff if (lprn) write (iout,'(80(1h*))') ! ! Determine the segments for which angles will be compared ! nn4 = nstart_sup+3 nne = min0(nend_sup,nres) if (lprn) write (iout,*) "nn4",nn4," nne",nne npart=npiece(jfrag,1) npiece_c=0 do i=1,npart ! write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. & icont(1,1).gt.ifrag(2,i,jfrag)) goto 11 jstart=1 do while (jstart.lt.ncont .and. & icont(1,jstart).lt.ifrag(1,i,jfrag)) ! write (iout,*) "jstart",jstart," icont",icont(1,jstart), ! & " ifrag",ifrag(1,i,jfrag) jstart=jstart+1 enddo ! write (iout,*) "jstart",jstart," icont",icont(1,jstart), ! & " ifrag",ifrag(1,i,jfrag) if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11 npiece_c=npiece_c+1 ic1=icont(1,jstart) ifrag_c(1,npiece_c)=icont(1,jstart) jend=ncont do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag)) ! write (iout,*) "jend",jend," icont",icont(1,jend), ! & " ifrag",ifrag(2,i,jfrag) jend=jend-1 enddo ! write (iout,*) "jend",jend," icont",icont(1,jend), ! & " ifrag",ifrag(2,i,jfrag) ic2=icont(1,jend) ifrag_c(2,npiece_c)=icont(1,jend)+1 ishift_c(npiece_c)=ishif1 ! write (iout,*) "1: i",i," jstart:",jstart," jend",jend, ! & " ic1",ic1," ic2",ic2, ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) 11 continue if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then idi=1 else idi=-1 endif ! write (iout,*) "idi",idi if (idi.eq.1) then if (icont(2,1).gt.ifrag(2,i,jfrag) .or. & icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12 jstart=1 do while (jstart.lt.ncont .and. & icont(2,jstart).lt.ifrag(1,i,jfrag)) ! write (iout,*) "jstart",jstart," icont",icont(2,jstart), ! & " ifrag",ifrag(1,i,jfrag) jstart=jstart+1 enddo ! write (iout,*) "jstart",jstart," icont",icont(2,jstart), ! & " ifrag",ifrag(1,i,jfrag) if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12 npiece_c=npiece_c+1 ic1=icont(2,jstart) ifrag_c(2,npiece_c)=icont(2,jstart)+1 jend=ncont do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag)) ! write (iout,*) "jend",jend," icont",icont(2,jend), ! & " ifrag",ifrag(2,i,jfrag) jend=jend-1 enddo ! write (iout,*) "jend",jend," icont",icont(2,jend), ! & " ifrag",ifrag(2,i,jfrag) else if (idi.eq.-1) then if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or. & icont(2,1).lt.ifrag(1,i,jfrag)) goto 12 jstart=ncont do while (jstart.gt.ncont .and. & icont(2,jstart).lt.ifrag(1,i,jfrag)) ! write (iout,*) "jstart",jstart," icont",icont(2,jstart), ! & " ifrag",ifrag(1,i,jfrag) jstart=jstart-1 enddo ! write (iout,*) "jstart",jstart," icont",icont(2,jstart), ! & " ifrag",ifrag(1,i,jfrag) if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12 npiece_c=npiece_c+1 ic1=icont(2,jstart) ifrag_c(2,npiece_c)=icont(2,jstart)+1 jend=1 do while (jend.lt.ncont .and. & icont(2,jend).gt.ifrag(2,i,jfrag)) ! write (iout,*) "jend",jend," icont",icont(2,jend), ! & " ifrag",ifrag(2,i,jfrag) jend=jend+1 enddo ! write (iout,*) "jend",jend," icont",icont(2,jend), ! & " ifrag",ifrag(2,i,jfrag) endif ic2=icont(2,jend) if (ic2.lt.ic1) then iic = ic1 ic1 = ic2 ic2 = iic endif ! write (iout,*) "2: i",i," ic1",ic1," ic2",ic2, ! & " jstart:",jstart," jend",jend, ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag) ifrag_c(1,npiece_c)=ic1 ifrag_c(2,npiece_c)=ic2+1 ishift_c(npiece_c)=ishif2 12 continue enddo if (lprn) then write (iout,*) "Before merge: npiece_c",npiece_c do i=1,npiece_c write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i) enddo endif ! ! Merge overlapping segments (e.g., avoid splitting helices) ! i=1 do while (i .lt. npiece_c) if (ishift_c(i).eq.ishift_c(i+1) .and. & ifrag_c(2,i).gt.ifrag_c(1,i+1)) then ifrag_c(2,i)=ifrag_c(2,i+1) do j=i+1,npiece_c ishift_c(j)=ishift_c(j+1) ifrag_c(1,j)=ifrag_c(1,j+1) ifrag_c(2,j)=ifrag_c(2,j+1) enddo npiece_c=npiece_c-1 else i=i+1 endif enddo if (lprn) then write (iout,*) "After merge: npiece_c",npiece_c do i=1,npiece_c write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i) enddo endif ! ! Compare angles ! angn=0.0d0 anorm=0 nn = 0 fract = 1.0d0 npart = npiece_c do i=1,npart ishifc=ishift_c(i) nbeg = ifrag_c(1,i) + 3 - ishifc if (nbeg.lt.nn4) nbeg=nn4 nend = ifrag_c(2,i) - ishifc + 1 if (nend.gt.nne) nend=nne if (nend.ge.nbeg) then nn = nn + nend - nbeg + 1 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,& " nn",nn," ishifc",ishifc if (lprn) write (iout,*) "angles" longest=0 ll = 0 do j=nbeg,nend ! deltang = pinorm(phi(j)-phi_ref(j+ishifc)) deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),& theta_ref(j+ishifc),phi(j),theta(j-1),theta(j)) if (dabs(deltang).gt.diffang_max) then if (ll.gt.longest) longest = ll ll = 0 else ll=ll+1 endif if (ll.gt.longest) longest = ll if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),& rad2deg*phi_ref(j+ishifc),rad2deg*deltang angn=angn+dabs(deltang) enddo longest=longest+3 ff = dfloat(longest)/dfloat(nend - nbeg + 4) if (lprn) write (iout,*)"segment",i," longest fragment within",& diffang_max*rad2deg,":",longest," fraction",ff if (ff.lt.fract) fract = ff endif enddo if (nn.gt.0) anorm = angn/nn if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract return end subroutine angnorm2 !------------------------------------------------------------------------- real(kind=8) function angnorm1(nang_pair,iang_pair,lprn) use geometry_data, only:phi,theta,rad2deg ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' logical :: lprn integer :: nang_pair,iang_pair(2,nres) real(kind=8) :: pinorm integer :: j,ia1,ia2 real(kind=8) :: angn,deltang angn=0.0d0 if (lprn) write (iout,'(80(1h*))') if (lprn) write (iout,*) "nang_pair",nang_pair if (lprn) write (iout,*) "angles" do j=1,nang_pair ia1 = iang_pair(1,j) ia2 = iang_pair(2,j) ! deltang = pinorm(phi(ia1)-phi_ref(ia2)) deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),& theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2)) if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),& rad2deg*phi_ref(ia2),rad2deg*deltang angn=angn+dabs(deltang) enddo if (lprn) & write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair angnorm1 = angn/nang_pair return end function angnorm1 !------------------------------------------------------------------------------ subroutine angnorm12(diff) use geometry_data, only:phi,theta,nstart_sup,nend_sup ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' real(kind=8) :: pinorm,diff integer :: nn4,nne,j diff=0.0d0 nn4 = nstart_sup+3 nne = min0(nend_sup,nres) ! do j=nn4-1,nne ! diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j))) ! enddo do j=nn4,nne ! diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j))) diff=diff+spherang(phi_ref(j),theta_ref(j-1),& theta_ref(j),phi(j),theta(j-1),theta(j)) enddo return end subroutine angnorm12 !-------------------------------------------------------------------------------- real(kind=8) function spherang(gam1,theta11,theta12,& gam2,theta21,theta22) ! implicit none use geometry, only:arcos real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,& x1,x2,xmed,f1,f2,fmed real(kind=8) :: tolx=1.0d-4, tolf=1.0d-4 real(kind=8) :: sumcos !el real(kind=8) :: pinorm,sumangp !arcos, integer :: it,maxit=100 ! Calculate the difference of the angles of two superposed 4-redidue fragments ! ! O P ! \ / ! O'--C--C ! \ ! P' ! ! The fragment O'-C-C-P' is rotated by angle fi about the C-C axis ! to achieve the minimum difference between the O'-C-O and P-C-P angles; ! the sum of these angles is the difference returned by the function. ! ! 4/28/04 AL ! If thetas match, take the difference of gamma and exit. if (dabs(theta11-theta12).lt.tolx & .and. dabs(theta21-theta22).lt.tolx) then spherang=dabs(pinorm(gam2-gam1)) return endif ! If the gammas are the same, take the difference of thetas and exit. x1=0.0d0 x2=0.5d0*pinorm(gam2-gam1) if (dabs(x2) .lt. tolx) then spherang=dabs(theta11-theta21)+dabs(theta12-theta22) return else if (x2.lt.0.0d0) then x1=x2 x2=0.0d0 endif ! Else apply regula falsi method to compute optimum overlap of the terminal Calphas f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1) f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2) do it=1,maxit xmed=x1-f1*(x2-x1)/(f2-f1) fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed) ! write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx) & .and. dabs(fmed).lt.tolf ) then x1=xmed f1=fmed goto 10 else if ( fmed*f1.lt.0.0d0 ) then x2=xmed f2=fmed else x1=xmed f1=fmed endif enddo 10 continue spherang=arcos(dcos(theta11)*dcos(theta12) & +dsin(theta11)*dsin(theta12)*dcos(x1))+ & arcos(dcos(theta21)*dcos(theta22)+ & dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1)) return end function spherang !-------------------------------------------------------------------------------- real(kind=8) function sumangp(gam1,theta11,theta12,gam2,& theta21,theta22,fi) ! implicit none real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,fi,& cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,& cosd2 ! derivarive of the sum of the difference of the angles of a 4-residue fragment. ! real(kind=8) :: arcos cost11=dcos(theta11) cost12=dcos(theta12) cost21=dcos(theta21) cost22=dcos(theta22) sint11=dsin(theta11) sint12=dsin(theta12) sint21=dsin(theta21) sint22=dsin(theta22) cosd1=cost11*cost12+sint11*sint12*dcos(fi) cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi) sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1) & +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2) return end function sumangp !----------------------------------------------------------------------------- ! contact.f !----------------------------------------------------------------------------- subroutine contact(lprint,ncont,icont,ist,ien) use calc_data use geometry_data, only:c,dc,dc_norm use energy_data, only:itype,maxcont ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'COMMON.CONTROL' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' ! include 'COMMON.CALC' ! include 'COMMON.CONTPAR' ! include 'COMMON.LOCAL' integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2 real(kind=8) :: csc !el,dist real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,& ddsc,ddla,ddlb integer :: ncont integer,dimension(2,maxcont) :: icont real(kind=8) :: u,v,a(3),b(3),dla,dlb logical :: lprint !el------- dla=0.0d0 dlb=0.0d0 !el------ ncont=0 kkk=3 if (lprint) then do i=1,nres write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),& c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),& dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i) enddo endif 110 format (a,'(',i3,')',9f8.3) do i=ist,ien-kkk iti=iabs(itype(i)) if (iti.le.0 .or. iti.gt.ntyp) cycle do j=i+kkk,ien itj=iabs(itype(j)) if (itj.le.0 .or. itj.gt.ntyp) cycle itypi=iti itypj=itj xj = c(1,nres+j)-c(1,nres+i) yj = c(2,nres+j)-c(2,nres+i) zj = c(3,nres+j)-c(3,nres+i) dxi = dc_norm(1,nres+i) dyi = dc_norm(2,nres+i) dzi = dc_norm(3,nres+i) dxj = dc_norm(1,nres+j) dyj = dc_norm(2,nres+j) dzj = dc_norm(3,nres+j) do k=1,3 a(k)=dc(k,nres+i) b(k)=dc(k,nres+j) enddo ! write (iout,*) (a(k),k=1,3),(b(k),k=1,3) if (icomparfunc.eq.1) then call contfunc(csc,iti,itj) else if (icomparfunc.eq.2) then call scdist(csc,iti,itj) else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then csc = dist(nres+i,nres+j) else if (icomparfunc.eq.4) then call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc) else write (*,*) "Error - Unknown sidechain contact function" write (iout,*) "Error - Unknown sidechain contact function" endif if (csc.lt.sc_cutoff(iti,itj)) then ! write(iout,*) "i",i," j",j," dla",dla,dsc(iti), ! & " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj), ! & dxi,dyi,dzi,dxi**2+dyi**2+dzi**2, ! & dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12, ! & xj,yj,zj ! write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2, ! & sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12, ! & chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw, ! & csc ncont=ncont+1 cscore(ncont)=csc icont(1,ncont)=i icont(2,ncont)=j omt1(ncont)=om1 omt2(ncont)=om2 omt12(ncont)=om12 ddsc(ncont)=1.0d0/rij ddla(ncont)=dla ddlb(ncont)=dlb endif enddo enddo if (lprint) then write (iout,'(a)') 'Contact map:' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') & i,restyp(it1),i1,restyp(it2),i2,cscore(i),& sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),& omt1(i),omt2(i),omt12(i) enddo endif return end subroutine contact #else !---------------------------------------------------------------------------- subroutine contact(lprint,ncont,icont) use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6) integer :: ncont,icont(2,maxcont) logical :: lprint integer :: kkk,i,j,i1,i2,it1,it2,iti,itj real(kind=8) :: rcomp ncont=0 kkk=3 ! print *,'nnt=',nnt,' nct=',nct do i=nnt+kkk,nct iti=iabs(itype(i)) do j=nnt,i-kkk itj=iabs(itype(j)) if (ipot.ne.4) then ! rcomp=sigmaii(iti,itj)+1.0D0 rcomp=facont*sigmaii(iti,itj) else ! rcomp=sigma(iti,itj)+1.0D0 rcomp=facont*sigma(iti,itj) endif ! rcomp=6.5D0 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j) if (dist(nres+i,nres+j).lt.rcomp) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j endif enddo enddo if (lprint) then write (iout,'(a)') 'Contact map:' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo endif return end subroutine contact #endif !---------------------------------------------------------------------------- real(kind=8) function contact_fract(ncont,ncont_ref,& icont,icont_ref) use energy_data, only:maxcont ! implicit none ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: i,j,nmatch integer :: ncont,ncont_ref integer,dimension(2,maxcont) :: icont,icont_ref nmatch=0 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont(1,i),i=1,ncont) ! write (iout,'(20i4)') (icont(2,i),i=1,ncont) do i=1,ncont do j=1,ncont_ref if (icont(1,i).eq.icont_ref(1,j) .and. & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1 enddo enddo ! print *,' nmatch=',nmatch ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref)) contact_fract=dfloat(nmatch)/dfloat(ncont_ref) return end function contact_fract #ifndef CLUSTER !------------------------------------------------------------------------------ subroutine pept_cont(lprint,ncont,icont) use geometry_data, only:c use energy_data, only:maxcont,nnt,nct,itype ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' integer :: ncont,icont(2,maxcont) integer :: i,j,k,kkk,i1,i2,it1,it2 logical :: lprint !el real(kind=8) :: dist real(kind=8) :: rcomp=5.5d0 ncont=0 kkk=0 print *,'Entering pept_cont: nnt=',nnt,' nct=',nct do i=nnt,nct-3 do k=1,3 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1)) enddo do j=i+2,nct-1 do k=1,3 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1)) enddo if (dist(2*nres+1,2*nres+2).lt.rcomp) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j endif enddo enddo if (lprint) then write (iout,'(a)') 'PP contact map:' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo endif return end subroutine pept_cont !----------------------------------------------------------------------------- ! cont_frag.f !----------------------------------------------------------------------------- subroutine contacts_between_fragments(lprint,is,ncont,icont,& ncont_interfrag,icont_interfrag) use energy_data, only:itype,maxcont ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.INTERACT' ! include 'COMMON.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.NAMES' integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),& icont_interfrag(2,maxcont,mmaxfrag) logical :: OK1,OK2,lprint integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2 ! Determine the contacts that occur within a fragment and between fragments. do i=1,nfrag(1) do j=1,i ind = icant(i,j) nc=0 ! write (iout,*) "i",i,(ifrag(1,k,i),ifrag(2,k,i) ! & ,k=1,npiece(i,1)) ! write (iout,*) "j",j,(ifrag(1,k,j),ifrag(2,k,j) ! & ,k=1,npiece(j,1)) ! write (iout,*) "ncont",ncont do k=1,ncont ic1=icont(1,k) ic2=icont(2,k) OK1=.false. l=0 do while (.not.OK1 .and. l.lt.npiece(j,1)) l=l+1 OK1=ic1.ge.ifrag(1,l,j)-is .and. & ic1.le.ifrag(2,l,j)+is enddo OK2=.false. l=0 do while (.not.OK2 .and. l.lt.npiece(i,1)) l=l+1 OK2=ic2.ge.ifrag(1,l,i)-is .and. & ic2.le.ifrag(2,l,i)+is enddo ! write(iout,*) "k",k," ic1",ic1," ic2",ic2," OK1",OK1, ! & " OK2",OK2 if (OK1.and.OK2) then nc=nc+1 icont_interfrag(1,nc,ind)=ic1 icont_interfrag(2,nc,ind)=ic2 ! write (iout,*) "nc",nc," ic1",ic1," ic2",ic2 endif enddo ncont_interfrag(ind)=nc ! do k=1,ncont_interfrag(ind) ! i1=icont_interfrag(1,k,ind) ! i2=icont_interfrag(2,k,ind) ! it1=itype(i1) ! it2=itype(i2) ! write (iout,'(i3,2x,a,i4,2x,a,i4)') ! & i,restyp(it1),i1,restyp(it2),i2 ! enddo enddo enddo if (lprint) then write (iout,*) "Contacts within fragments:" do i=1,nfrag(1) write (iout,*) "Fragment",i," (",(ifrag(1,k,i),& ifrag(2,k,i),k=1,npiece(i,1)),")" ind=icant(i,i) do k=1,ncont_interfrag(ind) i1=icont_interfrag(1,k,ind) i2=icont_interfrag(2,k,ind) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo enddo write (iout,*) write (iout,*) "Contacts between fragments:" do i=1,nfrag(1) do j=1,i-1 ind = icant(i,j) write (iout,*) "Fragments",i," (",(ifrag(1,k,i),& ifrag(2,k,i),k=1,npiece(i,1)),") and",j," (",& (ifrag(1,k,j),ifrag(2,k,j),k=1,npiece(j,1)),")" write (iout,*) "Number of contacts",& ncont_interfrag(ind) ind=icant(i,j) do k=1,ncont_interfrag(ind) i1=icont_interfrag(1,k,ind) i2=icont_interfrag(2,k,ind) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo enddo enddo endif return end subroutine contacts_between_fragments !----------------------------------------------------------------------------- ! contfunc.f !----------------------------------------------------------------------------- subroutine contfunc(cscore,itypi,itypj) ! ! This subroutine calculates the contact function based on ! the Gay-Berne potential of interaction. ! use calc_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTPAR' ! include 'COMMON.CALC' integer :: expon=6 integer :: itypi,itypj real(kind=8) :: cscore,sig0ij,rrij,sig,rij_shift,evdw,e2 ! sig0ij=sig_comp(itypi,itypj) chi1=chi_comp(itypi,itypj) chi2=chi_comp(itypj,itypi) chi12=chi1*chi2 chip1=chip_comp(itypi,itypj) chip2=chip_comp(itypj,itypi) chip12=chip1*chip2 rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) ! Calculate angle-dependent terms of the contact function erij(1)=xj*rij erij(2)=yj*rij erij(3)=zj*rij om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) om12=dxi*dxj+dyi*dyj+dzi*dzj chiom12=chi12*om12 ! print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2, ! & sig0ij, ! & rij,rrij,om1,om2,om12 ! Calculate eps1(om12) faceps1=1.0D0-om12*chiom12 faceps1_inv=1.0D0/faceps1 eps1=dsqrt(faceps1_inv) ! Following variable is eps1*deps1/dom12 eps1_om12=faceps1_inv*chiom12 ! Calculate sigma(om1,om2,om12) om1om2=om1*om2 chiom1=chi1*om1 chiom2=chi2*om2 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12 sigsq=1.0D0-facsig*faceps1_inv ! Calculate eps2 and its derivatives in om1, om2, and om12. chipom1=chip1*om1 chipom2=chip2*om2 chipom12=chip12*om12 facp=1.0D0-om12*chipom12 facp_inv=1.0D0/facp facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12 ! Following variable is the square root of eps2 eps2rt=1.0D0-facp1*facp_inv sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij if (rij_shift.le.0.0D0) then evdw=1.0D1 cscore = -dlog(evdw+1.0d-6) return endif rij_shift=1.0D0/rij_shift e2=(rij_shift*sig0ij)**expon evdw=dabs(eps1*eps2rt**2*e2) if (evdw.gt.1.0d1) evdw = 1.0d1 cscore = -dlog(evdw+1.0d-6) return end subroutine contfunc !------------------------------------------------------------------------------ subroutine scdist(cscore,itypi,itypj) ! ! This subroutine calculates the contact distance ! use calc_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTPAR' ! include 'COMMON.CALC' integer :: itypi,itypj real(kind=8) :: cscore,rrij chi1=chi_comp(itypi,itypj) chi2=chi_comp(itypj,itypi) chi12=chi1*chi2 rrij=xj*xj+yj*yj+zj*zj rij=dsqrt(rrij) ! Calculate angle-dependent terms of the contact function erij(1)=xj/rij erij(2)=yj/rij erij(3)=zj/rij om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) om12=dxi*dxj+dyi*dyj+dzi*dzj chiom12=chi12*om12 om1om2=om1*om2 chiom1=chi1*om1 chiom2=chi2*om2 cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12) return end subroutine scdist !------------------------------------------------------------------------------ ! elecont.f !------------------------------------------------------------------------------ subroutine elecont(lprint,ncont,icont,ist,ien) use geometry_data, only:c use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2 ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' ! include 'COMMON.LOCAL' logical :: lprint integer :: i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2 real(kind=8) :: rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,& xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,& vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,& eesij,ees,evdw,ene real(kind=8),dimension(2,2) :: elpp6c=reshape((/-0.2379d0,& -0.2056d0,-0.2056d0,-0.0610d0/),shape(elpp6c)) real(kind=8),dimension(2,2) :: elpp3c=reshape((/ 0.0503d0,& 0.0000d0, 0.0000d0, 0.0692d0/),shape(elpp3c)) real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc real(kind=8) :: elcutoff=-0.3d0 real(kind=8) :: elecutoff_14=-0.5d0 integer :: ncont,icont(2,maxcont) real(kind=8) :: econt(maxcont) ! ! Load the constants of peptide bond - peptide bond interactions. ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g. ! proline) - determined by averaging ECEPP energy. ! ! as of 7/06/91. ! ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ ! data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ !el data (elpp6c) /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/ !el data (elpp3c) / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/ !el data (elcutoff) /-0.3d0/ !el data (elecutoff_14) /-0.5d0/ ees=0.0d0 evdw=0.0d0 if (lprint) write (iout,'(a)') & "Constants of electrostatic interaction energy expression." do i=1,2 do j=1,2 rri=rpp(i,j)**6 appc(i,j)=epp(i,j)*rri*rri bppc(i,j)=-2.0*epp(i,j)*rri ael6c(i,j)=elpp6c(i,j)*4.2**6 ael3c(i,j)=elpp3c(i,j)*4.2**3 if (lprint) & write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),& ael3c(i,j) enddo enddo ncont=0 do 1 i=ist,ien-2 xi=c(1,i) yi=c(2,i) zi=c(3,i) dxi=c(1,i+1)-c(1,i) dyi=c(2,i+1)-c(2,i) dzi=c(3,i+1)-c(3,i) xmedi=xi+0.5*dxi ymedi=yi+0.5*dyi zmedi=zi+0.5*dzi do 4 j=i+2,ien-1 ind=ind+1 iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 if (iteli.eq.2 .and. itelj.eq.2 & .or.iteli.eq.0 .or.itelj.eq.0) goto 4 aaa=appc(iteli,itelj) bbb=bppc(iteli,itelj) ael6i=ael6c(iteli,itelj) ael3i=ael3c(iteli,itelj) dxj=c(1,j+1)-c(1,j) dyj=c(2,j+1)-c(2,j) dzj=c(3,j+1)-c(3,j) xj=c(1,j)+0.5*dxj-xmedi yj=c(2,j)+0.5*dyj-ymedi zj=c(3,j)+0.5*dzj-zmedi rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) r3ij=rrmij*rmij r6ij=r3ij*r3ij vrmij=vblinv*rmij cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij fac=cosa-3.0*cosb*cosg ev1=aaa*r6ij*r6ij ev2=bbb*r6ij fac3=ael6i*r6ij fac4=ael3i*r3ij evdwij=ev1+ev2 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg)) el2=fac4*fac eesij=el1+el2 if (j.gt.i+2 .and. eesij.le.elcutoff .or. & j.eq.i+2 .and. eesij.le.elecutoff_14) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j econt(ncont)=eesij endif ees=ees+eesij evdw=evdw+evdwij 4 continue 1 continue if (lprint) then write (iout,*) 'Total average electrostatic energy: ',ees write (iout,*) 'VDW energy between peptide-group centers: ',evdw write (iout,*) write (iout,*) 'Electrostatic contacts before pruning: ' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif ! For given residues keep only the contacts with the greatest energy. i=0 do while (i.lt.ncont) i=i+1 ene=econt(i) ic1=icont(1,i) ic2=icont(2,i) j=i do while (j.lt.ncont) j=j+1 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2, ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then if (ic1.eq.icont(1,j)) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)& .and. iabs(icont(1,k)-ic1).le.2 .and. & econt(k).lt.econt(j) ) goto 21 enddo else if (ic2.eq.icont(2,j) ) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)& .and. iabs(icont(2,k)-ic2).le.2 .and. & econt(k).lt.econt(j) ) goto 21 enddo endif ! Remove ith contact do k=i+1,ncont icont(1,k-1)=icont(1,k) icont(2,k-1)=icont(2,k) econt(k-1)=econt(k) enddo i=i-1 ncont=ncont-1 ! write (iout,*) "ncont",ncont ! do k=1,ncont ! write (iout,*) icont(1,k),icont(2,k) ! enddo goto 20 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) & then if (ic1.eq.icont(1,j)) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. & econt(k).lt.econt(i) ) goto 21 enddo else if (ic2.eq.icont(2,j) ) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. & econt(k).lt.econt(i) ) goto 21 enddo endif ! Remove jth contact do k=j+1,ncont icont(1,k-1)=icont(1,k) icont(2,k-1)=icont(2,k) econt(k-1)=econt(k) enddo ncont=ncont-1 ! write (iout,*) "ncont",ncont ! do k=1,ncont ! write (iout,*) icont(1,k),icont(2,k) ! enddo j=j-1 endif endif 21 continue enddo 20 continue enddo if (lprint) then write (iout,*) write (iout,*) 'Electrostatic contacts after pruning: ' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif return end subroutine elecont !------------------------------------------------------------------------------ ! match_contact.f !------------------------------------------------------------------------------ subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,& ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,& nc_frac,nc_req_set,istr,llocal,lprn) use energy_data, only:maxcont ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: ncont_ref,ncont,ishift,ishif2,nc_match integer,dimension(2,maxcont) :: icont_ref,icont !(2,maxcont) real(kind=8) :: nc_frac logical :: llocal,lprn integer :: ishif1,nc_match1_max,jfrag,n_shif1,n_shif2,& nc_req_set,istr,nc_match_max integer :: i,nc_req,nc_match1,is,js nc_match_max=0 do i=1,ncont_ref nc_match_max=nc_match_max+ & min0(icont_ref(2,i)-icont_ref(1,i)-1,3) enddo if (istr.eq.3) then nc_req=0 else if (nc_req_set.eq.0) then nc_req=nc_match_max*nc_frac else nc_req = dmin1(nc_match_max*nc_frac+0.5d0,& dfloat(nc_req_set)+1.0d-7) endif ! write (iout,*) "match_contact: nc_req:",nc_req ! write (iout,*) "nc_match_max",nc_match_max ! write (iout,*) "jfrag",jfrag," n_shif1",n_shif1, ! & " n_shif2",n_shif2 ! Match current contact map against reference contact map; exit, if at least ! half of the contacts match call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,& ncont,icont,jfrag,llocal,lprn) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",0,0," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. & nc_req.eq.0 .and. nc_match.eq.1) then ishif1=0 ishif2=0 return endif ! If sufficient matches are not found, try to shift contact maps up to three ! positions. if (n_shif1.gt.0) then do is=1,n_shif1 ! The following four tries help to find shifted beta-sheet patterns ! Shift "left" strand backward call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",-is,0," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. & nc_req.eq.0 .and. nc_match.eq.1) then ishif1=-is ishif2=0 return endif ! Shift "left" strand forward call ncont_match(nc_match,nc_match1,is,0,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",is,0," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. & nc_req.eq.0 .and. nc_match.eq.1) then ishif1=is ishif2=0 return endif enddo if (nc_req.eq.0) return ! Shift "right" strand backward do is=1,n_shif1 call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",0,-is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=0 ishif2=-is return endif ! Shift "right" strand upward call ncont_match(nc_match,nc_match1,0,is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",0,is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=0 ishif2=is return endif enddo ! is ! Now try to shift both residues in contacts. do is=1,n_shif1 do js=1,is if (js.ne.is) then call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",-is,-js," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=-is ishif2=-js return endif call ncont_match(nc_match,nc_match1,is,js,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",is,js," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=is ishif2=js return endif ! call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",-js,-is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=-js ishif2=-is return endif ! call ncont_match(nc_match,nc_match1,js,is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",js,is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=js ishif2=is return endif endif ! if (is+js.le.n_shif1) then call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",-is,js," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=-is ishif2=js return endif ! call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",js,-is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=js ishif2=-is return endif endif ! enddo !js enddo !is endif if (n_shif2.gt.0) then do is=1,n_shif2 call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",-is,-is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=-is ishif2=-is return endif call ncont_match(nc_match,nc_match1,is,is,ncont_ref,& icont_ref,ncont,icont,jfrag,llocal,lprn) if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 if (lprn .and. nc_match.gt.0) write (iout,*) & "Shift:",is,is," nc_match1",nc_match1,& " nc_match=",nc_match," req'd",nc_req if (nc_match.ge.nc_req) then ishif1=is ishif2=is return endif enddo endif ! If this point is reached, the contact maps are different. nc_match=0 ishif1=0 ishif2=0 return end subroutine match_contact !------------------------------------------------------------------------- subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,& ncont_ref,icont_ref,ncont,icont,jfrag,llocal,lprn) use energy_data, only:nnt,nct,maxcont ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.INTERACT' ! include 'COMMON.GEO' ! include 'COMMON.COMPAR' logical :: llocal,lprn integer ncont_ref,ncont,ishift,ishif2,nang_pair integer,dimension(2,maxcont) :: icont_ref,icont,icont_match !(2,maxcont) integer,dimension(2,nres) :: iang_pair !(2,maxres) integer :: nc_match,nc_match1,ishif1,jfrag integer :: i,j,ic1,ic2 real(kind=8) :: diffang,fract,rad2deg ! Compare the contact map against the reference contact map; they're stored ! in ICONT and ICONT_REF, respectively. The current contact map can be shifted. if (lprn) write (iout,'(80(1h*))') nc_match=0 nc_match1=0 ! Check the local structure by comparing dihedral angles. ! write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal if (llocal .and. ncont_ref.eq.0) then ! If there are no contacts just compare the dihedral angles and exit. call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,& lprn) if (lprn) write (iout,*) "diffang:",diffang*rad2deg,& " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag)) & then nc_match=1 else nc_match=0 endif return endif nang_pair=0 do i=1,ncont ic1=icont(1,i)+ishif1 ic2=icont(2,i)+ishif2 ! write (iout,*) "i",i," ic1",ic1," ic2",ic2 if (ic1.lt.nnt .or. ic2.gt.nct) goto 10 do j=1,ncont_ref if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3) nc_match1=nc_match1+1 icont_match(1,nc_match1)=ic1 icont_match(2,nc_match1)=ic2 ! call add_angpair(icont(1,i),icont_ref(1,j), ! & nang_pair,iang_pair) ! call add_angpair(icont(2,i),icont_ref(2,j), ! & nang_pair,iang_pair) if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),& " match",icont_ref(1,j),icont_ref(2,j),& " shifts",ishif1,ishif2 goto 10 endif enddo 10 continue enddo if (lprn) then write (iout,*) "nc_match",nc_match," nc_match1",nc_match1 write (iout,*) "icont_match" do i=1,nc_match1 write (iout,*) icont_match(1,i),icont_match(2,i) enddo endif if (llocal .and. nc_match.gt.0) then call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,& ang_cut1(jfrag),diffang,fract) if (lprn) write (iout,*) "diffang:",diffang*rad2deg,& " ang_cut:",ang_cut(jfrag)*rad2deg,& " ang_cut1",ang_cut1(jfrag)*rad2deg if (diffang.gt.ang_cut(jfrag) & .or. fract.lt.frac_min(jfrag)) nc_match=0 endif ! if (nc_match.gt.0) then ! diffang = angnorm1(nang_pair,iang_pair,lprn) ! if (diffang.gt.ang_cut(jfrag)) nc_match=0 ! endif if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,& " diffang",rad2deg*diffang," nc_match",nc_match return end subroutine ncont_match !------------------------------------------------------------------------------ subroutine match_secondary(jfrag,isecstr,nsec_match,lprn) ! This subroutine compares the secondary structure (isecstr) of fragment jfrag ! conformation considered to that of the reference conformation. ! Returns the number of equivalent residues (nsec_match). ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.PEPTCONT' ! include 'COMMON.COMPAR' logical :: lprn integer :: isecstr(nres) integer :: jfrag,nsec_match,npart,i,j npart = npiece(jfrag,1) nsec_match=0 if (lprn) then write (iout,*) "match_secondary jfrag",jfrag," ifrag",& (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart) write (iout,'(80i1)') (isec_ref(j),j=1,nres) write (iout,'(80i1)') (isecstr(j),j=1,nres) endif do i=1,npart do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag) ! The residue has equivalent conformational state to that of the reference ! structure, if: ! a) the conformational states are equal or ! b) the reference state is a coil and that of the conformation considered ! is a strand or ! c) the conformational state of the conformation considered is a strand ! and that of the reference conformation is a coil. ! 10/28/02 - case (b) deleted. if (isecstr(j).eq.isec_ref(j) .or. & ! & isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or. isec_ref(j).eq.0 .and. isecstr(j).eq.1) & nsec_match=nsec_match+1 enddo enddo return end subroutine match_secondary !------------------------------------------------------------------------------ ! odlodc.f !------------------------------------------------------------------------------ subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd) use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,& ! isccont_frag_ref ! implicit real*8 (a-h,o-z) real(kind=8),dimension(3) :: r1,r2,a,b,x,y real(kind=8) :: uu,vv,aa,bb,dd real(kind=8) :: ab,ar,br,det,dd1,dd2,dd3,dd4,dd5 !el odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 & !el + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2 ! print *,"r1",(r1(i),i=1,3) ! print *,"r2",(r2(i),i=1,3) ! print *,"a",(a(i),i=1,3) ! print *,"b",(b(i),i=1,3) aa = a(1)**2+a(2)**2+a(3)**2 bb = b(1)**2+b(2)**2+b(3)**2 ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3) ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3)) br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3)) det = aa*bb-ab**2 ! print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det uu = (-ar*bb+br*ab)/det vv = (br*aa-ar*ab)/det ! print *,u,v uu=dmin1(uu,1.0d0) uu=dmax1(uu,0.0d0) vv=dmin1(vv,1.0d0) vv=dmax1(vv,0.0d0) !el dd1 = odl(uu,vv) dd1 = odl(uu,vv,r1,r2,ar,br,ab,aa,bb) !el dd2 = odl(0.0d0,0.0d0) dd2 = odl(0.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb) !el dd3 = odl(0.0d0,1.0d0) dd3 = odl(0.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb) !el dd4 = odl(1.0d0,0.0d0) dd4 = odl(1.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb) !el dd5 = odl(1.0d0,1.0d0) dd5 = odl(1.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb) dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5)) if (dd.eq.dd2) then uu=0.0d0 vv=0.0d0 else if (dd.eq.dd3) then uu=0.0d0 vv=1.0d0 else if (dd.eq.dd4) then uu=1.0d0 vv=0.0d0 else if (dd.eq.dd5) then uu=1.0d0 vv=1.0d0 endif ! Control check ! do i=1,3 ! x(i)=r1(i)+u*a(i) ! y(i)=r2(i)+v*b(i) ! enddo ! dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2 ! dd1 = dsqrt(dd1) aa = dsqrt(aa) bb = dsqrt(bb) ! write (8,*) uu,vv,dd,dd1 ! print *,dd,dd1 return end subroutine odlodc !------------------------------------------------------------------------------ real(kind=8) function odl(u,v,r1,r2,ar,br,ab,aa,bb) real(kind=8),dimension(3) :: r1,r2 real(kind=8) :: aa,bb,u,v,ar,br,ab odl = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 & + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2 end function odl !------------------------------------------------------------------------------ ! proc_cont.f !------------------------------------------------------------------------------ subroutine proc_cont use geometry_data, only:rad2deg use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,& ! isccont_frag_ref ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.HEADER' ! include 'COMMON.CONTACTS1' ! include 'COMMON.PEPTCONT' ! include 'COMMON.GEO' integer :: i,j,k,ind,len_cut,ndigit,length_frag write (iout,*) "proc_cont: nlevel",nlevel if (nlevel.lt.0) then write (iout,*) "call define_fragments" call define_fragments else write (iout,*) "call secondary2" call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,& isec_ref) endif write (iout,'(80(1h=))') write (iout,*) "Electrostatic contacts" call contacts_between_fragments(.true.,0,ncont_pept_ref,& icont_pept_ref,ncont_frag_ref(1),icont_frag_ref(1,1,1)) write (iout,'(80(1h=))') write (iout,*) "Side chain contacts" call contacts_between_fragments(.true.,0,ncont_ref,& icont_ref,nsccont_frag_ref(1),isccont_frag_ref(1,1,1)) if (nlevel.lt.0) then do i=1,nfrag(1) ind=icant(i,i) len_cut=1000 if (istruct(i).le.1) then len_cut=max0(len_frag(i,1)*4/5,3) else if (istruct(i).eq.2 .or. istruct(i).eq.4) then len_cut=max0(len_frag(i,1)*2/5,3) endif write (iout,*) "i",i," istruct",istruct(i)," ncont_frag",& ncont_frag_ref(ind)," len_cut",len_cut,& " icont_single",icont_single," iloc_single",iloc_single iloc(i)=iloc_single if (iloc(i).gt.0) write (iout,*) & "Local structure used to compare structure of fragment",i,& " to native." if (istruct(i).ne.3 .and. istruct(i).ne.0 & .and. icont_single.gt.0 .and. & ncont_frag_ref(ind).ge.len_cut) then write (iout,*) "Electrostatic contacts used to compare",& " structure of fragment",i," to native." ielecont(i,1)=1 isccont(i,1)=0 else if (icont_single.gt.0 .and. nsccont_frag_ref(ind) & .ge.len_cut) then write (iout,*) "Side chain contacts used to compare",& " structure of fragment",i," to native." isccont(i,1)=1 ielecont(i,1)=0 else write (iout,*) "Contacts not used to compare",& " structure of fragment",i," to native." ielecont(i,1)=0 isccont(i,1)=0 nc_req_setf(i,1)=0 endif if (irms_single.gt.0 .or. isccont(i,1).eq.0 & .and. ielecont(i,1).eq.0) then write (iout,*) "RMSD used to compare",& " structure of fragment",i," to native." irms(i,1)=1 else write (iout,*) "RMSD not used to compare",& " structure of fragment",i," to native." irms(i,1)=0 endif enddo endif if (nlevel.lt.-1) then call define_pairs nlevel = -nlevel if (nlevel.gt.3) nlevel=3 if (nlevel.eq.3) then nfrag(3)=1 npiece(1,3)=nfrag(1) do i=1,nfrag(1) ipiece(i,1,3)=i enddo ielecont(1,3)=0 isccont(1,3)=0 irms(1,3)=1 n_shift(1,1,3)=0 n_shift(2,1,3)=0 endif else if (nlevel.eq.-1) then nlevel=1 endif isnfrag(1)=0 do i=1,nlevel isnfrag(i+1)=isnfrag(i)+nfrag(i) enddo ndigit=3*nfrag(1) do i=2,nlevel ndigit=ndigit+2*nfrag(i) enddo write (iout,*) "ndigit",ndigit if (.not.binary .and. ndigit.gt.30) then write (iout,*) "Highest class too large; switching to",& " binary representation." binary=.true. endif write (iout,*) "isnfrag",(isnfrag(i),i=1,nlevel+1) write(iout,*) "rmscut_base_up",rmscut_base_up,& " rmscut_base_low",rmscut_base_low," rmsup_lim",rmsup_lim do i=1,nlevel do j=1,nfrag(i) length_frag = 0 if (i.eq.1) then do k=1,npiece(j,i) length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1 enddo else do k=1,npiece(j,i) length_frag=length_frag+len_frag(ipiece(k,j,i),1) enddo endif len_frag(j,i)=length_frag rmscutfrag(1,j,i)=rmscut_base_up*length_frag rmscutfrag(2,j,i)=rmscut_base_low*length_frag if (rmscutfrag(1,j,i).lt.rmsup_lim) & rmscutfrag(1,j,i)=rmsup_lim if (rmscutfrag(1,j,i).gt.rmsupup_lim) & rmscutfrag(1,j,i)=rmsupup_lim enddo enddo write (iout,*) "Level",1," number of fragments:",nfrag(1) do j=1,nfrag(1) write (iout,*) npiece(j,1),(ifrag(1,k,j),ifrag(2,k,j),& k=1,npiece(j,1)),len_frag(j,1),rmscutfrag(1,j,1),& rmscutfrag(2,j,1),n_shift(1,j,1),n_shift(2,j,1),& ang_cut(j)*rad2deg,ang_cut1(j)*rad2deg,frac_min(j),& nc_fragm(j,1),nc_req_setf(j,1),istruct(j) enddo do i=2,nlevel write (iout,*) "Level",i," number of fragments:",nfrag(i) do j=1,nfrag(i) write (iout,*) npiece(j,i),(ipiece(k,j,i),& k=1,npiece(j,i)),len_frag(j,i),rmscutfrag(1,j,i),& rmscutfrag(2,j,i),n_shift(1,j,i),n_shift(2,j,i),& nc_fragm(j,i),nc_req_setf(j,i) enddo enddo return end subroutine proc_cont !------------------------------------------------------------------------------ ! define_pairs.f !------------------------------------------------------------------------------ subroutine define_pairs ! use energy_data, only:nsccont_frag_ref ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.COMPAR' ! include 'COMMON.FRAG' ! include 'COMMON.CHAIN' ! include 'COMMON.HEADER' ! include 'COMMON.GEO' ! include 'COMMON.CONTACTS1' ! include 'COMMON.PEPTCONT' integer :: j,k,i,length_frag,ind,ll1,ll2,len_cut do j=1,nfrag(1) length_frag = 0 do k=1,npiece(j,1) length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1 enddo len_frag(j,1)=length_frag write (iout,*) "Fragment",j," length",len_frag(j,1) enddo nfrag(2)=0 do i=1,nfrag(1) do j=i+1,nfrag(1) ind = icant(i,j) if (istruct(i).le.1 .or. istruct(j).le.1) then if (istruct(i).le.1) then ll1=len_frag(i,1) else ll1=len_frag(i,1)/2 endif if (istruct(j).le.1) then ll2=len_frag(j,1) else ll2=len_frag(j,1)/2 endif len_cut=max0(min0(ll1*2/3,ll2*4/5),3) else if (istruct(i).eq.2 .or. istruct(i).eq.4) then ll1=len_frag(i,1)/2 else ll1=len_frag(i,1) endif if (istruct(j).eq.2 .or. istruct(j).eq.4) then ll2=len_frag(j,1)/2 else ll2=len_frag(j,1) endif len_cut=max0(min0(ll1*4/5,ll2)*4/5,3) endif write (iout,*) "Fragments",i,j," structure",istruct(i),& istruct(j)," # contacts",& ncont_frag_ref(ind),nsccont_frag_ref(ind),& " lengths",len_frag(i,1),len_frag(j,1),& " ll1",ll1," ll2",ll2," len_cut",len_cut if ((istruct(i).eq.1 .or. istruct(j).eq.1) .and. & nsccont_frag_ref(ind).ge.len_cut ) then if (istruct(i).eq.1 .and. istruct(j).eq.1) then write (iout,*) "Adding pair of helices",i,j,& " based on SC contacts" else write (iout,*) "Adding helix+strand/sheet pair",i,j,& " based on SC contacts" endif nfrag(2)=nfrag(2)+1 if (icont_pair.gt.0) then write (iout,*) "# SC contacts will be used",& " in comparison." isccont(nfrag(2),2)=1 endif if (irms_pair.gt.0) then write (iout,*) "Fragment RMSD will be used",& " in comparison." irms(nfrag(2),2)=1 endif npiece(nfrag(2),2)=2 ipiece(1,nfrag(2),2)=i ipiece(2,nfrag(2),2)=j ielecont(nfrag(2),2)=0 n_shift(1,nfrag(2),2)=nshift_pair n_shift(2,nfrag(2),2)=nshift_pair nc_fragm(nfrag(2),2)=ncfrac_pair nc_req_setf(nfrag(2),2)=ncreq_pair else if ((istruct(i).ge.2 .and. istruct(i).le.4) & .and. (istruct(j).ge.2 .and. istruct(i).le.4) & .and. ncont_frag_ref(ind).ge.len_cut ) then nfrag(2)=nfrag(2)+1 write (iout,*) "Adding pair strands/sheets",i,j,& " based on pp contacts" if (icont_pair.gt.0) then write (iout,*) "# pp contacts will be used",& " in comparison." ielecont(nfrag(2),2)=1 endif if (irms_pair.gt.0) then write (iout,*) "Fragment RMSD will be used",& " in comparison." irms(nfrag(2),2)=1 endif npiece(nfrag(2),2)=2 ipiece(1,nfrag(2),2)=i ipiece(2,nfrag(2),2)=j ielecont(nfrag(2),2)=1 isccont(nfrag(2),2)=0 n_shift(1,nfrag(2),2)=nshift_pair n_shift(2,nfrag(2),2)=nshift_pair nc_fragm(nfrag(2),2)=ncfrac_bet nc_req_setf(nfrag(2),2)=ncreq_bet endif enddo enddo write (iout,*) "Pairs found" do i=1,nfrag(2) write (iout,*) ipiece(1,i,2),ipiece(2,i,2) enddo return end subroutine define_pairs !------------------------------------------------------------------------------ ! icant.f !------------------------------------------------------------------------------ INTEGER FUNCTION ICANT(I,J) integer :: i,j IF (I.GE.J) THEN ICANT=(I*(I-1))/2+J ELSE ICANT=(J*(J-1))/2+I ENDIF RETURN END FUNCTION ICANT !------------------------------------------------------------------------------ ! mysort.f !------------------------------------------------------------------------------ subroutine imysort(n, m, mm, x, y, z, z1, z2, z3, z4, z5, z6) ! implicit none integer :: n,m,mm integer :: x(m,mm,n),y(n),z(n),z1(2,n),z6(n),xmin,xtemp real(kind=8) :: z2(n),z3(n),z4(n),z5(n) real(kind=8) :: xxtemp integer :: i,j,k,imax do i=1,n xmin=x(1,1,i) imax=i do j=i+1,n if (x(1,1,j).lt.xmin) then imax=j xmin=x(1,1,j) endif enddo xxtemp=z2(imax) z2(imax)=z2(i) z2(i)=xxtemp xxtemp=z3(imax) z3(imax)=z3(i) z3(i)=xxtemp xxtemp=z4(imax) z4(imax)=z4(i) z4(i)=xxtemp xxtemp=z5(imax) z5(imax)=z5(i) z5(i)=xxtemp xtemp=y(imax) y(imax)=y(i) y(i)=xtemp xtemp=z(imax) z(imax)=z(i) z(i)=xtemp xtemp=z6(imax) z6(imax)=z6(i) z6(i)=xtemp do j=1,2 xtemp=z1(j,imax) z1(j,imax)=z1(j,i) z1(j,i)=xtemp enddo do j=1,m do k=1,mm xtemp=x(j,k,imax) x(j,k,imax)=x(j,k,i) x(j,k,i)=xtemp enddo enddo enddo return end subroutine imysort !------------------------------------------------------------------------------ ! qwolynes.f !------------------------------------------------------------------------------- real(kind=8) function qwolynes(ilevel,jfrag) use geometry_data, only:cref,nperm use control_data, only:symetr use energy_data, only:nnt,nct,itype ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.CONTROL' integer :: ilevel,jfrag,kkk integer :: i,j,jl,k,l,il,kl,nl,np,ip,kp integer :: nsep=3 real(kind=8),dimension(:),allocatable :: tempus !(maxperm) real(kind=8) :: maxiQ !dist, real(kind=8) :: qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM logical :: lprn=.false. real(kind=8) :: x !el sigm !el sigm(x)=0.25d0*x nperm=1 maxiQ=0 do i=1,symetr nperm=i*nperm enddo ! write (iout,*) "QWolyes: " jfrag",jfrag, ! & " ilevel",ilevel allocate(tempus(nperm)) do kkk=1,nperm qq = 0.0d0 if (ilevel.eq.0) then if (lprn) write (iout,*) "Q computed for whole molecule" nl=0 do il=nnt+nsep,nct do jl=nnt,il-nsep dij=0.0d0 dijCM=0.0d0 d0ij=0.0d0 d0ijCM=0.0d0 qqij=0.0d0 qqijCM=0.0d0 nl=nl+1 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ & (cref(2,jl,kkk)-cref(2,il,kkk))**2+ & (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) if (itype(il).ne.10 .or. itype(jl).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2) dijCM=dist(il+nres,jl+nres) qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2) endif qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "il",il," jl",jl,& " itype",itype(il),itype(jl) write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,& " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif enddo enddo qq = qq/nl if (lprn) write (iout,*) "nl",nl," qq",qq else if (ilevel.eq.1) then if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag nl=0 ! write (iout,*) "nlist_frag",nlist_frag(jfrag) do i=2,nlist_frag(jfrag) do j=1,i-1 il=list_frag(i,jfrag) jl=list_frag(j,jfrag) if (iabs(il-jl).gt.nsep) then dij=0.0d0 dijCM=0.0d0 d0ij=0.0d0 d0ijCM=0.0d0 qqij=0.0d0 qqijCM=0.0d0 nl=nl+1 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ & (cref(2,jl,kkk)-cref(2,il,kkk))**2+ & (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) if (itype(il).ne.10 .or. itype(jl).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2) dijCM=dist(il+nres,jl+nres) qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2) endif qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "i",i," j",j," il",il," jl",jl,& " itype",itype(il),itype(jl) write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,& " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif endif enddo enddo qq = qq/nl if (lprn) write (iout,*) "nl",nl," qq",qq else if (ilevel.eq.2) then np=npiece(jfrag,ilevel) nl=0 do i=2,np ip=ipiece(i,jfrag,ilevel) do j=1,nlist_frag(ip) il=list_frag(j,ip) do k=1,i-1 kp=ipiece(k,jfrag,ilevel) do l=1,nlist_frag(kp) kl=list_frag(l,kp) if (iabs(kl-il).gt.nsep) then nl=nl+1 dij=0.0d0 dijCM=0.0d0 d0ij=0.0d0 d0ijCM=0.0d0 qqij=0.0d0 qqijCM=0.0d0 d0ij=dsqrt((cref(1,kl,kkk)-cref(1,il,kkk))**2+ & (cref(2,kl,kkk)-cref(2,il,kkk))**2+ & (cref(3,kl,kkk)-cref(3,il,kkk))**2) dij=dist(il,kl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) if (itype(il).ne.10 .or. itype(kl).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ & (cref(2,kl+nres,kkk)-cref(2,il+nres,kkk))**2+ & (cref(3,kl+nres,kkk)-cref(3,il+nres,kkk))**2) dijCM=dist(il+nres,kl+nres) qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/ & (sigm(d0ijCM)))**2) endif qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "i",i," j",j," k",k," l",l," il",il,& " kl",kl," itype",itype(il),itype(kl) write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",& d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif endif enddo ! l enddo ! k enddo ! j enddo ! i qq = qq/nl if (lprn) write (iout,*) "nl",nl," qq",qq else write (iout,*)"Error: Q can be computed only for level 1 and 2." endif tempus(kkk)=qq enddo do kkk=1,nperm if (maxiQ.le.tempus(kkk)) maxiQ=tempus(kkk) enddo qwolynes=1.0d0-maxiQ deallocate(tempus) return end function qwolynes !------------------------------------------------------------------------------- real(kind=8) function sigm(x) real(kind=8) :: x sigm=0.25d0*x return end function sigm !------------------------------------------------------------------------------- subroutine fragment_list ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' logical :: lprn=.true. integer :: i,ilevel,j,k,jfrag do jfrag=1,nfrag(1) nlist_frag(jfrag)=0 do i=1,npiece(jfrag,1) if (lprn) write (iout,*) "jfrag=",jfrag,& "i=",i," fragment",ifrag(1,i,jfrag),& ifrag(2,i,jfrag) do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag) do k=1,nlist_frag(jfrag) if (list_frag(k,jfrag).eq.j) goto 10 enddo nlist_frag(jfrag)=nlist_frag(jfrag)+1 list_frag(nlist_frag(jfrag),jfrag)=j enddo 10 continue enddo enddo write (iout,*) "Fragment list" do j=1,nfrag(1) write (iout,*)"Fragment",j," list",(list_frag(k,j),& k=1,nlist_frag(j)) enddo return end subroutine fragment_list !------------------------------------------------------------------------------- real(kind=8) function rmscalc(ishif,i,j,jcon,lprn) use w_comm_local use control_data, only:symetr use geometry_data, only:nperm ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' real(kind=8) :: przes(3),obrot(3,3) !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2) !el logical :: iadded(nres) !el integer :: inumber(2,nres) !el common /ccc/ creff,cc,iadded,inumber logical :: lprn logical :: non_conv integer :: ishif,i,j,jcon,idup,kkk,l,k,kk real(kind=8) :: rminrms,rms if (lprn) then write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif write (iout,*) "npiece",npiece(j,i) call flush(iout) endif ! write (iout,*) "symetr",symetr ! call flush(iout) nperm=1 do idup=1,symetr nperm=nperm*idup enddo ! write (iout,*) "nperm",nperm ! call flush(iout) do kkk=1,nperm idup=0 do l=1,nres iadded(l)=.false. enddo ! write (iout,*) "kkk",kkk ! call flush(iout) do k=1,npiece(j,i) if (i.eq.1) then if (lprn) then write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",& ifrag(1,k,j),ifrag(2,k,j) call flush(iout) endif call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk) ! write (iout,*) "Exit cprep" ! call flush(iout) ! write (iout,*) "ii=",ii else kk = ipiece(k,j,i) ! write (iout,*) "kk",kk," npiece",npiece(kk,1) do l=1,npiece(kk,1) if (lprn) then write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,& " l=",l," adding fragment",& ifrag(1,l,kk),ifrag(2,l,kk) call flush(iout) endif call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk) ! write (iout,*) "After cprep" ! call flush(iout) enddo endif enddo enddo if (lprn) then write (iout,*) "tuszukaj" do kkk=1,nperm do k=1,idup write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),& inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3) enddo enddo call flush(iout) endif rminrms=1.0d10 do kkk=1,nperm call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv) if (non_conv) then print *,'Error: FITSQ non-convergent, jcon',jcon,i rms = 1.0d10 else if (rms.lt.-1.0d-6) then print *,'Error: rms^2 = ',rms,jcon,i rms = 1.0d10 else if (rms.ge.1.0d-6 .and. rms.lt.0) then rms = 0.0d0 endif ! write (iout,*) "rmsmin", rminrms, "rms", rms if (rms.le.rminrms) rminrms=rms enddo rmscalc = dsqrt(rminrms) ! write (iout, *) "analysys", rmscalc,anatemp return end function rmscalc !------------------------------------------------------------------------- subroutine cprep(if1,if2,ishif,idup,kwa) use w_comm_local use control_data, only:symetr use geometry_data, only:nperm,cref,c ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.CONTROL' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' real(kind=8) :: przes(3),obrot(3,3) !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2) !el logical :: iadded(nres) !el integer :: inumber(2,nres) integer :: iistrart,kwa,blar !el common /ccc/ creff,cc,iadded,inumber integer :: if1,if2,ishif,idup,kkk,l,m ! write (iout,*) "Calling cprep symetr",symetr," kwa",kwa nperm=1 do blar=1,symetr nperm=nperm*blar enddo ! write (iout,*) "nperm",nperm kkk=kwa ! ii=0 do l=if1,if2 ! write (iout,*) "l",l," iadded",iadded(l) ! call flush(iout) if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) & then idup=idup+1 iadded(l)=.true. inumber(1,idup)=l inumber(2,idup)=l+ishif do m=1,3 creff(m,idup)=cref(m,l,kkk) cc(m,idup)=c(m,l+ishif) enddo endif enddo return end subroutine cprep !------------------------------------------------------------------------- real(kind=8) function rmsnat(jcon) use control_data, only:symetr use geometry_data, only:nperm,cref,c use energy_data, only:itype ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' real(kind=8) :: przes(3),obrot(3,3),cc(3,2*nres),ccref(3,2*nres) logical :: non_conv integer :: ishif,i,j,resprzesun,jcon,kkk,nnsup real(kind=8) :: rminrms,rmsminsing,rms rminrms=10.0d10 rmsminsing=10d10 nperm=1 do i=1,symetr nperm=nperm*i enddo do kkk=1,nperm nnsup=0 do i=1,nres if (itype(i).ne.ntyp1) then nnsup=nnsup+1 do j=1,3 cc(j,nnsup)=c(j,i) ccref(j,nnsup)=cref(j,i,kkk) enddo endif enddo call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv) if (non_conv) then print *,'Error: FITSQ non-convergent, jcon',jcon,i rms=1.0d10 else if (rms.lt.-1.0d-6) then print *,'Error: rms^2 = ',rms,jcon,i rms = 1.0d10 else if (rms.ge.1.0d-6 .and. rms.lt.0) then rms=0.0d0 endif if (rms.le.rminrms) rminrms=rms ! write (iout,*) "kkk",kkk," rmsnat",rms , rminrms enddo rmsnat = dsqrt(rminrms) ! write (iout,*) "analysys",rmsnat, anatemp ! liczenie rmsdla pojedynczego lancucha return end function rmsnat !------------------------------------------------------------------------------- subroutine define_fragments use geometry_data, only:rad2deg use energy_data, only:itype use compare_data, only:nhfrag,nbfrag,bfrag,hfrag ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' ! include 'COMMON.FRAG' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.COMPAR' ! include 'COMMON.CHAIN' ! include 'COMMON.HEADER' ! include 'COMMON.GEO' ! include 'COMMON.CONTACTS' ! include 'COMMON.PEPTCONT' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' integer :: nstrand,istrand(2,nres/2) integer :: nhairp,ihairp(2,nres/5) character(len=16) :: strstr(4)=reshape((/'helix','hairpin',& 'strand','strand pair'/),shape(strstr)) integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,& 'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,& 'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,& ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,& ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,& ' MAXANG_HEL',angcut1_hel*rad2deg write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,& ' MAXANG_BET',angcut1_bet*rad2deg write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,& ' MAXANG_STRAND',angcut1_strand*rad2deg write (iout,*) 'FRAC_MIN',frac_min_set ! Find secondary structure elements (helices and beta-sheets) call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,& isec_ref) ! Define primary fragments. First include the helices. nhairp=0 nstrand=0 ! Merge helices ! AL 12/23/03 - to avoid splitting helices into very small fragments if (merge_helices) then write (iout,*) "Before merging helices: nhfrag",nhfrag do i=1,nhfrag write (2,*) hfrag(1,i),hfrag(2,i) enddo i=1 do while (i.lt.nhfrag) if (hfrag(1,i+1)-hfrag(2,i).le.1) then nhfrag=nhfrag-1 hfrag(2,i)=hfrag(2,i+1) do j=i+1,nhfrag hfrag(1,j)=hfrag(1,j+1) hfrag(2,j)=hfrag(2,j+1) enddo endif i=i+1 enddo write (iout,*) "After merging helices: nhfrag",nhfrag do i=1,nhfrag write (2,*) hfrag(1,i),hfrag(2,i) enddo endif nfrag(1)=nhfrag do i=1,nhfrag npiece(i,1)=1 ifrag(1,1,i)=hfrag(1,i) ifrag(2,1,i)=hfrag(2,i) n_shift(1,i,1)=0 n_shift(2,i,1)=nshift_hel ang_cut(i)=angcut_hel ang_cut1(i)=angcut1_hel frac_min(i)=frac_min_set nc_fragm(i,1)=ncfrac_hel nc_req_setf(i,1)=ncreq_hel istruct(i)=1 enddo write (iout,*) "isplit_bet",isplit_bet if (isplit_bet.gt.1) then ! Split beta-sheets into strands and store strands as primary fragments. call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp) do i=1,nstrand ii=i+nfrag(1) npiece(ii,1)=1 ifrag(1,1,ii)=istrand(1,i) ifrag(2,1,ii)=istrand(2,i) n_shift(1,ii,1)=nshift_strand n_shift(2,ii,1)=nshift_strand ang_cut(ii)=angcut_strand ang_cut1(ii)=angcut1_strand frac_min(ii)=frac_min_set nc_fragm(ii,1)=0 nc_req_setf(ii,1)=0 istruct(ii)=3 enddo nfrag(1)=nfrag(1)+nstrand else if (isplit_bet.eq.1) then ! Split only far beta-sheets; does not split hairpins. call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp) call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp) do i=1,nhairp ii=i+nfrag(1) npiece(ii,1)=1 ifrag(1,1,ii)=ihairp(1,i) ifrag(2,1,ii)=ihairp(2,i) n_shift(1,ii,1)=nshift_bet n_shift(2,ii,1)=nshift_bet ang_cut(ii)=angcut_bet ang_cut1(ii)=angcut1_bet frac_min(ii)=frac_min_set nc_fragm(ii,1)=ncfrac_bet nc_req_setf(ii,1)=ncreq_bet istruct(ii)=2 enddo nfrag(1)=nfrag(1)+nhairp do i=1,nstrand ii=i+nfrag(1) npiece(ii,1)=1 ifrag(1,1,ii)=istrand(1,i) ifrag(2,1,ii)=istrand(2,i) n_shift(1,ii,1)=nshift_strand n_shift(2,ii,1)=nshift_strand ang_cut(ii)=angcut_strand ang_cut1(ii)=angcut1_strand frac_min(ii)=frac_min_set nc_fragm(ii,1)=0 nc_req_setf(ii,1)=0 istruct(ii)=3 enddo nfrag(1)=nfrag(1)+nstrand else ! Do not split beta-sheets; each pair of strands is a primary element. call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp) do i=1,nhairp ii=i+nfrag(1) npiece(ii,1)=1 ifrag(1,1,ii)=ihairp(1,i) ifrag(2,1,ii)=ihairp(2,i) n_shift(1,ii,1)=nshift_bet n_shift(2,ii,1)=nshift_bet ang_cut(ii)=angcut_bet ang_cut1(ii)=angcut1_bet frac_min(ii)=frac_min_set nc_fragm(ii,1)=ncfrac_bet nc_req_setf(ii,1)=ncreq_bet istruct(ii)=2 enddo nfrag(1)=nfrag(1)+nhairp do i=1,nbfrag ii=i+nfrag(1) npiece(ii,1)=2 ifrag(1,1,ii)=bfrag(1,i) ifrag(2,1,ii)=bfrag(2,i) if (bfrag(3,i).lt.bfrag(4,i)) then ifrag(1,2,ii)=bfrag(3,i) ifrag(2,2,ii)=bfrag(4,i) else ifrag(1,2,ii)=bfrag(4,i) ifrag(2,2,ii)=bfrag(3,i) endif n_shift(1,ii,1)=nshift_bet n_shift(2,ii,1)=nshift_bet ang_cut(ii)=angcut_bet ang_cut1(ii)=angcut1_bet frac_min(ii)=frac_min_set nc_fragm(ii,1)=ncfrac_bet nc_req_setf(ii,1)=ncreq_bet istruct(ii)=4 enddo nfrag(1)=nfrag(1)+nbfrag endif write (iout,*) "The following primary fragments were found:" write (iout,*) "Helices:",nhfrag do i=1,nhfrag i1=ifrag(1,1,i) i2=ifrag(2,1,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo write (iout,*) "Hairpins:",nhairp do i=nhfrag+1,nhfrag+nhairp i1=ifrag(1,1,i) i2=ifrag(2,1,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') & i,restyp(it1),i1,restyp(it2),i2 enddo write (iout,*) "Far strand pairs:",nbfrag do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag i1=ifrag(1,1,i) i2=ifrag(2,1,i) it1=itype(i1) it2=itype(i2) i3=ifrag(1,2,i) i4=ifrag(2,2,i) it3=itype(i3) it4=itype(i4) write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2,& restyp(it3),i3,restyp(it4),i4 enddo write (iout,*) "Strands:",nstrand do i=nhfrag+nhairp+nbfrag+1,nfrag(1) i1=ifrag(1,1,i) i2=ifrag(2,1,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),& istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),& nc_fragm(1,1),nc_req_setf(1,1)) write (iout,*) "Fragments after sorting:" do i=1,nfrag(1) i1=ifrag(1,1,i) i2=ifrag(2,1,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,$)') & i,restyp(it1),i1,restyp(it2),i2 if (npiece(i,1).eq.1) then write (iout,'(2x,a)') strstr(istruct(i)) else i1=ifrag(1,2,i) i2=ifrag(2,2,i) it1=itype(i1) it2=itype(i2) write (iout,'(2x,a,i4,2x,a,i4,2x,a)') & restyp(it1),i1,restyp(it2),i2,strstr(istruct(i)) endif enddo return end subroutine define_fragments !------------------------------------------------------------------------------ subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' integer :: nbfrag,bfrag(4,nres/3) integer :: nhairp,ihairp(2,nres/5) integer :: i,j,k write (iout,*) "Entered find_and_remove_hairpins" write (iout,*) "nbfrag",nbfrag do i=1,nbfrag write (iout,*) i,(bfrag(k,i),k=1,4) enddo nhairp=0 i=1 do while (i.le.nbfrag) write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4) if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) & then write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i) nhairp=nhairp+1 ihairp(1,nhairp)=bfrag(1,i) ihairp(2,nhairp)=bfrag(3,i) nbfrag=nbfrag-1 do j=i,nbfrag do k=1,4 bfrag(k,j)=bfrag(k,j+1) enddo enddo else i=i+1 endif enddo write (iout,*) "After finding hairpins:" write (iout,*) "nhairp",nhairp do i=1,nhairp write (iout,*) i,ihairp(1,i),ihairp(2,i) enddo write (iout,*) "nbfrag",nbfrag do i=1,nbfrag write (iout,*) i,(bfrag(k,i),k=1,4) enddo return end subroutine find_and_remove_hairpins !------------------------------------------------------------------------------ subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' integer :: nbfrag,bfrag(4,nres/3) integer :: nstrand,istrand(2,nres/2) integer :: nhairp,ihairp(2,nres/5) logical :: found integer :: i,k write (iout,*) "Entered split_beta" write (iout,*) "nbfrag",nbfrag do i=1,nbfrag write (iout,*) i,(bfrag(k,i),k=1,4) enddo nstrand=0 do i=1,nbfrag write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i) call add_strand(nstrand,istrand,nhairp,ihairp,& bfrag(1,i),bfrag(2,i),found) if (bfrag(3,i).lt.bfrag(4,i)) then write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i) call add_strand(nstrand,istrand,nhairp,ihairp,& bfrag(3,i),bfrag(4,i),found) else write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i) call add_strand(nstrand,istrand,nhairp,ihairp,& bfrag(4,i),bfrag(3,i),found) endif enddo nbfrag=0 write (iout,*) "Strands found:",nstrand do i=1,nstrand write (iout,*) i,istrand(1,i),istrand(2,i) enddo return end subroutine split_beta !------------------------------------------------------------------------------ subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' ! include 'COMMON.IOUNITS' integer :: nstrand,istrand(2,nres/2) integer :: nhairp,ihairp(2,nres/5) logical :: found integer :: is1,is2,j,idelt found=.false. do j=1,nhairp idelt=(ihairp(2,j)-ihairp(1,j))/6 if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then write (iout,*) "strand",is1,is2," is part of hairpin",& ihairp(1,j),ihairp(2,j) return endif enddo do j=1,nstrand idelt=(istrand(2,j)-istrand(1,j))/3 if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) & then ! The strand already exists in the array; update its ends if necessary. write (iout,*) "strand",is1,is2," found at position",j,& ":",istrand(1,j),istrand(2,j) istrand(1,j)=min0(istrand(1,j),is1) istrand(2,j)=max0(istrand(2,j),is2) return endif enddo ! The strand has not been found; add it to the array. write (iout,*) "strand",is1,is2," added to the array." found=.true. nstrand=nstrand+1 istrand(1,nstrand)=is1 istrand(2,nstrand)=is2 return end subroutine add_strand !------------------------------------------------------------------------------ subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr) use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup use energy_data, only:itype,maxcont use compare_data, only:bfrag,hfrag,nbfrag,nhfrag use compare, only:freeres ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'COMMON.IOUNITS' ! include 'COMMON.FRAG' ! include 'COMMON.VAR' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.NAMES' ! include 'COMMON.INTERACT' integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres),& isecstr(nres) logical :: lprint,lprint_sec,not_done !el,freeres integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist integer :: nstrand,nbeta,nhelix,iii1,jjj1 real(kind=8) :: p1,p2 !rel external freeres character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec)) if (lprint) then write (iout,*) "entered secondary2",ncont write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup do i=1,ncont write (iout,*) icont(1,i),icont(2,i) enddo endif do i=1,nres isecstr(i)=0 enddo nbfrag=0 nhfrag=0 do i=1,nres isec(i,1)=0 isec(i,2)=0 nsec(i)=0 enddo ! finding parallel beta !d write (iout,*) '------- looking for parallel beta -----------' nbeta=0 nstrand=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (i1.ge.nstart_sup .and. i1.le.nend_sup & .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then !d write (iout,*) "parallel",i1,j1 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1 not_done=.true. do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. & freeres(i1,j1,nsec,isec)) goto 5 enddo not_done=.false. 5 continue !d write (iout,*) i1,j1,not_done enddo j1=j1-1 i1=i1-1 if (i1-ii1.gt.1) then ii1=max0(ii1-1,1) jj1=max0(jj1-1,1) nbeta=nbeta+1 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',& nbeta,ii1,i1,jj1,j1 nbfrag=nbfrag+1 bfrag(1,nbfrag)=ii1+1 bfrag(2,nbfrag)=i1+1 bfrag(3,nbfrag)=jj1+1 bfrag(4,nbfrag)=min0(j1+1,nres) do ij=ii1,i1 nsec(ij)=nsec(ij)+1 isec(ij,nsec(ij))=nbeta enddo do ij=jj1,j1 nsec(ij)=nsec(ij)+1 isec(ij,nsec(ij))=nbeta enddo if(lprint_sec) then nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" endif nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1 endif endif endif endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup enddo ! finding antiparallel beta !d write (iout,*) '--------- looking for antiparallel beta ---------' do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (freeres(i1,j1,nsec,isec)) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1 not_done=.true. do while (not_done) i1=i1+1 j1=j1-1 do j=1,ncont if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. & freeres(i1,j1,nsec,isec)) goto 6 enddo not_done=.false. 6 continue !d write (iout,*) i1,j1,not_done enddo i1=i1-1 j1=j1+1 if (i1-ii1.gt.1) then nbfrag=nbfrag+1 bfrag(1,nbfrag)=ii1 bfrag(2,nbfrag)=min0(i1+1,nres) bfrag(3,nbfrag)=min0(jj1+1,nres) bfrag(4,nbfrag)=j1 nbeta=nbeta+1 iii1=max0(ii1-1,1) do ij=iii1,i1 nsec(ij)=nsec(ij)+1 if (nsec(ij).le.2) then isec(ij,nsec(ij))=nbeta endif enddo jjj1=max0(j1-1,1) do ij=jjj1,jj1 nsec(ij)=nsec(ij)+1 if (nsec(ij).le.2) then isec(ij,nsec(ij))=nbeta endif enddo if (lprint_sec) then write (iout,'(a,i3,4i4)')'antiparallel beta',& nbeta,ii1-1,i1,jj1,j1-1 nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" endif nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2 endif endif endif enddo !d write (iout,*) "After beta:",nbfrag !d do i=1,nbfrag !d write (iout,*) (bfrag(j,i),j=1,4) !d enddo if (nstrand.gt.0.and.lprint_sec) then write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1" do i=2,nstrand if (i.le.9) then write(12,'(a9,i1,$)') " | strand",i else write(12,'(a9,i2,$)') " | strand",i endif enddo write(12,'(a1)') "'" endif ! finding alpha or 310 helix nhelix=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) p1=phi(i1+2)*rad2deg p2=0.0 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg if (j1.eq.i1+3 .and. & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2 ii1=i1 jj1=j1 if (nsec(ii1).eq.0) then not_done=.true. else not_done=.false. endif do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10 enddo not_done=.false. 10 continue p1=phi(i1+2)*rad2deg p2=phi(j1+2)*rad2deg if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) & not_done=.false. !d write (iout,*) i1,j1,not_done,p1,p2 enddo j1=j1+1 if (j1-ii1.gt.4) then nhelix=nhelix+1 !d write (iout,*)'helix',nhelix,ii1,j1 nhfrag=nhfrag+1 hfrag(1,nhfrag)=ii1 hfrag(2,nhfrag)=j1 do ij=ii1,j1 nsec(ij)=-1 enddo if (lprint_sec) then write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1 if (nhelix.le.9) then write(12,'(a17,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" else write(12,'(a17,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" endif endif endif endif enddo if (nhelix.gt.0.and.lprint_sec) then write(12,'(a26,$)') "DefPropRes 'helix' 'helix1" do i=2,nhelix if (nhelix.le.9) then write(12,'(a8,i1,$)') " | helix",i else write(12,'(a8,i2,$)') " | helix",i endif enddo write(12,'(a1)') "'" endif if (lprint_sec) then write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'" write(12,'(a20)') "XMacStand ribbon.mac" endif if (lprint) then write(iout,*) 'UNRES seq:',anatemp do j=1,nbfrag write(iout,*) 'beta ',(bfrag(i,j),i=1,4) enddo do j=1,nhfrag write(iout,*) 'helix ',(hfrag(i,j),i=1,2),anatemp enddo endif do j=1,nbfrag do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j)) isecstr(k)=1 enddo do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j)) isecstr(k)=1 enddo enddo do j=1,nhfrag do k=hfrag(1,j),hfrag(2,j) isecstr(k)=2 enddo enddo if (lprint) then write (iout,*) write (iout,*) "Secondary structure" do i=1,nres,80 ist=i ien=min0(i+79,nres) write (iout,*) write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10) write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien) write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien) enddo write (iout,*) endif return end subroutine secondary2 !------------------------------------------------- ! logical function freeres(i,j,nsec,isec) ! include 'DIMENSIONS' ! integer :: isec(nres,4),nsec(nres) ! integer :: i,j,k,l ! freeres=.false. ! ! if (nsec(i).gt.1.or.nsec(j).gt.1) return ! do k=1,nsec(i) ! do l=1,nsec(j) ! if (isec(i,k).eq.isec(j,l)) return ! enddo ! enddo ! freeres=.true. ! return ! end function freeres !------------------------------------------------- subroutine alloc_compar_arrays(nfrg,nlev) use energy_data, only:maxcont use w_comm_local integer :: nfrg,nlev !write(iout,*) "in alloc conpar arrays: nlevel=", nlevel," nfrag(1)=",nfrag(1) !------------------------ ! commom.contacts ! common /contacts/ allocate(nsccont_frag_ref(mmaxfrag)) !(mmaxfrag) !wham allocate(isccont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag) !wham !------------------------ ! COMMON.COMPAR ! common /compar/ allocate(rmsfrag(nfrg,nlev+1),nc_fragm(nfrg,nlev+1)) !(maxfrag,maxlevel) allocate(qfrag(nfrg,2)) !(maxfrag,2) allocate(rmscutfrag(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel) allocate(ang_cut(nfrg),ang_cut1(nfrg),frac_min(nfrg)) !(maxfrag) allocate(nc_req_setf(nfrg,nlev+1),npiece(nfrg,nlev+1),& ielecont(nfrg,nlev+1),isccont(nfrg,nlev+1),irms(nfrg,nlev+1),& ishifft(nfrg,nlev+1),len_frag(nfrg,nlev+1)) !(maxfrag,maxlevel) allocate(ncont_nat(2,nfrg,nlev+1)) allocate(n_shift(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel) ! allocate(nfrag(nlev)) !(maxlevel) allocate(isnfrag(nlev+2)) !(maxlevel+1) allocate(ifrag(2,maxpiece,nfrg)) !(2,maxpiece,maxfrag) allocate(ipiece(maxpiece,nfrg,2:nlev+1)) !(maxpiece,maxfrag,2:maxlevel) allocate(istruct(nfrg),iloc(nfrg),nlist_frag(nfrg)) !(maxfrag) allocate(iclass(nlev*nfrg,nlev+1)) !(maxlevel*maxfrag,maxlevel) allocate(list_frag(nres,nfrg)) !(maxres,maxfrag) !------------------------ ! COMMON.PEPTCONT ! common /peptcont/ ! integer,dimension(:,:),allocatable :: icont_pept_ref !(2,maxcont) allocate(ncont_frag_ref(mmaxfrag)) !(mmaxfrag) allocate(icont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag) ! integer,dimension(:),allocatable :: isec_ref !(maxres) !------------------------ ! module w_comm_local ! common /ccc/ allocate(creff(3,2*nres),cc(3,2*nres)) !(3,nres*2) allocate(iadded(nres)) !(nres) allocate(inumber(2,nres)) !(2,nres) !------------------------------------------------------------------------------- end subroutine alloc_compar_arrays #endif !------------------------------------------------------------------------------- end module conform_compar