small fix in .gitingnore
[unres.git] / source / wham / src / conf_compar.F
1       subroutine conf_compar(jcon,lprn,print_class)
2       implicit real*8 (a-h,o-z)
3 #ifdef MPI
4       include "mpif.h"
5 #endif
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8       include 'DIMENSIONS.COMPAR'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.CONTROL'
11       include 'COMMON.IOUNITS'
12       include 'COMMON.COMPAR'
13       include 'COMMON.CHAIN' 
14       include 'COMMON.INTERACT'
15       include 'COMMON.VAR'
16       include 'COMMON.PEPTCONT'
17       include 'COMMON.CONTACTS1'
18       include 'COMMON.HEADER'
19       include 'COMMON.FREE'
20       include 'COMMON.ENERGIES'
21 #ifdef MPI
22       include 'COMMON.MPI'
23 #endif
24       integer ilen
25       external ilen
26       logical lprn,print_class
27       integer ncont_frag(mmaxfrag),
28      & icont_frag(2,maxcont,mmaxfrag),ncontsc,
29      & icontsc(1,maxcont),nsccont_frag(mmaxfrag),
30      & isccont_frag(2,maxcont,mmaxfrag)
31       integer isecstr(maxres)
32       integer itemp(maxfrag)
33       character*4 liczba
34       double precision Epot
35 c      print *,"Enter conf_compar",jcon
36       call angnorm12(rmsang)
37 c Level 1: check secondary and supersecondary structure
38       call elecont(lprn,ncont,icont,nnt,nct)
39       call secondary2(lprn,.false.,ncont,icont,isecstr)
40       call contact(lprn,ncontsc,icontsc,nnt,nct)
41       if (lprn) write(iout,*) "Assigning electrostatic contacts"
42       call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,
43      &   icont_frag)
44       if (lprn) write(iout,*) "Assigning sidechain contacts"
45       call contacts_between_fragments(lprn,3,ncontsc,icontsc,
46      &   nsccont_frag,isccont_frag)
47       do i=1,nlevel
48         do j=1,isnfrag(nlevel+1)
49           iclass(j,i)=0
50         enddo
51       enddo
52       do j=1,nfrag(1)
53         ind = icant(j,j)
54         if (lprn) then
55           write (iout,'(80(1h=))') 
56           write (iout,*) "Level",1," fragment",j
57           write (iout,'(80(1h=))') 
58         endif
59         rmsfrag(j,1)=rmscalc(0,1,j,jcon,lprn)
60 c Compare electrostatic contacts in the current conf with that in the native
61 c structure.
62         if (lprn) write (iout,*) 
63      &    "Comparing electrostatic contact map and local structure" 
64         ncnat=ncont_frag_ref(ind)
65 c        write (iout,*) "before match_contact:",nc_fragm(j,1),
66 c     &   nc_req_setf(j,1)
67         call match_secondary(j,isecstr,nsec_match,lprn)
68         if (lprn) write (iout,*) "Fragment",j," nsec_match",
69      &    nsec_match," length",len_frag(j,1)," min_len",
70      &    frac_sec*len_frag(j,1)
71         if (nsec_match.lt.frac_sec*len_frag(j,1)) then
72           iclass(j,1)=0
73           if (lprn) write (iout,*) "Fragment",j,
74      &      " has incorrect secondary structure"
75         else
76           iclass(j,1)=1
77           if (lprn) write (iout,*) "Fragment",j,
78      &      " has correct secondary structure"
79         endif
80         if (ielecont(j,1).gt.0) then
81           call match_contact(ishif1,ishif2,nc_match,ncon_match,
82      &      ncont_frag_ref(ind),icont_frag_ref(1,1,ind),
83      &      ncont_frag(ind),icont_frag(1,1,ind),
84      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
85      &      nc_req_setf(j,1),istruct(j),.true.,lprn)
86         else if (isccont(j,1).gt.0) then
87           call match_contact(ishif1,ishif2,nc_match,ncon_match,
88      &      nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),
89      &      nsccont_frag(ind),isccont_frag(1,1,ind),
90      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
91      &      nc_req_setf(j,1),istruct(j),.true.,lprn)
92         else if (iloc(j).gt.0) then
93 c          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
94           call match_contact(ishif1,ishif2,nc_match,ncon_match,
95      &      0,icont_frag_ref(1,1,ind),
96      &      ncont_frag(ind),icont_frag(1,1,ind),
97      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
98      &      0,istruct(j),.true.,lprn)
99 c          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
100         else
101           ishif=0
102           nc_match=1
103         endif
104         if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
105         ishif=ishif1
106         qfrag(j,1)=qwolynes(1,j)
107         if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
108         if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match
109 c        write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1)
110         if (irms(j,1).gt.0) then
111           if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then
112             iclass_rms=2
113             ishifft_rms=0
114           else
115             ishiff=0
116             rms=1.0d2
117             iclass_rms=0
118             do while (rms.gt.rmscutfrag(1,j,1) .and. 
119      &         ishiff.lt.n_shift(1,j,1))
120               ishiff=ishiff+1
121               rms=rmscalc(-ishiff,1,j,jcon,lprn)
122 c              write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
123 c     &          " rms",rms," rmscut",rmscutfrag(1,j,1)
124               if (lprn) write (iout,*) "rms",rmsfrag(j,1) 
125               if (rms.gt.rmscutfrag(1,j,1)) then
126                 rms=rmscalc(ishiff,1,j,jcon,lprn)
127 c                write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
128 c     &           " rms",rms
129               endif
130               if (lprn) write (iout,*) "rms",rmsfrag(j,1) 
131             enddo
132 c            write (iout,*) "After loop: rms",rms,
133 c     &        " rmscut",rmscutfrag(1,j,1)
134 c            write (iout,*) "iclass_rms",iclass_rms
135             if (rms.le.rmscutfrag(1,j,1)) then
136               ishifft_rms=ishiff
137               rmsfrag(j,1)=rms
138               iclass_rms=1
139             endif
140 c            write (iout,*) "iclass_rms",iclass_rms
141           endif
142 c          write (iout,*) "ishif",ishif
143           if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
144         else
145           iclass_rms=1
146         endif
147 c        write (iout,*) "ishif",ishif," iclass",iclass(j,1),
148 c     &    " iclass_rms",iclass_rms
149         if (nc_match.gt.0 .and. iclass_rms.gt.0) then
150           if (ishif.eq.0) then
151             iclass(j,1)=iclass(j,1)+6
152           else
153             iclass(j,1)=iclass(j,1)+2
154           endif
155         endif
156         ncont_nat(1,j,1)=nc_match
157         ncont_nat(2,j,1)=ncon_match
158         ishifft(j,1)=ishif
159 c        write (iout,*) "iclass",iclass(j,1)
160       enddo
161 c Next levels: Check arrangements of elementary fragments.
162       do i=2,nlevel
163         do j=1,nfrag(i)
164         if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i))
165         if (lprn) then
166             write (iout,'(80(1h=))') 
167             write (iout,*) "Level",i," fragment",j
168             write (iout,'(80(1h=))') 
169         endif
170 c If an elementary fragment doesn't exist, don't check higher hierarchy levels.
171         do k=1,npiece(j,i)
172           ik=ipiece(k,j,i)
173           if (iclass(ik,1).eq.0) then
174             iclass(j,i)=0
175             goto 12
176           endif
177         enddo
178         if (i.eq.2 .and. ielecont(j,i).gt.0) then
179           iclass_con=0
180           ishifft_con=0
181           if (lprn) write (iout,*) 
182      &     "Comparing electrostatic contact map: fragments",
183      &      ipiece(1,j,i),ipiece(2,j,i)," ind",ind
184           call match_contact(ishif1,ishif2,nc_match,ncon_match,
185      &     ncont_frag_ref(ind),icont_frag_ref(1,1,ind),
186      &     ncont_frag(ind),icont_frag(1,1,ind),
187      &     j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),
188      &     nc_req_setf(j,i),2,.false.,lprn)
189           ishif=ishif1
190           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
191           if (nc_match.gt.0) then
192             if (ishif.eq.0) then
193               iclass_con=2
194             else
195               iclass_con=1
196             endif
197           endif
198           ncont_nat(1,j,i)=nc_match
199           ncont_nat(2,j,i)=ncon_match
200           ishifft_con=ishif
201         else if (i.eq.2 .and. isccont(j,i).gt.0) then
202           iclass_con=0
203           ishifft_con=0
204           if (lprn) write (iout,*) 
205      &     "Comparing sidechain contact map: fragments",
206      &    ipiece(1,j,i),ipiece(2,j,i)," ind",ind
207           call match_contact(ishif1,ishif2,nc_match,ncon_match,
208      &     nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),
209      &     nsccont_frag(ind),isccont_frag(1,1,ind),
210      &     j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),
211      &     nc_req_setf(j,i),2,.false.,lprn)
212           ishif=ishif1
213           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
214           if (nc_match.gt.0) then
215             if (ishif.eq.0) then
216               iclass_con=2
217             else
218               iclass_con=1
219             endif
220           endif
221           ncont_nat(1,j,i)=nc_match
222           ncont_nat(2,j,i)=ncon_match
223           ishifft_con=ishif
224         else if (i.eq.2) then
225           iclass_con=2
226           ishifft_con=0
227         endif
228         if (i.eq.2) qfrag(j,2)=qwolynes(2,j)
229         if (lprn) write (iout,*) 
230      &    "Comparing rms: fragments",
231      &     (ipiece(k,j,i),k=1,npiece(j,i))
232         rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn)
233         if (irms(j,i).gt.0) then
234           iclass_rms=0
235           ishifft_rms=0
236           if (lprn) write (iout,*) "rms",rmsfrag(j,i)
237 c          write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i),
238 c     &     " rmscutfrag",rmscutfrag(1,j,i)
239           if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then
240             iclass_rms=2
241             ishifft_rms=0
242           else
243             ishif=0
244             rms=1.0d2
245             do while (rms.gt.rmscutfrag(1,j,i) .and. 
246      &         ishif.lt.n_shift(1,j,i))
247               ishif=ishif+1
248               rms=rmscalc(-ishif,i,j,jcon,lprn)
249 c              print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
250               if (lprn) write (iout,*) "rms",rmsfrag(j,i) 
251               if (rms.gt.rmscutfrag(1,j,i)) then
252                 rms=rmscalc(ishif,i,j,jcon,lprn)
253 c                print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
254               endif
255               if (lprn) write (iout,*) "rms",rms
256             enddo
257             if (rms.le.rmscutfrag(1,j,i)) then
258               ishifft_rms=ishif
259               rmsfrag(j,i)=rms
260               iclass_rms=1
261             endif
262           endif
263         endif
264         if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and. 
265      &    isccont(j,i).eq.0 ) then
266           write (iout,*) "Error: no measure of comparison specified:",
267      &      " level",i," part",j
268           stop
269         endif
270         if (lprn) 
271      &  write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
272         if (i.eq.2) then
273           iclass(j,i) = min0(iclass_con,iclass_rms)
274           if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then
275             ishifft(j,i)=ishifft_rms
276           else
277             ishifft(j,i)=ishifft_con
278           endif
279         else if (i.gt.2) then
280           iclass(j,i) = iclass_rms
281           ishifft(j,i)= ishifft_rms
282         endif
283    12   continue
284         enddo
285       enddo
286       rms_nat=rmsnat(jcon)
287       qnat=qwolynes(0,0)
288 C Compute the structural class
289       iscor=0
290       IF (.NOT. BINARY) THEN
291       do i=1,nlevel
292         IF (I.EQ.1) THEN
293         do j=1,nfrag(i)
294           itemp(j)=iclass(j,i)
295         enddo
296         do kk=-1,1
297           do j=1,nfrag(i)
298             idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j
299             iex = 2**idig
300             im=mod(itemp(j),2)
301             itemp(j)=itemp(j)/2
302 c            write (iout,*) "i",i," j",j," idig",idig," iex",iex,
303 c     &        " iclass",iclass(j,i)," im",im
304             iscor=iscor+im*iex
305           enddo
306         enddo
307         ELSE
308         do j=1,nfrag(i)
309           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j
310           iex = 2**idig
311           if (iclass(j,i).gt.0) then
312             im=1
313           else
314             im=0
315           endif
316 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
317 c     &      " iclass",iclass(j,i)," im",im
318           iscor=iscor+im*iex
319         enddo
320         do j=1,nfrag(i)
321           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j
322           iex = 2**idig
323           if (iclass(j,i).gt.1) then
324             im=1
325           else
326             im=0
327           endif
328 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
329 c     &      " iclass",iclass(j,i)," im",im
330           iscor=iscor+im*iex
331         enddo
332         ENDIF
333       enddo
334       iscore=iscor
335       ENDIF
336       if (print_class) then
337 #ifdef MPI
338           write(istat,'(i6,$)') jcon+indstart(me)-1
339           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),
340      &     -entfac(jcon)
341 #else
342           write(istat,'(i6,$)') jcon
343           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),
344      &      -entfac(jcon)
345 #endif
346           write (istat,'(f8.3,2f6.3,$)') 
347      &      rms_nat,qnat,rmsang/(nres-3)
348           do j=1,nlevel
349             write(istat,'(1x,$,20(i3,$))')
350      &        (ncont_nat(1,k,j),k=1,nfrag(j))
351             if (j.lt.3) then
352               write(istat,'(1x,$,20(f5.1,f5.2$))')
353      &          (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j))
354             else
355               write(istat,'(1x,$,20(f5.1$))')
356      &          (rmsfrag(k,j),k=1,nfrag(j))
357             endif
358             write(istat,'(1x,$,20(i1,$))')
359      &        (iclass(k,j),k=1,nfrag(j))
360           enddo
361           if (binary) then
362             write (istat,'("  ",$)')
363             do j=1,nlevel
364               write (istat,'(100(i1,$))')(iclass(k,j),
365      &           k=1,nfrag(j))
366               if (j.lt.nlevel) write(iout,'(".",$)')
367             enddo
368             write (istat,*)
369           else
370             write (istat,'(i10)') iscore
371           endif
372       endif
373       RETURN
374       END