saxs and adam's corrections to multichain
[unres.git] / source / wham / src / angnorm.f
1       subroutine add_angpair(ici,icj,nang_pair,iang_pair)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.CHAIN'
6       integer ici,icj,nang_pair,iang_pair(2,maxres)
7       integer i,ian1,ian2
8 c      write (iout,*) "add_angpair: ici",ici," icj",icj,
9 c     &  " nang_pair",nang_pair
10       ian1=ici+2
11       if (ian1.lt.4 .or. ian1.gt.nres) return
12       ian2=icj+2
13 c      write (iout,*) "ian1",ian1," ian2",ian2
14       if (ian2.lt.4 .or. ian2.gt.nres) return
15       do i=1,nang_pair
16         if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
17       enddo
18       nang_pair=nang_pair+1
19       iang_pair(1,nang_pair)=ian1
20       iang_pair(2,nang_pair)=ian2
21       return
22       end
23 c-------------------------------------------------------------------------
24       subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,
25      &  lprn)
26       implicit real*8 (a-h,o-z)
27       include 'DIMENSIONS'
28       include 'DIMENSIONS.ZSCOPT'
29       include 'DIMENSIONS.COMPAR'
30       include 'COMMON.IOUNITS'
31       include 'COMMON.VAR'
32       include 'COMMON.COMPAR'
33       include 'COMMON.CHAIN'
34       include 'COMMON.GEO'
35       double precision pinorm,deltang
36       logical lprn
37       if (lprn) write (iout,'(80(1h*))')
38       angn=0.0d0
39       nn = 0
40       fract = 1.0d0
41       npart = npiece(jfrag,1)
42       nn4 = nstart_sup+3
43       nne = min0(nend_sup,nres)
44       if (lprn) write (iout,*) "nn4",nn4," nne",nne
45       do i=1,npart
46         nbeg = ifrag(1,i,jfrag) + 3 - ishif1
47         if (nbeg.lt.nn4) nbeg=nn4
48         nend = ifrag(2,i,jfrag) + 1 - ishif2
49         if (nend.gt.nne) nend=nne
50         if (nend.ge.nbeg) then
51         nn = nn + nend - nbeg + 1
52         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
53      &    " nn",nn," ishift1",ishif1," ishift2",ishif2
54         if (lprn) write (iout,*) "angles"
55         longest=0
56         ll = 0
57         do j=nbeg,nend
58 c          deltang = pinorm(phi(j)-phi_ref(j+ishif1))
59           deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),
60      &      theta_ref(j+ishif1),phi(j),theta(j-1),theta(j))
61           if (dabs(deltang).gt.diffang_max) then
62             if (ll.gt.longest) longest = ll
63             ll = 0
64           else
65             ll=ll+1
66           endif
67           if (ll.gt.longest) longest = ll
68           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
69      &     rad2deg*phi_ref(j+ishif1),rad2deg*deltang
70           angn=angn+dabs(deltang)
71         enddo
72         longest=longest+3
73         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
74         if (lprn) write (iout,*)"segment",i," longest fragment within",
75      &    diffang_max*rad2deg,":",longest," fraction",ff
76         if (ff.lt.fract) fract = ff
77         endif
78       enddo
79       if (nn.gt.0) then
80         angn = angn/nn
81       else
82         angn = dwapi
83       endif
84       if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,
85      &  " fract",fract
86       return
87       end
88 c-------------------------------------------------------------------------
89       subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,
90      &  diffang_max,anorm,fract)
91       implicit real*8 (a-h,o-z)
92       include 'DIMENSIONS'
93       include 'DIMENSIONS.ZSCOPT'
94       include 'DIMENSIONS.COMPAR'
95       include 'COMMON.IOUNITS'
96       include 'COMMON.VAR'
97       include 'COMMON.COMPAR'
98       include 'COMMON.CHAIN'
99       include 'COMMON.GEO'
100       integer ncont,icont(2,ncont),longest
101       double precision anorm,diffang_max,fract
102       integer npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
103       double precision pinorm
104       logical lprn
105       if (lprn) write (iout,'(80(1h*))')
106 c
107 c Determine the segments for which angles will be compared
108 c
109       nn4 = nstart_sup+3
110       nne = min0(nend_sup,nres)
111       if (lprn) write (iout,*) "nn4",nn4," nne",nne
112       npart=npiece(jfrag,1)
113       npiece_c=0
114       do i=1,npart
115 c        write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
116         if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. 
117      &    icont(1,1).gt.ifrag(2,i,jfrag)) goto 11
118         jstart=1
119         do while (jstart.lt.ncont .and. 
120      &   icont(1,jstart).lt.ifrag(1,i,jfrag))
121 c          write (iout,*) "jstart",jstart," icont",icont(1,jstart),
122 c     &     " ifrag",ifrag(1,i,jfrag)
123           jstart=jstart+1
124         enddo
125 c        write (iout,*) "jstart",jstart," icont",icont(1,jstart),
126 c     &   " ifrag",ifrag(1,i,jfrag)
127         if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11
128         npiece_c=npiece_c+1
129         ic1=icont(1,jstart)
130         ifrag_c(1,npiece_c)=icont(1,jstart)
131         jend=ncont
132         do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag))
133 c          write (iout,*) "jend",jend," icont",icont(1,jend),
134 c     &     " ifrag",ifrag(2,i,jfrag)
135           jend=jend-1
136         enddo
137 c        write (iout,*) "jend",jend," icont",icont(1,jend),
138 c     &   " ifrag",ifrag(2,i,jfrag)
139         ic2=icont(1,jend)
140         ifrag_c(2,npiece_c)=icont(1,jend)+1
141         ishift_c(npiece_c)=ishif1
142 c        write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
143 c     &    " ic1",ic1," ic2",ic2,
144 c     &    " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
145    11   continue
146         if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
147           idi=1
148         else
149           idi=-1
150         endif
151 c        write (iout,*) "idi",idi
152         if (idi.eq.1) then
153           if (icont(2,1).gt.ifrag(2,i,jfrag) .or. 
154      &      icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12
155           jstart=1
156           do while (jstart.lt.ncont .and. 
157      &     icont(2,jstart).lt.ifrag(1,i,jfrag))
158 c           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
159 c     &     " ifrag",ifrag(1,i,jfrag)
160             jstart=jstart+1
161           enddo
162 c          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
163 c     &     " ifrag",ifrag(1,i,jfrag)
164           if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
165           npiece_c=npiece_c+1
166           ic1=icont(2,jstart)
167           ifrag_c(2,npiece_c)=icont(2,jstart)+1
168           jend=ncont
169           do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag))
170 c            write (iout,*) "jend",jend," icont",icont(2,jend),
171 c     &     " ifrag",ifrag(2,i,jfrag)
172             jend=jend-1
173           enddo
174 c          write (iout,*) "jend",jend," icont",icont(2,jend),
175 c     &     " ifrag",ifrag(2,i,jfrag)
176         else if (idi.eq.-1) then
177           if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or.
178      &        icont(2,1).lt.ifrag(1,i,jfrag)) goto 12
179           jstart=ncont
180           do while (jstart.gt.ncont .and. 
181      &     icont(2,jstart).lt.ifrag(1,i,jfrag))
182 c           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
183 c     &     " ifrag",ifrag(1,i,jfrag)
184             jstart=jstart-1
185           enddo
186 c          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
187 c     &     " ifrag",ifrag(1,i,jfrag)
188           if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
189           npiece_c=npiece_c+1
190           ic1=icont(2,jstart)
191           ifrag_c(2,npiece_c)=icont(2,jstart)+1
192           jend=1
193           do while (jend.lt.ncont .and. 
194      &       icont(2,jend).gt.ifrag(2,i,jfrag))
195 c             write (iout,*) "jend",jend," icont",icont(2,jend),
196 c     &         " ifrag",ifrag(2,i,jfrag)
197             jend=jend+1
198           enddo
199 c          write (iout,*) "jend",jend," icont",icont(2,jend),
200 c     &     " ifrag",ifrag(2,i,jfrag)
201         endif
202         ic2=icont(2,jend)
203         if (ic2.lt.ic1) then
204           iic = ic1
205           ic1 = ic2
206           ic2 = iic
207         endif
208 c        write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
209 c     &    " jstart:",jstart," jend",jend,
210 c     &    " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
211         ifrag_c(1,npiece_c)=ic1
212         ifrag_c(2,npiece_c)=ic2+1
213         ishift_c(npiece_c)=ishif2
214    12   continue
215       enddo
216       if (lprn) then
217         write (iout,*) "Before merge: npiece_c",npiece_c
218         do i=1,npiece_c
219           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
220         enddo
221       endif
222 c
223 c Merge overlapping segments (e.g., avoid splitting helices)
224 c
225       i=1
226       do while (i .lt. npiece_c)
227         if (ishift_c(i).eq.ishift_c(i+1) .and. 
228      &     ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
229            ifrag_c(2,i)=ifrag_c(2,i+1)
230            do j=i+1,npiece_c
231              ishift_c(j)=ishift_c(j+1)
232              ifrag_c(1,j)=ifrag_c(1,j+1)
233              ifrag_c(2,j)=ifrag_c(2,j+1)
234            enddo
235            npiece_c=npiece_c-1
236         else
237           i=i+1
238         endif
239       enddo
240       if (lprn) then
241         write (iout,*) "After merge: npiece_c",npiece_c
242         do i=1,npiece_c
243           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
244         enddo
245       endif
246 c
247 c Compare angles
248 c
249       angn=0.0d0
250       anorm=0
251       nn = 0
252       fract = 1.0d0
253       npart = npiece_c
254       do i=1,npart
255         ishifc=ishift_c(i)
256         nbeg = ifrag_c(1,i) + 3 - ishifc
257         if (nbeg.lt.nn4) nbeg=nn4
258         nend = ifrag_c(2,i)  - ishifc + 1
259         if (nend.gt.nne) nend=nne
260         if (nend.ge.nbeg) then
261         nn = nn + nend - nbeg + 1
262         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
263      &    " nn",nn," ishifc",ishifc
264         if (lprn) write (iout,*) "angles"
265         longest=0
266         ll = 0
267         do j=nbeg,nend
268 c          deltang = pinorm(phi(j)-phi_ref(j+ishifc))
269           deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),
270      &      theta_ref(j+ishifc),phi(j),theta(j-1),theta(j))
271           if (dabs(deltang).gt.diffang_max) then
272             if (ll.gt.longest) longest = ll
273             ll = 0
274           else
275             ll=ll+1
276           endif
277           if (ll.gt.longest) longest = ll
278           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
279      &     rad2deg*phi_ref(j+ishifc),rad2deg*deltang
280           angn=angn+dabs(deltang)
281         enddo
282         longest=longest+3
283         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
284         if (lprn) write (iout,*)"segment",i," longest fragment within",
285      &    diffang_max*rad2deg,":",longest," fraction",ff
286         if (ff.lt.fract) fract = ff
287         endif
288       enddo
289       if (nn.gt.0) anorm = angn/nn
290       if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
291       return
292       end
293 c-------------------------------------------------------------------------
294       double precision function angnorm1(nang_pair,iang_pair,lprn)
295       implicit real*8 (a-h,o-z)
296       include 'DIMENSIONS'
297       include 'DIMENSIONS.ZSCOPT'
298       include 'DIMENSIONS.COMPAR'
299       include 'COMMON.IOUNITS'
300       include 'COMMON.VAR'
301       include 'COMMON.COMPAR'
302       include 'COMMON.CHAIN'
303       include 'COMMON.GEO'
304       logical lprn
305       integer nang_pair,iang_pair(2,maxres)
306       double precision pinorm
307       angn=0.0d0
308       if (lprn) write (iout,'(80(1h*))')
309       if (lprn) write (iout,*) "nang_pair",nang_pair
310       if (lprn) write (iout,*) "angles"
311       do j=1,nang_pair
312         ia1 = iang_pair(1,j)
313         ia2 = iang_pair(2,j)
314 c        deltang = pinorm(phi(ia1)-phi_ref(ia2))
315          deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),
316      &      theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2))
317         if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),
318      &   rad2deg*phi_ref(ia2),rad2deg*deltang
319         angn=angn+dabs(deltang)
320       enddo
321       if (lprn)
322      &write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
323       angnorm1 = angn/nang_pair
324       return
325       end
326 c------------------------------------------------------------------------------
327       subroutine angnorm12(diff)
328       implicit real*8 (a-h,o-z)
329       include 'DIMENSIONS'
330       include 'DIMENSIONS.ZSCOPT'
331       include 'DIMENSIONS.COMPAR'
332       include 'COMMON.IOUNITS'
333       include 'COMMON.VAR'
334       include 'COMMON.COMPAR'
335       include 'COMMON.CHAIN'
336       include 'COMMON.GEO'
337       double precision pinorm
338       diff=0.0d0
339       nn4 = nstart_sup+3
340       nne = min0(nend_sup,nres)
341 c      do j=nn4-1,nne
342 c        diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j)))
343 c      enddo
344       do j=nn4,nne 
345 c        diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j)))
346          diff=diff+spherang(phi_ref(j),theta_ref(j-1),
347      &      theta_ref(j),phi(j),theta(j-1),theta(j))
348       enddo
349       return
350       end
351 c--------------------------------------------------------------------------------
352       double precision function spherang(gam1,theta11,theta12,
353      &   gam2,theta21,theta22)
354       implicit none
355       double precision gam1,theta11,theta12,gam2,theta21,theta22,
356      &  x1,x2,xmed,f1,f2,fmed
357       double precision tolx /1.0d-4/, tolf /1.0d-4/
358       double precision sumcos
359       double precision arcos,pinorm,sumangp
360       integer it,maxit /100/
361 c Calculate the difference of the angles of two superposed 4-redidue fragments
362 c
363 c       O      P
364 c        \    /
365 c     O'--C--C       
366 c             \
367 c              P'
368 c
369 c The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
370 c to achieve the minimum difference between the O'-C-O and P-C-P angles;
371 c the sum of these angles is the difference returned by the function.
372 c
373 c 4/28/04 AL
374 c If thetas match, take the difference of gamma and exit.
375       if (dabs(theta11-theta12).lt.tolx 
376      & .and. dabs(theta21-theta22).lt.tolx) then
377         spherang=dabs(pinorm(gam2-gam1))
378         return
379       endif
380 c If the gammas are the same, take the difference of thetas and exit.
381       x1=0.0d0
382       x2=0.5d0*pinorm(gam2-gam1)
383       if (dabs(x2) .lt. tolx) then
384         spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
385         return
386       else if (x2.lt.0.0d0) then
387         x1=x2
388         x2=0.0d0
389       endif 
390 c Else apply regula falsi method to compute optimum overlap of the terminal Calphas
391       f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
392       f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
393       do it=1,maxit
394         xmed=x1-f1*(x2-x1)/(f2-f1)
395         fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
396 c        write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
397         if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx)
398      &       .and. dabs(fmed).lt.tolf ) then
399           x1=xmed
400           f1=fmed
401           goto 10
402         else if ( fmed*f1.lt.0.0d0 ) then
403           x2=xmed
404           f2=fmed
405         else
406           x1=xmed
407           f1=fmed
408         endif
409       enddo
410    10 continue
411       spherang=arcos(dcos(theta11)*dcos(theta12)
412      & +dsin(theta11)*dsin(theta12)*dcos(x1))+
413      & arcos(dcos(theta21)*dcos(theta22)+
414      & dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
415       return
416       end
417 c--------------------------------------------------------------------------------
418       double precision function sumangp(gam1,theta11,theta12,gam2,
419      & theta21,theta22,fi)
420       implicit none
421       double precision gam1,theta11,theta12,gam2,theta21,theta22,fi,
422      & cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,
423      & cosd2
424 c derivarive of the sum of the difference of the angles of a 4-residue fragment.
425       double precision arcos
426       cost11=dcos(theta11)
427       cost12=dcos(theta12)
428       cost21=dcos(theta21)
429       cost22=dcos(theta22)
430       sint11=dsin(theta11)
431       sint12=dsin(theta12)
432       sint21=dsin(theta21)
433       sint22=dsin(theta22)
434       cosd1=cost11*cost12+sint11*sint12*dcos(fi)
435       cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
436       sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1)
437      & +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)
438       return
439       end