Merge branch 'bartek' of mmka.chem.univ.gda.pl:unres into bartek
[unres.git] / source / unres / src_MD / contact.f
1       subroutine contact(lprint,ncont,icont,co)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.CHAIN'
6       include 'COMMON.INTERACT'
7       include 'COMMON.FFIELD'
8       include 'COMMON.NAMES'
9       real*8 facont /1.569D0/  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
10       integer ncont,icont(2,maxcont)
11       logical lprint
12       ncont=0
13       kkk=3
14       do i=nnt+kkk,nct
15         iti=itype(i)
16         do j=nnt,i-kkk
17           itj=itype(j)
18           if (ipot.ne.4) then
19 c           rcomp=sigmaii(iti,itj)+1.0D0
20             rcomp=facont*sigmaii(iti,itj)
21           else 
22 c           rcomp=sigma(iti,itj)+1.0D0
23             rcomp=facont*sigma(iti,itj)
24           endif
25 c         rcomp=6.5D0
26 c         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
27           if (dist(nres+i,nres+j).lt.rcomp) then
28             ncont=ncont+1
29             icont(1,ncont)=i
30             icont(2,ncont)=j
31           endif
32         enddo
33       enddo
34       if (lprint) then
35         write (iout,'(a)') 'Contact map:'
36         do i=1,ncont
37           i1=icont(1,i)
38           i2=icont(2,i)
39           it1=itype(i1)
40           it2=itype(i2)
41           write (iout,'(i3,2x,a,i4,2x,a,i4)') 
42      &     i,restyp(it1),i1,restyp(it2),i2 
43         enddo
44       endif
45       co = 0.0d0
46       do i=1,ncont
47         co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
48       enddo 
49       co = co / (nres*ncont)
50       return
51       end
52 c----------------------------------------------------------------------------
53       double precision function contact_fract(ncont,ncont_ref,
54      &                                     icont,icont_ref)
55       implicit real*8 (a-h,o-z)
56       include 'DIMENSIONS'
57       include 'COMMON.IOUNITS'
58       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
59       nmatch=0
60 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref 
61 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
62 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
63 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
64 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
65       do i=1,ncont
66         do j=1,ncont_ref
67           if (icont(1,i).eq.icont_ref(1,j) .and. 
68      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
69         enddo
70       enddo
71 c     print *,' nmatch=',nmatch
72 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
73       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
74       return
75       end
76 c----------------------------------------------------------------------------
77       double precision function contact_fract_nn(ncont,ncont_ref,
78      &                                     icont,icont_ref)
79       implicit real*8 (a-h,o-z)
80       include 'DIMENSIONS'
81       include 'COMMON.IOUNITS'
82       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
83       nmatch=0
84 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref 
85 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
86 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
87 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
88 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
89       do i=1,ncont
90         do j=1,ncont_ref
91           if (icont(1,i).eq.icont_ref(1,j) .and. 
92      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
93         enddo
94       enddo
95 c     print *,' nmatch=',nmatch
96 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
97       contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
98       return
99       end
100 c----------------------------------------------------------------------------
101       subroutine hairpin(lprint,nharp,iharp)
102       implicit real*8 (a-h,o-z)
103       include 'DIMENSIONS'
104       include 'COMMON.IOUNITS'
105       include 'COMMON.CHAIN'
106       include 'COMMON.INTERACT'
107       include 'COMMON.FFIELD'
108       include 'COMMON.NAMES'
109       integer ncont,icont(2,maxcont)
110       integer nharp,iharp(4,maxres/3)
111       logical lprint,not_done
112       real*8 rcomp /6.0d0/ 
113       ncont=0
114       kkk=0
115 c     print *,'nnt=',nnt,' nct=',nct
116       do i=nnt,nct-3
117         do k=1,3
118           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
119         enddo
120         do j=i+2,nct-1
121           do k=1,3
122             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
123           enddo
124           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
125             ncont=ncont+1
126             icont(1,ncont)=i
127             icont(2,ncont)=j
128           endif
129         enddo
130       enddo
131       if (lprint) then
132         write (iout,'(a)') 'PP contact map:'
133         do i=1,ncont
134           i1=icont(1,i)
135           i2=icont(2,i)
136           it1=itype(i1)
137           it2=itype(i2)
138           write (iout,'(i3,2x,a,i4,2x,a,i4)') 
139      &     i,restyp(it1),i1,restyp(it2),i2 
140         enddo
141       endif
142 c finding hairpins
143       nharp=0
144       do i=1,ncont
145         i1=icont(1,i)
146         j1=icont(2,i)
147         if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
148 c          write (iout,*) "found turn at ",i1,j1
149           ii1=i1
150           jj1=j1
151           not_done=.true.
152           do while (not_done)
153             i1=i1-1
154             j1=j1+1
155             do j=1,ncont
156               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
157             enddo
158             not_done=.false.
159   10        continue
160 c            write (iout,*) i1,j1,not_done  
161           enddo
162           i1=i1+1
163           j1=j1-1
164           if (j1-i1.gt.4) then
165             nharp=nharp+1
166             iharp(1,nharp)=i1
167             iharp(2,nharp)=j1
168             iharp(3,nharp)=ii1
169             iharp(4,nharp)=jj1 
170 c            write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
171           endif
172         endif
173       enddo
174 c      do i=1,nharp
175 c            write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
176 c      enddo
177       if (lprint) then
178       write (iout,*) "Hairpins:"
179       do i=1,nharp
180         i1=iharp(1,i)
181         j1=iharp(2,i)
182         ii1=iharp(3,i)
183         jj1=iharp(4,i)
184         write (iout,*)
185         write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1)
186         write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
187 c        do k=jj1,j1,-1
188 c         write (iout,'(a,i3,$)') restyp(itype(k)),k
189 c        enddo
190       enddo
191       endif
192       return
193       end
194 c----------------------------------------------------------------------------
195