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,lprn,print_class)
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
68 call angnorm12(rmsang)
69 ! Level 1: check secondary and supersecondary structure
70 call elecont(lprn,ncont,icont,nnt,nct)
72 write (iout,*) "elecont finished"
75 call secondary2(lprn,.false.,ncont,icont,isecstr)
77 write (iout,*) "secondary2 finished"
80 call contact(lprn,ncontsc,icontsc,nnt,nct)
82 write(iout,*) "Assigning electrostatic contacts"
85 call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,&
88 write(iout,*) "Assigning sidechain contacts"
91 call contacts_between_fragments(lprn,3,ncontsc,icontsc,&
92 nsccont_frag,isccont_frag)
94 write(iout,*) "--> After contacts_between_fragments"
98 do j=1,isnfrag(nlevel+1)
105 write (iout,'(80(1h=))')
106 write (iout,*) "Level",1," fragment",j
107 write (iout,'(80(1h=))')
110 rmsfrag(j,1)=rmscalc(0,1,j,jcon,lprn)
111 ! Compare electrostatic contacts in the current conf with that in the native
113 if (lprn) write (iout,*) &
114 "Comparing electrostatic contact map and local structure"
116 ncnat=ncont_frag_ref(ind)
117 ! write (iout,*) "before match_contact:",nc_fragm(j,1),
120 call match_secondary(j,isecstr,nsec_match,lprn)
121 if (lprn) write (iout,*) "Fragment",j," nsec_match",&
122 nsec_match," length",len_frag(j,1)," min_len",&
123 frac_sec*len_frag(j,1)
124 if (nsec_match.lt.frac_sec*len_frag(j,1)) then
126 if (lprn) write (iout,*) "Fragment",j,&
127 " has incorrect secondary structure"
130 if (lprn) write (iout,*) "Fragment",j,&
131 " has correct secondary structure"
133 if (ielecont(j,1).gt.0) then
134 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
135 ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
136 ncont_frag(ind),icont_frag(1,1,ind),&
137 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
138 nc_req_setf(j,1),istruct(j),.true.,lprn)
139 else if (isccont(j,1).gt.0) then
140 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
141 nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
142 nsccont_frag(ind),isccont_frag(1,1,ind),&
143 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
144 nc_req_setf(j,1),istruct(j),.true.,lprn)
145 else if (iloc(j).gt.0) then
146 ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
147 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
148 0,icont_frag_ref(1,1,ind),&
149 ncont_frag(ind),icont_frag(1,1,ind),&
150 j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
151 0,istruct(j),.true.,lprn)
152 ! write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
157 if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
159 qfrag(j,1)=qwolynes(1,j)
160 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
161 if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match
162 ! write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1)
163 if (irms(j,1).gt.0) then
164 if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then
171 do while (rms.gt.rmscutfrag(1,j,1) .and. &
172 ishiff.lt.n_shift(1,j,1))
174 rms=rmscalc(-ishiff,1,j,jcon,lprn)
175 ! write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
176 ! & " rms",rms," rmscut",rmscutfrag(1,j,1)
177 if (lprn) write (iout,*) "rms",rmsfrag(j,1)
178 if (rms.gt.rmscutfrag(1,j,1)) then
179 rms=rmscalc(ishiff,1,j,jcon,lprn)
180 ! write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
183 if (lprn) write (iout,*) "rms",rmsfrag(j,1)
185 ! write (iout,*) "After loop: rms",rms,
186 ! & " rmscut",rmscutfrag(1,j,1)
187 ! write (iout,*) "iclass_rms",iclass_rms
188 if (rms.le.rmscutfrag(1,j,1)) then
193 ! write (iout,*) "iclass_rms",iclass_rms
195 ! write (iout,*) "ishif",ishif
196 if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
200 ! write (iout,*) "ishif",ishif," iclass",iclass(j,1),
201 ! & " iclass_rms",iclass_rms
202 if (nc_match.gt.0 .and. iclass_rms.gt.0) then
204 iclass(j,1)=iclass(j,1)+6
206 iclass(j,1)=iclass(j,1)+2
209 ncont_nat(1,j,1)=nc_match
210 ncont_nat(2,j,1)=ncon_match
212 ! write (iout,*) "iclass",iclass(j,1)
214 ! Next levels: Check arrangements of elementary fragments.
217 if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i))
219 write (iout,'(80(1h=))')
220 write (iout,*) "Level",i," fragment",j
221 write (iout,'(80(1h=))')
223 ! If an elementary fragment doesn't exist, don't check higher hierarchy levels.
226 if (iclass(ik,1).eq.0) then
231 if (i.eq.2 .and. ielecont(j,i).gt.0) then
234 if (lprn) write (iout,*) &
235 "Comparing electrostatic contact map: fragments",&
236 ipiece(1,j,i),ipiece(2,j,i)," ind",ind
237 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
238 ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
239 ncont_frag(ind),icont_frag(1,1,ind),&
240 j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
241 nc_req_setf(j,i),2,.false.,lprn)
243 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
244 if (nc_match.gt.0) then
251 ncont_nat(1,j,i)=nc_match
252 ncont_nat(2,j,i)=ncon_match
254 else if (i.eq.2 .and. isccont(j,i).gt.0) then
257 if (lprn) write (iout,*) &
258 "Comparing sidechain contact map: fragments",&
259 ipiece(1,j,i),ipiece(2,j,i)," ind",ind
260 call match_contact(ishif1,ishif2,nc_match,ncon_match,&
261 nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
262 nsccont_frag(ind),isccont_frag(1,1,ind),&
263 j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
264 nc_req_setf(j,i),2,.false.,lprn)
266 if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
267 if (nc_match.gt.0) then
274 ncont_nat(1,j,i)=nc_match
275 ncont_nat(2,j,i)=ncon_match
277 else if (i.eq.2) then
281 if (i.eq.2) qfrag(j,2)=qwolynes(2,j)
282 if (lprn) write (iout,*) &
283 "Comparing rms: fragments",&
284 (ipiece(k,j,i),k=1,npiece(j,i))
285 rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn)
286 if (irms(j,i).gt.0) then
289 if (lprn) write (iout,*) "rms",rmsfrag(j,i)
290 ! write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i),
291 ! & " rmscutfrag",rmscutfrag(1,j,i)
292 if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then
298 do while (rms.gt.rmscutfrag(1,j,i) .and. &
299 ishif.lt.n_shift(1,j,i))
301 rms=rmscalc(-ishif,i,j,jcon,lprn)
302 ! print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
303 if (lprn) write (iout,*) "rms",rmsfrag(j,i)
304 if (rms.gt.rmscutfrag(1,j,i)) then
305 rms=rmscalc(ishif,i,j,jcon,lprn)
306 ! print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
308 if (lprn) write (iout,*) "rms",rms
310 if (rms.le.rmscutfrag(1,j,i)) then
317 if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and. &
318 isccont(j,i).eq.0 ) then
319 write (iout,*) "Error: no measure of comparison specified:",&
324 write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
326 iclass(j,i) = min0(iclass_con,iclass_rms)
327 if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then
328 ishifft(j,i)=ishifft_rms
330 ishifft(j,i)=ishifft_con
332 else if (i.gt.2) then
333 iclass(j,i) = iclass_rms
334 ishifft(j,i)= ishifft_rms
341 ! Compute the structural class
343 IF (.NOT. BINARY) THEN
351 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j
355 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
356 ! & " iclass",iclass(j,i)," im",im
362 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j
364 if (iclass(j,i).gt.0) then
369 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
370 ! & " iclass",iclass(j,i)," im",im
374 idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j
376 if (iclass(j,i).gt.1) then
381 ! write (iout,*) "i",i," j",j," idig",idig," iex",iex,
382 ! & " iclass",iclass(j,i)," im",im
389 if (print_class) then
391 write(istat,'(i6,$)') jcon+indstart(me)-1
392 write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
395 write(istat,'(i6,$)') jcon
396 write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
399 write (istat,'(f8.3,2f6.3,$)') &
400 rms_nat,qnat,rmsang/(nres-3)
402 write(istat,'(1x,$,20(i3,$))') &
403 (ncont_nat(1,k,j),k=1,nfrag(j))
405 write(istat,'(1x,$,20(f5.1,f5.2$))') &
406 (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j))
408 write(istat,'(1x,$,20(f5.1$))') &
409 (rmsfrag(k,j),k=1,nfrag(j))
411 write(istat,'(1x,$,20(i1,$))') &
412 (iclass(k,j),k=1,nfrag(j))
415 write (istat,'(" ",$)')
417 write (istat,'(100(i1,$))')(iclass(k,j),&
419 if (j.lt.nlevel) write(iout,'(".",$)')
423 write (istat,'(i10)') iscore
427 END subroutine conf_compar
428 !-----------------------------------------------------------------------------
430 !-----------------------------------------------------------------------------
431 subroutine add_angpair(ici,icj,nang_pair,iang_pair)
433 ! implicit real*8 (a-h,o-z)
434 ! include 'DIMENSIONS'
435 ! include 'COMMON.IOUNITS'
436 ! include 'COMMON.CHAIN'
437 integer :: ici,icj,nang_pair,iang_pair(2,nres)
438 integer :: i,ian1,ian2
439 ! write (iout,*) "add_angpair: ici",ici," icj",icj,
440 ! & " nang_pair",nang_pair
442 if (ian1.lt.4 .or. ian1.gt.nres) return
444 ! write (iout,*) "ian1",ian1," ian2",ian2
445 if (ian2.lt.4 .or. ian2.gt.nres) return
447 if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
449 nang_pair=nang_pair+1
450 iang_pair(1,nang_pair)=ian1
451 iang_pair(2,nang_pair)=ian2
453 end subroutine add_angpair
454 !-------------------------------------------------------------------------
455 subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,lprn)
457 use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
459 ! implicit real*8 (a-h,o-z)
460 ! include 'DIMENSIONS'
461 ! include 'DIMENSIONS.ZSCOPT'
462 ! include 'DIMENSIONS.COMPAR'
463 ! include 'COMMON.IOUNITS'
464 ! include 'COMMON.VAR'
465 ! include 'COMMON.COMPAR'
466 ! include 'COMMON.CHAIN'
467 ! include 'COMMON.GEO'
468 real(kind=8) :: pinorm,deltang
470 integer :: jfrag,ishif1,ishif2,nn,npart,nn4,nne
471 real(kind=8) :: diffang_max,angn,fract,ff
472 integer :: i,j,nbeg,nend,ll,longest
473 if (lprn) write (iout,'(80(1h*))')
477 npart = npiece(jfrag,1)
479 nne = min0(nend_sup,nres)
480 if (lprn) write (iout,*) "nn4",nn4," nne",nne
482 nbeg = ifrag(1,i,jfrag) + 3 - ishif1
483 if (nbeg.lt.nn4) nbeg=nn4
484 nend = ifrag(2,i,jfrag) + 1 - ishif2
485 if (nend.gt.nne) nend=nne
486 if (nend.ge.nbeg) then
487 nn = nn + nend - nbeg + 1
488 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
489 " nn",nn," ishift1",ishif1," ishift2",ishif2
490 if (lprn) write (iout,*) "angles"
494 ! deltang = pinorm(phi(j)-phi_ref(j+ishif1))
495 deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),&
496 theta_ref(j+ishif1),phi(j),theta(j-1),theta(j))
497 if (dabs(deltang).gt.diffang_max) then
498 if (ll.gt.longest) longest = ll
503 if (ll.gt.longest) longest = ll
504 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
505 rad2deg*phi_ref(j+ishif1),rad2deg*deltang
506 angn=angn+dabs(deltang)
509 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
510 if (lprn) write (iout,*)"segment",i," longest fragment within",&
511 diffang_max*rad2deg,":",longest," fraction",ff
512 if (ff.lt.fract) fract = ff
520 if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,&
523 end subroutine angnorm
524 !-------------------------------------------------------------------------
525 subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,&
526 diffang_max,anorm,fract)
528 use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
530 ! implicit real*8 (a-h,o-z)
531 ! include 'DIMENSIONS'
532 ! include 'DIMENSIONS.ZSCOPT'
533 ! include 'DIMENSIONS.COMPAR'
534 ! include 'COMMON.IOUNITS'
535 ! include 'COMMON.VAR'
536 ! include 'COMMON.COMPAR'
537 ! include 'COMMON.CHAIN'
538 ! include 'COMMON.GEO'
539 integer :: ncont,icont(2,ncont),longest
540 real(kind=8) :: anorm,diffang_max,fract
541 integer :: npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
542 real(kind=8) :: pinorm
544 integer :: jfrag,ishif1,ishif2
545 integer :: nn,nn4,nne,npart,i,j,jstart,jend,ic1,ic2,idi,iic
546 integer :: nbeg,nend,ll
547 real(kind=8) :: angn,ishifc,deltang,ff
549 if (lprn) write (iout,'(80(1h*))')
551 ! Determine the segments for which angles will be compared
554 nne = min0(nend_sup,nres)
555 if (lprn) write (iout,*) "nn4",nn4," nne",nne
556 npart=npiece(jfrag,1)
559 ! write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
560 if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. &
561 icont(1,1).gt.ifrag(2,i,jfrag)) goto 11
563 do while (jstart.lt.ncont .and. &
564 icont(1,jstart).lt.ifrag(1,i,jfrag))
565 ! write (iout,*) "jstart",jstart," icont",icont(1,jstart),
566 ! & " ifrag",ifrag(1,i,jfrag)
569 ! write (iout,*) "jstart",jstart," icont",icont(1,jstart),
570 ! & " ifrag",ifrag(1,i,jfrag)
571 if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11
574 ifrag_c(1,npiece_c)=icont(1,jstart)
576 do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag))
577 ! write (iout,*) "jend",jend," icont",icont(1,jend),
578 ! & " ifrag",ifrag(2,i,jfrag)
581 ! write (iout,*) "jend",jend," icont",icont(1,jend),
582 ! & " ifrag",ifrag(2,i,jfrag)
584 ifrag_c(2,npiece_c)=icont(1,jend)+1
585 ishift_c(npiece_c)=ishif1
586 ! write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
587 ! & " ic1",ic1," ic2",ic2,
588 ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
590 if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
595 ! write (iout,*) "idi",idi
597 if (icont(2,1).gt.ifrag(2,i,jfrag) .or. &
598 icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12
600 do while (jstart.lt.ncont .and. &
601 icont(2,jstart).lt.ifrag(1,i,jfrag))
602 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
603 ! & " ifrag",ifrag(1,i,jfrag)
606 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
607 ! & " ifrag",ifrag(1,i,jfrag)
608 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
611 ifrag_c(2,npiece_c)=icont(2,jstart)+1
613 do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag))
614 ! write (iout,*) "jend",jend," icont",icont(2,jend),
615 ! & " ifrag",ifrag(2,i,jfrag)
618 ! write (iout,*) "jend",jend," icont",icont(2,jend),
619 ! & " ifrag",ifrag(2,i,jfrag)
620 else if (idi.eq.-1) then
621 if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or. &
622 icont(2,1).lt.ifrag(1,i,jfrag)) goto 12
624 do while (jstart.gt.ncont .and. &
625 icont(2,jstart).lt.ifrag(1,i,jfrag))
626 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
627 ! & " ifrag",ifrag(1,i,jfrag)
630 ! write (iout,*) "jstart",jstart," icont",icont(2,jstart),
631 ! & " ifrag",ifrag(1,i,jfrag)
632 if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
635 ifrag_c(2,npiece_c)=icont(2,jstart)+1
637 do while (jend.lt.ncont .and. &
638 icont(2,jend).gt.ifrag(2,i,jfrag))
639 ! write (iout,*) "jend",jend," icont",icont(2,jend),
640 ! & " ifrag",ifrag(2,i,jfrag)
643 ! write (iout,*) "jend",jend," icont",icont(2,jend),
644 ! & " ifrag",ifrag(2,i,jfrag)
652 ! write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
653 ! & " jstart:",jstart," jend",jend,
654 ! & " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
655 ifrag_c(1,npiece_c)=ic1
656 ifrag_c(2,npiece_c)=ic2+1
657 ishift_c(npiece_c)=ishif2
661 write (iout,*) "Before merge: npiece_c",npiece_c
663 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
667 ! Merge overlapping segments (e.g., avoid splitting helices)
670 do while (i .lt. npiece_c)
671 if (ishift_c(i).eq.ishift_c(i+1) .and. &
672 ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
673 ifrag_c(2,i)=ifrag_c(2,i+1)
675 ishift_c(j)=ishift_c(j+1)
676 ifrag_c(1,j)=ifrag_c(1,j+1)
677 ifrag_c(2,j)=ifrag_c(2,j+1)
685 write (iout,*) "After merge: npiece_c",npiece_c
687 write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
700 nbeg = ifrag_c(1,i) + 3 - ishifc
701 if (nbeg.lt.nn4) nbeg=nn4
702 nend = ifrag_c(2,i) - ishifc + 1
703 if (nend.gt.nne) nend=nne
704 if (nend.ge.nbeg) then
705 nn = nn + nend - nbeg + 1
706 if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
707 " nn",nn," ishifc",ishifc
708 if (lprn) write (iout,*) "angles"
712 ! deltang = pinorm(phi(j)-phi_ref(j+ishifc))
713 deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),&
714 theta_ref(j+ishifc),phi(j),theta(j-1),theta(j))
715 if (dabs(deltang).gt.diffang_max) then
716 if (ll.gt.longest) longest = ll
721 if (ll.gt.longest) longest = ll
722 if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
723 rad2deg*phi_ref(j+ishifc),rad2deg*deltang
724 angn=angn+dabs(deltang)
727 ff = dfloat(longest)/dfloat(nend - nbeg + 4)
728 if (lprn) write (iout,*)"segment",i," longest fragment within",&
729 diffang_max*rad2deg,":",longest," fraction",ff
730 if (ff.lt.fract) fract = ff
733 if (nn.gt.0) anorm = angn/nn
734 if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
736 end subroutine angnorm2
737 !-------------------------------------------------------------------------
738 real(kind=8) function angnorm1(nang_pair,iang_pair,lprn)
740 use geometry_data, only:phi,theta,rad2deg
741 ! implicit real*8 (a-h,o-z)
742 ! include 'DIMENSIONS'
743 ! include 'DIMENSIONS.ZSCOPT'
744 ! include 'DIMENSIONS.COMPAR'
745 ! include 'COMMON.IOUNITS'
746 ! include 'COMMON.VAR'
747 ! include 'COMMON.COMPAR'
748 ! include 'COMMON.CHAIN'
749 ! include 'COMMON.GEO'
751 integer :: nang_pair,iang_pair(2,nres)
752 real(kind=8) :: pinorm
754 real(kind=8) :: angn,deltang
756 if (lprn) write (iout,'(80(1h*))')
757 if (lprn) write (iout,*) "nang_pair",nang_pair
758 if (lprn) write (iout,*) "angles"
762 ! deltang = pinorm(phi(ia1)-phi_ref(ia2))
763 deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),&
764 theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2))
765 if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),&
766 rad2deg*phi_ref(ia2),rad2deg*deltang
767 angn=angn+dabs(deltang)
770 write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
771 angnorm1 = angn/nang_pair
773 end function angnorm1
774 !------------------------------------------------------------------------------
775 subroutine angnorm12(diff)
777 use geometry_data, only:phi,theta,nstart_sup,nend_sup
778 ! implicit real*8 (a-h,o-z)
779 ! include 'DIMENSIONS'
780 ! include 'DIMENSIONS.ZSCOPT'
781 ! include 'DIMENSIONS.COMPAR'
782 ! include 'COMMON.IOUNITS'
783 ! include 'COMMON.VAR'
784 ! include 'COMMON.COMPAR'
785 ! include 'COMMON.CHAIN'
786 ! include 'COMMON.GEO'
787 real(kind=8) :: pinorm,diff
791 nne = min0(nend_sup,nres)
793 ! diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j)))
796 ! diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j)))
797 diff=diff+spherang(phi_ref(j),theta_ref(j-1),&
798 theta_ref(j),phi(j),theta(j-1),theta(j))
801 end subroutine angnorm12
802 !--------------------------------------------------------------------------------
803 real(kind=8) function spherang(gam1,theta11,theta12,&
804 gam2,theta21,theta22)
806 use geometry, only:arcos
807 real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,&
808 x1,x2,xmed,f1,f2,fmed
809 real(kind=8) :: tolx=1.0d-4, tolf=1.0d-4
810 real(kind=8) :: sumcos
811 !el real(kind=8) :: pinorm,sumangp !arcos,
812 integer :: it,maxit=100
813 ! Calculate the difference of the angles of two superposed 4-redidue fragments
821 ! The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
822 ! to achieve the minimum difference between the O'-C-O and P-C-P angles;
823 ! the sum of these angles is the difference returned by the function.
826 ! If thetas match, take the difference of gamma and exit.
827 if (dabs(theta11-theta12).lt.tolx &
828 .and. dabs(theta21-theta22).lt.tolx) then
829 spherang=dabs(pinorm(gam2-gam1))
832 ! If the gammas are the same, take the difference of thetas and exit.
834 x2=0.5d0*pinorm(gam2-gam1)
835 if (dabs(x2) .lt. tolx) then
836 spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
838 else if (x2.lt.0.0d0) then
842 ! Else apply regula falsi method to compute optimum overlap of the terminal Calphas
843 f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
844 f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
846 xmed=x1-f1*(x2-x1)/(f2-f1)
847 fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
848 ! write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
849 if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx) &
850 .and. dabs(fmed).lt.tolf ) then
854 else if ( fmed*f1.lt.0.0d0 ) then
863 spherang=arcos(dcos(theta11)*dcos(theta12) &
864 +dsin(theta11)*dsin(theta12)*dcos(x1))+ &
865 arcos(dcos(theta21)*dcos(theta22)+ &
866 dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
868 end function spherang
869 !--------------------------------------------------------------------------------
870 real(kind=8) function sumangp(gam1,theta11,theta12,gam2,&
873 real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,fi,&
874 cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,&
876 ! derivarive of the sum of the difference of the angles of a 4-residue fragment.
877 ! real(kind=8) :: arcos
886 cosd1=cost11*cost12+sint11*sint12*dcos(fi)
887 cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
888 sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1) &
889 +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)
892 !-----------------------------------------------------------------------------
894 !-----------------------------------------------------------------------------
895 subroutine contact(lprint,ncont,icont,ist,ien)
898 use geometry_data, only:c,dc,dc_norm
899 use energy_data, only:itype,maxcont
901 ! include 'DIMENSIONS'
902 ! include 'DIMENSIONS.ZSCOPT'
903 ! include 'COMMON.CONTROL'
904 ! include 'COMMON.IOUNITS'
905 ! include 'COMMON.CHAIN'
906 ! include 'COMMON.INTERACT'
907 ! include 'COMMON.FFIELD'
908 ! include 'COMMON.NAMES'
909 ! include 'COMMON.CALC'
910 ! include 'COMMON.CONTPAR'
911 ! include 'COMMON.LOCAL'
912 integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
913 real(kind=8) :: csc !el,dist
914 real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,&
917 integer,dimension(2,maxcont) :: icont
918 real(kind=8) :: u,v,a(3),b(3),dla,dlb
928 write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
929 c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),&
930 dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
933 110 format (a,'(',i3,')',9f8.3)
936 if (iti.le.0 .or. iti.gt.ntyp) cycle
939 if (itj.le.0 .or. itj.gt.ntyp) cycle
942 xj = c(1,nres+j)-c(1,nres+i)
943 yj = c(2,nres+j)-c(2,nres+i)
944 zj = c(3,nres+j)-c(3,nres+i)
945 dxi = dc_norm(1,nres+i)
946 dyi = dc_norm(2,nres+i)
947 dzi = dc_norm(3,nres+i)
948 dxj = dc_norm(1,nres+j)
949 dyj = dc_norm(2,nres+j)
950 dzj = dc_norm(3,nres+j)
955 ! write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
956 if (icomparfunc.eq.1) then
957 call contfunc(csc,iti,itj)
958 else if (icomparfunc.eq.2) then
959 call scdist(csc,iti,itj)
960 else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
961 csc = dist(nres+i,nres+j)
962 else if (icomparfunc.eq.4) then
963 call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
965 write (*,*) "Error - Unknown sidechain contact function"
966 write (iout,*) "Error - Unknown sidechain contact function"
968 if (csc.lt.sc_cutoff(iti,itj)) then
969 ! write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
970 ! & " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
971 ! & dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
972 ! & dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
974 ! write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
975 ! & sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
976 ! & chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
985 ddsc(ncont)=1.0d0/rij
992 write (iout,'(a)') 'Contact map:'
998 write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') &
999 i,restyp(it1),i1,restyp(it2),i2,cscore(i),&
1000 sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),&
1001 omt1(i),omt2(i),omt12(i)
1005 end subroutine contact
1007 !----------------------------------------------------------------------------
1008 subroutine contact(lprint,ncont,icont)
1010 use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii
1011 ! include 'DIMENSIONS'
1012 ! include 'COMMON.IOUNITS'
1013 ! include 'COMMON.CHAIN'
1014 ! include 'COMMON.INTERACT'
1015 ! include 'COMMON.FFIELD'
1016 ! include 'COMMON.NAMES'
1017 real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
1018 integer :: ncont,icont(2,maxcont)
1020 integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
1021 real(kind=8) :: rcomp
1024 ! print *,'nnt=',nnt,' nct=',nct
1030 ! rcomp=sigmaii(iti,itj)+1.0D0
1031 rcomp=facont*sigmaii(iti,itj)
1033 ! rcomp=sigma(iti,itj)+1.0D0
1034 rcomp=facont*sigma(iti,itj)
1037 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
1038 if (dist(nres+i,nres+j).lt.rcomp) then
1046 write (iout,'(a)') 'Contact map:'
1052 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1053 i,restyp(it1),i1,restyp(it2),i2
1057 end subroutine contact
1059 !----------------------------------------------------------------------------
1060 real(kind=8) function contact_fract(ncont,ncont_ref,&
1063 use energy_data, only:maxcont
1065 ! include 'DIMENSIONS'
1066 ! include 'COMMON.IOUNITS'
1067 integer :: i,j,nmatch
1068 integer :: ncont,ncont_ref
1069 integer,dimension(2,maxcont) :: icont,icont_ref
1071 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
1072 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
1073 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
1074 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
1075 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
1078 if (icont(1,i).eq.icont_ref(1,j) .and. &
1079 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
1082 ! print *,' nmatch=',nmatch
1083 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
1084 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
1086 end function contact_fract
1088 !------------------------------------------------------------------------------
1089 subroutine pept_cont(lprint,ncont,icont)
1091 use geometry_data, only:c
1092 use energy_data, only:maxcont,nnt,nct,itype
1094 ! include 'DIMENSIONS'
1095 ! include 'DIMENSIONS.ZSCOPT'
1096 ! include 'COMMON.IOUNITS'
1097 ! include 'COMMON.CHAIN'
1098 ! include 'COMMON.INTERACT'
1099 ! include 'COMMON.FFIELD'
1100 ! include 'COMMON.NAMES'
1101 integer :: ncont,icont(2,maxcont)
1102 integer :: i,j,k,kkk,i1,i2,it1,it2
1104 !el real(kind=8) :: dist
1105 real(kind=8) :: rcomp=5.5d0
1108 print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
1111 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
1115 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
1117 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
1125 write (iout,'(a)') 'PP contact map:'
1131 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1132 i,restyp(it1),i1,restyp(it2),i2
1136 end subroutine pept_cont
1137 !-----------------------------------------------------------------------------
1139 !-----------------------------------------------------------------------------
1140 subroutine contacts_between_fragments(lprint,is,ncont,icont,&
1141 ncont_interfrag,icont_interfrag)
1143 use energy_data, only:itype,maxcont
1144 ! include 'DIMENSIONS'
1145 ! include 'DIMENSIONS.ZSCOPT'
1146 ! include 'DIMENSIONS.COMPAR'
1147 ! include 'COMMON.INTERACT'
1148 ! include 'COMMON.COMPAR'
1149 ! include 'COMMON.IOUNITS'
1150 ! include 'COMMON.CHAIN'
1151 ! include 'COMMON.NAMES'
1152 integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),&
1153 icont_interfrag(2,maxcont,mmaxfrag)
1154 logical :: OK1,OK2,lprint
1155 integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2
1156 ! Determine the contacts that occur within a fragment and between fragments.
1161 ! write (iout,*) "i",i,(ifrag(1,k,i),ifrag(2,k,i)
1162 ! & ,k=1,npiece(i,1))
1163 ! write (iout,*) "j",j,(ifrag(1,k,j),ifrag(2,k,j)
1164 ! & ,k=1,npiece(j,1))
1165 ! write (iout,*) "ncont",ncont
1171 do while (.not.OK1 .and. l.lt.npiece(j,1))
1173 OK1=ic1.ge.ifrag(1,l,j)-is .and. &
1174 ic1.le.ifrag(2,l,j)+is
1178 do while (.not.OK2 .and. l.lt.npiece(i,1))
1180 OK2=ic2.ge.ifrag(1,l,i)-is .and. &
1181 ic2.le.ifrag(2,l,i)+is
1183 ! write(iout,*) "k",k," ic1",ic1," ic2",ic2," OK1",OK1,
1185 if (OK1.and.OK2) then
1187 icont_interfrag(1,nc,ind)=ic1
1188 icont_interfrag(2,nc,ind)=ic2
1189 ! write (iout,*) "nc",nc," ic1",ic1," ic2",ic2
1192 ncont_interfrag(ind)=nc
1193 ! do k=1,ncont_interfrag(ind)
1194 ! i1=icont_interfrag(1,k,ind)
1195 ! i2=icont_interfrag(2,k,ind)
1198 ! write (iout,'(i3,2x,a,i4,2x,a,i4)')
1199 ! & i,restyp(it1),i1,restyp(it2),i2
1204 write (iout,*) "Contacts within fragments:"
1206 write (iout,*) "Fragment",i," (",(ifrag(1,k,i),&
1207 ifrag(2,k,i),k=1,npiece(i,1)),")"
1209 do k=1,ncont_interfrag(ind)
1210 i1=icont_interfrag(1,k,ind)
1211 i2=icont_interfrag(2,k,ind)
1214 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1215 i,restyp(it1),i1,restyp(it2),i2
1219 write (iout,*) "Contacts between fragments:"
1223 write (iout,*) "Fragments",i," (",(ifrag(1,k,i),&
1224 ifrag(2,k,i),k=1,npiece(i,1)),") and",j," (",&
1225 (ifrag(1,k,j),ifrag(2,k,j),k=1,npiece(j,1)),")"
1226 write (iout,*) "Number of contacts",&
1227 ncont_interfrag(ind)
1229 do k=1,ncont_interfrag(ind)
1230 i1=icont_interfrag(1,k,ind)
1231 i2=icont_interfrag(2,k,ind)
1234 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1235 i,restyp(it1),i1,restyp(it2),i2
1241 end subroutine contacts_between_fragments
1242 !-----------------------------------------------------------------------------
1244 !-----------------------------------------------------------------------------
1245 subroutine contfunc(cscore,itypi,itypj)
1247 ! This subroutine calculates the contact function based on
1248 ! the Gay-Berne potential of interaction.
1251 ! implicit real*8 (a-h,o-z)
1252 ! include 'DIMENSIONS'
1253 ! include 'COMMON.CONTPAR'
1254 ! include 'COMMON.CALC'
1256 integer :: itypi,itypj
1257 real(kind=8) :: cscore,sig0ij,rrij,sig,rij_shift,evdw,e2
1259 sig0ij=sig_comp(itypi,itypj)
1260 chi1=chi_comp(itypi,itypj)
1261 chi2=chi_comp(itypj,itypi)
1263 chip1=chip_comp(itypi,itypj)
1264 chip2=chip_comp(itypj,itypi)
1266 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1268 ! Calculate angle-dependent terms of the contact function
1272 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1273 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1274 om12=dxi*dxj+dyi*dyj+dzi*dzj
1276 ! print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
1278 ! & rij,rrij,om1,om2,om12
1279 ! Calculate eps1(om12)
1280 faceps1=1.0D0-om12*chiom12
1281 faceps1_inv=1.0D0/faceps1
1282 eps1=dsqrt(faceps1_inv)
1283 ! Following variable is eps1*deps1/dom12
1284 eps1_om12=faceps1_inv*chiom12
1285 ! Calculate sigma(om1,om2,om12)
1289 facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1290 sigsq=1.0D0-facsig*faceps1_inv
1291 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1294 chipom12=chip12*om12
1295 facp=1.0D0-om12*chipom12
1297 facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1298 ! Following variable is the square root of eps2
1299 eps2rt=1.0D0-facp1*facp_inv
1301 sig=sig0ij*dsqrt(sigsq)
1302 rij_shift=1.0D0/rij-sig+sig0ij
1303 if (rij_shift.le.0.0D0) then
1305 cscore = -dlog(evdw+1.0d-6)
1308 rij_shift=1.0D0/rij_shift
1309 e2=(rij_shift*sig0ij)**expon
1310 evdw=dabs(eps1*eps2rt**2*e2)
1311 if (evdw.gt.1.0d1) evdw = 1.0d1
1312 cscore = -dlog(evdw+1.0d-6)
1314 end subroutine contfunc
1315 !------------------------------------------------------------------------------
1316 subroutine scdist(cscore,itypi,itypj)
1318 ! This subroutine calculates the contact distance
1321 ! implicit real*8 (a-h,o-z)
1322 ! include 'DIMENSIONS'
1323 ! include 'COMMON.CONTPAR'
1324 ! include 'COMMON.CALC'
1325 integer :: itypi,itypj
1326 real(kind=8) :: cscore,rrij
1328 chi1=chi_comp(itypi,itypj)
1329 chi2=chi_comp(itypj,itypi)
1331 rrij=xj*xj+yj*yj+zj*zj
1333 ! Calculate angle-dependent terms of the contact function
1337 om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1338 om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1339 om12=dxi*dxj+dyi*dyj+dzi*dzj
1344 cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)
1346 end subroutine scdist
1347 !------------------------------------------------------------------------------
1349 !------------------------------------------------------------------------------
1350 subroutine elecont(lprint,ncont,icont,ist,ien)
1352 use geometry_data, only:c
1353 use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2
1355 ! include 'DIMENSIONS'
1356 ! include 'DIMENSIONS.ZSCOPT'
1357 ! include 'DIMENSIONS.COMPAR'
1358 ! include 'COMMON.IOUNITS'
1359 ! include 'COMMON.CHAIN'
1360 ! include 'COMMON.INTERACT'
1361 ! include 'COMMON.FFIELD'
1362 ! include 'COMMON.NAMES'
1363 ! include 'COMMON.LOCAL'
1365 integer :: i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
1366 real(kind=8) :: rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,&
1367 xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,&
1368 vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,&
1370 real(kind=8),dimension(2,2) :: elpp6c=reshape((/-0.2379d0,&
1371 -0.2056d0,-0.2056d0,-0.0610d0/),shape(elpp6c))
1372 real(kind=8),dimension(2,2) :: elpp3c=reshape((/ 0.0503d0,&
1373 0.0000d0, 0.0000d0, 0.0692d0/),shape(elpp3c))
1374 real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc
1375 real(kind=8) :: elcutoff=-0.3d0
1376 real(kind=8) :: elecutoff_14=-0.5d0
1377 integer :: ncont,icont(2,maxcont)
1378 real(kind=8) :: econt(maxcont)
1380 ! Load the constants of peptide bond - peptide bond interactions.
1381 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
1382 ! proline) - determined by averaging ECEPP energy.
1386 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
1387 ! data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
1388 !el data (elpp6c) /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
1389 !el data (elpp3c) / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
1390 !el data (elcutoff) /-0.3d0/
1391 !el data (elecutoff_14) /-0.5d0/
1394 if (lprint) write (iout,'(a)') &
1395 "Constants of electrostatic interaction energy expression."
1399 appc(i,j)=epp(i,j)*rri*rri
1400 bppc(i,j)=-2.0*epp(i,j)*rri
1401 ael6c(i,j)=elpp6c(i,j)*4.2**6
1402 ael3c(i,j)=elpp3c(i,j)*4.2**3
1404 write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),&
1423 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1424 if (iteli.eq.2 .and. itelj.eq.2 &
1425 .or.iteli.eq.0 .or.itelj.eq.0) goto 4
1426 aaa=appc(iteli,itelj)
1427 bbb=bppc(iteli,itelj)
1428 ael6i=ael6c(iteli,itelj)
1429 ael3i=ael3c(iteli,itelj)
1433 xj=c(1,j)+0.5*dxj-xmedi
1434 yj=c(2,j)+0.5*dyj-ymedi
1435 zj=c(3,j)+0.5*dzj-zmedi
1436 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
1441 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
1442 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
1443 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
1444 fac=cosa-3.0*cosb*cosg
1450 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
1453 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
1454 j.eq.i+2 .and. eesij.le.elecutoff_14) then
1465 write (iout,*) 'Total average electrostatic energy: ',ees
1466 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
1468 write (iout,*) 'Electrostatic contacts before pruning: '
1474 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1475 i,restyp(it1),i1,restyp(it2),i2,econt(i)
1478 ! For given residues keep only the contacts with the greatest energy.
1480 do while (i.lt.ncont)
1486 do while (j.lt.ncont)
1488 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
1489 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
1490 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
1491 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
1492 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
1493 if (ic1.eq.icont(1,j)) then
1495 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)&
1496 .and. iabs(icont(1,k)-ic1).le.2 .and. &
1497 econt(k).lt.econt(j) ) goto 21
1499 else if (ic2.eq.icont(2,j) ) then
1501 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)&
1502 .and. iabs(icont(2,k)-ic2).le.2 .and. &
1503 econt(k).lt.econt(j) ) goto 21
1506 ! Remove ith contact
1508 icont(1,k-1)=icont(1,k)
1509 icont(2,k-1)=icont(2,k)
1514 ! write (iout,*) "ncont",ncont
1516 ! write (iout,*) icont(1,k),icont(2,k)
1519 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
1521 if (ic1.eq.icont(1,j)) then
1523 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
1524 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
1525 econt(k).lt.econt(i) ) goto 21
1527 else if (ic2.eq.icont(2,j) ) then
1529 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
1530 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
1531 econt(k).lt.econt(i) ) goto 21
1534 ! Remove jth contact
1536 icont(1,k-1)=icont(1,k)
1537 icont(2,k-1)=icont(2,k)
1541 ! write (iout,*) "ncont",ncont
1543 ! write (iout,*) icont(1,k),icont(2,k)
1554 write (iout,*) 'Electrostatic contacts after pruning: '
1560 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1561 i,restyp(it1),i1,restyp(it2),i2,econt(i)
1565 end subroutine elecont
1566 !------------------------------------------------------------------------------
1568 !------------------------------------------------------------------------------
1569 subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,&
1570 ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,&
1571 nc_frac,nc_req_set,istr,llocal,lprn)
1573 use energy_data, only:maxcont
1574 ! implicit real*8 (a-h,o-z)
1575 ! include 'DIMENSIONS'
1576 ! include 'COMMON.IOUNITS'
1577 integer :: ncont_ref,ncont,ishift,ishif2,nc_match
1578 integer,dimension(2,maxcont) :: icont_ref,icont !(2,maxcont)
1579 real(kind=8) :: nc_frac
1580 logical :: llocal,lprn
1581 integer :: ishif1,nc_match1_max,jfrag,n_shif1,n_shif2,&
1582 nc_req_set,istr,nc_match_max
1583 integer :: i,nc_req,nc_match1,is,js
1586 nc_match_max=nc_match_max+ &
1587 min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
1591 else if (nc_req_set.eq.0) then
1592 nc_req=nc_match_max*nc_frac
1594 nc_req = dmin1(nc_match_max*nc_frac+0.5d0,&
1595 dfloat(nc_req_set)+1.0d-7)
1597 ! write (iout,*) "match_contact: nc_req:",nc_req
1598 ! write (iout,*) "nc_match_max",nc_match_max
1599 ! write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
1600 ! & " n_shif2",n_shif2
1601 ! Match current contact map against reference contact map; exit, if at least
1602 ! half of the contacts match
1603 call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,&
1604 ncont,icont,jfrag,llocal,lprn)
1605 nc_match1_max=nc_match1
1606 if (lprn .and. nc_match.gt.0) write (iout,*) &
1607 "Shift:",0,0," nc_match1",nc_match1,&
1608 " nc_match=",nc_match," req'd",nc_req
1609 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1610 nc_req.eq.0 .and. nc_match.eq.1) then
1615 ! If sufficient matches are not found, try to shift contact maps up to three
1617 if (n_shif1.gt.0) then
1619 ! The following four tries help to find shifted beta-sheet patterns
1620 ! Shift "left" strand backward
1621 call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,&
1622 icont_ref,ncont,icont,jfrag,llocal,lprn)
1623 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1624 if (lprn .and. nc_match.gt.0) write (iout,*) &
1625 "Shift:",-is,0," nc_match1",nc_match1,&
1626 " nc_match=",nc_match," req'd",nc_req
1627 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1628 nc_req.eq.0 .and. nc_match.eq.1) then
1633 ! Shift "left" strand forward
1634 call ncont_match(nc_match,nc_match1,is,0,ncont_ref,&
1635 icont_ref,ncont,icont,jfrag,llocal,lprn)
1636 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1637 if (lprn .and. nc_match.gt.0) write (iout,*) &
1638 "Shift:",is,0," nc_match1",nc_match1,&
1639 " nc_match=",nc_match," req'd",nc_req
1640 if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1641 nc_req.eq.0 .and. nc_match.eq.1) then
1647 if (nc_req.eq.0) return
1648 ! Shift "right" strand backward
1650 call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,&
1651 icont_ref,ncont,icont,jfrag,llocal,lprn)
1652 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1653 if (lprn .and. nc_match.gt.0) write (iout,*) &
1654 "Shift:",0,-is," nc_match1",nc_match1,&
1655 " nc_match=",nc_match," req'd",nc_req
1656 if (nc_match.ge.nc_req) then
1661 ! Shift "right" strand upward
1662 call ncont_match(nc_match,nc_match1,0,is,ncont_ref,&
1663 icont_ref,ncont,icont,jfrag,llocal,lprn)
1664 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1665 if (lprn .and. nc_match.gt.0) write (iout,*) &
1666 "Shift:",0,is," nc_match1",nc_match1,&
1667 " nc_match=",nc_match," req'd",nc_req
1668 if (nc_match.ge.nc_req) then
1674 ! Now try to shift both residues in contacts.
1678 call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,&
1679 icont_ref,ncont,icont,jfrag,llocal,lprn)
1680 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1681 if (lprn .and. nc_match.gt.0) write (iout,*) &
1682 "Shift:",-is,-js," nc_match1",nc_match1,&
1683 " nc_match=",nc_match," req'd",nc_req
1684 if (nc_match.ge.nc_req) then
1689 call ncont_match(nc_match,nc_match1,is,js,ncont_ref,&
1690 icont_ref,ncont,icont,jfrag,llocal,lprn)
1691 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1692 if (lprn .and. nc_match.gt.0) write (iout,*) &
1693 "Shift:",is,js," nc_match1",nc_match1,&
1694 " nc_match=",nc_match," req'd",nc_req
1695 if (nc_match.ge.nc_req) then
1701 call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,&
1702 icont_ref,ncont,icont,jfrag,llocal,lprn)
1703 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1704 if (lprn .and. nc_match.gt.0) write (iout,*) &
1705 "Shift:",-js,-is," nc_match1",nc_match1,&
1706 " nc_match=",nc_match," req'd",nc_req
1707 if (nc_match.ge.nc_req) then
1713 call ncont_match(nc_match,nc_match1,js,is,ncont_ref,&
1714 icont_ref,ncont,icont,jfrag,llocal,lprn)
1715 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1716 if (lprn .and. nc_match.gt.0) write (iout,*) &
1717 "Shift:",js,is," nc_match1",nc_match1,&
1718 " nc_match=",nc_match," req'd",nc_req
1719 if (nc_match.ge.nc_req) then
1726 if (is+js.le.n_shif1) then
1727 call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,&
1728 icont_ref,ncont,icont,jfrag,llocal,lprn)
1729 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1730 if (lprn .and. nc_match.gt.0) write (iout,*) &
1731 "Shift:",-is,js," nc_match1",nc_match1,&
1732 " nc_match=",nc_match," req'd",nc_req
1733 if (nc_match.ge.nc_req) then
1739 call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,&
1740 icont_ref,ncont,icont,jfrag,llocal,lprn)
1741 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1742 if (lprn .and. nc_match.gt.0) write (iout,*) &
1743 "Shift:",js,-is," nc_match1",nc_match1,&
1744 " nc_match=",nc_match," req'd",nc_req
1745 if (nc_match.ge.nc_req) then
1756 if (n_shif2.gt.0) then
1758 call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,&
1759 icont_ref,ncont,icont,jfrag,llocal,lprn)
1760 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1761 if (lprn .and. nc_match.gt.0) write (iout,*) &
1762 "Shift:",-is,-is," nc_match1",nc_match1,&
1763 " nc_match=",nc_match," req'd",nc_req
1764 if (nc_match.ge.nc_req) then
1769 call ncont_match(nc_match,nc_match1,is,is,ncont_ref,&
1770 icont_ref,ncont,icont,jfrag,llocal,lprn)
1771 if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1772 if (lprn .and. nc_match.gt.0) write (iout,*) &
1773 "Shift:",is,is," nc_match1",nc_match1,&
1774 " nc_match=",nc_match," req'd",nc_req
1775 if (nc_match.ge.nc_req) then
1782 ! If this point is reached, the contact maps are different.
1787 end subroutine match_contact
1788 !-------------------------------------------------------------------------
1789 subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,&
1790 ncont_ref,icont_ref,ncont,icont,jfrag,llocal,lprn)
1792 use energy_data, only:nnt,nct,maxcont
1793 ! implicit real*8 (a-h,o-z)
1794 ! include 'DIMENSIONS'
1795 ! include 'DIMENSIONS.ZSCOPT'
1796 ! include 'DIMENSIONS.COMPAR'
1797 ! include 'COMMON.IOUNITS'
1798 ! include 'COMMON.INTERACT'
1799 ! include 'COMMON.GEO'
1800 ! include 'COMMON.COMPAR'
1801 logical :: llocal,lprn
1802 integer ncont_ref,ncont,ishift,ishif2,nang_pair
1803 integer,dimension(2,maxcont) :: icont_ref,icont,icont_match !(2,maxcont)
1804 integer,dimension(2,nres) :: iang_pair !(2,maxres)
1805 integer :: nc_match,nc_match1,ishif1,jfrag
1806 integer :: i,j,ic1,ic2
1807 real(kind=8) :: diffang,fract,rad2deg
1809 ! Compare the contact map against the reference contact map; they're stored
1810 ! in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
1811 if (lprn) write (iout,'(80(1h*))')
1814 ! Check the local structure by comparing dihedral angles.
1815 ! write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal
1816 if (llocal .and. ncont_ref.eq.0) then
1817 ! If there are no contacts just compare the dihedral angles and exit.
1818 call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,&
1820 if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1821 " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract
1822 if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag)) &
1832 ic1=icont(1,i)+ishif1
1833 ic2=icont(2,i)+ishif2
1834 ! write (iout,*) "i",i," ic1",ic1," ic2",ic2
1835 if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
1837 if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
1838 nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
1839 nc_match1=nc_match1+1
1840 icont_match(1,nc_match1)=ic1
1841 icont_match(2,nc_match1)=ic2
1842 ! call add_angpair(icont(1,i),icont_ref(1,j),
1843 ! & nang_pair,iang_pair)
1844 ! call add_angpair(icont(2,i),icont_ref(2,j),
1845 ! & nang_pair,iang_pair)
1846 if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),&
1847 " match",icont_ref(1,j),icont_ref(2,j),&
1848 " shifts",ishif1,ishif2
1855 write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
1856 write (iout,*) "icont_match"
1858 write (iout,*) icont_match(1,i),icont_match(2,i)
1861 if (llocal .and. nc_match.gt.0) then
1862 call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,&
1863 ang_cut1(jfrag),diffang,fract)
1864 if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1865 " ang_cut:",ang_cut(jfrag)*rad2deg,&
1866 " ang_cut1",ang_cut1(jfrag)*rad2deg
1867 if (diffang.gt.ang_cut(jfrag) &
1868 .or. fract.lt.frac_min(jfrag)) nc_match=0
1870 ! if (nc_match.gt.0) then
1871 ! diffang = angnorm1(nang_pair,iang_pair,lprn)
1872 ! if (diffang.gt.ang_cut(jfrag)) nc_match=0
1874 if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,&
1875 " diffang",rad2deg*diffang," nc_match",nc_match
1877 end subroutine ncont_match
1878 !------------------------------------------------------------------------------
1879 subroutine match_secondary(jfrag,isecstr,nsec_match,lprn)
1880 ! This subroutine compares the secondary structure (isecstr) of fragment jfrag
1881 ! conformation considered to that of the reference conformation.
1882 ! Returns the number of equivalent residues (nsec_match).
1883 ! implicit real*8 (a-h,o-z)
1884 ! include 'DIMENSIONS'
1885 ! include 'DIMENSIONS.ZSCOPT'
1886 ! include 'DIMENSIONS.COMPAR'
1887 ! include 'COMMON.IOUNITS'
1888 ! include 'COMMON.CHAIN'
1889 ! include 'COMMON.PEPTCONT'
1890 ! include 'COMMON.COMPAR'
1892 integer :: isecstr(nres)
1893 integer :: jfrag,nsec_match,npart,i,j
1894 npart = npiece(jfrag,1)
1897 write (iout,*) "match_secondary jfrag",jfrag," ifrag",&
1898 (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart)
1899 write (iout,'(80i1)') (isec_ref(j),j=1,nres)
1900 write (iout,'(80i1)') (isecstr(j),j=1,nres)
1903 do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
1904 ! The residue has equivalent conformational state to that of the reference
1906 ! a) the conformational states are equal or
1907 ! b) the reference state is a coil and that of the conformation considered
1909 ! c) the conformational state of the conformation considered is a strand
1910 ! and that of the reference conformation is a coil.
1911 ! 10/28/02 - case (b) deleted.
1912 if (isecstr(j).eq.isec_ref(j) .or. &
1913 ! & isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or.
1914 isec_ref(j).eq.0 .and. isecstr(j).eq.1) &
1915 nsec_match=nsec_match+1
1919 end subroutine match_secondary
1920 !------------------------------------------------------------------------------
1922 !------------------------------------------------------------------------------
1923 subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd)
1925 use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
1927 ! implicit real*8 (a-h,o-z)
1928 real(kind=8),dimension(3) :: r1,r2,a,b,x,y
1929 real(kind=8) :: uu,vv,aa,bb,dd
1930 real(kind=8) :: ab,ar,br,det,dd1,dd2,dd3,dd4,dd5
1931 !el odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
1932 !el + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
1933 ! print *,"r1",(r1(i),i=1,3)
1934 ! print *,"r2",(r2(i),i=1,3)
1935 ! print *,"a",(a(i),i=1,3)
1936 ! print *,"b",(b(i),i=1,3)
1937 aa = a(1)**2+a(2)**2+a(3)**2
1938 bb = b(1)**2+b(2)**2+b(3)**2
1939 ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
1940 ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3))
1941 br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3))
1943 ! print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det
1944 uu = (-ar*bb+br*ab)/det
1945 vv = (br*aa-ar*ab)/det
1951 !el dd1 = odl(uu,vv)
1952 dd1 = odl(uu,vv,r1,r2,ar,br,ab,aa,bb)
1953 !el dd2 = odl(0.0d0,0.0d0)
1954 dd2 = odl(0.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1955 !el dd3 = odl(0.0d0,1.0d0)
1956 dd3 = odl(0.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1957 !el dd4 = odl(1.0d0,0.0d0)
1958 dd4 = odl(1.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1959 !el dd5 = odl(1.0d0,1.0d0)
1960 dd5 = odl(1.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1961 dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5))
1965 else if (dd.eq.dd3) then
1968 else if (dd.eq.dd4) then
1971 else if (dd.eq.dd5) then
1980 ! dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2
1984 ! write (8,*) uu,vv,dd,dd1
1987 end subroutine odlodc
1988 !------------------------------------------------------------------------------
1989 real(kind=8) function odl(u,v,r1,r2,ar,br,ab,aa,bb)
1991 real(kind=8),dimension(3) :: r1,r2
1992 real(kind=8) :: aa,bb,u,v,ar,br,ab
1994 odl = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
1995 + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
1998 !------------------------------------------------------------------------------
2000 !------------------------------------------------------------------------------
2001 subroutine proc_cont
2003 use geometry_data, only:rad2deg
2004 use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
2006 ! implicit real*8 (a-h,o-z)
2007 ! include 'DIMENSIONS'
2008 ! include 'DIMENSIONS.ZSCOPT'
2009 ! include 'DIMENSIONS.COMPAR'
2010 ! include 'COMMON.IOUNITS'
2011 ! include 'COMMON.TIME1'
2012 ! include 'COMMON.SBRIDGE'
2013 ! include 'COMMON.CONTROL'
2014 ! include 'COMMON.COMPAR'
2015 ! include 'COMMON.CHAIN'
2016 ! include 'COMMON.HEADER'
2017 ! include 'COMMON.CONTACTS1'
2018 ! include 'COMMON.PEPTCONT'
2019 ! include 'COMMON.GEO'
2020 integer :: i,j,k,ind,len_cut,ndigit,length_frag
2022 write (iout,*) "proc_cont: nlevel",nlevel
2023 if (nlevel.lt.0) then
2024 write (iout,*) "call define_fragments"
2025 call define_fragments
2027 write (iout,*) "call secondary2"
2028 call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2031 write (iout,'(80(1h=))')
2032 write (iout,*) "Electrostatic contacts"
2033 call contacts_between_fragments(.true.,0,ncont_pept_ref,&
2034 icont_pept_ref,ncont_frag_ref(1),icont_frag_ref(1,1,1))
2035 write (iout,'(80(1h=))')
2036 write (iout,*) "Side chain contacts"
2037 call contacts_between_fragments(.true.,0,ncont_ref,&
2038 icont_ref,nsccont_frag_ref(1),isccont_frag_ref(1,1,1))
2039 if (nlevel.lt.0) then
2043 if (istruct(i).le.1) then
2044 len_cut=max0(len_frag(i,1)*4/5,3)
2045 else if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2046 len_cut=max0(len_frag(i,1)*2/5,3)
2048 write (iout,*) "i",i," istruct",istruct(i)," ncont_frag",&
2049 ncont_frag_ref(ind)," len_cut",len_cut,&
2050 " icont_single",icont_single," iloc_single",iloc_single
2052 if (iloc(i).gt.0) write (iout,*) &
2053 "Local structure used to compare structure of fragment",i,&
2055 if (istruct(i).ne.3 .and. istruct(i).ne.0 &
2056 .and. icont_single.gt.0 .and. &
2057 ncont_frag_ref(ind).ge.len_cut) then
2058 write (iout,*) "Electrostatic contacts used to compare",&
2059 " structure of fragment",i," to native."
2062 else if (icont_single.gt.0 .and. nsccont_frag_ref(ind) &
2064 write (iout,*) "Side chain contacts used to compare",&
2065 " structure of fragment",i," to native."
2069 write (iout,*) "Contacts not used to compare",&
2070 " structure of fragment",i," to native."
2075 if (irms_single.gt.0 .or. isccont(i,1).eq.0 &
2076 .and. ielecont(i,1).eq.0) then
2077 write (iout,*) "RMSD used to compare",&
2078 " structure of fragment",i," to native."
2081 write (iout,*) "RMSD not used to compare",&
2082 " structure of fragment",i," to native."
2087 if (nlevel.lt.-1) then
2090 if (nlevel.gt.3) nlevel=3
2091 if (nlevel.eq.3) then
2093 npiece(1,3)=nfrag(1)
2103 else if (nlevel.eq.-1) then
2108 isnfrag(i+1)=isnfrag(i)+nfrag(i)
2112 ndigit=ndigit+2*nfrag(i)
2114 write (iout,*) "ndigit",ndigit
2115 if (.not.binary .and. ndigit.gt.30) then
2116 write (iout,*) "Highest class too large; switching to",&
2117 " binary representation."
2120 write (iout,*) "isnfrag",(isnfrag(i),i=1,nlevel+1)
2121 write(iout,*) "rmscut_base_up",rmscut_base_up,&
2122 " rmscut_base_low",rmscut_base_low," rmsup_lim",rmsup_lim
2128 length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2132 length_frag=length_frag+len_frag(ipiece(k,j,i),1)
2135 len_frag(j,i)=length_frag
2136 rmscutfrag(1,j,i)=rmscut_base_up*length_frag
2137 rmscutfrag(2,j,i)=rmscut_base_low*length_frag
2138 if (rmscutfrag(1,j,i).lt.rmsup_lim) &
2139 rmscutfrag(1,j,i)=rmsup_lim
2140 if (rmscutfrag(1,j,i).gt.rmsupup_lim) &
2141 rmscutfrag(1,j,i)=rmsupup_lim
2144 write (iout,*) "Level",1," number of fragments:",nfrag(1)
2146 write (iout,*) npiece(j,1),(ifrag(1,k,j),ifrag(2,k,j),&
2147 k=1,npiece(j,1)),len_frag(j,1),rmscutfrag(1,j,1),&
2148 rmscutfrag(2,j,1),n_shift(1,j,1),n_shift(2,j,1),&
2149 ang_cut(j)*rad2deg,ang_cut1(j)*rad2deg,frac_min(j),&
2150 nc_fragm(j,1),nc_req_setf(j,1),istruct(j)
2153 write (iout,*) "Level",i," number of fragments:",nfrag(i)
2155 write (iout,*) npiece(j,i),(ipiece(k,j,i),&
2156 k=1,npiece(j,i)),len_frag(j,i),rmscutfrag(1,j,i),&
2157 rmscutfrag(2,j,i),n_shift(1,j,i),n_shift(2,j,i),&
2158 nc_fragm(j,i),nc_req_setf(j,i)
2162 end subroutine proc_cont
2163 !------------------------------------------------------------------------------
2165 !------------------------------------------------------------------------------
2166 subroutine define_pairs
2168 ! use energy_data, only:nsccont_frag_ref
2169 ! implicit real*8 (a-h,o-z)
2170 ! include 'DIMENSIONS'
2171 ! include 'DIMENSIONS.ZSCOPT'
2172 ! include 'DIMENSIONS.COMPAR'
2173 ! include 'COMMON.IOUNITS'
2174 ! include 'COMMON.TIME1'
2175 ! include 'COMMON.SBRIDGE'
2176 ! include 'COMMON.CONTROL'
2177 ! include 'COMMON.COMPAR'
2178 ! include 'COMMON.FRAG'
2179 ! include 'COMMON.CHAIN'
2180 ! include 'COMMON.HEADER'
2181 ! include 'COMMON.GEO'
2182 ! include 'COMMON.CONTACTS1'
2183 ! include 'COMMON.PEPTCONT'
2184 integer :: j,k,i,length_frag,ind,ll1,ll2,len_cut
2189 length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2191 len_frag(j,1)=length_frag
2192 write (iout,*) "Fragment",j," length",len_frag(j,1)
2198 if (istruct(i).le.1 .or. istruct(j).le.1) then
2199 if (istruct(i).le.1) then
2204 if (istruct(j).le.1) then
2209 len_cut=max0(min0(ll1*2/3,ll2*4/5),3)
2211 if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2216 if (istruct(j).eq.2 .or. istruct(j).eq.4) then
2221 len_cut=max0(min0(ll1*4/5,ll2)*4/5,3)
2223 write (iout,*) "Fragments",i,j," structure",istruct(i),&
2224 istruct(j)," # contacts",&
2225 ncont_frag_ref(ind),nsccont_frag_ref(ind),&
2226 " lengths",len_frag(i,1),len_frag(j,1),&
2227 " ll1",ll1," ll2",ll2," len_cut",len_cut
2228 if ((istruct(i).eq.1 .or. istruct(j).eq.1) .and. &
2229 nsccont_frag_ref(ind).ge.len_cut ) then
2230 if (istruct(i).eq.1 .and. istruct(j).eq.1) then
2231 write (iout,*) "Adding pair of helices",i,j,&
2232 " based on SC contacts"
2234 write (iout,*) "Adding helix+strand/sheet pair",i,j,&
2235 " based on SC contacts"
2238 if (icont_pair.gt.0) then
2239 write (iout,*) "# SC contacts will be used",&
2241 isccont(nfrag(2),2)=1
2243 if (irms_pair.gt.0) then
2244 write (iout,*) "Fragment RMSD will be used",&
2248 npiece(nfrag(2),2)=2
2249 ipiece(1,nfrag(2),2)=i
2250 ipiece(2,nfrag(2),2)=j
2251 ielecont(nfrag(2),2)=0
2252 n_shift(1,nfrag(2),2)=nshift_pair
2253 n_shift(2,nfrag(2),2)=nshift_pair
2254 nc_fragm(nfrag(2),2)=ncfrac_pair
2255 nc_req_setf(nfrag(2),2)=ncreq_pair
2256 else if ((istruct(i).ge.2 .and. istruct(i).le.4) &
2257 .and. (istruct(j).ge.2 .and. istruct(i).le.4) &
2258 .and. ncont_frag_ref(ind).ge.len_cut ) then
2260 write (iout,*) "Adding pair strands/sheets",i,j,&
2261 " based on pp contacts"
2262 if (icont_pair.gt.0) then
2263 write (iout,*) "# pp contacts will be used",&
2265 ielecont(nfrag(2),2)=1
2267 if (irms_pair.gt.0) then
2268 write (iout,*) "Fragment RMSD will be used",&
2272 npiece(nfrag(2),2)=2
2273 ipiece(1,nfrag(2),2)=i
2274 ipiece(2,nfrag(2),2)=j
2275 ielecont(nfrag(2),2)=1
2276 isccont(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_bet
2280 nc_req_setf(nfrag(2),2)=ncreq_bet
2284 write (iout,*) "Pairs found"
2286 write (iout,*) ipiece(1,i,2),ipiece(2,i,2)
2289 end subroutine define_pairs
2290 !------------------------------------------------------------------------------
2292 !------------------------------------------------------------------------------
2293 INTEGER FUNCTION ICANT(I,J)
2302 !------------------------------------------------------------------------------
2304 !------------------------------------------------------------------------------
2305 subroutine imysort(n, m, mm, x, y, z, z1, z2, z3, z4, z5, z6)
2308 integer :: x(m,mm,n),y(n),z(n),z1(2,n),z6(n),xmin,xtemp
2309 real(kind=8) :: z2(n),z3(n),z4(n),z5(n)
2310 real(kind=8) :: xxtemp
2311 integer :: i,j,k,imax
2316 if (x(1,1,j).lt.xmin) then
2350 x(j,k,imax)=x(j,k,i)
2356 end subroutine imysort
2357 !------------------------------------------------------------------------------
2359 !-------------------------------------------------------------------------------
2360 real(kind=8) function qwolynes(ilevel,jfrag)
2362 use geometry_data, only:cref,nperm
2363 use control_data, only:symetr
2364 use energy_data, only:nnt,nct,itype
2366 ! include 'DIMENSIONS'
2367 ! include 'DIMENSIONS.ZSCOPT'
2368 ! include 'DIMENSIONS.COMPAR'
2369 ! include 'COMMON.IOUNITS'
2370 ! include 'COMMON.COMPAR'
2371 ! include 'COMMON.CHAIN'
2372 ! include 'COMMON.INTERACT'
2373 ! include 'COMMON.CONTROL'
2374 integer :: ilevel,jfrag,kkk
2375 integer :: i,j,jl,k,l,il,kl,nl,np,ip,kp
2377 real(kind=8),dimension(:),allocatable :: tempus !(maxperm)
2378 real(kind=8) :: maxiQ !dist,
2379 real(kind=8) :: qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
2380 logical :: lprn=.false.
2381 real(kind=8) :: x !el sigm
2382 !el sigm(x)=0.25d0*x
2388 ! write (iout,*) "QWolyes: " jfrag",jfrag,
2389 ! & " ilevel",ilevel
2390 allocate(tempus(nperm))
2393 if (ilevel.eq.0) then
2394 if (lprn) write (iout,*) "Q computed for whole molecule"
2405 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2406 (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2407 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2409 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2410 if (itype(il).ne.10 .or. itype(jl).ne.10) then
2413 (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2414 (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2415 (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2416 dijCM=dist(il+nres,jl+nres)
2417 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2421 write (iout,*) "il",il," jl",jl,&
2422 " itype",itype(il),itype(jl)
2423 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2424 " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2429 if (lprn) write (iout,*) "nl",nl," qq",qq
2430 else if (ilevel.eq.1) then
2431 if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag
2433 ! write (iout,*) "nlist_frag",nlist_frag(jfrag)
2434 do i=2,nlist_frag(jfrag)
2436 il=list_frag(i,jfrag)
2437 jl=list_frag(j,jfrag)
2438 if (iabs(il-jl).gt.nsep) then
2446 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2447 (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2448 (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2450 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2451 if (itype(il).ne.10 .or. itype(jl).ne.10) then
2454 (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2455 (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2456 (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2457 dijCM=dist(il+nres,jl+nres)
2458 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2462 write (iout,*) "i",i," j",j," il",il," jl",jl,&
2463 " itype",itype(il),itype(jl)
2464 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2465 " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2471 if (lprn) write (iout,*) "nl",nl," qq",qq
2472 else if (ilevel.eq.2) then
2473 np=npiece(jfrag,ilevel)
2476 ip=ipiece(i,jfrag,ilevel)
2477 do j=1,nlist_frag(ip)
2480 kp=ipiece(k,jfrag,ilevel)
2481 do l=1,nlist_frag(kp)
2483 if (iabs(kl-il).gt.nsep) then
2491 d0ij=dsqrt((cref(1,kl,kkk)-cref(1,il,kkk))**2+ &
2492 (cref(2,kl,kkk)-cref(2,il,kkk))**2+ &
2493 (cref(3,kl,kkk)-cref(3,il,kkk))**2)
2495 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2496 if (itype(il).ne.10 .or. itype(kl).ne.10) then
2499 (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2500 (cref(2,kl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2501 (cref(3,kl+nres,kkk)-cref(3,il+nres,kkk))**2)
2502 dijCM=dist(il+nres,kl+nres)
2503 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/ &
2508 write (iout,*) "i",i," j",j," k",k," l",l," il",il,&
2509 " kl",kl," itype",itype(il),itype(kl)
2510 write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",&
2511 d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2519 if (lprn) write (iout,*) "nl",nl," qq",qq
2521 write (iout,*)"Error: Q can be computed only for level 1 and 2."
2526 if (maxiQ.le.tempus(kkk)) maxiQ=tempus(kkk)
2528 qwolynes=1.0d0-maxiQ
2531 end function qwolynes
2532 !-------------------------------------------------------------------------------
2533 real(kind=8) function sigm(x)
2538 !-------------------------------------------------------------------------------
2539 subroutine fragment_list
2541 ! include 'DIMENSIONS'
2542 ! include 'DIMENSIONS.ZSCOPT'
2543 ! include 'DIMENSIONS.COMPAR'
2544 ! include 'COMMON.IOUNITS'
2545 ! include 'COMMON.COMPAR'
2546 logical :: lprn=.true.
2547 integer :: i,ilevel,j,k,jfrag
2550 do i=1,npiece(jfrag,1)
2551 if (lprn) write (iout,*) "jfrag=",jfrag,&
2552 "i=",i," fragment",ifrag(1,i,jfrag),&
2554 do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
2555 do k=1,nlist_frag(jfrag)
2556 if (list_frag(k,jfrag).eq.j) goto 10
2558 nlist_frag(jfrag)=nlist_frag(jfrag)+1
2559 list_frag(nlist_frag(jfrag),jfrag)=j
2564 write (iout,*) "Fragment list"
2566 write (iout,*)"Fragment",j," list",(list_frag(k,j),&
2570 end subroutine fragment_list
2571 !-------------------------------------------------------------------------------
2572 real(kind=8) function rmscalc(ishif,i,j,jcon,lprn)
2575 use control_data, only:symetr
2576 use geometry_data, only:nperm
2577 ! implicit real*8 (a-h,o-z)
2578 ! include 'DIMENSIONS'
2579 ! include 'DIMENSIONS.ZSCOPT'
2580 ! include 'DIMENSIONS.COMPAR'
2581 ! include 'COMMON.IOUNITS'
2582 ! include 'COMMON.COMPAR'
2583 ! include 'COMMON.CHAIN'
2584 ! include 'COMMON.INTERACT'
2585 ! include 'COMMON.VAR'
2586 ! include 'COMMON.CONTROL'
2587 real(kind=8) :: przes(3),obrot(3,3)
2588 !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2589 !el logical :: iadded(nres)
2590 !el integer :: inumber(2,nres)
2591 !el common /ccc/ creff,cc,iadded,inumber
2594 integer :: ishif,i,j,jcon,idup,kkk,l,k,kk
2595 real(kind=8) :: rminrms,rms
2597 write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
2598 write (iout,*) "npiece",npiece(j,i)
2601 ! write (iout,*) "symetr",symetr
2607 ! write (iout,*) "nperm",nperm
2614 ! write (iout,*) "kkk",kkk
2619 write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",&
2620 ifrag(1,k,j),ifrag(2,k,j)
2623 call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
2624 ! write (iout,*) "Exit cprep"
2626 ! write (iout,*) "ii=",ii
2629 ! write (iout,*) "kk",kk," npiece",npiece(kk,1)
2632 write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,&
2633 " l=",l," adding fragment",&
2634 ifrag(1,l,kk),ifrag(2,l,kk)
2637 call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
2638 ! write (iout,*) "After cprep"
2645 write (iout,*) "tuszukaj"
2648 write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),&
2649 inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
2656 call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
2658 print *,'Error: FITSQ non-convergent, jcon',jcon,i
2660 else if (rms.lt.-1.0d-6) then
2661 print *,'Error: rms^2 = ',rms,jcon,i
2663 else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2666 ! write (iout,*) "rmsmin", rminrms, "rms", rms
2667 if (rms.le.rminrms) rminrms=rms
2669 rmscalc = dsqrt(rminrms)
2670 ! write (iout, *) "analysys", rmscalc,anatemp
2672 end function rmscalc
2673 !-------------------------------------------------------------------------
2674 subroutine cprep(if1,if2,ishif,idup,kwa)
2677 use control_data, only:symetr
2678 use geometry_data, only:nperm,cref,c
2679 ! implicit real*8 (a-h,o-z)
2680 ! include 'DIMENSIONS'
2681 ! include 'DIMENSIONS.ZSCOPT'
2682 ! include 'DIMENSIONS.COMPAR'
2683 ! include 'COMMON.CONTROL'
2684 ! include 'COMMON.IOUNITS'
2685 ! include 'COMMON.COMPAR'
2686 ! include 'COMMON.CHAIN'
2687 ! include 'COMMON.INTERACT'
2688 ! include 'COMMON.VAR'
2689 real(kind=8) :: przes(3),obrot(3,3)
2690 !el real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2691 !el logical :: iadded(nres)
2692 !el integer :: inumber(2,nres)
2693 integer :: iistrart,kwa,blar
2694 !el common /ccc/ creff,cc,iadded,inumber
2695 integer :: if1,if2,ishif,idup,kkk,l,m
2696 ! write (iout,*) "Calling cprep symetr",symetr," kwa",kwa
2701 ! write (iout,*) "nperm",nperm
2705 ! write (iout,*) "l",l," iadded",iadded(l)
2707 if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) &
2712 inumber(2,idup)=l+ishif
2714 creff(m,idup)=cref(m,l,kkk)
2715 cc(m,idup)=c(m,l+ishif)
2720 end subroutine cprep
2721 !-------------------------------------------------------------------------
2722 real(kind=8) function rmsnat(jcon)
2724 use control_data, only:symetr
2725 use geometry_data, only:nperm,cref,c
2726 use energy_data, only:itype
2727 ! implicit real*8 (a-h,o-z)
2728 ! include 'DIMENSIONS'
2729 ! include 'DIMENSIONS.ZSCOPT'
2730 ! include 'DIMENSIONS.COMPAR'
2731 ! include 'COMMON.IOUNITS'
2732 ! include 'COMMON.COMPAR'
2733 ! include 'COMMON.CHAIN'
2734 ! include 'COMMON.INTERACT'
2735 ! include 'COMMON.VAR'
2736 ! include 'COMMON.CONTROL'
2737 real(kind=8) :: przes(3),obrot(3,3),cc(3,2*nres),ccref(3,2*nres)
2739 integer :: ishif,i,j,resprzesun,jcon,kkk,nnsup
2740 real(kind=8) :: rminrms,rmsminsing,rms
2750 if (itype(i).ne.ntyp1) then
2754 ccref(j,nnsup)=cref(j,i,kkk)
2758 call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv)
2760 print *,'Error: FITSQ non-convergent, jcon',jcon,i
2762 else if (rms.lt.-1.0d-6) then
2763 print *,'Error: rms^2 = ',rms,jcon,i
2765 else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2768 if (rms.le.rminrms) rminrms=rms
2769 ! write (iout,*) "kkk",kkk," rmsnat",rms , rminrms
2771 rmsnat = dsqrt(rminrms)
2772 ! write (iout,*) "analysys",rmsnat, anatemp
2773 ! liczenie rmsdla pojedynczego lancucha
2776 !-------------------------------------------------------------------------------
2777 subroutine define_fragments
2779 use geometry_data, only:rad2deg
2780 use energy_data, only:itype
2781 use compare_data, only:nhfrag,nbfrag,bfrag,hfrag
2782 ! implicit real*8 (a-h,o-z)
2783 ! include 'DIMENSIONS'
2784 ! include 'DIMENSIONS.ZSCOPT'
2785 ! include 'DIMENSIONS.COMPAR'
2786 ! include 'COMMON.IOUNITS'
2787 ! include 'COMMON.TIME1'
2788 ! include 'COMMON.FRAG'
2789 ! include 'COMMON.SBRIDGE'
2790 ! include 'COMMON.CONTROL'
2791 ! include 'COMMON.COMPAR'
2792 ! include 'COMMON.CHAIN'
2793 ! include 'COMMON.HEADER'
2794 ! include 'COMMON.GEO'
2795 ! include 'COMMON.CONTACTS'
2796 ! include 'COMMON.PEPTCONT'
2797 ! include 'COMMON.INTERACT'
2798 ! include 'COMMON.NAMES'
2799 integer :: nstrand,istrand(2,nres/2)
2800 integer :: nhairp,ihairp(2,nres/5)
2801 character(len=16) :: strstr(4)=reshape((/'helix','hairpin',&
2802 'strand','strand pair'/),shape(strstr))
2803 integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4
2805 write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2806 'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2807 'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2808 ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet
2809 write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2810 ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2811 write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2812 ' MAXANG_HEL',angcut1_hel*rad2deg
2813 write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2814 ' MAXANG_BET',angcut1_bet*rad2deg
2815 write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2816 ' MAXANG_STRAND',angcut1_strand*rad2deg
2817 write (iout,*) 'FRAC_MIN',frac_min_set
2818 ! Find secondary structure elements (helices and beta-sheets)
2819 call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2821 ! Define primary fragments. First include the helices.
2825 ! AL 12/23/03 - to avoid splitting helices into very small fragments
2826 if (merge_helices) then
2827 write (iout,*) "Before merging helices: nhfrag",nhfrag
2829 write (2,*) hfrag(1,i),hfrag(2,i)
2832 do while (i.lt.nhfrag)
2833 if (hfrag(1,i+1)-hfrag(2,i).le.1) then
2835 hfrag(2,i)=hfrag(2,i+1)
2837 hfrag(1,j)=hfrag(1,j+1)
2838 hfrag(2,j)=hfrag(2,j+1)
2843 write (iout,*) "After merging helices: nhfrag",nhfrag
2845 write (2,*) hfrag(1,i),hfrag(2,i)
2851 ifrag(1,1,i)=hfrag(1,i)
2852 ifrag(2,1,i)=hfrag(2,i)
2854 n_shift(2,i,1)=nshift_hel
2855 ang_cut(i)=angcut_hel
2856 ang_cut1(i)=angcut1_hel
2857 frac_min(i)=frac_min_set
2858 nc_fragm(i,1)=ncfrac_hel
2859 nc_req_setf(i,1)=ncreq_hel
2862 write (iout,*) "isplit_bet",isplit_bet
2863 if (isplit_bet.gt.1) then
2864 ! Split beta-sheets into strands and store strands as primary fragments.
2865 call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2869 ifrag(1,1,ii)=istrand(1,i)
2870 ifrag(2,1,ii)=istrand(2,i)
2871 n_shift(1,ii,1)=nshift_strand
2872 n_shift(2,ii,1)=nshift_strand
2873 ang_cut(ii)=angcut_strand
2874 ang_cut1(ii)=angcut1_strand
2875 frac_min(ii)=frac_min_set
2880 nfrag(1)=nfrag(1)+nstrand
2881 else if (isplit_bet.eq.1) then
2882 ! Split only far beta-sheets; does not split hairpins.
2883 call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2884 call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2888 ifrag(1,1,ii)=ihairp(1,i)
2889 ifrag(2,1,ii)=ihairp(2,i)
2890 n_shift(1,ii,1)=nshift_bet
2891 n_shift(2,ii,1)=nshift_bet
2892 ang_cut(ii)=angcut_bet
2893 ang_cut1(ii)=angcut1_bet
2894 frac_min(ii)=frac_min_set
2895 nc_fragm(ii,1)=ncfrac_bet
2896 nc_req_setf(ii,1)=ncreq_bet
2899 nfrag(1)=nfrag(1)+nhairp
2903 ifrag(1,1,ii)=istrand(1,i)
2904 ifrag(2,1,ii)=istrand(2,i)
2905 n_shift(1,ii,1)=nshift_strand
2906 n_shift(2,ii,1)=nshift_strand
2907 ang_cut(ii)=angcut_strand
2908 ang_cut1(ii)=angcut1_strand
2909 frac_min(ii)=frac_min_set
2914 nfrag(1)=nfrag(1)+nstrand
2916 ! Do not split beta-sheets; each pair of strands is a primary element.
2917 call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2921 ifrag(1,1,ii)=ihairp(1,i)
2922 ifrag(2,1,ii)=ihairp(2,i)
2923 n_shift(1,ii,1)=nshift_bet
2924 n_shift(2,ii,1)=nshift_bet
2925 ang_cut(ii)=angcut_bet
2926 ang_cut1(ii)=angcut1_bet
2927 frac_min(ii)=frac_min_set
2928 nc_fragm(ii,1)=ncfrac_bet
2929 nc_req_setf(ii,1)=ncreq_bet
2932 nfrag(1)=nfrag(1)+nhairp
2936 ifrag(1,1,ii)=bfrag(1,i)
2937 ifrag(2,1,ii)=bfrag(2,i)
2938 if (bfrag(3,i).lt.bfrag(4,i)) then
2939 ifrag(1,2,ii)=bfrag(3,i)
2940 ifrag(2,2,ii)=bfrag(4,i)
2942 ifrag(1,2,ii)=bfrag(4,i)
2943 ifrag(2,2,ii)=bfrag(3,i)
2945 n_shift(1,ii,1)=nshift_bet
2946 n_shift(2,ii,1)=nshift_bet
2947 ang_cut(ii)=angcut_bet
2948 ang_cut1(ii)=angcut1_bet
2949 frac_min(ii)=frac_min_set
2950 nc_fragm(ii,1)=ncfrac_bet
2951 nc_req_setf(ii,1)=ncreq_bet
2954 nfrag(1)=nfrag(1)+nbfrag
2956 write (iout,*) "The following primary fragments were found:"
2957 write (iout,*) "Helices:",nhfrag
2963 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
2964 i,restyp(it1),i1,restyp(it2),i2
2966 write (iout,*) "Hairpins:",nhairp
2967 do i=nhfrag+1,nhfrag+nhairp
2972 write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') &
2973 i,restyp(it1),i1,restyp(it2),i2
2975 write (iout,*) "Far strand pairs:",nbfrag
2976 do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
2985 write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') &
2986 i,restyp(it1),i1,restyp(it2),i2,&
2987 restyp(it3),i3,restyp(it4),i4
2989 write (iout,*) "Strands:",nstrand
2990 do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
2995 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
2996 i,restyp(it1),i1,restyp(it2),i2
2998 call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),&
2999 istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),&
3000 nc_fragm(1,1),nc_req_setf(1,1))
3001 write (iout,*) "Fragments after sorting:"
3007 write (iout,'(i3,2x,a,i4,2x,a,i4,$)') &
3008 i,restyp(it1),i1,restyp(it2),i2
3009 if (npiece(i,1).eq.1) then
3010 write (iout,'(2x,a)') strstr(istruct(i))
3016 write (iout,'(2x,a,i4,2x,a,i4,2x,a)') &
3017 restyp(it1),i1,restyp(it2),i2,strstr(istruct(i))
3021 end subroutine define_fragments
3022 !------------------------------------------------------------------------------
3023 subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
3024 ! implicit real*8 (a-h,o-z)
3025 ! include 'DIMENSIONS'
3026 ! include 'DIMENSIONS.ZSCOPT'
3027 ! include 'DIMENSIONS.COMPAR'
3028 ! include 'COMMON.IOUNITS'
3029 integer :: nbfrag,bfrag(4,nres/3)
3030 integer :: nhairp,ihairp(2,nres/5)
3032 write (iout,*) "Entered find_and_remove_hairpins"
3033 write (iout,*) "nbfrag",nbfrag
3035 write (iout,*) i,(bfrag(k,i),k=1,4)
3039 do while (i.le.nbfrag)
3040 write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
3041 if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) &
3043 write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
3045 ihairp(1,nhairp)=bfrag(1,i)
3046 ihairp(2,nhairp)=bfrag(3,i)
3050 bfrag(k,j)=bfrag(k,j+1)
3057 write (iout,*) "After finding hairpins:"
3058 write (iout,*) "nhairp",nhairp
3060 write (iout,*) i,ihairp(1,i),ihairp(2,i)
3062 write (iout,*) "nbfrag",nbfrag
3064 write (iout,*) i,(bfrag(k,i),k=1,4)
3067 end subroutine find_and_remove_hairpins
3068 !------------------------------------------------------------------------------
3069 subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
3070 ! implicit real*8 (a-h,o-z)
3071 ! include 'DIMENSIONS'
3072 ! include 'DIMENSIONS.ZSCOPT'
3073 ! include 'DIMENSIONS.COMPAR'
3074 ! include 'COMMON.IOUNITS'
3075 integer :: nbfrag,bfrag(4,nres/3)
3076 integer :: nstrand,istrand(2,nres/2)
3077 integer :: nhairp,ihairp(2,nres/5)
3080 write (iout,*) "Entered split_beta"
3081 write (iout,*) "nbfrag",nbfrag
3083 write (iout,*) i,(bfrag(k,i),k=1,4)
3087 write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
3088 call add_strand(nstrand,istrand,nhairp,ihairp,&
3089 bfrag(1,i),bfrag(2,i),found)
3090 if (bfrag(3,i).lt.bfrag(4,i)) then
3091 write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
3092 call add_strand(nstrand,istrand,nhairp,ihairp,&
3093 bfrag(3,i),bfrag(4,i),found)
3095 write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
3096 call add_strand(nstrand,istrand,nhairp,ihairp,&
3097 bfrag(4,i),bfrag(3,i),found)
3101 write (iout,*) "Strands found:",nstrand
3103 write (iout,*) i,istrand(1,i),istrand(2,i)
3106 end subroutine split_beta
3107 !------------------------------------------------------------------------------
3108 subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
3109 ! implicit real*8 (a-h,o-z)
3110 ! include 'DIMENSIONS'
3111 ! include 'DIMENSIONS.ZSCOPT'
3112 ! include 'DIMENSIONS.COMPAR'
3113 ! include 'COMMON.IOUNITS'
3114 integer :: nstrand,istrand(2,nres/2)
3115 integer :: nhairp,ihairp(2,nres/5)
3117 integer :: is1,is2,j,idelt
3120 idelt=(ihairp(2,j)-ihairp(1,j))/6
3121 if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
3122 write (iout,*) "strand",is1,is2," is part of hairpin",&
3123 ihairp(1,j),ihairp(2,j)
3128 idelt=(istrand(2,j)-istrand(1,j))/3
3129 if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) &
3131 ! The strand already exists in the array; update its ends if necessary.
3132 write (iout,*) "strand",is1,is2," found at position",j,&
3133 ":",istrand(1,j),istrand(2,j)
3134 istrand(1,j)=min0(istrand(1,j),is1)
3135 istrand(2,j)=max0(istrand(2,j),is2)
3139 ! The strand has not been found; add it to the array.
3140 write (iout,*) "strand",is1,is2," added to the array."
3143 istrand(1,nstrand)=is1
3144 istrand(2,nstrand)=is2
3146 end subroutine add_strand
3147 !------------------------------------------------------------------------------
3148 subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
3150 use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup
3151 use energy_data, only:itype,maxcont
3152 use compare_data, only:bfrag,hfrag,nbfrag,nhfrag
3153 use compare, only:freeres
3154 ! implicit real*8 (a-h,o-z)
3155 ! include 'DIMENSIONS'
3156 ! include 'DIMENSIONS.ZSCOPT'
3157 ! include 'COMMON.IOUNITS'
3158 ! include 'COMMON.FRAG'
3159 ! include 'COMMON.VAR'
3160 ! include 'COMMON.GEO'
3161 ! include 'COMMON.CHAIN'
3162 ! include 'COMMON.NAMES'
3163 ! include 'COMMON.INTERACT'
3164 integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres),&
3166 logical :: lprint,lprint_sec,not_done !el,freeres
3167 integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
3168 integer :: nstrand,nbeta,nhelix,iii1,jjj1
3169 real(kind=8) :: p1,p2
3170 !rel external freeres
3171 character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec))
3173 write (iout,*) "entered secondary2",ncont
3174 write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup
3176 write (iout,*) icont(1,i),icont(2,i)
3190 ! finding parallel beta
3191 !d write (iout,*) '------- looking for parallel beta -----------'
3197 if (i1.ge.nstart_sup .and. i1.le.nend_sup &
3198 .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
3199 !d write (iout,*) "parallel",i1,j1
3200 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
3203 !d write (iout,*) i1,j1
3209 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
3210 freeres(i1,j1,nsec,isec)) goto 5
3214 !d write (iout,*) i1,j1,not_done
3218 if (i1-ii1.gt.1) then
3222 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
3226 bfrag(1,nbfrag)=ii1+1
3227 bfrag(2,nbfrag)=i1+1
3228 bfrag(3,nbfrag)=jj1+1
3229 bfrag(4,nbfrag)=min0(j1+1,nres)
3233 isec(ij,nsec(ij))=nbeta
3237 isec(ij,nsec(ij))=nbeta
3242 if (nbeta.le.9) then
3243 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3244 "DefPropRes 'strand",nstrand,&
3245 "' 'num = ",ii1-1,"..",i1-1,"'"
3247 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3248 "DefPropRes 'strand",nstrand,&
3249 "' 'num = ",ii1-1,"..",i1-1,"'"
3252 if (nbeta.le.9) then
3253 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3254 "DefPropRes 'strand",nstrand,&
3255 "' 'num = ",jj1-1,"..",j1-1,"'"
3257 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3258 "DefPropRes 'strand",nstrand,&
3259 "' 'num = ",jj1-1,"..",j1-1,"'"
3261 write(12,'(a8,4i4)') &
3262 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
3266 endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
3269 ! finding antiparallel beta
3270 !d write (iout,*) '--------- looking for antiparallel beta ---------'
3275 if (freeres(i1,j1,nsec,isec)) then
3278 !d write (iout,*) i1,j1
3285 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
3286 freeres(i1,j1,nsec,isec)) goto 6
3290 !d write (iout,*) i1,j1,not_done
3294 if (i1-ii1.gt.1) then
3298 bfrag(2,nbfrag)=min0(i1+1,nres)
3299 bfrag(3,nbfrag)=min0(jj1+1,nres)
3306 if (nsec(ij).le.2) then
3307 isec(ij,nsec(ij))=nbeta
3313 if (nsec(ij).le.2) then
3314 isec(ij,nsec(ij))=nbeta
3319 if (lprint_sec) then
3320 write (iout,'(a,i3,4i4)')'antiparallel beta',&
3321 nbeta,ii1-1,i1,jj1,j1-1
3323 if (nstrand.le.9) then
3324 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3325 "DefPropRes 'strand",nstrand,&
3326 "' 'num = ",ii1-2,"..",i1-1,"'"
3328 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3329 "DefPropRes 'strand",nstrand,&
3330 "' 'num = ",ii1-2,"..",i1-1,"'"
3333 if (nstrand.le.9) then
3334 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3335 "DefPropRes 'strand",nstrand,&
3336 "' 'num = ",j1-2,"..",jj1-1,"'"
3338 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3339 "DefPropRes 'strand",nstrand,&
3340 "' 'num = ",j1-2,"..",jj1-1,"'"
3342 write(12,'(a8,4i4)') &
3343 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
3349 !d write (iout,*) "After beta:",nbfrag
3351 !d write (iout,*) (bfrag(j,i),j=1,4)
3354 if (nstrand.gt.0.and.lprint_sec) then
3355 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
3358 write(12,'(a9,i1,$)') " | strand",i
3360 write(12,'(a9,i2,$)') " | strand",i
3363 write(12,'(a1)') "'"
3367 ! finding alpha or 310 helix
3373 p1=phi(i1+2)*rad2deg
3375 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
3378 if (j1.eq.i1+3 .and. &
3379 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
3380 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
3381 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
3382 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
3385 if (nsec(ii1).eq.0) then
3394 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
3398 p1=phi(i1+2)*rad2deg
3399 p2=phi(j1+2)*rad2deg
3400 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
3403 !d write (iout,*) i1,j1,not_done,p1,p2
3406 if (j1-ii1.gt.4) then
3408 !d write (iout,*)'helix',nhelix,ii1,j1
3417 if (lprint_sec) then
3418 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
3419 if (nhelix.le.9) then
3420 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
3421 "DefPropRes 'helix",nhelix,&
3422 "' 'num = ",ii1-1,"..",j1-2,"'"
3424 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
3425 "DefPropRes 'helix",nhelix,&
3426 "' 'num = ",ii1-1,"..",j1-2,"'"
3433 if (nhelix.gt.0.and.lprint_sec) then
3434 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
3436 if (nhelix.le.9) then
3437 write(12,'(a8,i1,$)') " | helix",i
3439 write(12,'(a8,i2,$)') " | helix",i
3442 write(12,'(a1)') "'"
3445 if (lprint_sec) then
3446 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
3447 write(12,'(a20)') "XMacStand ribbon.mac"
3452 write(iout,*) 'UNRES seq:',anatemp
3454 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
3458 write(iout,*) 'helix ',(hfrag(i,j),i=1,2),anatemp
3464 do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j))
3467 do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j))
3472 do k=hfrag(1,j),hfrag(2,j)
3478 write (iout,*) "Secondary structure"
3483 write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
3484 write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien)
3485 write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien)
3490 end subroutine secondary2
3491 !-------------------------------------------------
3492 ! logical function freeres(i,j,nsec,isec)
3493 ! include 'DIMENSIONS'
3494 ! integer :: isec(nres,4),nsec(nres)
3495 ! integer :: i,j,k,l
3498 ! if (nsec(i).gt.1.or.nsec(j).gt.1) return
3501 ! if (isec(i,k).eq.isec(j,l)) return
3506 ! end function freeres
3507 !-------------------------------------------------
3508 subroutine alloc_compar_arrays(nfrg,nlev)
3510 use energy_data, only:maxcont
3512 integer :: nfrg,nlev
3514 !write(iout,*) "in alloc conpar arrays: nlevel=", nlevel," nfrag(1)=",nfrag(1)
3515 !------------------------
3518 allocate(nsccont_frag_ref(mmaxfrag)) !(mmaxfrag) !wham
3519 allocate(isccont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag) !wham
3520 !------------------------
3523 allocate(rmsfrag(nfrg,nlev+1),nc_fragm(nfrg,nlev+1)) !(maxfrag,maxlevel)
3524 allocate(qfrag(nfrg,2)) !(maxfrag,2)
3525 allocate(rmscutfrag(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3526 allocate(ang_cut(nfrg),ang_cut1(nfrg),frac_min(nfrg)) !(maxfrag)
3527 allocate(nc_req_setf(nfrg,nlev+1),npiece(nfrg,nlev+1),&
3528 ielecont(nfrg,nlev+1),isccont(nfrg,nlev+1),irms(nfrg,nlev+1),&
3529 ishifft(nfrg,nlev+1),len_frag(nfrg,nlev+1)) !(maxfrag,maxlevel)
3530 allocate(ncont_nat(2,nfrg,nlev+1))
3531 allocate(n_shift(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3532 ! allocate(nfrag(nlev)) !(maxlevel)
3533 allocate(isnfrag(nlev+2)) !(maxlevel+1)
3534 allocate(ifrag(2,maxpiece,nfrg)) !(2,maxpiece,maxfrag)
3535 allocate(ipiece(maxpiece,nfrg,2:nlev+1)) !(maxpiece,maxfrag,2:maxlevel)
3536 allocate(istruct(nfrg),iloc(nfrg),nlist_frag(nfrg)) !(maxfrag)
3537 allocate(iclass(nlev*nfrg,nlev+1)) !(maxlevel*maxfrag,maxlevel)
3538 allocate(list_frag(nres,nfrg)) !(maxres,maxfrag)
3539 !------------------------
3542 ! integer,dimension(:,:),allocatable :: icont_pept_ref !(2,maxcont)
3543 allocate(ncont_frag_ref(mmaxfrag)) !(mmaxfrag)
3544 allocate(icont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag)
3545 ! integer,dimension(:),allocatable :: isec_ref !(maxres)
3546 !------------------------
3547 ! module w_comm_local
3549 allocate(creff(3,2*nres),cc(3,2*nres)) !(3,nres*2)
3550 allocate(iadded(nres)) !(nres)
3551 allocate(inumber(2,nres)) !(2,nres)
3554 !-------------------------------------------------------------------------------
3555 end subroutine alloc_compar_arrays
3557 !-------------------------------------------------------------------------------
3558 end module conform_compar