2 !-----------------------------------------------------------------------------
5 use geometry_data, only:nres
7 use geometry, only:dist
8 use regularize_, only:fitsq
16 ! include "COMMON.MPI"
19 !-----------------------------------------------------------------------------
22 !-----------------------------------------------------------------------------
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine conf_compar(jcon)
29 ! implicit real*8 (a-h,o-z)
30 use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
31 ! nsccont_frag_ref,isccont_frag_ref
35 ! include 'DIMENSIONS'
36 ! include 'DIMENSIONS.ZSCOPT'
37 ! include 'DIMENSIONS.COMPAR'
38 ! include 'DIMENSIONS.FREE'
39 ! include 'COMMON.CONTROL'
40 ! include 'COMMON.IOUNITS'
41 ! include 'COMMON.COMPAR'
42 ! include 'COMMON.CHAIN'
43 ! include 'COMMON.INTERACT'
44 ! include 'COMMON.VAR'
45 ! include 'COMMON.PEPTCONT'
46 ! include 'COMMON.CONTACTS1'
47 ! include 'COMMON.HEADER'
48 ! include 'COMMON.FREE'
49 ! include 'COMMON.ENERGIES'
51 ! include 'COMMON.MPI'
55 logical :: lprn,print_class
56 integer :: ncont_frag(mmaxfrag),&
57 icont_frag(2,maxcont,mmaxfrag),ncontsc,&
58 icontsc(1,maxcont),nsccont_frag(mmaxfrag),&
59 isccont_frag(2,maxcont,mmaxfrag)
60 integer :: isecstr(nres)
61 integer :: itemp(maxfrag)
62 character(len=4) :: liczba
63 real(kind=8) :: Epot,rms
64 integer :: jcon,i,j,ind,ncnat,nsec_match,ishift,ishif1,ishif2,&
65 nc_match,ncon_match,iclass_rms,ishifft_rms,ishiff,ishif
66 integer :: k,kk,iclass_con,iscor,ik,ishifft_con,idig,iex,im
67 ! print *,"Enter conf_compar",jcon
70 write (iout,*) "before anything"
73 ! call angnorm12(rmsang)
74 ! Level 1: check secondary and supersecondary structure
75 ! call elecont(lprn,ncont,icont,nnt,nct)
77 write (iout,*) "elecont finished"
80 call secondary2(lprn,.false.,ncont,icont,isecstr)
82 write (iout,*) "secondary2 finished"
85 call contact(lprn,ncontsc,icontsc,nnt,nct)
87 write(iout,*) "Assigning electrostatic contacts"
90 call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,&
93 write(iout,*) "Assigning sidechain contacts"
96 call contacts_between_fragments(lprn,3,ncontsc,icontsc,&
97 nsccont_frag,isccont_frag)
99 write(iout,*) "--> After contacts_between_fragments"
103 do j=1,isnfrag(nlevel+1)
110 write (iout,'(80(1h=))')
111 write (iout,*) "Level",1," fragment",j
112 write (iout,'(80(1h=))')
115 rmsfrag(j,1)=rmscalc(0,1,j,jcon,lprn)
116 ! Compare electrostatic contacts in the current conf with that in the native
118 if (lprn) write (iout,*) &
119 "Comparing electrostatic contact map and local structure"
121 ncnat=ncont_frag_ref(ind)
122 ! write (iout,*) "before match_contact:",nc_fragm(j,1),
125 call match_secondary(j,isecstr,nsec_match,lprn)
126 if (lprn) write (iout,*) "Fragment",j," nsec_match",&
127 nsec_match," length",len_frag(j,1)," min_len",&
128 frac_sec*len_frag(j,1)
129 if (nsec_match.lt.frac_sec*len_frag(j,1)) then
131 if (lprn) write (iout,*) "Fragment",j,&
132 " has incorrect secondary structure"
135 if (lprn) write (iout,*) "Fragment",j,&
136 " has correct secondary structure"
138 if (ielecont(j,1).gt.0) then
139 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
140 ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
141 ncont_frag(ind),icont_frag(1,1,ind),&
142 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
143 nc_req_setf(j,1),istruct(j),.true.,lprn)
144 else if (isccont(j,1).gt.0) then
145 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
146 nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
147 nsccont_frag(ind),isccont_frag(1,1,ind),&
148 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
149 nc_req_setf(j,1),istruct(j),.true.,lprn)
150 else if (iloc(j).gt.0) then
151 ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
152 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
153 0,icont_frag_ref(1,1,ind),&
154 ncont_frag(ind),icont_frag(1,1,ind),&
155 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
156 0,istruct(j),.true.,lprn)
157 ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
162 if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
164 qfrag(j,1)=qwolynes(1,j)
165 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
166 if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match
167 ! write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1)
168 if (irms(j,1).gt.0) then
169 if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then
176 do while (rms.gt.rmscutfrag(1,j,1) .and. &
177 ishiff.lt.n_shift(1,j,1))
179 rms=rmscalc(-ishiff,1,j,jcon,lprn)
180 ! write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
181 ! & " rms",rms," rmscut",rmscutfrag(1,j,1)
182 if (lprn) write (iout,*) "rms",rmsfrag(j,1)
183 if (rms.gt.rmscutfrag(1,j,1)) then
184 rms=rmscalc(ishiff,1,j,jcon,lprn)
185 ! write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
188 if (lprn) write (iout,*) "rms",rmsfrag(j,1)
190 ! write (iout,*) "After loop: rms",rms,
191 ! & " rmscut",rmscutfrag(1,j,1)
192 ! write (iout,*) "iclass_rms",iclass_rms
193 if (rms.le.rmscutfrag(1,j,1)) then
198 ! write (iout,*) "iclass_rms",iclass_rms
200 ! write (iout,*) "ishif",ishif
201 if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
205 ! write (iout,*) "ishif",ishif," iclass",iclass(j,1),
206 ! & " iclass_rms",iclass_rms
207 if (nc_match.gt.0 .and. iclass_rms.gt.0) then
209 iclass(j,1)=iclass(j,1)+6
211 iclass(j,1)=iclass(j,1)+2
214 ncont_nat(1,j,1)=nc_match
215 ncont_nat(2,j,1)=ncon_match
217 ! write (iout,*) "iclass",iclass(j,1)
219 ! Next levels: Check arrangements of elementary fragments.
222 if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i))
224 write (iout,'(80(1h=))')
225 write (iout,*) "Level",i," fragment",j
226 write (iout,'(80(1h=))')
228 ! If an elementary fragment doesn't exist, don't check higher hierarchy levels.
231 if (iclass(ik,1).eq.0) then
236 if (i.eq.2 .and. ielecont(j,i).gt.0) then
239 if (lprn) write (iout,*) &
240 "Comparing electrostatic contact map: fragments",&
241 ipiece(1,j,i),ipiece(2,j,i)," ind",ind
242 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
243 ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
244 ncont_frag(ind),icont_frag(1,1,ind),&
245 j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
246 nc_req_setf(j,i),2,.false.,lprn)
248 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
249 if (nc_match.gt.0) then
256 ncont_nat(1,j,i)=nc_match
257 ncont_nat(2,j,i)=ncon_match
259 else if (i.eq.2 .and. isccont(j,i).gt.0) then
262 if (lprn) write (iout,*) &
263 "Comparing sidechain contact map: fragments",&
264 ipiece(1,j,i),ipiece(2,j,i)," ind",ind
265 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
266 nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
267 nsccont_frag(ind),isccont_frag(1,1,ind),&
268 j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
269 nc_req_setf(j,i),2,.false.,lprn)
271 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
272 if (nc_match.gt.0) then
279 ncont_nat(1,j,i)=nc_match
280 ncont_nat(2,j,i)=ncon_match
282 else if (i.eq.2) then
286 if (i.eq.2) qfrag(j,2)=qwolynes(2,j)
287 if (lprn) write (iout,*) &
288 "Comparing rms: fragments",&
289 (ipiece(k,j,i),k=1,npiece(j,i))
290 rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn)
291 if (irms(j,i).gt.0) then
294 if (lprn) write (iout,*) "rms",rmsfrag(j,i)
295 ! write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i),
296 ! & " rmscutfrag",rmscutfrag(1,j,i)
297 if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then
303 do while (rms.gt.rmscutfrag(1,j,i) .and. &
304 ishif.lt.n_shift(1,j,i))
306 rms=rmscalc(-ishif,i,j,jcon,lprn)
307 ! print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
308 if (lprn) write (iout,*) "rms",rmsfrag(j,i)
309 if (rms.gt.rmscutfrag(1,j,i)) then
310 rms=rmscalc(ishif,i,j,jcon,lprn)
311 ! print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
313 if (lprn) write (iout,*) "rms",rms
315 if (rms.le.rmscutfrag(1,j,i)) then
322 if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and. &
323 isccont(j,i).eq.0 ) then
324 write (iout,*) "Error: no measure of comparison specified:",&
329 write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
331 iclass(j,i) = min0(iclass_con,iclass_rms)
332 if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then
333 ishifft(j,i)=ishifft_rms
335 ishifft(j,i)=ishifft_con
337 else if (i.gt.2) then
338 iclass(j,i) = iclass_rms
339 ishifft(j,i)= ishifft_rms
346 ! Compute the structural class
348 IF (.NOT. BINARY) THEN
356 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j
360 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
361 ! & " iclass",iclass(j,i)," im",im
367 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j
369 if (iclass(j,i).gt.0) then
374 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
375 ! & " iclass",iclass(j,i)," im",im
379 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j
381 if (iclass(j,i).gt.1) then
386 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
387 ! & " iclass",iclass(j,i)," im",im
394 if (print_class) then
396 write(istat,'(i6,$)') jcon+indstart(me)-1
397 write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
400 write(istat,'(i6,$)') jcon
401 write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
404 write (istat,'(f8.3,2f6.3,$)') &
405 rms_nat,qnat,rmsang/(nres-3)
407 write(istat,'(1x,$,20(i3,$))') &
408 (ncont_nat(1,k,j),k=1,nfrag(j))
410 write(istat,'(1x,$,20(f5.1,f5.2$))') &
411 (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j))
413 write(istat,'(1x,$,20(f5.1$))') &
414 (rmsfrag(k,j),k=1,nfrag(j))
416 write(istat,'(1x,$,20(i1,$))') &
417 (iclass(k,j),k=1,nfrag(j))
420 write (istat,'(" ",$)')
422 write (istat,'(100(i1,$))')(iclass(k,j),&
424 if (j.lt.nlevel) write(iout,'(".",$)')
428 write (istat,'(i10)') iscore
432 END subroutine conf_compar
433 !-----------------------------------------------------------------------------
435 !-----------------------------------------------------------------------------
436 subroutine add_angpair(ici,icj,nang_pair,iang_pair)
438 ! implicit real*8 (a-h,o-z)
439 ! include 'DIMENSIONS'
440 ! include 'COMMON.IOUNITS'
441 ! include 'COMMON.CHAIN'
442 integer :: ici,icj,nang_pair,iang_pair(2,nres)
443 integer :: i,ian1,ian2
444 ! write (iout,*) "add_angpair: ici",ici," icj",icj,
445 ! & " nang_pair",nang_pair
447 if (ian1.lt.4 .or. ian1.gt.nres) return
449 ! write (iout,*) "ian1",ian1," ian2",ian2
450 if (ian2.lt.4 .or. ian2.gt.nres) return
452 if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
454 nang_pair=nang_pair+1
455 iang_pair(1,nang_pair)=ian1
456 iang_pair(2,nang_pair)=ian2
458 end subroutine add_angpair
459 !-------------------------------------------------------------------------
460 subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,lprn)
462 use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
464 ! implicit real*8 (a-h,o-z)
465 ! include 'DIMENSIONS'
466 ! include 'DIMENSIONS.ZSCOPT'
467 ! include 'DIMENSIONS.COMPAR'
468 ! include 'COMMON.IOUNITS'
469 ! include 'COMMON.VAR'
470 ! include 'COMMON.COMPAR'
471 ! include 'COMMON.CHAIN'
472 ! include 'COMMON.GEO'
473 real(kind=8) :: pinorm,deltang
475 integer :: jfrag,ishif1,ishif2,nn,npart,nn4,nne
476 real(kind=8) :: diffang_max,angn,fract,ff
477 integer :: i,j,nbeg,nend,ll,longest
478 if (lprn) write (iout,'(80(1h*))')
482 npart = npiece(jfrag,1)
484 nne = min0(nend_sup,nres)
485 if (lprn) write (iout,*) "nn4",nn4," nne",nne
487 nbeg = ifrag(1,i,jfrag) + 3 - ishif1
488 if (nbeg.lt.nn4) nbeg=nn4
489 nend = ifrag(2,i,jfrag) + 1 - ishif2
490 if (nend.gt.nne) nend=nne
491 if (nend.ge.nbeg) then
492 nn = nn + nend - nbeg + 1
493 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
494 " nn",nn," ishift1",ishif1," ishift2",ishif2
495 if (lprn) write (iout,*) "angles"
499 write (iout,*) "ishif1",ishif1,j
500 ! deltang = pinorm(phi(j)-phi_ref(j+ishif1))
501 deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),&
502 theta_ref(j+ishif1),phi(j),theta(j-1),theta(j))
503 if (dabs(deltang).gt.diffang_max) then
504 if (ll.gt.longest) longest = ll
509 if (ll.gt.longest) longest = ll
510 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
511 rad2deg*phi_ref(j+ishif1),rad2deg*deltang
512 angn=angn+dabs(deltang)
515 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
516 if (lprn) write (iout,*)"segment",i," longest fragment within",&
517 diffang_max*rad2deg,":",longest," fraction",ff
518 if (ff.lt.fract) fract = ff
526 if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,&
529 end subroutine angnorm
530 !-------------------------------------------------------------------------
531 subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,&
532 diffang_max,anorm,fract)
534 use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
536 ! implicit real*8 (a-h,o-z)
537 ! include 'DIMENSIONS'
538 ! include 'DIMENSIONS.ZSCOPT'
539 ! include 'DIMENSIONS.COMPAR'
540 ! include 'COMMON.IOUNITS'
541 ! include 'COMMON.VAR'
542 ! include 'COMMON.COMPAR'
543 ! include 'COMMON.CHAIN'
544 ! include 'COMMON.GEO'
545 integer :: ncont,icont(2,ncont),longest
546 real(kind=8) :: anorm,diffang_max,fract
547 integer :: npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
548 real(kind=8) :: pinorm
550 integer :: jfrag,ishif1,ishif2
551 integer :: nn,nn4,nne,npart,i,j,jstart,jend,ic1,ic2,idi,iic
552 integer :: nbeg,nend,ll
553 real(kind=8) :: angn,ishifc,deltang,ff
555 if (lprn) write (iout,'(80(1h*))')
557 ! Determine the segments for which angles will be compared
560 nne = min0(nend_sup,nres)
561 if (lprn) write (iout,*) "nn4",nn4," nne",nne
562 npart=npiece(jfrag,1)
565 ! write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
566 if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. &
567 icont(1,1).gt.ifrag(2,i,jfrag)) goto 11
569 do while (jstart.lt.ncont .and. &
570 icont(1,jstart).lt.ifrag(1,i,jfrag))
571 ! write (iout,*) "jstart",jstart," icont",icont(1,jstart),
572 ! & " ifrag",ifrag(1,i,jfrag)
575 ! write (iout,*) "jstart",jstart," icont",icont(1,jstart),
576 ! & " ifrag",ifrag(1,i,jfrag)
577 if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11
580 ifrag_c(1,npiece_c)=icont(1,jstart)
582 do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag))
583 ! write (iout,*) "jend",jend," icont",icont(1,jend),
584 ! & " ifrag",ifrag(2,i,jfrag)
587 ! write (iout,*) "jend",jend," icont",icont(1,jend),
588 ! & " ifrag",ifrag(2,i,jfrag)
590 ifrag_c(2,npiece_c)=icont(1,jend)+1
591 ishift_c(npiece_c)=ishif1
592 ! write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
593 ! & " ic1",ic1," ic2",ic2,
594 ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
596 if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
601 ! write (iout,*) "idi",idi
603 if (icont(2,1).gt.ifrag(2,i,jfrag) .or. &
604 icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12
606 do while (jstart.lt.ncont .and. &
607 icont(2,jstart).lt.ifrag(1,i,jfrag))
608 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
609 ! & " ifrag",ifrag(1,i,jfrag)
612 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
613 ! & " ifrag",ifrag(1,i,jfrag)
614 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
617 ifrag_c(2,npiece_c)=icont(2,jstart)+1
619 do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag))
620 ! write (iout,*) "jend",jend," icont",icont(2,jend),
621 ! & " ifrag",ifrag(2,i,jfrag)
624 ! write (iout,*) "jend",jend," icont",icont(2,jend),
625 ! & " ifrag",ifrag(2,i,jfrag)
626 else if (idi.eq.-1) then
627 if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or. &
628 icont(2,1).lt.ifrag(1,i,jfrag)) goto 12
630 do while (jstart.gt.ncont .and. &
631 icont(2,jstart).lt.ifrag(1,i,jfrag))
632 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
633 ! & " ifrag",ifrag(1,i,jfrag)
636 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
637 ! & " ifrag",ifrag(1,i,jfrag)
638 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
641 ifrag_c(2,npiece_c)=icont(2,jstart)+1
643 do while (jend.lt.ncont .and. &
644 icont(2,jend).gt.ifrag(2,i,jfrag))
645 ! write (iout,*) "jend",jend," icont",icont(2,jend),
646 ! & " ifrag",ifrag(2,i,jfrag)
649 ! write (iout,*) "jend",jend," icont",icont(2,jend),
650 ! & " ifrag",ifrag(2,i,jfrag)
658 ! write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
659 ! & " jstart:",jstart," jend",jend,
660 ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
661 ifrag_c(1,npiece_c)=ic1
662 ifrag_c(2,npiece_c)=ic2+1
663 ishift_c(npiece_c)=ishif2
667 write (iout,*) "Before merge: npiece_c",npiece_c
669 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
673 ! Merge overlapping segments (e.g., avoid splitting helices)
676 do while (i .lt. npiece_c)
677 if (ishift_c(i).eq.ishift_c(i+1) .and. &
678 ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
679 ifrag_c(2,i)=ifrag_c(2,i+1)
681 ishift_c(j)=ishift_c(j+1)
682 ifrag_c(1,j)=ifrag_c(1,j+1)
683 ifrag_c(2,j)=ifrag_c(2,j+1)
691 write (iout,*) "After merge: npiece_c",npiece_c
693 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
706 nbeg = ifrag_c(1,i) + 3 - ishifc
707 if (nbeg.lt.nn4) nbeg=nn4
708 nend = ifrag_c(2,i) - ishifc + 1
709 if (nend.gt.nne) nend=nne
710 if (nend.ge.nbeg) then
711 nn = nn + nend - nbeg + 1
712 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
713 " nn",nn," ishifc",ishifc
714 if (lprn) write (iout,*) "angles"
717 do j=nbeg,nend-ishifc
719 ! deltang = pinorm(phi(j)-phi_ref(j+ishifc))
720 write (iout,*) "i=",j," nn",nn," ishifc",ishifc
721 deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),&
722 theta_ref(j+ishifc),phi(j),theta(j-1),theta(j))
723 if (dabs(deltang).gt.diffang_max) then
724 if (ll.gt.longest) longest = ll
729 if (ll.gt.longest) longest = ll
730 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
731 rad2deg*phi_ref(j+ishifc),rad2deg*deltang
732 angn=angn+dabs(deltang)
735 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
736 if (lprn) write (iout,*)"segment",i," longest fragment within",&
737 diffang_max*rad2deg,":",longest," fraction",ff
738 if (ff.lt.fract) fract = ff
741 if (nn.gt.0) anorm = angn/nn
742 if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
744 end subroutine angnorm2
745 !-------------------------------------------------------------------------
746 real(kind=8) function angnorm1(nang_pair,iang_pair,lprn)
748 use geometry_data, only:phi,theta,rad2deg
749 ! implicit real*8 (a-h,o-z)
750 ! include 'DIMENSIONS'
751 ! include 'DIMENSIONS.ZSCOPT'
752 ! include 'DIMENSIONS.COMPAR'
753 ! include 'COMMON.IOUNITS'
754 ! include 'COMMON.VAR'
755 ! include 'COMMON.COMPAR'
756 ! include 'COMMON.CHAIN'
757 ! include 'COMMON.GEO'
759 integer :: nang_pair,iang_pair(2,nres)
760 real(kind=8) :: pinorm
762 real(kind=8) :: angn,deltang
764 if (lprn) write (iout,'(80(1h*))')
765 if (lprn) write (iout,*) "nang_pair",nang_pair
766 if (lprn) write (iout,*) "angles"
770 ! deltang = pinorm(phi(ia1)-phi_ref(ia2))
771 deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),&
772 theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2))
773 if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),&
774 rad2deg*phi_ref(ia2),rad2deg*deltang
775 angn=angn+dabs(deltang)
778 write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
779 angnorm1 = angn/nang_pair
781 end function angnorm1
782 !------------------------------------------------------------------------------
783 subroutine angnorm12(diff)
785 use geometry_data, only:phi,theta,nstart_sup,nend_sup
786 ! implicit real*8 (a-h,o-z)
787 ! include 'DIMENSIONS'
788 ! include 'DIMENSIONS.ZSCOPT'
789 ! include 'DIMENSIONS.COMPAR'
790 ! include 'COMMON.IOUNITS'
791 ! include 'COMMON.VAR'
792 ! include 'COMMON.COMPAR'
793 ! include 'COMMON.CHAIN'
794 ! include 'COMMON.GEO'
795 real(kind=8) :: pinorm,diff
799 nne = min0(nend_sup,nres)
801 ! diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j)))
804 ! diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j)))
805 diff=diff+spherang(phi_ref(j),theta_ref(j-1),&
806 theta_ref(j),phi(j),theta(j-1),theta(j))
809 end subroutine angnorm12
810 !--------------------------------------------------------------------------------
811 real(kind=8) function spherang(gam1,theta11,theta12,&
812 gam2,theta21,theta22)
814 use geometry, only:arcos
815 real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,&
816 x1,x2,xmed,f1,f2,fmed
817 real(kind=8) :: tolx=1.0d-4, tolf=1.0d-4
818 real(kind=8) :: sumcos
819 !el real(kind=8) :: pinorm,sumangp !arcos,
820 integer :: it,maxit=100
821 ! Calculate the difference of the angles of two superposed 4-redidue fragments
829 ! The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
830 ! to achieve the minimum difference between the O'-C-O and P-C-P angles;
831 ! the sum of these angles is the difference returned by the function.
834 ! If thetas match, take the difference of gamma and exit.
835 if (dabs(theta11-theta12).lt.tolx &
836 .and. dabs(theta21-theta22).lt.tolx) then
837 spherang=dabs(pinorm(gam2-gam1))
840 ! If the gammas are the same, take the difference of thetas and exit.
842 x2=0.5d0*pinorm(gam2-gam1)
843 if (dabs(x2) .lt. tolx) then
844 spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
846 else if (x2.lt.0.0d0) then
850 ! Else apply regula falsi method to compute optimum overlap of the terminal Calphas
851 f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
852 f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
854 xmed=x1-f1*(x2-x1)/(f2-f1)
855 fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
856 ! write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
857 if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx) &
858 .and. dabs(fmed).lt.tolf ) then
862 else if ( fmed*f1.lt.0.0d0 ) then
871 spherang=arcos(dcos(theta11)*dcos(theta12) &
872 +dsin(theta11)*dsin(theta12)*dcos(x1))+ &
873 arcos(dcos(theta21)*dcos(theta22)+ &
874 dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
876 end function spherang
877 !--------------------------------------------------------------------------------
878 real(kind=8) function sumangp(gam1,theta11,theta12,gam2,&
881 real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,fi,&
882 cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,&
884 ! derivarive of the sum of the difference of the angles of a 4-residue fragment.
885 ! real(kind=8) :: arcos
894 cosd1=cost11*cost12+sint11*sint12*dcos(fi)
895 cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
896 sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1) &
897 +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)
900 !-----------------------------------------------------------------------------
902 !-----------------------------------------------------------------------------
903 subroutine contact(lprint,ncont,icont,ist,ien)
906 use geometry_data, only:c,dc,dc_norm
907 use energy_data, only:itype,maxcont,molnum
909 ! include 'DIMENSIONS'
910 ! include 'DIMENSIONS.ZSCOPT'
911 ! include 'COMMON.CONTROL'
912 ! include 'COMMON.IOUNITS'
913 ! include 'COMMON.CHAIN'
914 ! include 'COMMON.INTERACT'
915 ! include 'COMMON.FFIELD'
916 ! include 'COMMON.NAMES'
917 ! include 'COMMON.CALC'
918 ! include 'COMMON.CONTPAR'
919 ! include 'COMMON.LOCAL'
920 integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2,mnum,mnum2
921 real(kind=8) :: csc !el,dist
922 real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,&
924 integer :: ncont,mhum
925 integer,dimension(2,maxcont) :: icont
926 real(kind=8) :: u,v,a(3),b(3),dla,dlb
937 write (iout,110) restyp(itype(i,mnum),mnum),i,c(1,i),c(2,i),&
938 c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),&
939 dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
942 110 format (a,'(',i3,')',9f8.3)
945 iti=iabs(itype(i,mnum))
946 if (iti.le.0 .or. iti.gt.ntyp_molec(mnum)) cycle
949 itj=iabs(itype(j,mnum2))
950 if (itj.le.0 .or. itj.gt.ntyp_molec(mnum2)) cycle
953 xj = c(1,nres+j)-c(1,nres+i)
954 yj = c(2,nres+j)-c(2,nres+i)
955 zj = c(3,nres+j)-c(3,nres+i)
956 dxi = dc_norm(1,nres+i)
957 dyi = dc_norm(2,nres+i)
958 dzi = dc_norm(3,nres+i)
959 dxj = dc_norm(1,nres+j)
960 dyj = dc_norm(2,nres+j)
961 dzj = dc_norm(3,nres+j)
966 ! write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
967 if (icomparfunc.eq.1) then
968 call contfunc(csc,iti,itj)
969 else if (icomparfunc.eq.2) then
970 call scdist(csc,iti,itj)
971 else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
972 csc = dist(nres+i,nres+j)
973 else if (icomparfunc.eq.4) then
974 call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
976 write (*,*) "Error - Unknown sidechain contact function"
977 write (iout,*) "Error - Unknown sidechain contact function"
979 if (csc.lt.sc_cutoff(iti,itj)) then
980 ! write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
981 ! & " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
982 ! & dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
983 ! & dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
985 ! write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
986 ! & sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
987 ! & chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
996 ddsc(ncont)=1.0d0/rij
1003 write (iout,'(a)') 'Contact map:'
1011 if (mnum.eq.0) mnum=1
1012 print *,"CONTACT",i1,i2,mnum,it1,it2
1013 write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') &
1014 i,restyp(it1,mnum),i1,restyp(it2,mnum),i2,cscore(i),&
1015 sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),&
1016 omt1(i),omt2(i),omt12(i)
1020 end subroutine contact
1022 !----------------------------------------------------------------------------
1023 subroutine contact(lprint,ncont,icont)
1025 use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii,molnum
1026 ! include 'DIMENSIONS'
1027 ! include 'COMMON.IOUNITS'
1028 ! include 'COMMON.CHAIN'
1029 ! include 'COMMON.INTERACT'
1030 ! include 'COMMON.FFIELD'
1031 ! include 'COMMON.NAMES'
1032 real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
1033 integer :: ncont,icont(2,maxcont)
1035 integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
1036 real(kind=8) :: rcomp
1037 integer :: mnum,mnum2
1040 ! print *,'nnt=',nnt,' nct=',nct
1043 iti=iabs(itype(i,1))
1046 itj=iabs(itype(j,1))
1048 ! rcomp=sigmaii(iti,itj)+1.0D0
1049 rcomp=facont*sigmaii(iti,itj)
1051 ! rcomp=sigma(iti,itj)+1.0D0
1052 rcomp=facont*sigma(iti,itj)
1055 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
1056 if (dist(nres+i,nres+j).lt.rcomp) then
1064 write (iout,'(a)') 'Contact map:'
1071 write(iout,*) "test",i1,i2,it1,it2
1072 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1073 i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
1077 end subroutine contact
1079 !----------------------------------------------------------------------------
1080 real(kind=8) function contact_fract(ncont,ncont_ref,&
1083 use energy_data, only:maxcont
1085 ! include 'DIMENSIONS'
1086 ! include 'COMMON.IOUNITS'
1087 integer :: i,j,nmatch
1088 integer :: ncont,ncont_ref
1089 integer,dimension(2,maxcont) :: icont,icont_ref
1091 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
1092 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
1093 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
1094 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
1095 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
1098 if (icont(1,i).eq.icont_ref(1,j) .and. &
1099 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
1102 ! print *,' nmatch=',nmatch
1103 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
1104 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
1106 end function contact_fract
1108 !------------------------------------------------------------------------------
1109 subroutine pept_cont(lprint,ncont,icont)
1111 use geometry_data, only:c
1112 use energy_data, only:maxcont,nnt,nct,itype,molnum
1114 ! include 'DIMENSIONS'
1115 ! include 'DIMENSIONS.ZSCOPT'
1116 ! include 'COMMON.IOUNITS'
1117 ! include 'COMMON.CHAIN'
1118 ! include 'COMMON.INTERACT'
1119 ! include 'COMMON.FFIELD'
1120 ! include 'COMMON.NAMES'
1121 integer :: ncont,icont(2,maxcont)
1122 integer :: i,j,k,kkk,i1,i2,it1,it2,mnum
1124 !el real(kind=8) :: dist
1125 real(kind=8) :: rcomp=5.5d0
1128 print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
1131 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
1135 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
1137 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
1145 write (iout,'(a)') 'PP contact map:'
1152 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1153 i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
1157 end subroutine pept_cont
1158 !-----------------------------------------------------------------------------
1160 !-----------------------------------------------------------------------------
1161 subroutine contacts_between_fragments(lprint,is,ncont,icont,&
1162 ncont_interfrag,icont_interfrag)
1164 use energy_data, only:itype,maxcont,molnum
1165 ! include 'DIMENSIONS'
1166 ! include 'DIMENSIONS.ZSCOPT'
1167 ! include 'DIMENSIONS.COMPAR'
1168 ! include 'COMMON.INTERACT'
1169 ! include 'COMMON.COMPAR'
1170 ! include 'COMMON.IOUNITS'
1171 ! include 'COMMON.CHAIN'
1172 ! include 'COMMON.NAMES'
1173 integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),&
1174 icont_interfrag(2,maxcont,mmaxfrag)
1175 logical :: OK1,OK2,lprint
1176 integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2,mnum
1177 ! Determine the contacts that occur within a fragment and between fragments.
1182 ! write (iout,*) "i",i,(ifrag(1,k,i),ifrag(2,k,i)
1183 ! & ,k=1,npiece(i,1))
1184 ! write (iout,*) "j",j,(ifrag(1,k,j),ifrag(2,k,j)
1185 ! & ,k=1,npiece(j,1))
1186 ! write (iout,*) "ncont",ncont
1192 do while (.not.OK1 .and. l.lt.npiece(j,1))
1194 OK1=ic1.ge.ifrag(1,l,j)-is .and. &
1195 ic1.le.ifrag(2,l,j)+is
1199 do while (.not.OK2 .and. l.lt.npiece(i,1))
1201 OK2=ic2.ge.ifrag(1,l,i)-is .and. &
1202 ic2.le.ifrag(2,l,i)+is
1204 ! write(iout,*) "k",k," ic1",ic1," ic2",ic2," OK1",OK1,
1206 if (OK1.and.OK2) then
1208 icont_interfrag(1,nc,ind)=ic1
1209 icont_interfrag(2,nc,ind)=ic2
1210 ! write (iout,*) "nc",nc," ic1",ic1," ic2",ic2
1213 ncont_interfrag(ind)=nc
1214 ! do k=1,ncont_interfrag(ind)
1215 ! i1=icont_interfrag(1,k,ind)
1216 ! i2=icont_interfrag(2,k,ind)
1219 ! write (iout,'(i3,2x,a,i4,2x,a,i4)')
1220 ! & i,restyp(it1),i1,restyp(it2),i2
1225 write (iout,*) "Contacts within fragments:"
1227 write (iout,*) "Fragment",i," (",(ifrag(1,k,i),&
1228 ifrag(2,k,i),k=1,npiece(i,1)),")"
1230 do k=1,ncont_interfrag(ind)
1231 i1=icont_interfrag(1,k,ind)
1232 i2=icont_interfrag(2,k,ind)
1235 it2=itype(i2,molnum(i2))
1236 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1237 i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
1241 write (iout,*) "Contacts between fragments:"
1245 write (iout,*) "Fragments",i," (",(ifrag(1,k,i),&
1246 ifrag(2,k,i),k=1,npiece(i,1)),") and",j," (",&
1247 (ifrag(1,k,j),ifrag(2,k,j),k=1,npiece(j,1)),")"
1248 write (iout,*) "Number of contacts",&
1249 ncont_interfrag(ind)
1251 do k=1,ncont_interfrag(ind)
1252 i1=icont_interfrag(1,k,ind)
1253 i2=icont_interfrag(2,k,ind)
1256 it2=itype(i2,molnum(i2))
1257 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1258 i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
1264 end subroutine contacts_between_fragments
1265 !-----------------------------------------------------------------------------
1267 !-----------------------------------------------------------------------------
1268 subroutine contfunc(cscore,itypi,itypj)
1270 ! This subroutine calculates the contact function based on
1271 ! the Gay-Berne potential of interaction.
1274 ! implicit real*8 (a-h,o-z)
1275 ! include 'DIMENSIONS'
1276 ! include 'COMMON.CONTPAR'
1277 ! include 'COMMON.CALC'
1279 integer :: itypi,itypj
1280 real(kind=8) :: cscore,sig0ij,rrij,sig,rij_shift,evdw,e2
1282 sig0ij=sig_comp(itypi,itypj)
1283 chi1=chi_comp(itypi,itypj)
1284 chi2=chi_comp(itypj,itypi)
1286 chip1=chip_comp(itypi,itypj)
1287 chip2=chip_comp(itypj,itypi)
1289 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1291 ! Calculate angle-dependent terms of the contact function
1295 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1296 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1297 om12=dxi*dxj+dyi*dyj+dzi*dzj
1299 ! print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
1301 ! & rij,rrij,om1,om2,om12
1302 ! Calculate eps1(om12)
1303 faceps1=1.0D0-om12*chiom12
1304 faceps1_inv=1.0D0/faceps1
1305 eps1=dsqrt(faceps1_inv)
1306 ! Following variable is eps1*deps1/dom12
1307 eps1_om12=faceps1_inv*chiom12
1308 ! Calculate sigma(om1,om2,om12)
1312 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1313 sigsq=1.0D0-facsig*faceps1_inv
1314 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1317 chipom12=chip12*om12
1318 facp=1.0D0-om12*chipom12
1320 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1321 ! Following variable is the square root of eps2
1322 eps2rt=1.0D0-facp1*facp_inv
1324 sig=sig0ij*dsqrt(sigsq)
1325 rij_shift=1.0D0/rij-sig+sig0ij
1326 if (rij_shift.le.0.0D0) then
1328 cscore = -dlog(evdw+1.0d-6)
1331 rij_shift=1.0D0/rij_shift
1332 e2=(rij_shift*sig0ij)**expon
1333 evdw=dabs(eps1*eps2rt**2*e2)
1334 if (evdw.gt.1.0d1) evdw = 1.0d1
1335 cscore = -dlog(evdw+1.0d-6)
1337 end subroutine contfunc
1338 !------------------------------------------------------------------------------
1339 subroutine scdist(cscore,itypi,itypj)
1341 ! This subroutine calculates the contact distance
1344 ! implicit real*8 (a-h,o-z)
1345 ! include 'DIMENSIONS'
1346 ! include 'COMMON.CONTPAR'
1347 ! include 'COMMON.CALC'
1348 integer :: itypi,itypj
1349 real(kind=8) :: cscore,rrij
1351 chi1=chi_comp(itypi,itypj)
1352 chi2=chi_comp(itypj,itypi)
1354 rrij=xj*xj+yj*yj+zj*zj
1356 ! Calculate angle-dependent terms of the contact function
1360 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1361 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1362 om12=dxi*dxj+dyi*dyj+dzi*dzj
1367 cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)
1369 end subroutine scdist
1370 !------------------------------------------------------------------------------
1372 !------------------------------------------------------------------------------
1373 subroutine elecont(lprint,ncont,icont,ist,ien)
1375 use geometry_data, only:c
1376 use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2,molnum
1378 ! include 'DIMENSIONS'
1379 ! include 'DIMENSIONS.ZSCOPT'
1380 ! include 'DIMENSIONS.COMPAR'
1381 ! include 'COMMON.IOUNITS'
1382 ! include 'COMMON.CHAIN'
1383 ! include 'COMMON.INTERACT'
1384 ! include 'COMMON.FFIELD'
1385 ! include 'COMMON.NAMES'
1386 ! include 'COMMON.LOCAL'
1388 integer :: i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
1389 real(kind=8) :: rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,&
1390 xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,&
1391 vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,&
1393 real(kind=8),dimension(2,2) :: elpp6c=reshape((/-0.2379d0,&
1394 -0.2056d0,-0.2056d0,-0.0610d0/),shape(elpp6c))
1395 real(kind=8),dimension(2,2) :: elpp3c=reshape((/ 0.0503d0,&
1396 0.0000d0, 0.0000d0, 0.0692d0/),shape(elpp3c))
1397 real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc
1398 real(kind=8) :: elcutoff=-0.3d0
1399 real(kind=8) :: elecutoff_14=-0.5d0
1400 integer :: ncont,icont(2,maxcont),mnum
1401 real(kind=8) :: econt(maxcont)
1403 ! Load the constants of peptide bond - peptide bond interactions.
1404 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
1405 ! proline) - determined by averaging ECEPP energy.
1409 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
1410 ! data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
1411 !el data (elpp6c) /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
1412 !el data (elpp3c) / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
1413 !el data (elcutoff) /-0.3d0/
1414 !el data (elecutoff_14) /-0.5d0/
1417 if (lprint) write (iout,'(a)') &
1418 "Constants of electrostatic interaction energy expression."
1422 appc(i,j)=epp(i,j)*rri*rri
1423 bppc(i,j)=-2.0*epp(i,j)*rri
1424 ael6c(i,j)=elpp6c(i,j)*4.2**6
1425 ael3c(i,j)=elpp3c(i,j)*4.2**3
1427 write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),&
1446 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1447 if (iteli.eq.2 .and. itelj.eq.2 &
1448 .or.iteli.eq.0 .or.itelj.eq.0) goto 4
1449 aaa=appc(iteli,itelj)
1450 bbb=bppc(iteli,itelj)
1451 ael6i=ael6c(iteli,itelj)
1452 ael3i=ael3c(iteli,itelj)
1456 xj=c(1,j)+0.5*dxj-xmedi
1457 yj=c(2,j)+0.5*dyj-ymedi
1458 zj=c(3,j)+0.5*dzj-zmedi
1459 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
1464 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
1465 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
1466 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
1467 fac=cosa-3.0*cosb*cosg
1473 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
1476 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
1477 j.eq.i+2 .and. eesij.le.elecutoff_14) then
1488 write (iout,*) 'Total average electrostatic energy: ',ees
1489 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
1491 write (iout,*) 'Electrostatic contacts before pruning: '
1495 it1=itype(i1,molnum(i1))
1496 it2=itype(i2,molnum(i1))
1497 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1498 i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i1)),i2,econt(i)
1501 ! For given residues keep only the contacts with the greatest energy.
1503 do while (i.lt.ncont)
1509 do while (j.lt.ncont)
1511 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
1512 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
1513 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
1514 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
1515 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
1516 if (ic1.eq.icont(1,j)) then
1518 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)&
1519 .and. iabs(icont(1,k)-ic1).le.2 .and. &
1520 econt(k).lt.econt(j) ) goto 21
1522 else if (ic2.eq.icont(2,j) ) then
1524 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)&
1525 .and. iabs(icont(2,k)-ic2).le.2 .and. &
1526 econt(k).lt.econt(j) ) goto 21
1529 ! Remove ith contact
1531 icont(1,k-1)=icont(1,k)
1532 icont(2,k-1)=icont(2,k)
1537 ! write (iout,*) "ncont",ncont
1539 ! write (iout,*) icont(1,k),icont(2,k)
1542 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
1544 if (ic1.eq.icont(1,j)) then
1546 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
1547 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
1548 econt(k).lt.econt(i) ) goto 21
1550 else if (ic2.eq.icont(2,j) ) then
1552 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
1553 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
1554 econt(k).lt.econt(i) ) goto 21
1557 ! Remove jth contact
1559 icont(1,k-1)=icont(1,k)
1560 icont(2,k-1)=icont(2,k)
1564 ! write (iout,*) "ncont",ncont
1566 ! write (iout,*) icont(1,k),icont(2,k)
1577 write (iout,*) 'Electrostatic contacts after pruning: '
1584 it2=itype(i2,molnum(i2))
1585 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1586 i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2,econt(i)
1590 end subroutine elecont
1591 !------------------------------------------------------------------------------
1593 !------------------------------------------------------------------------------
1594 subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,&
1595 ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,&
1596 nc_frac,nc_req_set,istr,llocal,lprn)
1598 use energy_data, only:maxcont
1599 ! implicit real*8 (a-h,o-z)
1600 ! include 'DIMENSIONS'
1601 ! include 'COMMON.IOUNITS'
1602 integer :: ncont_ref,ncont,ishift,ishif2,nc_match
1603 integer,dimension(2,maxcont) :: icont_ref,icont !(2,maxcont)
1604 real(kind=8) :: nc_frac
1605 logical :: llocal,lprn
1606 integer :: ishif1,nc_match1_max,jfrag,n_shif1,n_shif2,&
1607 nc_req_set,istr,nc_match_max
1608 integer :: i,nc_req,nc_match1,is,js
1611 nc_match_max=nc_match_max+ &
1612 min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
1616 else if (nc_req_set.eq.0) then
1617 nc_req=nc_match_max*nc_frac
1619 nc_req = dmin1(nc_match_max*nc_frac+0.5d0,&
1620 dfloat(nc_req_set)+1.0d-7)
1622 ! write (iout,*) "match_contact: nc_req:",nc_req
1623 ! write (iout,*) "nc_match_max",nc_match_max
1624 ! write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
1625 ! & " n_shif2",n_shif2
1626 ! Match current contact map against reference contact map; exit, if at least
1627 ! half of the contacts match
1628 call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,&
1629 ncont,icont,jfrag,llocal,lprn)
1630 nc_match1_max=nc_match1
1631 if (lprn .and. nc_match.gt.0) write (iout,*) &
1632 "Shift:",0,0," nc_match1",nc_match1,&
1633 " nc_match=",nc_match," req'd",nc_req
1634 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1635 nc_req.eq.0 .and. nc_match.eq.1) then
1640 ! If sufficient matches are not found, try to shift contact maps up to three
1642 if (n_shif1.gt.0) then
1644 ! The following four tries help to find shifted beta-sheet patterns
1645 ! Shift "left" strand backward
1646 call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,&
1647 icont_ref,ncont,icont,jfrag,llocal,lprn)
1648 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1649 if (lprn .and. nc_match.gt.0) write (iout,*) &
1650 "Shift:",-is,0," nc_match1",nc_match1,&
1651 " nc_match=",nc_match," req'd",nc_req
1652 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1653 nc_req.eq.0 .and. nc_match.eq.1) then
1658 ! Shift "left" strand forward
1659 call ncont_match(nc_match,nc_match1,is,0,ncont_ref,&
1660 icont_ref,ncont,icont,jfrag,llocal,lprn)
1661 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1662 if (lprn .and. nc_match.gt.0) write (iout,*) &
1663 "Shift:",is,0," nc_match1",nc_match1,&
1664 " nc_match=",nc_match," req'd",nc_req
1665 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1666 nc_req.eq.0 .and. nc_match.eq.1) then
1672 if (nc_req.eq.0) return
1673 ! Shift "right" strand backward
1675 call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,&
1676 icont_ref,ncont,icont,jfrag,llocal,lprn)
1677 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1678 if (lprn .and. nc_match.gt.0) write (iout,*) &
1679 "Shift:",0,-is," nc_match1",nc_match1,&
1680 " nc_match=",nc_match," req'd",nc_req
1681 if (nc_match.ge.nc_req) then
1686 ! Shift "right" strand upward
1687 call ncont_match(nc_match,nc_match1,0,is,ncont_ref,&
1688 icont_ref,ncont,icont,jfrag,llocal,lprn)
1689 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1690 if (lprn .and. nc_match.gt.0) write (iout,*) &
1691 "Shift:",0,is," nc_match1",nc_match1,&
1692 " nc_match=",nc_match," req'd",nc_req
1693 if (nc_match.ge.nc_req) then
1699 ! Now try to shift both residues in contacts.
1703 call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,&
1704 icont_ref,ncont,icont,jfrag,llocal,lprn)
1705 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1706 if (lprn .and. nc_match.gt.0) write (iout,*) &
1707 "Shift:",-is,-js," nc_match1",nc_match1,&
1708 " nc_match=",nc_match," req'd",nc_req
1709 if (nc_match.ge.nc_req) then
1714 call ncont_match(nc_match,nc_match1,is,js,ncont_ref,&
1715 icont_ref,ncont,icont,jfrag,llocal,lprn)
1716 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1717 if (lprn .and. nc_match.gt.0) write (iout,*) &
1718 "Shift:",is,js," nc_match1",nc_match1,&
1719 " nc_match=",nc_match," req'd",nc_req
1720 if (nc_match.ge.nc_req) then
1726 call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,&
1727 icont_ref,ncont,icont,jfrag,llocal,lprn)
1728 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1729 if (lprn .and. nc_match.gt.0) write (iout,*) &
1730 "Shift:",-js,-is," nc_match1",nc_match1,&
1731 " nc_match=",nc_match," req'd",nc_req
1732 if (nc_match.ge.nc_req) then
1738 call ncont_match(nc_match,nc_match1,js,is,ncont_ref,&
1739 icont_ref,ncont,icont,jfrag,llocal,lprn)
1740 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1741 if (lprn .and. nc_match.gt.0) write (iout,*) &
1742 "Shift:",js,is," nc_match1",nc_match1,&
1743 " nc_match=",nc_match," req'd",nc_req
1744 if (nc_match.ge.nc_req) then
1751 if (is+js.le.n_shif1) then
1752 call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,&
1753 icont_ref,ncont,icont,jfrag,llocal,lprn)
1754 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1755 if (lprn .and. nc_match.gt.0) write (iout,*) &
1756 "Shift:",-is,js," nc_match1",nc_match1,&
1757 " nc_match=",nc_match," req'd",nc_req
1758 if (nc_match.ge.nc_req) then
1764 call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,&
1765 icont_ref,ncont,icont,jfrag,llocal,lprn)
1766 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1767 if (lprn .and. nc_match.gt.0) write (iout,*) &
1768 "Shift:",js,-is," nc_match1",nc_match1,&
1769 " nc_match=",nc_match," req'd",nc_req
1770 if (nc_match.ge.nc_req) then
1781 if (n_shif2.gt.0) then
1783 call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,&
1784 icont_ref,ncont,icont,jfrag,llocal,lprn)
1785 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1786 if (lprn .and. nc_match.gt.0) write (iout,*) &
1787 "Shift:",-is,-is," nc_match1",nc_match1,&
1788 " nc_match=",nc_match," req'd",nc_req
1789 if (nc_match.ge.nc_req) then
1794 call ncont_match(nc_match,nc_match1,is,is,ncont_ref,&
1795 icont_ref,ncont,icont,jfrag,llocal,lprn)
1796 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1797 if (lprn .and. nc_match.gt.0) write (iout,*) &
1798 "Shift:",is,is," nc_match1",nc_match1,&
1799 " nc_match=",nc_match," req'd",nc_req
1800 if (nc_match.ge.nc_req) then
1807 ! If this point is reached, the contact maps are different.
1812 end subroutine match_contact
1813 !-------------------------------------------------------------------------
1814 subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,&
1815 ncont_ref,icont_ref,ncont,icont,jfrag,llocal,lprn)
1817 use energy_data, only:nnt,nct,maxcont
1818 ! implicit real*8 (a-h,o-z)
1819 ! include 'DIMENSIONS'
1820 ! include 'DIMENSIONS.ZSCOPT'
1821 ! include 'DIMENSIONS.COMPAR'
1822 ! include 'COMMON.IOUNITS'
1823 ! include 'COMMON.INTERACT'
1824 ! include 'COMMON.GEO'
1825 ! include 'COMMON.COMPAR'
1826 logical :: llocal,lprn
1827 integer ncont_ref,ncont,ishift,ishif2,nang_pair
1828 integer,dimension(2,maxcont) :: icont_ref,icont,icont_match !(2,maxcont)
1829 integer,dimension(2,nres) :: iang_pair !(2,maxres)
1830 integer :: nc_match,nc_match1,ishif1,jfrag
1831 integer :: i,j,ic1,ic2
1832 real(kind=8) :: diffang,fract,rad2deg
1834 ! Compare the contact map against the reference contact map; they're stored
1835 ! in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
1836 if (lprn) write (iout,'(80(1h*))')
1839 ! Check the local structure by comparing dihedral angles.
1840 ! write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal
1841 if (llocal .and. ncont_ref.eq.0) then
1842 ! If there are no contacts just compare the dihedral angles and exit.
1843 call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,&
1845 if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1846 " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract
1847 if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag)) &
1857 ic1=icont(1,i)+ishif1
1858 ic2=icont(2,i)+ishif2
1859 ! write (iout,*) "i",i," ic1",ic1," ic2",ic2
1860 if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
1862 if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
1863 nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
1864 nc_match1=nc_match1+1
1865 icont_match(1,nc_match1)=ic1
1866 icont_match(2,nc_match1)=ic2
1867 ! call add_angpair(icont(1,i),icont_ref(1,j),
1868 ! & nang_pair,iang_pair)
1869 ! call add_angpair(icont(2,i),icont_ref(2,j),
1870 ! & nang_pair,iang_pair)
1871 if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),&
1872 " match",icont_ref(1,j),icont_ref(2,j),&
1873 " shifts",ishif1,ishif2
1880 write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
1881 write (iout,*) "icont_match"
1883 write (iout,*) icont_match(1,i),icont_match(2,i)
1886 if (llocal .and. nc_match.gt.0) then
1887 call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,&
1888 ang_cut1(jfrag),diffang,fract)
1889 if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1890 " ang_cut:",ang_cut(jfrag)*rad2deg,&
1891 " ang_cut1",ang_cut1(jfrag)*rad2deg
1892 if (diffang.gt.ang_cut(jfrag) &
1893 .or. fract.lt.frac_min(jfrag)) nc_match=0
1895 ! if (nc_match.gt.0) then
1896 ! diffang = angnorm1(nang_pair,iang_pair,lprn)
1897 ! if (diffang.gt.ang_cut(jfrag)) nc_match=0
1899 if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,&
1900 " diffang",rad2deg*diffang," nc_match",nc_match
1902 end subroutine ncont_match
1903 !------------------------------------------------------------------------------
1904 subroutine match_secondary(jfrag,isecstr,nsec_match,lprn)
1905 ! This subroutine compares the secondary structure (isecstr) of fragment jfrag
1906 ! conformation considered to that of the reference conformation.
1907 ! Returns the number of equivalent residues (nsec_match).
1908 ! implicit real*8 (a-h,o-z)
1909 ! include 'DIMENSIONS'
1910 ! include 'DIMENSIONS.ZSCOPT'
1911 ! include 'DIMENSIONS.COMPAR'
1912 ! include 'COMMON.IOUNITS'
1913 ! include 'COMMON.CHAIN'
1914 ! include 'COMMON.PEPTCONT'
1915 ! include 'COMMON.COMPAR'
1917 integer :: isecstr(nres)
1918 integer :: jfrag,nsec_match,npart,i,j
1919 npart = npiece(jfrag,1)
1922 write (iout,*) "match_secondary jfrag",jfrag," ifrag",&
1923 (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart)
1924 write (iout,'(80i1)') (isec_ref(j),j=1,nres)
1925 write (iout,'(80i1)') (isecstr(j),j=1,nres)
1928 do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
1929 ! The residue has equivalent conformational state to that of the reference
1931 ! a) the conformational states are equal or
1932 ! b) the reference state is a coil and that of the conformation considered
1934 ! c) the conformational state of the conformation considered is a strand
1935 ! and that of the reference conformation is a coil.
1936 ! 10/28/02 - case (b) deleted.
1937 if (isecstr(j).eq.isec_ref(j) .or. &
1938 ! & isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or.
1939 isec_ref(j).eq.0 .and. isecstr(j).eq.1) &
1940 nsec_match=nsec_match+1
1944 end subroutine match_secondary
1945 !------------------------------------------------------------------------------
1947 !------------------------------------------------------------------------------
1948 subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd)
1950 use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
1952 ! implicit real*8 (a-h,o-z)
1953 real(kind=8),dimension(3) :: r1,r2,a,b,x,y
1954 real(kind=8) :: uu,vv,aa,bb,dd
1955 real(kind=8) :: ab,ar,br,det,dd1,dd2,dd3,dd4,dd5
1956 !el odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
1957 !el + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
1958 ! print *,"r1",(r1(i),i=1,3)
1959 ! print *,"r2",(r2(i),i=1,3)
1960 ! print *,"a",(a(i),i=1,3)
1961 ! print *,"b",(b(i),i=1,3)
1962 aa = a(1)**2+a(2)**2+a(3)**2
1963 bb = b(1)**2+b(2)**2+b(3)**2
1964 ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
1965 ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3))
1966 br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3))
1968 ! print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det
1969 uu = (-ar*bb+br*ab)/det
1970 vv = (br*aa-ar*ab)/det
1976 !el dd1 = odl(uu,vv)
1977 dd1 = odl(uu,vv,r1,r2,ar,br,ab,aa,bb)
1978 !el dd2 = odl(0.0d0,0.0d0)
1979 dd2 = odl(0.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1980 !el dd3 = odl(0.0d0,1.0d0)
1981 dd3 = odl(0.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1982 !el dd4 = odl(1.0d0,0.0d0)
1983 dd4 = odl(1.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1984 !el dd5 = odl(1.0d0,1.0d0)
1985 dd5 = odl(1.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1986 dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5))
1990 else if (dd.eq.dd3) then
1993 else if (dd.eq.dd4) then
1996 else if (dd.eq.dd5) then
2005 ! dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2
2009 ! write (8,*) uu,vv,dd,dd1
2012 end subroutine odlodc
2013 !------------------------------------------------------------------------------
2014 real(kind=8) function odl(u,v,r1,r2,ar,br,ab,aa,bb)
2016 real(kind=8),dimension(3) :: r1,r2
2017 real(kind=8) :: aa,bb,u,v,ar,br,ab
2019 odl = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
2020 + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
2023 !------------------------------------------------------------------------------
2025 !------------------------------------------------------------------------------
2026 subroutine proc_cont
2028 use geometry_data, only:rad2deg
2029 use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
2031 ! implicit real*8 (a-h,o-z)
2032 ! include 'DIMENSIONS'
2033 ! include 'DIMENSIONS.ZSCOPT'
2034 ! include 'DIMENSIONS.COMPAR'
2035 ! include 'COMMON.IOUNITS'
2036 ! include 'COMMON.TIME1'
2037 ! include 'COMMON.SBRIDGE'
2038 ! include 'COMMON.CONTROL'
2039 ! include 'COMMON.COMPAR'
2040 ! include 'COMMON.CHAIN'
2041 ! include 'COMMON.HEADER'
2042 ! include 'COMMON.CONTACTS1'
2043 ! include 'COMMON.PEPTCONT'
2044 ! include 'COMMON.GEO'
2045 integer :: i,j,k,ind,len_cut,ndigit,length_frag
2047 write (iout,*) "proc_cont: nlevel",nlevel
2048 if (nlevel.lt.0) then
2049 write (iout,*) "call define_fragments"
2050 call define_fragments
2052 write (iout,*) "call secondary2"
2053 call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2056 write (iout,'(80(1h=))')
2057 write (iout,*) "Electrostatic contacts"
2058 call contacts_between_fragments(.true.,0,ncont_pept_ref,&
2059 icont_pept_ref,ncont_frag_ref(1),icont_frag_ref(1,1,1))
2060 write (iout,'(80(1h=))')
2061 write (iout,*) "Side chain contacts"
2062 call contacts_between_fragments(.true.,0,ncont_ref,&
2063 icont_ref,nsccont_frag_ref(1),isccont_frag_ref(1,1,1))
2064 if (nlevel.lt.0) then
2068 if (istruct(i).le.1) then
2069 len_cut=max0(len_frag(i,1)*4/5,3)
2070 else if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2071 len_cut=max0(len_frag(i,1)*2/5,3)
2073 write (iout,*) "i",i," istruct",istruct(i)," ncont_frag",&
2074 ncont_frag_ref(ind)," len_cut",len_cut,&
2075 " icont_single",icont_single," iloc_single",iloc_single
2077 if (iloc(i).gt.0) write (iout,*) &
2078 "Local structure used to compare structure of fragment",i,&
2080 if (istruct(i).ne.3 .and. istruct(i).ne.0 &
2081 .and. icont_single.gt.0 .and. &
2082 ncont_frag_ref(ind).ge.len_cut) then
2083 write (iout,*) "Electrostatic contacts used to compare",&
2084 " structure of fragment",i," to native."
2087 else if (icont_single.gt.0 .and. nsccont_frag_ref(ind) &
2089 write (iout,*) "Side chain contacts used to compare",&
2090 " structure of fragment",i," to native."
2094 write (iout,*) "Contacts not used to compare",&
2095 " structure of fragment",i," to native."
2100 if (irms_single.gt.0 .or. isccont(i,1).eq.0 &
2101 .and. ielecont(i,1).eq.0) then
2102 write (iout,*) "RMSD used to compare",&
2103 " structure of fragment",i," to native."
2106 write (iout,*) "RMSD not used to compare",&
2107 " structure of fragment",i," to native."
2112 if (nlevel.lt.-1) then
2115 if (nlevel.gt.3) nlevel=3
2116 if (nlevel.eq.3) then
2118 npiece(1,3)=nfrag(1)
2128 else if (nlevel.eq.-1) then
2133 isnfrag(i+1)=isnfrag(i)+nfrag(i)
2137 ndigit=ndigit+2*nfrag(i)
2139 write (iout,*) "ndigit",ndigit
2140 if (.not.binary .and. ndigit.gt.30) then
2141 write (iout,*) "Highest class too large; switching to",&
2142 " binary representation."
2145 write (iout,*) "isnfrag",(isnfrag(i),i=1,nlevel+1)
2146 write(iout,*) "rmscut_base_up",rmscut_base_up,&
2147 " rmscut_base_low",rmscut_base_low," rmsup_lim",rmsup_lim
2153 length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2157 length_frag=length_frag+len_frag(ipiece(k,j,i),1)
2160 len_frag(j,i)=length_frag
2161 rmscutfrag(1,j,i)=rmscut_base_up*length_frag
2162 rmscutfrag(2,j,i)=rmscut_base_low*length_frag
2163 if (rmscutfrag(1,j,i).lt.rmsup_lim) &
2164 rmscutfrag(1,j,i)=rmsup_lim
2165 if (rmscutfrag(1,j,i).gt.rmsupup_lim) &
2166 rmscutfrag(1,j,i)=rmsupup_lim
2169 write (iout,*) "Level",1," number of fragments:",nfrag(1)
2171 write (iout,*) npiece(j,1),(ifrag(1,k,j),ifrag(2,k,j),&
2172 k=1,npiece(j,1)),len_frag(j,1),rmscutfrag(1,j,1),&
2173 rmscutfrag(2,j,1),n_shift(1,j,1),n_shift(2,j,1),&
2174 ang_cut(j)*rad2deg,ang_cut1(j)*rad2deg,frac_min(j),&
2175 nc_fragm(j,1),nc_req_setf(j,1),istruct(j)
2178 write (iout,*) "Level",i," number of fragments:",nfrag(i)
2180 write (iout,*) npiece(j,i),(ipiece(k,j,i),&
2181 k=1,npiece(j,i)),len_frag(j,i),rmscutfrag(1,j,i),&
2182 rmscutfrag(2,j,i),n_shift(1,j,i),n_shift(2,j,i),&
2183 nc_fragm(j,i),nc_req_setf(j,i)
2187 end subroutine proc_cont
2188 !------------------------------------------------------------------------------
2190 !------------------------------------------------------------------------------
2191 subroutine define_pairs
2193 ! use energy_data, only:nsccont_frag_ref
2194 ! implicit real*8 (a-h,o-z)
2195 ! include 'DIMENSIONS'
2196 ! include 'DIMENSIONS.ZSCOPT'
2197 ! include 'DIMENSIONS.COMPAR'
2198 ! include 'COMMON.IOUNITS'
2199 ! include 'COMMON.TIME1'
2200 ! include 'COMMON.SBRIDGE'
2201 ! include 'COMMON.CONTROL'
2202 ! include 'COMMON.COMPAR'
2203 ! include 'COMMON.FRAG'
2204 ! include 'COMMON.CHAIN'
2205 ! include 'COMMON.HEADER'
2206 ! include 'COMMON.GEO'
2207 ! include 'COMMON.CONTACTS1'
2208 ! include 'COMMON.PEPTCONT'
2209 integer :: j,k,i,length_frag,ind,ll1,ll2,len_cut
2214 length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2216 len_frag(j,1)=length_frag
2217 write (iout,*) "Fragment",j," length",len_frag(j,1)
2223 if (istruct(i).le.1 .or. istruct(j).le.1) then
2224 if (istruct(i).le.1) then
2229 if (istruct(j).le.1) then
2234 len_cut=max0(min0(ll1*2/3,ll2*4/5),3)
2236 if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2241 if (istruct(j).eq.2 .or. istruct(j).eq.4) then
2246 len_cut=max0(min0(ll1*4/5,ll2)*4/5,3)
2248 write (iout,*) "Fragments",i,j," structure",istruct(i),&
2249 istruct(j)," # contacts",&
2250 ncont_frag_ref(ind),nsccont_frag_ref(ind),&
2251 " lengths",len_frag(i,1),len_frag(j,1),&
2252 " ll1",ll1," ll2",ll2," len_cut",len_cut
2253 if ((istruct(i).eq.1 .or. istruct(j).eq.1) .and. &
2254 nsccont_frag_ref(ind).ge.len_cut ) then
2255 if (istruct(i).eq.1 .and. istruct(j).eq.1) then
2256 write (iout,*) "Adding pair of helices",i,j,&
2257 " based on SC contacts"
2259 write (iout,*) "Adding helix+strand/sheet pair",i,j,&
2260 " based on SC contacts"
2263 if (icont_pair.gt.0) then
2264 write (iout,*) "# SC contacts will be used",&
2266 isccont(nfrag(2),2)=1
2268 if (irms_pair.gt.0) then
2269 write (iout,*) "Fragment RMSD will be used",&
2273 npiece(nfrag(2),2)=2
2274 ipiece(1,nfrag(2),2)=i
2275 ipiece(2,nfrag(2),2)=j
2276 ielecont(nfrag(2),2)=0
2277 n_shift(1,nfrag(2),2)=nshift_pair
2278 n_shift(2,nfrag(2),2)=nshift_pair
2279 nc_fragm(nfrag(2),2)=ncfrac_pair
2280 nc_req_setf(nfrag(2),2)=ncreq_pair
2281 else if ((istruct(i).ge.2 .and. istruct(i).le.4) &
2282 .and. (istruct(j).ge.2 .and. istruct(i).le.4) &
2283 .and. ncont_frag_ref(ind).ge.len_cut ) then
2285 write (iout,*) "Adding pair strands/sheets",i,j,&
2286 " based on pp contacts"
2287 if (icont_pair.gt.0) then
2288 write (iout,*) "# pp contacts will be used",&
2290 ielecont(nfrag(2),2)=1
2292 if (irms_pair.gt.0) then
2293 write (iout,*) "Fragment RMSD will be used",&
2297 npiece(nfrag(2),2)=2
2298 ipiece(1,nfrag(2),2)=i
2299 ipiece(2,nfrag(2),2)=j
2300 ielecont(nfrag(2),2)=1
2301 isccont(nfrag(2),2)=0
2302 n_shift(1,nfrag(2),2)=nshift_pair
2303 n_shift(2,nfrag(2),2)=nshift_pair
2304 nc_fragm(nfrag(2),2)=ncfrac_bet
2305 nc_req_setf(nfrag(2),2)=ncreq_bet
2309 write (iout,*) "Pairs found"
2311 write (iout,*) ipiece(1,i,2),ipiece(2,i,2)
2314 end subroutine define_pairs
2315 !------------------------------------------------------------------------------
2317 !------------------------------------------------------------------------------
2318 INTEGER FUNCTION ICANT(I,J)
2327 !------------------------------------------------------------------------------
2329 !------------------------------------------------------------------------------
2330 subroutine imysort(n, m, mm, x, y, z, z1, z2, z3, z4, z5, z6)
2333 integer :: x(m,mm,n),y(n),z(n),z1(2,n),z6(n),xmin,xtemp
2334 real(kind=8) :: z2(n),z3(n),z4(n),z5(n)
2335 real(kind=8) :: xxtemp
2336 integer :: i,j,k,imax
2341 if (x(1,1,j).lt.xmin) then
2375 x(j,k,imax)=x(j,k,i)
2381 end subroutine imysort
2382 !------------------------------------------------------------------------------
2384 !-------------------------------------------------------------------------------
2385 real(kind=8) function qwolynes(ilevel,jfrag)
2387 use geometry_data, only:cref,nperm
2388 use control_data, only:symetr
2389 use energy_data, only:nnt,nct,itype,molnum
2391 ! include 'DIMENSIONS'
2392 ! include 'DIMENSIONS.ZSCOPT'
2393 ! include 'DIMENSIONS.COMPAR'
2394 ! include 'COMMON.IOUNITS'
2395 ! include 'COMMON.COMPAR'
2396 ! include 'COMMON.CHAIN'
2397 ! include 'COMMON.INTERACT'
2398 ! include 'COMMON.CONTROL'
2399 integer :: ilevel,jfrag,kkk
2400 integer :: i,j,jl,k,l,il,kl,nl,np,ip,kp
2402 real(kind=8),dimension(:),allocatable :: tempus !(maxperm)
2403 real(kind=8) :: maxiQ !dist,
2404 real(kind=8) :: qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
2405 logical :: lprn=.false.
2406 real(kind=8) :: x !el sigm
2407 !el sigm(x)=0.25d0*x
2413 ! write (iout,*) "QWolyes: " jfrag",jfrag,
2414 ! & " ilevel",ilevel
2415 allocate(tempus(nperm))
2418 if (ilevel.eq.0) then
2419 if (lprn) write (iout,*) "Q computed for whole molecule"
2430 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2431 (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2432 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2434 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2435 if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
2438 (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2439 (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2440 (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2441 dijCM=dist(il+nres,jl+nres)
2442 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2446 write (iout,*) "il",il," jl",jl,&
2447 " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
2448 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2449 " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2454 if (lprn) write (iout,*) "nl",nl," qq",qq
2455 else if (ilevel.eq.1) then
2456 if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag
2458 ! write (iout,*) "nlist_frag",nlist_frag(jfrag)
2459 do i=2,nlist_frag(jfrag)
2461 il=list_frag(i,jfrag)
2462 jl=list_frag(j,jfrag)
2463 if (iabs(il-jl).gt.nsep) then
2471 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2472 (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2473 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2475 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2476 if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
2479 (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2480 (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2481 (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2482 dijCM=dist(il+nres,jl+nres)
2483 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2487 write (iout,*) "i",i," j",j," il",il," jl",jl,&
2488 " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
2489 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2490 " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2496 if (lprn) write (iout,*) "nl",nl," qq",qq
2497 else if (ilevel.eq.2) then
2498 np=npiece(jfrag,ilevel)
2501 ip=ipiece(i,jfrag,ilevel)
2502 do j=1,nlist_frag(ip)
2505 kp=ipiece(k,jfrag,ilevel)
2506 do l=1,nlist_frag(kp)
2508 if (iabs(kl-il).gt.nsep) then
2516 d0ij=dsqrt((cref(1,kl,kkk)-cref(1,il,kkk))**2+ &
2517 (cref(2,kl,kkk)-cref(2,il,kkk))**2+ &
2518 (cref(3,kl,kkk)-cref(3,il,kkk))**2)
2520 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2521 if (itype(il,molnum(il)).ne.10 .or. itype(kl,molnum(kl)).ne.10) then
2524 (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2525 (cref(2,kl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2526 (cref(3,kl+nres,kkk)-cref(3,il+nres,kkk))**2)
2527 dijCM=dist(il+nres,kl+nres)
2528 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/ &
2533 write (iout,*) "i",i," j",j," k",k," l",l," il",il,&
2534 " kl",kl," itype",itype(il,molnum(il)), &
2535 itype(kl,molnum(kl))
2536 write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",&
2537 d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2545 if (lprn) write (iout,*) "nl",nl," qq",qq
2547 write (iout,*)"Error: Q can be computed only for level 1 and 2."
2552 if (maxiQ.le.tempus(kkk)) maxiQ=tempus(kkk)
2554 qwolynes=1.0d0-maxiQ
2557 end function qwolynes
2558 !-------------------------------------------------------------------------------
2559 real(kind=8) function sigm(x)
2564 !-------------------------------------------------------------------------------
2565 subroutine fragment_list
2567 ! include 'DIMENSIONS'
2568 ! include 'DIMENSIONS.ZSCOPT'
2569 ! include 'DIMENSIONS.COMPAR'
2570 ! include 'COMMON.IOUNITS'
2571 ! include 'COMMON.COMPAR'
2572 logical :: lprn=.true.
2573 integer :: i,ilevel,j,k,jfrag
2576 do i=1,npiece(jfrag,1)
2577 if (lprn) write (iout,*) "jfrag=",jfrag,&
2578 "i=",i," fragment",ifrag(1,i,jfrag),&
2580 do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
2581 do k=1,nlist_frag(jfrag)
2582 if (list_frag(k,jfrag).eq.j) goto 10
2584 nlist_frag(jfrag)=nlist_frag(jfrag)+1
2585 list_frag(nlist_frag(jfrag),jfrag)=j
2590 write (iout,*) "Fragment list"
2592 write (iout,*)"Fragment",j," list",(list_frag(k,j),&
2596 end subroutine fragment_list
2597 !-------------------------------------------------------------------------------
2598 real(kind=8) function rmscalc(ishif,i,j,jcon,lprn)
2601 use control_data, only:symetr
2602 use geometry_data, only:nperm
2603 ! implicit real*8 (a-h,o-z)
2604 ! include 'DIMENSIONS'
2605 ! include 'DIMENSIONS.ZSCOPT'
2606 ! include 'DIMENSIONS.COMPAR'
2607 ! include 'COMMON.IOUNITS'
2608 ! include 'COMMON.COMPAR'
2609 ! include 'COMMON.CHAIN'
2610 ! include 'COMMON.INTERACT'
2611 ! include 'COMMON.VAR'
2612 ! include 'COMMON.CONTROL'
2613 real(kind=8) :: przes(3),obrot(3,3)
2614 !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2615 !el logical :: iadded(nres)
2616 !el integer :: inumber(2,nres)
2617 !el common /ccc/ creff,cc,iadded,inumber
2620 integer :: ishif,i,j,jcon,idup,kkk,l,k,kk
2621 real(kind=8) :: rminrms,rms
2623 write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
2624 write (iout,*) "npiece",npiece(j,i)
2627 ! write (iout,*) "symetr",symetr
2633 ! write (iout,*) "nperm",nperm
2640 ! write (iout,*) "kkk",kkk
2645 write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",&
2646 ifrag(1,k,j),ifrag(2,k,j)
2649 call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
2650 ! write (iout,*) "Exit cprep"
2652 ! write (iout,*) "ii=",ii
2655 ! write (iout,*) "kk",kk," npiece",npiece(kk,1)
2658 write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,&
2659 " l=",l," adding fragment",&
2660 ifrag(1,l,kk),ifrag(2,l,kk)
2663 call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
2664 ! write (iout,*) "After cprep"
2671 write (iout,*) "tuszukaj"
2674 write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),&
2675 inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
2682 call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
2684 print *,'Error: FITSQ non-convergent, jcon',jcon,i
2686 else if (rms.lt.-1.0d-6) then
2687 print *,'Error: rms^2 = ',rms,jcon,i
2689 else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2692 ! write (iout,*) "rmsmin", rminrms, "rms", rms
2693 if (rms.le.rminrms) rminrms=rms
2695 rmscalc = dsqrt(rminrms)
2696 ! write (iout, *) "analysys", rmscalc,anatemp
2698 end function rmscalc
2699 !-------------------------------------------------------------------------
2700 subroutine cprep(if1,if2,ishif,idup,kwa)
2703 use control_data, only:symetr
2704 use geometry_data, only:nperm,cref,c
2705 ! implicit real*8 (a-h,o-z)
2706 ! include 'DIMENSIONS'
2707 ! include 'DIMENSIONS.ZSCOPT'
2708 ! include 'DIMENSIONS.COMPAR'
2709 ! include 'COMMON.CONTROL'
2710 ! include 'COMMON.IOUNITS'
2711 ! include 'COMMON.COMPAR'
2712 ! include 'COMMON.CHAIN'
2713 ! include 'COMMON.INTERACT'
2714 ! include 'COMMON.VAR'
2715 real(kind=8) :: przes(3),obrot(3,3)
2716 !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2717 !el logical :: iadded(nres)
2718 !el integer :: inumber(2,nres)
2719 integer :: iistrart,kwa,blar
2720 !el common /ccc/ creff,cc,iadded,inumber
2721 integer :: if1,if2,ishif,idup,kkk,l,m
2722 ! write (iout,*) "Calling cprep symetr",symetr," kwa",kwa
2727 ! write (iout,*) "nperm",nperm
2731 ! write (iout,*) "l",l," iadded",iadded(l)
2733 if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) &
2738 inumber(2,idup)=l+ishif
2740 creff(m,idup)=cref(m,l,kkk)
2741 cc(m,idup)=c(m,l+ishif)
2746 end subroutine cprep
2747 !-------------------------------------------------------------------------
2748 real(kind=8) function rmsnat(jcon)
2750 use control_data, only:symetr
2751 use geometry_data, only:nperm,cref,c
2752 use energy_data, only:itype,molnum
2753 ! implicit real*8 (a-h,o-z)
2754 ! include 'DIMENSIONS'
2755 ! include 'DIMENSIONS.ZSCOPT'
2756 ! include 'DIMENSIONS.COMPAR'
2757 ! include 'COMMON.IOUNITS'
2758 ! include 'COMMON.COMPAR'
2759 ! include 'COMMON.CHAIN'
2760 ! include 'COMMON.INTERACT'
2761 ! include 'COMMON.VAR'
2762 ! include 'COMMON.CONTROL'
2763 real(kind=8) :: przes(3),obrot(3,3),cc(3,2*nres),ccref(3,2*nres)
2765 integer :: ishif,i,j,resprzesun,jcon,kkk,nnsup
2766 real(kind=8) :: rminrms,rmsminsing,rms
2776 if (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then
2780 ccref(j,nnsup)=cref(j,i,kkk)
2784 call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv)
2786 print *,'Error: FITSQ non-convergent, jcon',jcon,i
2788 else if (rms.lt.-1.0d-6) then
2789 print *,'Error: rms^2 = ',rms,jcon,i
2791 else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2794 if (rms.le.rminrms) rminrms=rms
2795 ! write (iout,*) "kkk",kkk," rmsnat",rms , rminrms
2797 rmsnat = dsqrt(rminrms)
2798 ! write (iout,*) "analysys",rmsnat, anatemp
2799 ! liczenie rmsdla pojedynczego lancucha
2802 !-------------------------------------------------------------------------------
2803 subroutine define_fragments
2805 use geometry_data, only:rad2deg
2806 use energy_data, only:itype,molnum
2807 use compare_data, only:nhfrag,nbfrag,bfrag,hfrag
2808 ! implicit real*8 (a-h,o-z)
2809 ! include 'DIMENSIONS'
2810 ! include 'DIMENSIONS.ZSCOPT'
2811 ! include 'DIMENSIONS.COMPAR'
2812 ! include 'COMMON.IOUNITS'
2813 ! include 'COMMON.TIME1'
2814 ! include 'COMMON.FRAG'
2815 ! include 'COMMON.SBRIDGE'
2816 ! include 'COMMON.CONTROL'
2817 ! include 'COMMON.COMPAR'
2818 ! include 'COMMON.CHAIN'
2819 ! include 'COMMON.HEADER'
2820 ! include 'COMMON.GEO'
2821 ! include 'COMMON.CONTACTS'
2822 ! include 'COMMON.PEPTCONT'
2823 ! include 'COMMON.INTERACT'
2824 ! include 'COMMON.NAMES'
2825 integer :: nstrand,istrand(2,nres/2)
2826 integer :: nhairp,ihairp(2,nres/5),mnum
2827 character(len=16) :: strstr(4)=reshape((/'helix','hairpin',&
2828 'strand','strand pair'/),shape(strstr))
2829 integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4
2831 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2832 'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2833 'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2834 ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet
2835 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2836 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2837 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2838 ' MAXANG_HEL',angcut1_hel*rad2deg
2839 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2840 ' MAXANG_BET',angcut1_bet*rad2deg
2841 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2842 ' MAXANG_STRAND',angcut1_strand*rad2deg
2843 write (iout,*) 'FRAC_MIN',frac_min_set
2844 ! Find secondary structure elements (helices and beta-sheets)
2845 call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2847 ! Define primary fragments. First include the helices.
2851 ! AL 12/23/03 - to avoid splitting helices into very small fragments
2852 if (merge_helices) then
2853 write (iout,*) "Before merging helices: nhfrag",nhfrag
2855 write (2,*) hfrag(1,i),hfrag(2,i)
2858 do while (i.lt.nhfrag)
2859 if (hfrag(1,i+1)-hfrag(2,i).le.1) then
2861 hfrag(2,i)=hfrag(2,i+1)
2863 hfrag(1,j)=hfrag(1,j+1)
2864 hfrag(2,j)=hfrag(2,j+1)
2869 write (iout,*) "After merging helices: nhfrag",nhfrag
2871 write (2,*) hfrag(1,i),hfrag(2,i)
2877 ifrag(1,1,i)=hfrag(1,i)
2878 ifrag(2,1,i)=hfrag(2,i)
2880 n_shift(2,i,1)=nshift_hel
2881 ang_cut(i)=angcut_hel
2882 ang_cut1(i)=angcut1_hel
2883 frac_min(i)=frac_min_set
2884 nc_fragm(i,1)=ncfrac_hel
2885 nc_req_setf(i,1)=ncreq_hel
2888 write (iout,*) "isplit_bet",isplit_bet
2889 if (isplit_bet.gt.1) then
2890 ! Split beta-sheets into strands and store strands as primary fragments.
2891 call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2895 ifrag(1,1,ii)=istrand(1,i)
2896 ifrag(2,1,ii)=istrand(2,i)
2897 n_shift(1,ii,1)=nshift_strand
2898 n_shift(2,ii,1)=nshift_strand
2899 ang_cut(ii)=angcut_strand
2900 ang_cut1(ii)=angcut1_strand
2901 frac_min(ii)=frac_min_set
2906 nfrag(1)=nfrag(1)+nstrand
2907 else if (isplit_bet.eq.1) then
2908 ! Split only far beta-sheets; does not split hairpins.
2909 call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2910 call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2914 ifrag(1,1,ii)=ihairp(1,i)
2915 ifrag(2,1,ii)=ihairp(2,i)
2916 n_shift(1,ii,1)=nshift_bet
2917 n_shift(2,ii,1)=nshift_bet
2918 ang_cut(ii)=angcut_bet
2919 ang_cut1(ii)=angcut1_bet
2920 frac_min(ii)=frac_min_set
2921 nc_fragm(ii,1)=ncfrac_bet
2922 nc_req_setf(ii,1)=ncreq_bet
2925 nfrag(1)=nfrag(1)+nhairp
2929 ifrag(1,1,ii)=istrand(1,i)
2930 ifrag(2,1,ii)=istrand(2,i)
2931 n_shift(1,ii,1)=nshift_strand
2932 n_shift(2,ii,1)=nshift_strand
2933 ang_cut(ii)=angcut_strand
2934 ang_cut1(ii)=angcut1_strand
2935 frac_min(ii)=frac_min_set
2940 nfrag(1)=nfrag(1)+nstrand
2942 ! Do not split beta-sheets; each pair of strands is a primary element.
2943 call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2947 ifrag(1,1,ii)=ihairp(1,i)
2948 ifrag(2,1,ii)=ihairp(2,i)
2949 n_shift(1,ii,1)=nshift_bet
2950 n_shift(2,ii,1)=nshift_bet
2951 ang_cut(ii)=angcut_bet
2952 ang_cut1(ii)=angcut1_bet
2953 frac_min(ii)=frac_min_set
2954 nc_fragm(ii,1)=ncfrac_bet
2955 nc_req_setf(ii,1)=ncreq_bet
2958 nfrag(1)=nfrag(1)+nhairp
2962 ifrag(1,1,ii)=bfrag(1,i)
2963 ifrag(2,1,ii)=bfrag(2,i)
2964 if (bfrag(3,i).lt.bfrag(4,i)) then
2965 ifrag(1,2,ii)=bfrag(3,i)
2966 ifrag(2,2,ii)=bfrag(4,i)
2968 ifrag(1,2,ii)=bfrag(4,i)
2969 ifrag(2,2,ii)=bfrag(3,i)
2971 n_shift(1,ii,1)=nshift_bet
2972 n_shift(2,ii,1)=nshift_bet
2973 ang_cut(ii)=angcut_bet
2974 ang_cut1(ii)=angcut1_bet
2975 frac_min(ii)=frac_min_set
2976 nc_fragm(ii,1)=ncfrac_bet
2977 nc_req_setf(ii,1)=ncreq_bet
2980 nfrag(1)=nfrag(1)+nbfrag
2982 write (iout,*) "The following primary fragments were found:"
2983 write (iout,*) "Helices:",nhfrag
2990 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
2991 i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
2993 write (iout,*) "Hairpins:",nhairp
2994 do i=nhfrag+1,nhfrag+nhairp
2997 it1=itype(i1,molnum(i1))
2998 it2=itype(i2,molnum(i2))
2999 write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') &
3000 i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2
3002 write (iout,*) "Far strand pairs:",nbfrag
3003 do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
3006 it1=itype(i1,molnum(i1))
3007 it2=itype(i2,molnum(i1))
3010 it3=itype(i3,molnum(i3))
3011 it4=itype(i4,molnum(i4))
3012 write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') &
3013 i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2,&
3014 restyp(it3,molnum(i3)),i3,restyp(it4,molnum(i4)),i4
3016 write (iout,*) "Strands:",nstrand
3017 do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
3023 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
3024 i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
3026 call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),&
3027 istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),&
3028 nc_fragm(1,1),nc_req_setf(1,1))
3029 write (iout,*) "Fragments after sorting:"
3036 write (iout,'(i3,2x,a,i4,2x,a,i4,$)') &
3037 i,restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2
3038 if (npiece(i,1).eq.1) then
3039 write (iout,'(2x,a)') strstr(istruct(i))
3045 write (iout,'(2x,a,i4,2x,a,i4,2x,a)') &
3046 restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2,strstr(istruct(i))
3050 end subroutine define_fragments
3051 !------------------------------------------------------------------------------
3052 subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
3053 ! implicit real*8 (a-h,o-z)
3054 ! include 'DIMENSIONS'
3055 ! include 'DIMENSIONS.ZSCOPT'
3056 ! include 'DIMENSIONS.COMPAR'
3057 ! include 'COMMON.IOUNITS'
3058 integer :: nbfrag,bfrag(4,nres/3)
3059 integer :: nhairp,ihairp(2,nres/5)
3061 write (iout,*) "Entered find_and_remove_hairpins"
3062 write (iout,*) "nbfrag",nbfrag
3064 write (iout,*) i,(bfrag(k,i),k=1,4)
3068 do while (i.le.nbfrag)
3069 write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
3070 if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) &
3072 write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
3074 ihairp(1,nhairp)=bfrag(1,i)
3075 ihairp(2,nhairp)=bfrag(3,i)
3079 bfrag(k,j)=bfrag(k,j+1)
3086 write (iout,*) "After finding hairpins:"
3087 write (iout,*) "nhairp",nhairp
3089 write (iout,*) i,ihairp(1,i),ihairp(2,i)
3091 write (iout,*) "nbfrag",nbfrag
3093 write (iout,*) i,(bfrag(k,i),k=1,4)
3096 end subroutine find_and_remove_hairpins
3097 !------------------------------------------------------------------------------
3098 subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
3099 ! implicit real*8 (a-h,o-z)
3100 ! include 'DIMENSIONS'
3101 ! include 'DIMENSIONS.ZSCOPT'
3102 ! include 'DIMENSIONS.COMPAR'
3103 ! include 'COMMON.IOUNITS'
3104 integer :: nbfrag,bfrag(4,nres/3)
3105 integer :: nstrand,istrand(2,nres/2)
3106 integer :: nhairp,ihairp(2,nres/5)
3109 write (iout,*) "Entered split_beta"
3110 write (iout,*) "nbfrag",nbfrag
3112 write (iout,*) i,(bfrag(k,i),k=1,4)
3116 write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
3117 call add_strand(nstrand,istrand,nhairp,ihairp,&
3118 bfrag(1,i),bfrag(2,i),found)
3119 if (bfrag(3,i).lt.bfrag(4,i)) then
3120 write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
3121 call add_strand(nstrand,istrand,nhairp,ihairp,&
3122 bfrag(3,i),bfrag(4,i),found)
3124 write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
3125 call add_strand(nstrand,istrand,nhairp,ihairp,&
3126 bfrag(4,i),bfrag(3,i),found)
3130 write (iout,*) "Strands found:",nstrand
3132 write (iout,*) i,istrand(1,i),istrand(2,i)
3135 end subroutine split_beta
3136 !------------------------------------------------------------------------------
3137 subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
3138 ! implicit real*8 (a-h,o-z)
3139 ! include 'DIMENSIONS'
3140 ! include 'DIMENSIONS.ZSCOPT'
3141 ! include 'DIMENSIONS.COMPAR'
3142 ! include 'COMMON.IOUNITS'
3143 integer :: nstrand,istrand(2,nres/2)
3144 integer :: nhairp,ihairp(2,nres/5)
3146 integer :: is1,is2,j,idelt
3149 idelt=(ihairp(2,j)-ihairp(1,j))/6
3150 if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
3151 write (iout,*) "strand",is1,is2," is part of hairpin",&
3152 ihairp(1,j),ihairp(2,j)
3157 idelt=(istrand(2,j)-istrand(1,j))/3
3158 if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) &
3160 ! The strand already exists in the array; update its ends if necessary.
3161 write (iout,*) "strand",is1,is2," found at position",j,&
3162 ":",istrand(1,j),istrand(2,j)
3163 istrand(1,j)=min0(istrand(1,j),is1)
3164 istrand(2,j)=max0(istrand(2,j),is2)
3168 ! The strand has not been found; add it to the array.
3169 write (iout,*) "strand",is1,is2," added to the array."
3172 istrand(1,nstrand)=is1
3173 istrand(2,nstrand)=is2
3175 end subroutine add_strand
3176 !------------------------------------------------------------------------------
3177 subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
3179 use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup
3180 use energy_data, only:itype,maxcont,molnum
3181 use compare_data, only:bfrag,hfrag,nbfrag,nhfrag
3182 use compare, only:freeres
3183 ! implicit real*8 (a-h,o-z)
3184 ! include 'DIMENSIONS'
3185 ! include 'DIMENSIONS.ZSCOPT'
3186 ! include 'COMMON.IOUNITS'
3187 ! include 'COMMON.FRAG'
3188 ! include 'COMMON.VAR'
3189 ! include 'COMMON.GEO'
3190 ! include 'COMMON.CHAIN'
3191 ! include 'COMMON.NAMES'
3192 ! include 'COMMON.INTERACT'
3193 integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres+1),&
3195 logical :: lprint,lprint_sec,not_done !el,freeres
3196 integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
3197 integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum
3198 real(kind=8) :: p1,p2
3199 !rel external freeres
3200 character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec))
3202 write (iout,*) "entered secondary2",ncont
3203 write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup
3205 write (iout,*) icont(1,i),icont(2,i)
3219 ! finding parallel beta
3220 !d write (iout,*) '------- looking for parallel beta -----------'
3226 if (i1.ge.nstart_sup .and. i1.le.nend_sup &
3227 .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
3228 !d write (iout,*) "parallel",i1,j1
3229 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
3232 !d write (iout,*) i1,j1
3238 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
3239 freeres(i1,j1,nsec,isec)) goto 5
3243 !d write (iout,*) i1,j1,not_done
3247 if (i1-ii1.gt.1) then
3251 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
3255 bfrag(1,nbfrag)=ii1+1
3256 bfrag(2,nbfrag)=i1+1
3257 bfrag(3,nbfrag)=jj1+1
3258 bfrag(4,nbfrag)=min0(j1+1,nres)
3262 isec(ij,nsec(ij))=nbeta
3266 isec(ij,nsec(ij))=nbeta
3271 if (nbeta.le.9) then
3272 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3273 "DefPropRes 'strand",nstrand,&
3274 "' 'num = ",ii1-1,"..",i1-1,"'"
3276 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3277 "DefPropRes 'strand",nstrand,&
3278 "' 'num = ",ii1-1,"..",i1-1,"'"
3281 if (nbeta.le.9) then
3282 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3283 "DefPropRes 'strand",nstrand,&
3284 "' 'num = ",jj1-1,"..",j1-1,"'"
3286 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3287 "DefPropRes 'strand",nstrand,&
3288 "' 'num = ",jj1-1,"..",j1-1,"'"
3290 write(12,'(a8,4i4)') &
3291 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
3295 endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
3298 ! finding antiparallel beta
3299 !d write (iout,*) '--------- looking for antiparallel beta ---------'
3304 if (freeres(i1,j1,nsec,isec)) then
3307 !d write (iout,*) i1,j1
3314 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
3315 freeres(i1,j1,nsec,isec)) goto 6
3319 !d write (iout,*) i1,j1,not_done
3323 if (i1-ii1.gt.1) then
3327 bfrag(2,nbfrag)=min0(i1+1,nres)
3328 bfrag(3,nbfrag)=min0(jj1+1,nres)
3335 if (nsec(ij).le.2) then
3336 isec(ij,nsec(ij))=nbeta
3342 if (nsec(ij).le.2) then
3343 isec(ij,nsec(ij))=nbeta
3348 if (lprint_sec) then
3349 write (iout,'(a,i3,4i4)')'antiparallel beta',&
3350 nbeta,ii1-1,i1,jj1,j1-1
3352 if (nstrand.le.9) then
3353 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3354 "DefPropRes 'strand",nstrand,&
3355 "' 'num = ",ii1-2,"..",i1-1,"'"
3357 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3358 "DefPropRes 'strand",nstrand,&
3359 "' 'num = ",ii1-2,"..",i1-1,"'"
3362 if (nstrand.le.9) then
3363 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3364 "DefPropRes 'strand",nstrand,&
3365 "' 'num = ",j1-2,"..",jj1-1,"'"
3367 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3368 "DefPropRes 'strand",nstrand,&
3369 "' 'num = ",j1-2,"..",jj1-1,"'"
3371 write(12,'(a8,4i4)') &
3372 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
3378 !d write (iout,*) "After beta:",nbfrag
3380 !d write (iout,*) (bfrag(j,i),j=1,4)
3383 if (nstrand.gt.0.and.lprint_sec) then
3384 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
3387 write(12,'(a9,i1,$)') " | strand",i
3389 write(12,'(a9,i2,$)') " | strand",i
3392 write(12,'(a1)') "'"
3396 ! finding alpha or 310 helix
3402 p1=phi(i1+2)*rad2deg
3404 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
3407 if (j1.eq.i1+3 .and. &
3408 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
3409 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
3410 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
3411 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
3414 if (nsec(ii1).eq.0) then
3423 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
3427 p1=phi(i1+2)*rad2deg
3428 p2=phi(j1+2)*rad2deg
3429 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
3432 !d write (iout,*) i1,j1,not_done,p1,p2
3435 if (j1-ii1.gt.4) then
3437 !d write (iout,*)'helix',nhelix,ii1,j1
3446 if (lprint_sec) then
3447 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
3448 if (nhelix.le.9) then
3449 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
3450 "DefPropRes 'helix",nhelix,&
3451 "' 'num = ",ii1-1,"..",j1-2,"'"
3453 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
3454 "DefPropRes 'helix",nhelix,&
3455 "' 'num = ",ii1-1,"..",j1-2,"'"
3462 if (nhelix.gt.0.and.lprint_sec) then
3463 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
3465 if (nhelix.le.9) then
3466 write(12,'(a8,i1,$)') " | helix",i
3468 write(12,'(a8,i2,$)') " | helix",i
3471 write(12,'(a1)') "'"
3474 if (lprint_sec) then
3475 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
3476 write(12,'(a20)') "XMacStand ribbon.mac"
3481 write(iout,*) 'UNRES seq:',anatemp
3483 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
3487 write(iout,*) 'helix ',(hfrag(i,j),i=1,2),anatemp
3493 do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j))
3496 do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j))
3501 do k=hfrag(1,j),hfrag(2,j)
3507 write (iout,*) "Secondary structure"
3513 write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
3514 write (iout,'(80a1)') (onelet(itype(k,mnum)),k=ist,ien)
3515 write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien)
3520 end subroutine secondary2
3521 !-------------------------------------------------
3522 ! logical function freeres(i,j,nsec,isec)
3523 ! include 'DIMENSIONS'
3524 ! integer :: isec(nres,4),nsec(nres)
3525 ! integer :: i,j,k,l
3528 ! if (nsec(i).gt.1.or.nsec(j).gt.1) return
3531 ! if (isec(i,k).eq.isec(j,l)) return
3536 ! end function freeres
3537 !-------------------------------------------------
3538 subroutine alloc_compar_arrays(nfrg,nlev)
3540 use energy_data, only:maxcont
3542 integer :: nfrg,nlev
3544 !write(iout,*) "in alloc conpar arrays: nlevel=", nlevel," nfrag(1)=",nfrag(1)
3545 !------------------------
3548 allocate(nsccont_frag_ref(mmaxfrag)) !(mmaxfrag) !wham
3549 allocate(isccont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag) !wham
3550 !------------------------
3553 allocate(rmsfrag(nfrg,nlev+1),nc_fragm(nfrg,nlev+1)) !(maxfrag,maxlevel)
3554 allocate(qfrag(nfrg,2)) !(maxfrag,2)
3555 allocate(rmscutfrag(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3556 allocate(ang_cut(nfrg),ang_cut1(nfrg),frac_min(nfrg)) !(maxfrag)
3557 allocate(nc_req_setf(nfrg,nlev+1),npiece(nfrg,nlev+1),&
3558 ielecont(nfrg,nlev+1),isccont(nfrg,nlev+1),irms(nfrg,nlev+1),&
3559 ishifft(nfrg,nlev+1),len_frag(nfrg,nlev+1)) !(maxfrag,maxlevel)
3560 allocate(ncont_nat(2,nfrg,nlev+1))
3561 allocate(n_shift(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3562 ! allocate(nfrag(nlev)) !(maxlevel)
3563 allocate(isnfrag(nlev+2)) !(maxlevel+1)
3564 allocate(ifrag(2,maxpiece,nfrg)) !(2,maxpiece,maxfrag)
3565 allocate(ipiece(maxpiece,nfrg,2:nlev+1)) !(maxpiece,maxfrag,2:maxlevel)
3566 allocate(istruct(nfrg),iloc(nfrg),nlist_frag(nfrg)) !(maxfrag)
3567 allocate(iclass(nlev*nfrg,nlev+1)) !(maxlevel*maxfrag,maxlevel)
3568 allocate(list_frag(nres,nfrg)) !(maxres,maxfrag)
3569 !------------------------
3572 ! integer,dimension(:,:),allocatable :: icont_pept_ref !(2,maxcont)
3573 allocate(ncont_frag_ref(mmaxfrag)) !(mmaxfrag)
3574 allocate(icont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag)
3575 ! integer,dimension(:),allocatable :: isec_ref !(maxres)
3576 !------------------------
3577 ! module w_comm_local
3579 allocate(creff(3,2*nres),cc(3,2*nres)) !(3,nres*2)
3580 allocate(iadded(nres)) !(nres)
3581 allocate(inumber(2,nres)) !(2,nres)
3584 !-------------------------------------------------------------------------------
3585 end subroutine alloc_compar_arrays
3587 !-------------------------------------------------------------------------------
3588 end module conform_compar