2 !-----------------------------------------------------------------------------
8 #if .not. defined WHAM_RUN && .not. defined CLUSTER
14 use control, only: hpb_partition
16 use minimm, only: sc_move, minimize
19 !-----------------------------------------------------------------------------
22 !-----------------------------------------------------------------------------
24 #if .not. defined WHAM_RUN && .not. defined CLUSTER
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine contact(lprint,ncont,icont,co)
30 use geometry, only:dist
31 ! implicit real*8 (a-h,o-z)
32 ! include 'DIMENSIONS'
33 ! include 'COMMON.IOUNITS'
34 ! include 'COMMON.CHAIN'
35 ! include 'COMMON.INTERACT'
36 ! include 'COMMON.FFIELD'
37 ! include 'COMMON.NAMES'
38 real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
40 integer,dimension(2,12*nres) :: icont!(2,12*nres) !(2,maxcont) (maxcont=12*maxres)
43 real(kind=8) :: co,rcomp
44 integer :: kkk,i,j,i1,i2,it1,it2,iti,itj,inum,jnum
49 iti=iabs(itype(i,molnum(i)))
50 if (molnum(i).lt.3) then
57 itj=iabs(itype(j,molnum(i)))
58 if (molnum(j).lt.3) then
64 ! rcomp=sigmaii(iti,itj)+1.0D0
65 rcomp=facont*sigmaii(iti,itj)
67 ! rcomp=sigma(iti,itj)+1.0D0
68 rcomp=facont*sigma(iti,itj)
71 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
72 if (dist(inum,jnum).lt.rcomp) then
80 write (iout,'(a)') 'Contact map:'
86 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
87 i,restyp(it1,1),i1,restyp(it2,1),i2
92 co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
94 co = co / (nres*ncont)
96 end subroutine contact
97 !-----------------------------------------------------------------------------
98 real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref)
100 ! implicit real*8 (a-h,o-z)
101 ! include 'DIMENSIONS'
102 ! include 'COMMON.IOUNITS'
103 integer :: ncont,ncont_ref
104 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
106 integer :: i,j,nmatch
108 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
109 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
110 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
111 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
112 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
115 if (icont(1,i).eq.icont_ref(1,j) .and. &
116 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
119 ! print *,' nmatch=',nmatch
120 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
121 contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
123 end function contact_fract
124 !-----------------------------------------------------------------------------
125 real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
127 ! implicit real*8 (a-h,o-z)
128 ! include 'DIMENSIONS'
129 ! include 'COMMON.IOUNITS'
130 integer :: ncont,ncont_ref
131 integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres)
133 integer :: i,j,nmatch
135 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref
136 ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
137 ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
138 ! write (iout,'(20i4)') (icont(1,i),i=1,ncont)
139 ! write (iout,'(20i4)') (icont(2,i),i=1,ncont)
142 if (icont(1,i).eq.icont_ref(1,j) .and. &
143 icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
146 ! print *,' nmatch=',nmatch
147 ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
148 contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
150 end function contact_fract_nn
151 !-----------------------------------------------------------------------------
152 subroutine hairpin(lprint,nharp,iharp)
154 use geometry, only:dist
155 ! implicit real*8 (a-h,o-z)
156 ! include 'DIMENSIONS'
157 ! include 'COMMON.IOUNITS'
158 ! include 'COMMON.CHAIN'
159 ! include 'COMMON.INTERACT'
160 ! include 'COMMON.FFIELD'
161 ! include 'COMMON.NAMES'
163 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
165 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
166 logical :: lprint,not_done
167 real(kind=8) :: rcomp=6.0d0
169 integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1
170 ! allocate(icont(2,12*nres))
174 ! print *,'nnt=',nnt,' nct=',nct
177 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
181 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
183 if (dist(2*nres+1,2*nres+2).lt.rcomp) then
191 write (iout,'(a)') 'PP contact map:'
197 write (iout,'(i3,2x,a,i4,2x,a,i4)') &
198 i,restyp(it1,1),i1,restyp(it2,1),i2
206 if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
207 ! write (iout,*) "found turn at ",i1,j1
215 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
219 ! write (iout,*) i1,j1,not_done
229 ! write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
234 ! write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
237 write (iout,*) "Hairpins:"
244 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
245 write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
247 ! write (iout,'(a,i3,$)') restyp(itype(k,1)),k
252 end subroutine hairpin
253 !-----------------------------------------------------------------------------
255 !-----------------------------------------------------------------------------
256 subroutine elecont(lprint,ncont,icont)
258 ! implicit real*8 (a-h,o-z)
259 ! include 'DIMENSIONS'
260 ! include 'COMMON.IOUNITS'
261 ! include 'COMMON.CHAIN'
262 ! include 'COMMON.INTERACT'
263 ! include 'COMMON.LOCAL'
264 ! include 'COMMON.FFIELD'
265 ! include 'COMMON.NAMES'
267 real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_
268 real(kind=8) :: ael6_i,ael3_i
269 real(kind=8),dimension(2,2) :: app_,bpp_,rpp_
271 integer,dimension(2,12*nres) :: icont !(2,12*nres)(2,maxcont) (maxcont=12*maxres)
272 real(kind=8),dimension(12*nres) :: econt !(maxcont)
274 integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2
275 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw
276 real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb
277 real(kind=8) :: xmedi,ymedi,zmedi
278 real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij
279 real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,&
280 evdwij,el1,el2,eesij,ene
282 ! Load the constants of peptide bond - peptide bond interactions.
283 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
284 ! proline) - determined by averaging ECEPP energy.
288 ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
289 data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
290 data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
291 data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
293 !el allocate(econt(12*nres)) !(maxcont)
296 elecutoff_14 = -0.5d0
297 if (lprint) write (iout,'(a)') &
298 "Constants of electrostatic interaction energy expression."
302 app_(i,j)=epp(i,j)*rri*rri
303 bpp_(i,j)=-2.0*epp(i,j)*rri
304 ael6_(i,j)=elpp_6(i,j)*4.2**6
305 ael3_(i,j)=elpp_3(i,j)*4.2**3
307 write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),&
314 ! print *, "nntt,nct",nnt,nct-2
316 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
327 if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
330 if (j.eq.i+2 .and. itelj.eq.2) iteli=2
331 if (iteli.eq.2 .and. itelj.eq.2) goto 4
332 aaa=app_(iteli,itelj)
333 bbb=bpp_(iteli,itelj)
334 ael6_i=ael6_(iteli,itelj)
335 ael3_i=ael3_(iteli,itelj)
339 xj=c(1,j)+0.5*dxj-xmedi
340 yj=c(2,j)+0.5*dyj-ymedi
341 zj=c(3,j)+0.5*dzj-zmedi
342 rrmij=1.0/(xj*xj+yj*yj+zj*zj)
347 cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
348 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
349 cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
350 fac=cosa-3.0*cosb*cosg
356 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
359 if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
360 j.eq.i+2 .and. eesij.le.elecutoff_14) then
371 write (iout,*) 'Total average electrostatic energy: ',ees
372 write (iout,*) 'VDW energy between peptide-group centers: ',evdw
374 write (iout,*) 'Electrostatic contacts before pruning: '
380 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
381 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
384 ! For given residues keep only the contacts with the greatest energy.
386 do while (i.lt.ncont)
392 do while (j.lt.ncont)
394 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
395 ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
396 ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
397 ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
398 if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
399 if (ic1.eq.icont(1,j)) then
401 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
402 .and. iabs(icont(1,k)-ic1).le.2 .and. &
403 econt(k).lt.econt(j) ) goto 21
405 else if (ic2.eq.icont(2,j) ) then
407 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
408 .and. iabs(icont(2,k)-ic2).le.2 .and. &
409 econt(k).lt.econt(j) ) goto 21
414 icont(1,k-1)=icont(1,k)
415 icont(2,k-1)=icont(2,k)
420 ! write (iout,*) "ncont",ncont
422 ! write (iout,*) icont(1,k),icont(2,k)
425 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
427 if (ic1.eq.icont(1,j)) then
429 if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
430 .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
431 econt(k).lt.econt(i) ) goto 21
433 else if (ic2.eq.icont(2,j) ) then
435 if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
436 .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
437 econt(k).lt.econt(i) ) goto 21
442 icont(1,k-1)=icont(1,k)
443 icont(2,k-1)=icont(2,k)
447 ! write (iout,*) "ncont",ncont
449 ! write (iout,*) icont(1,k),icont(2,k)
460 write (iout,*) 'Electrostatic contacts after pruning: '
466 write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
467 i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
471 end subroutine elecont
472 !-----------------------------------------------------------------------------
473 subroutine secondary2(lprint)
475 ! implicit real*8 (a-h,o-z)
476 ! include 'DIMENSIONS'
477 ! include 'COMMON.CHAIN'
478 ! include 'COMMON.IOUNITS'
479 ! include 'COMMON.DISTFIT'
480 ! include 'COMMON.VAR'
481 ! include 'COMMON.GEO'
482 ! include 'COMMON.CONTROL'
483 integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
485 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres)
486 integer,dimension(nres,0:4) :: isec !(maxres,4)
487 integer,dimension(nres) :: nsec !(maxres)
488 logical :: lprint,not_done !,freeres
489 real(kind=8) :: p1,p2
492 !el allocate(icont(2,12*nres),isec(nres,4),nsec(nres))
494 if(.not.dccart) call chainbuild_cart
495 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
496 !d call write_pdb(99,'sec structure',0d0)
506 call elecont(lprint,ncont,icont)
507 ! print *,"after elecont"
508 if (nres_molec(1).eq.0) return
510 ! finding parallel beta
511 !d write (iout,*) '------- looking for parallel beta -----------'
517 if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
520 !d write (iout,*) i1,j1
526 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
527 freeres(i1,j1,nsec,isec)) goto 5
531 !d write (iout,*) i1,j1,not_done
535 if (i1-ii1.gt.1) then
539 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
543 bfrag(1,nbfrag)=ii1+1
545 bfrag(3,nbfrag)=jj1+1
546 bfrag(4,nbfrag)=min0(j1+1,nres)
550 isec(ij,nsec(ij))=nbeta
554 isec(ij,nsec(ij))=nbeta
560 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
561 "DefPropRes 'strand",nstrand,&
562 "' 'num = ",ii1-1,"..",i1-1,"'"
564 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
565 "DefPropRes 'strand",nstrand,&
566 "' 'num = ",ii1-1,"..",i1-1,"'"
570 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
571 "DefPropRes 'strand",nstrand,&
572 "' 'num = ",jj1-1,"..",j1-1,"'"
574 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
575 "DefPropRes 'strand",nstrand,&
576 "' 'num = ",jj1-1,"..",j1-1,"'"
578 write(12,'(a8,4i4)') &
579 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
585 ! finding alpha or 310 helix
592 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
595 if (j1.eq.i1+3 .and. &
596 ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
597 ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
598 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
599 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
602 if (nsec(ii1).eq.0) then
611 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
617 if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
622 if (j1-ii1.gt.5) then
634 write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
635 if (nhelix.le.9) then
636 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
637 "DefPropRes 'helix",nhelix,&
638 "' 'num = ",ii1-1,"..",j1-2,"'"
640 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
641 "DefPropRes 'helix",nhelix,&
642 "' 'num = ",ii1-1,"..",j1-2,"'"
648 if (nhelix.gt.0.and.lprint) then
649 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
651 if (nhelix.le.9) then
652 write(12,'(a8,i1,$)') " | helix",i
654 write(12,'(a8,i2,$)') " | helix",i
661 ! finding antiparallel beta
662 !d write (iout,*) '--------- looking for antiparallel beta ---------'
667 if (freeres(i1,j1,nsec,isec)) then
670 !d write (iout,*) i1,j1
677 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
678 freeres(i1,j1,nsec,isec)) goto 6
682 !d write (iout,*) i1,j1,not_done
686 if (i1-ii1.gt.1) then
690 bfrag(2,nbfrag)=min0(i1+1,nres)
691 bfrag(3,nbfrag)=min0(jj1+1,nres)
698 if (nsec(ij).le.2) then
699 isec(ij,nsec(ij))=nbeta
705 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
706 isec(ij,nsec(ij))=nbeta
712 write (iout,'(a,i3,4i4)')'antiparallel beta',&
713 nbeta,ii1-1,i1,jj1,j1-1
715 if (nstrand.le.9) then
716 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
717 "DefPropRes 'strand",nstrand,&
718 "' 'num = ",ii1-2,"..",i1-1,"'"
720 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
721 "DefPropRes 'strand",nstrand,&
722 "' 'num = ",ii1-2,"..",i1-1,"'"
725 if (nstrand.le.9) then
726 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
727 "DefPropRes 'strand",nstrand,&
728 "' 'num = ",j1-2,"..",jj1-1,"'"
730 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
731 "DefPropRes 'strand",nstrand,&
732 "' 'num = ",j1-2,"..",jj1-1,"'"
734 write(12,'(a8,4i4)') &
735 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
741 if (nstrand.gt.0.and.lprint) then
742 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
745 write(12,'(a9,i1,$)') " | strand",i
747 write(12,'(a9,i2,$)') " | strand",i
756 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
757 write(12,'(a20)') "XMacStand ribbon.mac"
759 if (nres_molec(1).eq.0) return
760 write(iout,*) 'UNRES seq:'
762 write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
766 write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
771 end subroutine secondary2
773 !-----------------------------------------------------------------------------
774 logical function freeres(i,j,nsec,isec)
776 ! implicit real*8 (a-h,o-z)
777 ! include 'DIMENSIONS'
778 integer,dimension(nres,4) :: isec !(maxres,4)
779 integer,dimension(nres) :: nsec !(maxres)
786 if (nsec(i).lt.0.or.nsec(j).lt.0) return
788 if (nsec(i).gt.1.or.nsec(j).gt.1) return
791 if (isec(i,k).eq.isec(j,l)) return
797 !-----------------------------------------------------------------------------
799 !-----------------------------------------------------------------------------
800 logical function seq_comp(itypea,itypeb,length)
803 integer :: length,itypea(length),itypeb(length)
806 if (itypea(i).ne.itypeb(i)) then
813 end function seq_comp
815 !-----------------------------------------------------------------------------
817 !-----------------------------------------------------------------------------
818 subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
820 ! implicit real*8 (a-h,o-z)
821 ! include 'DIMENSIONS'
822 ! include 'COMMON.CHAIN'
823 ! include 'COMMON.CONTACTS'
824 ! include 'COMMON.IOUNITS'
825 real(kind=8) :: przes(3),obr(3,3)
826 logical :: non_conv,lprn
827 real(kind=8) :: rms,frac,frac_nn,co
828 ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
832 ! print *,"before contact"
833 !elte(iout,*) "rms_nacc before contact"
834 call contact(.false.,ncont,icont,co)
835 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
836 frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
837 if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
838 'RMS deviation from the reference structure:',rms,&
839 ' % of native contacts:',frac*100,&
840 ' % of nonnative contacts:',frac_nn*100,&
844 end subroutine rms_nac_nnc
845 !-----------------------------------------------------------------------------
846 subroutine rmsd(drms)
848 use regularize_, only:fitsq
849 ! implicit real*8 (a-h,o-z)
850 ! include 'DIMENSIONS'
854 ! include 'COMMON.CHAIN'
855 ! include 'COMMON.IOUNITS'
856 ! include 'COMMON.INTERACT'
857 ! include 'COMMON.CONTROL'
859 real(kind=8) :: przes(3),obrot(3,3)
860 real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
863 real(kind=8) :: drms,rminroz,roznica
864 integer :: i,j,iatom,kkk,iti,k,mnum
866 !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
874 ! print *,"nz_start",nz_start," nz_end",nz_end
875 ! if (symetr.le.1) then
877 ! do i=nz_start,nz_end
881 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
882 ! crefcopy(k,iatom,kkk)=cref(k,i,kkk)
884 ! if (iz_sc.eq.1.and.iti.ne.10) then
887 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
888 ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
895 ! print *,nz_start,nz_end,nstart_seq-nstart_sup
898 iti=itype(i+nstart_seq-nstart_sup,molnum(i))
899 if (iti.eq.ntyp1_molec(molnum(i))) cycle
901 mnum=molnum(i+nstart_seq-nstart_sup)
903 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
904 crefcopy(k,iatom)=cref(k,i,kkk)
906 if (iz_sc.eq.1.and.iti.ne.10.and.mnum.lt.4) then
909 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
910 crefcopy(k,iatom)=cref(k,nres+i,kkk)
919 ! write (iout,*) 'Ccopy and CREFcopy'
920 ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
921 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
922 ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
923 ! & (crefcopy(j,k),j=1,3),k=1,iatom)
925 ! ----- end diagnostics
927 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
928 przes,obrot,non_conv)
930 print *,'Problems in FITSQ!!! rmsd'
931 write (iout,*) 'Problems in FITSQ!!! rmsd'
932 print *,'Ccopy and CREFcopy'
933 write (iout,*) 'Ccopy and CREFcopy'
934 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
935 (crefcopy(j,k),j=1,3),k=1,iatom)
936 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
937 (crefcopy(j,k),j=1,3),k=1,iatom)
939 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
945 ! write (iout,*) "roznica", roznica,kkk
946 if (roznica.le.rminroz) rminroz=roznica
948 drms=dsqrt(dabs(rminroz))
950 ! write (iout,*) "nperm,symetr", nperm,symetr
951 ! ---- end diagnostics
954 !-----------------------------------------------------------------------------
955 subroutine rmsd_csa(drms)
957 use regularize_, only:fitsq
958 ! implicit real*8 (a-h,o-z)
959 ! include 'DIMENSIONS'
963 ! include 'COMMON.CHAIN'
964 ! include 'COMMON.IOUNITS'
965 ! include 'COMMON.INTERACT'
967 real(kind=8) :: przes(3),obrot(3,3)
968 real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
969 integer :: kkk,iatom,ierror,ierrcode
973 real(kind=8) :: drms,roznica
975 allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
983 ccopy(k,iatom)=c(k,i)
984 crefcopy(k,iatom)=crefjlee(k,i)
986 if (iz_sc.eq.1.and.iti.ne.10) then
989 ccopy(k,iatom)=c(k,nres+i)
990 crefcopy(k,iatom)=crefjlee(k,nres+i)
995 call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
996 przes,obrot,non_conv)
998 print *,'Problems in FITSQ!!! rmsd_csa'
999 write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
1000 print *,'Ccopy and CREFcopy'
1001 write (iout,*) 'Ccopy and CREFcopy'
1002 print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
1003 (crefcopy(j,k),j=1,3),k=1,iatom)
1004 write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
1005 (crefcopy(j,k),j=1,3),k=1,iatom)
1007 call mpi_abort(mpi_comm_world,ierror,ierrcode)
1012 drms=dsqrt(dabs(roznica))
1014 end subroutine rmsd_csa
1015 !-----------------------------------------------------------------------------
1017 !-----------------------------------------------------------------------------
1021 use geometry, only:pinorm
1022 use random, only:ran_number,iran_num
1023 ! implicit real*8 (a-h,o-z)
1024 ! include 'DIMENSIONS'
1026 ! include 'COMMON.GEO'
1027 ! include 'COMMON.VAR'
1028 ! include 'COMMON.INTERACT'
1029 ! include 'COMMON.IOUNITS'
1030 ! include 'COMMON.DISTFIT'
1031 ! include 'COMMON.SBRIDGE'
1032 ! include 'COMMON.CONTROL'
1033 ! include 'COMMON.FFIELD'
1034 ! include 'COMMON.MINIM'
1035 ! include 'COMMON.CHAIN'
1036 real(kind=8) :: time0,time1
1037 real(kind=8) :: energy(0:n_ene),ee
1038 real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres)
1039 integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1040 logical :: debug,accepted
1041 real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1044 !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1046 call geom_to_var(nvar,var1)
1051 write(iout,*) 'etot=',0,etot,rms
1052 call secondary2(.false.)
1054 call write_pdb(0,'first structure',etot)
1063 betbol=1.0D0/(1.9858D-3*temp)
1065 d=ran_number(-pi,pi)
1066 ! phi(jr)=pinorm(phi(jr)+d)
1071 write(iout,*) 'etot=',1,etot0,rms
1072 call write_pdb(1,'perturb structure',etot0)
1076 d=ran_number(-da,da)
1078 phi(jr)=pinorm(phi(jr)+d)
1083 if (etot.lt.etot0) then
1087 xxr=ran_number(0.0D0,1.0D0)
1088 xxh=betbol*(etot-etot0)
1089 if (xxh.lt.50.0D0) then
1091 if (xxh.gt.xxr) accepted=.true.
1095 ! print *,etot0,etot,accepted
1099 write(iout,*) 'etot=',i,etot,rms
1100 call write_pdb(i,'MC structure',etot)
1102 ! call geom_to_var(nvar,var1)
1103 call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1104 call geom_to_var(nvar,var)
1105 call minimize(etot,var,iretcode,nfun)
1106 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1107 call var_to_geom(nvar,var)
1110 write(iout,*) 'etot mcm=',i,etot,rms
1111 call write_pdb(i+1,'MCM structure',etot)
1112 call var_to_geom(nvar,var1)
1120 ! call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1121 ! call geom_to_var(nvar,var)
1124 ! call write_pdb(998 ,'sc min',etot)
1126 ! call minimize(etot,var,iretcode,nfun)
1127 ! write(iout,*)'------------------------------------------------'
1128 ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1130 ! call var_to_geom(nvar,var)
1132 ! call write_pdb(999,'full min',etot)
1136 !-----------------------------------------------------------------------------
1141 ! implicit real*8 (a-h,o-z)
1142 ! include 'DIMENSIONS'
1144 ! include 'COMMON.GEO'
1145 ! include 'COMMON.VAR'
1146 ! include 'COMMON.INTERACT'
1147 ! include 'COMMON.IOUNITS'
1148 ! include 'COMMON.DISTFIT'
1149 ! include 'COMMON.SBRIDGE'
1150 ! include 'COMMON.CONTROL'
1151 ! include 'COMMON.FFIELD'
1152 ! include 'COMMON.MINIM'
1153 ! include 'COMMON.CHAIN'
1154 real(kind=8) :: time0,time1
1155 real(kind=8) :: energy(0:n_ene),ee
1156 real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1160 integer :: i,ij,ieval,iretcode,nfun
1161 real(kind=8) :: etot
1163 allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres)
1165 call geom_to_var(nvar,var1)
1169 write(iout,*) nnt,nct,etot
1170 call write_pdb(1,'first structure',etot)
1171 call secondary2(.true.)
1180 call var_to_geom(nvar,var1)
1181 write(iout,*) 'N16 test',(jdata(i),i=1,5)
1182 call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1184 call geom_to_var(nvar,var)
1190 call minimize(etot,var,iretcode,nfun)
1191 write(iout,*)'------------------------------------------------'
1192 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1198 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
1199 nfun/(time1-time0),' eval/s'
1201 call var_to_geom(nvar,var)
1203 call write_pdb(ij*100+99,'full min',etot)
1210 end subroutine test_n16
1212 !-----------------------------------------------------------------------------
1213 subroutine test_local
1215 ! implicit real*8 (a-h,o-z)
1216 ! include 'DIMENSIONS'
1217 ! include 'COMMON.GEO'
1218 ! include 'COMMON.VAR'
1219 ! include 'COMMON.INTERACT'
1220 ! include 'COMMON.IOUNITS'
1221 real(kind=8) :: time0,time1
1222 real(kind=8) :: energy(0:n_ene),ee
1223 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1225 real(kind=8) :: etot
1227 ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres)
1229 ! call geom_to_var(nvar,varia)
1230 call write_pdb(1,'first structure',0d0)
1234 write(iout,*) nnt,nct,etot
1236 write(iout,*) 'calling sc_move'
1237 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1238 write(iout,*) nft_sc,etot
1239 call write_pdb(2,'second structure',etot)
1241 write(iout,*) 'calling local_move'
1242 call local_move_init(.false.)
1243 call local_move(24,29,20d0,50d0)
1245 call write_pdb(3,'third structure',etot)
1247 write(iout,*) 'calling sc_move'
1248 call sc_move(24,29,5,10d0,nft_sc,etot)
1249 write(iout,*) nft_sc,etot
1250 call write_pdb(2,'last structure',etot)
1253 end subroutine test_local
1254 !-----------------------------------------------------------------------------
1257 ! implicit real*8 (a-h,o-z)
1258 ! include 'DIMENSIONS'
1259 ! include 'COMMON.GEO'
1260 ! include 'COMMON.VAR'
1261 ! include 'COMMON.INTERACT'
1262 ! include 'COMMON.IOUNITS'
1263 real(kind=8) :: time0,time1,etot
1264 real(kind=8) :: energy(0:n_ene),ee
1265 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1269 ! call geom_to_var(nvar,varia)
1270 call write_pdb(1,'first structure',0d0)
1274 write(iout,*) nnt,nct,etot
1276 write(iout,*) 'calling sc_move'
1278 call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1279 write(iout,*) nft_sc,etot
1280 call write_pdb(2,'second structure',etot)
1282 write(iout,*) 'calling sc_move 2nd time'
1284 call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1285 write(iout,*) nft_sc,etot
1286 call write_pdb(3,'last structure',etot)
1288 end subroutine test_sc
1289 !-----------------------------------------------------------------------------
1290 subroutine bgrow(bstrand,nbstrand,in,ind,new)
1292 ! implicit real*8 (a-h,o-z)
1293 ! include 'DIMENSIONS'
1294 ! include 'COMMON.CHAIN'
1295 integer,dimension(nres/3,6) :: bstrand !(maxres/3,6)
1298 integer :: nbstrand,in,ind,new,ishift,i
1300 ishift=iabs(bstrand(in,ind+4)-new)
1302 print *,'bgrow',bstrand(in,ind+4),new,ishift
1307 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1309 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1310 if (bstrand(i,5).lt.bstrand(i,6)) then
1311 bstrand(i,5)=bstrand(i,5)-ishift
1313 bstrand(i,5)=bstrand(i,5)+ishift
1318 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1320 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1321 if (bstrand(i,6).lt.bstrand(i,5)) then
1322 bstrand(i,6)=bstrand(i,6)-ishift
1324 bstrand(i,6)=bstrand(i,6)+ishift
1331 end subroutine bgrow
1332 !-----------------------------------------------------------------------------
1335 use geometry, only:dist
1337 ! implicit real*8 (a-h,o-z)
1338 ! include 'DIMENSIONS'
1340 ! include 'COMMON.GEO'
1341 ! include 'COMMON.CHAIN'
1342 ! include 'COMMON.IOUNITS'
1343 ! include 'COMMON.VAR'
1344 ! include 'COMMON.CONTROL'
1345 ! include 'COMMON.SBRIDGE'
1346 ! include 'COMMON.FFIELD'
1347 ! include 'COMMON.MINIM'
1349 ! include 'COMMON.DISTFIT'
1350 integer :: if(20,nres),nif,ifa(20)
1351 integer :: ibc(0:nres,0:nres),istrand(20)
1352 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1353 integer :: itmp(20,nres)
1354 real(kind=8) :: time0,time1
1355 real(kind=8) :: energy(0:n_ene),ee
1356 real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres)
1358 logical :: debug,ltest,usedbfrag(nres/3)
1359 character(len=50) :: linia
1361 integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1362 integer :: bstrand(nres/3,6),nbstrand
1363 real(kind=8) :: etot
1364 integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1365 in,ind,ifun,nfun,iretcode
1366 !------------------------
1369 !------------------------
1376 call geom_to_var(nvar,vorg)
1377 call secondary2(debug)
1379 if (nbfrag.le.1) return
1382 usedbfrag(i)=.false.
1386 nbetasheet=nbetasheet+1
1388 bstrand(1,1)=bfrag(1,1)
1389 bstrand(1,2)=bfrag(2,1)
1390 bstrand(1,3)=nbetasheet
1392 bstrand(1,5)=bfrag(1,1)
1393 bstrand(1,6)=bfrag(2,1)
1394 do i=bfrag(1,1),bfrag(2,1)
1395 betasheet(i)=nbetasheet
1399 bstrand(2,1)=bfrag(3,1)
1400 bstrand(2,2)=bfrag(4,1)
1401 bstrand(2,3)=nbetasheet
1402 bstrand(2,5)=bfrag(3,1)
1403 bstrand(2,6)=bfrag(4,1)
1405 if (bfrag(3,1).le.bfrag(4,1)) then
1407 do i=bfrag(3,1),bfrag(4,1)
1408 betasheet(i)=nbetasheet
1413 do i=bfrag(4,1),bfrag(3,1)
1414 betasheet(i)=nbetasheet
1421 do while (iused_nbfrag.ne.nbfrag)
1425 IF (.not.usedbfrag(j)) THEN
1427 write (*,*) j,(bfrag(i,j),i=1,4)
1429 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1431 write (*,*) '------------------'
1434 if (bfrag(3,j).le.bfrag(4,j)) then
1435 do i=bfrag(3,j),bfrag(4,j)
1436 if(betasheet(i).eq.nbetasheet) then
1438 do k=bfrag(3,j),bfrag(4,j)
1439 betasheet(k)=nbetasheet
1444 iused_nbfrag=iused_nbfrag+1
1445 do k=bfrag(1,j),bfrag(2,j)
1446 betasheet(k)=nbetasheet
1447 ibetasheet(k)=nbstrand
1449 if (bstrand(in,4).lt.0) then
1450 bstrand(nbstrand,1)=bfrag(2,j)
1451 bstrand(nbstrand,2)=bfrag(1,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).lt.bfrag(4,j)) then
1457 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1459 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1460 (bstrand(in,5)-bfrag(4,j))
1462 if(bstrand(in,2).gt.bfrag(3,j)) then
1463 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1465 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1466 (-bstrand(in,6)+bfrag(3,j))
1469 bstrand(nbstrand,1)=bfrag(1,j)
1470 bstrand(nbstrand,2)=bfrag(2,j)
1471 bstrand(nbstrand,3)=nbetasheet
1472 bstrand(nbstrand,4)=nbstrand
1473 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1474 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1475 if(bstrand(in,1).gt.bfrag(3,j)) then
1476 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1478 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1479 (-bstrand(in,5)+bfrag(3,j))
1481 if(bstrand(in,2).lt.bfrag(4,j)) then
1482 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1484 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1485 (bstrand(in,6)-bfrag(4,j))
1490 if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1491 in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1492 do k=bfrag(1,j),bfrag(2,j)
1493 betasheet(k)=nbetasheet
1498 iused_nbfrag=iused_nbfrag+1
1499 do k=bfrag(3,1),bfrag(4,1)
1500 betasheet(k)=nbetasheet
1501 ibetasheet(k)=nbstrand
1503 if (bstrand(in,4).lt.0) then
1504 bstrand(nbstrand,1)=bfrag(4,j)
1505 bstrand(nbstrand,2)=bfrag(3,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).lt.bfrag(2,j)) then
1511 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1513 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1514 (bstrand(in,5)-bfrag(2,j))
1516 if(bstrand(in,2).gt.bfrag(1,j)) then
1517 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1519 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1520 (-bstrand(in,6)+bfrag(1,j))
1523 bstrand(nbstrand,1)=bfrag(3,j)
1524 bstrand(nbstrand,2)=bfrag(4,j)
1525 bstrand(nbstrand,3)=nbetasheet
1526 bstrand(nbstrand,4)=nbstrand
1527 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1528 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1529 if(bstrand(in,1).gt.bfrag(1,j)) then
1530 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1532 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1533 (-bstrand(in,5)+bfrag(1,j))
1535 if(bstrand(in,2).lt.bfrag(2,j)) then
1536 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1538 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1539 (bstrand(in,6)-bfrag(2,j))
1546 do i=bfrag(4,j),bfrag(3,j)
1547 if(betasheet(i).eq.nbetasheet) then
1549 do k=bfrag(4,j),bfrag(3,j)
1550 betasheet(k)=nbetasheet
1555 iused_nbfrag=iused_nbfrag+1
1556 do k=bfrag(1,j),bfrag(2,j)
1557 betasheet(k)=nbetasheet
1558 ibetasheet(k)=nbstrand
1560 if (bstrand(in,4).lt.0) then
1561 bstrand(nbstrand,1)=bfrag(1,j)
1562 bstrand(nbstrand,2)=bfrag(2,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).lt.bfrag(3,j)) then
1568 call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1570 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1571 (bstrand(in,5)-bfrag(3,j))
1573 if(bstrand(in,2).gt.bfrag(4,j)) then
1574 call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1576 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1577 (-bstrand(in,6)+bfrag(4,j))
1580 bstrand(nbstrand,1)=bfrag(2,j)
1581 bstrand(nbstrand,2)=bfrag(1,j)
1582 bstrand(nbstrand,3)=nbetasheet
1583 bstrand(nbstrand,4)=-nbstrand
1584 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1585 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1586 if(bstrand(in,1).gt.bfrag(4,j)) then
1587 call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1589 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1590 (-bstrand(in,5)+bfrag(4,j))
1592 if(bstrand(in,2).lt.bfrag(3,j)) then
1593 call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1595 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1596 (bstrand(in,6)-bfrag(3,j))
1601 if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1602 in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1603 do k=bfrag(1,j),bfrag(2,j)
1604 betasheet(k)=nbetasheet
1609 iused_nbfrag=iused_nbfrag+1
1610 do k=bfrag(4,j),bfrag(3,j)
1611 betasheet(k)=nbetasheet
1612 ibetasheet(k)=nbstrand
1614 if (bstrand(in,4).lt.0) then
1615 bstrand(nbstrand,1)=bfrag(4,j)
1616 bstrand(nbstrand,2)=bfrag(3,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).lt.bfrag(2,j)) then
1622 call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1624 bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1625 (bstrand(in,5)-bfrag(2,j))
1627 if(bstrand(in,2).gt.bfrag(1,j)) then
1628 call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1630 bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1631 (-bstrand(in,6)+bfrag(1,j))
1634 bstrand(nbstrand,1)=bfrag(3,j)
1635 bstrand(nbstrand,2)=bfrag(4,j)
1636 bstrand(nbstrand,3)=nbetasheet
1637 bstrand(nbstrand,4)=-nbstrand
1638 bstrand(nbstrand,5)=bstrand(nbstrand,1)
1639 bstrand(nbstrand,6)=bstrand(nbstrand,2)
1640 if(bstrand(in,1).gt.bfrag(1,j)) then
1641 call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1643 bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1644 (-bstrand(in,5)+bfrag(1,j))
1646 if(bstrand(in,2).lt.bfrag(2,j)) then
1647 call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1649 bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1650 (bstrand(in,6)-bfrag(2,j))
1664 do while (usedbfrag(j))
1669 nbetasheet=nbetasheet+1
1670 bstrand(nbstrand,1)=bfrag(1,j)
1671 bstrand(nbstrand,2)=bfrag(2,j)
1672 bstrand(nbstrand,3)=nbetasheet
1673 bstrand(nbstrand,5)=bfrag(1,j)
1674 bstrand(nbstrand,6)=bfrag(2,j)
1676 bstrand(nbstrand,4)=nbstrand
1677 do i=bfrag(1,j),bfrag(2,j)
1678 betasheet(i)=nbetasheet
1679 ibetasheet(i)=nbstrand
1683 bstrand(nbstrand,1)=bfrag(3,j)
1684 bstrand(nbstrand,2)=bfrag(4,j)
1685 bstrand(nbstrand,3)=nbetasheet
1686 bstrand(nbstrand,5)=bfrag(3,j)
1687 bstrand(nbstrand,6)=bfrag(4,j)
1689 if (bfrag(3,j).le.bfrag(4,j)) then
1690 bstrand(nbstrand,4)=nbstrand
1691 do i=bfrag(3,j),bfrag(4,j)
1692 betasheet(i)=nbetasheet
1693 ibetasheet(i)=nbstrand
1696 bstrand(nbstrand,4)=-nbstrand
1697 do i=bfrag(4,j),bfrag(3,j)
1698 betasheet(i)=nbetasheet
1699 ibetasheet(i)=nbstrand
1703 iused_nbfrag=iused_nbfrag+1
1709 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1716 if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1720 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1723 !------------------------
1727 if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1728 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1730 ifb(nifb,1)=bstrand(i,4)
1731 ifb(nifb,2)=bstrand(j,4)
1738 write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1744 write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1746 nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1748 if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1749 nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1755 if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1756 if (if(j,i).gt.0) then
1757 if(betasheet(if(j,i)).eq.0 .or. &
1758 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1763 write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1766 ! read (inp,*) (ifa(i),i=1,4)
1768 ! read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1772 !------------------------
1777 !ccccccccccccccccccccccccccccccccc
1779 !ccccccccccccccccccccccccccccccccc
1783 istrand(is-j+1)=int(ii/is**(is-j))
1784 ii=ii-istrand(is-j+1)*is**(is-j)
1788 istrand(k)=istrand(k)+1
1789 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1793 if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1794 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1803 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1804 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1805 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1806 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1812 if (mod(isa,2).eq.0) then
1814 if (istrand(k).eq.1) ltest=.false.
1817 do k=(isa+1)/2+1,isa
1818 if (istrand(k).eq.1) ltest=.false.
1822 IF (ltest.and.lifb0.eq.1) THEN
1825 call var_to_geom(nvar,vorg)
1827 write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1828 write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1829 write (linia,'(10i3)') (istrand(k),k=1,isa)
1839 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1841 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1845 itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1851 write(*,*) (itmp(j,i),j=1,4)
1855 ! ifa(1),ifa(2),ifa(3),ifa(4)
1856 ! if(1,i),if(2,i),if(3,i),if(4,i)
1861 ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1862 ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1863 -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1864 -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1872 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1874 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1878 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1880 ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1883 !------------------------
1886 ! freeze sec.elements
1896 do i=bfrag(1,j),bfrag(2,j)
1901 if (bfrag(3,j).le.bfrag(4,j)) then
1902 do i=bfrag(3,j),bfrag(4,j)
1908 do i=bfrag(4,j),bfrag(3,j)
1916 do i=hfrag(1,j),hfrag(2,j)
1924 !------------------------
1925 ! generate constrains
1933 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1941 else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1949 else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1957 else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1965 else if ( ibc(i,j).gt.0 ) then
1966 d0(ind)=DIST(i,ibc(i,j))
1973 else if ( ibc(j,i).gt.0 ) then
1974 d0(ind)=DIST(ibc(j,i),j)
1988 !d--------------------------
1990 write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1991 ibc(jhpb(i),ihpb(i)),' --',&
1992 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1999 call contact_cp_min(varia,ifun,iconf,linia,debug)
2004 call minimize(etot,varia,iretcode,nfun)
2005 write(iout,*)'------------------------------------------------'
2006 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2012 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2013 nfun/(time1-time0),' eval/s'
2015 write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
2016 call var_to_geom(nvar,varia)
2018 call write_pdb(900+iconf,linia,etot)
2023 call enerprint(energy)
2025 !d call briefout(0,etot)
2026 !d call secondary2(.true.)
2030 !ccccccccccccccccccccccccccccccccccc
2033 !ccccccccccccccccccccccccccccccccccc
2036 10 write (iout,'(a)') 'Error reading test structure.'
2038 end subroutine test11
2039 !-----------------------------------------------------------------------------
2042 use geometry, only:dist
2044 ! implicit real*8 (a-h,o-z)
2045 ! include 'DIMENSIONS'
2047 ! include 'COMMON.GEO'
2048 ! include 'COMMON.CHAIN'
2049 ! include 'COMMON.IOUNITS'
2050 ! include 'COMMON.VAR'
2051 ! include 'COMMON.CONTROL'
2052 ! include 'COMMON.SBRIDGE'
2053 ! include 'COMMON.FFIELD'
2054 ! include 'COMMON.MINIM'
2056 ! include 'COMMON.DISTFIT'
2057 integer :: if(3,nres),nif
2058 integer :: ibc(nres,nres),istrand(20)
2059 integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2060 real(kind=8) :: time0,time1
2061 real(kind=8) :: energy(0:n_ene),ee
2062 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2064 logical :: debug,ltest
2065 character(len=50) :: linia
2066 integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2067 real(kind=8) :: etot
2070 read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2073 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2077 !------------------------
2078 call secondary2(debug)
2079 !------------------------
2087 ! freeze sec.elements and store indexes for beta constrains
2097 do i=bfrag(1,j),bfrag(2,j)
2102 if (bfrag(3,j).le.bfrag(4,j)) then
2103 do i=bfrag(3,j),bfrag(4,j)
2107 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2110 do i=bfrag(4,j),bfrag(3,j)
2114 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2119 do i=hfrag(1,j),hfrag(2,j)
2128 ! ---------------- test --------------
2130 if (ibc(if(1,i),if(2,i)).eq.-1) then
2131 ibc(if(1,i),if(2,i))=if(3,i)
2132 ibc(if(1,i),if(3,i))=if(2,i)
2133 else if (ibc(if(2,i),if(1,i)).eq.-1) then
2134 ibc(if(2,i),if(1,i))=0
2135 ibc(if(1,i),if(2,i))=if(3,i)
2136 ibc(if(1,i),if(3,i))=if(2,i)
2138 ibc(if(1,i),if(2,i))=if(3,i)
2139 ibc(if(1,i),if(3,i))=if(2,i)
2145 if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j)
2148 !------------------------
2154 if ( ibc(i,j).eq.-1 ) then
2162 else if ( ibc(i,j).gt.0 ) then
2163 d0(ind)=DIST(i,ibc(i,j))
2170 else if ( ibc(j,i).gt.0 ) then
2171 d0(ind)=DIST(ibc(j,i),j)
2185 !d--------------------------
2186 write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2187 ibc(jhpb(i),ihpb(i)),' --',&
2188 ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2196 call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2201 call minimize(etot,varia,iretcode,nfun)
2202 write(iout,*)'------------------------------------------------'
2203 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2209 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2210 nfun/(time1-time0),' eval/s'
2213 call var_to_geom(nvar,varia)
2215 call write_pdb(999,'full min',etot)
2220 call enerprint(energy)
2222 call briefout(0,etot)
2223 call secondary2(.true.)
2226 10 write (iout,'(a)') 'Error reading test structure.'
2228 end subroutine test3
2229 !-----------------------------------------------------------------------------
2233 ! implicit real*8 (a-h,o-z)
2234 ! include 'DIMENSIONS'
2236 ! include 'COMMON.GEO'
2237 ! include 'COMMON.CHAIN'
2238 ! include 'COMMON.IOUNITS'
2239 ! include 'COMMON.VAR'
2240 ! include 'COMMON.CONTROL'
2241 ! include 'COMMON.SBRIDGE'
2242 ! include 'COMMON.FFIELD'
2243 ! include 'COMMON.MINIM'
2245 ! include 'COMMON.DISTFIT'
2246 integer :: if(2,2),ind
2247 integer :: iff(nres)
2248 real(kind=8) :: time0,time1
2249 real(kind=8) :: energy(0:n_ene),ee
2250 real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2251 theta1,phi1,alph1,omeg1 !(maxres)
2252 real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres)
2254 integer :: i,j,nn,ifun,iretcode,nfun
2255 real(kind=8) :: etot
2258 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2259 write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2260 read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2261 read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2262 read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2263 read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2265 theta2(i)=deg2rad*theta2(i)
2266 phi2(i)=deg2rad*phi2(i)
2267 alph2(i)=deg2rad*alph2(i)
2268 omeg2(i)=deg2rad*omeg2(i)
2282 !------------------------
2287 do i=if(j,1),if(j,2)
2293 call geom_to_var(nvar,varia)
2294 call write_pdb(1,'first structure',0d0)
2296 call secondary(.true.)
2298 call secondary2(.true.)
2301 if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2302 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2303 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2306 if (bfrag(3,j).lt.bfrag(4,j)) then
2307 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2308 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2309 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2311 write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2312 "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2313 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2326 call geom_to_var(nvar,varia2)
2327 call write_pdb(2,'second structure',0d0)
2331 !-------------------------------------------------------
2334 call contact_cp(varia,varia2,iff,ifun,7)
2339 call minimize(etot,varia,iretcode,nfun)
2340 write(iout,*)'------------------------------------------------'
2341 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2347 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
2348 nfun/(time1-time0),' eval/s'
2351 call var_to_geom(nvar,varia)
2353 call write_pdb(999,'full min',etot)
2358 call enerprint(energy)
2360 call briefout(0,etot)
2363 10 write (iout,'(a)') 'Error reading test structure.'
2365 end subroutine test__
2366 !-----------------------------------------------------------------------------
2367 subroutine secondary(lprint)
2369 ! implicit real*8 (a-h,o-z)
2370 ! include 'DIMENSIONS'
2371 ! include 'COMMON.CHAIN'
2372 ! include 'COMMON.IOUNITS'
2373 ! include 'COMMON.DISTFIT'
2375 integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2376 logical :: lprint,not_done
2377 real(kind=4) :: dcont(nres*nres/2),d
2378 real(kind=4) :: rcomp = 7.0
2379 real(kind=4) :: rbeta = 5.2
2380 real(kind=4) :: ralfa = 5.2
2381 real(kind=4) :: r310 = 6.6
2382 real(kind=8),dimension(3) :: xpi,xpj
2383 integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2385 call chainbuild_cart
2386 !d call write_pdb(99,'sec structure',0d0)
2398 xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2402 xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2404 !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2405 !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2406 !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2407 !d print *,'CA',i,j,d
2408 d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2409 (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2410 (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2411 if ( d.lt.rcomp*rcomp) then
2415 dcont(ncont)=sqrt(d)
2421 write (iout,'(a)') '#PP contact map distances:'
2423 write (iout,'(3i4,f10.5)') &
2424 i,icont(1,i),icont(2,i),dcont(i)
2428 ! finding parallel beta
2429 !d write (iout,*) '------- looking for parallel beta -----------'
2435 if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2436 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2437 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2438 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2439 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2440 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2444 !d write (iout,*) i1,j1,dcont(i)
2450 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2451 .and. dcont(j).le.rbeta .and. &
2452 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2453 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2454 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2455 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2456 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2461 !d write (iout,*) i1,j1,dcont(j),not_done
2465 if (i1-ii1.gt.1) then
2469 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2478 isec(ij,1)=isec(ij,1)+1
2479 isec(ij,1+isec(ij,1))=nbeta
2482 isec(ij,1)=isec(ij,1)+1
2483 isec(ij,1+isec(ij,1))=nbeta
2488 if (nbeta.le.9) then
2489 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2490 "DefPropRes 'strand",nstrand,&
2491 "' 'num = ",ii1-1,"..",i1-1,"'"
2493 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2494 "DefPropRes 'strand",nstrand,&
2495 "' 'num = ",ii1-1,"..",i1-1,"'"
2498 if (nbeta.le.9) then
2499 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2500 "DefPropRes 'strand",nstrand,&
2501 "' 'num = ",jj1-1,"..",j1-1,"'"
2503 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2504 "DefPropRes 'strand",nstrand,&
2505 "' 'num = ",jj1-1,"..",j1-1,"'"
2507 write(12,'(a8,4i4)') &
2508 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2514 ! finding antiparallel beta
2515 !d write (iout,*) '--------- looking for antiparallel beta ---------'
2520 if (dcont(i).le.rbeta.and. &
2521 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2522 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2523 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2524 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2525 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2529 !d write (iout,*) i1,j1,dcont(i)
2536 if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
2537 isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2538 (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2539 (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2540 (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2541 (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2542 .and. dcont(j).le.rbeta ) goto 6
2546 !d write (iout,*) i1,j1,dcont(j),not_done
2550 if (i1-ii1.gt.1) then
2551 if(lprint)write (iout,*)'antiparallel beta',&
2552 nbeta,ii1-1,i1,jj1,j1-1
2555 bfrag(1,nbfrag)=max0(ii1-1,1)
2558 bfrag(4,nbfrag)=max0(j1-1,1)
2563 isec(ij,1)=isec(ij,1)+1
2564 isec(ij,1+isec(ij,1))=nbeta
2568 isec(ij,1)=isec(ij,1)+1
2569 isec(ij,1+isec(ij,1))=nbeta
2575 if (nstrand.le.9) then
2576 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2577 "DefPropRes 'strand",nstrand,&
2578 "' 'num = ",ii1-2,"..",i1-1,"'"
2580 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2581 "DefPropRes 'strand",nstrand,&
2582 "' 'num = ",ii1-2,"..",i1-1,"'"
2585 if (nstrand.le.9) then
2586 write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2587 "DefPropRes 'strand",nstrand,&
2588 "' 'num = ",j1-2,"..",jj1-1,"'"
2590 write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2591 "DefPropRes 'strand",nstrand,&
2592 "' 'num = ",j1-2,"..",jj1-1,"'"
2594 write(12,'(a8,4i4)') &
2595 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2601 if (nstrand.gt.0.and.lprint) then
2602 write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2605 write(12,'(a9,i1,$)') " | strand",i
2607 write(12,'(a9,i2,$)') " | strand",i
2610 write(12,'(a1)') "'"
2614 ! finding alpha or 310 helix
2620 if (j1.eq.i1+3.and.dcont(i).le.r310 &
2621 .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2622 !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2623 !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2626 if (isec(ii1,1).eq.0) then
2635 if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2639 !d write (iout,*) i1,j1,not_done
2642 if (j1-ii1.gt.4) then
2644 !d write (iout,*)'helix',nhelix,ii1,j1
2648 hfrag(2,nhfrag)=max0(j1-1,1)
2654 write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2655 if (nhelix.le.9) then
2656 write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2657 "DefPropRes 'helix",nhelix,&
2658 "' 'num = ",ii1-1,"..",j1-2,"'"
2660 write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2661 "DefPropRes 'helix",nhelix,&
2662 "' 'num = ",ii1-1,"..",j1-2,"'"
2669 if (nhelix.gt.0.and.lprint) then
2670 write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2672 if (nhelix.le.9) then
2673 write(12,'(a8,i1,$)') " | helix",i
2675 write(12,'(a8,i2,$)') " | helix",i
2678 write(12,'(a1)') "'"
2682 write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2683 write(12,'(a20)') "XMacStand ribbon.mac"
2687 end subroutine secondary
2688 !-----------------------------------------------------------------------------
2689 subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2691 use geometry, only:dist
2693 ! implicit real*8 (a-h,o-z)
2694 ! include 'DIMENSIONS'
2696 ! include 'COMMON.SBRIDGE'
2697 ! include 'COMMON.FFIELD'
2698 ! include 'COMMON.IOUNITS'
2699 ! include 'COMMON.DISTFIT'
2700 ! include 'COMMON.VAR'
2701 ! include 'COMMON.CHAIN'
2702 ! include 'COMMON.MINIM'
2704 character(len=50) :: linia
2706 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2707 real(kind=8) :: time0,time1
2708 integer :: iff(nres),ieval
2709 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2712 integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2713 real(kind=8) :: wstrain0,etot
2715 maxres22=nres*(nres+1)/2
2717 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2718 call var_to_geom(nvar,var)
2725 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2747 call var_to_geom(nvar,var2)
2750 if ( iff(i).eq.1 ) then
2759 !d call write_pdb(3,'combined structure',0d0)
2760 !d time0=MPI_WTIME()
2763 NY=((NRES-4)*(NRES-5))/2
2764 call distfit(.true.,200)
2766 !d time1=MPI_WTIME()
2767 !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2777 call geom_to_var(nvar,var)
2778 !d time0=MPI_WTIME()
2779 call minimize(etot,var,iretcode,nfun)
2780 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2782 !d time1=MPI_WTIME()
2783 !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2784 !d & nfun/(time1-time0),' SOFT eval/s'
2785 call var_to_geom(nvar,var)
2791 if (iff(1).eq.1) then
2797 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2802 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2808 if (iff(nres).eq.1) then
2814 !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2815 !d & "select",ij(1),"-",ij(2),
2816 !d & ",",ij(3),"-",ij(4)
2817 !d call write_pdb(in_pdb,linia,etot)
2823 !d time0=MPI_WTIME()
2824 call minimize(etot,var,iretcode,nfun)
2825 !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2828 !d time1=MPI_WTIME()
2829 !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2830 !d & nfun/(time1-time0),' eval/s'
2831 !d call var_to_geom(nvar,var)
2833 !d call write_pdb(6,'dist structure',etot)
2842 end subroutine contact_cp2
2843 !-----------------------------------------------------------------------------
2844 subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2846 use geometry, only:dist
2848 ! implicit real*8 (a-h,o-z)
2849 ! include 'DIMENSIONS'
2850 ! include 'COMMON.SBRIDGE'
2851 ! include 'COMMON.FFIELD'
2852 ! include 'COMMON.IOUNITS'
2853 ! include 'COMMON.DISTFIT'
2854 ! include 'COMMON.VAR'
2855 ! include 'COMMON.CHAIN'
2856 ! include 'COMMON.MINIM'
2858 character(len=50) :: linia
2860 real(kind=8) :: energy(0:n_ene)
2861 real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres)
2862 real(kind=8) :: time0,time1
2863 integer :: iff(nres),ieval
2864 real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres)
2868 integer :: in_pdb,i,j,ind,iwsk
2872 if (ieval.eq.-1) debug=.true.
2876 ! store selected dist. constrains from 1st structure
2879 ! Intercept NaNs in the coordinates
2880 ! write(iout,*) (var(i),i=1,nvar)
2885 if (x_sum.ne.x_sum) then
2886 write(iout,*)" *** contact_cp : Found NaN in coordinates"
2888 print *," *** contact_cp : Found NaN in coordinates"
2894 call var_to_geom(nvar,var)
2901 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2924 ! freeze sec.elements from 2nd structure
2932 call var_to_geom(nvar,var2)
2933 call secondary2(debug)
2935 do i=bfrag(1,j),bfrag(2,j)
2940 if (bfrag(3,j).le.bfrag(4,j)) then
2941 do i=bfrag(3,j),bfrag(4,j)
2947 do i=bfrag(4,j),bfrag(3,j)
2955 do i=hfrag(1,j),hfrag(2,j)
2964 ! copy selected res from 1st to 2nd structure
2968 if ( iff(i).eq.1 ) then
2978 ! prepare description in linia variable
2982 if (iff(1).eq.1) then
2988 if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2993 if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2999 if (iff(nres).eq.1) then
3004 write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3005 "SELECT",ij(1)-1,"-",ij(2)-1,&
3006 ",",ij(3)-1,"-",ij(4)-1
3012 call contact_cp_min(var,ieval,in_pdb,linia,debug)
3015 end subroutine contact_cp
3016 !-----------------------------------------------------------------------------
3017 subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3021 ! input : theta,phi,alph,omeg,in_pdb,linia,debug
3022 ! output : var,ieval
3024 ! implicit real*8 (a-h,o-z)
3025 ! include 'DIMENSIONS'
3027 ! include 'COMMON.SBRIDGE'
3028 ! include 'COMMON.FFIELD'
3029 ! include 'COMMON.IOUNITS'
3030 ! include 'COMMON.DISTFIT'
3031 ! include 'COMMON.VAR'
3032 ! include 'COMMON.CHAIN'
3033 ! include 'COMMON.MINIM'
3035 character(len=50) :: linia
3037 real(kind=8) :: energy(0:n_ene)
3038 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3039 real(kind=8) :: time0,time1
3040 integer :: ieval,info(3)
3041 logical :: debug,fail,reduce,change !check_var,
3044 integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3046 real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3047 wtor_d01,wstrain0,etot
3049 write(iout,'(a20,i6,a20)') &
3050 '------------------',in_pdb,'-------------------'
3054 call write_pdb(1000+in_pdb,'combined structure',0d0)
3061 ! run optimization of distances
3063 ! uses d0(),w() and mask() for frozen 2D
3065 !test---------------------------------------------
3067 !test NY=((NRES-4)*(NRES-5))/2
3068 !test call distfit(debug,5000)
3100 call geom_to_var(nvar,var)
3101 !de change=reduce(var)
3102 if (check_var(var,info)) then
3103 write(iout,*) 'cp_min error in input'
3104 print *,'cp_min error in input'
3108 !d call etotal(energy(0))
3109 !d call enerprint(energy(0))
3113 !dtest call minimize(etot,var,iretcode,nfun)
3114 !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3117 !d call etotal(energy(0))
3118 !d call enerprint(energy(0))
3137 !test--------------------------------------------------
3143 write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3144 call write_pdb(2000+in_pdb,'distfit structure',0d0)
3152 ! run soft pot. optimization
3154 ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3156 ! mask_phi(),mask_theta(),mask_side(),mask_r
3162 !de change=reduce(var)
3163 !de if (check_var(var,info)) write(iout,*) 'error before soft'
3167 call minimize(etot,var,iretcode,nfun)
3169 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3173 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3174 nfun/(time1-time0),' SOFT eval/s'
3177 call var_to_geom(nvar,var)
3179 call write_pdb(3000+in_pdb,'soft structure',etot)
3183 ! run full UNRES optimization with constrains and frozen 2D
3184 ! the same variables as soft pot. optimizatio
3190 ! check overlaps before calling full UNRES minim
3192 call var_to_geom(nvar,var)
3196 write(iout,*) 'N7 ',energy(0)
3197 if (energy(0).ne.energy(0)) then
3198 write(iout,*) 'N7 error - gives NaN',energy(0)
3202 if (energy(1).eq.1.0d20) then
3203 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3204 call overlap_sc(fail)
3208 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3220 !dte time0=MPI_WTIME()
3221 !de change=reduce(var)
3222 !de if (check_var(var,info)) then
3223 !de write(iout,*) 'error before mask dist'
3224 !de call var_to_geom(nvar,var)
3226 !de call write_pdb(10000+in_pdb,'before mask dist',etot)
3228 !dte call minimize(etot,var,iretcode,nfun)
3229 !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3230 !dte & ' eval ',nfun
3231 !dte ieval=ieval+nfun
3233 !dte time1=MPI_WTIME()
3234 !dte write (iout,'(a,f6.2,f8.2,a)')
3235 !dte & ' Time for mask dist min.',time1-time0,
3236 !dte & nfun/(time1-time0),' eval/s'
3237 !dte call flush(iout)
3240 call var_to_geom(nvar,var)
3242 call write_pdb(4000+in_pdb,'mask dist',etot)
3245 ! switch off freezing of 2D and
3246 ! run full UNRES optimization with constrains
3252 !de change=reduce(var)
3253 !de if (check_var(var,info)) then
3254 !de write(iout,*) 'error before dist'
3255 !de call var_to_geom(nvar,var)
3257 !de call write_pdb(11000+in_pdb,'before dist',etot)
3260 call minimize(etot,var,iretcode,nfun)
3262 !de change=reduce(var)
3263 !de if (check_var(var,info)) then
3264 !de write(iout,*) 'error after dist',ico
3265 !de call var_to_geom(nvar,var)
3267 !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3269 write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3275 write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3276 nfun/(time1-time0),' eval/s'
3278 !de call etotal(energy(0))
3279 !de write(iout,*) 'N7 after dist',energy(0)
3283 call var_to_geom(nvar,var)
3285 call write_pdb(in_pdb,linia,etot)
3297 end subroutine contact_cp_min
3298 !-----------------------------------------------------------------------------
3301 use geometry, only:dist
3303 ! implicit real*8 (a-h,o-z)
3304 ! include 'DIMENSIONS'
3306 ! include 'COMMON.GEO'
3307 ! include 'COMMON.CHAIN'
3308 ! include 'COMMON.IOUNITS'
3309 ! include 'COMMON.VAR'
3310 ! include 'COMMON.CONTROL'
3311 ! include 'COMMON.SBRIDGE'
3312 ! include 'COMMON.FFIELD'
3313 ! include 'COMMON.MINIM'
3314 ! include 'COMMON.INTERACT'
3316 ! include 'COMMON.DISTFIT'
3317 integer :: iff(nres)
3318 real(kind=8) :: time0,time1
3319 real(kind=8) :: energy(0:n_ene),ee
3320 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3322 logical :: debug,ltest,fail
3323 character(len=50) :: linia
3324 integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3326 real(kind=8) :: wstrain0,wang0,etot
3332 !------------------------
3334 ! freeze sec.elements
3344 do i=bfrag(1,j),bfrag(2,j)
3349 if (bfrag(3,j).le.bfrag(4,j)) then
3350 do i=bfrag(3,j),bfrag(4,j)
3356 do i=bfrag(4,j),bfrag(3,j)
3364 do i=hfrag(1,j),hfrag(2,j)
3376 ! store dist. constrains
3380 if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3385 dhpb(nhpb)=DIST(i,j)
3393 call write_pdb(100+in_pdb,'input reg. structure',0d0)
3403 ! run soft pot. optimization
3409 call geom_to_var(nvar,var)
3413 call minimize(etot,var,iretcode,nfun)
3415 write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3419 write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3420 nfun/(time1-time0),' SOFT eval/s'
3422 call var_to_geom(nvar,var)
3424 call write_pdb(300+in_pdb,'soft structure',etot)
3427 ! run full UNRES optimization with constrains and frozen 2D
3428 ! the same variables as soft pot. optimizatio
3437 call minimize(etot,var,iretcode,nfun)
3438 write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3445 write (iout,'(a,f6.2,f8.2,a)') &
3446 ' Time for mask dist min.',time1-time0,&
3447 nfun/(time1-time0),' eval/s'
3449 call var_to_geom(nvar,var)
3451 call write_pdb(400+in_pdb,'mask & dist',etot)
3454 ! switch off constrains and
3455 ! run full UNRES optimization with frozen 2D
3470 call minimize(etot,var,iretcode,nfun)
3471 write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3477 write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,&
3478 nfun/(time1-time0),' eval/s'
3482 call var_to_geom(nvar,var)
3484 call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3491 ! run full UNRES optimization with constrains and NO frozen 2D
3501 wstrain=wstrain0/ico
3506 call minimize(etot,var,iretcode,nfun)
3507 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3508 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3515 write (iout,'(a,f6.2,f8.2,a)') &
3516 ' Time for dist min.',time1-time0,&
3517 nfun/(time1-time0),' eval/s'
3519 call var_to_geom(nvar,var)
3521 call write_pdb(600+in_pdb+ico,'dist cons',etot)
3538 call minimize(etot,var,iretcode,nfun)
3539 write(iout,*)'------------------------------------------------'
3540 write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3546 write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,&
3547 nfun/(time1-time0),' eval/s'
3550 call var_to_geom(nvar,var)
3552 call write_pdb(999,'full min',etot)
3556 end subroutine softreg
3557 !-----------------------------------------------------------------------------
3558 subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3560 use geometry, only:dist
3562 ! implicit real*8 (a-h,o-z)
3563 ! include 'DIMENSIONS'
3565 ! include 'COMMON.GEO'
3566 ! include 'COMMON.VAR'
3567 ! include 'COMMON.INTERACT'
3568 ! include 'COMMON.IOUNITS'
3569 ! include 'COMMON.DISTFIT'
3570 ! include 'COMMON.SBRIDGE'
3571 ! include 'COMMON.CONTROL'
3572 ! include 'COMMON.FFIELD'
3573 ! include 'COMMON.MINIM'
3574 ! include 'COMMON.CHAIN'
3575 real(kind=8) :: time0,time1
3576 real(kind=8) :: energy(0:n_ene),ee
3577 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3578 integer :: jdata(5),isec(nres)
3581 integer :: i1,i2,i3,i4,i5,ieval,ij
3582 integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3583 real(kind=8) :: etot,wscloc0,wstrain0
3591 call secondary2(.false.)
3597 do i=bfrag(1,j),bfrag(2,j)
3600 do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3605 do i=hfrag(1,j),hfrag(2,j)
3611 ! cut strands at the ends
3613 if (jdata(2)-jdata(1).gt.3) then
3616 if (jdata(3).lt.jdata(4)) then
3626 !v call etotal(energy(0))
3628 !v write(iout,*) nnt,nct,etot
3629 !v call write_pdb(ij*100,'first structure',etot)
3630 !v write(iout,*) 'N16 test',(jdata(i),i=1,5)
3632 !------------------------
3633 ! generate constrains
3636 if(ishift.eq.0) ishift=-2
3639 do i=jdata(1),jdata(2)
3641 if(jdata(4).gt.jdata(3))then
3642 do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3644 !d print *,i,j,j+ishift
3649 dhpb(nhpb)=DIST(i,j+ishift)
3652 do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3654 !d print *,i,j,j+ishift
3659 dhpb(nhpb)=DIST(i,j+ishift)
3666 if(isec(i).gt.0.or.isec(j).gt.0) then
3672 dhpb(nhpb)=DIST(i,j)
3679 call geom_to_var(nvar,var)
3686 wstrain=wstrain0/ico
3688 !v time0=MPI_WTIME()
3689 call minimize(etot,var,iretcode,nfun)
3690 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3691 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3694 !v time1=MPI_WTIME()
3695 !v write (iout,'(a,f6.2,f8.2,a)')
3696 !v & ' Time for dist min.',time1-time0,
3697 !v & nfun/(time1-time0),' eval/s'
3698 !v call var_to_geom(nvar,var)
3700 !v call write_pdb(ij*100+ico,'dist cons',etot)
3712 call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3715 !v call etotal(energy(0))
3717 !v call write_pdb(ij*100+10,'sc_move',etot)
3719 !d print *,nft_sc,etot
3722 end subroutine beta_slide
3723 !-----------------------------------------------------------------------------
3724 subroutine beta_zip(i1,i2,ieval,ij)
3727 ! implicit real*8 (a-h,o-z)
3728 ! include 'DIMENSIONS'
3730 ! include 'COMMON.GEO'
3731 ! include 'COMMON.VAR'
3732 ! include 'COMMON.INTERACT'
3733 ! include 'COMMON.IOUNITS'
3734 ! include 'COMMON.DISTFIT'
3735 ! include 'COMMON.SBRIDGE'
3736 ! include 'COMMON.CONTROL'
3737 ! include 'COMMON.FFIELD'
3738 ! include 'COMMON.MINIM'
3739 ! include 'COMMON.CHAIN'
3740 real(kind=8) :: time0,time1
3741 real(kind=8) :: energy(0:n_ene),ee
3742 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
3743 character(len=10) :: test
3745 integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3746 real(kind=8) :: etot,wstrain0
3748 !v call etotal(energy(0))
3750 !v write(test,'(2i5)') i1,i2
3751 !v call write_pdb(ij*100,test,etot)
3752 !v write(iout,*) 'N17 test',i1,i2,etot,ij
3755 ! generate constrains
3766 call geom_to_var(nvar,var)
3772 wstrain=wstrain0/ico
3773 !v time0=MPI_WTIME()
3774 call minimize(etot,var,iretcode,nfun)
3775 write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3776 ' SUMSL DIST',wstrain,' return code is',iretcode,&
3779 !v time1=MPI_WTIME()
3780 !v write (iout,'(a,f6.2,f8.2,a)')
3781 !v & ' Time for dist min.',time1-time0,
3782 !v & nfun/(time1-time0),' eval/s'
3783 ! do not comment the next line
3784 call var_to_geom(nvar,var)
3786 !v call write_pdb(ij*100+ico,'dist cons',etot)
3794 !v call etotal(energy(0))
3796 !v write(iout,*) 'N17 test end',i1,i2,etot,ij
3799 end subroutine beta_zip
3800 !-----------------------------------------------------------------------------
3802 !-----------------------------------------------------------------------------
3803 subroutine thread_seq
3805 use geometry, only:dist
3806 use random, only:iran_num
3807 use control, only:tcpu
3808 use regularize_, only:regularize
3809 use mcm_data, only: nsave_part,nacc_tot
3810 ! Thread the sequence through a database of known structures
3811 ! implicit real*8 (a-h,o-z)
3812 ! include 'DIMENSIONS'
3813 use MPI_data !include 'COMMON.INFO'
3818 ! include 'COMMON.CONTROL'
3819 ! include 'COMMON.CHAIN'
3820 ! include 'COMMON.DBASE'
3821 ! include 'COMMON.INTERACT'
3822 ! include 'COMMON.VAR'
3823 ! include 'COMMON.THREAD'
3824 ! include 'COMMON.FFIELD'
3825 ! include 'COMMON.SBRIDGE'
3826 ! include 'COMMON.HEADER'
3827 ! include 'COMMON.IOUNITS'
3828 ! include 'COMMON.TIME1'
3829 ! include 'COMMON.CONTACTS'
3830 ! include 'COMMON.MCM'
3831 ! include 'COMMON.NAMES'
3833 integer :: ThreadId,ThreadType,Kwita
3835 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3836 real(kind=8) :: przes(3),obr(3,3)
3837 real(kind=8) :: time_for_thread
3838 logical :: found_pattern,non_conv
3839 character(len=32) :: head_pdb
3840 real(kind=8) :: energia(0:n_ene)
3841 integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3843 real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3845 n_ene_comp=nprint_ene
3850 if (me.eq.king) then
3869 ave_time_for_thread=0.0D0
3870 max_time_for_thread=0.0D0
3871 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3872 nthread=nexcl+nthread
3873 do ithread=1,nthread
3874 found_pattern=.false.
3876 do while (.not.found_pattern)
3878 if (itrial.gt.1000) then
3879 write (iout,'(/a/)') 'Too many attempts to find pattern.'
3882 call recv_stop_sig(Kwita)
3883 call send_stop_sig(-3)
3887 ! Find long enough chain in the database
3889 nres_t=nres_base(1,ii)
3890 ! Select the starting position to thread.
3891 print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3892 nres_t,' nres0=',nres0
3893 if (nres_t.ge.nres0) then
3894 ist=iran_num(0,nres_t-nres0)
3896 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3897 if (Kwita.lt.0) then
3898 write (iout,*) 'Stop signal received. Terminating.'
3899 write (*,*) 'Stop signal received. Terminating.'
3901 write (*,*) 'ithread=',ithread,' nthread=',nthread
3904 call pattern_receive
3907 if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3909 found_pattern=.true.
3911 ! If this point is reached, the pattern has not yet been examined.
3913 ! print *,'found_pattern:',found_pattern
3919 if (Kwita.eq.0) call recv_stop_sig(Kwita)
3920 if (Kwita.lt.0) then
3921 write (iout,*) 'Stop signal received. Terminating.'
3923 write (*,*) 'ithread=',ithread,' nthread=',nthread
3929 ipatt(2,ithread)=ist
3931 write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3932 'Processor:',me,' Attempt:',ithread,&
3933 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3934 ' start at res.',ist+1
3935 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3936 ' Attempt:',ithread,&
3937 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3938 ' start at res.',ist+1
3940 write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3941 'Attempt:',ithread,&
3942 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3943 ' start at res.',ist+1
3944 write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3945 'Attempt:',ithread,&
3946 ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3947 ' start at res.',ist+1
3950 ! Copy coordinates from the database.
3954 c(j,i)=cart_base(j,i+ist,ii)
3957 !d write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
3959 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3961 !d write (iout,'(a,f10.5)')
3962 !d & 'Initial RMS deviation from reference structure:',rms
3963 if (itype(nres,1).eq.ntyp1) then
3965 dcj=c(j,nres-2)-c(j,nres-3)
3966 c(j,nres)=c(j,nres-1)+dcj
3967 c(j,2*nres)=c(j,nres)
3970 if (itype(1,1).eq.ntyp1) then
3977 call int_from_cart(.false.,.false.)
3978 !d print *,'Exit INT_FROM_CART.'
3979 !d print *,'nhpb=',nhpb
3984 ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3986 ! stop 'End generate'
3987 ! Generate SC conformations.
3991 !d print *,'Processor:',me,': exit GEN_SIDE.'
3993 !d print *,'Exit GEN_SIDE.'
3995 ! Calculate initial energy.
3997 call etotal(energia)
4000 ener0(i,ithread)=energia(i)
4002 ener0(n_ene_comp+1,ithread)=energia(0)
4004 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4005 ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
4007 ener0(n_ene_comp+2,ithread)=rms
4008 ener0(n_ene_comp+4,ithread)=frac
4009 ener0(n_ene_comp+5,ithread)=frac_nn
4011 ener0(n_ene_comp+3,ithread)=0.0d0
4014 print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
4016 print*,'ithread=',ithread,' Start REGULARIZE.'
4019 call regularize(nct-nnt+1,etot,rms,&
4020 cart_base(1,ist+nnt,ipattern),iretcode)
4022 time_for_thread=curr_tim1-curr_tim
4023 ave_time_for_thread= &
4024 ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4025 if (time_for_thread.gt.max_time_for_thread) &
4026 max_time_for_thread=time_for_thread
4028 print *,'Processor',me,': Exit REGULARIZE.'
4029 if (WhatsUp.eq.2) then
4031 'Sufficient number of confs. collected. Terminating.'
4034 else if (WhatsUp.eq.-1) then
4036 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4037 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4038 call send_stop_sig(-2)
4040 else if (WhatsUp.eq.-2) then
4042 write (iout,*) 'Timeup signal received. Terminating.'
4044 else if (WhatsUp.eq.-3) then
4046 write (iout,*) 'Error stop signal received. Terminating.'
4050 print *,'Exit REGULARIZE.'
4051 if (iretcode.eq.11) then
4052 write (iout,'(/a/)') &
4053 '******* Allocated time exceeded in SUMSL. The program will stop.'
4058 head_pdb=titel(:24)//':'//str_nam(ipattern)
4059 if (outpdb) call pdbout(etot,head_pdb,ipdb)
4060 if (outmol2) call mol2out(etot,head_pdb)
4062 call briefout(ithread,etot)
4064 link_end=min0(link_end,nss)
4065 write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4067 call etotal(energia)
4068 ! call enerprint(energia(0))
4071 !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4072 !d write (iout,'(a,f10.5)')
4073 !d & 'RMS deviation from reference structure:',dsqrt(rms)
4075 ener(i,ithread)=energia(i)
4077 ener(n_ene_comp+1,ithread)=energia(0)
4078 ener(n_ene_comp+3,ithread)=rms
4080 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4081 ener(n_ene_comp+2,ithread)=rms
4082 ener(n_ene_comp+4,ithread)=frac
4083 ener(n_ene_comp+5,ithread)=frac_nn
4085 call write_stat_thread(ithread,ipattern,ist)
4086 ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4087 ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4088 ! & (ener(k,ithread),k=12,14)
4090 if (me.eq.king) then
4092 call pattern_receive
4093 call receive_MCM_info
4094 if (nacc_tot.ge.nthread) then
4096 'Sufficient number of conformations collected nacc_tot=',&
4097 nacc_tot,'. Stopping other processors and terminating.'
4099 'Sufficient number of conformations collected nacc_tot=',&
4100 nacc_tot,'. Stopping other processors and terminating.'
4101 call recv_stop_sig(Kwita)
4102 if (Kwita.eq.0) call send_stop_sig(-1)
4107 call send_MCM_info(2)
4110 if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4111 write (iout,'(/2a)') &
4112 '********** There would be not enough time for another thread. ',&
4113 'The program will stop.'
4115 '********** There would be not enough time for another thread. ',&
4116 'The program will stop.'
4117 write (iout,'(a,1pe14.4/)') &
4118 'Elapsed time for last threading step: ',time_for_thread
4121 call recv_stop_sig(Kwita)
4122 call send_stop_sig(-2)
4127 write (iout,'(a,1pe14.4)') &
4128 'Elapsed time for this threading step: ',time_for_thread
4131 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4132 if (Kwita.lt.0) then
4133 write (iout,*) 'Stop signal received. Terminating.'
4134 write (*,*) 'Stop signal received. Terminating.'
4136 write (*,*) 'nthread=',nthread,' ithread=',ithread
4142 call send_stop_sig(-1)
4146 ! Any messages left for me?
4147 call pattern_receive
4148 if (Kwita.eq.0) call recv_stop_sig(Kwita)
4150 call write_thread_summary
4152 if (king.eq.king) then
4154 do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4157 call recv_stop_sig(Kwita)
4158 call receive_MCM_info
4161 call receive_thread_results(iproc)
4163 call write_thread_summary
4165 call send_thread_results
4169 end subroutine thread_seq
4170 !-----------------------------------------------------------------------------
4173 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4175 use random, only:iran_num
4176 ! implicit real*8 (a-h,o-z)
4177 ! include 'DIMENSIONS'
4178 ! include 'COMMON.CHAIN'
4179 ! include 'COMMON.DBASE'
4180 ! include 'COMMON.INTERACT'
4181 ! include 'COMMON.VAR'
4182 ! include 'COMMON.THREAD'
4183 ! include 'COMMON.FFIELD'
4184 ! include 'COMMON.SBRIDGE'
4185 ! include 'COMMON.HEADER'
4186 ! include 'COMMON.GEO'
4187 ! include 'COMMON.IOUNITS'
4188 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
4189 !el integer :: icall
4190 !el common /srutu/ icall
4191 real(kind=8) :: energia(0:n_ene)
4192 logical :: glycine,fail
4193 integer :: i,maxsample,link_end0,ind_sc,isample
4194 real(kind=8) :: alph0,omeg0,e1,e0
4198 link_end=min0(link_end,nss)
4200 if (itype(i,1).ne.10) then
4201 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
4202 call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1)
4206 call etotal(energia)
4208 do isample=1,maxsample
4209 ! Choose a non-glycine side chain.
4212 ind_sc=iran_num(nnt,nct)
4213 glycine=(itype(ind_sc,1).eq.10)
4217 call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
4218 omeg(ind_sc),fail,1)
4220 call etotal(energia)
4221 !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4222 !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4223 !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4234 end subroutine sc_conf
4235 !-----------------------------------------------------------------------------
4237 !-----------------------------------------------------------------------------
4238 logical function check_var(var,info)
4242 ! implicit real*8 (a-h,o-z)
4243 ! include 'DIMENSIONS'
4244 ! include 'COMMON.VAR'
4245 ! include 'COMMON.IOUNITS'
4246 ! include 'COMMON.GEO'
4247 ! include 'COMMON.SETUP'
4248 real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
4249 integer,dimension(3) :: info
4253 do i=nphi+ntheta+1,nphi+ntheta+nside
4254 ! Check the side chain "valence" angles alpha
4255 if (var(i).lt.1.0d-7) then
4256 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4257 write (iout,*) 'Processor',me,'received bad variables!!!!'
4258 write (iout,*) 'Variables'
4259 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4260 write (iout,*) 'Continuing calculations at this point',&
4261 ' could destroy the results obtained so far... ABORTING!!!!!!'
4262 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4263 'valence angle alpha',i-nphi-ntheta,var(i),&
4264 'n it',info(1),info(2),'mv ',info(3)
4265 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4266 write (*,*) 'Processor',me,'received bad variables!!!!'
4267 write (*,*) 'Variables'
4268 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4269 write (*,*) 'Continuing calculations at this point',&
4270 ' could destroy the results obtained so far... ABORTING!!!!!!'
4271 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4272 'valence angle alpha',i-nphi-ntheta,var(i),&
4273 'n it',info(1),info(2),'mv ',info(3)
4278 ! Check the backbone "valence" angles theta
4279 do i=nphi+1,nphi+ntheta
4280 if (var(i).lt.1.0d-7) then
4281 write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4282 write (iout,*) 'Processor',me,'received bad variables!!!!'
4283 write (iout,*) 'Variables'
4284 write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4285 write (iout,*) 'Continuing calculations at this point',&
4286 ' could destroy the results obtained so far... ABORTING!!!!!!'
4287 write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4288 'valence angle theta',i-nphi,var(i),&
4289 'n it',info(1),info(2),'mv ',info(3)
4290 write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4291 write (*,*) 'Processor',me,'received bad variables!!!!'
4292 write (*,*) 'Variables'
4293 write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4294 write (*,*) 'Continuing calculations at this point',&
4295 ' could destroy the results obtained so far... ABORTING!!!!!!'
4296 write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4297 'valence angle theta',i-nphi,var(i),&
4298 'n it',info(1),info(2),'mv ',info(3)
4304 end function check_var
4305 !-----------------------------------------------------------------------------
4307 !-----------------------------------------------------------------------------
4308 subroutine distfit(debug,maxit)
4310 use geometry_data, only: phi
4313 ! implicit real*8 (a-h,o-z)
4314 ! include 'DIMENSIONS'
4315 ! include 'COMMON.CHAIN'
4316 ! include 'COMMON.VAR'
4317 ! include 'COMMON.IOUNITS'
4318 ! include 'COMMON.DISTFIT'
4319 integer :: i,maxit,MAXMAR,IT,IMAR
4320 real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres)
4321 logical :: debug,sing
4322 real(kind=8) :: TOL,RL,F0,AIN,F1
4324 !input------------------------------------
4326 ! NY=((NRES-4)*(NRES-5))/2
4327 !input------------------------------------
4333 CALL TRANSFER(NRES,phi,phiold)
4337 !d WRITE (IOUT,*) 'DISTFIT: F0=',F0
4353 CALL TRANSFER(NX,XX,X)
4354 CALL BANACH(NX,NRES,H,X,sing)
4359 IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4361 WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4362 WRITE (IOUT,*) 'IT=',it,'F=',F0
4367 phi(I)=phiold(I)+mask(i)*X(I-3)
4372 !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4374 CALL TRANSFER(NRES,phi,phiold)
4377 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4379 WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4380 WRITE (IOUT,*) 'IT=',it,'F=',F1
4386 WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4387 WRITE (IOUT,*) 'IT=',it,'F=',F0
4388 CALL TRANSFER(NRES,phiold,phi)
4391 !d write (iout,*) "it",it," imar",imar," f0",f0
4393 WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4395 end subroutine distfit
4396 !-----------------------------------------------------------------------------
4397 real(kind=8) function RDIF()
4400 use geometry, only: dist
4401 ! implicit real*8 (a-h,o-z)
4402 ! include 'DIMENSIONS'
4403 ! include 'COMMON.CHAIN'
4404 ! include 'COMMON.DISTFIT'
4406 real(kind=8) :: suma,DIJ
4415 if (w(ind).ne.0.0) then
4417 suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4419 ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4427 !-----------------------------------------------------------------------------
4432 use geometry, only:dist
4433 ! implicit real*8 (a-h,o-z)
4434 ! include 'DIMENSIONS'
4435 ! include 'COMMON.CHAIN'
4436 ! include 'COMMON.DISTFIT'
4437 ! include 'COMMON.GEO'
4438 integer :: i,j,k,l,I1,I2,IND
4439 real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4452 R13(K)=C(K,J)-C(K,I1)
4456 R24(L)=C(L,K)-C(L,I2)
4458 IND=((J-1)*(2*NRES-J-6))/2+K-3
4459 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4460 PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4461 PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4462 DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4467 end subroutine RDERIV
4468 !-----------------------------------------------------------------------------
4472 ! implicit real*8 (a-h,o-z)
4473 ! include 'DIMENSIONS'
4474 ! include 'COMMON.CHAIN'
4475 ! include 'COMMON.DISTFIT'
4477 real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4485 XI=XI+BKIWK*(D0(K)-DDD(K))
4493 HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4500 end subroutine HEVAL
4501 !-----------------------------------------------------------------------------
4502 subroutine VEC(I,J,U)
4504 use geometry_data, only: C
4505 ! Find the unit vector from atom (I) to atom (J). Store in U.
4507 ! implicit real*8 (a-h,o-z)
4508 ! include 'DIMENSIONS'
4509 ! include 'COMMON.CHAIN'
4511 real(kind=8),DIMENSION(3) :: U
4512 real(kind=8) :: ANORM,UK
4526 !-----------------------------------------------------------------------------
4527 subroutine TRANSFER(N,X1,X2)
4529 ! implicit real*8 (a-h,o-z)
4530 ! include 'DIMENSIONS'
4532 real(kind=8),DIMENSION(N) :: X1,X2
4536 end subroutine TRANSFER
4537 !-----------------------------------------------------------------------------
4538 !-----------------------------------------------------------------------------
4539 subroutine alloc_compare_arrays
4541 maxres22=nres*(nres+1)/2
4543 ! common /struct/ in io_common: read_threadbase
4544 ! allocate(cart_base !(3,maxres_base,maxseq)
4545 ! allocate(nres_base !(3,maxseq)
4546 ! allocate(str_nam !(maxseq)
4548 ! COMMON /c_frag/ in io_conf: readpdb
4549 if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4550 if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4552 if(.not.allocated(w)) allocate(w(maxres22),d0(maxres22)) !(maxres22)
4554 !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4555 if (.not.allocated(DDD)) allocate(DDD(maxres22)) !(maxres22)
4556 if (.not.allocated(H)) allocate(H(nres,nres)) !(MAXRES,MAXRES)
4557 if (.not.allocated(XX)) allocate(XX(nres)) !(MAXRES)
4559 if (.not.allocated(mask)) allocate(mask(nres)) !(maxres)
4562 if (.not.allocated(iexam))allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4564 if (.not.allocated(ener0)) allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4567 end subroutine alloc_compare_arrays
4568 !-----------------------------------------------------------------------------
4570 !-----------------------------------------------------------------------------