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
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),i1,restyp(it2),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),i1,restyp(it2),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)),k,k=i1,ii1)
234 write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
236 ! write (iout,'(a,i3,$)') restyp(itype(k)),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),&
304 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
315 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
318 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
319 if (iteli.eq.2 .and. itelj.eq.2) goto 4
320 aaa=app_(iteli,itelj)
321 bbb=bpp_(iteli,itelj)
322 ael6_i=ael6_(iteli,itelj)
323 ael3_i=ael3_(iteli,itelj)
327 xj=c(1,j)+0.5*dxj-xmedi
328 yj=c(2,j)+0.5*dyj-ymedi
329 zj=c(3,j)+0.5*dzj-zmedi
330 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
335 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
336 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
337 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
338 fac=cosa-3.0*cosb*cosg
344 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
347 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
348 j.eq.i+2 .and. eesij.le.elecutoff_14) then
359 write (iout,*) 'Total average electrostatic energy: ',ees
360 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
362 write (iout,*) 'Electrostatic contacts before pruning: '
368 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
369 i,restyp(it1),i1,restyp(it2),i2,econt(i)
372 ! For given residues keep only the contacts with the greatest energy.
374 do while (i.lt.ncont)
380 do while (j.lt.ncont)
382 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
383 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
384 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
385 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
386 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
387 if (ic1.eq.icont(1,j)) then
389 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
390 .and. iabs(icont(1,k)-ic1).le.2 .and. &
391 econt(k).lt.econt(j) ) goto 21
393 else if (ic2.eq.icont(2,j) ) then
395 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
396 .and. iabs(icont(2,k)-ic2).le.2 .and. &
397 econt(k).lt.econt(j) ) goto 21
402 icont(1,k-1)=icont(1,k)
403 icont(2,k-1)=icont(2,k)
408 ! write (iout,*) "ncont",ncont
410 ! write (iout,*) icont(1,k),icont(2,k)
413 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
415 if (ic1.eq.icont(1,j)) then
417 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
418 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
419 econt(k).lt.econt(i) ) goto 21
421 else if (ic2.eq.icont(2,j) ) then
423 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
424 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
425 econt(k).lt.econt(i) ) goto 21
430 icont(1,k-1)=icont(1,k)
431 icont(2,k-1)=icont(2,k)
435 ! write (iout,*) "ncont",ncont
437 ! write (iout,*) icont(1,k),icont(2,k)
448 write (iout,*) 'Electrostatic contacts after pruning: '
454 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
455 i,restyp(it1),i1,restyp(it2),i2,econt(i)
459 end subroutine elecont
460 !-----------------------------------------------------------------------------
461 subroutine secondary2(lprint)
463 ! implicit real*8 (a-h,o-z)
464 ! include 'DIMENSIONS'
465 ! include 'COMMON.CHAIN'
466 ! include 'COMMON.IOUNITS'
467 ! include 'COMMON.DISTFIT'
468 ! include 'COMMON.VAR'
469 ! include 'COMMON.GEO'
470 ! include 'COMMON.CONTROL'
471 integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
473 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
474 integer,dimension(nres,4) :: isec !(maxres,4)
475 integer,dimension(nres) :: nsec !(maxres)
476 logical :: lprint,not_done !,freeres
477 real(kind=8) :: p1,p2
480 !el allocate(icont(2,12*nres),isec(nres,4),nsec(nres))
482 if(.not.dccart) call chainbuild_cart
483 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
484 !d call write_pdb(99,'sec structure',0d0)
494 call elecont(lprint,ncont,icont)
496 ! finding parallel beta
497 !d write (iout,*) '------- looking for parallel beta -----------'
503 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
506 !d write (iout,*) i1,j1
512 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
513 freeres(i1,j1,nsec,isec)) goto 5
517 !d write (iout,*) i1,j1,not_done
521 if (i1-ii1.gt.1) then
525 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
529 bfrag(1,nbfrag)=ii1+1
531 bfrag(3,nbfrag)=jj1+1
532 bfrag(4,nbfrag)=min0(j1+1,nres)
536 isec(ij,nsec(ij))=nbeta
540 isec(ij,nsec(ij))=nbeta
546 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
547 "DefPropRes 'strand",nstrand,&
548 "' 'num = ",ii1-1,"..",i1-1,"'"
550 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
551 "DefPropRes 'strand",nstrand,&
552 "' 'num = ",ii1-1,"..",i1-1,"'"
556 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
557 "DefPropRes 'strand",nstrand,&
558 "' 'num = ",jj1-1,"..",j1-1,"'"
560 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
561 "DefPropRes 'strand",nstrand,&
562 "' 'num = ",jj1-1,"..",j1-1,"'"
564 write(12,'(a8,4i4)') &
565 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
571 ! finding alpha or 310 helix
578 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
581 if (j1.eq.i1+3 .and. &
582 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
583 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
584 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
585 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
588 if (nsec(ii1).eq.0) then
597 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
603 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
608 if (j1-ii1.gt.5) then
620 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
621 if (nhelix.le.9) then
622 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
623 "DefPropRes 'helix",nhelix,&
624 "' 'num = ",ii1-1,"..",j1-2,"'"
626 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
627 "DefPropRes 'helix",nhelix,&
628 "' 'num = ",ii1-1,"..",j1-2,"'"
634 if (nhelix.gt.0.and.lprint) then
635 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
637 if (nhelix.le.9) then
638 write(12,'(a8,i1,$)') " | helix",i
640 write(12,'(a8,i2,$)') " | helix",i
647 ! finding antiparallel beta
648 !d write (iout,*) '--------- looking for antiparallel beta ---------'
653 if (freeres(i1,j1,nsec,isec)) then
656 !d write (iout,*) i1,j1
663 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
664 freeres(i1,j1,nsec,isec)) goto 6
668 !d write (iout,*) i1,j1,not_done
672 if (i1-ii1.gt.1) then
676 bfrag(2,nbfrag)=min0(i1+1,nres)
677 bfrag(3,nbfrag)=min0(jj1+1,nres)
684 if (nsec(ij).le.2) then
685 isec(ij,nsec(ij))=nbeta
691 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
692 isec(ij,nsec(ij))=nbeta
698 write (iout,'(a,i3,4i4)')'antiparallel beta',&
699 nbeta,ii1-1,i1,jj1,j1-1
701 if (nstrand.le.9) then
702 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
703 "DefPropRes 'strand",nstrand,&
704 "' 'num = ",ii1-2,"..",i1-1,"'"
706 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
707 "DefPropRes 'strand",nstrand,&
708 "' 'num = ",ii1-2,"..",i1-1,"'"
711 if (nstrand.le.9) then
712 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
713 "DefPropRes 'strand",nstrand,&
714 "' 'num = ",j1-2,"..",jj1-1,"'"
716 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
717 "DefPropRes 'strand",nstrand,&
718 "' 'num = ",j1-2,"..",jj1-1,"'"
720 write(12,'(a8,4i4)') &
721 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
727 if (nstrand.gt.0.and.lprint) then
728 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
731 write(12,'(a9,i1,$)') " | strand",i
733 write(12,'(a9,i2,$)') " | strand",i
742 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
743 write(12,'(a20)') "XMacStand ribbon.mac"
746 write(iout,*) 'UNRES seq:'
748 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
752 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
757 end subroutine secondary2
759 !-----------------------------------------------------------------------------
760 logical function freeres(i,j,nsec,isec)
762 ! implicit real*8 (a-h,o-z)
763 ! include 'DIMENSIONS'
764 integer,dimension(nres,4) :: isec !(maxres,4)
765 integer,dimension(nres) :: nsec !(maxres)
772 if (nsec(i).lt.0.or.nsec(j).lt.0) return
774 if (nsec(i).gt.1.or.nsec(j).gt.1) return
777 if (isec(i,k).eq.isec(j,l)) return
783 !-----------------------------------------------------------------------------
785 !-----------------------------------------------------------------------------
786 logical function seq_comp(itypea,itypeb,length)
789 integer :: length,itypea(length),itypeb(length)
792 if (itypea(i).ne.itypeb(i)) then
799 end function seq_comp
801 !-----------------------------------------------------------------------------
803 !-----------------------------------------------------------------------------
804 subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
806 ! implicit real*8 (a-h,o-z)
807 ! include 'DIMENSIONS'
808 ! include 'COMMON.CHAIN'
809 ! include 'COMMON.CONTACTS'
810 ! include 'COMMON.IOUNITS'
811 real(kind=8) :: przes(3),obr(3,3)
812 logical :: non_conv,lprn
813 real(kind=8) :: rms,frac,frac_nn,co
814 ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
818 !elte(iout,*) "rms_nacc before contact"
819 call contact(.false.,ncont,icont,co)
820 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
821 frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
822 if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
823 'RMS deviation from the reference structure:',rms,&
824 ' % of native contacts:',frac*100,&
825 ' % of nonnative contacts:',frac_nn*100,&
829 end subroutine rms_nac_nnc
830 !-----------------------------------------------------------------------------
831 subroutine rmsd(drms)
833 use regularize_, only:fitsq
834 ! implicit real*8 (a-h,o-z)
835 ! include 'DIMENSIONS'
839 ! include 'COMMON.CHAIN'
840 ! include 'COMMON.IOUNITS'
841 ! include 'COMMON.INTERACT'
842 ! include 'COMMON.CONTROL'
844 real(kind=8) :: przes(3),obrot(3,3)
845 real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
848 real(kind=8) :: drms,rminroz,roznica
849 integer :: i,j,iatom,kkk,iti,k
851 !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
859 ! print *,"nz_start",nz_start," nz_end",nz_end
860 ! if (symetr.le.1) then
862 ! do i=nz_start,nz_end
866 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
867 ! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
869 ! if (iz_sc.eq.1.and.iti.ne.10) then
872 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
873 ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
884 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
885 crefcopy(k,iatom)=cref(k,i,kkk)
887 if (iz_sc.eq.1.and.iti.ne.10) then
890 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
891 crefcopy(k,iatom)=cref(k,nres+i,kkk)
900 ! write (iout,*) 'Ccopy and CREFcopy'
901 ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
902 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
903 ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
904 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
906 ! ----- end diagnostics
908 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
909 przes,obrot,non_conv)
911 print *,'Problems in FITSQ!!! rmsd'
912 write (iout,*) 'Problems in FITSQ!!! rmsd'
913 print *,'Ccopy and CREFcopy'
914 write (iout,*) 'Ccopy and CREFcopy'
915 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
916 (crefcopy(j,k),j=1,3),k=1,iatom)
917 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
918 (crefcopy(j,k),j=1,3),k=1,iatom)
920 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
926 ! write (iout,*) "roznica", roznica,kkk
927 if (roznica.le.rminroz) rminroz=roznica
929 drms=dsqrt(dabs(rminroz))
931 ! write (iout,*) "nperm,symetr", nperm,symetr
932 ! ---- end diagnostics
935 !-----------------------------------------------------------------------------
936 subroutine rmsd_csa(drms)
938 use regularize_, only:fitsq
939 ! implicit real*8 (a-h,o-z)
940 ! include 'DIMENSIONS'
944 ! include 'COMMON.CHAIN'
945 ! include 'COMMON.IOUNITS'
946 ! include 'COMMON.INTERACT'
948 real(kind=8) :: przes(3),obrot(3,3)
949 real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
950 integer :: kkk,iatom,ierror,ierrcode
954 real(kind=8) :: drms,roznica
956 allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
964 ccopy(k,iatom)=c(k,i)
965 crefcopy(k,iatom)=crefjlee(k,i)
967 if (iz_sc.eq.1.and.iti.ne.10) then
970 ccopy(k,iatom)=c(k,nres+i)
971 crefcopy(k,iatom)=crefjlee(k,nres+i)
976 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
977 przes,obrot,non_conv)
979 print *,'Problems in FITSQ!!! rmsd_csa'
980 write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
981 print *,'Ccopy and CREFcopy'
982 write (iout,*) 'Ccopy and CREFcopy'
983 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
984 (crefcopy(j,k),j=1,3),k=1,iatom)
985 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
986 (crefcopy(j,k),j=1,3),k=1,iatom)
988 call mpi_abort(mpi_comm_world,ierror,ierrcode)
993 drms=dsqrt(dabs(roznica))
995 end subroutine rmsd_csa
996 !-----------------------------------------------------------------------------
998 !-----------------------------------------------------------------------------
1002 use geometry, only:pinorm
1003 use random, only:ran_number,iran_num
1004 ! implicit real*8 (a-h,o-z)
1005 ! include 'DIMENSIONS'
1007 ! include 'COMMON.GEO'
1008 ! include 'COMMON.VAR'
1009 ! include 'COMMON.INTERACT'
1010 ! include 'COMMON.IOUNITS'
1011 ! include 'COMMON.DISTFIT'
1012 ! include 'COMMON.SBRIDGE'
1013 ! include 'COMMON.CONTROL'
1014 ! include 'COMMON.FFIELD'
1015 ! include 'COMMON.MINIM'
1016 ! include 'COMMON.CHAIN'
1017 real(kind=8) :: time0,time1
1018 real(kind=8) :: energy(0:n_ene),ee
1019 real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres)
1020 integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1021 logical :: debug,accepted
1022 real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1025 !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1027 call geom_to_var(nvar,var1)
1032 write(iout,*) 'etot=',0,etot,rms
1033 call secondary2(.false.)
1035 call write_pdb(0,'first structure',etot)
1044 betbol=1.0D0/(1.9858D-3*temp)
1046 d=ran_number(-pi,pi)
1047 ! phi(jr)=pinorm(phi(jr)+d)
1052 write(iout,*) 'etot=',1,etot0,rms
1053 call write_pdb(1,'perturb structure',etot0)
1057 d=ran_number(-da,da)
1059 phi(jr)=pinorm(phi(jr)+d)
1064 if (etot.lt.etot0) then
1068 xxr=ran_number(0.0D0,1.0D0)
1069 xxh=betbol*(etot-etot0)
1070 if (xxh.lt.50.0D0) then
1072 if (xxh.gt.xxr) accepted=.true.
1076 ! print *,etot0,etot,accepted
1080 write(iout,*) 'etot=',i,etot,rms
1081 call write_pdb(i,'MC structure',etot)
1083 ! call geom_to_var(nvar,var1)
1084 call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1085 call geom_to_var(nvar,var)
1086 call minimize(etot,var,iretcode,nfun)
1087 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1088 call var_to_geom(nvar,var)
1091 write(iout,*) 'etot mcm=',i,etot,rms
1092 call write_pdb(i+1,'MCM structure',etot)
1093 call var_to_geom(nvar,var1)
1101 ! call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1102 ! call geom_to_var(nvar,var)
1105 ! call write_pdb(998 ,'sc min',etot)
1107 ! call minimize(etot,var,iretcode,nfun)
1108 ! write(iout,*)'------------------------------------------------'
1109 ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1111 ! call var_to_geom(nvar,var)
1113 ! call write_pdb(999,'full min',etot)
1117 !-----------------------------------------------------------------------------
1122 ! implicit real*8 (a-h,o-z)
1123 ! include 'DIMENSIONS'
1125 ! include 'COMMON.GEO'
1126 ! include 'COMMON.VAR'
1127 ! include 'COMMON.INTERACT'
1128 ! include 'COMMON.IOUNITS'
1129 ! include 'COMMON.DISTFIT'
1130 ! include 'COMMON.SBRIDGE'
1131 ! include 'COMMON.CONTROL'
1132 ! include 'COMMON.FFIELD'
1133 ! include 'COMMON.MINIM'
1134 ! include 'COMMON.CHAIN'
1135 real(kind=8) :: time0,time1
1136 real(kind=8) :: energy(0:n_ene),ee
1137 real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1141 integer :: i,ij,ieval,iretcode,nfun
1142 real(kind=8) :: etot
1144 allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1146 call geom_to_var(nvar,var1)
1150 write(iout,*) nnt,nct,etot
1151 call write_pdb(1,'first structure',etot)
1152 call secondary2(.true.)
1161 call var_to_geom(nvar,var1)
1162 write(iout,*) 'N16 test',(jdata(i),i=1,5)
1163 call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1165 call geom_to_var(nvar,var)
1171 call minimize(etot,var,iretcode,nfun)
1172 write(iout,*)'------------------------------------------------'
1173 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1179 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1180 nfun/(time1-time0),' eval/s'
1182 call var_to_geom(nvar,var)
1184 call write_pdb(ij*100+99,'full min',etot)
1191 end subroutine test_n16
1193 !-----------------------------------------------------------------------------
1194 subroutine test_local
1196 ! implicit real*8 (a-h,o-z)
1197 ! include 'DIMENSIONS'
1198 ! include 'COMMON.GEO'
1199 ! include 'COMMON.VAR'
1200 ! include 'COMMON.INTERACT'
1201 ! include 'COMMON.IOUNITS'
1202 real(kind=8) :: time0,time1
1203 real(kind=8) :: energy(0:n_ene),ee
1204 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1206 real(kind=8) :: etot
1208 ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres)
1210 ! call geom_to_var(nvar,varia)
1211 call write_pdb(1,'first structure',0d0)
1215 write(iout,*) nnt,nct,etot
1217 write(iout,*) 'calling sc_move'
1218 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1219 write(iout,*) nft_sc,etot
1220 call write_pdb(2,'second structure',etot)
1222 write(iout,*) 'calling local_move'
1223 call local_move_init(.false.)
1224 call local_move(24,29,20d0,50d0)
1226 call write_pdb(3,'third structure',etot)
1228 write(iout,*) 'calling sc_move'
1229 call sc_move(24,29,5,10d0,nft_sc,etot)
1230 write(iout,*) nft_sc,etot
1231 call write_pdb(2,'last structure',etot)
1234 end subroutine test_local
1235 !-----------------------------------------------------------------------------
1238 ! implicit real*8 (a-h,o-z)
1239 ! include 'DIMENSIONS'
1240 ! include 'COMMON.GEO'
1241 ! include 'COMMON.VAR'
1242 ! include 'COMMON.INTERACT'
1243 ! include 'COMMON.IOUNITS'
1244 real(kind=8) :: time0,time1,etot
1245 real(kind=8) :: energy(0:n_ene),ee
1246 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1250 ! call geom_to_var(nvar,varia)
1251 call write_pdb(1,'first structure',0d0)
1255 write(iout,*) nnt,nct,etot
1257 write(iout,*) 'calling sc_move'
1259 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1260 write(iout,*) nft_sc,etot
1261 call write_pdb(2,'second structure',etot)
1263 write(iout,*) 'calling sc_move 2nd time'
1265 call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1266 write(iout,*) nft_sc,etot
1267 call write_pdb(3,'last structure',etot)
1269 end subroutine test_sc
1270 !-----------------------------------------------------------------------------
1271 subroutine bgrow(bstrand,nbstrand,in,ind,new)
1273 ! implicit real*8 (a-h,o-z)
1274 ! include 'DIMENSIONS'
1275 ! include 'COMMON.CHAIN'
1276 integer,dimension(nres/3,6) :: bstrand !(maxres/3,6)
1279 integer :: nbstrand,in,ind,new,ishift,i
1281 ishift=iabs(bstrand(in,ind+4)-new)
1283 print *,'bgrow',bstrand(in,ind+4),new,ishift
1288 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1290 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1291 if (bstrand(i,5).lt.bstrand(i,6)) then
1292 bstrand(i,5)=bstrand(i,5)-ishift
1294 bstrand(i,5)=bstrand(i,5)+ishift
1299 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1301 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1302 if (bstrand(i,6).lt.bstrand(i,5)) then
1303 bstrand(i,6)=bstrand(i,6)-ishift
1305 bstrand(i,6)=bstrand(i,6)+ishift
1312 end subroutine bgrow
1313 !-----------------------------------------------------------------------------
1316 use geometry, only:dist
1318 ! implicit real*8 (a-h,o-z)
1319 ! include 'DIMENSIONS'
1321 ! include 'COMMON.GEO'
1322 ! include 'COMMON.CHAIN'
1323 ! include 'COMMON.IOUNITS'
1324 ! include 'COMMON.VAR'
1325 ! include 'COMMON.CONTROL'
1326 ! include 'COMMON.SBRIDGE'
1327 ! include 'COMMON.FFIELD'
1328 ! include 'COMMON.MINIM'
1330 ! include 'COMMON.DISTFIT'
1331 integer :: if(20,nres),nif,ifa(20)
1332 integer :: ibc(0:nres,0:nres),istrand(20)
1333 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1334 integer :: itmp(20,nres)
1335 real(kind=8) :: time0,time1
1336 real(kind=8) :: energy(0:n_ene),ee
1337 real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres)
1339 logical :: debug,ltest,usedbfrag(nres/3)
1340 character(len=50) :: linia
1342 integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1343 integer :: bstrand(nres/3,6),nbstrand
1344 real(kind=8) :: etot
1345 integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1346 in,ind,ifun,nfun,iretcode
1347 !------------------------
1350 !------------------------
1357 call geom_to_var(nvar,vorg)
1358 call secondary2(debug)
1360 if (nbfrag.le.1) return
1363 usedbfrag(i)=.false.
1367 nbetasheet=nbetasheet+1
1369 bstrand(1,1)=bfrag(1,1)
1370 bstrand(1,2)=bfrag(2,1)
1371 bstrand(1,3)=nbetasheet
1373 bstrand(1,5)=bfrag(1,1)
1374 bstrand(1,6)=bfrag(2,1)
1375 do i=bfrag(1,1),bfrag(2,1)
1376 betasheet(i)=nbetasheet
1380 bstrand(2,1)=bfrag(3,1)
1381 bstrand(2,2)=bfrag(4,1)
1382 bstrand(2,3)=nbetasheet
1383 bstrand(2,5)=bfrag(3,1)
1384 bstrand(2,6)=bfrag(4,1)
1386 if (bfrag(3,1).le.bfrag(4,1)) then
1388 do i=bfrag(3,1),bfrag(4,1)
1389 betasheet(i)=nbetasheet
1394 do i=bfrag(4,1),bfrag(3,1)
1395 betasheet(i)=nbetasheet
1402 do while (iused_nbfrag.ne.nbfrag)
1406 IF (.not.usedbfrag(j)) THEN
1408 write (*,*) j,(bfrag(i,j),i=1,4)
1410 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1412 write (*,*) '------------------'
1415 if (bfrag(3,j).le.bfrag(4,j)) then
1416 do i=bfrag(3,j),bfrag(4,j)
1417 if(betasheet(i).eq.nbetasheet) then
1419 do k=bfrag(3,j),bfrag(4,j)
1420 betasheet(k)=nbetasheet
1425 iused_nbfrag=iused_nbfrag+1
1426 do k=bfrag(1,j),bfrag(2,j)
1427 betasheet(k)=nbetasheet
1428 ibetasheet(k)=nbstrand
1430 if (bstrand(in,4).lt.0) then
1431 bstrand(nbstrand,1)=bfrag(2,j)
1432 bstrand(nbstrand,2)=bfrag(1,j)
1433 bstrand(nbstrand,3)=nbetasheet
1434 bstrand(nbstrand,4)=-nbstrand
1435 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1436 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1437 if(bstrand(in,1).lt.bfrag(4,j)) then
1438 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1440 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1441 (bstrand(in,5)-bfrag(4,j))
1443 if(bstrand(in,2).gt.bfrag(3,j)) then
1444 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1446 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1447 (-bstrand(in,6)+bfrag(3,j))
1450 bstrand(nbstrand,1)=bfrag(1,j)
1451 bstrand(nbstrand,2)=bfrag(2,j)
1452 bstrand(nbstrand,3)=nbetasheet
1453 bstrand(nbstrand,4)=nbstrand
1454 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1455 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1456 if(bstrand(in,1).gt.bfrag(3,j)) then
1457 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1459 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1460 (-bstrand(in,5)+bfrag(3,j))
1462 if(bstrand(in,2).lt.bfrag(4,j)) then
1463 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1465 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1466 (bstrand(in,6)-bfrag(4,j))
1471 if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1472 in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1473 do k=bfrag(1,j),bfrag(2,j)
1474 betasheet(k)=nbetasheet
1479 iused_nbfrag=iused_nbfrag+1
1480 do k=bfrag(3,1),bfrag(4,1)
1481 betasheet(k)=nbetasheet
1482 ibetasheet(k)=nbstrand
1484 if (bstrand(in,4).lt.0) then
1485 bstrand(nbstrand,1)=bfrag(4,j)
1486 bstrand(nbstrand,2)=bfrag(3,j)
1487 bstrand(nbstrand,3)=nbetasheet
1488 bstrand(nbstrand,4)=-nbstrand
1489 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1490 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1491 if(bstrand(in,1).lt.bfrag(2,j)) then
1492 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1494 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1495 (bstrand(in,5)-bfrag(2,j))
1497 if(bstrand(in,2).gt.bfrag(1,j)) then
1498 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1500 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1501 (-bstrand(in,6)+bfrag(1,j))
1504 bstrand(nbstrand,1)=bfrag(3,j)
1505 bstrand(nbstrand,2)=bfrag(4,j)
1506 bstrand(nbstrand,3)=nbetasheet
1507 bstrand(nbstrand,4)=nbstrand
1508 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1509 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1510 if(bstrand(in,1).gt.bfrag(1,j)) then
1511 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1513 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1514 (-bstrand(in,5)+bfrag(1,j))
1516 if(bstrand(in,2).lt.bfrag(2,j)) then
1517 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1519 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1520 (bstrand(in,6)-bfrag(2,j))
1527 do i=bfrag(4,j),bfrag(3,j)
1528 if(betasheet(i).eq.nbetasheet) then
1530 do k=bfrag(4,j),bfrag(3,j)
1531 betasheet(k)=nbetasheet
1536 iused_nbfrag=iused_nbfrag+1
1537 do k=bfrag(1,j),bfrag(2,j)
1538 betasheet(k)=nbetasheet
1539 ibetasheet(k)=nbstrand
1541 if (bstrand(in,4).lt.0) then
1542 bstrand(nbstrand,1)=bfrag(1,j)
1543 bstrand(nbstrand,2)=bfrag(2,j)
1544 bstrand(nbstrand,3)=nbetasheet
1545 bstrand(nbstrand,4)=nbstrand
1546 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1547 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1548 if(bstrand(in,1).lt.bfrag(3,j)) then
1549 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1551 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1552 (bstrand(in,5)-bfrag(3,j))
1554 if(bstrand(in,2).gt.bfrag(4,j)) then
1555 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1557 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1558 (-bstrand(in,6)+bfrag(4,j))
1561 bstrand(nbstrand,1)=bfrag(2,j)
1562 bstrand(nbstrand,2)=bfrag(1,j)
1563 bstrand(nbstrand,3)=nbetasheet
1564 bstrand(nbstrand,4)=-nbstrand
1565 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1566 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1567 if(bstrand(in,1).gt.bfrag(4,j)) then
1568 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1570 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1571 (-bstrand(in,5)+bfrag(4,j))
1573 if(bstrand(in,2).lt.bfrag(3,j)) then
1574 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1576 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1577 (bstrand(in,6)-bfrag(3,j))
1582 if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1583 in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1584 do k=bfrag(1,j),bfrag(2,j)
1585 betasheet(k)=nbetasheet
1590 iused_nbfrag=iused_nbfrag+1
1591 do k=bfrag(4,j),bfrag(3,j)
1592 betasheet(k)=nbetasheet
1593 ibetasheet(k)=nbstrand
1595 if (bstrand(in,4).lt.0) then
1596 bstrand(nbstrand,1)=bfrag(4,j)
1597 bstrand(nbstrand,2)=bfrag(3,j)
1598 bstrand(nbstrand,3)=nbetasheet
1599 bstrand(nbstrand,4)=nbstrand
1600 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1601 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1602 if(bstrand(in,1).lt.bfrag(2,j)) then
1603 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1605 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1606 (bstrand(in,5)-bfrag(2,j))
1608 if(bstrand(in,2).gt.bfrag(1,j)) then
1609 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1611 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1612 (-bstrand(in,6)+bfrag(1,j))
1615 bstrand(nbstrand,1)=bfrag(3,j)
1616 bstrand(nbstrand,2)=bfrag(4,j)
1617 bstrand(nbstrand,3)=nbetasheet
1618 bstrand(nbstrand,4)=-nbstrand
1619 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1620 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1621 if(bstrand(in,1).gt.bfrag(1,j)) then
1622 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1624 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1625 (-bstrand(in,5)+bfrag(1,j))
1627 if(bstrand(in,2).lt.bfrag(2,j)) then
1628 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1630 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1631 (bstrand(in,6)-bfrag(2,j))
1645 do while (usedbfrag(j))
1650 nbetasheet=nbetasheet+1
1651 bstrand(nbstrand,1)=bfrag(1,j)
1652 bstrand(nbstrand,2)=bfrag(2,j)
1653 bstrand(nbstrand,3)=nbetasheet
1654 bstrand(nbstrand,5)=bfrag(1,j)
1655 bstrand(nbstrand,6)=bfrag(2,j)
1657 bstrand(nbstrand,4)=nbstrand
1658 do i=bfrag(1,j),bfrag(2,j)
1659 betasheet(i)=nbetasheet
1660 ibetasheet(i)=nbstrand
1664 bstrand(nbstrand,1)=bfrag(3,j)
1665 bstrand(nbstrand,2)=bfrag(4,j)
1666 bstrand(nbstrand,3)=nbetasheet
1667 bstrand(nbstrand,5)=bfrag(3,j)
1668 bstrand(nbstrand,6)=bfrag(4,j)
1670 if (bfrag(3,j).le.bfrag(4,j)) then
1671 bstrand(nbstrand,4)=nbstrand
1672 do i=bfrag(3,j),bfrag(4,j)
1673 betasheet(i)=nbetasheet
1674 ibetasheet(i)=nbstrand
1677 bstrand(nbstrand,4)=-nbstrand
1678 do i=bfrag(4,j),bfrag(3,j)
1679 betasheet(i)=nbetasheet
1680 ibetasheet(i)=nbstrand
1684 iused_nbfrag=iused_nbfrag+1
1690 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1697 if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1701 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1704 !------------------------
1708 if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1709 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1711 ifb(nifb,1)=bstrand(i,4)
1712 ifb(nifb,2)=bstrand(j,4)
1719 write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1725 write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1727 nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1729 if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1730 nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1736 if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1737 if (if(j,i).gt.0) then
1738 if(betasheet(if(j,i)).eq.0 .or. &
1739 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1744 write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1747 ! read (inp,*) (ifa(i),i=1,4)
1749 ! read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1753 !------------------------
1758 !ccccccccccccccccccccccccccccccccc
1760 !ccccccccccccccccccccccccccccccccc
1764 istrand(is-j+1)=int(ii/is**(is-j))
1765 ii=ii-istrand(is-j+1)*is**(is-j)
1769 istrand(k)=istrand(k)+1
1770 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1774 if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1775 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1784 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1785 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1786 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1787 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1793 if (mod(isa,2).eq.0) then
1795 if (istrand(k).eq.1) ltest=.false.
1798 do k=(isa+1)/2+1,isa
1799 if (istrand(k).eq.1) ltest=.false.
1803 IF (ltest.and.lifb0.eq.1) THEN
1806 call var_to_geom(nvar,vorg)
1808 write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1809 write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1810 write (linia,'(10i3)') (istrand(k),k=1,isa)
1820 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1822 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1826 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1832 write(*,*) (itmp(j,i),j=1,4)
1836 ! ifa(1),ifa(2),ifa(3),ifa(4)
1837 ! if(1,i),if(2,i),if(3,i),if(4,i)
1842 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1843 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1844 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1845 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1853 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1855 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1859 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1861 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1864 !------------------------
1867 ! freeze sec.elements
1877 do i=bfrag(1,j),bfrag(2,j)
1882 if (bfrag(3,j).le.bfrag(4,j)) then
1883 do i=bfrag(3,j),bfrag(4,j)
1889 do i=bfrag(4,j),bfrag(3,j)
1897 do i=hfrag(1,j),hfrag(2,j)
1905 !------------------------
1906 ! generate constrains
1914 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1922 else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1930 else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1938 else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1946 else if ( ibc(i,j).gt.0 ) then
1947 d0(ind)=DIST(i,ibc(i,j))
1954 else if ( ibc(j,i).gt.0 ) then
1955 d0(ind)=DIST(ibc(j,i),j)
1969 !d--------------------------
1971 write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1972 ibc(jhpb(i),ihpb(i)),' --',&
1973 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1980 call contact_cp_min(varia,ifun,iconf,linia,debug)
1985 call minimize(etot,varia,iretcode,nfun)
1986 write(iout,*)'------------------------------------------------'
1987 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1993 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1994 nfun/(time1-time0),' eval/s'
1996 write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
1997 call var_to_geom(nvar,varia)
1999 call write_pdb(900+iconf,linia,etot)
2004 call enerprint(energy)
2006 !d call briefout(0,etot)
2007 !d call secondary2(.true.)
2011 !ccccccccccccccccccccccccccccccccccc
2014 !ccccccccccccccccccccccccccccccccccc
2017 10 write (iout,'(a)') 'Error reading test structure.'
2019 end subroutine test11
2020 !-----------------------------------------------------------------------------
2023 use geometry, only:dist
2025 ! implicit real*8 (a-h,o-z)
2026 ! include 'DIMENSIONS'
2028 ! include 'COMMON.GEO'
2029 ! include 'COMMON.CHAIN'
2030 ! include 'COMMON.IOUNITS'
2031 ! include 'COMMON.VAR'
2032 ! include 'COMMON.CONTROL'
2033 ! include 'COMMON.SBRIDGE'
2034 ! include 'COMMON.FFIELD'
2035 ! include 'COMMON.MINIM'
2037 ! include 'COMMON.DISTFIT'
2038 integer :: if(3,nres),nif
2039 integer :: ibc(nres,nres),istrand(20)
2040 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2041 real(kind=8) :: time0,time1
2042 real(kind=8) :: energy(0:n_ene),ee
2043 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2045 logical :: debug,ltest
2046 character(len=50) :: linia
2047 integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2048 real(kind=8) :: etot
2051 read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2054 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2058 !------------------------
2059 call secondary2(debug)
2060 !------------------------
2068 ! freeze sec.elements and store indexes for beta constrains
2078 do i=bfrag(1,j),bfrag(2,j)
2083 if (bfrag(3,j).le.bfrag(4,j)) then
2084 do i=bfrag(3,j),bfrag(4,j)
2088 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2091 do i=bfrag(4,j),bfrag(3,j)
2095 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2100 do i=hfrag(1,j),hfrag(2,j)
2109 ! ---------------- test --------------
2111 if (ibc(if(1,i),if(2,i)).eq.-1) then
2112 ibc(if(1,i),if(2,i))=if(3,i)
2113 ibc(if(1,i),if(3,i))=if(2,i)
2114 else if (ibc(if(2,i),if(1,i)).eq.-1) then
2115 ibc(if(2,i),if(1,i))=0
2116 ibc(if(1,i),if(2,i))=if(3,i)
2117 ibc(if(1,i),if(3,i))=if(2,i)
2119 ibc(if(1,i),if(2,i))=if(3,i)
2120 ibc(if(1,i),if(3,i))=if(2,i)
2126 if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j)
2129 !------------------------
2135 if ( ibc(i,j).eq.-1 ) then
2143 else if ( ibc(i,j).gt.0 ) then
2144 d0(ind)=DIST(i,ibc(i,j))
2151 else if ( ibc(j,i).gt.0 ) then
2152 d0(ind)=DIST(ibc(j,i),j)
2166 !d--------------------------
2167 write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2168 ibc(jhpb(i),ihpb(i)),' --',&
2169 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2177 call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2182 call minimize(etot,varia,iretcode,nfun)
2183 write(iout,*)'------------------------------------------------'
2184 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2190 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2191 nfun/(time1-time0),' eval/s'
2194 call var_to_geom(nvar,varia)
2196 call write_pdb(999,'full min',etot)
2201 call enerprint(energy)
2203 call briefout(0,etot)
2204 call secondary2(.true.)
2207 10 write (iout,'(a)') 'Error reading test structure.'
2209 end subroutine test3
2210 !-----------------------------------------------------------------------------
2214 ! implicit real*8 (a-h,o-z)
2215 ! include 'DIMENSIONS'
2217 ! include 'COMMON.GEO'
2218 ! include 'COMMON.CHAIN'
2219 ! include 'COMMON.IOUNITS'
2220 ! include 'COMMON.VAR'
2221 ! include 'COMMON.CONTROL'
2222 ! include 'COMMON.SBRIDGE'
2223 ! include 'COMMON.FFIELD'
2224 ! include 'COMMON.MINIM'
2226 ! include 'COMMON.DISTFIT'
2227 integer :: if(2,2),ind
2228 integer :: iff(nres)
2229 real(kind=8) :: time0,time1
2230 real(kind=8) :: energy(0:n_ene),ee
2231 real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2232 theta1,phi1,alph1,omeg1 !(maxres)
2233 real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres)
2235 integer :: i,j,nn,ifun,iretcode,nfun
2236 real(kind=8) :: etot
2239 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2240 write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2241 read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2242 read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2243 read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2244 read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2246 theta2(i)=deg2rad*theta2(i)
2247 phi2(i)=deg2rad*phi2(i)
2248 alph2(i)=deg2rad*alph2(i)
2249 omeg2(i)=deg2rad*omeg2(i)
2263 !------------------------
2268 do i=if(j,1),if(j,2)
2274 call geom_to_var(nvar,varia)
2275 call write_pdb(1,'first structure',0d0)
2277 call secondary(.true.)
2279 call secondary2(.true.)
2282 if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2283 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2284 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2287 if (bfrag(3,j).lt.bfrag(4,j)) then
2288 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2289 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2290 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2292 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2293 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2294 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2307 call geom_to_var(nvar,varia2)
2308 call write_pdb(2,'second structure',0d0)
2312 !-------------------------------------------------------
2315 call contact_cp(varia,varia2,iff,ifun,7)
2320 call minimize(etot,varia,iretcode,nfun)
2321 write(iout,*)'------------------------------------------------'
2322 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2328 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2329 nfun/(time1-time0),' eval/s'
2332 call var_to_geom(nvar,varia)
2334 call write_pdb(999,'full min',etot)
2339 call enerprint(energy)
2341 call briefout(0,etot)
2344 10 write (iout,'(a)') 'Error reading test structure.'
2346 end subroutine test__
2347 !-----------------------------------------------------------------------------
2348 subroutine secondary(lprint)
2350 ! implicit real*8 (a-h,o-z)
2351 ! include 'DIMENSIONS'
2352 ! include 'COMMON.CHAIN'
2353 ! include 'COMMON.IOUNITS'
2354 ! include 'COMMON.DISTFIT'
2356 integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2357 logical :: lprint,not_done
2358 real(kind=4) :: dcont(nres*nres/2),d
2359 real(kind=4) :: rcomp = 7.0
2360 real(kind=4) :: rbeta = 5.2
2361 real(kind=4) :: ralfa = 5.2
2362 real(kind=4) :: r310 = 6.6
2363 real(kind=8),dimension(3) :: xpi,xpj
2364 integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2366 call chainbuild_cart
2367 !d call write_pdb(99,'sec structure',0d0)
2379 xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2383 xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2385 !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2386 !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2387 !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2388 !d print *,'CA',i,j,d
2389 d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2390 (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2391 (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2392 if ( d.lt.rcomp*rcomp) then
2396 dcont(ncont)=sqrt(d)
2402 write (iout,'(a)') '#PP contact map distances:'
2404 write (iout,'(3i4,f10.5)') &
2405 i,icont(1,i),icont(2,i),dcont(i)
2409 ! finding parallel beta
2410 !d write (iout,*) '------- looking for parallel beta -----------'
2416 if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2417 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2418 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2419 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2420 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2421 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2425 !d write (iout,*) i1,j1,dcont(i)
2431 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2432 .and. dcont(j).le.rbeta .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) &
2442 !d write (iout,*) i1,j1,dcont(j),not_done
2446 if (i1-ii1.gt.1) then
2450 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2459 isec(ij,1)=isec(ij,1)+1
2460 isec(ij,1+isec(ij,1))=nbeta
2463 isec(ij,1)=isec(ij,1)+1
2464 isec(ij,1+isec(ij,1))=nbeta
2469 if (nbeta.le.9) then
2470 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2471 "DefPropRes 'strand",nstrand,&
2472 "' 'num = ",ii1-1,"..",i1-1,"'"
2474 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2475 "DefPropRes 'strand",nstrand,&
2476 "' 'num = ",ii1-1,"..",i1-1,"'"
2479 if (nbeta.le.9) then
2480 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2481 "DefPropRes 'strand",nstrand,&
2482 "' 'num = ",jj1-1,"..",j1-1,"'"
2484 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2485 "DefPropRes 'strand",nstrand,&
2486 "' 'num = ",jj1-1,"..",j1-1,"'"
2488 write(12,'(a8,4i4)') &
2489 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2495 ! finding antiparallel beta
2496 !d write (iout,*) '--------- looking for antiparallel beta ---------'
2501 if (dcont(i).le.rbeta.and. &
2502 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2503 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2504 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2505 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2506 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2510 !d write (iout,*) i1,j1,dcont(i)
2517 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .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) &
2523 .and. dcont(j).le.rbeta ) goto 6
2527 !d write (iout,*) i1,j1,dcont(j),not_done
2531 if (i1-ii1.gt.1) then
2532 if(lprint)write (iout,*)'antiparallel beta',&
2533 nbeta,ii1-1,i1,jj1,j1-1
2536 bfrag(1,nbfrag)=max0(ii1-1,1)
2539 bfrag(4,nbfrag)=max0(j1-1,1)
2544 isec(ij,1)=isec(ij,1)+1
2545 isec(ij,1+isec(ij,1))=nbeta
2549 isec(ij,1)=isec(ij,1)+1
2550 isec(ij,1+isec(ij,1))=nbeta
2556 if (nstrand.le.9) then
2557 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2558 "DefPropRes 'strand",nstrand,&
2559 "' 'num = ",ii1-2,"..",i1-1,"'"
2561 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2562 "DefPropRes 'strand",nstrand,&
2563 "' 'num = ",ii1-2,"..",i1-1,"'"
2566 if (nstrand.le.9) then
2567 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2568 "DefPropRes 'strand",nstrand,&
2569 "' 'num = ",j1-2,"..",jj1-1,"'"
2571 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2572 "DefPropRes 'strand",nstrand,&
2573 "' 'num = ",j1-2,"..",jj1-1,"'"
2575 write(12,'(a8,4i4)') &
2576 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2582 if (nstrand.gt.0.and.lprint) then
2583 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2586 write(12,'(a9,i1,$)') " | strand",i
2588 write(12,'(a9,i2,$)') " | strand",i
2591 write(12,'(a1)') "'"
2595 ! finding alpha or 310 helix
2601 if (j1.eq.i1+3.and.dcont(i).le.r310 &
2602 .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2603 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2604 !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2607 if (isec(ii1,1).eq.0) then
2616 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2620 !d write (iout,*) i1,j1,not_done
2623 if (j1-ii1.gt.4) then
2625 !d write (iout,*)'helix',nhelix,ii1,j1
2629 hfrag(2,nhfrag)=max0(j1-1,1)
2635 write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2636 if (nhelix.le.9) then
2637 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2638 "DefPropRes 'helix",nhelix,&
2639 "' 'num = ",ii1-1,"..",j1-2,"'"
2641 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2642 "DefPropRes 'helix",nhelix,&
2643 "' 'num = ",ii1-1,"..",j1-2,"'"
2650 if (nhelix.gt.0.and.lprint) then
2651 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2653 if (nhelix.le.9) then
2654 write(12,'(a8,i1,$)') " | helix",i
2656 write(12,'(a8,i2,$)') " | helix",i
2659 write(12,'(a1)') "'"
2663 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2664 write(12,'(a20)') "XMacStand ribbon.mac"
2668 end subroutine secondary
2669 !-----------------------------------------------------------------------------
2670 subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2672 use geometry, only:dist
2674 ! implicit real*8 (a-h,o-z)
2675 ! include 'DIMENSIONS'
2677 ! include 'COMMON.SBRIDGE'
2678 ! include 'COMMON.FFIELD'
2679 ! include 'COMMON.IOUNITS'
2680 ! include 'COMMON.DISTFIT'
2681 ! include 'COMMON.VAR'
2682 ! include 'COMMON.CHAIN'
2683 ! include 'COMMON.MINIM'
2685 character(len=50) :: linia
2687 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2688 real(kind=8) :: time0,time1
2689 integer :: iff(nres),ieval
2690 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2693 integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2694 real(kind=8) :: wstrain0,etot
2696 maxres22=nres*(nres+1)/2
2698 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2699 call var_to_geom(nvar,var)
2706 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2728 call var_to_geom(nvar,var2)
2731 if ( iff(i).eq.1 ) then
2740 !d call write_pdb(3,'combined structure',0d0)
2741 !d time0=MPI_WTIME()
2744 NY=((NRES-4)*(NRES-5))/2
2745 call distfit(.true.,200)
2747 !d time1=MPI_WTIME()
2748 !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2758 call geom_to_var(nvar,var)
2759 !d time0=MPI_WTIME()
2760 call minimize(etot,var,iretcode,nfun)
2761 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2763 !d time1=MPI_WTIME()
2764 !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2765 !d & nfun/(time1-time0),' SOFT eval/s'
2766 call var_to_geom(nvar,var)
2772 if (iff(1).eq.1) then
2778 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2783 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2789 if (iff(nres).eq.1) then
2795 !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2796 !d & "select",ij(1),"-",ij(2),
2797 !d & ",",ij(3),"-",ij(4)
2798 !d call write_pdb(in_pdb,linia,etot)
2804 !d time0=MPI_WTIME()
2805 call minimize(etot,var,iretcode,nfun)
2806 !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2809 !d time1=MPI_WTIME()
2810 !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2811 !d & nfun/(time1-time0),' eval/s'
2812 !d call var_to_geom(nvar,var)
2814 !d call write_pdb(6,'dist structure',etot)
2823 end subroutine contact_cp2
2824 !-----------------------------------------------------------------------------
2825 subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2827 use geometry, only:dist
2829 ! implicit real*8 (a-h,o-z)
2830 ! include 'DIMENSIONS'
2831 ! include 'COMMON.SBRIDGE'
2832 ! include 'COMMON.FFIELD'
2833 ! include 'COMMON.IOUNITS'
2834 ! include 'COMMON.DISTFIT'
2835 ! include 'COMMON.VAR'
2836 ! include 'COMMON.CHAIN'
2837 ! include 'COMMON.MINIM'
2839 character(len=50) :: linia
2841 real(kind=8) :: energy(0:n_ene)
2842 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2843 real(kind=8) :: time0,time1
2844 integer :: iff(nres),ieval
2845 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2849 integer :: in_pdb,i,j,ind,iwsk
2853 if (ieval.eq.-1) debug=.true.
2857 ! store selected dist. constrains from 1st structure
2860 ! Intercept NaNs in the coordinates
2861 ! write(iout,*) (var(i),i=1,nvar)
2866 if (x_sum.ne.x_sum) then
2867 write(iout,*)" *** contact_cp : Found NaN in coordinates"
2869 print *," *** contact_cp : Found NaN in coordinates"
2875 call var_to_geom(nvar,var)
2882 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2905 ! freeze sec.elements from 2nd structure
2913 call var_to_geom(nvar,var2)
2914 call secondary2(debug)
2916 do i=bfrag(1,j),bfrag(2,j)
2921 if (bfrag(3,j).le.bfrag(4,j)) then
2922 do i=bfrag(3,j),bfrag(4,j)
2928 do i=bfrag(4,j),bfrag(3,j)
2936 do i=hfrag(1,j),hfrag(2,j)
2945 ! copy selected res from 1st to 2nd structure
2949 if ( iff(i).eq.1 ) then
2959 ! prepare description in linia variable
2963 if (iff(1).eq.1) then
2969 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2974 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2980 if (iff(nres).eq.1) then
2985 write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2986 "SELECT",ij(1)-1,"-",ij(2)-1,&
2987 ",",ij(3)-1,"-",ij(4)-1
2993 call contact_cp_min(var,ieval,in_pdb,linia,debug)
2996 end subroutine contact_cp
2997 !-----------------------------------------------------------------------------
2998 subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3002 ! input : theta,phi,alph,omeg,in_pdb,linia,debug
3003 ! output : var,ieval
3005 ! implicit real*8 (a-h,o-z)
3006 ! include 'DIMENSIONS'
3008 ! include 'COMMON.SBRIDGE'
3009 ! include 'COMMON.FFIELD'
3010 ! include 'COMMON.IOUNITS'
3011 ! include 'COMMON.DISTFIT'
3012 ! include 'COMMON.VAR'
3013 ! include 'COMMON.CHAIN'
3014 ! include 'COMMON.MINIM'
3016 character(len=50) :: linia
3018 real(kind=8) :: energy(0:n_ene)
3019 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3020 real(kind=8) :: time0,time1
3021 integer :: ieval,info(3)
3022 logical :: debug,fail,reduce,change !check_var,
3025 integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3027 real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3028 wtor_d01,wstrain0,etot
3030 write(iout,'(a20,i6,a20)') &
3031 '------------------',in_pdb,'-------------------'
3035 call write_pdb(1000+in_pdb,'combined structure',0d0)
3042 ! run optimization of distances
3044 ! uses d0(),w() and mask() for frozen 2D
3046 !test---------------------------------------------
3048 !test NY=((NRES-4)*(NRES-5))/2
3049 !test call distfit(debug,5000)
3081 call geom_to_var(nvar,var)
3082 !de change=reduce(var)
3083 if (check_var(var,info)) then
3084 write(iout,*) 'cp_min error in input'
3085 print *,'cp_min error in input'
3089 !d call etotal(energy(0))
3090 !d call enerprint(energy(0))
3094 !dtest call minimize(etot,var,iretcode,nfun)
3095 !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3098 !d call etotal(energy(0))
3099 !d call enerprint(energy(0))
3118 !test--------------------------------------------------
3124 write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3125 call write_pdb(2000+in_pdb,'distfit structure',0d0)
3133 ! run soft pot. optimization
3135 ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3137 ! mask_phi(),mask_theta(),mask_side(),mask_r
3143 !de change=reduce(var)
3144 !de if (check_var(var,info)) write(iout,*) 'error before soft'
3148 call minimize(etot,var,iretcode,nfun)
3150 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3154 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3155 nfun/(time1-time0),' SOFT eval/s'
3158 call var_to_geom(nvar,var)
3160 call write_pdb(3000+in_pdb,'soft structure',etot)
3164 ! run full UNRES optimization with constrains and frozen 2D
3165 ! the same variables as soft pot. optimizatio
3171 ! check overlaps before calling full UNRES minim
3173 call var_to_geom(nvar,var)
3177 write(iout,*) 'N7 ',energy(0)
3178 if (energy(0).ne.energy(0)) then
3179 write(iout,*) 'N7 error - gives NaN',energy(0)
3183 if (energy(1).eq.1.0d20) then
3184 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3185 call overlap_sc(fail)
3189 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3201 !dte time0=MPI_WTIME()
3202 !de change=reduce(var)
3203 !de if (check_var(var,info)) then
3204 !de write(iout,*) 'error before mask dist'
3205 !de call var_to_geom(nvar,var)
3207 !de call write_pdb(10000+in_pdb,'before mask dist',etot)
3209 !dte call minimize(etot,var,iretcode,nfun)
3210 !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3211 !dte & ' eval ',nfun
3212 !dte ieval=ieval+nfun
3214 !dte time1=MPI_WTIME()
3215 !dte write (iout,'(a,f6.2,f8.2,a)')
3216 !dte & ' Time for mask dist min.',time1-time0,
3217 !dte & nfun/(time1-time0),' eval/s'
3218 !dte call flush(iout)
3221 call var_to_geom(nvar,var)
3223 call write_pdb(4000+in_pdb,'mask dist',etot)
3226 ! switch off freezing of 2D and
3227 ! run full UNRES optimization with constrains
3233 !de change=reduce(var)
3234 !de if (check_var(var,info)) then
3235 !de write(iout,*) 'error before dist'
3236 !de call var_to_geom(nvar,var)
3238 !de call write_pdb(11000+in_pdb,'before dist',etot)
3241 call minimize(etot,var,iretcode,nfun)
3243 !de change=reduce(var)
3244 !de if (check_var(var,info)) then
3245 !de write(iout,*) 'error after dist',ico
3246 !de call var_to_geom(nvar,var)
3248 !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3250 write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3256 write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3257 nfun/(time1-time0),' eval/s'
3259 !de call etotal(energy(0))
3260 !de write(iout,*) 'N7 after dist',energy(0)
3264 call var_to_geom(nvar,var)
3266 call write_pdb(in_pdb,linia,etot)
3278 end subroutine contact_cp_min
3279 !-----------------------------------------------------------------------------
3282 use geometry, only:dist
3284 ! implicit real*8 (a-h,o-z)
3285 ! include 'DIMENSIONS'
3287 ! include 'COMMON.GEO'
3288 ! include 'COMMON.CHAIN'
3289 ! include 'COMMON.IOUNITS'
3290 ! include 'COMMON.VAR'
3291 ! include 'COMMON.CONTROL'
3292 ! include 'COMMON.SBRIDGE'
3293 ! include 'COMMON.FFIELD'
3294 ! include 'COMMON.MINIM'
3295 ! include 'COMMON.INTERACT'
3297 ! include 'COMMON.DISTFIT'
3298 integer :: iff(nres)
3299 real(kind=8) :: time0,time1
3300 real(kind=8) :: energy(0:n_ene),ee
3301 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3303 logical :: debug,ltest,fail
3304 character(len=50) :: linia
3305 integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3307 real(kind=8) :: wstrain0,wang0,etot
3313 !------------------------
3315 ! freeze sec.elements
3325 do i=bfrag(1,j),bfrag(2,j)
3330 if (bfrag(3,j).le.bfrag(4,j)) then
3331 do i=bfrag(3,j),bfrag(4,j)
3337 do i=bfrag(4,j),bfrag(3,j)
3345 do i=hfrag(1,j),hfrag(2,j)
3357 ! store dist. constrains
3361 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3366 dhpb(nhpb)=DIST(i,j)
3374 call write_pdb(100+in_pdb,'input reg. structure',0d0)
3384 ! run soft pot. optimization
3390 call geom_to_var(nvar,var)
3394 call minimize(etot,var,iretcode,nfun)
3396 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3400 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3401 nfun/(time1-time0),' SOFT eval/s'
3403 call var_to_geom(nvar,var)
3405 call write_pdb(300+in_pdb,'soft structure',etot)
3408 ! run full UNRES optimization with constrains and frozen 2D
3409 ! the same variables as soft pot. optimizatio
3418 call minimize(etot,var,iretcode,nfun)
3419 write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3426 write (iout,'(a,f6.2,f8.2,a)') &
3427 ' Time for mask dist min.',time1-time0,&
3428 nfun/(time1-time0),' eval/s'
3430 call var_to_geom(nvar,var)
3432 call write_pdb(400+in_pdb,'mask & dist',etot)
3435 ! switch off constrains and
3436 ! run full UNRES optimization with frozen 2D
3451 call minimize(etot,var,iretcode,nfun)
3452 write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3458 write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,&
3459 nfun/(time1-time0),' eval/s'
3463 call var_to_geom(nvar,var)
3465 call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3472 ! run full UNRES optimization with constrains and NO frozen 2D
3482 wstrain=wstrain0/ico
3487 call minimize(etot,var,iretcode,nfun)
3488 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3489 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3496 write (iout,'(a,f6.2,f8.2,a)') &
3497 ' Time for dist min.',time1-time0,&
3498 nfun/(time1-time0),' eval/s'
3500 call var_to_geom(nvar,var)
3502 call write_pdb(600+in_pdb+ico,'dist cons',etot)
3519 call minimize(etot,var,iretcode,nfun)
3520 write(iout,*)'------------------------------------------------'
3521 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3527 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
3528 nfun/(time1-time0),' eval/s'
3531 call var_to_geom(nvar,var)
3533 call write_pdb(999,'full min',etot)
3537 end subroutine softreg
3538 !-----------------------------------------------------------------------------
3539 subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3541 use geometry, only:dist
3543 ! implicit real*8 (a-h,o-z)
3544 ! include 'DIMENSIONS'
3546 ! include 'COMMON.GEO'
3547 ! include 'COMMON.VAR'
3548 ! include 'COMMON.INTERACT'
3549 ! include 'COMMON.IOUNITS'
3550 ! include 'COMMON.DISTFIT'
3551 ! include 'COMMON.SBRIDGE'
3552 ! include 'COMMON.CONTROL'
3553 ! include 'COMMON.FFIELD'
3554 ! include 'COMMON.MINIM'
3555 ! include 'COMMON.CHAIN'
3556 real(kind=8) :: time0,time1
3557 real(kind=8) :: energy(0:n_ene),ee
3558 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3559 integer :: jdata(5),isec(nres)
3562 integer :: i1,i2,i3,i4,i5,ieval,ij
3563 integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3564 real(kind=8) :: etot,wscloc0,wstrain0
3572 call secondary2(.false.)
3578 do i=bfrag(1,j),bfrag(2,j)
3581 do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3586 do i=hfrag(1,j),hfrag(2,j)
3592 ! cut strands at the ends
3594 if (jdata(2)-jdata(1).gt.3) then
3597 if (jdata(3).lt.jdata(4)) then
3607 !v call etotal(energy(0))
3609 !v write(iout,*) nnt,nct,etot
3610 !v call write_pdb(ij*100,'first structure',etot)
3611 !v write(iout,*) 'N16 test',(jdata(i),i=1,5)
3613 !------------------------
3614 ! generate constrains
3617 if(ishift.eq.0) ishift=-2
3620 do i=jdata(1),jdata(2)
3622 if(jdata(4).gt.jdata(3))then
3623 do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3625 !d print *,i,j,j+ishift
3630 dhpb(nhpb)=DIST(i,j+ishift)
3633 do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3635 !d print *,i,j,j+ishift
3640 dhpb(nhpb)=DIST(i,j+ishift)
3647 if(isec(i).gt.0.or.isec(j).gt.0) then
3653 dhpb(nhpb)=DIST(i,j)
3660 call geom_to_var(nvar,var)
3667 wstrain=wstrain0/ico
3669 !v time0=MPI_WTIME()
3670 call minimize(etot,var,iretcode,nfun)
3671 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3672 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3675 !v time1=MPI_WTIME()
3676 !v write (iout,'(a,f6.2,f8.2,a)')
3677 !v & ' Time for dist min.',time1-time0,
3678 !v & nfun/(time1-time0),' eval/s'
3679 !v call var_to_geom(nvar,var)
3681 !v call write_pdb(ij*100+ico,'dist cons',etot)
3693 call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3696 !v call etotal(energy(0))
3698 !v call write_pdb(ij*100+10,'sc_move',etot)
3700 !d print *,nft_sc,etot
3703 end subroutine beta_slide
3704 !-----------------------------------------------------------------------------
3705 subroutine beta_zip(i1,i2,ieval,ij)
3708 ! implicit real*8 (a-h,o-z)
3709 ! include 'DIMENSIONS'
3711 ! include 'COMMON.GEO'
3712 ! include 'COMMON.VAR'
3713 ! include 'COMMON.INTERACT'
3714 ! include 'COMMON.IOUNITS'
3715 ! include 'COMMON.DISTFIT'
3716 ! include 'COMMON.SBRIDGE'
3717 ! include 'COMMON.CONTROL'
3718 ! include 'COMMON.FFIELD'
3719 ! include 'COMMON.MINIM'
3720 ! include 'COMMON.CHAIN'
3721 real(kind=8) :: time0,time1
3722 real(kind=8) :: energy(0:n_ene),ee
3723 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3724 character(len=10) :: test
3726 integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3727 real(kind=8) :: etot,wstrain0
3729 !v call etotal(energy(0))
3731 !v write(test,'(2i5)') i1,i2
3732 !v call write_pdb(ij*100,test,etot)
3733 !v write(iout,*) 'N17 test',i1,i2,etot,ij
3736 ! generate constrains
3747 call geom_to_var(nvar,var)
3753 wstrain=wstrain0/ico
3754 !v time0=MPI_WTIME()
3755 call minimize(etot,var,iretcode,nfun)
3756 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3757 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3760 !v time1=MPI_WTIME()
3761 !v write (iout,'(a,f6.2,f8.2,a)')
3762 !v & ' Time for dist min.',time1-time0,
3763 !v & nfun/(time1-time0),' eval/s'
3764 ! do not comment the next line
3765 call var_to_geom(nvar,var)
3767 !v call write_pdb(ij*100+ico,'dist cons',etot)
3775 !v call etotal(energy(0))
3777 !v write(iout,*) 'N17 test end',i1,i2,etot,ij
3780 end subroutine beta_zip
3781 !-----------------------------------------------------------------------------
3783 !-----------------------------------------------------------------------------
3784 subroutine thread_seq
3786 use geometry, only:dist
3787 use random, only:iran_num
3788 use control, only:tcpu
3789 use regularize_, only:regularize
3790 use mcm_data, only: nsave_part,nacc_tot
3791 ! Thread the sequence through a database of known structures
3792 ! implicit real*8 (a-h,o-z)
3793 ! include 'DIMENSIONS'
3794 use MPI_data !include 'COMMON.INFO'
3799 ! include 'COMMON.CONTROL'
3800 ! include 'COMMON.CHAIN'
3801 ! include 'COMMON.DBASE'
3802 ! include 'COMMON.INTERACT'
3803 ! include 'COMMON.VAR'
3804 ! include 'COMMON.THREAD'
3805 ! include 'COMMON.FFIELD'
3806 ! include 'COMMON.SBRIDGE'
3807 ! include 'COMMON.HEADER'
3808 ! include 'COMMON.IOUNITS'
3809 ! include 'COMMON.TIME1'
3810 ! include 'COMMON.CONTACTS'
3811 ! include 'COMMON.MCM'
3812 ! include 'COMMON.NAMES'
3814 integer :: ThreadId,ThreadType,Kwita
3816 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3817 real(kind=8) :: przes(3),obr(3,3)
3818 real(kind=8) :: time_for_thread
3819 logical :: found_pattern,non_conv
3820 character(len=32) :: head_pdb
3821 real(kind=8) :: energia(0:n_ene)
3822 integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3824 real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3826 n_ene_comp=nprint_ene
3831 if (me.eq.king) then
3850 ave_time_for_thread=0.0D0
3851 max_time_for_thread=0.0D0
3852 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3853 nthread=nexcl+nthread
3854 do ithread=1,nthread
3855 found_pattern=.false.
3857 do while (.not.found_pattern)
3859 if (itrial.gt.1000) then
3860 write (iout,'(/a/)') 'Too many attempts to find pattern.'
3863 call recv_stop_sig(Kwita)
3864 call send_stop_sig(-3)
3868 ! Find long enough chain in the database
3870 nres_t=nres_base(1,ii)
3871 ! Select the starting position to thread.
3872 print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3873 nres_t,' nres0=',nres0
3874 if (nres_t.ge.nres0) then
3875 ist=iran_num(0,nres_t-nres0)
3877 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3878 if (Kwita.lt.0) then
3879 write (iout,*) 'Stop signal received. Terminating.'
3880 write (*,*) 'Stop signal received. Terminating.'
3882 write (*,*) 'ithread=',ithread,' nthread=',nthread
3885 call pattern_receive
3888 if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3890 found_pattern=.true.
3892 ! If this point is reached, the pattern has not yet been examined.
3894 ! print *,'found_pattern:',found_pattern
3900 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3901 if (Kwita.lt.0) then
3902 write (iout,*) 'Stop signal received. Terminating.'
3904 write (*,*) 'ithread=',ithread,' nthread=',nthread
3910 ipatt(2,ithread)=ist
3912 write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3913 'Processor:',me,' Attempt:',ithread,&
3914 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3915 ' start at res.',ist+1
3916 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3917 ' Attempt:',ithread,&
3918 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3919 ' start at res.',ist+1
3921 write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3922 'Attempt:',ithread,&
3923 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3924 ' start at res.',ist+1
3925 write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3926 'Attempt:',ithread,&
3927 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3928 ' start at res.',ist+1
3931 ! Copy coordinates from the database.
3935 c(j,i)=cart_base(j,i+ist,ii)
3938 !d write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3)
3940 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3942 !d write (iout,'(a,f10.5)')
3943 !d & 'Initial RMS deviation from reference structure:',rms
3944 if (itype(nres).eq.ntyp1) then
3946 dcj=c(j,nres-2)-c(j,nres-3)
3947 c(j,nres)=c(j,nres-1)+dcj
3948 c(j,2*nres)=c(j,nres)
3951 if (itype(1).eq.ntyp1) then
3958 call int_from_cart(.false.,.false.)
3959 !d print *,'Exit INT_FROM_CART.'
3960 !d print *,'nhpb=',nhpb
3965 ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3967 ! stop 'End generate'
3968 ! Generate SC conformations.
3972 !d print *,'Processor:',me,': exit GEN_SIDE.'
3974 !d print *,'Exit GEN_SIDE.'
3976 ! Calculate initial energy.
3978 call etotal(energia)
3981 ener0(i,ithread)=energia(i)
3983 ener0(n_ene_comp+1,ithread)=energia(0)
3985 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
3986 ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
3988 ener0(n_ene_comp+2,ithread)=rms
3989 ener0(n_ene_comp+4,ithread)=frac
3990 ener0(n_ene_comp+5,ithread)=frac_nn
3992 ener0(n_ene_comp+3,ithread)=0.0d0
3995 print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
3997 print*,'ithread=',ithread,' Start REGULARIZE.'
4000 call regularize(nct-nnt+1,etot,rms,&
4001 cart_base(1,ist+nnt,ipattern),iretcode)
4003 time_for_thread=curr_tim1-curr_tim
4004 ave_time_for_thread= &
4005 ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4006 if (time_for_thread.gt.max_time_for_thread) &
4007 max_time_for_thread=time_for_thread
4009 print *,'Processor',me,': Exit REGULARIZE.'
4010 if (WhatsUp.eq.2) then
4012 'Sufficient number of confs. collected. Terminating.'
4015 else if (WhatsUp.eq.-1) then
4017 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4018 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4019 call send_stop_sig(-2)
4021 else if (WhatsUp.eq.-2) then
4023 write (iout,*) 'Timeup signal received. Terminating.'
4025 else if (WhatsUp.eq.-3) then
4027 write (iout,*) 'Error stop signal received. Terminating.'
4031 print *,'Exit REGULARIZE.'
4032 if (iretcode.eq.11) then
4033 write (iout,'(/a/)') &
4034 '******* Allocated time exceeded in SUMSL. The program will stop.'
4039 head_pdb=titel(:24)//':'//str_nam(ipattern)
4040 if (outpdb) call pdbout(etot,head_pdb,ipdb)
4041 if (outmol2) call mol2out(etot,head_pdb)
4043 call briefout(ithread,etot)
4045 link_end=min0(link_end,nss)
4046 write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4048 call etotal(energia)
4049 ! call enerprint(energia(0))
4052 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4053 !d write (iout,'(a,f10.5)')
4054 !d & 'RMS deviation from reference structure:',dsqrt(rms)
4056 ener(i,ithread)=energia(i)
4058 ener(n_ene_comp+1,ithread)=energia(0)
4059 ener(n_ene_comp+3,ithread)=rms
4061 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4062 ener(n_ene_comp+2,ithread)=rms
4063 ener(n_ene_comp+4,ithread)=frac
4064 ener(n_ene_comp+5,ithread)=frac_nn
4066 call write_stat_thread(ithread,ipattern,ist)
4067 ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4068 ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4069 ! & (ener(k,ithread),k=12,14)
4071 if (me.eq.king) then
4073 call pattern_receive
4074 call receive_MCM_info
4075 if (nacc_tot.ge.nthread) then
4077 'Sufficient number of conformations collected nacc_tot=',&
4078 nacc_tot,'. Stopping other processors and terminating.'
4080 'Sufficient number of conformations collected nacc_tot=',&
4081 nacc_tot,'. Stopping other processors and terminating.'
4082 call recv_stop_sig(Kwita)
4083 if (Kwita.eq.0) call send_stop_sig(-1)
4088 call send_MCM_info(2)
4091 if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4092 write (iout,'(/2a)') &
4093 '********** There would be not enough time for another thread. ',&
4094 'The program will stop.'
4096 '********** There would be not enough time for another thread. ',&
4097 'The program will stop.'
4098 write (iout,'(a,1pe14.4/)') &
4099 'Elapsed time for last threading step: ',time_for_thread
4102 call recv_stop_sig(Kwita)
4103 call send_stop_sig(-2)
4108 write (iout,'(a,1pe14.4)') &
4109 'Elapsed time for this threading step: ',time_for_thread
4112 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4113 if (Kwita.lt.0) then
4114 write (iout,*) 'Stop signal received. Terminating.'
4115 write (*,*) 'Stop signal received. Terminating.'
4117 write (*,*) 'nthread=',nthread,' ithread=',ithread
4123 call send_stop_sig(-1)
4127 ! Any messages left for me?
4128 call pattern_receive
4129 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4131 call write_thread_summary
4133 if (king.eq.king) then
4135 do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4138 call recv_stop_sig(Kwita)
4139 call receive_MCM_info
4142 call receive_thread_results(iproc)
4144 call write_thread_summary
4146 call send_thread_results
4150 end subroutine thread_seq
4151 !-----------------------------------------------------------------------------
4154 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4156 use random, only:iran_num
4157 ! implicit real*8 (a-h,o-z)
4158 ! include 'DIMENSIONS'
4159 ! include 'COMMON.CHAIN'
4160 ! include 'COMMON.DBASE'
4161 ! include 'COMMON.INTERACT'
4162 ! include 'COMMON.VAR'
4163 ! include 'COMMON.THREAD'
4164 ! include 'COMMON.FFIELD'
4165 ! include 'COMMON.SBRIDGE'
4166 ! include 'COMMON.HEADER'
4167 ! include 'COMMON.GEO'
4168 ! include 'COMMON.IOUNITS'
4169 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
4170 !el integer :: icall
4171 !el common /srutu/ icall
4172 real(kind=8) :: energia(0:n_ene)
4173 logical :: glycine,fail
4174 integer :: i,maxsample,link_end0,ind_sc,isample
4175 real(kind=8) :: alph0,omeg0,e1,e0
4179 link_end=min0(link_end,nss)
4181 if (itype(i).ne.10) then
4182 !d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
4183 call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
4187 call etotal(energia)
4189 do isample=1,maxsample
4190 ! Choose a non-glycine side chain.
4193 ind_sc=iran_num(nnt,nct)
4194 glycine=(itype(ind_sc).eq.10)
4198 call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),&
4201 call etotal(energia)
4202 !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4203 !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4204 !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4215 end subroutine sc_conf
4216 !-----------------------------------------------------------------------------
4218 !-----------------------------------------------------------------------------
4219 logical function check_var(var,info)
4223 ! implicit real*8 (a-h,o-z)
4224 ! include 'DIMENSIONS'
4225 ! include 'COMMON.VAR'
4226 ! include 'COMMON.IOUNITS'
4227 ! include 'COMMON.GEO'
4228 ! include 'COMMON.SETUP'
4229 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4230 integer,dimension(3) :: info
4234 do i=nphi+ntheta+1,nphi+ntheta+nside
4235 ! Check the side chain "valence" angles alpha
4236 if (var(i).lt.1.0d-7) then
4237 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4238 write (iout,*) 'Processor',me,'received bad variables!!!!'
4239 write (iout,*) 'Variables'
4240 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4241 write (iout,*) 'Continuing calculations at this point',&
4242 ' could destroy the results obtained so far... ABORTING!!!!!!'
4243 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4244 'valence angle alpha',i-nphi-ntheta,var(i),&
4245 'n it',info(1),info(2),'mv ',info(3)
4246 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4247 write (*,*) 'Processor',me,'received bad variables!!!!'
4248 write (*,*) 'Variables'
4249 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4250 write (*,*) 'Continuing calculations at this point',&
4251 ' could destroy the results obtained so far... ABORTING!!!!!!'
4252 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4253 'valence angle alpha',i-nphi-ntheta,var(i),&
4254 'n it',info(1),info(2),'mv ',info(3)
4259 ! Check the backbone "valence" angles theta
4260 do i=nphi+1,nphi+ntheta
4261 if (var(i).lt.1.0d-7) then
4262 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4263 write (iout,*) 'Processor',me,'received bad variables!!!!'
4264 write (iout,*) 'Variables'
4265 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4266 write (iout,*) 'Continuing calculations at this point',&
4267 ' could destroy the results obtained so far... ABORTING!!!!!!'
4268 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4269 'valence angle theta',i-nphi,var(i),&
4270 'n it',info(1),info(2),'mv ',info(3)
4271 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4272 write (*,*) 'Processor',me,'received bad variables!!!!'
4273 write (*,*) 'Variables'
4274 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4275 write (*,*) 'Continuing calculations at this point',&
4276 ' could destroy the results obtained so far... ABORTING!!!!!!'
4277 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4278 'valence angle theta',i-nphi,var(i),&
4279 'n it',info(1),info(2),'mv ',info(3)
4285 end function check_var
4286 !-----------------------------------------------------------------------------
4288 !-----------------------------------------------------------------------------
4289 subroutine distfit(debug,maxit)
4291 use geometry_data, only: phi
4294 ! implicit real*8 (a-h,o-z)
4295 ! include 'DIMENSIONS'
4296 ! include 'COMMON.CHAIN'
4297 ! include 'COMMON.VAR'
4298 ! include 'COMMON.IOUNITS'
4299 ! include 'COMMON.DISTFIT'
4300 integer :: i,maxit,MAXMAR,IT,IMAR
4301 real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres)
4302 logical :: debug,sing
4303 real(kind=8) :: TOL,RL,F0,AIN,F1
4305 !input------------------------------------
4307 ! NY=((NRES-4)*(NRES-5))/2
4308 !input------------------------------------
4314 CALL TRANSFER(NRES,phi,phiold)
4318 !d WRITE (IOUT,*) 'DISTFIT: F0=',F0
4334 CALL TRANSFER(NX,XX,X)
4335 CALL BANACH(NX,NRES,H,X,sing)
4340 IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4342 WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4343 WRITE (IOUT,*) 'IT=',it,'F=',F0
4348 phi(I)=phiold(I)+mask(i)*X(I-3)
4353 !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4355 CALL TRANSFER(NRES,phi,phiold)
4358 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4360 WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4361 WRITE (IOUT,*) 'IT=',it,'F=',F1
4367 WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4368 WRITE (IOUT,*) 'IT=',it,'F=',F0
4369 CALL TRANSFER(NRES,phiold,phi)
4372 !d write (iout,*) "it",it," imar",imar," f0",f0
4374 WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4376 end subroutine distfit
4377 !-----------------------------------------------------------------------------
4378 real(kind=8) function RDIF()
4381 use geometry, only: dist
4382 ! implicit real*8 (a-h,o-z)
4383 ! include 'DIMENSIONS'
4384 ! include 'COMMON.CHAIN'
4385 ! include 'COMMON.DISTFIT'
4387 real(kind=8) :: suma,DIJ
4396 if (w(ind).ne.0.0) then
4398 suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4400 ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4408 !-----------------------------------------------------------------------------
4413 use geometry, only:dist
4414 ! implicit real*8 (a-h,o-z)
4415 ! include 'DIMENSIONS'
4416 ! include 'COMMON.CHAIN'
4417 ! include 'COMMON.DISTFIT'
4418 ! include 'COMMON.GEO'
4419 integer :: i,j,k,l,I1,I2,IND
4420 real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4433 R13(K)=C(K,J)-C(K,I1)
4437 R24(L)=C(L,K)-C(L,I2)
4439 IND=((J-1)*(2*NRES-J-6))/2+K-3
4440 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4441 PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4442 PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4443 DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4448 end subroutine RDERIV
4449 !-----------------------------------------------------------------------------
4453 ! implicit real*8 (a-h,o-z)
4454 ! include 'DIMENSIONS'
4455 ! include 'COMMON.CHAIN'
4456 ! include 'COMMON.DISTFIT'
4458 real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4466 XI=XI+BKIWK*(D0(K)-DDD(K))
4474 HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4481 end subroutine HEVAL
4482 !-----------------------------------------------------------------------------
4483 subroutine VEC(I,J,U)
4485 use geometry_data, only: C
4486 ! Find the unit vector from atom (I) to atom (J). Store in U.
4488 ! implicit real*8 (a-h,o-z)
4489 ! include 'DIMENSIONS'
4490 ! include 'COMMON.CHAIN'
4492 real(kind=8),DIMENSION(3) :: U
4493 real(kind=8) :: ANORM,UK
4507 !-----------------------------------------------------------------------------
4508 subroutine TRANSFER(N,X1,X2)
4510 ! implicit real*8 (a-h,o-z)
4511 ! include 'DIMENSIONS'
4513 real(kind=8),DIMENSION(N) :: X1,X2
4517 end subroutine TRANSFER
4518 !-----------------------------------------------------------------------------
4519 !-----------------------------------------------------------------------------
4520 subroutine alloc_compare_arrays
4522 maxres22=nres*(nres+1)/2
4524 ! common /struct/ in io_common: read_threadbase
4525 ! allocate(cart_base !(3,maxres_base,maxseq)
4526 ! allocate(nres_base !(3,maxseq)
4527 ! allocate(str_nam !(maxseq)
4529 ! COMMON /c_frag/ in io_conf: readpdb
4530 if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4531 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4533 allocate(w(maxres22),d0(maxres22)) !(maxres22)
4535 !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4536 allocate(DDD(maxres22)) !(maxres22)
4537 allocate(H(nres,nres)) !(MAXRES,MAXRES)
4538 allocate(XX(nres)) !(MAXRES)
4540 allocate(mask(nres)) !(maxres)
4543 allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4545 allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4548 end subroutine alloc_compare_arrays
4549 !-----------------------------------------------------------------------------
4551 !-----------------------------------------------------------------------------