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