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
49 iti=iabs(itype(i,molnum(i)))
51 itj=iabs(itype(j,molnum(i)))
53 ! rcomp=sigmaii(iti,itj)+1.0D0
54 rcomp=facont*sigmaii(iti,itj)
56 ! rcomp=sigma(iti,itj)+1.0D0
57 rcomp=facont*sigma(iti,itj)
60 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
61 if (dist(nres+i,nres+j).lt.rcomp) then
69 write (iout,'(a)') 'Contact map:'
75 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
76 i,restyp(it1,1),i1,restyp(it2,1),i2
81 co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
83 co = co / (nres*ncont)
85 end subroutine contact
86 !-----------------------------------------------------------------------------
87 real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref)
89 ! implicit real*8 (a-h,o-z)
90 ! include 'DIMENSIONS'
91 ! include 'COMMON.IOUNITS'
92 integer :: ncont,ncont_ref
93 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
97 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
98 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
99 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
100 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
101 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
104 if (icont(1,i).eq.icont_ref(1,j) .and. &
105 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
108 ! print *,' nmatch=',nmatch
109 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
110 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
112 end function contact_fract
113 !-----------------------------------------------------------------------------
114 real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
116 ! implicit real*8 (a-h,o-z)
117 ! include 'DIMENSIONS'
118 ! include 'COMMON.IOUNITS'
119 integer :: ncont,ncont_ref
120 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
122 integer :: i,j,nmatch
124 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
125 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
126 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
127 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
128 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
131 if (icont(1,i).eq.icont_ref(1,j) .and. &
132 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
135 ! print *,' nmatch=',nmatch
136 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
137 contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
139 end function contact_fract_nn
140 !-----------------------------------------------------------------------------
141 subroutine hairpin(lprint,nharp,iharp)
143 use geometry, only:dist
144 ! implicit real*8 (a-h,o-z)
145 ! include 'DIMENSIONS'
146 ! include 'COMMON.IOUNITS'
147 ! include 'COMMON.CHAIN'
148 ! include 'COMMON.INTERACT'
149 ! include 'COMMON.FFIELD'
150 ! include 'COMMON.NAMES'
152 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
154 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
155 logical :: lprint,not_done
156 real(kind=8) :: rcomp=6.0d0
158 integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1
159 ! allocate(icont(2,12*nres))
163 ! print *,'nnt=',nnt,' nct=',nct
166 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
170 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
172 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
180 write (iout,'(a)') 'PP contact map:'
186 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
187 i,restyp(it1,1),i1,restyp(it2,1),i2
195 if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
196 ! write (iout,*) "found turn at ",i1,j1
204 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
208 ! write (iout,*) i1,j1,not_done
218 ! write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
223 ! write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
226 write (iout,*) "Hairpins:"
233 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
234 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
236 ! write (iout,'(a,i3,$)') restyp(itype(k,1)),k
241 end subroutine hairpin
242 !-----------------------------------------------------------------------------
244 !-----------------------------------------------------------------------------
245 subroutine elecont(lprint,ncont,icont)
247 ! implicit real*8 (a-h,o-z)
248 ! include 'DIMENSIONS'
249 ! include 'COMMON.IOUNITS'
250 ! include 'COMMON.CHAIN'
251 ! include 'COMMON.INTERACT'
252 ! include 'COMMON.LOCAL'
253 ! include 'COMMON.FFIELD'
254 ! include 'COMMON.NAMES'
256 real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_
257 real(kind=8) :: ael6_i,ael3_i
258 real(kind=8),dimension(2,2) :: app_,bpp_,rpp_
260 integer,dimension(2,12*nres) :: icont !(2,12*nres)(2,maxcont) (maxcont=12*maxres)
261 real(kind=8),dimension(12*nres) :: econt !(maxcont)
263 integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2
264 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw
265 real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb
266 real(kind=8) :: xmedi,ymedi,zmedi
267 real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij
268 real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,&
269 evdwij,el1,el2,eesij,ene
271 ! Load the constants of peptide bond - peptide bond interactions.
272 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
273 ! proline) - determined by averaging ECEPP energy.
277 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
278 data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
279 data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
280 data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
282 !el allocate(econt(12*nres)) !(maxcont)
285 elecutoff_14 = -0.5d0
286 if (lprint) write (iout,'(a)') &
287 "Constants of electrostatic interaction energy expression."
291 app_(i,j)=epp(i,j)*rri*rri
292 bpp_(i,j)=-2.0*epp(i,j)*rri
293 ael6_(i,j)=elpp_6(i,j)*4.2**6
294 ael3_(i,j)=elpp_3(i,j)*4.2**3
296 write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),&
303 print *, "nntt,nct",nnt,nct-2
305 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
316 if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
319 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
320 if (iteli.eq.2 .and. itelj.eq.2) goto 4
321 aaa=app_(iteli,itelj)
322 bbb=bpp_(iteli,itelj)
323 ael6_i=ael6_(iteli,itelj)
324 ael3_i=ael3_(iteli,itelj)
328 xj=c(1,j)+0.5*dxj-xmedi
329 yj=c(2,j)+0.5*dyj-ymedi
330 zj=c(3,j)+0.5*dzj-zmedi
331 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
336 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
337 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
338 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
339 fac=cosa-3.0*cosb*cosg
345 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
348 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
349 j.eq.i+2 .and. eesij.le.elecutoff_14) then
360 write (iout,*) 'Total average electrostatic energy: ',ees
361 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
363 write (iout,*) 'Electrostatic contacts before pruning: '
369 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
370 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
373 ! For given residues keep only the contacts with the greatest energy.
375 do while (i.lt.ncont)
381 do while (j.lt.ncont)
383 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
384 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
385 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
386 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
387 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
388 if (ic1.eq.icont(1,j)) then
390 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
391 .and. iabs(icont(1,k)-ic1).le.2 .and. &
392 econt(k).lt.econt(j) ) goto 21
394 else if (ic2.eq.icont(2,j) ) then
396 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
397 .and. iabs(icont(2,k)-ic2).le.2 .and. &
398 econt(k).lt.econt(j) ) goto 21
403 icont(1,k-1)=icont(1,k)
404 icont(2,k-1)=icont(2,k)
409 ! write (iout,*) "ncont",ncont
411 ! write (iout,*) icont(1,k),icont(2,k)
414 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
416 if (ic1.eq.icont(1,j)) then
418 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
419 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
420 econt(k).lt.econt(i) ) goto 21
422 else if (ic2.eq.icont(2,j) ) then
424 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
425 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
426 econt(k).lt.econt(i) ) goto 21
431 icont(1,k-1)=icont(1,k)
432 icont(2,k-1)=icont(2,k)
436 ! write (iout,*) "ncont",ncont
438 ! write (iout,*) icont(1,k),icont(2,k)
449 write (iout,*) 'Electrostatic contacts after pruning: '
455 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
456 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
460 end subroutine elecont
461 !-----------------------------------------------------------------------------
462 subroutine secondary2(lprint)
464 ! implicit real*8 (a-h,o-z)
465 ! include 'DIMENSIONS'
466 ! include 'COMMON.CHAIN'
467 ! include 'COMMON.IOUNITS'
468 ! include 'COMMON.DISTFIT'
469 ! include 'COMMON.VAR'
470 ! include 'COMMON.GEO'
471 ! include 'COMMON.CONTROL'
472 integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
474 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
475 integer,dimension(nres,0:4) :: isec !(maxres,4)
476 integer,dimension(nres) :: nsec !(maxres)
477 logical :: lprint,not_done !,freeres
478 real(kind=8) :: p1,p2
481 !el allocate(icont(2,12*nres),isec(nres,4),nsec(nres))
483 if(.not.dccart) call chainbuild_cart
484 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
485 !d call write_pdb(99,'sec structure',0d0)
495 call elecont(lprint,ncont,icont)
496 print *,"after elecont"
497 if (nres_molec(1).eq.0) return
499 ! finding parallel beta
500 !d write (iout,*) '------- looking for parallel beta -----------'
506 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
509 !d write (iout,*) i1,j1
515 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
516 freeres(i1,j1,nsec,isec)) goto 5
520 !d write (iout,*) i1,j1,not_done
524 if (i1-ii1.gt.1) then
528 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
532 bfrag(1,nbfrag)=ii1+1
534 bfrag(3,nbfrag)=jj1+1
535 bfrag(4,nbfrag)=min0(j1+1,nres)
539 isec(ij,nsec(ij))=nbeta
543 isec(ij,nsec(ij))=nbeta
549 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
550 "DefPropRes 'strand",nstrand,&
551 "' 'num = ",ii1-1,"..",i1-1,"'"
553 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
554 "DefPropRes 'strand",nstrand,&
555 "' 'num = ",ii1-1,"..",i1-1,"'"
559 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
560 "DefPropRes 'strand",nstrand,&
561 "' 'num = ",jj1-1,"..",j1-1,"'"
563 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
564 "DefPropRes 'strand",nstrand,&
565 "' 'num = ",jj1-1,"..",j1-1,"'"
567 write(12,'(a8,4i4)') &
568 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
574 ! finding alpha or 310 helix
581 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
584 if (j1.eq.i1+3 .and. &
585 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
586 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
587 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
588 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
591 if (nsec(ii1).eq.0) then
600 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
606 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
611 if (j1-ii1.gt.5) then
623 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
624 if (nhelix.le.9) then
625 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
626 "DefPropRes 'helix",nhelix,&
627 "' 'num = ",ii1-1,"..",j1-2,"'"
629 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
630 "DefPropRes 'helix",nhelix,&
631 "' 'num = ",ii1-1,"..",j1-2,"'"
637 if (nhelix.gt.0.and.lprint) then
638 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
640 if (nhelix.le.9) then
641 write(12,'(a8,i1,$)') " | helix",i
643 write(12,'(a8,i2,$)') " | helix",i
650 ! finding antiparallel beta
651 !d write (iout,*) '--------- looking for antiparallel beta ---------'
656 if (freeres(i1,j1,nsec,isec)) then
659 !d write (iout,*) i1,j1
666 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
667 freeres(i1,j1,nsec,isec)) goto 6
671 !d write (iout,*) i1,j1,not_done
675 if (i1-ii1.gt.1) then
679 bfrag(2,nbfrag)=min0(i1+1,nres)
680 bfrag(3,nbfrag)=min0(jj1+1,nres)
687 if (nsec(ij).le.2) then
688 isec(ij,nsec(ij))=nbeta
694 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
695 isec(ij,nsec(ij))=nbeta
701 write (iout,'(a,i3,4i4)')'antiparallel beta',&
702 nbeta,ii1-1,i1,jj1,j1-1
704 if (nstrand.le.9) then
705 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
706 "DefPropRes 'strand",nstrand,&
707 "' 'num = ",ii1-2,"..",i1-1,"'"
709 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
710 "DefPropRes 'strand",nstrand,&
711 "' 'num = ",ii1-2,"..",i1-1,"'"
714 if (nstrand.le.9) then
715 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
716 "DefPropRes 'strand",nstrand,&
717 "' 'num = ",j1-2,"..",jj1-1,"'"
719 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
720 "DefPropRes 'strand",nstrand,&
721 "' 'num = ",j1-2,"..",jj1-1,"'"
723 write(12,'(a8,4i4)') &
724 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
730 if (nstrand.gt.0.and.lprint) then
731 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
734 write(12,'(a9,i1,$)') " | strand",i
736 write(12,'(a9,i2,$)') " | strand",i
745 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
746 write(12,'(a20)') "XMacStand ribbon.mac"
748 if (nres_molec(1).eq.0) return
749 write(iout,*) 'UNRES seq:'
751 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
755 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
760 end subroutine secondary2
762 !-----------------------------------------------------------------------------
763 logical function freeres(i,j,nsec,isec)
765 ! implicit real*8 (a-h,o-z)
766 ! include 'DIMENSIONS'
767 integer,dimension(nres,4) :: isec !(maxres,4)
768 integer,dimension(nres) :: nsec !(maxres)
775 if (nsec(i).lt.0.or.nsec(j).lt.0) return
777 if (nsec(i).gt.1.or.nsec(j).gt.1) return
780 if (isec(i,k).eq.isec(j,l)) return
786 !-----------------------------------------------------------------------------
788 !-----------------------------------------------------------------------------
789 logical function seq_comp(itypea,itypeb,length)
792 integer :: length,itypea(length),itypeb(length)
795 if (itypea(i).ne.itypeb(i)) then
802 end function seq_comp
804 !-----------------------------------------------------------------------------
806 !-----------------------------------------------------------------------------
807 subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
809 ! implicit real*8 (a-h,o-z)
810 ! include 'DIMENSIONS'
811 ! include 'COMMON.CHAIN'
812 ! include 'COMMON.CONTACTS'
813 ! include 'COMMON.IOUNITS'
814 real(kind=8) :: przes(3),obr(3,3)
815 logical :: non_conv,lprn
816 real(kind=8) :: rms,frac,frac_nn,co
817 ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
821 ! print *,"before contact"
822 !elte(iout,*) "rms_nacc before contact"
823 call contact(.false.,ncont,icont,co)
824 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
825 frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
826 if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
827 'RMS deviation from the reference structure:',rms,&
828 ' % of native contacts:',frac*100,&
829 ' % of nonnative contacts:',frac_nn*100,&
833 end subroutine rms_nac_nnc
834 !-----------------------------------------------------------------------------
835 subroutine rmsd(drms)
837 use regularize_, only:fitsq
838 ! implicit real*8 (a-h,o-z)
839 ! include 'DIMENSIONS'
843 ! include 'COMMON.CHAIN'
844 ! include 'COMMON.IOUNITS'
845 ! include 'COMMON.INTERACT'
846 ! include 'COMMON.CONTROL'
848 real(kind=8) :: przes(3),obrot(3,3)
849 real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
852 real(kind=8) :: drms,rminroz,roznica
853 integer :: i,j,iatom,kkk,iti,k
855 !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
863 ! print *,"nz_start",nz_start," nz_end",nz_end
864 ! if (symetr.le.1) then
866 ! do i=nz_start,nz_end
870 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
871 ! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
873 ! if (iz_sc.eq.1.and.iti.ne.10) then
876 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
877 ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
884 print *,nz_start,nz_end,nstart_seq-nstart_sup
889 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
890 crefcopy(k,iatom)=cref(k,i,kkk)
892 if (iz_sc.eq.1.and.iti.ne.10) then
895 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
896 crefcopy(k,iatom)=cref(k,nres+i,kkk)
905 ! write (iout,*) 'Ccopy and CREFcopy'
906 ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
907 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
908 ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
909 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
911 ! ----- end diagnostics
913 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
914 przes,obrot,non_conv)
916 print *,'Problems in FITSQ!!! rmsd'
917 write (iout,*) 'Problems in FITSQ!!! rmsd'
918 print *,'Ccopy and CREFcopy'
919 write (iout,*) 'Ccopy and CREFcopy'
920 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
921 (crefcopy(j,k),j=1,3),k=1,iatom)
922 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
923 (crefcopy(j,k),j=1,3),k=1,iatom)
925 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
931 ! write (iout,*) "roznica", roznica,kkk
932 if (roznica.le.rminroz) rminroz=roznica
934 drms=dsqrt(dabs(rminroz))
936 ! write (iout,*) "nperm,symetr", nperm,symetr
937 ! ---- end diagnostics
940 !-----------------------------------------------------------------------------
941 subroutine rmsd_csa(drms)
943 use regularize_, only:fitsq
944 ! implicit real*8 (a-h,o-z)
945 ! include 'DIMENSIONS'
949 ! include 'COMMON.CHAIN'
950 ! include 'COMMON.IOUNITS'
951 ! include 'COMMON.INTERACT'
953 real(kind=8) :: przes(3),obrot(3,3)
954 real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
955 integer :: kkk,iatom,ierror,ierrcode
959 real(kind=8) :: drms,roznica
961 allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
969 ccopy(k,iatom)=c(k,i)
970 crefcopy(k,iatom)=crefjlee(k,i)
972 if (iz_sc.eq.1.and.iti.ne.10) then
975 ccopy(k,iatom)=c(k,nres+i)
976 crefcopy(k,iatom)=crefjlee(k,nres+i)
981 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
982 przes,obrot,non_conv)
984 print *,'Problems in FITSQ!!! rmsd_csa'
985 write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
986 print *,'Ccopy and CREFcopy'
987 write (iout,*) 'Ccopy and CREFcopy'
988 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
989 (crefcopy(j,k),j=1,3),k=1,iatom)
990 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
991 (crefcopy(j,k),j=1,3),k=1,iatom)
993 call mpi_abort(mpi_comm_world,ierror,ierrcode)
998 drms=dsqrt(dabs(roznica))
1000 end subroutine rmsd_csa
1001 !-----------------------------------------------------------------------------
1003 !-----------------------------------------------------------------------------
1007 use geometry, only:pinorm
1008 use random, only:ran_number,iran_num
1009 ! implicit real*8 (a-h,o-z)
1010 ! include 'DIMENSIONS'
1012 ! include 'COMMON.GEO'
1013 ! include 'COMMON.VAR'
1014 ! include 'COMMON.INTERACT'
1015 ! include 'COMMON.IOUNITS'
1016 ! include 'COMMON.DISTFIT'
1017 ! include 'COMMON.SBRIDGE'
1018 ! include 'COMMON.CONTROL'
1019 ! include 'COMMON.FFIELD'
1020 ! include 'COMMON.MINIM'
1021 ! include 'COMMON.CHAIN'
1022 real(kind=8) :: time0,time1
1023 real(kind=8) :: energy(0:n_ene),ee
1024 real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres)
1025 integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1026 logical :: debug,accepted
1027 real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1030 !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1032 call geom_to_var(nvar,var1)
1037 write(iout,*) 'etot=',0,etot,rms
1038 call secondary2(.false.)
1040 call write_pdb(0,'first structure',etot)
1049 betbol=1.0D0/(1.9858D-3*temp)
1051 d=ran_number(-pi,pi)
1052 ! phi(jr)=pinorm(phi(jr)+d)
1057 write(iout,*) 'etot=',1,etot0,rms
1058 call write_pdb(1,'perturb structure',etot0)
1062 d=ran_number(-da,da)
1064 phi(jr)=pinorm(phi(jr)+d)
1069 if (etot.lt.etot0) then
1073 xxr=ran_number(0.0D0,1.0D0)
1074 xxh=betbol*(etot-etot0)
1075 if (xxh.lt.50.0D0) then
1077 if (xxh.gt.xxr) accepted=.true.
1081 ! print *,etot0,etot,accepted
1085 write(iout,*) 'etot=',i,etot,rms
1086 call write_pdb(i,'MC structure',etot)
1088 ! call geom_to_var(nvar,var1)
1089 call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1090 call geom_to_var(nvar,var)
1091 call minimize(etot,var,iretcode,nfun)
1092 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1093 call var_to_geom(nvar,var)
1096 write(iout,*) 'etot mcm=',i,etot,rms
1097 call write_pdb(i+1,'MCM structure',etot)
1098 call var_to_geom(nvar,var1)
1106 ! call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1107 ! call geom_to_var(nvar,var)
1110 ! call write_pdb(998 ,'sc min',etot)
1112 ! call minimize(etot,var,iretcode,nfun)
1113 ! write(iout,*)'------------------------------------------------'
1114 ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1116 ! call var_to_geom(nvar,var)
1118 ! call write_pdb(999,'full min',etot)
1122 !-----------------------------------------------------------------------------
1127 ! implicit real*8 (a-h,o-z)
1128 ! include 'DIMENSIONS'
1130 ! include 'COMMON.GEO'
1131 ! include 'COMMON.VAR'
1132 ! include 'COMMON.INTERACT'
1133 ! include 'COMMON.IOUNITS'
1134 ! include 'COMMON.DISTFIT'
1135 ! include 'COMMON.SBRIDGE'
1136 ! include 'COMMON.CONTROL'
1137 ! include 'COMMON.FFIELD'
1138 ! include 'COMMON.MINIM'
1139 ! include 'COMMON.CHAIN'
1140 real(kind=8) :: time0,time1
1141 real(kind=8) :: energy(0:n_ene),ee
1142 real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1146 integer :: i,ij,ieval,iretcode,nfun
1147 real(kind=8) :: etot
1149 allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1151 call geom_to_var(nvar,var1)
1155 write(iout,*) nnt,nct,etot
1156 call write_pdb(1,'first structure',etot)
1157 call secondary2(.true.)
1166 call var_to_geom(nvar,var1)
1167 write(iout,*) 'N16 test',(jdata(i),i=1,5)
1168 call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1170 call geom_to_var(nvar,var)
1176 call minimize(etot,var,iretcode,nfun)
1177 write(iout,*)'------------------------------------------------'
1178 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1184 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1185 nfun/(time1-time0),' eval/s'
1187 call var_to_geom(nvar,var)
1189 call write_pdb(ij*100+99,'full min',etot)
1196 end subroutine test_n16
1198 !-----------------------------------------------------------------------------
1199 subroutine test_local
1201 ! implicit real*8 (a-h,o-z)
1202 ! include 'DIMENSIONS'
1203 ! include 'COMMON.GEO'
1204 ! include 'COMMON.VAR'
1205 ! include 'COMMON.INTERACT'
1206 ! include 'COMMON.IOUNITS'
1207 real(kind=8) :: time0,time1
1208 real(kind=8) :: energy(0:n_ene),ee
1209 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1211 real(kind=8) :: etot
1213 ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres)
1215 ! call geom_to_var(nvar,varia)
1216 call write_pdb(1,'first structure',0d0)
1220 write(iout,*) nnt,nct,etot
1222 write(iout,*) 'calling sc_move'
1223 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1224 write(iout,*) nft_sc,etot
1225 call write_pdb(2,'second structure',etot)
1227 write(iout,*) 'calling local_move'
1228 call local_move_init(.false.)
1229 call local_move(24,29,20d0,50d0)
1231 call write_pdb(3,'third structure',etot)
1233 write(iout,*) 'calling sc_move'
1234 call sc_move(24,29,5,10d0,nft_sc,etot)
1235 write(iout,*) nft_sc,etot
1236 call write_pdb(2,'last structure',etot)
1239 end subroutine test_local
1240 !-----------------------------------------------------------------------------
1243 ! implicit real*8 (a-h,o-z)
1244 ! include 'DIMENSIONS'
1245 ! include 'COMMON.GEO'
1246 ! include 'COMMON.VAR'
1247 ! include 'COMMON.INTERACT'
1248 ! include 'COMMON.IOUNITS'
1249 real(kind=8) :: time0,time1,etot
1250 real(kind=8) :: energy(0:n_ene),ee
1251 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1255 ! call geom_to_var(nvar,varia)
1256 call write_pdb(1,'first structure',0d0)
1260 write(iout,*) nnt,nct,etot
1262 write(iout,*) 'calling sc_move'
1264 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1265 write(iout,*) nft_sc,etot
1266 call write_pdb(2,'second structure',etot)
1268 write(iout,*) 'calling sc_move 2nd time'
1270 call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1271 write(iout,*) nft_sc,etot
1272 call write_pdb(3,'last structure',etot)
1274 end subroutine test_sc
1275 !-----------------------------------------------------------------------------
1276 subroutine bgrow(bstrand,nbstrand,in,ind,new)
1278 ! implicit real*8 (a-h,o-z)
1279 ! include 'DIMENSIONS'
1280 ! include 'COMMON.CHAIN'
1281 integer,dimension(nres/3,6) :: bstrand !(maxres/3,6)
1284 integer :: nbstrand,in,ind,new,ishift,i
1286 ishift=iabs(bstrand(in,ind+4)-new)
1288 print *,'bgrow',bstrand(in,ind+4),new,ishift
1293 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1295 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1296 if (bstrand(i,5).lt.bstrand(i,6)) then
1297 bstrand(i,5)=bstrand(i,5)-ishift
1299 bstrand(i,5)=bstrand(i,5)+ishift
1304 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1306 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1307 if (bstrand(i,6).lt.bstrand(i,5)) then
1308 bstrand(i,6)=bstrand(i,6)-ishift
1310 bstrand(i,6)=bstrand(i,6)+ishift
1317 end subroutine bgrow
1318 !-----------------------------------------------------------------------------
1321 use geometry, only:dist
1323 ! implicit real*8 (a-h,o-z)
1324 ! include 'DIMENSIONS'
1326 ! include 'COMMON.GEO'
1327 ! include 'COMMON.CHAIN'
1328 ! include 'COMMON.IOUNITS'
1329 ! include 'COMMON.VAR'
1330 ! include 'COMMON.CONTROL'
1331 ! include 'COMMON.SBRIDGE'
1332 ! include 'COMMON.FFIELD'
1333 ! include 'COMMON.MINIM'
1335 ! include 'COMMON.DISTFIT'
1336 integer :: if(20,nres),nif,ifa(20)
1337 integer :: ibc(0:nres,0:nres),istrand(20)
1338 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1339 integer :: itmp(20,nres)
1340 real(kind=8) :: time0,time1
1341 real(kind=8) :: energy(0:n_ene),ee
1342 real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres)
1344 logical :: debug,ltest,usedbfrag(nres/3)
1345 character(len=50) :: linia
1347 integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1348 integer :: bstrand(nres/3,6),nbstrand
1349 real(kind=8) :: etot
1350 integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1351 in,ind,ifun,nfun,iretcode
1352 !------------------------
1355 !------------------------
1362 call geom_to_var(nvar,vorg)
1363 call secondary2(debug)
1365 if (nbfrag.le.1) return
1368 usedbfrag(i)=.false.
1372 nbetasheet=nbetasheet+1
1374 bstrand(1,1)=bfrag(1,1)
1375 bstrand(1,2)=bfrag(2,1)
1376 bstrand(1,3)=nbetasheet
1378 bstrand(1,5)=bfrag(1,1)
1379 bstrand(1,6)=bfrag(2,1)
1380 do i=bfrag(1,1),bfrag(2,1)
1381 betasheet(i)=nbetasheet
1385 bstrand(2,1)=bfrag(3,1)
1386 bstrand(2,2)=bfrag(4,1)
1387 bstrand(2,3)=nbetasheet
1388 bstrand(2,5)=bfrag(3,1)
1389 bstrand(2,6)=bfrag(4,1)
1391 if (bfrag(3,1).le.bfrag(4,1)) then
1393 do i=bfrag(3,1),bfrag(4,1)
1394 betasheet(i)=nbetasheet
1399 do i=bfrag(4,1),bfrag(3,1)
1400 betasheet(i)=nbetasheet
1407 do while (iused_nbfrag.ne.nbfrag)
1411 IF (.not.usedbfrag(j)) THEN
1413 write (*,*) j,(bfrag(i,j),i=1,4)
1415 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1417 write (*,*) '------------------'
1420 if (bfrag(3,j).le.bfrag(4,j)) then
1421 do i=bfrag(3,j),bfrag(4,j)
1422 if(betasheet(i).eq.nbetasheet) then
1424 do k=bfrag(3,j),bfrag(4,j)
1425 betasheet(k)=nbetasheet
1430 iused_nbfrag=iused_nbfrag+1
1431 do k=bfrag(1,j),bfrag(2,j)
1432 betasheet(k)=nbetasheet
1433 ibetasheet(k)=nbstrand
1435 if (bstrand(in,4).lt.0) then
1436 bstrand(nbstrand,1)=bfrag(2,j)
1437 bstrand(nbstrand,2)=bfrag(1,j)
1438 bstrand(nbstrand,3)=nbetasheet
1439 bstrand(nbstrand,4)=-nbstrand
1440 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1441 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1442 if(bstrand(in,1).lt.bfrag(4,j)) then
1443 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1445 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1446 (bstrand(in,5)-bfrag(4,j))
1448 if(bstrand(in,2).gt.bfrag(3,j)) then
1449 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1451 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1452 (-bstrand(in,6)+bfrag(3,j))
1455 bstrand(nbstrand,1)=bfrag(1,j)
1456 bstrand(nbstrand,2)=bfrag(2,j)
1457 bstrand(nbstrand,3)=nbetasheet
1458 bstrand(nbstrand,4)=nbstrand
1459 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1460 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1461 if(bstrand(in,1).gt.bfrag(3,j)) then
1462 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1464 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1465 (-bstrand(in,5)+bfrag(3,j))
1467 if(bstrand(in,2).lt.bfrag(4,j)) then
1468 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1470 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1471 (bstrand(in,6)-bfrag(4,j))
1476 if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1477 in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1478 do k=bfrag(1,j),bfrag(2,j)
1479 betasheet(k)=nbetasheet
1484 iused_nbfrag=iused_nbfrag+1
1485 do k=bfrag(3,1),bfrag(4,1)
1486 betasheet(k)=nbetasheet
1487 ibetasheet(k)=nbstrand
1489 if (bstrand(in,4).lt.0) then
1490 bstrand(nbstrand,1)=bfrag(4,j)
1491 bstrand(nbstrand,2)=bfrag(3,j)
1492 bstrand(nbstrand,3)=nbetasheet
1493 bstrand(nbstrand,4)=-nbstrand
1494 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1495 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1496 if(bstrand(in,1).lt.bfrag(2,j)) then
1497 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1499 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1500 (bstrand(in,5)-bfrag(2,j))
1502 if(bstrand(in,2).gt.bfrag(1,j)) then
1503 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1505 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1506 (-bstrand(in,6)+bfrag(1,j))
1509 bstrand(nbstrand,1)=bfrag(3,j)
1510 bstrand(nbstrand,2)=bfrag(4,j)
1511 bstrand(nbstrand,3)=nbetasheet
1512 bstrand(nbstrand,4)=nbstrand
1513 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1514 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1515 if(bstrand(in,1).gt.bfrag(1,j)) then
1516 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1518 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1519 (-bstrand(in,5)+bfrag(1,j))
1521 if(bstrand(in,2).lt.bfrag(2,j)) then
1522 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1524 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1525 (bstrand(in,6)-bfrag(2,j))
1532 do i=bfrag(4,j),bfrag(3,j)
1533 if(betasheet(i).eq.nbetasheet) then
1535 do k=bfrag(4,j),bfrag(3,j)
1536 betasheet(k)=nbetasheet
1541 iused_nbfrag=iused_nbfrag+1
1542 do k=bfrag(1,j),bfrag(2,j)
1543 betasheet(k)=nbetasheet
1544 ibetasheet(k)=nbstrand
1546 if (bstrand(in,4).lt.0) then
1547 bstrand(nbstrand,1)=bfrag(1,j)
1548 bstrand(nbstrand,2)=bfrag(2,j)
1549 bstrand(nbstrand,3)=nbetasheet
1550 bstrand(nbstrand,4)=nbstrand
1551 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1552 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1553 if(bstrand(in,1).lt.bfrag(3,j)) then
1554 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1556 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1557 (bstrand(in,5)-bfrag(3,j))
1559 if(bstrand(in,2).gt.bfrag(4,j)) then
1560 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1562 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1563 (-bstrand(in,6)+bfrag(4,j))
1566 bstrand(nbstrand,1)=bfrag(2,j)
1567 bstrand(nbstrand,2)=bfrag(1,j)
1568 bstrand(nbstrand,3)=nbetasheet
1569 bstrand(nbstrand,4)=-nbstrand
1570 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1571 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1572 if(bstrand(in,1).gt.bfrag(4,j)) then
1573 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1575 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1576 (-bstrand(in,5)+bfrag(4,j))
1578 if(bstrand(in,2).lt.bfrag(3,j)) then
1579 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1581 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1582 (bstrand(in,6)-bfrag(3,j))
1587 if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1588 in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1589 do k=bfrag(1,j),bfrag(2,j)
1590 betasheet(k)=nbetasheet
1595 iused_nbfrag=iused_nbfrag+1
1596 do k=bfrag(4,j),bfrag(3,j)
1597 betasheet(k)=nbetasheet
1598 ibetasheet(k)=nbstrand
1600 if (bstrand(in,4).lt.0) then
1601 bstrand(nbstrand,1)=bfrag(4,j)
1602 bstrand(nbstrand,2)=bfrag(3,j)
1603 bstrand(nbstrand,3)=nbetasheet
1604 bstrand(nbstrand,4)=nbstrand
1605 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1606 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1607 if(bstrand(in,1).lt.bfrag(2,j)) then
1608 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1610 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1611 (bstrand(in,5)-bfrag(2,j))
1613 if(bstrand(in,2).gt.bfrag(1,j)) then
1614 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1616 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1617 (-bstrand(in,6)+bfrag(1,j))
1620 bstrand(nbstrand,1)=bfrag(3,j)
1621 bstrand(nbstrand,2)=bfrag(4,j)
1622 bstrand(nbstrand,3)=nbetasheet
1623 bstrand(nbstrand,4)=-nbstrand
1624 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1625 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1626 if(bstrand(in,1).gt.bfrag(1,j)) then
1627 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1629 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1630 (-bstrand(in,5)+bfrag(1,j))
1632 if(bstrand(in,2).lt.bfrag(2,j)) then
1633 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1635 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1636 (bstrand(in,6)-bfrag(2,j))
1650 do while (usedbfrag(j))
1655 nbetasheet=nbetasheet+1
1656 bstrand(nbstrand,1)=bfrag(1,j)
1657 bstrand(nbstrand,2)=bfrag(2,j)
1658 bstrand(nbstrand,3)=nbetasheet
1659 bstrand(nbstrand,5)=bfrag(1,j)
1660 bstrand(nbstrand,6)=bfrag(2,j)
1662 bstrand(nbstrand,4)=nbstrand
1663 do i=bfrag(1,j),bfrag(2,j)
1664 betasheet(i)=nbetasheet
1665 ibetasheet(i)=nbstrand
1669 bstrand(nbstrand,1)=bfrag(3,j)
1670 bstrand(nbstrand,2)=bfrag(4,j)
1671 bstrand(nbstrand,3)=nbetasheet
1672 bstrand(nbstrand,5)=bfrag(3,j)
1673 bstrand(nbstrand,6)=bfrag(4,j)
1675 if (bfrag(3,j).le.bfrag(4,j)) then
1676 bstrand(nbstrand,4)=nbstrand
1677 do i=bfrag(3,j),bfrag(4,j)
1678 betasheet(i)=nbetasheet
1679 ibetasheet(i)=nbstrand
1682 bstrand(nbstrand,4)=-nbstrand
1683 do i=bfrag(4,j),bfrag(3,j)
1684 betasheet(i)=nbetasheet
1685 ibetasheet(i)=nbstrand
1689 iused_nbfrag=iused_nbfrag+1
1695 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1702 if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1706 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1709 !------------------------
1713 if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1714 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1716 ifb(nifb,1)=bstrand(i,4)
1717 ifb(nifb,2)=bstrand(j,4)
1724 write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1730 write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1732 nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1734 if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1735 nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1741 if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1742 if (if(j,i).gt.0) then
1743 if(betasheet(if(j,i)).eq.0 .or. &
1744 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1749 write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1752 ! read (inp,*) (ifa(i),i=1,4)
1754 ! read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1758 !------------------------
1763 !ccccccccccccccccccccccccccccccccc
1765 !ccccccccccccccccccccccccccccccccc
1769 istrand(is-j+1)=int(ii/is**(is-j))
1770 ii=ii-istrand(is-j+1)*is**(is-j)
1774 istrand(k)=istrand(k)+1
1775 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1779 if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1780 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1789 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1790 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1791 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1792 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1798 if (mod(isa,2).eq.0) then
1800 if (istrand(k).eq.1) ltest=.false.
1803 do k=(isa+1)/2+1,isa
1804 if (istrand(k).eq.1) ltest=.false.
1808 IF (ltest.and.lifb0.eq.1) THEN
1811 call var_to_geom(nvar,vorg)
1813 write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1814 write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1815 write (linia,'(10i3)') (istrand(k),k=1,isa)
1825 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1827 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1831 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1837 write(*,*) (itmp(j,i),j=1,4)
1841 ! ifa(1),ifa(2),ifa(3),ifa(4)
1842 ! if(1,i),if(2,i),if(3,i),if(4,i)
1847 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1848 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1849 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1850 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1858 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1860 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1864 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1866 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1869 !------------------------
1872 ! freeze sec.elements
1882 do i=bfrag(1,j),bfrag(2,j)
1887 if (bfrag(3,j).le.bfrag(4,j)) then
1888 do i=bfrag(3,j),bfrag(4,j)
1894 do i=bfrag(4,j),bfrag(3,j)
1902 do i=hfrag(1,j),hfrag(2,j)
1910 !------------------------
1911 ! generate constrains
1919 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1927 else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1935 else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1943 else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1951 else if ( ibc(i,j).gt.0 ) then
1952 d0(ind)=DIST(i,ibc(i,j))
1959 else if ( ibc(j,i).gt.0 ) then
1960 d0(ind)=DIST(ibc(j,i),j)
1974 !d--------------------------
1976 write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1977 ibc(jhpb(i),ihpb(i)),' --',&
1978 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1985 call contact_cp_min(varia,ifun,iconf,linia,debug)
1990 call minimize(etot,varia,iretcode,nfun)
1991 write(iout,*)'------------------------------------------------'
1992 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1998 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1999 nfun/(time1-time0),' eval/s'
2001 write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
2002 call var_to_geom(nvar,varia)
2004 call write_pdb(900+iconf,linia,etot)
2009 call enerprint(energy)
2011 !d call briefout(0,etot)
2012 !d call secondary2(.true.)
2016 !ccccccccccccccccccccccccccccccccccc
2019 !ccccccccccccccccccccccccccccccccccc
2022 10 write (iout,'(a)') 'Error reading test structure.'
2024 end subroutine test11
2025 !-----------------------------------------------------------------------------
2028 use geometry, only:dist
2030 ! implicit real*8 (a-h,o-z)
2031 ! include 'DIMENSIONS'
2033 ! include 'COMMON.GEO'
2034 ! include 'COMMON.CHAIN'
2035 ! include 'COMMON.IOUNITS'
2036 ! include 'COMMON.VAR'
2037 ! include 'COMMON.CONTROL'
2038 ! include 'COMMON.SBRIDGE'
2039 ! include 'COMMON.FFIELD'
2040 ! include 'COMMON.MINIM'
2042 ! include 'COMMON.DISTFIT'
2043 integer :: if(3,nres),nif
2044 integer :: ibc(nres,nres),istrand(20)
2045 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2046 real(kind=8) :: time0,time1
2047 real(kind=8) :: energy(0:n_ene),ee
2048 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2050 logical :: debug,ltest
2051 character(len=50) :: linia
2052 integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2053 real(kind=8) :: etot
2056 read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2059 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2063 !------------------------
2064 call secondary2(debug)
2065 !------------------------
2073 ! freeze sec.elements and store indexes for beta constrains
2083 do i=bfrag(1,j),bfrag(2,j)
2088 if (bfrag(3,j).le.bfrag(4,j)) then
2089 do i=bfrag(3,j),bfrag(4,j)
2093 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2096 do i=bfrag(4,j),bfrag(3,j)
2100 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2105 do i=hfrag(1,j),hfrag(2,j)
2114 ! ---------------- test --------------
2116 if (ibc(if(1,i),if(2,i)).eq.-1) then
2117 ibc(if(1,i),if(2,i))=if(3,i)
2118 ibc(if(1,i),if(3,i))=if(2,i)
2119 else if (ibc(if(2,i),if(1,i)).eq.-1) then
2120 ibc(if(2,i),if(1,i))=0
2121 ibc(if(1,i),if(2,i))=if(3,i)
2122 ibc(if(1,i),if(3,i))=if(2,i)
2124 ibc(if(1,i),if(2,i))=if(3,i)
2125 ibc(if(1,i),if(3,i))=if(2,i)
2131 if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j)
2134 !------------------------
2140 if ( ibc(i,j).eq.-1 ) then
2148 else if ( ibc(i,j).gt.0 ) then
2149 d0(ind)=DIST(i,ibc(i,j))
2156 else if ( ibc(j,i).gt.0 ) then
2157 d0(ind)=DIST(ibc(j,i),j)
2171 !d--------------------------
2172 write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2173 ibc(jhpb(i),ihpb(i)),' --',&
2174 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2182 call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2187 call minimize(etot,varia,iretcode,nfun)
2188 write(iout,*)'------------------------------------------------'
2189 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2195 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2196 nfun/(time1-time0),' eval/s'
2199 call var_to_geom(nvar,varia)
2201 call write_pdb(999,'full min',etot)
2206 call enerprint(energy)
2208 call briefout(0,etot)
2209 call secondary2(.true.)
2212 10 write (iout,'(a)') 'Error reading test structure.'
2214 end subroutine test3
2215 !-----------------------------------------------------------------------------
2219 ! implicit real*8 (a-h,o-z)
2220 ! include 'DIMENSIONS'
2222 ! include 'COMMON.GEO'
2223 ! include 'COMMON.CHAIN'
2224 ! include 'COMMON.IOUNITS'
2225 ! include 'COMMON.VAR'
2226 ! include 'COMMON.CONTROL'
2227 ! include 'COMMON.SBRIDGE'
2228 ! include 'COMMON.FFIELD'
2229 ! include 'COMMON.MINIM'
2231 ! include 'COMMON.DISTFIT'
2232 integer :: if(2,2),ind
2233 integer :: iff(nres)
2234 real(kind=8) :: time0,time1
2235 real(kind=8) :: energy(0:n_ene),ee
2236 real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2237 theta1,phi1,alph1,omeg1 !(maxres)
2238 real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres)
2240 integer :: i,j,nn,ifun,iretcode,nfun
2241 real(kind=8) :: etot
2244 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2245 write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2246 read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2247 read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2248 read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2249 read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2251 theta2(i)=deg2rad*theta2(i)
2252 phi2(i)=deg2rad*phi2(i)
2253 alph2(i)=deg2rad*alph2(i)
2254 omeg2(i)=deg2rad*omeg2(i)
2268 !------------------------
2273 do i=if(j,1),if(j,2)
2279 call geom_to_var(nvar,varia)
2280 call write_pdb(1,'first structure',0d0)
2282 call secondary(.true.)
2284 call secondary2(.true.)
2287 if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2288 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2289 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2292 if (bfrag(3,j).lt.bfrag(4,j)) then
2293 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2294 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2295 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2297 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2298 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2299 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2312 call geom_to_var(nvar,varia2)
2313 call write_pdb(2,'second structure',0d0)
2317 !-------------------------------------------------------
2320 call contact_cp(varia,varia2,iff,ifun,7)
2325 call minimize(etot,varia,iretcode,nfun)
2326 write(iout,*)'------------------------------------------------'
2327 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2333 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2334 nfun/(time1-time0),' eval/s'
2337 call var_to_geom(nvar,varia)
2339 call write_pdb(999,'full min',etot)
2344 call enerprint(energy)
2346 call briefout(0,etot)
2349 10 write (iout,'(a)') 'Error reading test structure.'
2351 end subroutine test__
2352 !-----------------------------------------------------------------------------
2353 subroutine secondary(lprint)
2355 ! implicit real*8 (a-h,o-z)
2356 ! include 'DIMENSIONS'
2357 ! include 'COMMON.CHAIN'
2358 ! include 'COMMON.IOUNITS'
2359 ! include 'COMMON.DISTFIT'
2361 integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2362 logical :: lprint,not_done
2363 real(kind=4) :: dcont(nres*nres/2),d
2364 real(kind=4) :: rcomp = 7.0
2365 real(kind=4) :: rbeta = 5.2
2366 real(kind=4) :: ralfa = 5.2
2367 real(kind=4) :: r310 = 6.6
2368 real(kind=8),dimension(3) :: xpi,xpj
2369 integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2371 call chainbuild_cart
2372 !d call write_pdb(99,'sec structure',0d0)
2384 xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2388 xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2390 !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2391 !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2392 !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2393 !d print *,'CA',i,j,d
2394 d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2395 (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2396 (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2397 if ( d.lt.rcomp*rcomp) then
2401 dcont(ncont)=sqrt(d)
2407 write (iout,'(a)') '#PP contact map distances:'
2409 write (iout,'(3i4,f10.5)') &
2410 i,icont(1,i),icont(2,i),dcont(i)
2414 ! finding parallel beta
2415 !d write (iout,*) '------- looking for parallel beta -----------'
2421 if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2422 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2423 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2424 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2425 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2426 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2430 !d write (iout,*) i1,j1,dcont(i)
2436 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2437 .and. dcont(j).le.rbeta .and. &
2438 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2439 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2440 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2441 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2442 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2447 !d write (iout,*) i1,j1,dcont(j),not_done
2451 if (i1-ii1.gt.1) then
2455 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2464 isec(ij,1)=isec(ij,1)+1
2465 isec(ij,1+isec(ij,1))=nbeta
2468 isec(ij,1)=isec(ij,1)+1
2469 isec(ij,1+isec(ij,1))=nbeta
2474 if (nbeta.le.9) then
2475 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2476 "DefPropRes 'strand",nstrand,&
2477 "' 'num = ",ii1-1,"..",i1-1,"'"
2479 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2480 "DefPropRes 'strand",nstrand,&
2481 "' 'num = ",ii1-1,"..",i1-1,"'"
2484 if (nbeta.le.9) then
2485 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2486 "DefPropRes 'strand",nstrand,&
2487 "' 'num = ",jj1-1,"..",j1-1,"'"
2489 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2490 "DefPropRes 'strand",nstrand,&
2491 "' 'num = ",jj1-1,"..",j1-1,"'"
2493 write(12,'(a8,4i4)') &
2494 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2500 ! finding antiparallel beta
2501 !d write (iout,*) '--------- looking for antiparallel beta ---------'
2506 if (dcont(i).le.rbeta.and. &
2507 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2508 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2509 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2510 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2511 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2515 !d write (iout,*) i1,j1,dcont(i)
2522 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
2523 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2524 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2525 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2526 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2527 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2528 .and. dcont(j).le.rbeta ) goto 6
2532 !d write (iout,*) i1,j1,dcont(j),not_done
2536 if (i1-ii1.gt.1) then
2537 if(lprint)write (iout,*)'antiparallel beta',&
2538 nbeta,ii1-1,i1,jj1,j1-1
2541 bfrag(1,nbfrag)=max0(ii1-1,1)
2544 bfrag(4,nbfrag)=max0(j1-1,1)
2549 isec(ij,1)=isec(ij,1)+1
2550 isec(ij,1+isec(ij,1))=nbeta
2554 isec(ij,1)=isec(ij,1)+1
2555 isec(ij,1+isec(ij,1))=nbeta
2561 if (nstrand.le.9) then
2562 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2563 "DefPropRes 'strand",nstrand,&
2564 "' 'num = ",ii1-2,"..",i1-1,"'"
2566 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2567 "DefPropRes 'strand",nstrand,&
2568 "' 'num = ",ii1-2,"..",i1-1,"'"
2571 if (nstrand.le.9) then
2572 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2573 "DefPropRes 'strand",nstrand,&
2574 "' 'num = ",j1-2,"..",jj1-1,"'"
2576 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2577 "DefPropRes 'strand",nstrand,&
2578 "' 'num = ",j1-2,"..",jj1-1,"'"
2580 write(12,'(a8,4i4)') &
2581 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2587 if (nstrand.gt.0.and.lprint) then
2588 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2591 write(12,'(a9,i1,$)') " | strand",i
2593 write(12,'(a9,i2,$)') " | strand",i
2596 write(12,'(a1)') "'"
2600 ! finding alpha or 310 helix
2606 if (j1.eq.i1+3.and.dcont(i).le.r310 &
2607 .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2608 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2609 !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2612 if (isec(ii1,1).eq.0) then
2621 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2625 !d write (iout,*) i1,j1,not_done
2628 if (j1-ii1.gt.4) then
2630 !d write (iout,*)'helix',nhelix,ii1,j1
2634 hfrag(2,nhfrag)=max0(j1-1,1)
2640 write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2641 if (nhelix.le.9) then
2642 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2643 "DefPropRes 'helix",nhelix,&
2644 "' 'num = ",ii1-1,"..",j1-2,"'"
2646 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2647 "DefPropRes 'helix",nhelix,&
2648 "' 'num = ",ii1-1,"..",j1-2,"'"
2655 if (nhelix.gt.0.and.lprint) then
2656 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2658 if (nhelix.le.9) then
2659 write(12,'(a8,i1,$)') " | helix",i
2661 write(12,'(a8,i2,$)') " | helix",i
2664 write(12,'(a1)') "'"
2668 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2669 write(12,'(a20)') "XMacStand ribbon.mac"
2673 end subroutine secondary
2674 !-----------------------------------------------------------------------------
2675 subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2677 use geometry, only:dist
2679 ! implicit real*8 (a-h,o-z)
2680 ! include 'DIMENSIONS'
2682 ! include 'COMMON.SBRIDGE'
2683 ! include 'COMMON.FFIELD'
2684 ! include 'COMMON.IOUNITS'
2685 ! include 'COMMON.DISTFIT'
2686 ! include 'COMMON.VAR'
2687 ! include 'COMMON.CHAIN'
2688 ! include 'COMMON.MINIM'
2690 character(len=50) :: linia
2692 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2693 real(kind=8) :: time0,time1
2694 integer :: iff(nres),ieval
2695 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2698 integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2699 real(kind=8) :: wstrain0,etot
2701 maxres22=nres*(nres+1)/2
2703 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2704 call var_to_geom(nvar,var)
2711 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2733 call var_to_geom(nvar,var2)
2736 if ( iff(i).eq.1 ) then
2745 !d call write_pdb(3,'combined structure',0d0)
2746 !d time0=MPI_WTIME()
2749 NY=((NRES-4)*(NRES-5))/2
2750 call distfit(.true.,200)
2752 !d time1=MPI_WTIME()
2753 !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2763 call geom_to_var(nvar,var)
2764 !d time0=MPI_WTIME()
2765 call minimize(etot,var,iretcode,nfun)
2766 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2768 !d time1=MPI_WTIME()
2769 !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2770 !d & nfun/(time1-time0),' SOFT eval/s'
2771 call var_to_geom(nvar,var)
2777 if (iff(1).eq.1) then
2783 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2788 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2794 if (iff(nres).eq.1) then
2800 !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2801 !d & "select",ij(1),"-",ij(2),
2802 !d & ",",ij(3),"-",ij(4)
2803 !d call write_pdb(in_pdb,linia,etot)
2809 !d time0=MPI_WTIME()
2810 call minimize(etot,var,iretcode,nfun)
2811 !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2814 !d time1=MPI_WTIME()
2815 !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2816 !d & nfun/(time1-time0),' eval/s'
2817 !d call var_to_geom(nvar,var)
2819 !d call write_pdb(6,'dist structure',etot)
2828 end subroutine contact_cp2
2829 !-----------------------------------------------------------------------------
2830 subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2832 use geometry, only:dist
2834 ! implicit real*8 (a-h,o-z)
2835 ! include 'DIMENSIONS'
2836 ! include 'COMMON.SBRIDGE'
2837 ! include 'COMMON.FFIELD'
2838 ! include 'COMMON.IOUNITS'
2839 ! include 'COMMON.DISTFIT'
2840 ! include 'COMMON.VAR'
2841 ! include 'COMMON.CHAIN'
2842 ! include 'COMMON.MINIM'
2844 character(len=50) :: linia
2846 real(kind=8) :: energy(0:n_ene)
2847 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2848 real(kind=8) :: time0,time1
2849 integer :: iff(nres),ieval
2850 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2854 integer :: in_pdb,i,j,ind,iwsk
2858 if (ieval.eq.-1) debug=.true.
2862 ! store selected dist. constrains from 1st structure
2865 ! Intercept NaNs in the coordinates
2866 ! write(iout,*) (var(i),i=1,nvar)
2871 if (x_sum.ne.x_sum) then
2872 write(iout,*)" *** contact_cp : Found NaN in coordinates"
2874 print *," *** contact_cp : Found NaN in coordinates"
2880 call var_to_geom(nvar,var)
2887 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2910 ! freeze sec.elements from 2nd structure
2918 call var_to_geom(nvar,var2)
2919 call secondary2(debug)
2921 do i=bfrag(1,j),bfrag(2,j)
2926 if (bfrag(3,j).le.bfrag(4,j)) then
2927 do i=bfrag(3,j),bfrag(4,j)
2933 do i=bfrag(4,j),bfrag(3,j)
2941 do i=hfrag(1,j),hfrag(2,j)
2950 ! copy selected res from 1st to 2nd structure
2954 if ( iff(i).eq.1 ) then
2964 ! prepare description in linia variable
2968 if (iff(1).eq.1) then
2974 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2979 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2985 if (iff(nres).eq.1) then
2990 write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2991 "SELECT",ij(1)-1,"-",ij(2)-1,&
2992 ",",ij(3)-1,"-",ij(4)-1
2998 call contact_cp_min(var,ieval,in_pdb,linia,debug)
3001 end subroutine contact_cp
3002 !-----------------------------------------------------------------------------
3003 subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3007 ! input : theta,phi,alph,omeg,in_pdb,linia,debug
3008 ! output : var,ieval
3010 ! implicit real*8 (a-h,o-z)
3011 ! include 'DIMENSIONS'
3013 ! include 'COMMON.SBRIDGE'
3014 ! include 'COMMON.FFIELD'
3015 ! include 'COMMON.IOUNITS'
3016 ! include 'COMMON.DISTFIT'
3017 ! include 'COMMON.VAR'
3018 ! include 'COMMON.CHAIN'
3019 ! include 'COMMON.MINIM'
3021 character(len=50) :: linia
3023 real(kind=8) :: energy(0:n_ene)
3024 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3025 real(kind=8) :: time0,time1
3026 integer :: ieval,info(3)
3027 logical :: debug,fail,reduce,change !check_var,
3030 integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3032 real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3033 wtor_d01,wstrain0,etot
3035 write(iout,'(a20,i6,a20)') &
3036 '------------------',in_pdb,'-------------------'
3040 call write_pdb(1000+in_pdb,'combined structure',0d0)
3047 ! run optimization of distances
3049 ! uses d0(),w() and mask() for frozen 2D
3051 !test---------------------------------------------
3053 !test NY=((NRES-4)*(NRES-5))/2
3054 !test call distfit(debug,5000)
3086 call geom_to_var(nvar,var)
3087 !de change=reduce(var)
3088 if (check_var(var,info)) then
3089 write(iout,*) 'cp_min error in input'
3090 print *,'cp_min error in input'
3094 !d call etotal(energy(0))
3095 !d call enerprint(energy(0))
3099 !dtest call minimize(etot,var,iretcode,nfun)
3100 !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3103 !d call etotal(energy(0))
3104 !d call enerprint(energy(0))
3123 !test--------------------------------------------------
3129 write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3130 call write_pdb(2000+in_pdb,'distfit structure',0d0)
3138 ! run soft pot. optimization
3140 ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3142 ! mask_phi(),mask_theta(),mask_side(),mask_r
3148 !de change=reduce(var)
3149 !de if (check_var(var,info)) write(iout,*) 'error before soft'
3153 call minimize(etot,var,iretcode,nfun)
3155 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3159 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3160 nfun/(time1-time0),' SOFT eval/s'
3163 call var_to_geom(nvar,var)
3165 call write_pdb(3000+in_pdb,'soft structure',etot)
3169 ! run full UNRES optimization with constrains and frozen 2D
3170 ! the same variables as soft pot. optimizatio
3176 ! check overlaps before calling full UNRES minim
3178 call var_to_geom(nvar,var)
3182 write(iout,*) 'N7 ',energy(0)
3183 if (energy(0).ne.energy(0)) then
3184 write(iout,*) 'N7 error - gives NaN',energy(0)
3188 if (energy(1).eq.1.0d20) then
3189 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3190 call overlap_sc(fail)
3194 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3206 !dte time0=MPI_WTIME()
3207 !de change=reduce(var)
3208 !de if (check_var(var,info)) then
3209 !de write(iout,*) 'error before mask dist'
3210 !de call var_to_geom(nvar,var)
3212 !de call write_pdb(10000+in_pdb,'before mask dist',etot)
3214 !dte call minimize(etot,var,iretcode,nfun)
3215 !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3216 !dte & ' eval ',nfun
3217 !dte ieval=ieval+nfun
3219 !dte time1=MPI_WTIME()
3220 !dte write (iout,'(a,f6.2,f8.2,a)')
3221 !dte & ' Time for mask dist min.',time1-time0,
3222 !dte & nfun/(time1-time0),' eval/s'
3223 !dte call flush(iout)
3226 call var_to_geom(nvar,var)
3228 call write_pdb(4000+in_pdb,'mask dist',etot)
3231 ! switch off freezing of 2D and
3232 ! run full UNRES optimization with constrains
3238 !de change=reduce(var)
3239 !de if (check_var(var,info)) then
3240 !de write(iout,*) 'error before dist'
3241 !de call var_to_geom(nvar,var)
3243 !de call write_pdb(11000+in_pdb,'before dist',etot)
3246 call minimize(etot,var,iretcode,nfun)
3248 !de change=reduce(var)
3249 !de if (check_var(var,info)) then
3250 !de write(iout,*) 'error after dist',ico
3251 !de call var_to_geom(nvar,var)
3253 !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3255 write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3261 write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3262 nfun/(time1-time0),' eval/s'
3264 !de call etotal(energy(0))
3265 !de write(iout,*) 'N7 after dist',energy(0)
3269 call var_to_geom(nvar,var)
3271 call write_pdb(in_pdb,linia,etot)
3283 end subroutine contact_cp_min
3284 !-----------------------------------------------------------------------------
3287 use geometry, only:dist
3289 ! implicit real*8 (a-h,o-z)
3290 ! include 'DIMENSIONS'
3292 ! include 'COMMON.GEO'
3293 ! include 'COMMON.CHAIN'
3294 ! include 'COMMON.IOUNITS'
3295 ! include 'COMMON.VAR'
3296 ! include 'COMMON.CONTROL'
3297 ! include 'COMMON.SBRIDGE'
3298 ! include 'COMMON.FFIELD'
3299 ! include 'COMMON.MINIM'
3300 ! include 'COMMON.INTERACT'
3302 ! include 'COMMON.DISTFIT'
3303 integer :: iff(nres)
3304 real(kind=8) :: time0,time1
3305 real(kind=8) :: energy(0:n_ene),ee
3306 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3308 logical :: debug,ltest,fail
3309 character(len=50) :: linia
3310 integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3312 real(kind=8) :: wstrain0,wang0,etot
3318 !------------------------
3320 ! freeze sec.elements
3330 do i=bfrag(1,j),bfrag(2,j)
3335 if (bfrag(3,j).le.bfrag(4,j)) then
3336 do i=bfrag(3,j),bfrag(4,j)
3342 do i=bfrag(4,j),bfrag(3,j)
3350 do i=hfrag(1,j),hfrag(2,j)
3362 ! store dist. constrains
3366 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3371 dhpb(nhpb)=DIST(i,j)
3379 call write_pdb(100+in_pdb,'input reg. structure',0d0)
3389 ! run soft pot. optimization
3395 call geom_to_var(nvar,var)
3399 call minimize(etot,var,iretcode,nfun)
3401 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3405 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3406 nfun/(time1-time0),' SOFT eval/s'
3408 call var_to_geom(nvar,var)
3410 call write_pdb(300+in_pdb,'soft structure',etot)
3413 ! run full UNRES optimization with constrains and frozen 2D
3414 ! the same variables as soft pot. optimizatio
3423 call minimize(etot,var,iretcode,nfun)
3424 write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3431 write (iout,'(a,f6.2,f8.2,a)') &
3432 ' Time for mask dist min.',time1-time0,&
3433 nfun/(time1-time0),' eval/s'
3435 call var_to_geom(nvar,var)
3437 call write_pdb(400+in_pdb,'mask & dist',etot)
3440 ! switch off constrains and
3441 ! run full UNRES optimization with frozen 2D
3456 call minimize(etot,var,iretcode,nfun)
3457 write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3463 write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,&
3464 nfun/(time1-time0),' eval/s'
3468 call var_to_geom(nvar,var)
3470 call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3477 ! run full UNRES optimization with constrains and NO frozen 2D
3487 wstrain=wstrain0/ico
3492 call minimize(etot,var,iretcode,nfun)
3493 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3494 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3501 write (iout,'(a,f6.2,f8.2,a)') &
3502 ' Time for dist min.',time1-time0,&
3503 nfun/(time1-time0),' eval/s'
3505 call var_to_geom(nvar,var)
3507 call write_pdb(600+in_pdb+ico,'dist cons',etot)
3524 call minimize(etot,var,iretcode,nfun)
3525 write(iout,*)'------------------------------------------------'
3526 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3532 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
3533 nfun/(time1-time0),' eval/s'
3536 call var_to_geom(nvar,var)
3538 call write_pdb(999,'full min',etot)
3542 end subroutine softreg
3543 !-----------------------------------------------------------------------------
3544 subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3546 use geometry, only:dist
3548 ! implicit real*8 (a-h,o-z)
3549 ! include 'DIMENSIONS'
3551 ! include 'COMMON.GEO'
3552 ! include 'COMMON.VAR'
3553 ! include 'COMMON.INTERACT'
3554 ! include 'COMMON.IOUNITS'
3555 ! include 'COMMON.DISTFIT'
3556 ! include 'COMMON.SBRIDGE'
3557 ! include 'COMMON.CONTROL'
3558 ! include 'COMMON.FFIELD'
3559 ! include 'COMMON.MINIM'
3560 ! include 'COMMON.CHAIN'
3561 real(kind=8) :: time0,time1
3562 real(kind=8) :: energy(0:n_ene),ee
3563 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3564 integer :: jdata(5),isec(nres)
3567 integer :: i1,i2,i3,i4,i5,ieval,ij
3568 integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3569 real(kind=8) :: etot,wscloc0,wstrain0
3577 call secondary2(.false.)
3583 do i=bfrag(1,j),bfrag(2,j)
3586 do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3591 do i=hfrag(1,j),hfrag(2,j)
3597 ! cut strands at the ends
3599 if (jdata(2)-jdata(1).gt.3) then
3602 if (jdata(3).lt.jdata(4)) then
3612 !v call etotal(energy(0))
3614 !v write(iout,*) nnt,nct,etot
3615 !v call write_pdb(ij*100,'first structure',etot)
3616 !v write(iout,*) 'N16 test',(jdata(i),i=1,5)
3618 !------------------------
3619 ! generate constrains
3622 if(ishift.eq.0) ishift=-2
3625 do i=jdata(1),jdata(2)
3627 if(jdata(4).gt.jdata(3))then
3628 do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3630 !d print *,i,j,j+ishift
3635 dhpb(nhpb)=DIST(i,j+ishift)
3638 do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3640 !d print *,i,j,j+ishift
3645 dhpb(nhpb)=DIST(i,j+ishift)
3652 if(isec(i).gt.0.or.isec(j).gt.0) then
3658 dhpb(nhpb)=DIST(i,j)
3665 call geom_to_var(nvar,var)
3672 wstrain=wstrain0/ico
3674 !v time0=MPI_WTIME()
3675 call minimize(etot,var,iretcode,nfun)
3676 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3677 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3680 !v time1=MPI_WTIME()
3681 !v write (iout,'(a,f6.2,f8.2,a)')
3682 !v & ' Time for dist min.',time1-time0,
3683 !v & nfun/(time1-time0),' eval/s'
3684 !v call var_to_geom(nvar,var)
3686 !v call write_pdb(ij*100+ico,'dist cons',etot)
3698 call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3701 !v call etotal(energy(0))
3703 !v call write_pdb(ij*100+10,'sc_move',etot)
3705 !d print *,nft_sc,etot
3708 end subroutine beta_slide
3709 !-----------------------------------------------------------------------------
3710 subroutine beta_zip(i1,i2,ieval,ij)
3713 ! implicit real*8 (a-h,o-z)
3714 ! include 'DIMENSIONS'
3716 ! include 'COMMON.GEO'
3717 ! include 'COMMON.VAR'
3718 ! include 'COMMON.INTERACT'
3719 ! include 'COMMON.IOUNITS'
3720 ! include 'COMMON.DISTFIT'
3721 ! include 'COMMON.SBRIDGE'
3722 ! include 'COMMON.CONTROL'
3723 ! include 'COMMON.FFIELD'
3724 ! include 'COMMON.MINIM'
3725 ! include 'COMMON.CHAIN'
3726 real(kind=8) :: time0,time1
3727 real(kind=8) :: energy(0:n_ene),ee
3728 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3729 character(len=10) :: test
3731 integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3732 real(kind=8) :: etot,wstrain0
3734 !v call etotal(energy(0))
3736 !v write(test,'(2i5)') i1,i2
3737 !v call write_pdb(ij*100,test,etot)
3738 !v write(iout,*) 'N17 test',i1,i2,etot,ij
3741 ! generate constrains
3752 call geom_to_var(nvar,var)
3758 wstrain=wstrain0/ico
3759 !v time0=MPI_WTIME()
3760 call minimize(etot,var,iretcode,nfun)
3761 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3762 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3765 !v time1=MPI_WTIME()
3766 !v write (iout,'(a,f6.2,f8.2,a)')
3767 !v & ' Time for dist min.',time1-time0,
3768 !v & nfun/(time1-time0),' eval/s'
3769 ! do not comment the next line
3770 call var_to_geom(nvar,var)
3772 !v call write_pdb(ij*100+ico,'dist cons',etot)
3780 !v call etotal(energy(0))
3782 !v write(iout,*) 'N17 test end',i1,i2,etot,ij
3785 end subroutine beta_zip
3786 !-----------------------------------------------------------------------------
3788 !-----------------------------------------------------------------------------
3789 subroutine thread_seq
3791 use geometry, only:dist
3792 use random, only:iran_num
3793 use control, only:tcpu
3794 use regularize_, only:regularize
3795 use mcm_data, only: nsave_part,nacc_tot
3796 ! Thread the sequence through a database of known structures
3797 ! implicit real*8 (a-h,o-z)
3798 ! include 'DIMENSIONS'
3799 use MPI_data !include 'COMMON.INFO'
3804 ! include 'COMMON.CONTROL'
3805 ! include 'COMMON.CHAIN'
3806 ! include 'COMMON.DBASE'
3807 ! include 'COMMON.INTERACT'
3808 ! include 'COMMON.VAR'
3809 ! include 'COMMON.THREAD'
3810 ! include 'COMMON.FFIELD'
3811 ! include 'COMMON.SBRIDGE'
3812 ! include 'COMMON.HEADER'
3813 ! include 'COMMON.IOUNITS'
3814 ! include 'COMMON.TIME1'
3815 ! include 'COMMON.CONTACTS'
3816 ! include 'COMMON.MCM'
3817 ! include 'COMMON.NAMES'
3819 integer :: ThreadId,ThreadType,Kwita
3821 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3822 real(kind=8) :: przes(3),obr(3,3)
3823 real(kind=8) :: time_for_thread
3824 logical :: found_pattern,non_conv
3825 character(len=32) :: head_pdb
3826 real(kind=8) :: energia(0:n_ene)
3827 integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3829 real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3831 n_ene_comp=nprint_ene
3836 if (me.eq.king) then
3855 ave_time_for_thread=0.0D0
3856 max_time_for_thread=0.0D0
3857 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3858 nthread=nexcl+nthread
3859 do ithread=1,nthread
3860 found_pattern=.false.
3862 do while (.not.found_pattern)
3864 if (itrial.gt.1000) then
3865 write (iout,'(/a/)') 'Too many attempts to find pattern.'
3868 call recv_stop_sig(Kwita)
3869 call send_stop_sig(-3)
3873 ! Find long enough chain in the database
3875 nres_t=nres_base(1,ii)
3876 ! Select the starting position to thread.
3877 print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3878 nres_t,' nres0=',nres0
3879 if (nres_t.ge.nres0) then
3880 ist=iran_num(0,nres_t-nres0)
3882 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3883 if (Kwita.lt.0) then
3884 write (iout,*) 'Stop signal received. Terminating.'
3885 write (*,*) 'Stop signal received. Terminating.'
3887 write (*,*) 'ithread=',ithread,' nthread=',nthread
3890 call pattern_receive
3893 if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3895 found_pattern=.true.
3897 ! If this point is reached, the pattern has not yet been examined.
3899 ! print *,'found_pattern:',found_pattern
3905 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3906 if (Kwita.lt.0) then
3907 write (iout,*) 'Stop signal received. Terminating.'
3909 write (*,*) 'ithread=',ithread,' nthread=',nthread
3915 ipatt(2,ithread)=ist
3917 write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3918 'Processor:',me,' Attempt:',ithread,&
3919 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3920 ' start at res.',ist+1
3921 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3922 ' Attempt:',ithread,&
3923 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3924 ' start at res.',ist+1
3926 write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3927 'Attempt:',ithread,&
3928 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3929 ' start at res.',ist+1
3930 write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3931 'Attempt:',ithread,&
3932 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3933 ' start at res.',ist+1
3936 ! Copy coordinates from the database.
3940 c(j,i)=cart_base(j,i+ist,ii)
3943 !d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
3945 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3947 !d write (iout,'(a,f10.5)')
3948 !d & 'Initial RMS deviation from reference structure:',rms
3949 if (itype(nres,1).eq.ntyp1) then
3951 dcj=c(j,nres-2)-c(j,nres-3)
3952 c(j,nres)=c(j,nres-1)+dcj
3953 c(j,2*nres)=c(j,nres)
3956 if (itype(1,1).eq.ntyp1) then
3963 call int_from_cart(.false.,.false.)
3964 !d print *,'Exit INT_FROM_CART.'
3965 !d print *,'nhpb=',nhpb
3970 ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3972 ! stop 'End generate'
3973 ! Generate SC conformations.
3977 !d print *,'Processor:',me,': exit GEN_SIDE.'
3979 !d print *,'Exit GEN_SIDE.'
3981 ! Calculate initial energy.
3983 call etotal(energia)
3986 ener0(i,ithread)=energia(i)
3988 ener0(n_ene_comp+1,ithread)=energia(0)
3990 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
3991 ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
3993 ener0(n_ene_comp+2,ithread)=rms
3994 ener0(n_ene_comp+4,ithread)=frac
3995 ener0(n_ene_comp+5,ithread)=frac_nn
3997 ener0(n_ene_comp+3,ithread)=0.0d0
4000 print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
4002 print*,'ithread=',ithread,' Start REGULARIZE.'
4005 call regularize(nct-nnt+1,etot,rms,&
4006 cart_base(1,ist+nnt,ipattern),iretcode)
4008 time_for_thread=curr_tim1-curr_tim
4009 ave_time_for_thread= &
4010 ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4011 if (time_for_thread.gt.max_time_for_thread) &
4012 max_time_for_thread=time_for_thread
4014 print *,'Processor',me,': Exit REGULARIZE.'
4015 if (WhatsUp.eq.2) then
4017 'Sufficient number of confs. collected. Terminating.'
4020 else if (WhatsUp.eq.-1) then
4022 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4023 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4024 call send_stop_sig(-2)
4026 else if (WhatsUp.eq.-2) then
4028 write (iout,*) 'Timeup signal received. Terminating.'
4030 else if (WhatsUp.eq.-3) then
4032 write (iout,*) 'Error stop signal received. Terminating.'
4036 print *,'Exit REGULARIZE.'
4037 if (iretcode.eq.11) then
4038 write (iout,'(/a/)') &
4039 '******* Allocated time exceeded in SUMSL. The program will stop.'
4044 head_pdb=titel(:24)//':'//str_nam(ipattern)
4045 if (outpdb) call pdbout(etot,head_pdb,ipdb)
4046 if (outmol2) call mol2out(etot,head_pdb)
4048 call briefout(ithread,etot)
4050 link_end=min0(link_end,nss)
4051 write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4053 call etotal(energia)
4054 ! call enerprint(energia(0))
4057 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4058 !d write (iout,'(a,f10.5)')
4059 !d & 'RMS deviation from reference structure:',dsqrt(rms)
4061 ener(i,ithread)=energia(i)
4063 ener(n_ene_comp+1,ithread)=energia(0)
4064 ener(n_ene_comp+3,ithread)=rms
4066 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4067 ener(n_ene_comp+2,ithread)=rms
4068 ener(n_ene_comp+4,ithread)=frac
4069 ener(n_ene_comp+5,ithread)=frac_nn
4071 call write_stat_thread(ithread,ipattern,ist)
4072 ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4073 ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4074 ! & (ener(k,ithread),k=12,14)
4076 if (me.eq.king) then
4078 call pattern_receive
4079 call receive_MCM_info
4080 if (nacc_tot.ge.nthread) then
4082 'Sufficient number of conformations collected nacc_tot=',&
4083 nacc_tot,'. Stopping other processors and terminating.'
4085 'Sufficient number of conformations collected nacc_tot=',&
4086 nacc_tot,'. Stopping other processors and terminating.'
4087 call recv_stop_sig(Kwita)
4088 if (Kwita.eq.0) call send_stop_sig(-1)
4093 call send_MCM_info(2)
4096 if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4097 write (iout,'(/2a)') &
4098 '********** There would be not enough time for another thread. ',&
4099 'The program will stop.'
4101 '********** There would be not enough time for another thread. ',&
4102 'The program will stop.'
4103 write (iout,'(a,1pe14.4/)') &
4104 'Elapsed time for last threading step: ',time_for_thread
4107 call recv_stop_sig(Kwita)
4108 call send_stop_sig(-2)
4113 write (iout,'(a,1pe14.4)') &
4114 'Elapsed time for this threading step: ',time_for_thread
4117 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4118 if (Kwita.lt.0) then
4119 write (iout,*) 'Stop signal received. Terminating.'
4120 write (*,*) 'Stop signal received. Terminating.'
4122 write (*,*) 'nthread=',nthread,' ithread=',ithread
4128 call send_stop_sig(-1)
4132 ! Any messages left for me?
4133 call pattern_receive
4134 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4136 call write_thread_summary
4138 if (king.eq.king) then
4140 do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4143 call recv_stop_sig(Kwita)
4144 call receive_MCM_info
4147 call receive_thread_results(iproc)
4149 call write_thread_summary
4151 call send_thread_results
4155 end subroutine thread_seq
4156 !-----------------------------------------------------------------------------
4159 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4161 use random, only:iran_num
4162 ! implicit real*8 (a-h,o-z)
4163 ! include 'DIMENSIONS'
4164 ! include 'COMMON.CHAIN'
4165 ! include 'COMMON.DBASE'
4166 ! include 'COMMON.INTERACT'
4167 ! include 'COMMON.VAR'
4168 ! include 'COMMON.THREAD'
4169 ! include 'COMMON.FFIELD'
4170 ! include 'COMMON.SBRIDGE'
4171 ! include 'COMMON.HEADER'
4172 ! include 'COMMON.GEO'
4173 ! include 'COMMON.IOUNITS'
4174 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
4175 !el integer :: icall
4176 !el common /srutu/ icall
4177 real(kind=8) :: energia(0:n_ene)
4178 logical :: glycine,fail
4179 integer :: i,maxsample,link_end0,ind_sc,isample
4180 real(kind=8) :: alph0,omeg0,e1,e0
4184 link_end=min0(link_end,nss)
4186 if (itype(i,1).ne.10) then
4187 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
4188 call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
4192 call etotal(energia)
4194 do isample=1,maxsample
4195 ! Choose a non-glycine side chain.
4198 ind_sc=iran_num(nnt,nct)
4199 glycine=(itype(ind_sc,1).eq.10)
4203 call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
4206 call etotal(energia)
4207 !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4208 !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4209 !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4220 end subroutine sc_conf
4221 !-----------------------------------------------------------------------------
4223 !-----------------------------------------------------------------------------
4224 logical function check_var(var,info)
4228 ! implicit real*8 (a-h,o-z)
4229 ! include 'DIMENSIONS'
4230 ! include 'COMMON.VAR'
4231 ! include 'COMMON.IOUNITS'
4232 ! include 'COMMON.GEO'
4233 ! include 'COMMON.SETUP'
4234 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4235 integer,dimension(3) :: info
4239 do i=nphi+ntheta+1,nphi+ntheta+nside
4240 ! Check the side chain "valence" angles alpha
4241 if (var(i).lt.1.0d-7) then
4242 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4243 write (iout,*) 'Processor',me,'received bad variables!!!!'
4244 write (iout,*) 'Variables'
4245 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4246 write (iout,*) 'Continuing calculations at this point',&
4247 ' could destroy the results obtained so far... ABORTING!!!!!!'
4248 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4249 'valence angle alpha',i-nphi-ntheta,var(i),&
4250 'n it',info(1),info(2),'mv ',info(3)
4251 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4252 write (*,*) 'Processor',me,'received bad variables!!!!'
4253 write (*,*) 'Variables'
4254 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4255 write (*,*) 'Continuing calculations at this point',&
4256 ' could destroy the results obtained so far... ABORTING!!!!!!'
4257 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4258 'valence angle alpha',i-nphi-ntheta,var(i),&
4259 'n it',info(1),info(2),'mv ',info(3)
4264 ! Check the backbone "valence" angles theta
4265 do i=nphi+1,nphi+ntheta
4266 if (var(i).lt.1.0d-7) then
4267 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4268 write (iout,*) 'Processor',me,'received bad variables!!!!'
4269 write (iout,*) 'Variables'
4270 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4271 write (iout,*) 'Continuing calculations at this point',&
4272 ' could destroy the results obtained so far... ABORTING!!!!!!'
4273 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4274 'valence angle theta',i-nphi,var(i),&
4275 'n it',info(1),info(2),'mv ',info(3)
4276 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4277 write (*,*) 'Processor',me,'received bad variables!!!!'
4278 write (*,*) 'Variables'
4279 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4280 write (*,*) 'Continuing calculations at this point',&
4281 ' could destroy the results obtained so far... ABORTING!!!!!!'
4282 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4283 'valence angle theta',i-nphi,var(i),&
4284 'n it',info(1),info(2),'mv ',info(3)
4290 end function check_var
4291 !-----------------------------------------------------------------------------
4293 !-----------------------------------------------------------------------------
4294 subroutine distfit(debug,maxit)
4296 use geometry_data, only: phi
4299 ! implicit real*8 (a-h,o-z)
4300 ! include 'DIMENSIONS'
4301 ! include 'COMMON.CHAIN'
4302 ! include 'COMMON.VAR'
4303 ! include 'COMMON.IOUNITS'
4304 ! include 'COMMON.DISTFIT'
4305 integer :: i,maxit,MAXMAR,IT,IMAR
4306 real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres)
4307 logical :: debug,sing
4308 real(kind=8) :: TOL,RL,F0,AIN,F1
4310 !input------------------------------------
4312 ! NY=((NRES-4)*(NRES-5))/2
4313 !input------------------------------------
4319 CALL TRANSFER(NRES,phi,phiold)
4323 !d WRITE (IOUT,*) 'DISTFIT: F0=',F0
4339 CALL TRANSFER(NX,XX,X)
4340 CALL BANACH(NX,NRES,H,X,sing)
4345 IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4347 WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4348 WRITE (IOUT,*) 'IT=',it,'F=',F0
4353 phi(I)=phiold(I)+mask(i)*X(I-3)
4358 !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4360 CALL TRANSFER(NRES,phi,phiold)
4363 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4365 WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4366 WRITE (IOUT,*) 'IT=',it,'F=',F1
4372 WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4373 WRITE (IOUT,*) 'IT=',it,'F=',F0
4374 CALL TRANSFER(NRES,phiold,phi)
4377 !d write (iout,*) "it",it," imar",imar," f0",f0
4379 WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4381 end subroutine distfit
4382 !-----------------------------------------------------------------------------
4383 real(kind=8) function RDIF()
4386 use geometry, only: dist
4387 ! implicit real*8 (a-h,o-z)
4388 ! include 'DIMENSIONS'
4389 ! include 'COMMON.CHAIN'
4390 ! include 'COMMON.DISTFIT'
4392 real(kind=8) :: suma,DIJ
4401 if (w(ind).ne.0.0) then
4403 suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4405 ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4413 !-----------------------------------------------------------------------------
4418 use geometry, only:dist
4419 ! implicit real*8 (a-h,o-z)
4420 ! include 'DIMENSIONS'
4421 ! include 'COMMON.CHAIN'
4422 ! include 'COMMON.DISTFIT'
4423 ! include 'COMMON.GEO'
4424 integer :: i,j,k,l,I1,I2,IND
4425 real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4438 R13(K)=C(K,J)-C(K,I1)
4442 R24(L)=C(L,K)-C(L,I2)
4444 IND=((J-1)*(2*NRES-J-6))/2+K-3
4445 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4446 PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4447 PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4448 DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4453 end subroutine RDERIV
4454 !-----------------------------------------------------------------------------
4458 ! implicit real*8 (a-h,o-z)
4459 ! include 'DIMENSIONS'
4460 ! include 'COMMON.CHAIN'
4461 ! include 'COMMON.DISTFIT'
4463 real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4471 XI=XI+BKIWK*(D0(K)-DDD(K))
4479 HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4486 end subroutine HEVAL
4487 !-----------------------------------------------------------------------------
4488 subroutine VEC(I,J,U)
4490 use geometry_data, only: C
4491 ! Find the unit vector from atom (I) to atom (J). Store in U.
4493 ! implicit real*8 (a-h,o-z)
4494 ! include 'DIMENSIONS'
4495 ! include 'COMMON.CHAIN'
4497 real(kind=8),DIMENSION(3) :: U
4498 real(kind=8) :: ANORM,UK
4512 !-----------------------------------------------------------------------------
4513 subroutine TRANSFER(N,X1,X2)
4515 ! implicit real*8 (a-h,o-z)
4516 ! include 'DIMENSIONS'
4518 real(kind=8),DIMENSION(N) :: X1,X2
4522 end subroutine TRANSFER
4523 !-----------------------------------------------------------------------------
4524 !-----------------------------------------------------------------------------
4525 subroutine alloc_compare_arrays
4527 maxres22=nres*(nres+1)/2
4529 ! common /struct/ in io_common: read_threadbase
4530 ! allocate(cart_base !(3,maxres_base,maxseq)
4531 ! allocate(nres_base !(3,maxseq)
4532 ! allocate(str_nam !(maxseq)
4534 ! COMMON /c_frag/ in io_conf: readpdb
4535 if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4536 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4538 allocate(w(maxres22),d0(maxres22)) !(maxres22)
4540 !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4541 allocate(DDD(maxres22)) !(maxres22)
4542 allocate(H(nres,nres)) !(MAXRES,MAXRES)
4543 allocate(XX(nres)) !(MAXRES)
4545 allocate(mask(nres)) !(maxres)
4548 allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4550 allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4553 end subroutine alloc_compare_arrays
4554 !-----------------------------------------------------------------------------
4556 !-----------------------------------------------------------------------------