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,100*nres) :: icont!(2,100*nres) !(2,maxcont) (maxcont=12*maxres)
43 real(kind=8) :: co,rcomp
44 integer :: kkk,i,j,i1,i2,it1,it2,iti,itj,inum,jnum
48 do i=nnt_molec(1)+kkk,nct_molec(1)
49 iti=iabs(itype(i,molnum(i)))
50 if (molnum(i).lt.3) then
57 itj=iabs(itype(j,molnum(i)))
58 if (molnum(j).lt.3) then
64 ! rcomp=sigmaii(iti,itj)+1.0D0
65 rcomp=facont*sigmaii(iti,itj)
67 ! rcomp=sigma(iti,itj)+1.0D0
68 rcomp=facont*sigma(iti,itj)
71 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
72 if (dist(inum,jnum).lt.rcomp) then
74 if (ncont.gt.nres*100) ncont=nres*100
81 write (iout,'(a)') 'Contact map:'
87 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
88 i,restyp(it1,1),i1,restyp(it2,1),i2
93 co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
95 co = co / (nres*ncont)
97 end subroutine contact
98 !-----------------------------------------------------------------------------
99 real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref)
101 ! implicit real*8 (a-h,o-z)
102 ! include 'DIMENSIONS'
103 ! include 'COMMON.IOUNITS'
104 integer :: ncont,ncont_ref
105 integer,dimension(2,100*nres) :: icont,icont_ref !(2,100*nres) (2,maxcont) (maxcont=12*maxres)
107 integer :: i,j,nmatch
109 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
110 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
111 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
112 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
113 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
116 if (icont(1,i).eq.icont_ref(1,j) .and. &
117 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
120 ! print *,' nmatch=',nmatch
121 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
122 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
124 end function contact_fract
125 !-----------------------------------------------------------------------------
126 real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
128 ! implicit real*8 (a-h,o-z)
129 ! include 'DIMENSIONS'
130 ! include 'COMMON.IOUNITS'
131 integer :: ncont,ncont_ref
132 integer,dimension(2,100*nres) :: icont,icont_ref !(2,100*nres) (2,maxcont) (maxcont=12*maxres)
134 integer :: i,j,nmatch
136 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
137 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
138 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
139 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
140 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
143 if (icont(1,i).eq.icont_ref(1,j) .and. &
144 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
147 ! print *,' nmatch=',nmatch
148 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
149 contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
151 end function contact_fract_nn
152 !-----------------------------------------------------------------------------
153 subroutine hairpin(lprint,nharp,iharp)
155 use geometry, only:dist
156 ! implicit real*8 (a-h,o-z)
157 ! include 'DIMENSIONS'
158 ! include 'COMMON.IOUNITS'
159 ! include 'COMMON.CHAIN'
160 ! include 'COMMON.INTERACT'
161 ! include 'COMMON.FFIELD'
162 ! include 'COMMON.NAMES'
164 integer,dimension(2,100*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
166 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
167 logical :: lprint,not_done
168 real(kind=8) :: rcomp=6.0d0
170 integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1
171 ! allocate(icont(2,100*nres))
175 ! print *,'nnt=',nnt,' nct=',nct
176 do i=nnt_molec(1),nct_molec(1)-3
178 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
182 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
184 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
192 write (iout,'(a)') 'PP contact map:'
198 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
199 i,restyp(it1,1),i1,restyp(it2,1),i2
207 if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
208 ! write (iout,*) "found turn at ",i1,j1
216 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
220 ! write (iout,*) i1,j1,not_done
230 ! write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
235 ! write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
238 write (iout,*) "Hairpins:"
245 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
246 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
248 ! write (iout,'(a,i3,$)') restyp(itype(k,1)),k
253 end subroutine hairpin
254 !-----------------------------------------------------------------------------
256 !-----------------------------------------------------------------------------
257 subroutine elecont(lprint,ncont,icont)
259 ! implicit real*8 (a-h,o-z)
260 ! include 'DIMENSIONS'
261 ! include 'COMMON.IOUNITS'
262 ! include 'COMMON.CHAIN'
263 ! include 'COMMON.INTERACT'
264 ! include 'COMMON.LOCAL'
265 ! include 'COMMON.FFIELD'
266 ! include 'COMMON.NAMES'
268 real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_
269 real(kind=8) :: ael6_i,ael3_i
270 real(kind=8),dimension(2,2) :: app_,bpp_,rpp_
272 integer,dimension(2,100*nres) :: icont !(2,100*nres)(2,maxcont) (maxcont=12*maxres)
273 real(kind=8),dimension(100*nres) :: econt !(maxcont)
275 integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2
276 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw
277 real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb
278 real(kind=8) :: xmedi,ymedi,zmedi
279 real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij
280 real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,&
281 evdwij,el1,el2,eesij,ene
283 ! Load the constants of peptide bond - peptide bond interactions.
284 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
285 ! proline) - determined by averaging ECEPP energy.
289 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
290 data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
291 data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
292 data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
294 !el allocate(econt(100*nres)) !(maxcont)
297 elecutoff_14 = -0.5d0
298 if (lprint) write (iout,'(a)') &
299 "Constants of electrostatic interaction energy expression."
303 app_(i,j)=epp(i,j)*rri*rri
304 bpp_(i,j)=-2.0*epp(i,j)*rri
305 ael6_(i,j)=elpp_6(i,j)*4.2**6
306 ael3_(i,j)=elpp_3(i,j)*4.2**3
308 write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),&
315 ! print *, "nntt,nct",nnt,nct-2
316 do 1 i=nnt_molec(1),nct_molec(1)-2
317 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
328 if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
331 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
332 if (iteli.eq.2 .and. itelj.eq.2) goto 4
333 aaa=app_(iteli,itelj)
334 bbb=bpp_(iteli,itelj)
335 ael6_i=ael6_(iteli,itelj)
336 ael3_i=ael3_(iteli,itelj)
340 xj=c(1,j)+0.5*dxj-xmedi
341 yj=c(2,j)+0.5*dyj-ymedi
342 zj=c(3,j)+0.5*dzj-zmedi
343 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
348 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
349 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
350 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
351 fac=cosa-3.0*cosb*cosg
357 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
360 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
361 j.eq.i+2 .and. eesij.le.elecutoff_14) then
372 write (iout,*) 'Total average electrostatic energy: ',ees
373 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
375 write (iout,*) 'Electrostatic contacts before pruning: '
381 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
382 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
385 ! For given residues keep only the contacts with the greatest energy.
387 do while (i.lt.ncont)
393 do while (j.lt.ncont)
395 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
396 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
397 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
398 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
399 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
400 if (ic1.eq.icont(1,j)) then
402 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
403 .and. iabs(icont(1,k)-ic1).le.2 .and. &
404 econt(k).lt.econt(j) ) goto 21
406 else if (ic2.eq.icont(2,j) ) then
408 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
409 .and. iabs(icont(2,k)-ic2).le.2 .and. &
410 econt(k).lt.econt(j) ) goto 21
415 icont(1,k-1)=icont(1,k)
416 icont(2,k-1)=icont(2,k)
421 ! write (iout,*) "ncont",ncont
423 ! write (iout,*) icont(1,k),icont(2,k)
426 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
428 if (ic1.eq.icont(1,j)) then
430 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
431 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
432 econt(k).lt.econt(i) ) goto 21
434 else if (ic2.eq.icont(2,j) ) then
436 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
437 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
438 econt(k).lt.econt(i) ) goto 21
443 icont(1,k-1)=icont(1,k)
444 icont(2,k-1)=icont(2,k)
448 ! write (iout,*) "ncont",ncont
450 ! write (iout,*) icont(1,k),icont(2,k)
461 write (iout,*) 'Electrostatic contacts after pruning: '
467 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
468 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
472 end subroutine elecont
473 !-----------------------------------------------------------------------------
474 subroutine secondary2(lprint)
476 ! implicit real*8 (a-h,o-z)
477 ! include 'DIMENSIONS'
478 ! include 'COMMON.CHAIN'
479 ! include 'COMMON.IOUNITS'
480 ! include 'COMMON.DISTFIT'
481 ! include 'COMMON.VAR'
482 ! include 'COMMON.GEO'
483 ! include 'COMMON.CONTROL'
484 integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
486 integer,dimension(2,100*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
487 integer,dimension(nres,0:4) :: isec !(maxres,4)
488 integer,dimension(nres) :: nsec !(maxres)
489 logical :: lprint,not_done !,freeres
490 real(kind=8) :: p1,p2
493 !el allocate(icont(2,100*nres),isec(nres,4),nsec(nres))
495 if(.not.dccart) call chainbuild_cart
496 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
497 !d call write_pdb(99,'sec structure',0d0)
507 call elecont(lprint,ncont,icont)
508 ! print *,"after elecont"
509 if (nres_molec(1).eq.0) return
511 ! finding parallel beta
512 !d write (iout,*) '------- looking for parallel beta -----------'
518 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
521 !d write (iout,*) i1,j1
527 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
528 freeres(i1,j1,nsec,isec)) goto 5
532 !d write (iout,*) i1,j1,not_done
536 if (i1-ii1.gt.1) then
540 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
544 bfrag(1,nbfrag)=ii1+1
546 bfrag(3,nbfrag)=jj1+1
547 bfrag(4,nbfrag)=min0(j1+1,nres)
551 isec(ij,nsec(ij))=nbeta
555 isec(ij,nsec(ij))=nbeta
561 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
562 "DefPropRes 'strand",nstrand,&
563 "' 'num = ",ii1-1,"..",i1-1,"'"
565 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
566 "DefPropRes 'strand",nstrand,&
567 "' 'num = ",ii1-1,"..",i1-1,"'"
571 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
572 "DefPropRes 'strand",nstrand,&
573 "' 'num = ",jj1-1,"..",j1-1,"'"
575 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
576 "DefPropRes 'strand",nstrand,&
577 "' 'num = ",jj1-1,"..",j1-1,"'"
579 write(12,'(a8,4i4)') &
580 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
586 ! finding alpha or 310 helix
593 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
596 if (j1.eq.i1+3 .and. &
597 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
598 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
599 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
600 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
603 if (nsec(ii1).eq.0) then
612 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
618 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
623 if (j1-ii1.gt.5) then
635 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
636 if (nhelix.le.9) then
637 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
638 "DefPropRes 'helix",nhelix,&
639 "' 'num = ",ii1-1,"..",j1-2,"'"
641 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
642 "DefPropRes 'helix",nhelix,&
643 "' 'num = ",ii1-1,"..",j1-2,"'"
649 if (nhelix.gt.0.and.lprint) then
650 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
652 if (nhelix.le.9) then
653 write(12,'(a8,i1,$)') " | helix",i
655 write(12,'(a8,i2,$)') " | helix",i
662 ! finding antiparallel beta
663 !d write (iout,*) '--------- looking for antiparallel beta ---------'
668 if (freeres(i1,j1,nsec,isec)) then
671 !d write (iout,*) i1,j1
678 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
679 freeres(i1,j1,nsec,isec)) goto 6
683 !d write (iout,*) i1,j1,not_done
687 if (i1-ii1.gt.1) then
691 bfrag(2,nbfrag)=min0(i1+1,nres)
692 bfrag(3,nbfrag)=min0(jj1+1,nres)
699 if (nsec(ij).le.2) then
700 isec(ij,nsec(ij))=nbeta
706 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
707 isec(ij,nsec(ij))=nbeta
713 write (iout,'(a,i3,4i4)')'antiparallel beta',&
714 nbeta,ii1-1,i1,jj1,j1-1
716 if (nstrand.le.9) then
717 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
718 "DefPropRes 'strand",nstrand,&
719 "' 'num = ",ii1-2,"..",i1-1,"'"
721 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
722 "DefPropRes 'strand",nstrand,&
723 "' 'num = ",ii1-2,"..",i1-1,"'"
726 if (nstrand.le.9) then
727 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
728 "DefPropRes 'strand",nstrand,&
729 "' 'num = ",j1-2,"..",jj1-1,"'"
731 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
732 "DefPropRes 'strand",nstrand,&
733 "' 'num = ",j1-2,"..",jj1-1,"'"
735 write(12,'(a8,4i4)') &
736 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
742 if (nstrand.gt.0.and.lprint) then
743 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
746 write(12,'(a9,i1,$)') " | strand",i
748 write(12,'(a9,i2,$)') " | strand",i
757 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
758 write(12,'(a20)') "XMacStand ribbon.mac"
760 if (nres_molec(1).eq.0) return
761 write(iout,*) 'UNRES seq:'
763 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
767 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
772 end subroutine secondary2
774 !-----------------------------------------------------------------------------
775 logical function freeres(i,j,nsec,isec)
777 ! implicit real*8 (a-h,o-z)
778 ! include 'DIMENSIONS'
779 integer,dimension(nres,4) :: isec !(maxres,4)
780 integer,dimension(nres) :: nsec !(maxres)
787 if (nsec(i).lt.0.or.nsec(j).lt.0) return
789 if (nsec(i).gt.1.or.nsec(j).gt.1) return
792 if (isec(i,k).eq.isec(j,l)) return
798 !-----------------------------------------------------------------------------
800 !-----------------------------------------------------------------------------
801 logical function seq_comp(itypea,itypeb,length)
804 integer :: length,itypea(length),itypeb(length)
807 if (itypea(i).ne.itypeb(i)) then
814 end function seq_comp
816 !-----------------------------------------------------------------------------
818 !-----------------------------------------------------------------------------
819 subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
821 ! implicit real*8 (a-h,o-z)
822 ! include 'DIMENSIONS'
823 ! include 'COMMON.CHAIN'
824 ! include 'COMMON.CONTACTS'
825 ! include 'COMMON.IOUNITS'
826 real(kind=8) :: przes(3),obr(3,3)
827 logical :: non_conv,lprn
828 real(kind=8) :: rms,frac,frac_nn,co
829 ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
833 ! print *,"before contact"
834 !elte(iout,*) "rms_nacc before contact"
835 call contact(.false.,ncont,icont,co)
836 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
837 frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
838 if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
839 'RMS deviation from the reference structure:',rms,&
840 ' % of native contacts:',frac*100,&
841 ' % of nonnative contacts:',frac_nn*100,&
845 end subroutine rms_nac_nnc
846 !-----------------------------------------------------------------------------
847 subroutine rmsd(drms)
849 use regularize_, only:fitsq
850 ! implicit real*8 (a-h,o-z)
851 ! include 'DIMENSIONS'
855 ! include 'COMMON.CHAIN'
856 ! include 'COMMON.IOUNITS'
857 ! include 'COMMON.INTERACT'
858 ! include 'COMMON.CONTROL'
860 real(kind=8) :: przes(3),obrot(3,3)
861 real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
864 real(kind=8) :: drms,rminroz,roznica
865 integer :: i,j,iatom,kkk,iti,k,mnum
867 !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
875 ! print *,"nz_start",nz_start," nz_end",nz_end
876 ! if (symetr.le.1) then
878 ! do i=nz_start,nz_end
882 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
883 ! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
885 ! if (iz_sc.eq.1.and.iti.ne.10) then
888 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
889 ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
896 ! print *,nz_start,nz_end,nstart_seq-nstart_sup
899 iti=itype(i+nstart_seq-nstart_sup,molnum(i))
900 if (iti.eq.ntyp1_molec(molnum(i))) cycle
902 mnum=molnum(i+nstart_seq-nstart_sup)
904 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
905 crefcopy(k,iatom)=cref(k,i,kkk)
907 if (iz_sc.eq.1.and.iti.ne.10.and.mnum.lt.4) then
910 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
911 crefcopy(k,iatom)=cref(k,nres+i,kkk)
920 ! write (iout,*) 'Ccopy and CREFcopy'
921 ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
922 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
923 ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
924 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
926 ! ----- end diagnostics
928 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
929 przes,obrot,non_conv)
931 print *,'Problems in FITSQ!!! rmsd'
932 write (iout,*) 'Problems in FITSQ!!! rmsd'
933 print *,'Ccopy and CREFcopy'
934 write (iout,*) 'Ccopy and CREFcopy'
935 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
936 (crefcopy(j,k),j=1,3),k=1,iatom)
937 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
938 (crefcopy(j,k),j=1,3),k=1,iatom)
940 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
946 ! write (iout,*) "roznica", roznica,kkk
947 if (roznica.le.rminroz) rminroz=roznica
949 drms=dsqrt(dabs(rminroz))
951 ! write (iout,*) "nperm,symetr", nperm,symetr
952 ! ---- end diagnostics
955 !-----------------------------------------------------------------------------
956 subroutine rmsd_csa(drms)
958 use regularize_, only:fitsq
959 ! implicit real*8 (a-h,o-z)
960 ! include 'DIMENSIONS'
964 ! include 'COMMON.CHAIN'
965 ! include 'COMMON.IOUNITS'
966 ! include 'COMMON.INTERACT'
968 real(kind=8) :: przes(3),obrot(3,3)
969 real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
970 integer :: kkk,iatom,ierror,ierrcode
974 real(kind=8) :: drms,roznica
976 allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
984 ccopy(k,iatom)=c(k,i)
985 crefcopy(k,iatom)=crefjlee(k,i)
987 if (iz_sc.eq.1.and.iti.ne.10) then
990 ccopy(k,iatom)=c(k,nres+i)
991 crefcopy(k,iatom)=crefjlee(k,nres+i)
996 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
997 przes,obrot,non_conv)
999 print *,'Problems in FITSQ!!! rmsd_csa'
1000 write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
1001 print *,'Ccopy and CREFcopy'
1002 write (iout,*) 'Ccopy and CREFcopy'
1003 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
1004 (crefcopy(j,k),j=1,3),k=1,iatom)
1005 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
1006 (crefcopy(j,k),j=1,3),k=1,iatom)
1008 call mpi_abort(mpi_comm_world,ierror,ierrcode)
1013 drms=dsqrt(dabs(roznica))
1015 end subroutine rmsd_csa
1016 !-----------------------------------------------------------------------------
1018 !-----------------------------------------------------------------------------
1022 use geometry, only:pinorm
1023 use random, only:ran_number,iran_num
1024 ! implicit real*8 (a-h,o-z)
1025 ! include 'DIMENSIONS'
1027 ! include 'COMMON.GEO'
1028 ! include 'COMMON.VAR'
1029 ! include 'COMMON.INTERACT'
1030 ! include 'COMMON.IOUNITS'
1031 ! include 'COMMON.DISTFIT'
1032 ! include 'COMMON.SBRIDGE'
1033 ! include 'COMMON.CONTROL'
1034 ! include 'COMMON.FFIELD'
1035 ! include 'COMMON.MINIM'
1036 ! include 'COMMON.CHAIN'
1037 real(kind=8) :: time0,time1
1038 real(kind=8) :: energy(0:n_ene),ee
1039 real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres)
1040 integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1041 logical :: debug,accepted
1042 real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1045 !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1047 call geom_to_var(nvar,var1)
1052 write(iout,*) 'etot=',0,etot,rms
1053 call secondary2(.false.)
1055 call write_pdb(0,'first structure',etot)
1064 betbol=1.0D0/(1.9858D-3*temp)
1066 d=ran_number(-pi,pi)
1067 ! phi(jr)=pinorm(phi(jr)+d)
1072 write(iout,*) 'etot=',1,etot0,rms
1073 call write_pdb(1,'perturb structure',etot0)
1077 d=ran_number(-da,da)
1079 phi(jr)=pinorm(phi(jr)+d)
1084 if (etot.lt.etot0) then
1088 xxr=ran_number(0.0D0,1.0D0)
1089 xxh=betbol*(etot-etot0)
1090 if (xxh.lt.50.0D0) then
1092 if (xxh.gt.xxr) accepted=.true.
1096 ! print *,etot0,etot,accepted
1100 write(iout,*) 'etot=',i,etot,rms
1101 call write_pdb(i,'MC structure',etot)
1103 ! call geom_to_var(nvar,var1)
1104 call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1105 call geom_to_var(nvar,var)
1106 call minimize(etot,var,iretcode,nfun)
1107 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1108 call var_to_geom(nvar,var)
1111 write(iout,*) 'etot mcm=',i,etot,rms
1112 call write_pdb(i+1,'MCM structure',etot)
1113 call var_to_geom(nvar,var1)
1121 ! call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1122 ! call geom_to_var(nvar,var)
1125 ! call write_pdb(998 ,'sc min',etot)
1127 ! call minimize(etot,var,iretcode,nfun)
1128 ! write(iout,*)'------------------------------------------------'
1129 ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1131 ! call var_to_geom(nvar,var)
1133 ! call write_pdb(999,'full min',etot)
1137 !-----------------------------------------------------------------------------
1142 ! implicit real*8 (a-h,o-z)
1143 ! include 'DIMENSIONS'
1145 ! include 'COMMON.GEO'
1146 ! include 'COMMON.VAR'
1147 ! include 'COMMON.INTERACT'
1148 ! include 'COMMON.IOUNITS'
1149 ! include 'COMMON.DISTFIT'
1150 ! include 'COMMON.SBRIDGE'
1151 ! include 'COMMON.CONTROL'
1152 ! include 'COMMON.FFIELD'
1153 ! include 'COMMON.MINIM'
1154 ! include 'COMMON.CHAIN'
1155 real(kind=8) :: time0,time1
1156 real(kind=8) :: energy(0:n_ene),ee
1157 real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1161 integer :: i,ij,ieval,iretcode,nfun
1162 real(kind=8) :: etot
1164 allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1166 call geom_to_var(nvar,var1)
1170 write(iout,*) nnt,nct,etot
1171 call write_pdb(1,'first structure',etot)
1172 call secondary2(.true.)
1181 call var_to_geom(nvar,var1)
1182 write(iout,*) 'N16 test',(jdata(i),i=1,5)
1183 call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1185 call geom_to_var(nvar,var)
1191 call minimize(etot,var,iretcode,nfun)
1192 write(iout,*)'------------------------------------------------'
1193 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1199 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1200 nfun/(time1-time0),' eval/s'
1202 call var_to_geom(nvar,var)
1204 call write_pdb(ij*100+99,'full min',etot)
1211 end subroutine test_n16
1213 !-----------------------------------------------------------------------------
1214 subroutine test_local
1216 ! implicit real*8 (a-h,o-z)
1217 ! include 'DIMENSIONS'
1218 ! include 'COMMON.GEO'
1219 ! include 'COMMON.VAR'
1220 ! include 'COMMON.INTERACT'
1221 ! include 'COMMON.IOUNITS'
1222 real(kind=8) :: time0,time1
1223 real(kind=8) :: energy(0:n_ene),ee
1224 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1226 real(kind=8) :: etot
1228 ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres)
1230 ! call geom_to_var(nvar,varia)
1231 call write_pdb(1,'first structure',0d0)
1235 write(iout,*) nnt,nct,etot
1237 write(iout,*) 'calling sc_move'
1238 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1239 write(iout,*) nft_sc,etot
1240 call write_pdb(2,'second structure',etot)
1242 write(iout,*) 'calling local_move'
1243 call local_move_init(.false.)
1244 call local_move(24,29,20d0,50d0)
1246 call write_pdb(3,'third structure',etot)
1248 write(iout,*) 'calling sc_move'
1249 call sc_move(24,29,5,10d0,nft_sc,etot)
1250 write(iout,*) nft_sc,etot
1251 call write_pdb(2,'last structure',etot)
1254 end subroutine test_local
1255 !-----------------------------------------------------------------------------
1258 ! implicit real*8 (a-h,o-z)
1259 ! include 'DIMENSIONS'
1260 ! include 'COMMON.GEO'
1261 ! include 'COMMON.VAR'
1262 ! include 'COMMON.INTERACT'
1263 ! include 'COMMON.IOUNITS'
1264 real(kind=8) :: time0,time1,etot
1265 real(kind=8) :: energy(0:n_ene),ee
1266 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1270 ! call geom_to_var(nvar,varia)
1271 call write_pdb(1,'first structure',0d0)
1275 write(iout,*) nnt,nct,etot
1277 write(iout,*) 'calling sc_move'
1279 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1280 write(iout,*) nft_sc,etot
1281 call write_pdb(2,'second structure',etot)
1283 write(iout,*) 'calling sc_move 2nd time'
1285 call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1286 write(iout,*) nft_sc,etot
1287 call write_pdb(3,'last structure',etot)
1289 end subroutine test_sc
1290 !-----------------------------------------------------------------------------
1291 subroutine bgrow(bstrand,nbstrand,in,ind,new)
1293 ! implicit real*8 (a-h,o-z)
1294 ! include 'DIMENSIONS'
1295 ! include 'COMMON.CHAIN'
1296 integer,dimension(nres/3,6) :: bstrand !(maxres/3,6)
1299 integer :: nbstrand,in,ind,new,ishift,i
1301 ishift=iabs(bstrand(in,ind+4)-new)
1303 print *,'bgrow',bstrand(in,ind+4),new,ishift
1308 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1310 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1311 if (bstrand(i,5).lt.bstrand(i,6)) then
1312 bstrand(i,5)=bstrand(i,5)-ishift
1314 bstrand(i,5)=bstrand(i,5)+ishift
1319 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1321 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1322 if (bstrand(i,6).lt.bstrand(i,5)) then
1323 bstrand(i,6)=bstrand(i,6)-ishift
1325 bstrand(i,6)=bstrand(i,6)+ishift
1332 end subroutine bgrow
1333 !-----------------------------------------------------------------------------
1336 use geometry, only:dist
1338 ! implicit real*8 (a-h,o-z)
1339 ! include 'DIMENSIONS'
1341 ! include 'COMMON.GEO'
1342 ! include 'COMMON.CHAIN'
1343 ! include 'COMMON.IOUNITS'
1344 ! include 'COMMON.VAR'
1345 ! include 'COMMON.CONTROL'
1346 ! include 'COMMON.SBRIDGE'
1347 ! include 'COMMON.FFIELD'
1348 ! include 'COMMON.MINIM'
1350 ! include 'COMMON.DISTFIT'
1351 integer :: if(20,nres),nif,ifa(20)
1352 integer :: ibc(0:nres,0:nres),istrand(20)
1353 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1354 integer :: itmp(20,nres)
1355 real(kind=8) :: time0,time1
1356 real(kind=8) :: energy(0:n_ene),ee
1357 real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres)
1359 logical :: debug,ltest,usedbfrag(nres/3)
1360 character(len=50) :: linia
1362 integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1363 integer :: bstrand(nres/3,6),nbstrand
1364 real(kind=8) :: etot
1365 integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1366 in,ind,ifun,nfun,iretcode
1367 !------------------------
1370 !------------------------
1377 call geom_to_var(nvar,vorg)
1378 call secondary2(debug)
1380 if (nbfrag.le.1) return
1383 usedbfrag(i)=.false.
1387 nbetasheet=nbetasheet+1
1389 bstrand(1,1)=bfrag(1,1)
1390 bstrand(1,2)=bfrag(2,1)
1391 bstrand(1,3)=nbetasheet
1393 bstrand(1,5)=bfrag(1,1)
1394 bstrand(1,6)=bfrag(2,1)
1395 do i=bfrag(1,1),bfrag(2,1)
1396 betasheet(i)=nbetasheet
1400 bstrand(2,1)=bfrag(3,1)
1401 bstrand(2,2)=bfrag(4,1)
1402 bstrand(2,3)=nbetasheet
1403 bstrand(2,5)=bfrag(3,1)
1404 bstrand(2,6)=bfrag(4,1)
1406 if (bfrag(3,1).le.bfrag(4,1)) then
1408 do i=bfrag(3,1),bfrag(4,1)
1409 betasheet(i)=nbetasheet
1414 do i=bfrag(4,1),bfrag(3,1)
1415 betasheet(i)=nbetasheet
1422 do while (iused_nbfrag.ne.nbfrag)
1426 IF (.not.usedbfrag(j)) THEN
1428 write (*,*) j,(bfrag(i,j),i=1,4)
1430 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1432 write (*,*) '------------------'
1435 if (bfrag(3,j).le.bfrag(4,j)) then
1436 do i=bfrag(3,j),bfrag(4,j)
1437 if(betasheet(i).eq.nbetasheet) then
1439 do k=bfrag(3,j),bfrag(4,j)
1440 betasheet(k)=nbetasheet
1445 iused_nbfrag=iused_nbfrag+1
1446 do k=bfrag(1,j),bfrag(2,j)
1447 betasheet(k)=nbetasheet
1448 ibetasheet(k)=nbstrand
1450 if (bstrand(in,4).lt.0) then
1451 bstrand(nbstrand,1)=bfrag(2,j)
1452 bstrand(nbstrand,2)=bfrag(1,j)
1453 bstrand(nbstrand,3)=nbetasheet
1454 bstrand(nbstrand,4)=-nbstrand
1455 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1456 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1457 if(bstrand(in,1).lt.bfrag(4,j)) then
1458 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1460 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1461 (bstrand(in,5)-bfrag(4,j))
1463 if(bstrand(in,2).gt.bfrag(3,j)) then
1464 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1466 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1467 (-bstrand(in,6)+bfrag(3,j))
1470 bstrand(nbstrand,1)=bfrag(1,j)
1471 bstrand(nbstrand,2)=bfrag(2,j)
1472 bstrand(nbstrand,3)=nbetasheet
1473 bstrand(nbstrand,4)=nbstrand
1474 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1475 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1476 if(bstrand(in,1).gt.bfrag(3,j)) then
1477 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1479 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1480 (-bstrand(in,5)+bfrag(3,j))
1482 if(bstrand(in,2).lt.bfrag(4,j)) then
1483 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1485 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1486 (bstrand(in,6)-bfrag(4,j))
1491 if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1492 in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1493 do k=bfrag(1,j),bfrag(2,j)
1494 betasheet(k)=nbetasheet
1499 iused_nbfrag=iused_nbfrag+1
1500 do k=bfrag(3,1),bfrag(4,1)
1501 betasheet(k)=nbetasheet
1502 ibetasheet(k)=nbstrand
1504 if (bstrand(in,4).lt.0) then
1505 bstrand(nbstrand,1)=bfrag(4,j)
1506 bstrand(nbstrand,2)=bfrag(3,j)
1507 bstrand(nbstrand,3)=nbetasheet
1508 bstrand(nbstrand,4)=-nbstrand
1509 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1510 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1511 if(bstrand(in,1).lt.bfrag(2,j)) then
1512 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1514 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1515 (bstrand(in,5)-bfrag(2,j))
1517 if(bstrand(in,2).gt.bfrag(1,j)) then
1518 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1520 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1521 (-bstrand(in,6)+bfrag(1,j))
1524 bstrand(nbstrand,1)=bfrag(3,j)
1525 bstrand(nbstrand,2)=bfrag(4,j)
1526 bstrand(nbstrand,3)=nbetasheet
1527 bstrand(nbstrand,4)=nbstrand
1528 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1529 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1530 if(bstrand(in,1).gt.bfrag(1,j)) then
1531 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1533 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1534 (-bstrand(in,5)+bfrag(1,j))
1536 if(bstrand(in,2).lt.bfrag(2,j)) then
1537 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1539 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1540 (bstrand(in,6)-bfrag(2,j))
1547 do i=bfrag(4,j),bfrag(3,j)
1548 if(betasheet(i).eq.nbetasheet) then
1550 do k=bfrag(4,j),bfrag(3,j)
1551 betasheet(k)=nbetasheet
1556 iused_nbfrag=iused_nbfrag+1
1557 do k=bfrag(1,j),bfrag(2,j)
1558 betasheet(k)=nbetasheet
1559 ibetasheet(k)=nbstrand
1561 if (bstrand(in,4).lt.0) then
1562 bstrand(nbstrand,1)=bfrag(1,j)
1563 bstrand(nbstrand,2)=bfrag(2,j)
1564 bstrand(nbstrand,3)=nbetasheet
1565 bstrand(nbstrand,4)=nbstrand
1566 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1567 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1568 if(bstrand(in,1).lt.bfrag(3,j)) then
1569 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1571 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1572 (bstrand(in,5)-bfrag(3,j))
1574 if(bstrand(in,2).gt.bfrag(4,j)) then
1575 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1577 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1578 (-bstrand(in,6)+bfrag(4,j))
1581 bstrand(nbstrand,1)=bfrag(2,j)
1582 bstrand(nbstrand,2)=bfrag(1,j)
1583 bstrand(nbstrand,3)=nbetasheet
1584 bstrand(nbstrand,4)=-nbstrand
1585 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1586 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1587 if(bstrand(in,1).gt.bfrag(4,j)) then
1588 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1590 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1591 (-bstrand(in,5)+bfrag(4,j))
1593 if(bstrand(in,2).lt.bfrag(3,j)) then
1594 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1596 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1597 (bstrand(in,6)-bfrag(3,j))
1602 if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1603 in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1604 do k=bfrag(1,j),bfrag(2,j)
1605 betasheet(k)=nbetasheet
1610 iused_nbfrag=iused_nbfrag+1
1611 do k=bfrag(4,j),bfrag(3,j)
1612 betasheet(k)=nbetasheet
1613 ibetasheet(k)=nbstrand
1615 if (bstrand(in,4).lt.0) then
1616 bstrand(nbstrand,1)=bfrag(4,j)
1617 bstrand(nbstrand,2)=bfrag(3,j)
1618 bstrand(nbstrand,3)=nbetasheet
1619 bstrand(nbstrand,4)=nbstrand
1620 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1621 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1622 if(bstrand(in,1).lt.bfrag(2,j)) then
1623 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1625 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1626 (bstrand(in,5)-bfrag(2,j))
1628 if(bstrand(in,2).gt.bfrag(1,j)) then
1629 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1631 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1632 (-bstrand(in,6)+bfrag(1,j))
1635 bstrand(nbstrand,1)=bfrag(3,j)
1636 bstrand(nbstrand,2)=bfrag(4,j)
1637 bstrand(nbstrand,3)=nbetasheet
1638 bstrand(nbstrand,4)=-nbstrand
1639 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1640 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1641 if(bstrand(in,1).gt.bfrag(1,j)) then
1642 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1644 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1645 (-bstrand(in,5)+bfrag(1,j))
1647 if(bstrand(in,2).lt.bfrag(2,j)) then
1648 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1650 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1651 (bstrand(in,6)-bfrag(2,j))
1665 do while (usedbfrag(j))
1670 nbetasheet=nbetasheet+1
1671 bstrand(nbstrand,1)=bfrag(1,j)
1672 bstrand(nbstrand,2)=bfrag(2,j)
1673 bstrand(nbstrand,3)=nbetasheet
1674 bstrand(nbstrand,5)=bfrag(1,j)
1675 bstrand(nbstrand,6)=bfrag(2,j)
1677 bstrand(nbstrand,4)=nbstrand
1678 do i=bfrag(1,j),bfrag(2,j)
1679 betasheet(i)=nbetasheet
1680 ibetasheet(i)=nbstrand
1684 bstrand(nbstrand,1)=bfrag(3,j)
1685 bstrand(nbstrand,2)=bfrag(4,j)
1686 bstrand(nbstrand,3)=nbetasheet
1687 bstrand(nbstrand,5)=bfrag(3,j)
1688 bstrand(nbstrand,6)=bfrag(4,j)
1690 if (bfrag(3,j).le.bfrag(4,j)) then
1691 bstrand(nbstrand,4)=nbstrand
1692 do i=bfrag(3,j),bfrag(4,j)
1693 betasheet(i)=nbetasheet
1694 ibetasheet(i)=nbstrand
1697 bstrand(nbstrand,4)=-nbstrand
1698 do i=bfrag(4,j),bfrag(3,j)
1699 betasheet(i)=nbetasheet
1700 ibetasheet(i)=nbstrand
1704 iused_nbfrag=iused_nbfrag+1
1710 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1717 if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1721 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1724 !------------------------
1728 if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1729 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1731 ifb(nifb,1)=bstrand(i,4)
1732 ifb(nifb,2)=bstrand(j,4)
1739 write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1745 write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1747 nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1749 if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1750 nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1756 if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1757 if (if(j,i).gt.0) then
1758 if(betasheet(if(j,i)).eq.0 .or. &
1759 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1764 write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1767 ! read (inp,*) (ifa(i),i=1,4)
1769 ! read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1773 !------------------------
1778 !ccccccccccccccccccccccccccccccccc
1780 !ccccccccccccccccccccccccccccccccc
1784 istrand(is-j+1)=int(ii/is**(is-j))
1785 ii=ii-istrand(is-j+1)*is**(is-j)
1789 istrand(k)=istrand(k)+1
1790 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1794 if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1795 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1804 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1805 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1806 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1807 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1813 if (mod(isa,2).eq.0) then
1815 if (istrand(k).eq.1) ltest=.false.
1818 do k=(isa+1)/2+1,isa
1819 if (istrand(k).eq.1) ltest=.false.
1823 IF (ltest.and.lifb0.eq.1) THEN
1826 call var_to_geom(nvar,vorg)
1828 write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1829 write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1830 write (linia,'(10i3)') (istrand(k),k=1,isa)
1840 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1842 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1846 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1852 write(*,*) (itmp(j,i),j=1,4)
1856 ! ifa(1),ifa(2),ifa(3),ifa(4)
1857 ! if(1,i),if(2,i),if(3,i),if(4,i)
1862 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1863 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1864 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1865 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1873 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1875 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1879 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1881 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1884 !------------------------
1887 ! freeze sec.elements
1897 do i=bfrag(1,j),bfrag(2,j)
1902 if (bfrag(3,j).le.bfrag(4,j)) then
1903 do i=bfrag(3,j),bfrag(4,j)
1909 do i=bfrag(4,j),bfrag(3,j)
1917 do i=hfrag(1,j),hfrag(2,j)
1925 !------------------------
1926 ! generate constrains
1934 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1942 else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1950 else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1958 else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1966 else if ( ibc(i,j).gt.0 ) then
1967 d0(ind)=DIST(i,ibc(i,j))
1974 else if ( ibc(j,i).gt.0 ) then
1975 d0(ind)=DIST(ibc(j,i),j)
1989 !d--------------------------
1991 write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1992 ibc(jhpb(i),ihpb(i)),' --',&
1993 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2000 call contact_cp_min(varia,ifun,iconf,linia,debug)
2005 call minimize(etot,varia,iretcode,nfun)
2006 write(iout,*)'------------------------------------------------'
2007 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2013 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2014 nfun/(time1-time0),' eval/s'
2016 write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
2017 call var_to_geom(nvar,varia)
2019 call write_pdb(900+iconf,linia,etot)
2024 call enerprint(energy)
2026 !d call briefout(0,etot)
2027 !d call secondary2(.true.)
2031 !ccccccccccccccccccccccccccccccccccc
2034 !ccccccccccccccccccccccccccccccccccc
2037 10 write (iout,'(a)') 'Error reading test structure.'
2039 end subroutine test11
2040 !-----------------------------------------------------------------------------
2043 use geometry, only:dist
2045 ! implicit real*8 (a-h,o-z)
2046 ! include 'DIMENSIONS'
2048 ! include 'COMMON.GEO'
2049 ! include 'COMMON.CHAIN'
2050 ! include 'COMMON.IOUNITS'
2051 ! include 'COMMON.VAR'
2052 ! include 'COMMON.CONTROL'
2053 ! include 'COMMON.SBRIDGE'
2054 ! include 'COMMON.FFIELD'
2055 ! include 'COMMON.MINIM'
2057 ! include 'COMMON.DISTFIT'
2058 integer :: if(3,nres),nif
2059 integer :: ibc(nres,nres),istrand(20)
2060 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2061 real(kind=8) :: time0,time1
2062 real(kind=8) :: energy(0:n_ene),ee
2063 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2065 logical :: debug,ltest
2066 character(len=50) :: linia
2067 integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2068 real(kind=8) :: etot
2071 read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2074 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2078 !------------------------
2079 call secondary2(debug)
2080 !------------------------
2088 ! freeze sec.elements and store indexes for beta constrains
2098 do i=bfrag(1,j),bfrag(2,j)
2103 if (bfrag(3,j).le.bfrag(4,j)) then
2104 do i=bfrag(3,j),bfrag(4,j)
2108 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2111 do i=bfrag(4,j),bfrag(3,j)
2115 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2120 do i=hfrag(1,j),hfrag(2,j)
2129 ! ---------------- test --------------
2131 if (ibc(if(1,i),if(2,i)).eq.-1) then
2132 ibc(if(1,i),if(2,i))=if(3,i)
2133 ibc(if(1,i),if(3,i))=if(2,i)
2134 else if (ibc(if(2,i),if(1,i)).eq.-1) then
2135 ibc(if(2,i),if(1,i))=0
2136 ibc(if(1,i),if(2,i))=if(3,i)
2137 ibc(if(1,i),if(3,i))=if(2,i)
2139 ibc(if(1,i),if(2,i))=if(3,i)
2140 ibc(if(1,i),if(3,i))=if(2,i)
2146 if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j)
2149 !------------------------
2155 if ( ibc(i,j).eq.-1 ) then
2163 else if ( ibc(i,j).gt.0 ) then
2164 d0(ind)=DIST(i,ibc(i,j))
2171 else if ( ibc(j,i).gt.0 ) then
2172 d0(ind)=DIST(ibc(j,i),j)
2186 !d--------------------------
2187 write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2188 ibc(jhpb(i),ihpb(i)),' --',&
2189 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2197 call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2202 call minimize(etot,varia,iretcode,nfun)
2203 write(iout,*)'------------------------------------------------'
2204 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2210 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2211 nfun/(time1-time0),' eval/s'
2214 call var_to_geom(nvar,varia)
2216 call write_pdb(999,'full min',etot)
2221 call enerprint(energy)
2223 call briefout(0,etot)
2224 call secondary2(.true.)
2227 10 write (iout,'(a)') 'Error reading test structure.'
2229 end subroutine test3
2230 !-----------------------------------------------------------------------------
2234 ! implicit real*8 (a-h,o-z)
2235 ! include 'DIMENSIONS'
2237 ! include 'COMMON.GEO'
2238 ! include 'COMMON.CHAIN'
2239 ! include 'COMMON.IOUNITS'
2240 ! include 'COMMON.VAR'
2241 ! include 'COMMON.CONTROL'
2242 ! include 'COMMON.SBRIDGE'
2243 ! include 'COMMON.FFIELD'
2244 ! include 'COMMON.MINIM'
2246 ! include 'COMMON.DISTFIT'
2247 integer :: if(2,2),ind
2248 integer :: iff(nres)
2249 real(kind=8) :: time0,time1
2250 real(kind=8) :: energy(0:n_ene),ee
2251 real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2252 theta1,phi1,alph1,omeg1 !(maxres)
2253 real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres)
2255 integer :: i,j,nn,ifun,iretcode,nfun
2256 real(kind=8) :: etot
2259 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2260 write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2261 read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2262 read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2263 read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2264 read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2266 theta2(i)=deg2rad*theta2(i)
2267 phi2(i)=deg2rad*phi2(i)
2268 alph2(i)=deg2rad*alph2(i)
2269 omeg2(i)=deg2rad*omeg2(i)
2283 !------------------------
2288 do i=if(j,1),if(j,2)
2294 call geom_to_var(nvar,varia)
2295 call write_pdb(1,'first structure',0d0)
2297 call secondary(.true.)
2299 call secondary2(.true.)
2302 if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2303 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2304 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2307 if (bfrag(3,j).lt.bfrag(4,j)) then
2308 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2309 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2310 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2312 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2313 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2314 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2327 call geom_to_var(nvar,varia2)
2328 call write_pdb(2,'second structure',0d0)
2332 !-------------------------------------------------------
2335 call contact_cp(varia,varia2,iff,ifun,7)
2340 call minimize(etot,varia,iretcode,nfun)
2341 write(iout,*)'------------------------------------------------'
2342 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2348 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2349 nfun/(time1-time0),' eval/s'
2352 call var_to_geom(nvar,varia)
2354 call write_pdb(999,'full min',etot)
2359 call enerprint(energy)
2361 call briefout(0,etot)
2364 10 write (iout,'(a)') 'Error reading test structure.'
2366 end subroutine test__
2367 !-----------------------------------------------------------------------------
2368 subroutine secondary(lprint)
2370 ! implicit real*8 (a-h,o-z)
2371 ! include 'DIMENSIONS'
2372 ! include 'COMMON.CHAIN'
2373 ! include 'COMMON.IOUNITS'
2374 ! include 'COMMON.DISTFIT'
2376 integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2377 logical :: lprint,not_done
2378 real(kind=4) :: dcont(nres*nres/2),d
2379 real(kind=4) :: rcomp = 7.0
2380 real(kind=4) :: rbeta = 5.2
2381 real(kind=4) :: ralfa = 5.2
2382 real(kind=4) :: r310 = 6.6
2383 real(kind=8),dimension(3) :: xpi,xpj
2384 integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2386 call chainbuild_cart
2387 !d call write_pdb(99,'sec structure',0d0)
2399 xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2403 xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2405 !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2406 !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2407 !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2408 !d print *,'CA',i,j,d
2409 d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2410 (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2411 (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2412 if ( d.lt.rcomp*rcomp) then
2416 dcont(ncont)=sqrt(d)
2422 write (iout,'(a)') '#PP contact map distances:'
2424 write (iout,'(3i4,f10.5)') &
2425 i,icont(1,i),icont(2,i),dcont(i)
2429 ! finding parallel beta
2430 !d write (iout,*) '------- looking for parallel beta -----------'
2436 if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2437 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2438 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2439 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2440 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2441 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2445 !d write (iout,*) i1,j1,dcont(i)
2451 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2452 .and. dcont(j).le.rbeta .and. &
2453 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2454 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2455 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2456 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2457 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2462 !d write (iout,*) i1,j1,dcont(j),not_done
2466 if (i1-ii1.gt.1) then
2470 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2479 isec(ij,1)=isec(ij,1)+1
2480 isec(ij,1+isec(ij,1))=nbeta
2483 isec(ij,1)=isec(ij,1)+1
2484 isec(ij,1+isec(ij,1))=nbeta
2489 if (nbeta.le.9) then
2490 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2491 "DefPropRes 'strand",nstrand,&
2492 "' 'num = ",ii1-1,"..",i1-1,"'"
2494 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2495 "DefPropRes 'strand",nstrand,&
2496 "' 'num = ",ii1-1,"..",i1-1,"'"
2499 if (nbeta.le.9) then
2500 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2501 "DefPropRes 'strand",nstrand,&
2502 "' 'num = ",jj1-1,"..",j1-1,"'"
2504 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2505 "DefPropRes 'strand",nstrand,&
2506 "' 'num = ",jj1-1,"..",j1-1,"'"
2508 write(12,'(a8,4i4)') &
2509 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2515 ! finding antiparallel beta
2516 !d write (iout,*) '--------- looking for antiparallel beta ---------'
2521 if (dcont(i).le.rbeta.and. &
2522 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2523 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2524 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2525 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2526 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2530 !d write (iout,*) i1,j1,dcont(i)
2537 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
2538 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2539 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2540 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2541 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2542 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2543 .and. dcont(j).le.rbeta ) goto 6
2547 !d write (iout,*) i1,j1,dcont(j),not_done
2551 if (i1-ii1.gt.1) then
2552 if(lprint)write (iout,*)'antiparallel beta',&
2553 nbeta,ii1-1,i1,jj1,j1-1
2556 bfrag(1,nbfrag)=max0(ii1-1,1)
2559 bfrag(4,nbfrag)=max0(j1-1,1)
2564 isec(ij,1)=isec(ij,1)+1
2565 isec(ij,1+isec(ij,1))=nbeta
2569 isec(ij,1)=isec(ij,1)+1
2570 isec(ij,1+isec(ij,1))=nbeta
2576 if (nstrand.le.9) then
2577 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2578 "DefPropRes 'strand",nstrand,&
2579 "' 'num = ",ii1-2,"..",i1-1,"'"
2581 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2582 "DefPropRes 'strand",nstrand,&
2583 "' 'num = ",ii1-2,"..",i1-1,"'"
2586 if (nstrand.le.9) then
2587 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2588 "DefPropRes 'strand",nstrand,&
2589 "' 'num = ",j1-2,"..",jj1-1,"'"
2591 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2592 "DefPropRes 'strand",nstrand,&
2593 "' 'num = ",j1-2,"..",jj1-1,"'"
2595 write(12,'(a8,4i4)') &
2596 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2602 if (nstrand.gt.0.and.lprint) then
2603 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2606 write(12,'(a9,i1,$)') " | strand",i
2608 write(12,'(a9,i2,$)') " | strand",i
2611 write(12,'(a1)') "'"
2615 ! finding alpha or 310 helix
2621 if (j1.eq.i1+3.and.dcont(i).le.r310 &
2622 .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2623 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2624 !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2627 if (isec(ii1,1).eq.0) then
2636 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2640 !d write (iout,*) i1,j1,not_done
2643 if (j1-ii1.gt.4) then
2645 !d write (iout,*)'helix',nhelix,ii1,j1
2649 hfrag(2,nhfrag)=max0(j1-1,1)
2655 write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2656 if (nhelix.le.9) then
2657 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2658 "DefPropRes 'helix",nhelix,&
2659 "' 'num = ",ii1-1,"..",j1-2,"'"
2661 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2662 "DefPropRes 'helix",nhelix,&
2663 "' 'num = ",ii1-1,"..",j1-2,"'"
2670 if (nhelix.gt.0.and.lprint) then
2671 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2673 if (nhelix.le.9) then
2674 write(12,'(a8,i1,$)') " | helix",i
2676 write(12,'(a8,i2,$)') " | helix",i
2679 write(12,'(a1)') "'"
2683 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2684 write(12,'(a20)') "XMacStand ribbon.mac"
2688 end subroutine secondary
2689 !-----------------------------------------------------------------------------
2690 subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2692 use geometry, only:dist
2694 ! implicit real*8 (a-h,o-z)
2695 ! include 'DIMENSIONS'
2697 ! include 'COMMON.SBRIDGE'
2698 ! include 'COMMON.FFIELD'
2699 ! include 'COMMON.IOUNITS'
2700 ! include 'COMMON.DISTFIT'
2701 ! include 'COMMON.VAR'
2702 ! include 'COMMON.CHAIN'
2703 ! include 'COMMON.MINIM'
2705 character(len=50) :: linia
2707 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2708 real(kind=8) :: time0,time1
2709 integer :: iff(nres),ieval
2710 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2713 integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2714 real(kind=8) :: wstrain0,etot
2716 maxres22=nres*(nres+1)/2
2718 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2719 call var_to_geom(nvar,var)
2726 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2748 call var_to_geom(nvar,var2)
2751 if ( iff(i).eq.1 ) then
2760 !d call write_pdb(3,'combined structure',0d0)
2761 !d time0=MPI_WTIME()
2764 NY=((NRES-4)*(NRES-5))/2
2765 call distfit(.true.,200)
2767 !d time1=MPI_WTIME()
2768 !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2778 call geom_to_var(nvar,var)
2779 !d time0=MPI_WTIME()
2780 call minimize(etot,var,iretcode,nfun)
2781 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2783 !d time1=MPI_WTIME()
2784 !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2785 !d & nfun/(time1-time0),' SOFT eval/s'
2786 call var_to_geom(nvar,var)
2792 if (iff(1).eq.1) then
2798 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2803 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2809 if (iff(nres).eq.1) then
2815 !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2816 !d & "select",ij(1),"-",ij(2),
2817 !d & ",",ij(3),"-",ij(4)
2818 !d call write_pdb(in_pdb,linia,etot)
2824 !d time0=MPI_WTIME()
2825 call minimize(etot,var,iretcode,nfun)
2826 !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2829 !d time1=MPI_WTIME()
2830 !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2831 !d & nfun/(time1-time0),' eval/s'
2832 !d call var_to_geom(nvar,var)
2834 !d call write_pdb(6,'dist structure',etot)
2843 end subroutine contact_cp2
2844 !-----------------------------------------------------------------------------
2845 subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2847 use geometry, only:dist
2849 ! implicit real*8 (a-h,o-z)
2850 ! include 'DIMENSIONS'
2851 ! include 'COMMON.SBRIDGE'
2852 ! include 'COMMON.FFIELD'
2853 ! include 'COMMON.IOUNITS'
2854 ! include 'COMMON.DISTFIT'
2855 ! include 'COMMON.VAR'
2856 ! include 'COMMON.CHAIN'
2857 ! include 'COMMON.MINIM'
2859 character(len=50) :: linia
2861 real(kind=8) :: energy(0:n_ene)
2862 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2863 real(kind=8) :: time0,time1
2864 integer :: iff(nres),ieval
2865 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2869 integer :: in_pdb,i,j,ind,iwsk
2873 if (ieval.eq.-1) debug=.true.
2877 ! store selected dist. constrains from 1st structure
2880 ! Intercept NaNs in the coordinates
2881 ! write(iout,*) (var(i),i=1,nvar)
2886 if (x_sum.ne.x_sum) then
2887 write(iout,*)" *** contact_cp : Found NaN in coordinates"
2889 print *," *** contact_cp : Found NaN in coordinates"
2895 call var_to_geom(nvar,var)
2902 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2925 ! freeze sec.elements from 2nd structure
2933 call var_to_geom(nvar,var2)
2934 call secondary2(debug)
2936 do i=bfrag(1,j),bfrag(2,j)
2941 if (bfrag(3,j).le.bfrag(4,j)) then
2942 do i=bfrag(3,j),bfrag(4,j)
2948 do i=bfrag(4,j),bfrag(3,j)
2956 do i=hfrag(1,j),hfrag(2,j)
2965 ! copy selected res from 1st to 2nd structure
2969 if ( iff(i).eq.1 ) then
2979 ! prepare description in linia variable
2983 if (iff(1).eq.1) then
2989 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2994 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
3000 if (iff(nres).eq.1) then
3005 write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3006 "SELECT",ij(1)-1,"-",ij(2)-1,&
3007 ",",ij(3)-1,"-",ij(4)-1
3013 call contact_cp_min(var,ieval,in_pdb,linia,debug)
3016 end subroutine contact_cp
3017 !-----------------------------------------------------------------------------
3018 subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3022 ! input : theta,phi,alph,omeg,in_pdb,linia,debug
3023 ! output : var,ieval
3025 ! implicit real*8 (a-h,o-z)
3026 ! include 'DIMENSIONS'
3028 ! include 'COMMON.SBRIDGE'
3029 ! include 'COMMON.FFIELD'
3030 ! include 'COMMON.IOUNITS'
3031 ! include 'COMMON.DISTFIT'
3032 ! include 'COMMON.VAR'
3033 ! include 'COMMON.CHAIN'
3034 ! include 'COMMON.MINIM'
3036 character(len=50) :: linia
3038 real(kind=8) :: energy(0:n_ene)
3039 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3040 real(kind=8) :: time0,time1
3041 integer :: ieval,info(3)
3042 logical :: debug,fail,reduce,change !check_var,
3045 integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3047 real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3048 wtor_d01,wstrain0,etot
3050 write(iout,'(a20,i6,a20)') &
3051 '------------------',in_pdb,'-------------------'
3055 call write_pdb(1000+in_pdb,'combined structure',0d0)
3062 ! run optimization of distances
3064 ! uses d0(),w() and mask() for frozen 2D
3066 !test---------------------------------------------
3068 !test NY=((NRES-4)*(NRES-5))/2
3069 !test call distfit(debug,5000)
3101 call geom_to_var(nvar,var)
3102 !de change=reduce(var)
3103 if (check_var(var,info)) then
3104 write(iout,*) 'cp_min error in input'
3105 print *,'cp_min error in input'
3109 !d call etotal(energy(0))
3110 !d call enerprint(energy(0))
3114 !dtest call minimize(etot,var,iretcode,nfun)
3115 !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3118 !d call etotal(energy(0))
3119 !d call enerprint(energy(0))
3138 !test--------------------------------------------------
3144 write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3145 call write_pdb(2000+in_pdb,'distfit structure',0d0)
3153 ! run soft pot. optimization
3155 ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3157 ! mask_phi(),mask_theta(),mask_side(),mask_r
3163 !de change=reduce(var)
3164 !de if (check_var(var,info)) write(iout,*) 'error before soft'
3168 call minimize(etot,var,iretcode,nfun)
3170 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3174 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3175 nfun/(time1-time0),' SOFT eval/s'
3178 call var_to_geom(nvar,var)
3180 call write_pdb(3000+in_pdb,'soft structure',etot)
3184 ! run full UNRES optimization with constrains and frozen 2D
3185 ! the same variables as soft pot. optimizatio
3191 ! check overlaps before calling full UNRES minim
3193 call var_to_geom(nvar,var)
3197 write(iout,*) 'N7 ',energy(0)
3198 if (energy(0).ne.energy(0)) then
3199 write(iout,*) 'N7 error - gives NaN',energy(0)
3203 if (energy(1).eq.1.0d20) then
3204 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3205 call overlap_sc(fail)
3209 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3221 !dte time0=MPI_WTIME()
3222 !de change=reduce(var)
3223 !de if (check_var(var,info)) then
3224 !de write(iout,*) 'error before mask dist'
3225 !de call var_to_geom(nvar,var)
3227 !de call write_pdb(10000+in_pdb,'before mask dist',etot)
3229 !dte call minimize(etot,var,iretcode,nfun)
3230 !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3231 !dte & ' eval ',nfun
3232 !dte ieval=ieval+nfun
3234 !dte time1=MPI_WTIME()
3235 !dte write (iout,'(a,f6.2,f8.2,a)')
3236 !dte & ' Time for mask dist min.',time1-time0,
3237 !dte & nfun/(time1-time0),' eval/s'
3238 !dte call flush(iout)
3241 call var_to_geom(nvar,var)
3243 call write_pdb(4000+in_pdb,'mask dist',etot)
3246 ! switch off freezing of 2D and
3247 ! run full UNRES optimization with constrains
3253 !de change=reduce(var)
3254 !de if (check_var(var,info)) then
3255 !de write(iout,*) 'error before dist'
3256 !de call var_to_geom(nvar,var)
3258 !de call write_pdb(11000+in_pdb,'before dist',etot)
3261 call minimize(etot,var,iretcode,nfun)
3263 !de change=reduce(var)
3264 !de if (check_var(var,info)) then
3265 !de write(iout,*) 'error after dist',ico
3266 !de call var_to_geom(nvar,var)
3268 !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3270 write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3276 write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3277 nfun/(time1-time0),' eval/s'
3279 !de call etotal(energy(0))
3280 !de write(iout,*) 'N7 after dist',energy(0)
3284 call var_to_geom(nvar,var)
3286 call write_pdb(in_pdb,linia,etot)
3298 end subroutine contact_cp_min
3299 !-----------------------------------------------------------------------------
3302 use geometry, only:dist
3304 ! implicit real*8 (a-h,o-z)
3305 ! include 'DIMENSIONS'
3307 ! include 'COMMON.GEO'
3308 ! include 'COMMON.CHAIN'
3309 ! include 'COMMON.IOUNITS'
3310 ! include 'COMMON.VAR'
3311 ! include 'COMMON.CONTROL'
3312 ! include 'COMMON.SBRIDGE'
3313 ! include 'COMMON.FFIELD'
3314 ! include 'COMMON.MINIM'
3315 ! include 'COMMON.INTERACT'
3317 ! include 'COMMON.DISTFIT'
3318 integer :: iff(nres)
3319 real(kind=8) :: time0,time1
3320 real(kind=8) :: energy(0:n_ene),ee
3321 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3323 logical :: debug,ltest,fail
3324 character(len=50) :: linia
3325 integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3327 real(kind=8) :: wstrain0,wang0,etot
3333 !------------------------
3335 ! freeze sec.elements
3345 do i=bfrag(1,j),bfrag(2,j)
3350 if (bfrag(3,j).le.bfrag(4,j)) then
3351 do i=bfrag(3,j),bfrag(4,j)
3357 do i=bfrag(4,j),bfrag(3,j)
3365 do i=hfrag(1,j),hfrag(2,j)
3377 ! store dist. constrains
3381 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3386 dhpb(nhpb)=DIST(i,j)
3394 call write_pdb(100+in_pdb,'input reg. structure',0d0)
3404 ! run soft pot. optimization
3410 call geom_to_var(nvar,var)
3414 call minimize(etot,var,iretcode,nfun)
3416 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3420 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3421 nfun/(time1-time0),' SOFT eval/s'
3423 call var_to_geom(nvar,var)
3425 call write_pdb(300+in_pdb,'soft structure',etot)
3428 ! run full UNRES optimization with constrains and frozen 2D
3429 ! the same variables as soft pot. optimizatio
3438 call minimize(etot,var,iretcode,nfun)
3439 write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3446 write (iout,'(a,f6.2,f8.2,a)') &
3447 ' Time for mask dist min.',time1-time0,&
3448 nfun/(time1-time0),' eval/s'
3450 call var_to_geom(nvar,var)
3452 call write_pdb(400+in_pdb,'mask & dist',etot)
3455 ! switch off constrains and
3456 ! run full UNRES optimization with frozen 2D
3471 call minimize(etot,var,iretcode,nfun)
3472 write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3478 write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,&
3479 nfun/(time1-time0),' eval/s'
3483 call var_to_geom(nvar,var)
3485 call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3492 ! run full UNRES optimization with constrains and NO frozen 2D
3502 wstrain=wstrain0/ico
3507 call minimize(etot,var,iretcode,nfun)
3508 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3509 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3516 write (iout,'(a,f6.2,f8.2,a)') &
3517 ' Time for dist min.',time1-time0,&
3518 nfun/(time1-time0),' eval/s'
3520 call var_to_geom(nvar,var)
3522 call write_pdb(600+in_pdb+ico,'dist cons',etot)
3539 call minimize(etot,var,iretcode,nfun)
3540 write(iout,*)'------------------------------------------------'
3541 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3547 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
3548 nfun/(time1-time0),' eval/s'
3551 call var_to_geom(nvar,var)
3553 call write_pdb(999,'full min',etot)
3557 end subroutine softreg
3558 !-----------------------------------------------------------------------------
3559 subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3561 use geometry, only:dist
3563 ! implicit real*8 (a-h,o-z)
3564 ! include 'DIMENSIONS'
3566 ! include 'COMMON.GEO'
3567 ! include 'COMMON.VAR'
3568 ! include 'COMMON.INTERACT'
3569 ! include 'COMMON.IOUNITS'
3570 ! include 'COMMON.DISTFIT'
3571 ! include 'COMMON.SBRIDGE'
3572 ! include 'COMMON.CONTROL'
3573 ! include 'COMMON.FFIELD'
3574 ! include 'COMMON.MINIM'
3575 ! include 'COMMON.CHAIN'
3576 real(kind=8) :: time0,time1
3577 real(kind=8) :: energy(0:n_ene),ee
3578 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3579 integer :: jdata(5),isec(nres)
3582 integer :: i1,i2,i3,i4,i5,ieval,ij
3583 integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3584 real(kind=8) :: etot,wscloc0,wstrain0
3592 call secondary2(.false.)
3598 do i=bfrag(1,j),bfrag(2,j)
3601 do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3606 do i=hfrag(1,j),hfrag(2,j)
3612 ! cut strands at the ends
3614 if (jdata(2)-jdata(1).gt.3) then
3617 if (jdata(3).lt.jdata(4)) then
3627 !v call etotal(energy(0))
3629 !v write(iout,*) nnt,nct,etot
3630 !v call write_pdb(ij*100,'first structure',etot)
3631 !v write(iout,*) 'N16 test',(jdata(i),i=1,5)
3633 !------------------------
3634 ! generate constrains
3637 if(ishift.eq.0) ishift=-2
3640 do i=jdata(1),jdata(2)
3642 if(jdata(4).gt.jdata(3))then
3643 do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3645 !d print *,i,j,j+ishift
3650 dhpb(nhpb)=DIST(i,j+ishift)
3653 do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3655 !d print *,i,j,j+ishift
3660 dhpb(nhpb)=DIST(i,j+ishift)
3667 if(isec(i).gt.0.or.isec(j).gt.0) then
3673 dhpb(nhpb)=DIST(i,j)
3680 call geom_to_var(nvar,var)
3687 wstrain=wstrain0/ico
3689 !v time0=MPI_WTIME()
3690 call minimize(etot,var,iretcode,nfun)
3691 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3692 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3695 !v time1=MPI_WTIME()
3696 !v write (iout,'(a,f6.2,f8.2,a)')
3697 !v & ' Time for dist min.',time1-time0,
3698 !v & nfun/(time1-time0),' eval/s'
3699 !v call var_to_geom(nvar,var)
3701 !v call write_pdb(ij*100+ico,'dist cons',etot)
3713 call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3716 !v call etotal(energy(0))
3718 !v call write_pdb(ij*100+10,'sc_move',etot)
3720 !d print *,nft_sc,etot
3723 end subroutine beta_slide
3724 !-----------------------------------------------------------------------------
3725 subroutine beta_zip(i1,i2,ieval,ij)
3728 ! implicit real*8 (a-h,o-z)
3729 ! include 'DIMENSIONS'
3731 ! include 'COMMON.GEO'
3732 ! include 'COMMON.VAR'
3733 ! include 'COMMON.INTERACT'
3734 ! include 'COMMON.IOUNITS'
3735 ! include 'COMMON.DISTFIT'
3736 ! include 'COMMON.SBRIDGE'
3737 ! include 'COMMON.CONTROL'
3738 ! include 'COMMON.FFIELD'
3739 ! include 'COMMON.MINIM'
3740 ! include 'COMMON.CHAIN'
3741 real(kind=8) :: time0,time1
3742 real(kind=8) :: energy(0:n_ene),ee
3743 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3744 character(len=10) :: test
3746 integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3747 real(kind=8) :: etot,wstrain0
3749 !v call etotal(energy(0))
3751 !v write(test,'(2i5)') i1,i2
3752 !v call write_pdb(ij*100,test,etot)
3753 !v write(iout,*) 'N17 test',i1,i2,etot,ij
3756 ! generate constrains
3767 call geom_to_var(nvar,var)
3773 wstrain=wstrain0/ico
3774 !v time0=MPI_WTIME()
3775 call minimize(etot,var,iretcode,nfun)
3776 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3777 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3780 !v time1=MPI_WTIME()
3781 !v write (iout,'(a,f6.2,f8.2,a)')
3782 !v & ' Time for dist min.',time1-time0,
3783 !v & nfun/(time1-time0),' eval/s'
3784 ! do not comment the next line
3785 call var_to_geom(nvar,var)
3787 !v call write_pdb(ij*100+ico,'dist cons',etot)
3795 !v call etotal(energy(0))
3797 !v write(iout,*) 'N17 test end',i1,i2,etot,ij
3800 end subroutine beta_zip
3801 !-----------------------------------------------------------------------------
3803 !-----------------------------------------------------------------------------
3804 subroutine thread_seq
3806 use geometry, only:dist
3807 use random, only:iran_num
3808 use control, only:tcpu
3809 use regularize_, only:regularize
3810 use mcm_data, only: nsave_part,nacc_tot
3811 ! Thread the sequence through a database of known structures
3812 ! implicit real*8 (a-h,o-z)
3813 ! include 'DIMENSIONS'
3814 use MPI_data !include 'COMMON.INFO'
3819 ! include 'COMMON.CONTROL'
3820 ! include 'COMMON.CHAIN'
3821 ! include 'COMMON.DBASE'
3822 ! include 'COMMON.INTERACT'
3823 ! include 'COMMON.VAR'
3824 ! include 'COMMON.THREAD'
3825 ! include 'COMMON.FFIELD'
3826 ! include 'COMMON.SBRIDGE'
3827 ! include 'COMMON.HEADER'
3828 ! include 'COMMON.IOUNITS'
3829 ! include 'COMMON.TIME1'
3830 ! include 'COMMON.CONTACTS'
3831 ! include 'COMMON.MCM'
3832 ! include 'COMMON.NAMES'
3834 integer :: ThreadId,ThreadType,Kwita
3836 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3837 real(kind=8) :: przes(3),obr(3,3)
3838 real(kind=8) :: time_for_thread
3839 logical :: found_pattern,non_conv
3840 character(len=32) :: head_pdb
3841 real(kind=8) :: energia(0:n_ene)
3842 integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3844 real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3846 n_ene_comp=nprint_ene
3851 if (me.eq.king) then
3870 ave_time_for_thread=0.0D0
3871 max_time_for_thread=0.0D0
3872 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3873 nthread=nexcl+nthread
3874 do ithread=1,nthread
3875 found_pattern=.false.
3877 do while (.not.found_pattern)
3879 if (itrial.gt.1000) then
3880 write (iout,'(/a/)') 'Too many attempts to find pattern.'
3883 call recv_stop_sig(Kwita)
3884 call send_stop_sig(-3)
3888 ! Find long enough chain in the database
3890 nres_t=nres_base(1,ii)
3891 ! Select the starting position to thread.
3892 print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3893 nres_t,' nres0=',nres0
3894 if (nres_t.ge.nres0) then
3895 ist=iran_num(0,nres_t-nres0)
3897 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3898 if (Kwita.lt.0) then
3899 write (iout,*) 'Stop signal received. Terminating.'
3900 write (*,*) 'Stop signal received. Terminating.'
3902 write (*,*) 'ithread=',ithread,' nthread=',nthread
3905 call pattern_receive
3908 if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3910 found_pattern=.true.
3912 ! If this point is reached, the pattern has not yet been examined.
3914 ! print *,'found_pattern:',found_pattern
3920 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3921 if (Kwita.lt.0) then
3922 write (iout,*) 'Stop signal received. Terminating.'
3924 write (*,*) 'ithread=',ithread,' nthread=',nthread
3930 ipatt(2,ithread)=ist
3932 write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3933 'Processor:',me,' Attempt:',ithread,&
3934 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3935 ' start at res.',ist+1
3936 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3937 ' Attempt:',ithread,&
3938 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3939 ' start at res.',ist+1
3941 write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3942 'Attempt:',ithread,&
3943 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3944 ' start at res.',ist+1
3945 write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3946 'Attempt:',ithread,&
3947 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3948 ' start at res.',ist+1
3951 ! Copy coordinates from the database.
3955 c(j,i)=cart_base(j,i+ist,ii)
3958 !d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
3960 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3962 !d write (iout,'(a,f10.5)')
3963 !d & 'Initial RMS deviation from reference structure:',rms
3964 if (itype(nres,1).eq.ntyp1) then
3966 dcj=c(j,nres-2)-c(j,nres-3)
3967 c(j,nres)=c(j,nres-1)+dcj
3968 c(j,2*nres)=c(j,nres)
3971 if (itype(1,1).eq.ntyp1) then
3978 call int_from_cart(.false.,.false.)
3979 !d print *,'Exit INT_FROM_CART.'
3980 !d print *,'nhpb=',nhpb
3985 ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3987 ! stop 'End generate'
3988 ! Generate SC conformations.
3992 !d print *,'Processor:',me,': exit GEN_SIDE.'
3994 !d print *,'Exit GEN_SIDE.'
3996 ! Calculate initial energy.
3998 call etotal(energia)
4001 ener0(i,ithread)=energia(i)
4003 ener0(n_ene_comp+1,ithread)=energia(0)
4005 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4006 ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
4008 ener0(n_ene_comp+2,ithread)=rms
4009 ener0(n_ene_comp+4,ithread)=frac
4010 ener0(n_ene_comp+5,ithread)=frac_nn
4012 ener0(n_ene_comp+3,ithread)=0.0d0
4015 print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
4017 print*,'ithread=',ithread,' Start REGULARIZE.'
4020 call regularize(nct-nnt+1,etot,rms,&
4021 cart_base(1,ist+nnt,ipattern),iretcode)
4023 time_for_thread=curr_tim1-curr_tim
4024 ave_time_for_thread= &
4025 ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4026 if (time_for_thread.gt.max_time_for_thread) &
4027 max_time_for_thread=time_for_thread
4029 print *,'Processor',me,': Exit REGULARIZE.'
4030 if (WhatsUp.eq.2) then
4032 'Sufficient number of confs. collected. Terminating.'
4035 else if (WhatsUp.eq.-1) then
4037 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4038 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4039 call send_stop_sig(-2)
4041 else if (WhatsUp.eq.-2) then
4043 write (iout,*) 'Timeup signal received. Terminating.'
4045 else if (WhatsUp.eq.-3) then
4047 write (iout,*) 'Error stop signal received. Terminating.'
4051 print *,'Exit REGULARIZE.'
4052 if (iretcode.eq.11) then
4053 write (iout,'(/a/)') &
4054 '******* Allocated time exceeded in SUMSL. The program will stop.'
4059 head_pdb=titel(:24)//':'//str_nam(ipattern)
4060 if (outpdb) call pdbout(etot,head_pdb,ipdb)
4061 if (outmol2) call mol2out(etot,head_pdb)
4063 call briefout(ithread,etot)
4065 link_end=min0(link_end,nss)
4066 write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4068 call etotal(energia)
4069 ! call enerprint(energia(0))
4072 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4073 !d write (iout,'(a,f10.5)')
4074 !d & 'RMS deviation from reference structure:',dsqrt(rms)
4076 ener(i,ithread)=energia(i)
4078 ener(n_ene_comp+1,ithread)=energia(0)
4079 ener(n_ene_comp+3,ithread)=rms
4081 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4082 ener(n_ene_comp+2,ithread)=rms
4083 ener(n_ene_comp+4,ithread)=frac
4084 ener(n_ene_comp+5,ithread)=frac_nn
4086 call write_stat_thread(ithread,ipattern,ist)
4087 ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4088 ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4089 ! & (ener(k,ithread),k=12,14)
4091 if (me.eq.king) then
4093 call pattern_receive
4094 call receive_MCM_info
4095 if (nacc_tot.ge.nthread) then
4097 'Sufficient number of conformations collected nacc_tot=',&
4098 nacc_tot,'. Stopping other processors and terminating.'
4100 'Sufficient number of conformations collected nacc_tot=',&
4101 nacc_tot,'. Stopping other processors and terminating.'
4102 call recv_stop_sig(Kwita)
4103 if (Kwita.eq.0) call send_stop_sig(-1)
4108 call send_MCM_info(2)
4111 if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4112 write (iout,'(/2a)') &
4113 '********** There would be not enough time for another thread. ',&
4114 'The program will stop.'
4116 '********** There would be not enough time for another thread. ',&
4117 'The program will stop.'
4118 write (iout,'(a,1pe14.4/)') &
4119 'Elapsed time for last threading step: ',time_for_thread
4122 call recv_stop_sig(Kwita)
4123 call send_stop_sig(-2)
4128 write (iout,'(a,1pe14.4)') &
4129 'Elapsed time for this threading step: ',time_for_thread
4132 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4133 if (Kwita.lt.0) then
4134 write (iout,*) 'Stop signal received. Terminating.'
4135 write (*,*) 'Stop signal received. Terminating.'
4137 write (*,*) 'nthread=',nthread,' ithread=',ithread
4143 call send_stop_sig(-1)
4147 ! Any messages left for me?
4148 call pattern_receive
4149 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4151 call write_thread_summary
4153 if (king.eq.king) then
4155 do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4158 call recv_stop_sig(Kwita)
4159 call receive_MCM_info
4162 call receive_thread_results(iproc)
4164 call write_thread_summary
4166 call send_thread_results
4170 end subroutine thread_seq
4171 !-----------------------------------------------------------------------------
4174 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4176 use random, only:iran_num
4177 ! implicit real*8 (a-h,o-z)
4178 ! include 'DIMENSIONS'
4179 ! include 'COMMON.CHAIN'
4180 ! include 'COMMON.DBASE'
4181 ! include 'COMMON.INTERACT'
4182 ! include 'COMMON.VAR'
4183 ! include 'COMMON.THREAD'
4184 ! include 'COMMON.FFIELD'
4185 ! include 'COMMON.SBRIDGE'
4186 ! include 'COMMON.HEADER'
4187 ! include 'COMMON.GEO'
4188 ! include 'COMMON.IOUNITS'
4189 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
4190 !el integer :: icall
4191 !el common /srutu/ icall
4192 real(kind=8) :: energia(0:n_ene)
4193 logical :: glycine,fail
4194 integer :: i,maxsample,link_end0,ind_sc,isample
4195 real(kind=8) :: alph0,omeg0,e1,e0
4199 link_end=min0(link_end,nss)
4201 if (itype(i,1).ne.10) then
4202 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
4203 call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1)
4207 call etotal(energia)
4209 do isample=1,maxsample
4210 ! Choose a non-glycine side chain.
4213 ind_sc=iran_num(nnt,nct)
4214 glycine=(itype(ind_sc,1).eq.10)
4218 call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
4219 omeg(ind_sc),fail,1)
4221 call etotal(energia)
4222 !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4223 !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4224 !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4235 end subroutine sc_conf
4236 !-----------------------------------------------------------------------------
4238 !-----------------------------------------------------------------------------
4239 logical function check_var(var,info)
4243 ! implicit real*8 (a-h,o-z)
4244 ! include 'DIMENSIONS'
4245 ! include 'COMMON.VAR'
4246 ! include 'COMMON.IOUNITS'
4247 ! include 'COMMON.GEO'
4248 ! include 'COMMON.SETUP'
4249 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4250 integer,dimension(3) :: info
4254 do i=nphi+ntheta+1,nphi+ntheta+nside
4255 ! Check the side chain "valence" angles alpha
4256 if (var(i).lt.1.0d-7) then
4257 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4258 write (iout,*) 'Processor',me,'received bad variables!!!!'
4259 write (iout,*) 'Variables'
4260 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4261 write (iout,*) 'Continuing calculations at this point',&
4262 ' could destroy the results obtained so far... ABORTING!!!!!!'
4263 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4264 'valence angle alpha',i-nphi-ntheta,var(i),&
4265 'n it',info(1),info(2),'mv ',info(3)
4266 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4267 write (*,*) 'Processor',me,'received bad variables!!!!'
4268 write (*,*) 'Variables'
4269 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4270 write (*,*) 'Continuing calculations at this point',&
4271 ' could destroy the results obtained so far... ABORTING!!!!!!'
4272 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4273 'valence angle alpha',i-nphi-ntheta,var(i),&
4274 'n it',info(1),info(2),'mv ',info(3)
4279 ! Check the backbone "valence" angles theta
4280 do i=nphi+1,nphi+ntheta
4281 if (var(i).lt.1.0d-7) then
4282 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4283 write (iout,*) 'Processor',me,'received bad variables!!!!'
4284 write (iout,*) 'Variables'
4285 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4286 write (iout,*) 'Continuing calculations at this point',&
4287 ' could destroy the results obtained so far... ABORTING!!!!!!'
4288 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4289 'valence angle theta',i-nphi,var(i),&
4290 'n it',info(1),info(2),'mv ',info(3)
4291 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4292 write (*,*) 'Processor',me,'received bad variables!!!!'
4293 write (*,*) 'Variables'
4294 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4295 write (*,*) 'Continuing calculations at this point',&
4296 ' could destroy the results obtained so far... ABORTING!!!!!!'
4297 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4298 'valence angle theta',i-nphi,var(i),&
4299 'n it',info(1),info(2),'mv ',info(3)
4305 end function check_var
4306 !-----------------------------------------------------------------------------
4308 !-----------------------------------------------------------------------------
4309 subroutine distfit(debug,maxit)
4311 use geometry_data, only: phi
4314 ! implicit real*8 (a-h,o-z)
4315 ! include 'DIMENSIONS'
4316 ! include 'COMMON.CHAIN'
4317 ! include 'COMMON.VAR'
4318 ! include 'COMMON.IOUNITS'
4319 ! include 'COMMON.DISTFIT'
4320 integer :: i,maxit,MAXMAR,IT,IMAR
4321 real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres)
4322 logical :: debug,sing
4323 real(kind=8) :: TOL,RL,F0,AIN,F1
4325 !input------------------------------------
4327 ! NY=((NRES-4)*(NRES-5))/2
4328 !input------------------------------------
4334 CALL TRANSFER(NRES,phi,phiold)
4338 !d WRITE (IOUT,*) 'DISTFIT: F0=',F0
4354 CALL TRANSFER(NX,XX,X)
4355 CALL BANACH(NX,NRES,H,X,sing)
4360 IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4362 WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4363 WRITE (IOUT,*) 'IT=',it,'F=',F0
4368 phi(I)=phiold(I)+mask(i)*X(I-3)
4373 !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4375 CALL TRANSFER(NRES,phi,phiold)
4378 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4380 WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4381 WRITE (IOUT,*) 'IT=',it,'F=',F1
4387 WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4388 WRITE (IOUT,*) 'IT=',it,'F=',F0
4389 CALL TRANSFER(NRES,phiold,phi)
4392 !d write (iout,*) "it",it," imar",imar," f0",f0
4394 WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4396 end subroutine distfit
4397 !-----------------------------------------------------------------------------
4398 real(kind=8) function RDIF()
4401 use geometry, only: dist
4402 ! implicit real*8 (a-h,o-z)
4403 ! include 'DIMENSIONS'
4404 ! include 'COMMON.CHAIN'
4405 ! include 'COMMON.DISTFIT'
4407 real(kind=8) :: suma,DIJ
4416 if (w(ind).ne.0.0) then
4418 suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4420 ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4428 !-----------------------------------------------------------------------------
4433 use geometry, only:dist
4434 ! implicit real*8 (a-h,o-z)
4435 ! include 'DIMENSIONS'
4436 ! include 'COMMON.CHAIN'
4437 ! include 'COMMON.DISTFIT'
4438 ! include 'COMMON.GEO'
4439 integer :: i,j,k,l,I1,I2,IND
4440 real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4453 R13(K)=C(K,J)-C(K,I1)
4457 R24(L)=C(L,K)-C(L,I2)
4459 IND=((J-1)*(2*NRES-J-6))/2+K-3
4460 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4461 PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4462 PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4463 DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4468 end subroutine RDERIV
4469 !-----------------------------------------------------------------------------
4473 ! implicit real*8 (a-h,o-z)
4474 ! include 'DIMENSIONS'
4475 ! include 'COMMON.CHAIN'
4476 ! include 'COMMON.DISTFIT'
4478 real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4486 XI=XI+BKIWK*(D0(K)-DDD(K))
4494 HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4501 end subroutine HEVAL
4502 !-----------------------------------------------------------------------------
4503 subroutine VEC(I,J,U)
4505 use geometry_data, only: C
4506 ! Find the unit vector from atom (I) to atom (J). Store in U.
4508 ! implicit real*8 (a-h,o-z)
4509 ! include 'DIMENSIONS'
4510 ! include 'COMMON.CHAIN'
4512 real(kind=8),DIMENSION(3) :: U
4513 real(kind=8) :: ANORM,UK
4527 !-----------------------------------------------------------------------------
4528 subroutine TRANSFER(N,X1,X2)
4530 ! implicit real*8 (a-h,o-z)
4531 ! include 'DIMENSIONS'
4533 real(kind=8),DIMENSION(N) :: X1,X2
4537 end subroutine TRANSFER
4538 !-----------------------------------------------------------------------------
4539 !-----------------------------------------------------------------------------
4540 subroutine alloc_compare_arrays
4542 maxres22=nres*(nres+1)/2
4544 ! common /struct/ in io_common: read_threadbase
4545 ! allocate(cart_base !(3,maxres_base,maxseq)
4546 ! allocate(nres_base !(3,maxseq)
4547 ! allocate(str_nam !(maxseq)
4549 ! COMMON /c_frag/ in io_conf: readpdb
4550 if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4551 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4554 if(.not.allocated(w)) allocate(w(maxres22),d0(maxres22)) !(maxres22)
4556 !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4557 if (.not.allocated(DDD)) allocate(DDD(maxres22)) !(maxres22)
4558 if (.not.allocated(H)) allocate(H(nres,nres)) !(MAXRES,MAXRES)
4560 if (.not.allocated(XX)) allocate(XX(nres)) !(MAXRES)
4562 if (.not.allocated(mask)) allocate(mask(nres)) !(maxres)
4565 if (.not.allocated(iexam))allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4567 if (.not.allocated(ener0)) allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4570 end subroutine alloc_compare_arrays
4571 !-----------------------------------------------------------------------------
4573 !-----------------------------------------------------------------------------