1 subroutine add_angpair(ici,icj,nang_pair,iang_pair)
2 implicit real*8 (a-h,o-z)
4 include 'COMMON.IOUNITS'
6 integer ici,icj,nang_pair,iang_pair(2,maxres)
8 c write (iout,*) "add_angpair: ici",ici," icj",icj,
9 c & " nang_pair",nang_pair
11 if (ian1.lt.4 .or. ian1.gt.nres) return
13 c write (iout,*) "ian1",ian1," ian2",ian2
14 if (ian2.lt.4 .or. ian2.gt.nres) return
16 if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
19 iang_pair(1,nang_pair)=ian1
20 iang_pair(2,nang_pair)=ian2
23 c-------------------------------------------------------------------------
24 subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,
26 implicit real*8 (a-h,o-z)
28 include 'DIMENSIONS.ZSCOPT'
29 include 'DIMENSIONS.COMPAR'
30 include 'COMMON.IOUNITS'
32 include 'COMMON.COMPAR'
33 include 'COMMON.CHAIN'
35 double precision pinorm,deltang
38 write (iout,'(80(1h*))')
39 write (iout,*) "angnorm"
40 write (iout,'(80(1h*))')
45 npart = npiece(jfrag,1)
47 nne = min0(nend_sup,nres)
48 if (lprn) write (iout,*) "nn4",nn4," nne",nne
50 nbeg = ifrag(1,i,jfrag) + 3 - ishif1
51 if (nbeg.lt.nn4) nbeg=nn4
52 nend = ifrag(2,i,jfrag) + 1 - ishif2
53 if (nend.gt.nne) nend=nne
54 if (nend.ge.nbeg) then
55 nn = nn + nend - nbeg + 1
57 write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
58 & " nn",nn," ishift1",ishif1," ishift2",ishif2
59 write (iout,*) "angles"
65 c deltang = pinorm(phi(j)-phi_ref(j+ishif1))
66 deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),
67 & theta_ref(j+ishif1),phi(iperm(j,ipermmin)),
68 & theta(iperm(j-1,ipermmin)),theta(iperm(j,ipermmin)))
69 if (dabs(deltang).gt.diffang_max) then
70 if (ll.gt.longest) longest = ll
75 if (ll.gt.longest) longest = ll
76 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
77 & rad2deg*phi_ref(j+ishif1),rad2deg*deltang
78 angn=angn+dabs(deltang)
81 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
82 if (lprn) write (iout,*)"segment",i," longest fragment within",
83 & diffang_max*rad2deg,":",longest," fraction",ff
84 if (ff.lt.fract) fract = ff
92 if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,
96 c-------------------------------------------------------------------------
97 subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,
98 & diffang_max,anorm,fract,ipermmin)
99 implicit real*8 (a-h,o-z)
101 include 'DIMENSIONS.ZSCOPT'
102 include 'DIMENSIONS.COMPAR'
103 include 'COMMON.IOUNITS'
105 include 'COMMON.COMPAR'
106 include 'COMMON.CHAIN'
108 integer ncont,icont(2,ncont),longest
109 double precision anorm,diffang_max,fract
110 integer npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
111 double precision pinorm
112 integer iperm,ipermmin
114 if (lprn) write (iout,'(80(1h*))')
116 c Determine the segments for which angles will be compared
119 nne = min0(nend_sup,nres)
120 if (lprn) write (iout,*) "nn4",nn4," nne",nne
121 npart=npiece(jfrag,1)
124 c write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
125 if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or.
126 & icont(1,1).gt.ifrag(2,i,jfrag)) goto 11
128 do while (jstart.lt.ncont .and.
129 & icont(1,jstart).lt.ifrag(1,i,jfrag))
130 c write (iout,*) "jstart",jstart," icont",icont(1,jstart),
131 c & " ifrag",ifrag(1,i,jfrag)
134 c write (iout,*) "jstart",jstart," icont",icont(1,jstart),
135 c & " ifrag",ifrag(1,i,jfrag)
136 if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11
139 ifrag_c(1,npiece_c)=icont(1,jstart)
141 do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag))
142 c write (iout,*) "jend",jend," icont",icont(1,jend),
143 c & " ifrag",ifrag(2,i,jfrag)
146 c write (iout,*) "jend",jend," icont",icont(1,jend),
147 c & " ifrag",ifrag(2,i,jfrag)
149 ifrag_c(2,npiece_c)=icont(1,jend)+1
150 ishift_c(npiece_c)=ishif1
151 c write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
152 c & " ic1",ic1," ic2",ic2,
153 c & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
155 if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
160 c write (iout,*) "idi",idi
162 if (icont(2,1).gt.ifrag(2,i,jfrag) .or.
163 & icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12
165 do while (jstart.lt.ncont .and.
166 & icont(2,jstart).lt.ifrag(1,i,jfrag))
167 c write (iout,*) "jstart",jstart," icont",icont(2,jstart),
168 c & " ifrag",ifrag(1,i,jfrag)
171 c write (iout,*) "jstart",jstart," icont",icont(2,jstart),
172 c & " ifrag",ifrag(1,i,jfrag)
173 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
176 ifrag_c(2,npiece_c)=icont(2,jstart)+1
178 do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag))
179 c write (iout,*) "jend",jend," icont",icont(2,jend),
180 c & " ifrag",ifrag(2,i,jfrag)
183 c write (iout,*) "jend",jend," icont",icont(2,jend),
184 c & " ifrag",ifrag(2,i,jfrag)
185 else if (idi.eq.-1) then
186 if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or.
187 & icont(2,1).lt.ifrag(1,i,jfrag)) goto 12
189 do while (jstart.gt.ncont .and.
190 & icont(2,jstart).lt.ifrag(1,i,jfrag))
191 c write (iout,*) "jstart",jstart," icont",icont(2,jstart),
192 c & " ifrag",ifrag(1,i,jfrag)
195 c write (iout,*) "jstart",jstart," icont",icont(2,jstart),
196 c & " ifrag",ifrag(1,i,jfrag)
197 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
200 ifrag_c(2,npiece_c)=icont(2,jstart)+1
202 do while (jend.lt.ncont .and.
203 & icont(2,jend).gt.ifrag(2,i,jfrag))
204 c write (iout,*) "jend",jend," icont",icont(2,jend),
205 c & " ifrag",ifrag(2,i,jfrag)
208 c write (iout,*) "jend",jend," icont",icont(2,jend),
209 c & " ifrag",ifrag(2,i,jfrag)
217 c write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
218 c & " jstart:",jstart," jend",jend,
219 c & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
220 ifrag_c(1,npiece_c)=ic1
221 ifrag_c(2,npiece_c)=ic2+1
222 ishift_c(npiece_c)=ishif2
226 write (iout,*) "Before merge: npiece_c",npiece_c
228 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
232 c Merge overlapping segments (e.g., avoid splitting helices)
235 do while (i .lt. npiece_c)
236 if (ishift_c(i).eq.ishift_c(i+1) .and.
237 & ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
238 ifrag_c(2,i)=ifrag_c(2,i+1)
240 ishift_c(j)=ishift_c(j+1)
241 ifrag_c(1,j)=ifrag_c(1,j+1)
242 ifrag_c(2,j)=ifrag_c(2,j+1)
250 write (iout,*) "After merge: npiece_c",npiece_c
252 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
265 nbeg = ifrag_c(1,i) + 3 - ishifc
266 if (nbeg.lt.nn4) nbeg=nn4
267 nend = ifrag_c(2,i) - ishifc + 1
268 if (nend.gt.nne) nend=nne
269 if (nend.ge.nbeg) then
270 nn = nn + nend - nbeg + 1
271 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,
272 & " nn",nn," ishifc",ishifc
273 if (lprn) write (iout,*) "angles"
277 c deltang = pinorm(phi(j)-phi_ref(j+ishifc))
278 deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),
279 & theta_ref(j+ishifc),phi(iperm(j,ipermmin)),
280 & theta(iperm(j-1,ipermmin)),theta(iperm(j,ipermmin)))
281 if (dabs(deltang).gt.diffang_max) then
282 if (ll.gt.longest) longest = ll
287 if (ll.gt.longest) longest = ll
288 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),
289 & rad2deg*phi_ref(j+ishifc),rad2deg*deltang
290 angn=angn+dabs(deltang)
293 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
294 if (lprn) write (iout,*)"segment",i," longest fragment within",
295 & diffang_max*rad2deg,":",longest," fraction",ff
296 if (ff.lt.fract) fract = ff
299 if (nn.gt.0) anorm = angn/nn
300 if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
303 c-------------------------------------------------------------------------
304 double precision function angnorm1(nang_pair,iang_pair,ipermmin,
306 implicit real*8 (a-h,o-z)
308 include 'DIMENSIONS.ZSCOPT'
309 include 'DIMENSIONS.COMPAR'
310 include 'COMMON.IOUNITS'
312 include 'COMMON.COMPAR'
313 include 'COMMON.CHAIN'
316 integer nang_pair,iang_pair(2,maxres)
317 double precision pinorm
318 integer iperm,ipermmin
320 if (lprn) write (iout,'(80(1h*))')
321 if (lprn) write (iout,*) "nang_pair",nang_pair
322 if (lprn) write (iout,*) "angles"
326 c deltang = pinorm(phi(ia1)-phi_ref(ia2))
327 deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),
328 & theta_ref(ia2),phi(iperm(ia2,ipemmin)),
329 & theta(iperm(ia2-1,ipermmin)),theta(iperm(ia2-1,ipermmin)))
330 if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),
331 & rad2deg*phi_ref(ia2),rad2deg*deltang
332 angn=angn+dabs(deltang)
335 &write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
336 angnorm1 = angn/nang_pair
339 c------------------------------------------------------------------------------
340 subroutine angnorm12(diff,ipermmin)
341 implicit real*8 (a-h,o-z)
343 include 'DIMENSIONS.ZSCOPT'
344 include 'DIMENSIONS.COMPAR'
345 include 'COMMON.IOUNITS'
347 include 'COMMON.COMPAR'
348 include 'COMMON.CHAIN'
350 double precision pinorm
351 integer ipermm,ipermmin
354 nne = min0(nend_sup,nres)
356 c diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j)))
359 c diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j)))
360 diff=diff+spherang(phi_ref(j),theta_ref(j-1),
361 & theta_ref(j),phi(iperm(j,ipermmin)),theta(iperm(j-1,ipermmin)),
362 & theta(iperm(j,ipermmin)))
366 c--------------------------------------------------------------------------------
367 double precision function spherang(gam1,theta11,theta12,
368 & gam2,theta21,theta22)
370 double precision gam1,theta11,theta12,gam2,theta21,theta22,
371 & x1,x2,xmed,f1,f2,fmed
372 double precision tolx /1.0d-4/, tolf /1.0d-4/
373 double precision sumcos
374 double precision arcos,pinorm,sumangp
375 integer it,maxit /100/
376 c Calculate the difference of the angles of two superposed 4-redidue fragments
384 c The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
385 c to achieve the minimum difference between the O'-C-O and P-C-P angles;
386 c the sum of these angles is the difference returned by the function.
389 c If thetas match, take the difference of gamma and exit.
390 if (dabs(theta11-theta12).lt.tolx
391 & .and. dabs(theta21-theta22).lt.tolx) then
392 spherang=dabs(pinorm(gam2-gam1))
395 c If the gammas are the same, take the difference of thetas and exit.
397 x2=0.5d0*pinorm(gam2-gam1)
398 if (dabs(x2) .lt. tolx) then
399 spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
401 else if (x2.lt.0.0d0) then
405 c Else apply regula falsi method to compute optimum overlap of the terminal Calphas
406 f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
407 f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
409 xmed=x1-f1*(x2-x1)/(f2-f1)
410 fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
411 c write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
412 if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx)
413 & .and. dabs(fmed).lt.tolf ) then
417 else if ( fmed*f1.lt.0.0d0 ) then
426 spherang=arcos(dcos(theta11)*dcos(theta12)
427 & +dsin(theta11)*dsin(theta12)*dcos(x1))+
428 & arcos(dcos(theta21)*dcos(theta22)+
429 & dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
432 c--------------------------------------------------------------------------------
433 double precision function sumangp(gam1,theta11,theta12,gam2,
434 & theta21,theta22,fi)
436 double precision gam1,theta11,theta12,gam2,theta21,theta22,fi,
437 & cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,
439 c derivarive of the sum of the difference of the angles of a 4-residue fragment.
440 double precision arcos
449 cosd1=cost11*cost12+sint11*sint12*dcos(fi)
450 cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
451 sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1)
452 & +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)