+++ /dev/null
- 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
-
projects / unres4.git / blobdiff