2 !-----------------------------------------------------------------------------
8 #if .not. defined WHAM_RUN && .not. defined CLUSTER
14 use control, only: hpb_partition
16 use minimm, only: sc_move, minimize
19 !-----------------------------------------------------------------------------
22 !-----------------------------------------------------------------------------
24 #if .not. defined WHAM_RUN && .not. defined CLUSTER
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine contact(lprint,ncont,icont,co)
30 use geometry, only:dist
31 ! implicit real*8 (a-h,o-z)
32 ! include 'DIMENSIONS'
33 ! include 'COMMON.IOUNITS'
34 ! include 'COMMON.CHAIN'
35 ! include 'COMMON.INTERACT'
36 ! include 'COMMON.FFIELD'
37 ! include 'COMMON.NAMES'
38 real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
40 integer,dimension(2,12*nres) :: icont!(2,12*nres) !(2,maxcont) (maxcont=12*maxres)
43 real(kind=8) :: co,rcomp
44 integer :: kkk,i,j,i1,i2,it1,it2,iti,itj,inum,jnum
49 iti=iabs(itype(i,molnum(i)))
50 if (molnum(i).lt.3) then
57 itj=iabs(itype(j,molnum(i)))
58 if (molnum(j).lt.3) then
64 ! rcomp=sigmaii(iti,itj)+1.0D0
65 rcomp=facont*sigmaii(iti,itj)
67 ! rcomp=sigma(iti,itj)+1.0D0
68 rcomp=facont*sigma(iti,itj)
71 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
72 if (dist(inum,jnum).lt.rcomp) then
80 write (iout,'(a)') 'Contact map:'
86 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
87 i,restyp(it1,1),i1,restyp(it2,1),i2
92 co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
94 co = co / (nres*ncont)
96 end subroutine contact
97 !-----------------------------------------------------------------------------
98 real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref)
100 ! implicit real*8 (a-h,o-z)
101 ! include 'DIMENSIONS'
102 ! include 'COMMON.IOUNITS'
103 integer :: ncont,ncont_ref
104 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
106 integer :: i,j,nmatch
108 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
109 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
110 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
111 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
112 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
115 if (icont(1,i).eq.icont_ref(1,j) .and. &
116 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
119 ! print *,' nmatch=',nmatch
120 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
121 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
123 end function contact_fract
124 !-----------------------------------------------------------------------------
125 real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
127 ! implicit real*8 (a-h,o-z)
128 ! include 'DIMENSIONS'
129 ! include 'COMMON.IOUNITS'
130 integer :: ncont,ncont_ref
131 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
133 integer :: i,j,nmatch
135 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
136 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
137 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
138 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
139 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
142 if (icont(1,i).eq.icont_ref(1,j) .and. &
143 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
146 ! print *,' nmatch=',nmatch
147 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
148 contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
150 end function contact_fract_nn
151 !-----------------------------------------------------------------------------
152 subroutine hairpin(lprint,nharp,iharp)
154 use geometry, only:dist
155 ! implicit real*8 (a-h,o-z)
156 ! include 'DIMENSIONS'
157 ! include 'COMMON.IOUNITS'
158 ! include 'COMMON.CHAIN'
159 ! include 'COMMON.INTERACT'
160 ! include 'COMMON.FFIELD'
161 ! include 'COMMON.NAMES'
163 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
165 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
166 logical :: lprint,not_done
167 real(kind=8) :: rcomp=6.0d0
169 integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1
170 ! allocate(icont(2,12*nres))
174 ! print *,'nnt=',nnt,' nct=',nct
177 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
181 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
183 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
191 write (iout,'(a)') 'PP contact map:'
197 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
198 i,restyp(it1,1),i1,restyp(it2,1),i2
206 if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
207 ! write (iout,*) "found turn at ",i1,j1
215 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
219 ! write (iout,*) i1,j1,not_done
229 ! write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
234 ! write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
237 write (iout,*) "Hairpins:"
244 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
245 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
247 ! write (iout,'(a,i3,$)') restyp(itype(k,1)),k
252 end subroutine hairpin
253 !-----------------------------------------------------------------------------
255 !-----------------------------------------------------------------------------
256 subroutine elecont(lprint,ncont,icont)
258 ! implicit real*8 (a-h,o-z)
259 ! include 'DIMENSIONS'
260 ! include 'COMMON.IOUNITS'
261 ! include 'COMMON.CHAIN'
262 ! include 'COMMON.INTERACT'
263 ! include 'COMMON.LOCAL'
264 ! include 'COMMON.FFIELD'
265 ! include 'COMMON.NAMES'
267 real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_
268 real(kind=8) :: ael6_i,ael3_i
269 real(kind=8),dimension(2,2) :: app_,bpp_,rpp_
271 integer,dimension(2,12*nres) :: icont !(2,12*nres)(2,maxcont) (maxcont=12*maxres)
272 real(kind=8),dimension(12*nres) :: econt !(maxcont)
274 integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2
275 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw
276 real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb
277 real(kind=8) :: xmedi,ymedi,zmedi
278 real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij
279 real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,&
280 evdwij,el1,el2,eesij,ene
282 ! Load the constants of peptide bond - peptide bond interactions.
283 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
284 ! proline) - determined by averaging ECEPP energy.
288 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
289 data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
290 data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
291 data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
293 !el allocate(econt(12*nres)) !(maxcont)
296 elecutoff_14 = -0.5d0
297 if (lprint) write (iout,'(a)') &
298 "Constants of electrostatic interaction energy expression."
302 app_(i,j)=epp(i,j)*rri*rri
303 bpp_(i,j)=-2.0*epp(i,j)*rri
304 ael6_(i,j)=elpp_6(i,j)*4.2**6
305 ael3_(i,j)=elpp_3(i,j)*4.2**3
307 write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),&
314 print *, "nntt,nct",nnt,nct-2
316 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
327 if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
330 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
331 if (iteli.eq.2 .and. itelj.eq.2) goto 4
332 aaa=app_(iteli,itelj)
333 bbb=bpp_(iteli,itelj)
334 ael6_i=ael6_(iteli,itelj)
335 ael3_i=ael3_(iteli,itelj)
339 xj=c(1,j)+0.5*dxj-xmedi
340 yj=c(2,j)+0.5*dyj-ymedi
341 zj=c(3,j)+0.5*dzj-zmedi
342 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
347 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
348 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
349 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
350 fac=cosa-3.0*cosb*cosg
356 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
359 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
360 j.eq.i+2 .and. eesij.le.elecutoff_14) then
371 write (iout,*) 'Total average electrostatic energy: ',ees
372 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
374 write (iout,*) 'Electrostatic contacts before pruning: '
380 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
381 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
384 ! For given residues keep only the contacts with the greatest energy.
386 do while (i.lt.ncont)
392 do while (j.lt.ncont)
394 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
395 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
396 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
397 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
398 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
399 if (ic1.eq.icont(1,j)) then
401 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
402 .and. iabs(icont(1,k)-ic1).le.2 .and. &
403 econt(k).lt.econt(j) ) goto 21
405 else if (ic2.eq.icont(2,j) ) then
407 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
408 .and. iabs(icont(2,k)-ic2).le.2 .and. &
409 econt(k).lt.econt(j) ) goto 21
414 icont(1,k-1)=icont(1,k)
415 icont(2,k-1)=icont(2,k)
420 ! write (iout,*) "ncont",ncont
422 ! write (iout,*) icont(1,k),icont(2,k)
425 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
427 if (ic1.eq.icont(1,j)) then
429 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
430 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
431 econt(k).lt.econt(i) ) goto 21
433 else if (ic2.eq.icont(2,j) ) then
435 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
436 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
437 econt(k).lt.econt(i) ) goto 21
442 icont(1,k-1)=icont(1,k)
443 icont(2,k-1)=icont(2,k)
447 ! write (iout,*) "ncont",ncont
449 ! write (iout,*) icont(1,k),icont(2,k)
460 write (iout,*) 'Electrostatic contacts after pruning: '
466 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
467 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
471 end subroutine elecont
472 !-----------------------------------------------------------------------------
473 subroutine secondary2(lprint)
475 ! implicit real*8 (a-h,o-z)
476 ! include 'DIMENSIONS'
477 ! include 'COMMON.CHAIN'
478 ! include 'COMMON.IOUNITS'
479 ! include 'COMMON.DISTFIT'
480 ! include 'COMMON.VAR'
481 ! include 'COMMON.GEO'
482 ! include 'COMMON.CONTROL'
483 integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
485 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
486 integer,dimension(nres,0:4) :: isec !(maxres,4)
487 integer,dimension(nres) :: nsec !(maxres)
488 logical :: lprint,not_done !,freeres
489 real(kind=8) :: p1,p2
492 !el allocate(icont(2,12*nres),isec(nres,4),nsec(nres))
494 if(.not.dccart) call chainbuild_cart
495 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
496 !d call write_pdb(99,'sec structure',0d0)
506 call elecont(lprint,ncont,icont)
507 print *,"after elecont"
508 if (nres_molec(1).eq.0) return
510 ! finding parallel beta
511 !d write (iout,*) '------- looking for parallel beta -----------'
517 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
520 !d write (iout,*) i1,j1
526 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
527 freeres(i1,j1,nsec,isec)) goto 5
531 !d write (iout,*) i1,j1,not_done
535 if (i1-ii1.gt.1) then
539 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
543 bfrag(1,nbfrag)=ii1+1
545 bfrag(3,nbfrag)=jj1+1
546 bfrag(4,nbfrag)=min0(j1+1,nres)
550 isec(ij,nsec(ij))=nbeta
554 isec(ij,nsec(ij))=nbeta
560 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
561 "DefPropRes 'strand",nstrand,&
562 "' 'num = ",ii1-1,"..",i1-1,"'"
564 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
565 "DefPropRes 'strand",nstrand,&
566 "' 'num = ",ii1-1,"..",i1-1,"'"
570 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
571 "DefPropRes 'strand",nstrand,&
572 "' 'num = ",jj1-1,"..",j1-1,"'"
574 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
575 "DefPropRes 'strand",nstrand,&
576 "' 'num = ",jj1-1,"..",j1-1,"'"
578 write(12,'(a8,4i4)') &
579 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
585 ! finding alpha or 310 helix
592 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
595 if (j1.eq.i1+3 .and. &
596 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
597 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
598 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
599 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
602 if (nsec(ii1).eq.0) then
611 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
617 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
622 if (j1-ii1.gt.5) then
634 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
635 if (nhelix.le.9) then
636 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
637 "DefPropRes 'helix",nhelix,&
638 "' 'num = ",ii1-1,"..",j1-2,"'"
640 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
641 "DefPropRes 'helix",nhelix,&
642 "' 'num = ",ii1-1,"..",j1-2,"'"
648 if (nhelix.gt.0.and.lprint) then
649 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
651 if (nhelix.le.9) then
652 write(12,'(a8,i1,$)') " | helix",i
654 write(12,'(a8,i2,$)') " | helix",i
661 ! finding antiparallel beta
662 !d write (iout,*) '--------- looking for antiparallel beta ---------'
667 if (freeres(i1,j1,nsec,isec)) then
670 !d write (iout,*) i1,j1
677 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
678 freeres(i1,j1,nsec,isec)) goto 6
682 !d write (iout,*) i1,j1,not_done
686 if (i1-ii1.gt.1) then
690 bfrag(2,nbfrag)=min0(i1+1,nres)
691 bfrag(3,nbfrag)=min0(jj1+1,nres)
698 if (nsec(ij).le.2) then
699 isec(ij,nsec(ij))=nbeta
705 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
706 isec(ij,nsec(ij))=nbeta
712 write (iout,'(a,i3,4i4)')'antiparallel beta',&
713 nbeta,ii1-1,i1,jj1,j1-1
715 if (nstrand.le.9) then
716 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
717 "DefPropRes 'strand",nstrand,&
718 "' 'num = ",ii1-2,"..",i1-1,"'"
720 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
721 "DefPropRes 'strand",nstrand,&
722 "' 'num = ",ii1-2,"..",i1-1,"'"
725 if (nstrand.le.9) then
726 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
727 "DefPropRes 'strand",nstrand,&
728 "' 'num = ",j1-2,"..",jj1-1,"'"
730 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
731 "DefPropRes 'strand",nstrand,&
732 "' 'num = ",j1-2,"..",jj1-1,"'"
734 write(12,'(a8,4i4)') &
735 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
741 if (nstrand.gt.0.and.lprint) then
742 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
745 write(12,'(a9,i1,$)') " | strand",i
747 write(12,'(a9,i2,$)') " | strand",i
756 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
757 write(12,'(a20)') "XMacStand ribbon.mac"
759 if (nres_molec(1).eq.0) return
760 write(iout,*) 'UNRES seq:'
762 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
766 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
771 end subroutine secondary2
773 !-----------------------------------------------------------------------------
774 logical function freeres(i,j,nsec,isec)
776 ! implicit real*8 (a-h,o-z)
777 ! include 'DIMENSIONS'
778 integer,dimension(nres,4) :: isec !(maxres,4)
779 integer,dimension(nres) :: nsec !(maxres)
786 if (nsec(i).lt.0.or.nsec(j).lt.0) return
788 if (nsec(i).gt.1.or.nsec(j).gt.1) return
791 if (isec(i,k).eq.isec(j,l)) return
797 !-----------------------------------------------------------------------------
799 !-----------------------------------------------------------------------------
800 logical function seq_comp(itypea,itypeb,length)
803 integer :: length,itypea(length),itypeb(length)
806 if (itypea(i).ne.itypeb(i)) then
813 end function seq_comp
815 !-----------------------------------------------------------------------------
817 !-----------------------------------------------------------------------------
818 subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
820 ! implicit real*8 (a-h,o-z)
821 ! include 'DIMENSIONS'
822 ! include 'COMMON.CHAIN'
823 ! include 'COMMON.CONTACTS'
824 ! include 'COMMON.IOUNITS'
825 real(kind=8) :: przes(3),obr(3,3)
826 logical :: non_conv,lprn
827 real(kind=8) :: rms,frac,frac_nn,co
828 ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
832 ! print *,"before contact"
833 !elte(iout,*) "rms_nacc before contact"
834 call contact(.false.,ncont,icont,co)
835 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
836 frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
837 if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
838 'RMS deviation from the reference structure:',rms,&
839 ' % of native contacts:',frac*100,&
840 ' % of nonnative contacts:',frac_nn*100,&
844 end subroutine rms_nac_nnc
845 !-----------------------------------------------------------------------------
846 subroutine rmsd(drms)
848 use regularize_, only:fitsq
849 ! implicit real*8 (a-h,o-z)
850 ! include 'DIMENSIONS'
854 ! include 'COMMON.CHAIN'
855 ! include 'COMMON.IOUNITS'
856 ! include 'COMMON.INTERACT'
857 ! include 'COMMON.CONTROL'
859 real(kind=8) :: przes(3),obrot(3,3)
860 real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
863 real(kind=8) :: drms,rminroz,roznica
864 integer :: i,j,iatom,kkk,iti,k
866 !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
874 ! print *,"nz_start",nz_start," nz_end",nz_end
875 ! if (symetr.le.1) then
877 ! do i=nz_start,nz_end
881 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
882 ! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
884 ! if (iz_sc.eq.1.and.iti.ne.10) then
887 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
888 ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
895 print *,nz_start,nz_end,nstart_seq-nstart_sup
900 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
901 crefcopy(k,iatom)=cref(k,i,kkk)
903 if (iz_sc.eq.1.and.iti.ne.10) then
906 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
907 crefcopy(k,iatom)=cref(k,nres+i,kkk)
916 ! write (iout,*) 'Ccopy and CREFcopy'
917 ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
918 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
919 ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
920 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
922 ! ----- end diagnostics
924 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
925 przes,obrot,non_conv)
927 print *,'Problems in FITSQ!!! rmsd'
928 write (iout,*) 'Problems in FITSQ!!! rmsd'
929 print *,'Ccopy and CREFcopy'
930 write (iout,*) 'Ccopy and CREFcopy'
931 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
932 (crefcopy(j,k),j=1,3),k=1,iatom)
933 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
934 (crefcopy(j,k),j=1,3),k=1,iatom)
936 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
942 ! write (iout,*) "roznica", roznica,kkk
943 if (roznica.le.rminroz) rminroz=roznica
945 drms=dsqrt(dabs(rminroz))
947 ! write (iout,*) "nperm,symetr", nperm,symetr
948 ! ---- end diagnostics
951 !-----------------------------------------------------------------------------
952 subroutine rmsd_csa(drms)
954 use regularize_, only:fitsq
955 ! implicit real*8 (a-h,o-z)
956 ! include 'DIMENSIONS'
960 ! include 'COMMON.CHAIN'
961 ! include 'COMMON.IOUNITS'
962 ! include 'COMMON.INTERACT'
964 real(kind=8) :: przes(3),obrot(3,3)
965 real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
966 integer :: kkk,iatom,ierror,ierrcode
970 real(kind=8) :: drms,roznica
972 allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
980 ccopy(k,iatom)=c(k,i)
981 crefcopy(k,iatom)=crefjlee(k,i)
983 if (iz_sc.eq.1.and.iti.ne.10) then
986 ccopy(k,iatom)=c(k,nres+i)
987 crefcopy(k,iatom)=crefjlee(k,nres+i)
992 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
993 przes,obrot,non_conv)
995 print *,'Problems in FITSQ!!! rmsd_csa'
996 write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
997 print *,'Ccopy and CREFcopy'
998 write (iout,*) 'Ccopy and CREFcopy'
999 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
1000 (crefcopy(j,k),j=1,3),k=1,iatom)
1001 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
1002 (crefcopy(j,k),j=1,3),k=1,iatom)
1004 call mpi_abort(mpi_comm_world,ierror,ierrcode)
1009 drms=dsqrt(dabs(roznica))
1011 end subroutine rmsd_csa
1012 !-----------------------------------------------------------------------------
1014 !-----------------------------------------------------------------------------
1018 use geometry, only:pinorm
1019 use random, only:ran_number,iran_num
1020 ! implicit real*8 (a-h,o-z)
1021 ! include 'DIMENSIONS'
1023 ! include 'COMMON.GEO'
1024 ! include 'COMMON.VAR'
1025 ! include 'COMMON.INTERACT'
1026 ! include 'COMMON.IOUNITS'
1027 ! include 'COMMON.DISTFIT'
1028 ! include 'COMMON.SBRIDGE'
1029 ! include 'COMMON.CONTROL'
1030 ! include 'COMMON.FFIELD'
1031 ! include 'COMMON.MINIM'
1032 ! include 'COMMON.CHAIN'
1033 real(kind=8) :: time0,time1
1034 real(kind=8) :: energy(0:n_ene),ee
1035 real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres)
1036 integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1037 logical :: debug,accepted
1038 real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1041 !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1043 call geom_to_var(nvar,var1)
1048 write(iout,*) 'etot=',0,etot,rms
1049 call secondary2(.false.)
1051 call write_pdb(0,'first structure',etot)
1060 betbol=1.0D0/(1.9858D-3*temp)
1062 d=ran_number(-pi,pi)
1063 ! phi(jr)=pinorm(phi(jr)+d)
1068 write(iout,*) 'etot=',1,etot0,rms
1069 call write_pdb(1,'perturb structure',etot0)
1073 d=ran_number(-da,da)
1075 phi(jr)=pinorm(phi(jr)+d)
1080 if (etot.lt.etot0) then
1084 xxr=ran_number(0.0D0,1.0D0)
1085 xxh=betbol*(etot-etot0)
1086 if (xxh.lt.50.0D0) then
1088 if (xxh.gt.xxr) accepted=.true.
1092 ! print *,etot0,etot,accepted
1096 write(iout,*) 'etot=',i,etot,rms
1097 call write_pdb(i,'MC structure',etot)
1099 ! call geom_to_var(nvar,var1)
1100 call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1101 call geom_to_var(nvar,var)
1102 call minimize(etot,var,iretcode,nfun)
1103 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1104 call var_to_geom(nvar,var)
1107 write(iout,*) 'etot mcm=',i,etot,rms
1108 call write_pdb(i+1,'MCM structure',etot)
1109 call var_to_geom(nvar,var1)
1117 ! call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1118 ! call geom_to_var(nvar,var)
1121 ! call write_pdb(998 ,'sc min',etot)
1123 ! call minimize(etot,var,iretcode,nfun)
1124 ! write(iout,*)'------------------------------------------------'
1125 ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1127 ! call var_to_geom(nvar,var)
1129 ! call write_pdb(999,'full min',etot)
1133 !-----------------------------------------------------------------------------
1138 ! implicit real*8 (a-h,o-z)
1139 ! include 'DIMENSIONS'
1141 ! include 'COMMON.GEO'
1142 ! include 'COMMON.VAR'
1143 ! include 'COMMON.INTERACT'
1144 ! include 'COMMON.IOUNITS'
1145 ! include 'COMMON.DISTFIT'
1146 ! include 'COMMON.SBRIDGE'
1147 ! include 'COMMON.CONTROL'
1148 ! include 'COMMON.FFIELD'
1149 ! include 'COMMON.MINIM'
1150 ! include 'COMMON.CHAIN'
1151 real(kind=8) :: time0,time1
1152 real(kind=8) :: energy(0:n_ene),ee
1153 real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1157 integer :: i,ij,ieval,iretcode,nfun
1158 real(kind=8) :: etot
1160 allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1162 call geom_to_var(nvar,var1)
1166 write(iout,*) nnt,nct,etot
1167 call write_pdb(1,'first structure',etot)
1168 call secondary2(.true.)
1177 call var_to_geom(nvar,var1)
1178 write(iout,*) 'N16 test',(jdata(i),i=1,5)
1179 call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1181 call geom_to_var(nvar,var)
1187 call minimize(etot,var,iretcode,nfun)
1188 write(iout,*)'------------------------------------------------'
1189 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1195 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1196 nfun/(time1-time0),' eval/s'
1198 call var_to_geom(nvar,var)
1200 call write_pdb(ij*100+99,'full min',etot)
1207 end subroutine test_n16
1209 !-----------------------------------------------------------------------------
1210 subroutine test_local
1212 ! implicit real*8 (a-h,o-z)
1213 ! include 'DIMENSIONS'
1214 ! include 'COMMON.GEO'
1215 ! include 'COMMON.VAR'
1216 ! include 'COMMON.INTERACT'
1217 ! include 'COMMON.IOUNITS'
1218 real(kind=8) :: time0,time1
1219 real(kind=8) :: energy(0:n_ene),ee
1220 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1222 real(kind=8) :: etot
1224 ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres)
1226 ! call geom_to_var(nvar,varia)
1227 call write_pdb(1,'first structure',0d0)
1231 write(iout,*) nnt,nct,etot
1233 write(iout,*) 'calling sc_move'
1234 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1235 write(iout,*) nft_sc,etot
1236 call write_pdb(2,'second structure',etot)
1238 write(iout,*) 'calling local_move'
1239 call local_move_init(.false.)
1240 call local_move(24,29,20d0,50d0)
1242 call write_pdb(3,'third structure',etot)
1244 write(iout,*) 'calling sc_move'
1245 call sc_move(24,29,5,10d0,nft_sc,etot)
1246 write(iout,*) nft_sc,etot
1247 call write_pdb(2,'last structure',etot)
1250 end subroutine test_local
1251 !-----------------------------------------------------------------------------
1254 ! implicit real*8 (a-h,o-z)
1255 ! include 'DIMENSIONS'
1256 ! include 'COMMON.GEO'
1257 ! include 'COMMON.VAR'
1258 ! include 'COMMON.INTERACT'
1259 ! include 'COMMON.IOUNITS'
1260 real(kind=8) :: time0,time1,etot
1261 real(kind=8) :: energy(0:n_ene),ee
1262 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1266 ! call geom_to_var(nvar,varia)
1267 call write_pdb(1,'first structure',0d0)
1271 write(iout,*) nnt,nct,etot
1273 write(iout,*) 'calling sc_move'
1275 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1276 write(iout,*) nft_sc,etot
1277 call write_pdb(2,'second structure',etot)
1279 write(iout,*) 'calling sc_move 2nd time'
1281 call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1282 write(iout,*) nft_sc,etot
1283 call write_pdb(3,'last structure',etot)
1285 end subroutine test_sc
1286 !-----------------------------------------------------------------------------
1287 subroutine bgrow(bstrand,nbstrand,in,ind,new)
1289 ! implicit real*8 (a-h,o-z)
1290 ! include 'DIMENSIONS'
1291 ! include 'COMMON.CHAIN'
1292 integer,dimension(nres/3,6) :: bstrand !(maxres/3,6)
1295 integer :: nbstrand,in,ind,new,ishift,i
1297 ishift=iabs(bstrand(in,ind+4)-new)
1299 print *,'bgrow',bstrand(in,ind+4),new,ishift
1304 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1306 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1307 if (bstrand(i,5).lt.bstrand(i,6)) then
1308 bstrand(i,5)=bstrand(i,5)-ishift
1310 bstrand(i,5)=bstrand(i,5)+ishift
1315 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1317 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1318 if (bstrand(i,6).lt.bstrand(i,5)) then
1319 bstrand(i,6)=bstrand(i,6)-ishift
1321 bstrand(i,6)=bstrand(i,6)+ishift
1328 end subroutine bgrow
1329 !-----------------------------------------------------------------------------
1332 use geometry, only:dist
1334 ! implicit real*8 (a-h,o-z)
1335 ! include 'DIMENSIONS'
1337 ! include 'COMMON.GEO'
1338 ! include 'COMMON.CHAIN'
1339 ! include 'COMMON.IOUNITS'
1340 ! include 'COMMON.VAR'
1341 ! include 'COMMON.CONTROL'
1342 ! include 'COMMON.SBRIDGE'
1343 ! include 'COMMON.FFIELD'
1344 ! include 'COMMON.MINIM'
1346 ! include 'COMMON.DISTFIT'
1347 integer :: if(20,nres),nif,ifa(20)
1348 integer :: ibc(0:nres,0:nres),istrand(20)
1349 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1350 integer :: itmp(20,nres)
1351 real(kind=8) :: time0,time1
1352 real(kind=8) :: energy(0:n_ene),ee
1353 real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres)
1355 logical :: debug,ltest,usedbfrag(nres/3)
1356 character(len=50) :: linia
1358 integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1359 integer :: bstrand(nres/3,6),nbstrand
1360 real(kind=8) :: etot
1361 integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1362 in,ind,ifun,nfun,iretcode
1363 !------------------------
1366 !------------------------
1373 call geom_to_var(nvar,vorg)
1374 call secondary2(debug)
1376 if (nbfrag.le.1) return
1379 usedbfrag(i)=.false.
1383 nbetasheet=nbetasheet+1
1385 bstrand(1,1)=bfrag(1,1)
1386 bstrand(1,2)=bfrag(2,1)
1387 bstrand(1,3)=nbetasheet
1389 bstrand(1,5)=bfrag(1,1)
1390 bstrand(1,6)=bfrag(2,1)
1391 do i=bfrag(1,1),bfrag(2,1)
1392 betasheet(i)=nbetasheet
1396 bstrand(2,1)=bfrag(3,1)
1397 bstrand(2,2)=bfrag(4,1)
1398 bstrand(2,3)=nbetasheet
1399 bstrand(2,5)=bfrag(3,1)
1400 bstrand(2,6)=bfrag(4,1)
1402 if (bfrag(3,1).le.bfrag(4,1)) then
1404 do i=bfrag(3,1),bfrag(4,1)
1405 betasheet(i)=nbetasheet
1410 do i=bfrag(4,1),bfrag(3,1)
1411 betasheet(i)=nbetasheet
1418 do while (iused_nbfrag.ne.nbfrag)
1422 IF (.not.usedbfrag(j)) THEN
1424 write (*,*) j,(bfrag(i,j),i=1,4)
1426 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1428 write (*,*) '------------------'
1431 if (bfrag(3,j).le.bfrag(4,j)) then
1432 do i=bfrag(3,j),bfrag(4,j)
1433 if(betasheet(i).eq.nbetasheet) then
1435 do k=bfrag(3,j),bfrag(4,j)
1436 betasheet(k)=nbetasheet
1441 iused_nbfrag=iused_nbfrag+1
1442 do k=bfrag(1,j),bfrag(2,j)
1443 betasheet(k)=nbetasheet
1444 ibetasheet(k)=nbstrand
1446 if (bstrand(in,4).lt.0) then
1447 bstrand(nbstrand,1)=bfrag(2,j)
1448 bstrand(nbstrand,2)=bfrag(1,j)
1449 bstrand(nbstrand,3)=nbetasheet
1450 bstrand(nbstrand,4)=-nbstrand
1451 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1452 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1453 if(bstrand(in,1).lt.bfrag(4,j)) then
1454 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1456 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1457 (bstrand(in,5)-bfrag(4,j))
1459 if(bstrand(in,2).gt.bfrag(3,j)) then
1460 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1462 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1463 (-bstrand(in,6)+bfrag(3,j))
1466 bstrand(nbstrand,1)=bfrag(1,j)
1467 bstrand(nbstrand,2)=bfrag(2,j)
1468 bstrand(nbstrand,3)=nbetasheet
1469 bstrand(nbstrand,4)=nbstrand
1470 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1471 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1472 if(bstrand(in,1).gt.bfrag(3,j)) then
1473 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1475 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1476 (-bstrand(in,5)+bfrag(3,j))
1478 if(bstrand(in,2).lt.bfrag(4,j)) then
1479 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1481 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1482 (bstrand(in,6)-bfrag(4,j))
1487 if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1488 in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1489 do k=bfrag(1,j),bfrag(2,j)
1490 betasheet(k)=nbetasheet
1495 iused_nbfrag=iused_nbfrag+1
1496 do k=bfrag(3,1),bfrag(4,1)
1497 betasheet(k)=nbetasheet
1498 ibetasheet(k)=nbstrand
1500 if (bstrand(in,4).lt.0) then
1501 bstrand(nbstrand,1)=bfrag(4,j)
1502 bstrand(nbstrand,2)=bfrag(3,j)
1503 bstrand(nbstrand,3)=nbetasheet
1504 bstrand(nbstrand,4)=-nbstrand
1505 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1506 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1507 if(bstrand(in,1).lt.bfrag(2,j)) then
1508 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1510 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1511 (bstrand(in,5)-bfrag(2,j))
1513 if(bstrand(in,2).gt.bfrag(1,j)) then
1514 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1516 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1517 (-bstrand(in,6)+bfrag(1,j))
1520 bstrand(nbstrand,1)=bfrag(3,j)
1521 bstrand(nbstrand,2)=bfrag(4,j)
1522 bstrand(nbstrand,3)=nbetasheet
1523 bstrand(nbstrand,4)=nbstrand
1524 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1525 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1526 if(bstrand(in,1).gt.bfrag(1,j)) then
1527 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1529 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1530 (-bstrand(in,5)+bfrag(1,j))
1532 if(bstrand(in,2).lt.bfrag(2,j)) then
1533 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1535 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1536 (bstrand(in,6)-bfrag(2,j))
1543 do i=bfrag(4,j),bfrag(3,j)
1544 if(betasheet(i).eq.nbetasheet) then
1546 do k=bfrag(4,j),bfrag(3,j)
1547 betasheet(k)=nbetasheet
1552 iused_nbfrag=iused_nbfrag+1
1553 do k=bfrag(1,j),bfrag(2,j)
1554 betasheet(k)=nbetasheet
1555 ibetasheet(k)=nbstrand
1557 if (bstrand(in,4).lt.0) then
1558 bstrand(nbstrand,1)=bfrag(1,j)
1559 bstrand(nbstrand,2)=bfrag(2,j)
1560 bstrand(nbstrand,3)=nbetasheet
1561 bstrand(nbstrand,4)=nbstrand
1562 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1563 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1564 if(bstrand(in,1).lt.bfrag(3,j)) then
1565 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1567 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1568 (bstrand(in,5)-bfrag(3,j))
1570 if(bstrand(in,2).gt.bfrag(4,j)) then
1571 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1573 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1574 (-bstrand(in,6)+bfrag(4,j))
1577 bstrand(nbstrand,1)=bfrag(2,j)
1578 bstrand(nbstrand,2)=bfrag(1,j)
1579 bstrand(nbstrand,3)=nbetasheet
1580 bstrand(nbstrand,4)=-nbstrand
1581 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1582 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1583 if(bstrand(in,1).gt.bfrag(4,j)) then
1584 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1586 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1587 (-bstrand(in,5)+bfrag(4,j))
1589 if(bstrand(in,2).lt.bfrag(3,j)) then
1590 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1592 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1593 (bstrand(in,6)-bfrag(3,j))
1598 if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1599 in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1600 do k=bfrag(1,j),bfrag(2,j)
1601 betasheet(k)=nbetasheet
1606 iused_nbfrag=iused_nbfrag+1
1607 do k=bfrag(4,j),bfrag(3,j)
1608 betasheet(k)=nbetasheet
1609 ibetasheet(k)=nbstrand
1611 if (bstrand(in,4).lt.0) then
1612 bstrand(nbstrand,1)=bfrag(4,j)
1613 bstrand(nbstrand,2)=bfrag(3,j)
1614 bstrand(nbstrand,3)=nbetasheet
1615 bstrand(nbstrand,4)=nbstrand
1616 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1617 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1618 if(bstrand(in,1).lt.bfrag(2,j)) then
1619 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1621 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1622 (bstrand(in,5)-bfrag(2,j))
1624 if(bstrand(in,2).gt.bfrag(1,j)) then
1625 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1627 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1628 (-bstrand(in,6)+bfrag(1,j))
1631 bstrand(nbstrand,1)=bfrag(3,j)
1632 bstrand(nbstrand,2)=bfrag(4,j)
1633 bstrand(nbstrand,3)=nbetasheet
1634 bstrand(nbstrand,4)=-nbstrand
1635 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1636 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1637 if(bstrand(in,1).gt.bfrag(1,j)) then
1638 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1640 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1641 (-bstrand(in,5)+bfrag(1,j))
1643 if(bstrand(in,2).lt.bfrag(2,j)) then
1644 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1646 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1647 (bstrand(in,6)-bfrag(2,j))
1661 do while (usedbfrag(j))
1666 nbetasheet=nbetasheet+1
1667 bstrand(nbstrand,1)=bfrag(1,j)
1668 bstrand(nbstrand,2)=bfrag(2,j)
1669 bstrand(nbstrand,3)=nbetasheet
1670 bstrand(nbstrand,5)=bfrag(1,j)
1671 bstrand(nbstrand,6)=bfrag(2,j)
1673 bstrand(nbstrand,4)=nbstrand
1674 do i=bfrag(1,j),bfrag(2,j)
1675 betasheet(i)=nbetasheet
1676 ibetasheet(i)=nbstrand
1680 bstrand(nbstrand,1)=bfrag(3,j)
1681 bstrand(nbstrand,2)=bfrag(4,j)
1682 bstrand(nbstrand,3)=nbetasheet
1683 bstrand(nbstrand,5)=bfrag(3,j)
1684 bstrand(nbstrand,6)=bfrag(4,j)
1686 if (bfrag(3,j).le.bfrag(4,j)) then
1687 bstrand(nbstrand,4)=nbstrand
1688 do i=bfrag(3,j),bfrag(4,j)
1689 betasheet(i)=nbetasheet
1690 ibetasheet(i)=nbstrand
1693 bstrand(nbstrand,4)=-nbstrand
1694 do i=bfrag(4,j),bfrag(3,j)
1695 betasheet(i)=nbetasheet
1696 ibetasheet(i)=nbstrand
1700 iused_nbfrag=iused_nbfrag+1
1706 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1713 if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1717 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1720 !------------------------
1724 if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1725 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1727 ifb(nifb,1)=bstrand(i,4)
1728 ifb(nifb,2)=bstrand(j,4)
1735 write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1741 write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1743 nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1745 if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1746 nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1752 if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1753 if (if(j,i).gt.0) then
1754 if(betasheet(if(j,i)).eq.0 .or. &
1755 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1760 write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1763 ! read (inp,*) (ifa(i),i=1,4)
1765 ! read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1769 !------------------------
1774 !ccccccccccccccccccccccccccccccccc
1776 !ccccccccccccccccccccccccccccccccc
1780 istrand(is-j+1)=int(ii/is**(is-j))
1781 ii=ii-istrand(is-j+1)*is**(is-j)
1785 istrand(k)=istrand(k)+1
1786 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1790 if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1791 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1800 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1801 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1802 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1803 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1809 if (mod(isa,2).eq.0) then
1811 if (istrand(k).eq.1) ltest=.false.
1814 do k=(isa+1)/2+1,isa
1815 if (istrand(k).eq.1) ltest=.false.
1819 IF (ltest.and.lifb0.eq.1) THEN
1822 call var_to_geom(nvar,vorg)
1824 write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1825 write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1826 write (linia,'(10i3)') (istrand(k),k=1,isa)
1836 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1838 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1842 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1848 write(*,*) (itmp(j,i),j=1,4)
1852 ! ifa(1),ifa(2),ifa(3),ifa(4)
1853 ! if(1,i),if(2,i),if(3,i),if(4,i)
1858 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1859 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1860 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1861 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1869 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1871 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1875 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1877 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1880 !------------------------
1883 ! freeze sec.elements
1893 do i=bfrag(1,j),bfrag(2,j)
1898 if (bfrag(3,j).le.bfrag(4,j)) then
1899 do i=bfrag(3,j),bfrag(4,j)
1905 do i=bfrag(4,j),bfrag(3,j)
1913 do i=hfrag(1,j),hfrag(2,j)
1921 !------------------------
1922 ! generate constrains
1930 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1938 else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1946 else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1954 else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1962 else if ( ibc(i,j).gt.0 ) then
1963 d0(ind)=DIST(i,ibc(i,j))
1970 else if ( ibc(j,i).gt.0 ) then
1971 d0(ind)=DIST(ibc(j,i),j)
1985 !d--------------------------
1987 write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1988 ibc(jhpb(i),ihpb(i)),' --',&
1989 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1996 call contact_cp_min(varia,ifun,iconf,linia,debug)
2001 call minimize(etot,varia,iretcode,nfun)
2002 write(iout,*)'------------------------------------------------'
2003 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2009 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2010 nfun/(time1-time0),' eval/s'
2012 write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
2013 call var_to_geom(nvar,varia)
2015 call write_pdb(900+iconf,linia,etot)
2020 call enerprint(energy)
2022 !d call briefout(0,etot)
2023 !d call secondary2(.true.)
2027 !ccccccccccccccccccccccccccccccccccc
2030 !ccccccccccccccccccccccccccccccccccc
2033 10 write (iout,'(a)') 'Error reading test structure.'
2035 end subroutine test11
2036 !-----------------------------------------------------------------------------
2039 use geometry, only:dist
2041 ! implicit real*8 (a-h,o-z)
2042 ! include 'DIMENSIONS'
2044 ! include 'COMMON.GEO'
2045 ! include 'COMMON.CHAIN'
2046 ! include 'COMMON.IOUNITS'
2047 ! include 'COMMON.VAR'
2048 ! include 'COMMON.CONTROL'
2049 ! include 'COMMON.SBRIDGE'
2050 ! include 'COMMON.FFIELD'
2051 ! include 'COMMON.MINIM'
2053 ! include 'COMMON.DISTFIT'
2054 integer :: if(3,nres),nif
2055 integer :: ibc(nres,nres),istrand(20)
2056 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2057 real(kind=8) :: time0,time1
2058 real(kind=8) :: energy(0:n_ene),ee
2059 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2061 logical :: debug,ltest
2062 character(len=50) :: linia
2063 integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2064 real(kind=8) :: etot
2067 read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2070 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2074 !------------------------
2075 call secondary2(debug)
2076 !------------------------
2084 ! freeze sec.elements and store indexes for beta constrains
2094 do i=bfrag(1,j),bfrag(2,j)
2099 if (bfrag(3,j).le.bfrag(4,j)) then
2100 do i=bfrag(3,j),bfrag(4,j)
2104 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2107 do i=bfrag(4,j),bfrag(3,j)
2111 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2116 do i=hfrag(1,j),hfrag(2,j)
2125 ! ---------------- test --------------
2127 if (ibc(if(1,i),if(2,i)).eq.-1) then
2128 ibc(if(1,i),if(2,i))=if(3,i)
2129 ibc(if(1,i),if(3,i))=if(2,i)
2130 else if (ibc(if(2,i),if(1,i)).eq.-1) then
2131 ibc(if(2,i),if(1,i))=0
2132 ibc(if(1,i),if(2,i))=if(3,i)
2133 ibc(if(1,i),if(3,i))=if(2,i)
2135 ibc(if(1,i),if(2,i))=if(3,i)
2136 ibc(if(1,i),if(3,i))=if(2,i)
2142 if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j)
2145 !------------------------
2151 if ( ibc(i,j).eq.-1 ) then
2159 else if ( ibc(i,j).gt.0 ) then
2160 d0(ind)=DIST(i,ibc(i,j))
2167 else if ( ibc(j,i).gt.0 ) then
2168 d0(ind)=DIST(ibc(j,i),j)
2182 !d--------------------------
2183 write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2184 ibc(jhpb(i),ihpb(i)),' --',&
2185 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2193 call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2198 call minimize(etot,varia,iretcode,nfun)
2199 write(iout,*)'------------------------------------------------'
2200 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2206 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2207 nfun/(time1-time0),' eval/s'
2210 call var_to_geom(nvar,varia)
2212 call write_pdb(999,'full min',etot)
2217 call enerprint(energy)
2219 call briefout(0,etot)
2220 call secondary2(.true.)
2223 10 write (iout,'(a)') 'Error reading test structure.'
2225 end subroutine test3
2226 !-----------------------------------------------------------------------------
2230 ! implicit real*8 (a-h,o-z)
2231 ! include 'DIMENSIONS'
2233 ! include 'COMMON.GEO'
2234 ! include 'COMMON.CHAIN'
2235 ! include 'COMMON.IOUNITS'
2236 ! include 'COMMON.VAR'
2237 ! include 'COMMON.CONTROL'
2238 ! include 'COMMON.SBRIDGE'
2239 ! include 'COMMON.FFIELD'
2240 ! include 'COMMON.MINIM'
2242 ! include 'COMMON.DISTFIT'
2243 integer :: if(2,2),ind
2244 integer :: iff(nres)
2245 real(kind=8) :: time0,time1
2246 real(kind=8) :: energy(0:n_ene),ee
2247 real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2248 theta1,phi1,alph1,omeg1 !(maxres)
2249 real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres)
2251 integer :: i,j,nn,ifun,iretcode,nfun
2252 real(kind=8) :: etot
2255 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2256 write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2257 read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2258 read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2259 read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2260 read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2262 theta2(i)=deg2rad*theta2(i)
2263 phi2(i)=deg2rad*phi2(i)
2264 alph2(i)=deg2rad*alph2(i)
2265 omeg2(i)=deg2rad*omeg2(i)
2279 !------------------------
2284 do i=if(j,1),if(j,2)
2290 call geom_to_var(nvar,varia)
2291 call write_pdb(1,'first structure',0d0)
2293 call secondary(.true.)
2295 call secondary2(.true.)
2298 if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2299 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2300 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2303 if (bfrag(3,j).lt.bfrag(4,j)) then
2304 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2305 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2306 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2308 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2309 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2310 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2323 call geom_to_var(nvar,varia2)
2324 call write_pdb(2,'second structure',0d0)
2328 !-------------------------------------------------------
2331 call contact_cp(varia,varia2,iff,ifun,7)
2336 call minimize(etot,varia,iretcode,nfun)
2337 write(iout,*)'------------------------------------------------'
2338 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2344 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2345 nfun/(time1-time0),' eval/s'
2348 call var_to_geom(nvar,varia)
2350 call write_pdb(999,'full min',etot)
2355 call enerprint(energy)
2357 call briefout(0,etot)
2360 10 write (iout,'(a)') 'Error reading test structure.'
2362 end subroutine test__
2363 !-----------------------------------------------------------------------------
2364 subroutine secondary(lprint)
2366 ! implicit real*8 (a-h,o-z)
2367 ! include 'DIMENSIONS'
2368 ! include 'COMMON.CHAIN'
2369 ! include 'COMMON.IOUNITS'
2370 ! include 'COMMON.DISTFIT'
2372 integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2373 logical :: lprint,not_done
2374 real(kind=4) :: dcont(nres*nres/2),d
2375 real(kind=4) :: rcomp = 7.0
2376 real(kind=4) :: rbeta = 5.2
2377 real(kind=4) :: ralfa = 5.2
2378 real(kind=4) :: r310 = 6.6
2379 real(kind=8),dimension(3) :: xpi,xpj
2380 integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2382 call chainbuild_cart
2383 !d call write_pdb(99,'sec structure',0d0)
2395 xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2399 xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2401 !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2402 !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2403 !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2404 !d print *,'CA',i,j,d
2405 d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2406 (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2407 (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2408 if ( d.lt.rcomp*rcomp) then
2412 dcont(ncont)=sqrt(d)
2418 write (iout,'(a)') '#PP contact map distances:'
2420 write (iout,'(3i4,f10.5)') &
2421 i,icont(1,i),icont(2,i),dcont(i)
2425 ! finding parallel beta
2426 !d write (iout,*) '------- looking for parallel beta -----------'
2432 if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2433 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2434 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2435 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2436 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2437 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2441 !d write (iout,*) i1,j1,dcont(i)
2447 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2448 .and. dcont(j).le.rbeta .and. &
2449 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2450 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2451 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2452 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2453 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2458 !d write (iout,*) i1,j1,dcont(j),not_done
2462 if (i1-ii1.gt.1) then
2466 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2475 isec(ij,1)=isec(ij,1)+1
2476 isec(ij,1+isec(ij,1))=nbeta
2479 isec(ij,1)=isec(ij,1)+1
2480 isec(ij,1+isec(ij,1))=nbeta
2485 if (nbeta.le.9) then
2486 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2487 "DefPropRes 'strand",nstrand,&
2488 "' 'num = ",ii1-1,"..",i1-1,"'"
2490 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2491 "DefPropRes 'strand",nstrand,&
2492 "' 'num = ",ii1-1,"..",i1-1,"'"
2495 if (nbeta.le.9) then
2496 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2497 "DefPropRes 'strand",nstrand,&
2498 "' 'num = ",jj1-1,"..",j1-1,"'"
2500 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2501 "DefPropRes 'strand",nstrand,&
2502 "' 'num = ",jj1-1,"..",j1-1,"'"
2504 write(12,'(a8,4i4)') &
2505 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2511 ! finding antiparallel beta
2512 !d write (iout,*) '--------- looking for antiparallel beta ---------'
2517 if (dcont(i).le.rbeta.and. &
2518 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2519 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2520 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2521 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2522 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2526 !d write (iout,*) i1,j1,dcont(i)
2533 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
2534 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2535 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2536 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2537 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2538 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2539 .and. dcont(j).le.rbeta ) goto 6
2543 !d write (iout,*) i1,j1,dcont(j),not_done
2547 if (i1-ii1.gt.1) then
2548 if(lprint)write (iout,*)'antiparallel beta',&
2549 nbeta,ii1-1,i1,jj1,j1-1
2552 bfrag(1,nbfrag)=max0(ii1-1,1)
2555 bfrag(4,nbfrag)=max0(j1-1,1)
2560 isec(ij,1)=isec(ij,1)+1
2561 isec(ij,1+isec(ij,1))=nbeta
2565 isec(ij,1)=isec(ij,1)+1
2566 isec(ij,1+isec(ij,1))=nbeta
2572 if (nstrand.le.9) then
2573 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2574 "DefPropRes 'strand",nstrand,&
2575 "' 'num = ",ii1-2,"..",i1-1,"'"
2577 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2578 "DefPropRes 'strand",nstrand,&
2579 "' 'num = ",ii1-2,"..",i1-1,"'"
2582 if (nstrand.le.9) then
2583 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2584 "DefPropRes 'strand",nstrand,&
2585 "' 'num = ",j1-2,"..",jj1-1,"'"
2587 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2588 "DefPropRes 'strand",nstrand,&
2589 "' 'num = ",j1-2,"..",jj1-1,"'"
2591 write(12,'(a8,4i4)') &
2592 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2598 if (nstrand.gt.0.and.lprint) then
2599 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2602 write(12,'(a9,i1,$)') " | strand",i
2604 write(12,'(a9,i2,$)') " | strand",i
2607 write(12,'(a1)') "'"
2611 ! finding alpha or 310 helix
2617 if (j1.eq.i1+3.and.dcont(i).le.r310 &
2618 .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2619 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2620 !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2623 if (isec(ii1,1).eq.0) then
2632 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2636 !d write (iout,*) i1,j1,not_done
2639 if (j1-ii1.gt.4) then
2641 !d write (iout,*)'helix',nhelix,ii1,j1
2645 hfrag(2,nhfrag)=max0(j1-1,1)
2651 write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2652 if (nhelix.le.9) then
2653 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2654 "DefPropRes 'helix",nhelix,&
2655 "' 'num = ",ii1-1,"..",j1-2,"'"
2657 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2658 "DefPropRes 'helix",nhelix,&
2659 "' 'num = ",ii1-1,"..",j1-2,"'"
2666 if (nhelix.gt.0.and.lprint) then
2667 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2669 if (nhelix.le.9) then
2670 write(12,'(a8,i1,$)') " | helix",i
2672 write(12,'(a8,i2,$)') " | helix",i
2675 write(12,'(a1)') "'"
2679 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2680 write(12,'(a20)') "XMacStand ribbon.mac"
2684 end subroutine secondary
2685 !-----------------------------------------------------------------------------
2686 subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2688 use geometry, only:dist
2690 ! implicit real*8 (a-h,o-z)
2691 ! include 'DIMENSIONS'
2693 ! include 'COMMON.SBRIDGE'
2694 ! include 'COMMON.FFIELD'
2695 ! include 'COMMON.IOUNITS'
2696 ! include 'COMMON.DISTFIT'
2697 ! include 'COMMON.VAR'
2698 ! include 'COMMON.CHAIN'
2699 ! include 'COMMON.MINIM'
2701 character(len=50) :: linia
2703 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2704 real(kind=8) :: time0,time1
2705 integer :: iff(nres),ieval
2706 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2709 integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2710 real(kind=8) :: wstrain0,etot
2712 maxres22=nres*(nres+1)/2
2714 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2715 call var_to_geom(nvar,var)
2722 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2744 call var_to_geom(nvar,var2)
2747 if ( iff(i).eq.1 ) then
2756 !d call write_pdb(3,'combined structure',0d0)
2757 !d time0=MPI_WTIME()
2760 NY=((NRES-4)*(NRES-5))/2
2761 call distfit(.true.,200)
2763 !d time1=MPI_WTIME()
2764 !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2774 call geom_to_var(nvar,var)
2775 !d time0=MPI_WTIME()
2776 call minimize(etot,var,iretcode,nfun)
2777 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2779 !d time1=MPI_WTIME()
2780 !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2781 !d & nfun/(time1-time0),' SOFT eval/s'
2782 call var_to_geom(nvar,var)
2788 if (iff(1).eq.1) then
2794 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2799 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2805 if (iff(nres).eq.1) then
2811 !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2812 !d & "select",ij(1),"-",ij(2),
2813 !d & ",",ij(3),"-",ij(4)
2814 !d call write_pdb(in_pdb,linia,etot)
2820 !d time0=MPI_WTIME()
2821 call minimize(etot,var,iretcode,nfun)
2822 !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2825 !d time1=MPI_WTIME()
2826 !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2827 !d & nfun/(time1-time0),' eval/s'
2828 !d call var_to_geom(nvar,var)
2830 !d call write_pdb(6,'dist structure',etot)
2839 end subroutine contact_cp2
2840 !-----------------------------------------------------------------------------
2841 subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2843 use geometry, only:dist
2845 ! implicit real*8 (a-h,o-z)
2846 ! include 'DIMENSIONS'
2847 ! include 'COMMON.SBRIDGE'
2848 ! include 'COMMON.FFIELD'
2849 ! include 'COMMON.IOUNITS'
2850 ! include 'COMMON.DISTFIT'
2851 ! include 'COMMON.VAR'
2852 ! include 'COMMON.CHAIN'
2853 ! include 'COMMON.MINIM'
2855 character(len=50) :: linia
2857 real(kind=8) :: energy(0:n_ene)
2858 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2859 real(kind=8) :: time0,time1
2860 integer :: iff(nres),ieval
2861 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2865 integer :: in_pdb,i,j,ind,iwsk
2869 if (ieval.eq.-1) debug=.true.
2873 ! store selected dist. constrains from 1st structure
2876 ! Intercept NaNs in the coordinates
2877 ! write(iout,*) (var(i),i=1,nvar)
2882 if (x_sum.ne.x_sum) then
2883 write(iout,*)" *** contact_cp : Found NaN in coordinates"
2885 print *," *** contact_cp : Found NaN in coordinates"
2891 call var_to_geom(nvar,var)
2898 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2921 ! freeze sec.elements from 2nd structure
2929 call var_to_geom(nvar,var2)
2930 call secondary2(debug)
2932 do i=bfrag(1,j),bfrag(2,j)
2937 if (bfrag(3,j).le.bfrag(4,j)) then
2938 do i=bfrag(3,j),bfrag(4,j)
2944 do i=bfrag(4,j),bfrag(3,j)
2952 do i=hfrag(1,j),hfrag(2,j)
2961 ! copy selected res from 1st to 2nd structure
2965 if ( iff(i).eq.1 ) then
2975 ! prepare description in linia variable
2979 if (iff(1).eq.1) then
2985 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2990 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2996 if (iff(nres).eq.1) then
3001 write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3002 "SELECT",ij(1)-1,"-",ij(2)-1,&
3003 ",",ij(3)-1,"-",ij(4)-1
3009 call contact_cp_min(var,ieval,in_pdb,linia,debug)
3012 end subroutine contact_cp
3013 !-----------------------------------------------------------------------------
3014 subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3018 ! input : theta,phi,alph,omeg,in_pdb,linia,debug
3019 ! output : var,ieval
3021 ! implicit real*8 (a-h,o-z)
3022 ! include 'DIMENSIONS'
3024 ! include 'COMMON.SBRIDGE'
3025 ! include 'COMMON.FFIELD'
3026 ! include 'COMMON.IOUNITS'
3027 ! include 'COMMON.DISTFIT'
3028 ! include 'COMMON.VAR'
3029 ! include 'COMMON.CHAIN'
3030 ! include 'COMMON.MINIM'
3032 character(len=50) :: linia
3034 real(kind=8) :: energy(0:n_ene)
3035 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3036 real(kind=8) :: time0,time1
3037 integer :: ieval,info(3)
3038 logical :: debug,fail,reduce,change !check_var,
3041 integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3043 real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3044 wtor_d01,wstrain0,etot
3046 write(iout,'(a20,i6,a20)') &
3047 '------------------',in_pdb,'-------------------'
3051 call write_pdb(1000+in_pdb,'combined structure',0d0)
3058 ! run optimization of distances
3060 ! uses d0(),w() and mask() for frozen 2D
3062 !test---------------------------------------------
3064 !test NY=((NRES-4)*(NRES-5))/2
3065 !test call distfit(debug,5000)
3097 call geom_to_var(nvar,var)
3098 !de change=reduce(var)
3099 if (check_var(var,info)) then
3100 write(iout,*) 'cp_min error in input'
3101 print *,'cp_min error in input'
3105 !d call etotal(energy(0))
3106 !d call enerprint(energy(0))
3110 !dtest call minimize(etot,var,iretcode,nfun)
3111 !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3114 !d call etotal(energy(0))
3115 !d call enerprint(energy(0))
3134 !test--------------------------------------------------
3140 write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3141 call write_pdb(2000+in_pdb,'distfit structure',0d0)
3149 ! run soft pot. optimization
3151 ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3153 ! mask_phi(),mask_theta(),mask_side(),mask_r
3159 !de change=reduce(var)
3160 !de if (check_var(var,info)) write(iout,*) 'error before soft'
3164 call minimize(etot,var,iretcode,nfun)
3166 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3170 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3171 nfun/(time1-time0),' SOFT eval/s'
3174 call var_to_geom(nvar,var)
3176 call write_pdb(3000+in_pdb,'soft structure',etot)
3180 ! run full UNRES optimization with constrains and frozen 2D
3181 ! the same variables as soft pot. optimizatio
3187 ! check overlaps before calling full UNRES minim
3189 call var_to_geom(nvar,var)
3193 write(iout,*) 'N7 ',energy(0)
3194 if (energy(0).ne.energy(0)) then
3195 write(iout,*) 'N7 error - gives NaN',energy(0)
3199 if (energy(1).eq.1.0d20) then
3200 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3201 call overlap_sc(fail)
3205 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3217 !dte time0=MPI_WTIME()
3218 !de change=reduce(var)
3219 !de if (check_var(var,info)) then
3220 !de write(iout,*) 'error before mask dist'
3221 !de call var_to_geom(nvar,var)
3223 !de call write_pdb(10000+in_pdb,'before mask dist',etot)
3225 !dte call minimize(etot,var,iretcode,nfun)
3226 !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3227 !dte & ' eval ',nfun
3228 !dte ieval=ieval+nfun
3230 !dte time1=MPI_WTIME()
3231 !dte write (iout,'(a,f6.2,f8.2,a)')
3232 !dte & ' Time for mask dist min.',time1-time0,
3233 !dte & nfun/(time1-time0),' eval/s'
3234 !dte call flush(iout)
3237 call var_to_geom(nvar,var)
3239 call write_pdb(4000+in_pdb,'mask dist',etot)
3242 ! switch off freezing of 2D and
3243 ! run full UNRES optimization with constrains
3249 !de change=reduce(var)
3250 !de if (check_var(var,info)) then
3251 !de write(iout,*) 'error before dist'
3252 !de call var_to_geom(nvar,var)
3254 !de call write_pdb(11000+in_pdb,'before dist',etot)
3257 call minimize(etot,var,iretcode,nfun)
3259 !de change=reduce(var)
3260 !de if (check_var(var,info)) then
3261 !de write(iout,*) 'error after dist',ico
3262 !de call var_to_geom(nvar,var)
3264 !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3266 write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3272 write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3273 nfun/(time1-time0),' eval/s'
3275 !de call etotal(energy(0))
3276 !de write(iout,*) 'N7 after dist',energy(0)
3280 call var_to_geom(nvar,var)
3282 call write_pdb(in_pdb,linia,etot)
3294 end subroutine contact_cp_min
3295 !-----------------------------------------------------------------------------
3298 use geometry, only:dist
3300 ! implicit real*8 (a-h,o-z)
3301 ! include 'DIMENSIONS'
3303 ! include 'COMMON.GEO'
3304 ! include 'COMMON.CHAIN'
3305 ! include 'COMMON.IOUNITS'
3306 ! include 'COMMON.VAR'
3307 ! include 'COMMON.CONTROL'
3308 ! include 'COMMON.SBRIDGE'
3309 ! include 'COMMON.FFIELD'
3310 ! include 'COMMON.MINIM'
3311 ! include 'COMMON.INTERACT'
3313 ! include 'COMMON.DISTFIT'
3314 integer :: iff(nres)
3315 real(kind=8) :: time0,time1
3316 real(kind=8) :: energy(0:n_ene),ee
3317 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3319 logical :: debug,ltest,fail
3320 character(len=50) :: linia
3321 integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3323 real(kind=8) :: wstrain0,wang0,etot
3329 !------------------------
3331 ! freeze sec.elements
3341 do i=bfrag(1,j),bfrag(2,j)
3346 if (bfrag(3,j).le.bfrag(4,j)) then
3347 do i=bfrag(3,j),bfrag(4,j)
3353 do i=bfrag(4,j),bfrag(3,j)
3361 do i=hfrag(1,j),hfrag(2,j)
3373 ! store dist. constrains
3377 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3382 dhpb(nhpb)=DIST(i,j)
3390 call write_pdb(100+in_pdb,'input reg. structure',0d0)
3400 ! run soft pot. optimization
3406 call geom_to_var(nvar,var)
3410 call minimize(etot,var,iretcode,nfun)
3412 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3416 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3417 nfun/(time1-time0),' SOFT eval/s'
3419 call var_to_geom(nvar,var)
3421 call write_pdb(300+in_pdb,'soft structure',etot)
3424 ! run full UNRES optimization with constrains and frozen 2D
3425 ! the same variables as soft pot. optimizatio
3434 call minimize(etot,var,iretcode,nfun)
3435 write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3442 write (iout,'(a,f6.2,f8.2,a)') &
3443 ' Time for mask dist min.',time1-time0,&
3444 nfun/(time1-time0),' eval/s'
3446 call var_to_geom(nvar,var)
3448 call write_pdb(400+in_pdb,'mask & dist',etot)
3451 ! switch off constrains and
3452 ! run full UNRES optimization with frozen 2D
3467 call minimize(etot,var,iretcode,nfun)
3468 write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3474 write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,&
3475 nfun/(time1-time0),' eval/s'
3479 call var_to_geom(nvar,var)
3481 call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3488 ! run full UNRES optimization with constrains and NO frozen 2D
3498 wstrain=wstrain0/ico
3503 call minimize(etot,var,iretcode,nfun)
3504 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3505 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3512 write (iout,'(a,f6.2,f8.2,a)') &
3513 ' Time for dist min.',time1-time0,&
3514 nfun/(time1-time0),' eval/s'
3516 call var_to_geom(nvar,var)
3518 call write_pdb(600+in_pdb+ico,'dist cons',etot)
3535 call minimize(etot,var,iretcode,nfun)
3536 write(iout,*)'------------------------------------------------'
3537 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3543 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
3544 nfun/(time1-time0),' eval/s'
3547 call var_to_geom(nvar,var)
3549 call write_pdb(999,'full min',etot)
3553 end subroutine softreg
3554 !-----------------------------------------------------------------------------
3555 subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3557 use geometry, only:dist
3559 ! implicit real*8 (a-h,o-z)
3560 ! include 'DIMENSIONS'
3562 ! include 'COMMON.GEO'
3563 ! include 'COMMON.VAR'
3564 ! include 'COMMON.INTERACT'
3565 ! include 'COMMON.IOUNITS'
3566 ! include 'COMMON.DISTFIT'
3567 ! include 'COMMON.SBRIDGE'
3568 ! include 'COMMON.CONTROL'
3569 ! include 'COMMON.FFIELD'
3570 ! include 'COMMON.MINIM'
3571 ! include 'COMMON.CHAIN'
3572 real(kind=8) :: time0,time1
3573 real(kind=8) :: energy(0:n_ene),ee
3574 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3575 integer :: jdata(5),isec(nres)
3578 integer :: i1,i2,i3,i4,i5,ieval,ij
3579 integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3580 real(kind=8) :: etot,wscloc0,wstrain0
3588 call secondary2(.false.)
3594 do i=bfrag(1,j),bfrag(2,j)
3597 do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3602 do i=hfrag(1,j),hfrag(2,j)
3608 ! cut strands at the ends
3610 if (jdata(2)-jdata(1).gt.3) then
3613 if (jdata(3).lt.jdata(4)) then
3623 !v call etotal(energy(0))
3625 !v write(iout,*) nnt,nct,etot
3626 !v call write_pdb(ij*100,'first structure',etot)
3627 !v write(iout,*) 'N16 test',(jdata(i),i=1,5)
3629 !------------------------
3630 ! generate constrains
3633 if(ishift.eq.0) ishift=-2
3636 do i=jdata(1),jdata(2)
3638 if(jdata(4).gt.jdata(3))then
3639 do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3641 !d print *,i,j,j+ishift
3646 dhpb(nhpb)=DIST(i,j+ishift)
3649 do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3651 !d print *,i,j,j+ishift
3656 dhpb(nhpb)=DIST(i,j+ishift)
3663 if(isec(i).gt.0.or.isec(j).gt.0) then
3669 dhpb(nhpb)=DIST(i,j)
3676 call geom_to_var(nvar,var)
3683 wstrain=wstrain0/ico
3685 !v time0=MPI_WTIME()
3686 call minimize(etot,var,iretcode,nfun)
3687 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3688 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3691 !v time1=MPI_WTIME()
3692 !v write (iout,'(a,f6.2,f8.2,a)')
3693 !v & ' Time for dist min.',time1-time0,
3694 !v & nfun/(time1-time0),' eval/s'
3695 !v call var_to_geom(nvar,var)
3697 !v call write_pdb(ij*100+ico,'dist cons',etot)
3709 call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3712 !v call etotal(energy(0))
3714 !v call write_pdb(ij*100+10,'sc_move',etot)
3716 !d print *,nft_sc,etot
3719 end subroutine beta_slide
3720 !-----------------------------------------------------------------------------
3721 subroutine beta_zip(i1,i2,ieval,ij)
3724 ! implicit real*8 (a-h,o-z)
3725 ! include 'DIMENSIONS'
3727 ! include 'COMMON.GEO'
3728 ! include 'COMMON.VAR'
3729 ! include 'COMMON.INTERACT'
3730 ! include 'COMMON.IOUNITS'
3731 ! include 'COMMON.DISTFIT'
3732 ! include 'COMMON.SBRIDGE'
3733 ! include 'COMMON.CONTROL'
3734 ! include 'COMMON.FFIELD'
3735 ! include 'COMMON.MINIM'
3736 ! include 'COMMON.CHAIN'
3737 real(kind=8) :: time0,time1
3738 real(kind=8) :: energy(0:n_ene),ee
3739 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3740 character(len=10) :: test
3742 integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3743 real(kind=8) :: etot,wstrain0
3745 !v call etotal(energy(0))
3747 !v write(test,'(2i5)') i1,i2
3748 !v call write_pdb(ij*100,test,etot)
3749 !v write(iout,*) 'N17 test',i1,i2,etot,ij
3752 ! generate constrains
3763 call geom_to_var(nvar,var)
3769 wstrain=wstrain0/ico
3770 !v time0=MPI_WTIME()
3771 call minimize(etot,var,iretcode,nfun)
3772 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3773 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3776 !v time1=MPI_WTIME()
3777 !v write (iout,'(a,f6.2,f8.2,a)')
3778 !v & ' Time for dist min.',time1-time0,
3779 !v & nfun/(time1-time0),' eval/s'
3780 ! do not comment the next line
3781 call var_to_geom(nvar,var)
3783 !v call write_pdb(ij*100+ico,'dist cons',etot)
3791 !v call etotal(energy(0))
3793 !v write(iout,*) 'N17 test end',i1,i2,etot,ij
3796 end subroutine beta_zip
3797 !-----------------------------------------------------------------------------
3799 !-----------------------------------------------------------------------------
3800 subroutine thread_seq
3802 use geometry, only:dist
3803 use random, only:iran_num
3804 use control, only:tcpu
3805 use regularize_, only:regularize
3806 use mcm_data, only: nsave_part,nacc_tot
3807 ! Thread the sequence through a database of known structures
3808 ! implicit real*8 (a-h,o-z)
3809 ! include 'DIMENSIONS'
3810 use MPI_data !include 'COMMON.INFO'
3815 ! include 'COMMON.CONTROL'
3816 ! include 'COMMON.CHAIN'
3817 ! include 'COMMON.DBASE'
3818 ! include 'COMMON.INTERACT'
3819 ! include 'COMMON.VAR'
3820 ! include 'COMMON.THREAD'
3821 ! include 'COMMON.FFIELD'
3822 ! include 'COMMON.SBRIDGE'
3823 ! include 'COMMON.HEADER'
3824 ! include 'COMMON.IOUNITS'
3825 ! include 'COMMON.TIME1'
3826 ! include 'COMMON.CONTACTS'
3827 ! include 'COMMON.MCM'
3828 ! include 'COMMON.NAMES'
3830 integer :: ThreadId,ThreadType,Kwita
3832 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3833 real(kind=8) :: przes(3),obr(3,3)
3834 real(kind=8) :: time_for_thread
3835 logical :: found_pattern,non_conv
3836 character(len=32) :: head_pdb
3837 real(kind=8) :: energia(0:n_ene)
3838 integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3840 real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3842 n_ene_comp=nprint_ene
3847 if (me.eq.king) then
3866 ave_time_for_thread=0.0D0
3867 max_time_for_thread=0.0D0
3868 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3869 nthread=nexcl+nthread
3870 do ithread=1,nthread
3871 found_pattern=.false.
3873 do while (.not.found_pattern)
3875 if (itrial.gt.1000) then
3876 write (iout,'(/a/)') 'Too many attempts to find pattern.'
3879 call recv_stop_sig(Kwita)
3880 call send_stop_sig(-3)
3884 ! Find long enough chain in the database
3886 nres_t=nres_base(1,ii)
3887 ! Select the starting position to thread.
3888 print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3889 nres_t,' nres0=',nres0
3890 if (nres_t.ge.nres0) then
3891 ist=iran_num(0,nres_t-nres0)
3893 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3894 if (Kwita.lt.0) then
3895 write (iout,*) 'Stop signal received. Terminating.'
3896 write (*,*) 'Stop signal received. Terminating.'
3898 write (*,*) 'ithread=',ithread,' nthread=',nthread
3901 call pattern_receive
3904 if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3906 found_pattern=.true.
3908 ! If this point is reached, the pattern has not yet been examined.
3910 ! print *,'found_pattern:',found_pattern
3916 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3917 if (Kwita.lt.0) then
3918 write (iout,*) 'Stop signal received. Terminating.'
3920 write (*,*) 'ithread=',ithread,' nthread=',nthread
3926 ipatt(2,ithread)=ist
3928 write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3929 'Processor:',me,' Attempt:',ithread,&
3930 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3931 ' start at res.',ist+1
3932 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3933 ' Attempt:',ithread,&
3934 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3935 ' start at res.',ist+1
3937 write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3938 'Attempt:',ithread,&
3939 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3940 ' start at res.',ist+1
3941 write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3942 'Attempt:',ithread,&
3943 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3944 ' start at res.',ist+1
3947 ! Copy coordinates from the database.
3951 c(j,i)=cart_base(j,i+ist,ii)
3954 !d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
3956 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3958 !d write (iout,'(a,f10.5)')
3959 !d & 'Initial RMS deviation from reference structure:',rms
3960 if (itype(nres,1).eq.ntyp1) then
3962 dcj=c(j,nres-2)-c(j,nres-3)
3963 c(j,nres)=c(j,nres-1)+dcj
3964 c(j,2*nres)=c(j,nres)
3967 if (itype(1,1).eq.ntyp1) then
3974 call int_from_cart(.false.,.false.)
3975 !d print *,'Exit INT_FROM_CART.'
3976 !d print *,'nhpb=',nhpb
3981 ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3983 ! stop 'End generate'
3984 ! Generate SC conformations.
3988 !d print *,'Processor:',me,': exit GEN_SIDE.'
3990 !d print *,'Exit GEN_SIDE.'
3992 ! Calculate initial energy.
3994 call etotal(energia)
3997 ener0(i,ithread)=energia(i)
3999 ener0(n_ene_comp+1,ithread)=energia(0)
4001 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4002 ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
4004 ener0(n_ene_comp+2,ithread)=rms
4005 ener0(n_ene_comp+4,ithread)=frac
4006 ener0(n_ene_comp+5,ithread)=frac_nn
4008 ener0(n_ene_comp+3,ithread)=0.0d0
4011 print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
4013 print*,'ithread=',ithread,' Start REGULARIZE.'
4016 call regularize(nct-nnt+1,etot,rms,&
4017 cart_base(1,ist+nnt,ipattern),iretcode)
4019 time_for_thread=curr_tim1-curr_tim
4020 ave_time_for_thread= &
4021 ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4022 if (time_for_thread.gt.max_time_for_thread) &
4023 max_time_for_thread=time_for_thread
4025 print *,'Processor',me,': Exit REGULARIZE.'
4026 if (WhatsUp.eq.2) then
4028 'Sufficient number of confs. collected. Terminating.'
4031 else if (WhatsUp.eq.-1) then
4033 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4034 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4035 call send_stop_sig(-2)
4037 else if (WhatsUp.eq.-2) then
4039 write (iout,*) 'Timeup signal received. Terminating.'
4041 else if (WhatsUp.eq.-3) then
4043 write (iout,*) 'Error stop signal received. Terminating.'
4047 print *,'Exit REGULARIZE.'
4048 if (iretcode.eq.11) then
4049 write (iout,'(/a/)') &
4050 '******* Allocated time exceeded in SUMSL. The program will stop.'
4055 head_pdb=titel(:24)//':'//str_nam(ipattern)
4056 if (outpdb) call pdbout(etot,head_pdb,ipdb)
4057 if (outmol2) call mol2out(etot,head_pdb)
4059 call briefout(ithread,etot)
4061 link_end=min0(link_end,nss)
4062 write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4064 call etotal(energia)
4065 ! call enerprint(energia(0))
4068 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4069 !d write (iout,'(a,f10.5)')
4070 !d & 'RMS deviation from reference structure:',dsqrt(rms)
4072 ener(i,ithread)=energia(i)
4074 ener(n_ene_comp+1,ithread)=energia(0)
4075 ener(n_ene_comp+3,ithread)=rms
4077 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4078 ener(n_ene_comp+2,ithread)=rms
4079 ener(n_ene_comp+4,ithread)=frac
4080 ener(n_ene_comp+5,ithread)=frac_nn
4082 call write_stat_thread(ithread,ipattern,ist)
4083 ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4084 ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4085 ! & (ener(k,ithread),k=12,14)
4087 if (me.eq.king) then
4089 call pattern_receive
4090 call receive_MCM_info
4091 if (nacc_tot.ge.nthread) then
4093 'Sufficient number of conformations collected nacc_tot=',&
4094 nacc_tot,'. Stopping other processors and terminating.'
4096 'Sufficient number of conformations collected nacc_tot=',&
4097 nacc_tot,'. Stopping other processors and terminating.'
4098 call recv_stop_sig(Kwita)
4099 if (Kwita.eq.0) call send_stop_sig(-1)
4104 call send_MCM_info(2)
4107 if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4108 write (iout,'(/2a)') &
4109 '********** There would be not enough time for another thread. ',&
4110 'The program will stop.'
4112 '********** There would be not enough time for another thread. ',&
4113 'The program will stop.'
4114 write (iout,'(a,1pe14.4/)') &
4115 'Elapsed time for last threading step: ',time_for_thread
4118 call recv_stop_sig(Kwita)
4119 call send_stop_sig(-2)
4124 write (iout,'(a,1pe14.4)') &
4125 'Elapsed time for this threading step: ',time_for_thread
4128 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4129 if (Kwita.lt.0) then
4130 write (iout,*) 'Stop signal received. Terminating.'
4131 write (*,*) 'Stop signal received. Terminating.'
4133 write (*,*) 'nthread=',nthread,' ithread=',ithread
4139 call send_stop_sig(-1)
4143 ! Any messages left for me?
4144 call pattern_receive
4145 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4147 call write_thread_summary
4149 if (king.eq.king) then
4151 do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4154 call recv_stop_sig(Kwita)
4155 call receive_MCM_info
4158 call receive_thread_results(iproc)
4160 call write_thread_summary
4162 call send_thread_results
4166 end subroutine thread_seq
4167 !-----------------------------------------------------------------------------
4170 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4172 use random, only:iran_num
4173 ! implicit real*8 (a-h,o-z)
4174 ! include 'DIMENSIONS'
4175 ! include 'COMMON.CHAIN'
4176 ! include 'COMMON.DBASE'
4177 ! include 'COMMON.INTERACT'
4178 ! include 'COMMON.VAR'
4179 ! include 'COMMON.THREAD'
4180 ! include 'COMMON.FFIELD'
4181 ! include 'COMMON.SBRIDGE'
4182 ! include 'COMMON.HEADER'
4183 ! include 'COMMON.GEO'
4184 ! include 'COMMON.IOUNITS'
4185 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
4186 !el integer :: icall
4187 !el common /srutu/ icall
4188 real(kind=8) :: energia(0:n_ene)
4189 logical :: glycine,fail
4190 integer :: i,maxsample,link_end0,ind_sc,isample
4191 real(kind=8) :: alph0,omeg0,e1,e0
4195 link_end=min0(link_end,nss)
4197 if (itype(i,1).ne.10) then
4198 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
4199 call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1)
4203 call etotal(energia)
4205 do isample=1,maxsample
4206 ! Choose a non-glycine side chain.
4209 ind_sc=iran_num(nnt,nct)
4210 glycine=(itype(ind_sc,1).eq.10)
4214 call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
4215 omeg(ind_sc),fail,1)
4217 call etotal(energia)
4218 !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4219 !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4220 !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4231 end subroutine sc_conf
4232 !-----------------------------------------------------------------------------
4234 !-----------------------------------------------------------------------------
4235 logical function check_var(var,info)
4239 ! implicit real*8 (a-h,o-z)
4240 ! include 'DIMENSIONS'
4241 ! include 'COMMON.VAR'
4242 ! include 'COMMON.IOUNITS'
4243 ! include 'COMMON.GEO'
4244 ! include 'COMMON.SETUP'
4245 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4246 integer,dimension(3) :: info
4250 do i=nphi+ntheta+1,nphi+ntheta+nside
4251 ! Check the side chain "valence" angles alpha
4252 if (var(i).lt.1.0d-7) then
4253 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4254 write (iout,*) 'Processor',me,'received bad variables!!!!'
4255 write (iout,*) 'Variables'
4256 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4257 write (iout,*) 'Continuing calculations at this point',&
4258 ' could destroy the results obtained so far... ABORTING!!!!!!'
4259 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4260 'valence angle alpha',i-nphi-ntheta,var(i),&
4261 'n it',info(1),info(2),'mv ',info(3)
4262 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4263 write (*,*) 'Processor',me,'received bad variables!!!!'
4264 write (*,*) 'Variables'
4265 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4266 write (*,*) 'Continuing calculations at this point',&
4267 ' could destroy the results obtained so far... ABORTING!!!!!!'
4268 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4269 'valence angle alpha',i-nphi-ntheta,var(i),&
4270 'n it',info(1),info(2),'mv ',info(3)
4275 ! Check the backbone "valence" angles theta
4276 do i=nphi+1,nphi+ntheta
4277 if (var(i).lt.1.0d-7) then
4278 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4279 write (iout,*) 'Processor',me,'received bad variables!!!!'
4280 write (iout,*) 'Variables'
4281 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4282 write (iout,*) 'Continuing calculations at this point',&
4283 ' could destroy the results obtained so far... ABORTING!!!!!!'
4284 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4285 'valence angle theta',i-nphi,var(i),&
4286 'n it',info(1),info(2),'mv ',info(3)
4287 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4288 write (*,*) 'Processor',me,'received bad variables!!!!'
4289 write (*,*) 'Variables'
4290 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4291 write (*,*) 'Continuing calculations at this point',&
4292 ' could destroy the results obtained so far... ABORTING!!!!!!'
4293 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4294 'valence angle theta',i-nphi,var(i),&
4295 'n it',info(1),info(2),'mv ',info(3)
4301 end function check_var
4302 !-----------------------------------------------------------------------------
4304 !-----------------------------------------------------------------------------
4305 subroutine distfit(debug,maxit)
4307 use geometry_data, only: phi
4310 ! implicit real*8 (a-h,o-z)
4311 ! include 'DIMENSIONS'
4312 ! include 'COMMON.CHAIN'
4313 ! include 'COMMON.VAR'
4314 ! include 'COMMON.IOUNITS'
4315 ! include 'COMMON.DISTFIT'
4316 integer :: i,maxit,MAXMAR,IT,IMAR
4317 real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres)
4318 logical :: debug,sing
4319 real(kind=8) :: TOL,RL,F0,AIN,F1
4321 !input------------------------------------
4323 ! NY=((NRES-4)*(NRES-5))/2
4324 !input------------------------------------
4330 CALL TRANSFER(NRES,phi,phiold)
4334 !d WRITE (IOUT,*) 'DISTFIT: F0=',F0
4350 CALL TRANSFER(NX,XX,X)
4351 CALL BANACH(NX,NRES,H,X,sing)
4356 IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4358 WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4359 WRITE (IOUT,*) 'IT=',it,'F=',F0
4364 phi(I)=phiold(I)+mask(i)*X(I-3)
4369 !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4371 CALL TRANSFER(NRES,phi,phiold)
4374 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4376 WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4377 WRITE (IOUT,*) 'IT=',it,'F=',F1
4383 WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4384 WRITE (IOUT,*) 'IT=',it,'F=',F0
4385 CALL TRANSFER(NRES,phiold,phi)
4388 !d write (iout,*) "it",it," imar",imar," f0",f0
4390 WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4392 end subroutine distfit
4393 !-----------------------------------------------------------------------------
4394 real(kind=8) function RDIF()
4397 use geometry, only: dist
4398 ! implicit real*8 (a-h,o-z)
4399 ! include 'DIMENSIONS'
4400 ! include 'COMMON.CHAIN'
4401 ! include 'COMMON.DISTFIT'
4403 real(kind=8) :: suma,DIJ
4412 if (w(ind).ne.0.0) then
4414 suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4416 ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4424 !-----------------------------------------------------------------------------
4429 use geometry, only:dist
4430 ! implicit real*8 (a-h,o-z)
4431 ! include 'DIMENSIONS'
4432 ! include 'COMMON.CHAIN'
4433 ! include 'COMMON.DISTFIT'
4434 ! include 'COMMON.GEO'
4435 integer :: i,j,k,l,I1,I2,IND
4436 real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4449 R13(K)=C(K,J)-C(K,I1)
4453 R24(L)=C(L,K)-C(L,I2)
4455 IND=((J-1)*(2*NRES-J-6))/2+K-3
4456 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4457 PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4458 PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4459 DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4464 end subroutine RDERIV
4465 !-----------------------------------------------------------------------------
4469 ! implicit real*8 (a-h,o-z)
4470 ! include 'DIMENSIONS'
4471 ! include 'COMMON.CHAIN'
4472 ! include 'COMMON.DISTFIT'
4474 real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4482 XI=XI+BKIWK*(D0(K)-DDD(K))
4490 HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4497 end subroutine HEVAL
4498 !-----------------------------------------------------------------------------
4499 subroutine VEC(I,J,U)
4501 use geometry_data, only: C
4502 ! Find the unit vector from atom (I) to atom (J). Store in U.
4504 ! implicit real*8 (a-h,o-z)
4505 ! include 'DIMENSIONS'
4506 ! include 'COMMON.CHAIN'
4508 real(kind=8),DIMENSION(3) :: U
4509 real(kind=8) :: ANORM,UK
4523 !-----------------------------------------------------------------------------
4524 subroutine TRANSFER(N,X1,X2)
4526 ! implicit real*8 (a-h,o-z)
4527 ! include 'DIMENSIONS'
4529 real(kind=8),DIMENSION(N) :: X1,X2
4533 end subroutine TRANSFER
4534 !-----------------------------------------------------------------------------
4535 !-----------------------------------------------------------------------------
4536 subroutine alloc_compare_arrays
4538 maxres22=nres*(nres+1)/2
4540 ! common /struct/ in io_common: read_threadbase
4541 ! allocate(cart_base !(3,maxres_base,maxseq)
4542 ! allocate(nres_base !(3,maxseq)
4543 ! allocate(str_nam !(maxseq)
4545 ! COMMON /c_frag/ in io_conf: readpdb
4546 if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4547 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4549 allocate(w(maxres22),d0(maxres22)) !(maxres22)
4551 !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4552 allocate(DDD(maxres22)) !(maxres22)
4553 allocate(H(nres,nres)) !(MAXRES,MAXRES)
4554 allocate(XX(nres)) !(MAXRES)
4556 allocate(mask(nres)) !(maxres)
4559 allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4561 allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4564 end subroutine alloc_compare_arrays
4565 !-----------------------------------------------------------------------------
4567 !-----------------------------------------------------------------------------