restoring read2sigma code after wrong merge
[unres.git] / source / wham / src / contact.f
1       subroutine contact(lprint,ncont,icont,ist,ien)
2       implicit none
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.ZSCOPT'
5       include 'COMMON.CONTROL'
6       include 'COMMON.IOUNITS'
7       include 'COMMON.CHAIN'
8       include 'COMMON.INTERACT'
9       include 'COMMON.FFIELD'
10       include 'COMMON.NAMES'
11       include 'COMMON.CALC'
12       include 'COMMON.CONTPAR'
13       include 'COMMON.LOCAL'
14       integer ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
15       real*8 csc,dist
16       real*8 cscore(maxcont),omt1(maxcont),omt2(maxcont),omt12(maxcont),
17      & ddsc(maxcont),ddla(maxcont),ddlb(maxcont)
18       integer ncont,icont(2,maxcont)
19       real*8 u,v,a(3),b(3),dla,dlb
20       logical lprint
21       ncont=0
22       kkk=3
23       if (lprint) then
24       do i=1,nres
25         write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),
26      &    c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),
27      &    dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
28       enddo
29       endif
30   110 format (a,'(',i3,')',9f8.3)
31       do i=ist,ien-kkk
32         iti=itype(i)
33         do j=i+kkk,ien
34           itj=itype(j)
35           itypi=iti
36           itypj=itj
37           xj = c(1,nres+j)-c(1,nres+i)    
38           yj = c(2,nres+j)-c(2,nres+i)    
39           zj = c(3,nres+j)-c(3,nres+i)    
40           dxi = dc_norm(1,nres+i)
41           dyi = dc_norm(2,nres+i)
42           dzi = dc_norm(3,nres+i)
43           dxj = dc_norm(1,nres+j)
44           dyj = dc_norm(2,nres+j)
45           dzj = dc_norm(3,nres+j)
46           do k=1,3
47             a(k)=dc(k,nres+i)
48             b(k)=dc(k,nres+j)
49           enddo
50 c          write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
51           if (icomparfunc.eq.1) then
52             call contfunc(csc,iti,itj)
53           else if (icomparfunc.eq.2) then
54             call scdist(csc,iti,itj)
55           else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
56             csc = dist(nres+i,nres+j)
57           else if (icomparfunc.eq.4) then
58             call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
59           else
60             write (*,*) "Error - Unknown sidechain contact function"
61             write (iout,*) "Error - Unknown sidechain contact function"
62           endif
63           if (csc.lt.sc_cutoff(iti,itj)) then
64 c            write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
65 c     &      " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
66 c     &      dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
67 c     &      dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
68 c     &      xj,yj,zj
69 c            write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
70 c     &       sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
71 c     &       chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
72 c     &       csc
73             ncont=ncont+1
74             cscore(ncont)=csc
75             icont(1,ncont)=i
76             icont(2,ncont)=j
77             omt1(ncont)=om1
78             omt2(ncont)=om2
79             omt12(ncont)=om12
80             ddsc(ncont)=1.0d0/rij
81             ddla(ncont)=dla
82             ddlb(ncont)=dlb
83           endif
84         enddo
85       enddo
86       if (lprint) then
87         write (iout,'(a)') 'Contact map:'
88         do i=1,ncont
89           i1=icont(1,i)
90           i2=icont(2,i)
91           it1=itype(i1)
92           it2=itype(i2)
93           write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') 
94      &     i,restyp(it1),i1,restyp(it2),i2,cscore(i),
95      &     sc_cutoff(it1,it2),ddsc(i),ddla(i),ddlb(i),
96      &     omt1(i),omt2(i),omt12(i)
97         enddo
98       endif
99       return
100       end
101 c----------------------------------------------------------------------------
102       double precision function contact_fract(ncont,ncont_ref,
103      &                                     icont,icont_ref)
104       implicit none
105       include 'DIMENSIONS'
106       include 'COMMON.IOUNITS'
107       integer i,j,nmatch
108       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
109       nmatch=0
110 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref 
111 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
112 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
113 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
114 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
115       do i=1,ncont
116         do j=1,ncont_ref
117           if (icont(1,i).eq.icont_ref(1,j) .and. 
118      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
119         enddo
120       enddo
121 c     print *,' nmatch=',nmatch
122 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
123       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
124       return
125       end
126 c------------------------------------------------------------------------------
127       subroutine pept_cont(lprint,ncont,icont)
128       implicit none
129       include 'DIMENSIONS'
130       include 'DIMENSIONS.ZSCOPT'
131       include 'COMMON.IOUNITS'
132       include 'COMMON.CHAIN'
133       include 'COMMON.INTERACT'
134       include 'COMMON.FFIELD'
135       include 'COMMON.NAMES'
136       integer ncont,icont(2,maxcont)
137       integer i,j,k,kkk,i1,i2,it1,it2
138       logical lprint
139       real*8 dist
140       real*8 rcomp /5.5d0/
141       ncont=0
142       kkk=0
143       print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
144       do i=nnt,nct-3
145         do k=1,3
146           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
147         enddo
148         do j=i+2,nct-1
149           do k=1,3
150             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
151           enddo
152           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
153             ncont=ncont+1
154             icont(1,ncont)=i
155             icont(2,ncont)=j
156           endif
157         enddo
158       enddo
159       if (lprint) then
160         write (iout,'(a)') 'PP contact map:'
161         do i=1,ncont
162           i1=icont(1,i)
163           i2=icont(2,i)
164           it1=itype(i1)
165           it2=itype(i2)
166           write (iout,'(i3,2x,a,i4,2x,a,i4)')
167      &     i,restyp(it1),i1,restyp(it2),i2
168         enddo
169       endif
170       return
171       end