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