module compare !----------------------------------------------------------------------------- use io_units use names use geometry_data use energy_data use control_data #if !defined(WHAM_RUN) && !defined(CLUSTER) use compare_data use io_base use io_config use geometry use energy use control, only: hpb_partition use minim_data use minimm, only: sc_move, minimize #endif implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains #if !defined(WHAM_RUN) && !defined(CLUSTER) !----------------------------------------------------------------------------- ! contact.f !----------------------------------------------------------------------------- subroutine contact(lprint,ncont,icont,co) use geometry, only:dist ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' real(kind=8) :: facont=1.569D0 ! facont = (2/(1-sqrt(1-1/4)))**(1/6) integer :: ncont integer,dimension(2,12*nres) :: icont!(2,12*nres) !(2,maxcont) (maxcont=12*maxres) logical :: lprint !el local variables real(kind=8) :: co,rcomp integer :: kkk,i,j,i1,i2,it1,it2,iti,itj ncont=0 kkk=3 do i=nnt+kkk,nct iti=iabs(itype(i)) do j=nnt,i-kkk itj=iabs(itype(j)) if (ipot.ne.4) then ! rcomp=sigmaii(iti,itj)+1.0D0 rcomp=facont*sigmaii(iti,itj) else ! rcomp=sigma(iti,itj)+1.0D0 rcomp=facont*sigma(iti,itj) endif ! rcomp=6.5D0 ! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j) if (dist(nres+i,nres+j).lt.rcomp) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j endif enddo enddo if (lprint) then write (iout,'(a)') 'Contact map:' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo endif co = 0.0d0 do i=1,ncont co = co + dfloat(iabs(icont(1,i)-icont(2,i))) enddo co = co / (nres*ncont) return end subroutine contact !----------------------------------------------------------------------------- real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: ncont,ncont_ref integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres) !el local variables integer :: i,j,nmatch nmatch=0 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont(1,i),i=1,ncont) ! write (iout,'(20i4)') (icont(2,i),i=1,ncont) do i=1,ncont do j=1,ncont_ref if (icont(1,i).eq.icont_ref(1,j) .and. & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1 enddo enddo ! print *,' nmatch=',nmatch ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref)) contact_fract=dfloat(nmatch)/dfloat(ncont_ref) return end function contact_fract !----------------------------------------------------------------------------- real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' integer :: ncont,ncont_ref integer,dimension(2,12*nres) :: icont,icont_ref !(2,12*nres) (2,maxcont) (maxcont=12*maxres) !el local variables integer :: i,j,nmatch nmatch=0 ! print *,'ncont=',ncont,' ncont_ref=',ncont_ref ! write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref) ! write (iout,'(20i4)') (icont(1,i),i=1,ncont) ! write (iout,'(20i4)') (icont(2,i),i=1,ncont) do i=1,ncont do j=1,ncont_ref if (icont(1,i).eq.icont_ref(1,j) .and. & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1 enddo enddo ! print *,' nmatch=',nmatch ! contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref)) contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont) return end function contact_fract_nn !----------------------------------------------------------------------------- subroutine hairpin(lprint,nharp,iharp) use geometry, only:dist ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' integer :: ncont integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres) integer :: nharp integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) logical :: lprint,not_done real(kind=8) :: rcomp=6.0d0 !el local variables integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1 ! allocate(icont(2,12*nres)) ncont=0 kkk=0 ! print *,'nnt=',nnt,' nct=',nct do i=nnt,nct-3 do k=1,3 c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1)) enddo do j=i+2,nct-1 do k=1,3 c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1)) enddo if (dist(2*nres+1,2*nres+2).lt.rcomp) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j endif enddo enddo if (lprint) then write (iout,'(a)') 'PP contact map:' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4)') & i,restyp(it1),i1,restyp(it2),i2 enddo endif ! finding hairpins nharp=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then ! write (iout,*) "found turn at ",i1,j1 ii1=i1 jj1=j1 not_done=.true. do while (not_done) i1=i1-1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10 enddo not_done=.false. 10 continue ! write (iout,*) i1,j1,not_done enddo i1=i1+1 j1=j1-1 if (j1-i1.gt.4) then nharp=nharp+1 iharp(1,nharp)=i1 iharp(2,nharp)=j1 iharp(3,nharp)=ii1 iharp(4,nharp)=jj1 ! write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4) endif endif enddo ! do i=1,nharp ! write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4) ! enddo if (lprint) then write (iout,*) "Hairpins:" do i=1,nharp i1=iharp(1,i) j1=iharp(2,i) ii1=iharp(3,i) jj1=iharp(4,i) write (iout,*) write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1) write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1) ! do k=jj1,j1,-1 ! write (iout,'(a,i3,$)') restyp(itype(k)),k ! enddo enddo endif return end subroutine hairpin !----------------------------------------------------------------------------- ! elecont.f !----------------------------------------------------------------------------- subroutine elecont(lprint,ncont,icont) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.LOCAL' ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' logical :: lprint real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_ real(kind=8) :: ael6_i,ael3_i real(kind=8),dimension(2,2) :: app_,bpp_,rpp_ integer :: ncont integer,dimension(2,12*nres) :: icont !(2,12*nres)(2,maxcont) (maxcont=12*maxres) real(kind=8),dimension(12*nres) :: econt !(maxcont) !el local variables integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2 real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb real(kind=8) :: xmedi,ymedi,zmedi real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,& evdwij,el1,el2,eesij,ene ! ! Load the constants of peptide bond - peptide bond interactions. ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g. ! proline) - determined by averaging ECEPP energy. ! ! as of 7/06/91. ! ! data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/ data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/ !el allocate(econt(12*nres)) !(maxcont) elcutoff = -0.3d0 elecutoff_14 = -0.5d0 if (lprint) write (iout,'(a)') & "Constants of electrostatic interaction energy expression." do i=1,2 do j=1,2 rri=rpp_(i,j)**6 app_(i,j)=epp(i,j)*rri*rri bpp_(i,j)=-2.0*epp(i,j)*rri ael6_(i,j)=elpp_6(i,j)*4.2**6 ael3_(i,j)=elpp_3(i,j)*4.2**3 if (lprint) & write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),& ael3_(i,j) enddo enddo ncont=0 ees=0.0 evdw=0.0 do 1 i=nnt,nct-2 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1 xi=c(1,i) yi=c(2,i) zi=c(3,i) dxi=c(1,i+1)-c(1,i) dyi=c(2,i+1)-c(2,i) dzi=c(3,i+1)-c(3,i) xmedi=xi+0.5*dxi ymedi=yi+0.5*dyi zmedi=zi+0.5*dzi do 4 j=i+2,nct-1 if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4 iteli=itel(i) itelj=itel(j) if (j.eq.i+2 .and. itelj.eq.2) iteli=2 if (iteli.eq.2 .and. itelj.eq.2) goto 4 aaa=app_(iteli,itelj) bbb=bpp_(iteli,itelj) ael6_i=ael6_(iteli,itelj) ael3_i=ael3_(iteli,itelj) dxj=c(1,j+1)-c(1,j) dyj=c(2,j+1)-c(2,j) dzj=c(3,j+1)-c(3,j) xj=c(1,j)+0.5*dxj-xmedi yj=c(2,j)+0.5*dyj-ymedi zj=c(3,j)+0.5*dzj-zmedi rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) r3ij=rrmij*rmij r6ij=r3ij*r3ij vrmij=vblinv*rmij cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2 cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij fac=cosa-3.0*cosb*cosg ev1=aaa*r6ij*r6ij ev2=bbb*r6ij fac3=ael6_i*r6ij fac4=ael3_i*r3ij evdwij=ev1+ev2 el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg)) el2=fac4*fac eesij=el1+el2 if (j.gt.i+2 .and. eesij.le.elcutoff .or. & j.eq.i+2 .and. eesij.le.elecutoff_14) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j econt(ncont)=eesij endif ees=ees+eesij evdw=evdw+evdwij 4 continue 1 continue if (lprint) then write (iout,*) 'Total average electrostatic energy: ',ees write (iout,*) 'VDW energy between peptide-group centers: ',evdw write (iout,*) write (iout,*) 'Electrostatic contacts before pruning: ' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif ! For given residues keep only the contacts with the greatest energy. i=0 do while (i.lt.ncont) i=i+1 ene=econt(i) ic1=icont(1,i) ic2=icont(2,i) j=i do while (j.lt.ncont) j=j+1 if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then ! write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2, ! & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then if (ic1.eq.icont(1,j)) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) & .and. iabs(icont(1,k)-ic1).le.2 .and. & econt(k).lt.econt(j) ) goto 21 enddo else if (ic2.eq.icont(2,j) ) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) & .and. iabs(icont(2,k)-ic2).le.2 .and. & econt(k).lt.econt(j) ) goto 21 enddo endif ! Remove ith contact do k=i+1,ncont icont(1,k-1)=icont(1,k) icont(2,k-1)=icont(2,k) econt(k-1)=econt(k) enddo i=i-1 ncont=ncont-1 ! write (iout,*) "ncont",ncont ! do k=1,ncont ! write (iout,*) icont(1,k),icont(2,k) ! enddo goto 20 else if (econt(j).gt.ene .and. ic2.ne.ic1+2) & then if (ic1.eq.icont(1,j)) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. & econt(k).lt.econt(i) ) goto 21 enddo else if (ic2.eq.icont(2,j) ) then do k=1,ncont if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. & econt(k).lt.econt(i) ) goto 21 enddo endif ! Remove jth contact do k=j+1,ncont icont(1,k-1)=icont(1,k) icont(2,k-1)=icont(2,k) econt(k-1)=econt(k) enddo ncont=ncont-1 ! write (iout,*) "ncont",ncont ! do k=1,ncont ! write (iout,*) icont(1,k),icont(2,k) ! enddo j=j-1 endif endif 21 continue enddo 20 continue enddo if (lprint) then write (iout,*) write (iout,*) 'Electrostatic contacts after pruning: ' do i=1,ncont i1=icont(1,i) i2=icont(2,i) it1=itype(i1) it2=itype(i2) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif return end subroutine elecont !----------------------------------------------------------------------------- subroutine secondary2(lprint) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.VAR' ! include 'COMMON.GEO' ! include 'COMMON.CONTROL' integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,& iii1,jjj1 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres) integer,dimension(nres,4) :: isec !(maxres,4) integer,dimension(nres) :: nsec !(maxres) logical :: lprint,not_done !,freeres real(kind=8) :: p1,p2 !el external freeres !el allocate(icont(2,12*nres),isec(nres,4),nsec(nres)) if(.not.dccart) call chainbuild if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3) !d call write_pdb(99,'sec structure',0d0) ncont=0 nbfrag=0 nhfrag=0 do i=1,nres isec(i,1)=0 isec(i,2)=0 nsec(i)=0 enddo call elecont(lprint,ncont,icont) ! finding parallel beta !d write (iout,*) '------- looking for parallel beta -----------' nbeta=0 nstrand=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1 not_done=.true. do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. & freeres(i1,j1,nsec,isec)) goto 5 enddo not_done=.false. 5 continue !d write (iout,*) i1,j1,not_done enddo j1=j1-1 i1=i1-1 if (i1-ii1.gt.1) then ii1=max0(ii1-1,1) jj1=max0(jj1-1,1) nbeta=nbeta+1 if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',& nbeta,ii1,i1,jj1,j1 nbfrag=nbfrag+1 bfrag(1,nbfrag)=ii1+1 bfrag(2,nbfrag)=i1+1 bfrag(3,nbfrag)=jj1+1 bfrag(4,nbfrag)=min0(j1+1,nres) do ij=ii1,i1 nsec(ij)=nsec(ij)+1 isec(ij,nsec(ij))=nbeta enddo do ij=jj1,j1 nsec(ij)=nsec(ij)+1 isec(ij,nsec(ij))=nbeta enddo if(lprint) then nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" endif nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1 endif endif endif enddo ! finding alpha or 310 helix nhelix=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) p1=phi(i1+2)*rad2deg p2=0.0 if (j1+2.le.nres) p2=phi(j1+2)*rad2deg if (j1.eq.i1+3 .and. & ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. & ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2 !o if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2 ii1=i1 jj1=j1 if (nsec(ii1).eq.0) then not_done=.true. else not_done=.false. endif do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10 enddo not_done=.false. 10 continue p1=phi(i1+2)*rad2deg p2=phi(j1+2)*rad2deg if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) & not_done=.false. !d enddo j1=j1+1 if (j1-ii1.gt.5) then nhelix=nhelix+1 !d nhfrag=nhfrag+1 hfrag(1,nhfrag)=ii1 hfrag(2,nhfrag)=j1 do ij=ii1,j1 nsec(ij)=-1 enddo if (lprint) then write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1 if (nhelix.le.9) then write(12,'(a17,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" else write(12,'(a17,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" endif endif endif endif enddo if (nhelix.gt.0.and.lprint) then write(12,'(a26,$)') "DefPropRes 'helix' 'helix1" do i=2,nhelix if (nhelix.le.9) then write(12,'(a8,i1,$)') " | helix",i else write(12,'(a8,i2,$)') " | helix",i endif enddo write(12,'(a1)') "'" endif ! finding antiparallel beta !d write (iout,*) '--------- looking for antiparallel beta ---------' do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (freeres(i1,j1,nsec,isec)) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1 not_done=.true. do while (not_done) i1=i1+1 j1=j1-1 do j=1,ncont if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. & freeres(i1,j1,nsec,isec)) goto 6 enddo not_done=.false. 6 continue !d write (iout,*) i1,j1,not_done enddo i1=i1-1 j1=j1+1 if (i1-ii1.gt.1) then nbfrag=nbfrag+1 bfrag(1,nbfrag)=ii1 bfrag(2,nbfrag)=min0(i1+1,nres) bfrag(3,nbfrag)=min0(jj1+1,nres) bfrag(4,nbfrag)=j1 nbeta=nbeta+1 iii1=max0(ii1-1,1) do ij=iii1,i1 nsec(ij)=nsec(ij)+1 if (nsec(ij).le.2) then isec(ij,nsec(ij))=nbeta endif enddo jjj1=max0(j1-1,1) do ij=jjj1,jj1 nsec(ij)=nsec(ij)+1 if (nsec(ij).le.2 .and. nsec(ij).gt.0) then isec(ij,nsec(ij))=nbeta endif enddo if (lprint) then write (iout,'(a,i3,4i4)')'antiparallel beta',& nbeta,ii1-1,i1,jj1,j1-1 nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" endif nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2 endif endif endif enddo if (nstrand.gt.0.and.lprint) then write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1" do i=2,nstrand if (i.le.9) then write(12,'(a9,i1,$)') " | strand",i else write(12,'(a9,i2,$)') " | strand",i endif enddo write(12,'(a1)') "'" endif if (lprint) then write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'" write(12,'(a20)') "XMacStand ribbon.mac" write(iout,*) 'UNRES seq:' do j=1,nbfrag write(iout,*) 'beta ',(bfrag(i,j),i=1,4) enddo do j=1,nhfrag write(iout,*) 'helix ',(hfrag(i,j),i=1,2) enddo endif return end subroutine secondary2 #endif !----------------------------------------------------------------------------- logical function freeres(i,j,nsec,isec) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' integer,dimension(nres,4) :: isec !(maxres,4) integer,dimension(nres) :: nsec !(maxres) !el local variables integer :: i,j,k,l freeres=.false. #ifndef WHAM_RUN if (nsec(i).lt.0.or.nsec(j).lt.0) return #endif if (nsec(i).gt.1.or.nsec(j).gt.1) return do k=1,nsec(i) do l=1,nsec(j) if (isec(i,k).eq.isec(j,l)) return enddo enddo freeres=.true. return end function freeres !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) !el implicit none integer :: length,itypea(length),itypeb(length) integer :: i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end function seq_comp #ifndef WHAM_RUN !----------------------------------------------------------------------------- ! rmsd.F !----------------------------------------------------------------------------- subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.CONTACTS' ! include 'COMMON.IOUNITS' real(kind=8) :: przes(3),obr(3,3) logical :: non_conv,lprn real(kind=8) :: rms,frac,frac_nn,co ! call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes, ! & obr,non_conv) ! rms=dsqrt(rms) call rmsd(rms) !elte(iout,*) "rms_nacc before contact" call contact(.false.,ncont,icont,co) frac=contact_fract(ncont,ncont_ref,icont,icont_ref) frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref) if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') & 'RMS deviation from the reference structure:',rms,& ' % of native contacts:',frac*100,& ' % of nonnative contacts:',frac_nn*100,& ' contact order:',co return end subroutine rms_nac_nnc !----------------------------------------------------------------------------- subroutine rmsd(drms) use regularize_, only:fitsq ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.INTERACT' ! include 'COMMON.CONTROL' logical :: non_conv real(kind=8) :: przes(3),obrot(3,3) real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres !el local variables real(kind=8) :: drms,rminroz,roznica integer :: i,j,iatom,kkk,iti,k !el allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres nperm=1 do i=1,symetr nperm=nperm*i enddo iatom=0 rminroz=100d2 ! print *,"nz_start",nz_start," nz_end",nz_end ! if (symetr.le.1) then do kkk=1,nperm ! do i=nz_start,nz_end ! iatom=iatom+1 ! iti=itype(i) ! do k=1,3 ! ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup) ! crefcopy(k,iatom,kkk)=cref(k,i,kkk) ! enddo ! if (iz_sc.eq.1.and.iti.ne.10) then ! iatom=iatom+1 ! do k=1,3 ! ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup) ! crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk) ! enddo ! endif ! enddo ! else ! do kkk=1,nperm iatom=0 do i=nz_start,nz_end iatom=iatom+1 iti=itype(i) do k=1,3 ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup) crefcopy(k,iatom)=cref(k,i,kkk) enddo if (iz_sc.eq.1.and.iti.ne.10) then iatom=iatom+1 do k=1,3 ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup) crefcopy(k,iatom)=cref(k,nres+i,kkk) enddo endif enddo ! enddo ! endif ! ----- diagnostics ! do kkk=1,nperm ! write (iout,*) 'Ccopy and CREFcopy' ! print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3), ! & (crefcopy(j,k),j=1,3),k=1,iatom) ! write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3), ! & (crefcopy(j,k),j=1,3),k=1,iatom) ! enddo ! ----- end diagnostics ! do kkk=1,nperm call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,& przes,obrot,non_conv) if (non_conv) then print *,'Problems in FITSQ!!! rmsd' write (iout,*) 'Problems in FITSQ!!! rmsd' print *,'Ccopy and CREFcopy' write (iout,*) 'Ccopy and CREFcopy' print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),& (crefcopy(j,k),j=1,3),k=1,iatom) write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),& (crefcopy(j,k),j=1,3),k=1,iatom) #ifdef MPI ! call mpi_abort(mpi_comm_world,ierror,ierrcode) roznica=100.0d10 #else stop #endif endif ! write (iout,*) "roznica", roznica,kkk if (roznica.le.rminroz) rminroz=roznica enddo drms=dsqrt(dabs(rminroz)) ! ---- diagnostics ! write (iout,*) "nperm,symetr", nperm,symetr ! ---- end diagnostics return end subroutine rmsd !----------------------------------------------------------------------------- subroutine rmsd_csa(drms) use regularize_, only:fitsq ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.INTERACT' logical :: non_conv real(kind=8) :: przes(3),obrot(3,3) real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres integer :: kkk,iatom,ierror,ierrcode !el local variables integer ::i,j,k,iti real(kind=8) :: drms,roznica allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres kkk=1 iatom=0 do i=nz_start,nz_end iatom=iatom+1 iti=itype(i) do k=1,3 ccopy(k,iatom)=c(k,i) crefcopy(k,iatom)=crefjlee(k,i) enddo if (iz_sc.eq.1.and.iti.ne.10) then iatom=iatom+1 do k=1,3 ccopy(k,iatom)=c(k,nres+i) crefcopy(k,iatom)=crefjlee(k,nres+i) enddo endif enddo call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,& przes,obrot,non_conv) if (non_conv) then print *,'Problems in FITSQ!!! rmsd_csa' write (iout,*) 'Problems in FITSQ!!! rmsd_csa' print *,'Ccopy and CREFcopy' write (iout,*) 'Ccopy and CREFcopy' print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),& (crefcopy(j,k),j=1,3),k=1,iatom) write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),& (crefcopy(j,k),j=1,3),k=1,iatom) #ifdef MPI call mpi_abort(mpi_comm_world,ierror,ierrcode) #else stop #endif endif drms=dsqrt(dabs(roznica)) return end subroutine rmsd_csa !----------------------------------------------------------------------------- ! test.F !----------------------------------------------------------------------------- subroutine test !el use minim use geometry, only:pinorm use random, only:ran_number,iran_num ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! include 'COMMON.CHAIN' real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: var,var1 !(maxvar) (maxvar=6*maxres) integer :: j1,j2,jr,i,iretcode,nfun,nft_sc logical :: debug,accepted real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,& xxr,xxh debug=.true. !el allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres) call geom_to_var(nvar,var1) call chainbuild call etotal(energy) etot=energy(0) call rmsd(rms) write(iout,*) 'etot=',0,etot,rms call secondary2(.false.) call write_pdb(0,'first structure',etot) j1=13 j2=21 da=180.0*deg2rad temp=3000.0d0 betbol=1.0D0/(1.9858D-3*temp) jr=iran_num(j1,j2) d=ran_number(-pi,pi) ! phi(jr)=pinorm(phi(jr)+d) call chainbuild call etotal(energy) etot0=energy(0) call rmsd(rms) write(iout,*) 'etot=',1,etot0,rms call write_pdb(1,'perturb structure',etot0) do i=2,500,2 jr=iran_num(j1,j2) d=ran_number(-da,da) phiold=phi(jr) phi(jr)=pinorm(phi(jr)+d) call chainbuild call etotal(energy) etot=energy(0) if (etot.lt.etot0) then accepted=.true. else accepted=.false. xxr=ran_number(0.0D0,1.0D0) xxh=betbol*(etot-etot0) if (xxh.lt.50.0D0) then xxh=dexp(-xxh) if (xxh.gt.xxr) accepted=.true. endif endif accepted=.true. ! print *,etot0,etot,accepted if (accepted) then etot0=etot call rmsd(rms) write(iout,*) 'etot=',i,etot,rms call write_pdb(i,'MC structure',etot) ! minimize ! call geom_to_var(nvar,var1) call sc_move(2,nres-1,1,10d0,nft_sc,etot) call geom_to_var(nvar,var) call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun call var_to_geom(nvar,var) call chainbuild call rmsd(rms) write(iout,*) 'etot mcm=',i,etot,rms call write_pdb(i+1,'MCM structure',etot) call var_to_geom(nvar,var1) ! -------- else phi(jr)=phiold endif enddo ! minimize ! call sc_move(2,nres-1,1,10d0,nft_sc,etot) ! call geom_to_var(nvar,var) ! ! call chainbuild ! call write_pdb(998 ,'sc min',etot) ! ! call minimize(etot,var,iretcode,nfun) ! write(iout,*)'------------------------------------------------' ! write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun ! ! call var_to_geom(nvar,var) ! call chainbuild ! call write_pdb(999,'full min',etot) return end subroutine test !----------------------------------------------------------------------------- !el#ifdef MPI subroutine test_n16 !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! include 'COMMON.CHAIN' real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres) integer :: jdata(5) logical :: debug !el local variables integer :: i,ij,ieval,iretcode,nfun real(kind=8) :: etot debug=.true. allocate(var(6*nres),var1(6*nres)) !(maxvar) (maxvar=6*maxres) ! call geom_to_var(nvar,var1) call chainbuild call etotal(energy) etot=energy(0) write(iout,*) nnt,nct,etot call write_pdb(1,'first structure',etot) call secondary2(.true.) do i=1,4 jdata(i)=bfrag(i,2) enddo DO ij=1,4 ieval=0 jdata(5)=ij call var_to_geom(nvar,var1) write(iout,*) 'N16 test',(jdata(i),i=1,5) call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), & ieval,ij) call geom_to_var(nvar,var) if (minim) then !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,& '+ DIST eval',ieval !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,& nfun/(time1-time0),' eval/s' call var_to_geom(nvar,var) call chainbuild call write_pdb(ij*100+99,'full min',etot) endif ENDDO return end subroutine test_n16 !el#endif !----------------------------------------------------------------------------- subroutine test_local ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) integer :: nft_sc real(kind=8) :: etot ! ! allocate(varia(6*nres)) !(maxvar) (maxvar=6*maxres) call chainbuild ! call geom_to_var(nvar,varia) call write_pdb(1,'first structure',0d0) call etotal(energy) etot=energy(0) write(iout,*) nnt,nct,etot write(iout,*) 'calling sc_move' call sc_move(nnt,nct,5,10d0,nft_sc,etot) write(iout,*) nft_sc,etot call write_pdb(2,'second structure',etot) write(iout,*) 'calling local_move' call local_move_init(.false.) call local_move(24,29,20d0,50d0) call chainbuild call write_pdb(3,'third structure',etot) write(iout,*) 'calling sc_move' call sc_move(24,29,5,10d0,nft_sc,etot) write(iout,*) nft_sc,etot call write_pdb(2,'last structure',etot) return end subroutine test_local !----------------------------------------------------------------------------- subroutine test_sc ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' real(kind=8) :: time0,time1,etot real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) integer :: nft_sc ! call chainbuild ! call geom_to_var(nvar,varia) call write_pdb(1,'first structure',0d0) call etotal(energy) etot=energy(0) write(iout,*) nnt,nct,etot write(iout,*) 'calling sc_move' call sc_move(nnt,nct,5,10d0,nft_sc,etot) write(iout,*) nft_sc,etot call write_pdb(2,'second structure',etot) write(iout,*) 'calling sc_move 2nd time' call sc_move(nnt,nct,5,1d0,nft_sc,etot) write(iout,*) nft_sc,etot call write_pdb(3,'last structure',etot) return end subroutine test_sc !----------------------------------------------------------------------------- subroutine bgrow(bstrand,nbstrand,in,ind,new) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' integer,dimension(nres/3,6) :: bstrand !(maxres/3,6) !el local variables integer :: nbstrand,in,ind,new,ishift,i ishift=iabs(bstrand(in,ind+4)-new) print *,'bgrow',bstrand(in,ind+4),new,ishift bstrand(in,ind)=new if(ind.eq.1)then bstrand(nbstrand,5)=bstrand(nbstrand,1) do i=1,nbstrand-1 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN if (bstrand(i,5).lt.bstrand(i,6)) then bstrand(i,5)=bstrand(i,5)-ishift else bstrand(i,5)=bstrand(i,5)+ishift endif ENDIF enddo else bstrand(nbstrand,6)=bstrand(nbstrand,2) do i=1,nbstrand-1 IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN if (bstrand(i,6).lt.bstrand(i,5)) then bstrand(i,6)=bstrand(i,6)-ishift else bstrand(i,6)=bstrand(i,6)+ishift endif ENDIF enddo endif return end subroutine bgrow !----------------------------------------------------------------------------- subroutine test11 use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! ! include 'COMMON.DISTFIT' integer :: if(20,nres),nif,ifa(20) integer :: ibc(0:nres,0:nres),istrand(20) integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0 integer :: itmp(20,nres) real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: varia,vorg !(maxvar) (maxvar=6*maxres) ! logical :: debug,ltest,usedbfrag(nres/3) character(len=50) :: linia ! integer :: betasheet(nres),ibetasheet(nres),nbetasheet integer :: bstrand(nres/3,6),nbstrand real(kind=8) :: etot integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,& in,ind,ifun,nfun,iretcode !------------------------ debug=.true. !------------------------ nbstrand=0 nbetasheet=0 do i=1,nres betasheet(i)=0 ibetasheet(i)=0 enddo call geom_to_var(nvar,vorg) call secondary2(debug) if (nbfrag.le.1) return do i=1,nbfrag usedbfrag(i)=.false. enddo nbetasheet=nbetasheet+1 nbstrand=2 bstrand(1,1)=bfrag(1,1) bstrand(1,2)=bfrag(2,1) bstrand(1,3)=nbetasheet bstrand(1,4)=1 bstrand(1,5)=bfrag(1,1) bstrand(1,6)=bfrag(2,1) do i=bfrag(1,1),bfrag(2,1) betasheet(i)=nbetasheet ibetasheet(i)=1 enddo ! bstrand(2,1)=bfrag(3,1) bstrand(2,2)=bfrag(4,1) bstrand(2,3)=nbetasheet bstrand(2,5)=bfrag(3,1) bstrand(2,6)=bfrag(4,1) if (bfrag(3,1).le.bfrag(4,1)) then bstrand(2,4)=2 do i=bfrag(3,1),bfrag(4,1) betasheet(i)=nbetasheet ibetasheet(i)=2 enddo else bstrand(2,4)=-2 do i=bfrag(4,1),bfrag(3,1) betasheet(i)=nbetasheet ibetasheet(i)=2 enddo endif iused_nbfrag=1 do while (iused_nbfrag.ne.nbfrag) do j=2,nbfrag IF (.not.usedbfrag(j)) THEN write (*,*) j,(bfrag(i,j),i=1,4) do jk=6,1,-1 write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand) enddo write (*,*) '------------------' if (bfrag(3,j).le.bfrag(4,j)) then do i=bfrag(3,j),bfrag(4,j) if(betasheet(i).eq.nbetasheet) then in=ibetasheet(i) do k=bfrag(3,j),bfrag(4,j) betasheet(k)=nbetasheet ibetasheet(k)=in enddo nbstrand=nbstrand+1 usedbfrag(j)=.true. iused_nbfrag=iused_nbfrag+1 do k=bfrag(1,j),bfrag(2,j) betasheet(k)=nbetasheet ibetasheet(k)=nbstrand enddo if (bstrand(in,4).lt.0) then bstrand(nbstrand,1)=bfrag(2,j) bstrand(nbstrand,2)=bfrag(1,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=-nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).lt.bfrag(4,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(4,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)+& (bstrand(in,5)-bfrag(4,j)) endif if(bstrand(in,2).gt.bfrag(3,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(3,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)-& (-bstrand(in,6)+bfrag(3,j)) endif else bstrand(nbstrand,1)=bfrag(1,j) bstrand(nbstrand,2)=bfrag(2,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).gt.bfrag(3,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(3,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)-& (-bstrand(in,5)+bfrag(3,j)) endif if(bstrand(in,2).lt.bfrag(4,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(4,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)+& (bstrand(in,6)-bfrag(4,j)) endif endif goto 11 endif if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then in=ibetasheet(bfrag(1,j)+i-bfrag(3,j)) do k=bfrag(1,j),bfrag(2,j) betasheet(k)=nbetasheet ibetasheet(k)=in enddo nbstrand=nbstrand+1 usedbfrag(j)=.true. iused_nbfrag=iused_nbfrag+1 do k=bfrag(3,1),bfrag(4,1) betasheet(k)=nbetasheet ibetasheet(k)=nbstrand enddo if (bstrand(in,4).lt.0) then bstrand(nbstrand,1)=bfrag(4,j) bstrand(nbstrand,2)=bfrag(3,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=-nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).lt.bfrag(2,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(2,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)+& (bstrand(in,5)-bfrag(2,j)) endif if(bstrand(in,2).gt.bfrag(1,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(1,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)-& (-bstrand(in,6)+bfrag(1,j)) endif else bstrand(nbstrand,1)=bfrag(3,j) bstrand(nbstrand,2)=bfrag(4,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).gt.bfrag(1,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(1,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)-& (-bstrand(in,5)+bfrag(1,j)) endif if(bstrand(in,2).lt.bfrag(2,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(2,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)+& (bstrand(in,6)-bfrag(2,j)) endif endif goto 11 endif enddo else do i=bfrag(4,j),bfrag(3,j) if(betasheet(i).eq.nbetasheet) then in=ibetasheet(i) do k=bfrag(4,j),bfrag(3,j) betasheet(k)=nbetasheet ibetasheet(k)=in enddo nbstrand=nbstrand+1 usedbfrag(j)=.true. iused_nbfrag=iused_nbfrag+1 do k=bfrag(1,j),bfrag(2,j) betasheet(k)=nbetasheet ibetasheet(k)=nbstrand enddo if (bstrand(in,4).lt.0) then bstrand(nbstrand,1)=bfrag(1,j) bstrand(nbstrand,2)=bfrag(2,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).lt.bfrag(3,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(3,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)-& (bstrand(in,5)-bfrag(3,j)) endif if(bstrand(in,2).gt.bfrag(4,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(4,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)+& (-bstrand(in,6)+bfrag(4,j)) endif else bstrand(nbstrand,1)=bfrag(2,j) bstrand(nbstrand,2)=bfrag(1,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=-nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).gt.bfrag(4,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(4,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)+& (-bstrand(in,5)+bfrag(4,j)) endif if(bstrand(in,2).lt.bfrag(3,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(3,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)-& (bstrand(in,6)-bfrag(3,j)) endif endif goto 11 endif if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then in=ibetasheet(bfrag(2,j)-i+bfrag(4,j)) do k=bfrag(1,j),bfrag(2,j) betasheet(k)=nbetasheet ibetasheet(k)=in enddo nbstrand=nbstrand+1 usedbfrag(j)=.true. iused_nbfrag=iused_nbfrag+1 do k=bfrag(4,j),bfrag(3,j) betasheet(k)=nbetasheet ibetasheet(k)=nbstrand enddo if (bstrand(in,4).lt.0) then bstrand(nbstrand,1)=bfrag(4,j) bstrand(nbstrand,2)=bfrag(3,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).lt.bfrag(2,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(2,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)-& (bstrand(in,5)-bfrag(2,j)) endif if(bstrand(in,2).gt.bfrag(1,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(1,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)+& (-bstrand(in,6)+bfrag(1,j)) endif else bstrand(nbstrand,1)=bfrag(3,j) bstrand(nbstrand,2)=bfrag(4,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,4)=-nbstrand bstrand(nbstrand,5)=bstrand(nbstrand,1) bstrand(nbstrand,6)=bstrand(nbstrand,2) if(bstrand(in,1).gt.bfrag(1,j)) then call bgrow(bstrand,nbstrand,in,1,bfrag(1,j)) else bstrand(nbstrand,5)=bstrand(nbstrand,5)+& (-bstrand(in,5)+bfrag(1,j)) endif if(bstrand(in,2).lt.bfrag(2,j)) then call bgrow(bstrand,nbstrand,in,2,bfrag(2,j)) else bstrand(nbstrand,6)=bstrand(nbstrand,6)-& (bstrand(in,6)-bfrag(2,j)) endif endif goto 11 endif enddo endif ENDIF enddo j=2 do while (usedbfrag(j)) j=j+1 enddo nbstrand=nbstrand+1 nbetasheet=nbetasheet+1 bstrand(nbstrand,1)=bfrag(1,j) bstrand(nbstrand,2)=bfrag(2,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,5)=bfrag(1,j) bstrand(nbstrand,6)=bfrag(2,j) bstrand(nbstrand,4)=nbstrand do i=bfrag(1,j),bfrag(2,j) betasheet(i)=nbetasheet ibetasheet(i)=nbstrand enddo ! nbstrand=nbstrand+1 bstrand(nbstrand,1)=bfrag(3,j) bstrand(nbstrand,2)=bfrag(4,j) bstrand(nbstrand,3)=nbetasheet bstrand(nbstrand,5)=bfrag(3,j) bstrand(nbstrand,6)=bfrag(4,j) if (bfrag(3,j).le.bfrag(4,j)) then bstrand(nbstrand,4)=nbstrand do i=bfrag(3,j),bfrag(4,j) betasheet(i)=nbetasheet ibetasheet(i)=nbstrand enddo else bstrand(nbstrand,4)=-nbstrand do i=bfrag(4,j),bfrag(3,j) betasheet(i)=nbetasheet ibetasheet(i)=nbstrand enddo endif iused_nbfrag=iused_nbfrag+1 usedbfrag(j)=.true. 11 continue do jk=6,1,-1 write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand) enddo enddo do i=1,nres if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i) enddo write(*,*) do j=6,1,-1 write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand) enddo !------------------------ nifb=0 do i=1,nbstrand do j=i+1,nbstrand if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. & iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then nifb=nifb+1 ifb(nifb,1)=bstrand(i,4) ifb(nifb,2)=bstrand(j,4) endif enddo enddo write(*,*) do i=1,nifb write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2) enddo do i=1,nbstrand ifa(i)=bstrand(i,4) enddo write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand) nif=iabs(bstrand(1,6)-bstrand(1,5))+1 do j=2,nbstrand if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) & nif=iabs(bstrand(j,6)-bstrand(j,5))+1 enddo write(*,*) nif do i=1,nif do j=1,nbstrand if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6)) if (if(j,i).gt.0) then if(betasheet(if(j,i)).eq.0 .or. & ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0 else if(j,i)=0 endif enddo write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand) enddo ! read (inp,*) (ifa(i),i=1,4) ! do i=1,nres ! read (inp,*,err=20,end=20) (if(j,i),j=1,4) ! enddo ! 20 nif=i-1 stop !------------------------ isa=4 is=2*isa-1 iconf=0 !ccccccccccccccccccccccccccccccccc DO ig=1,is**isa-1 !ccccccccccccccccccccccccccccccccc ii=ig do j=1,is istrand(is-j+1)=int(ii/is**(is-j)) ii=ii-istrand(is-j+1)*is**(is-j) enddo ltest=.true. do k=1,isa istrand(k)=istrand(k)+1 if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1 enddo do k=1,isa do l=1,isa if(istrand(k).eq.istrand(l).and.k.ne.l.or. & istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false. enddo enddo lifb0=1 do m=1,nifb lifb(m)=0 do k=1,isa-1 if( & ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. & ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. & -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. & -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) & lifb(m)=1 enddo lifb0=lifb0*lifb(m) enddo if (mod(isa,2).eq.0) then do k=isa/2+1,isa if (istrand(k).eq.1) ltest=.false. enddo else do k=(isa+1)/2+1,isa if (istrand(k).eq.1) ltest=.false. enddo endif IF (ltest.and.lifb0.eq.1) THEN iconf=iconf+1 call var_to_geom(nvar,vorg) write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa) write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa) write (linia,'(10i3)') (istrand(k),k=1,isa) do i=1,nres do j=1,nres ibc(i,j)=0 enddo enddo do i=1,4 if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then do j=1,nif itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j) enddo else do j=1,nif itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1) enddo endif enddo do i=1,nif write(*,*) (itmp(j,i),j=1,4) enddo do i=1,nif ! ifa(1),ifa(2),ifa(3),ifa(4) ! if(1,i),if(2,i),if(3,i),if(4,i) do k=1,isa-1 ltest=.false. do m=1,nifb if( & ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. & ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. & -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. & -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) & then ltest=.true. goto 110 endif enddo 110 continue if (ltest) then ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1 else ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2 endif ! if (k.lt.3) & ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3 if (k.lt.2) & ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4 enddo enddo !------------------------ ! ! freeze sec.elements ! do i=1,nres mask(i)=1 mask_phi(i)=1 mask_theta(i)=1 mask_side(i)=1 enddo do j=1,nbfrag do i=bfrag(1,j),bfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo if (bfrag(3,j).le.bfrag(4,j)) then do i=bfrag(3,j),bfrag(4,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo else do i=bfrag(4,j),bfrag(3,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo endif enddo do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo enddo mask_r=.true. !------------------------ ! generate constrains ! nhpb0=nhpb call chainbuild ind=0 do i=1,nres-3 do j=i+3,nres ind=ind+1 if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then d0(ind)=DIST(i,j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then d0(ind)=5.0 w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then d0(ind)=11.0 w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then d0(ind)=16.0 w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(i,j).gt.0 ) then d0(ind)=DIST(i,ibc(i,j)) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(j,i).gt.0 ) then d0(ind)=DIST(ibc(j,i),j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else w(ind)=0.0 endif ddd(ind)=d0(ind) enddo enddo call hpb_partition !d-------------------------- write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),& ibc(jhpb(i),ihpb(i)),' --',& ihpb(i),jhpb(i),dhpb(i),i=1,nhpb) !d nhpb=0 !d goto 901 ! ! !el#ifdef MPI call contact_cp_min(varia,ifun,iconf,linia,debug) if (minim) then !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,& '+ DIST eval',ifun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,& nfun/(time1-time0),' eval/s' write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa) call var_to_geom(nvar,varia) call chainbuild call write_pdb(900+iconf,linia,etot) endif !el#endif call etotal(energy) etot=energy(0) call enerprint(energy) !d call intout !d call briefout(0,etot) !d call secondary2(.true.) 901 CONTINUE !test return !ccccccccccccccccccccccccccccccccccc ENDIF ENDDO !ccccccccccccccccccccccccccccccccccc return 10 write (iout,'(a)') 'Error reading test structure.' return end subroutine test11 !----------------------------------------------------------------------------- subroutine test3 use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! ! include 'COMMON.DISTFIT' integer :: if(3,nres),nif integer :: ibc(nres,nres),istrand(20) integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0 real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) ! logical :: debug,ltest character(len=50) :: linia integer :: ieval,i,j,ind,in_pdb,nfun,iretcode real(kind=8) :: etot ! do i=1,nres read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i) enddo 20 nif=i-1 write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),& i=1,nif) !------------------------ call secondary2(debug) !------------------------ do i=1,nres do j=1,nres ibc(i,j)=0 enddo enddo ! ! freeze sec.elements and store indexes for beta constrains ! do i=1,nres mask(i)=1 mask_phi(i)=1 mask_theta(i)=1 mask_side(i)=1 enddo do j=1,nbfrag do i=bfrag(1,j),bfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo if (bfrag(3,j).le.bfrag(4,j)) then do i=bfrag(3,j),bfrag(4,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1 enddo else do i=bfrag(4,j),bfrag(3,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1 enddo endif enddo do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo enddo mask_r=.true. ! ---------------- test -------------- do i=1,nif if (ibc(if(1,i),if(2,i)).eq.-1) then ibc(if(1,i),if(2,i))=if(3,i) ibc(if(1,i),if(3,i))=if(2,i) else if (ibc(if(2,i),if(1,i)).eq.-1) then ibc(if(2,i),if(1,i))=0 ibc(if(1,i),if(2,i))=if(3,i) ibc(if(1,i),if(3,i))=if(2,i) else ibc(if(1,i),if(2,i))=if(3,i) ibc(if(1,i),if(3,i))=if(2,i) endif enddo do i=1,nres do j=1,nres if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j) enddo enddo !------------------------ call chainbuild ind=0 do i=1,nres-3 do j=i+3,nres ind=ind+1 if ( ibc(i,j).eq.-1 ) then d0(ind)=DIST(i,j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(i,j).gt.0 ) then d0(ind)=DIST(i,ibc(i,j)) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else if ( ibc(j,i).gt.0 ) then d0(ind)=DIST(ibc(j,i),j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else w(ind)=0.0 endif enddo enddo call hpb_partition !d-------------------------- write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),& ibc(jhpb(i),ihpb(i)),' --',& ihpb(i),jhpb(i),dhpb(i),i=1,nhpb) linia='dist' debug=.true. in_pdb=7 ! !el#ifdef MPI call contact_cp_min(varia,ieval,in_pdb,linia,debug) if (minim) then !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,& '+ DIST eval',ieval !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,& nfun/(time1-time0),' eval/s' call var_to_geom(nvar,varia) call chainbuild call write_pdb(999,'full min',etot) endif !el#endif call etotal(energy) etot=energy(0) call enerprint(energy) call intout call briefout(0,etot) call secondary2(.true.) return 10 write (iout,'(a)') 'Error reading test structure.' return end subroutine test3 !----------------------------------------------------------------------------- subroutine test__ !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! ! include 'COMMON.DISTFIT' integer :: if(2,2),ind integer :: iff(nres) real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,& theta1,phi1,alph1,omeg1 !(maxres) real(kind=8),dimension(6*nres) :: varia,varia2 !(maxvar) (maxvar=6*maxres) ! integer :: i,j,nn,ifun,iretcode,nfun real(kind=8) :: etot nn=0 read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2) write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2) read (inp,*,err=10,end=10) (theta2(i),i=3,nres) read (inp,*,err=10,end=10) (phi2(i),i=4,nres) read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1) read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1) do i=1,nres theta2(i)=deg2rad*theta2(i) phi2(i)=deg2rad*phi2(i) alph2(i)=deg2rad*alph2(i) omeg2(i)=deg2rad*omeg2(i) enddo do i=1,nres theta1(i)=theta(i) phi1(i)=phi(i) alph1(i)=alph(i) omeg1(i)=omeg(i) enddo do i=1,nres mask(i)=1 enddo !------------------------ do i=1,nres iff(i)=0 enddo do j=1,2 do i=if(j,1),if(j,2) iff(i)=1 enddo enddo call chainbuild call geom_to_var(nvar,varia) call write_pdb(1,'first structure',0d0) call secondary(.true.) call secondary2(.true.) do j=1,nbfrag if ( (bfrag(3,j).lt.bfrag(4,j) .or. & bfrag(4,j)-bfrag(2,j).gt.4) .and. & bfrag(2,j)-bfrag(1,j).gt.3 ) then nn=nn+1 if (bfrag(3,j).lt.bfrag(4,j)) then write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,& ",",bfrag(3,j)-1,"-",bfrag(4,j)-1 else write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,& ",",bfrag(4,j)-1,"-",bfrag(3,j)-1 endif endif enddo do i=1,nres theta(i)=theta2(i) phi(i)=phi2(i) alph(i)=alph2(i) omeg(i)=omeg2(i) enddo call chainbuild call geom_to_var(nvar,varia2) call write_pdb(2,'second structure',0d0) !------------------------------------------------------- !el#ifdef MPI ifun=-1 call contact_cp(varia,varia2,iff,ifun,7) if (minim) then !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,& '+ DIST eval',ifun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,& nfun/(time1-time0),' eval/s' call var_to_geom(nvar,varia) call chainbuild call write_pdb(999,'full min',etot) endif !el#endif call etotal(energy) etot=energy(0) call enerprint(energy) call intout call briefout(0,etot) return 10 write (iout,'(a)') 'Error reading test structure.' return end subroutine test__ !----------------------------------------------------------------------------- subroutine secondary(lprint) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' integer :: ncont,icont(2,nres*nres/2),isec(nres,3) logical :: lprint,not_done real(kind=4) :: dcont(nres*nres/2),d real(kind=4) :: rcomp = 7.0 real(kind=4) :: rbeta = 5.2 real(kind=4) :: ralfa = 5.2 real(kind=4) :: r310 = 6.6 real(kind=8),dimension(3) :: xpi,xpj integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,& nhelix call chainbuild !d call write_pdb(99,'sec structure',0d0) ncont=0 nbfrag=0 nhfrag=0 do i=1,nres isec(i,1)=0 isec(i,2)=0 isec(i,3)=0 enddo do i=2,nres-3 do k=1,3 xpi(k)=0.5d0*(c(k,i-1)+c(k,i)) enddo do j=i+2,nres do k=1,3 xpj(k)=0.5d0*(c(k,j-1)+c(k,j)) enddo !d d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) + !d & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) + !d & (c(3,i)-c(3,j))*(c(3,i)-c(3,j)) !d print *,'CA',i,j,d d = (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + & (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + & (xpi(3)-xpj(3))*(xpi(3)-xpj(3)) if ( d.lt.rcomp*rcomp) then ncont=ncont+1 icont(1,ncont)=i icont(2,ncont)=j dcont(ncont)=sqrt(d) endif enddo enddo if (lprint) then write (iout,*) write (iout,'(a)') '#PP contact map distances:' do i=1,ncont write (iout,'(3i4,f10.5)') & i,icont(1,i),icont(2,i),dcont(i) enddo endif ! finding parallel beta !d write (iout,*) '------- looking for parallel beta -----------' nbeta=0 nstrand=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. & isec(i1,1).le.1.and.isec(j1,1).le.1.and. & (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. & (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. & (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. & (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) & ) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1,dcont(i) not_done=.true. do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) & .and. dcont(j).le.rbeta .and. & isec(i1,1).le.1.and.isec(j1,1).le.1.and. & (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. & (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. & (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. & (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) & ) goto 5 enddo not_done=.false. 5 continue !d write (iout,*) i1,j1,dcont(j),not_done enddo j1=j1-1 i1=i1-1 if (i1-ii1.gt.1) then ii1=max0(ii1-1,1) jj1=max0(jj1-1,1) nbeta=nbeta+1 if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1 nbfrag=nbfrag+1 bfrag(1,nbfrag)=ii1 bfrag(2,nbfrag)=i1 bfrag(3,nbfrag)=jj1 bfrag(4,nbfrag)=j1 do ij=ii1,i1 isec(ij,1)=isec(ij,1)+1 isec(ij,1+isec(ij,1))=nbeta enddo do ij=jj1,j1 isec(ij,1)=isec(ij,1)+1 isec(ij,1+isec(ij,1))=nbeta enddo if(lprint) then nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-1,"..",i1-1,"'" endif nstrand=nstrand+1 if (nbeta.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",jj1-1,"..",j1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-1,i1-1,jj1-1,j1-1 endif endif endif enddo ! finding antiparallel beta !d write (iout,*) '--------- looking for antiparallel beta ---------' do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (dcont(i).le.rbeta.and. & isec(i1,1).le.1.and.isec(j1,1).le.1.and. & (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. & (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. & (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. & (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) & ) then ii1=i1 jj1=j1 !d write (iout,*) i1,j1,dcont(i) not_done=.true. do while (not_done) i1=i1+1 j1=j1-1 do j=1,ncont if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. & isec(i1,1).le.1.and.isec(j1,1).le.1.and. & (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. & (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. & (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. & (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) & .and. dcont(j).le.rbeta ) goto 6 enddo not_done=.false. 6 continue !d write (iout,*) i1,j1,dcont(j),not_done enddo i1=i1-1 j1=j1+1 if (i1-ii1.gt.1) then if(lprint)write (iout,*)'antiparallel beta',& nbeta,ii1-1,i1,jj1,j1-1 nbfrag=nbfrag+1 bfrag(1,nbfrag)=max0(ii1-1,1) bfrag(2,nbfrag)=i1 bfrag(3,nbfrag)=jj1 bfrag(4,nbfrag)=max0(j1-1,1) nbeta=nbeta+1 iii1=max0(ii1-1,1) do ij=iii1,i1 isec(ij,1)=isec(ij,1)+1 isec(ij,1+isec(ij,1))=nbeta enddo jjj1=max0(j1-1,1) do ij=jjj1,jj1 isec(ij,1)=isec(ij,1)+1 isec(ij,1+isec(ij,1))=nbeta enddo if (lprint) then nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",ii1-2,"..",i1-1,"'" endif nstrand=nstrand+1 if (nstrand.le.9) then write(12,'(a18,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" else write(12,'(a18,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'strand",nstrand,& "' 'num = ",j1-2,"..",jj1-1,"'" endif write(12,'(a8,4i4)') & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2 endif endif endif enddo if (nstrand.gt.0.and.lprint) then write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1" do i=2,nstrand if (i.le.9) then write(12,'(a9,i1,$)') " | strand",i else write(12,'(a9,i2,$)') " | strand",i endif enddo write(12,'(a1)') "'" endif ! finding alpha or 310 helix nhelix=0 do i=1,ncont i1=icont(1,i) j1=icont(2,i) if (j1.eq.i1+3.and.dcont(i).le.r310 & .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then !d if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i) !d if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i) ii1=i1 jj1=j1 if (isec(ii1,1).eq.0) then not_done=.true. else not_done=.false. endif do while (not_done) i1=i1+1 j1=j1+1 do j=1,ncont if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10 enddo not_done=.false. 10 continue !d write (iout,*) i1,j1,not_done enddo j1=j1-1 if (j1-ii1.gt.4) then nhelix=nhelix+1 !d write (iout,*)'helix',nhelix,ii1,j1 nhfrag=nhfrag+1 hfrag(1,nhfrag)=ii1 hfrag(2,nhfrag)=max0(j1-1,1) do ij=ii1,j1 isec(ij,1)=-1 enddo if (lprint) then write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2 if (nhelix.le.9) then write(12,'(a17,i1,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" else write(12,'(a17,i2,a9,i3,a2,i3,a1)') & "DefPropRes 'helix",nhelix,& "' 'num = ",ii1-1,"..",j1-2,"'" endif endif endif endif enddo if (nhelix.gt.0.and.lprint) then write(12,'(a26,$)') "DefPropRes 'helix' 'helix1" do i=2,nhelix if (nhelix.le.9) then write(12,'(a8,i1,$)') " | helix",i else write(12,'(a8,i2,$)') " | helix",i endif enddo write(12,'(a1)') "'" endif if (lprint) then write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'" write(12,'(a20)') "XMacStand ribbon.mac" endif return end subroutine secondary !----------------------------------------------------------------------------- subroutine contact_cp2(var,var2,iff,ieval,in_pdb) use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.MINIM' character(len=50) :: linia integer :: nf,ij(4) real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres) real(kind=8) :: time0,time1 integer :: iff(nres),ieval real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres) !el local variables integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode real(kind=8) :: wstrain0,etot integer :: maxres22 maxres22=nres*(nres+1)/2 if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES) call var_to_geom(nvar,var) call chainbuild nhpb0=nhpb ind=0 do i=1,nres-3 do j=i+3,nres ind=ind+1 if ( iff(i).eq.1.and.iff(j).eq.1 ) then d0(ind)=DIST(i,j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else w(ind)=0.0 endif enddo enddo call hpb_partition do i=1,nres theta1(i)=theta(i) phi1(i)=phi(i) alph1(i)=alph(i) omeg1(i)=omeg(i) enddo call var_to_geom(nvar,var2) do i=1,nres if ( iff(i).eq.1 ) then theta(i)=theta1(i) phi(i)=phi1(i) alph(i)=alph1(i) omeg(i)=omeg1(i) endif enddo call chainbuild !d call write_pdb(3,'combined structure',0d0) !d time0=MPI_WTIME() NX=NRES-3 NY=((NRES-4)*(NRES-5))/2 call distfit(.true.,200) !d time1=MPI_WTIME() !d write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec' ipot0=ipot maxmin0=maxmin maxfun0=maxfun wstrain0=wstrain ipot=6 maxmin=2000 maxfun=5000 call geom_to_var(nvar,var) !d time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun !d time1=MPI_WTIME() !d write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0, !d & nfun/(time1-time0),' SOFT eval/s' call var_to_geom(nvar,var) call chainbuild iwsk=0 nf=0 if (iff(1).eq.1) then iwsk=1 nf=nf+1 ij(nf)=0 endif do i=2,nres if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then iwsk=1 nf=nf+1 ij(nf)=i endif if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then iwsk=0 nf=nf+1 ij(nf)=i-1 endif enddo if (iff(nres).eq.1) then nf=nf+1 ij(nf)=nres endif !d write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') !d & "select",ij(1),"-",ij(2), !d & ",",ij(3),"-",ij(4) !d call write_pdb(in_pdb,linia,etot) ipot=ipot0 maxmin=maxmin0 maxfun=maxfun0 !d time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) !d write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun ieval=nfun !d time1=MPI_WTIME() !d write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0, !d & nfun/(time1-time0),' eval/s' !d call var_to_geom(nvar,var) !d call chainbuild !d call write_pdb(6,'dist structure',etot) nhpb= nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 return end subroutine contact_cp2 !----------------------------------------------------------------------------- subroutine contact_cp(var,var2,iff,ieval,in_pdb) use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.MINIM' character(len=50) :: linia integer :: nf,ij(4) real(kind=8) :: energy(0:n_ene) real(kind=8),dimension(6*nres) :: var,var2 !(maxvar) (maxvar=6*maxres) real(kind=8) :: time0,time1 integer :: iff(nres),ieval real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1 !(maxres) logical :: debug !el local variables integer :: in_pdb,i,j,ind,iwsk debug=.false. ! debug=.true. if (ieval.eq.-1) debug=.true. ! ! store selected dist. constrains from 1st structure ! #ifdef OSF ! Intercept NaNs in the coordinates ! write(iout,*) (var(i),i=1,nvar) x_sum=0.D0 do i=1,nvar x_sum=x_sum+var(i) enddo if (x_sum.ne.x_sum) then write(iout,*)" *** contact_cp : Found NaN in coordinates" call flush(iout) print *," *** contact_cp : Found NaN in coordinates" return endif #endif call var_to_geom(nvar,var) call chainbuild nhpb0=nhpb ind=0 do i=1,nres-3 do j=i+3,nres ind=ind+1 if ( iff(i).eq.1.and.iff(j).eq.1 ) then d0(ind)=DIST(i,j) w(ind)=10.0 nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=10.0 dhpb(nhpb)=d0(ind) else w(ind)=0.0 endif enddo enddo call hpb_partition do i=1,nres theta1(i)=theta(i) phi1(i)=phi(i) alph1(i)=alph(i) omeg1(i)=omeg(i) enddo ! ! freeze sec.elements from 2nd structure ! do i=1,nres mask_phi(i)=1 mask_theta(i)=1 mask_side(i)=1 enddo call var_to_geom(nvar,var2) call secondary2(debug) do j=1,nbfrag do i=bfrag(1,j),bfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo if (bfrag(3,j).le.bfrag(4,j)) then do i=bfrag(3,j),bfrag(4,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo else do i=bfrag(4,j),bfrag(3,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo endif enddo do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) mask(i)=0 mask_phi(i)=0 mask_theta(i)=0 enddo enddo mask_r=.true. ! ! copy selected res from 1st to 2nd structure ! do i=1,nres if ( iff(i).eq.1 ) then theta(i)=theta1(i) phi(i)=phi1(i) alph(i)=alph1(i) omeg(i)=omeg1(i) endif enddo if(debug) then ! ! prepare description in linia variable ! iwsk=0 nf=0 if (iff(1).eq.1) then iwsk=1 nf=nf+1 ij(nf)=1 endif do i=2,nres if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then iwsk=1 nf=nf+1 ij(nf)=i endif if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then iwsk=0 nf=nf+1 ij(nf)=i-1 endif enddo if (iff(nres).eq.1) then nf=nf+1 ij(nf)=nres endif write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') & "SELECT",ij(1)-1,"-",ij(2)-1,& ",",ij(3)-1,"-",ij(4)-1 endif ! ! run optimization ! call contact_cp_min(var,ieval,in_pdb,linia,debug) return end subroutine contact_cp !----------------------------------------------------------------------------- subroutine contact_cp_min(var,ieval,in_pdb,linia,debug) !el use minim ! ! input : theta,phi,alph,omeg,in_pdb,linia,debug ! output : var,ieval ! ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.MINIM' character(len=50) :: linia integer :: nf,ij(4) real(kind=8) :: energy(0:n_ene) real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) real(kind=8) :: time0,time1 integer :: ieval,info(3) logical :: debug,fail,reduce,change !check_var, !el local variables integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,& iretcode,nfun real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,& wtor_d01,wstrain0,etot write(iout,'(a20,i6,a20)') & '------------------',in_pdb,'-------------------' !el#ifdef MPI if (debug) then call chainbuild call write_pdb(1000+in_pdb,'combined structure',0d0) !el#ifdef MPI time0=MPI_WTIME() !el#endif endif !el#endif ! ! run optimization of distances ! ! uses d0(),w() and mask() for frozen 2D ! !test--------------------------------------------- !test NX=NRES-3 !test NY=((NRES-4)*(NRES-5))/2 !test call distfit(debug,5000) do i=1,nres mask_side(i)=0 enddo ipot01=ipot maxmin01=maxmin maxfun01=maxfun ! wstrain01=wstrain wsc01=wsc wscp01=wscp welec01=welec wvdwpp01=wvdwpp ! wang01=wang wscloc01=wscloc wtor01=wtor wtor_d01=wtor_d ipot=6 maxmin=2000 maxfun=4000 ! wstrain=1.0 wsc=0.0 wscp=0.0 welec=0.0 wvdwpp=0.0 ! wang=0.0 wscloc=0.0 wtor=0.0 wtor_d=0.0 call geom_to_var(nvar,var) !de change=reduce(var) if (check_var(var,info)) then write(iout,*) 'cp_min error in input' print *,'cp_min error in input' return endif !d call etotal(energy(0)) !d call enerprint(energy(0)) !d call check_eint !el#ifdef MPI time0=MPI_WTIME() !dtest call minimize(etot,var,iretcode,nfun) !dtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun time1=MPI_WTIME() !el#endif !d call etotal(energy(0)) !d call enerprint(energy(0)) !d call check_eint do i=1,nres mask_side(i)=1 enddo ipot=ipot01 maxmin=maxmin01 maxfun=maxfun01 ! wstrain=wstrain01 wsc=wsc01 wscp=wscp01 welec=welec01 wvdwpp=wvdwpp01 ! wang=wang01 wscloc=wscloc01 wtor=wtor01 wtor_d=wtor_d01 !test-------------------------------------------------- if(debug) then !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec' call write_pdb(2000+in_pdb,'distfit structure',0d0) endif ipot0=ipot maxmin0=maxmin maxfun0=maxfun wstrain0=wstrain ! ! run soft pot. optimization ! with constrains: ! nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition ! and frozen 2D: ! mask_phi(),mask_theta(),mask_side(),mask_r ! ipot=6 maxmin=2000 maxfun=4000 !el#ifdef MPI !de change=reduce(var) !de if (check_var(var,info)) write(iout,*) 'error before soft' !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,& nfun/(time1-time0),' SOFT eval/s' if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(3000+in_pdb,'soft structure',etot) endif !el#endif ! ! run full UNRES optimization with constrains and frozen 2D ! the same variables as soft pot. optimizatio ! ipot=ipot0 maxmin=maxmin0 maxfun=maxfun0 ! ! check overlaps before calling full UNRES minim ! call var_to_geom(nvar,var) call chainbuild call etotal(energy) #ifdef OSF write(iout,*) 'N7 ',energy(0) if (energy(0).ne.energy(0)) then write(iout,*) 'N7 error - gives NaN',energy(0) endif #endif ieval=1 if (energy(1).eq.1.0d20) then write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1) call overlap_sc(fail) if(.not.fail) then call etotal(energy) ieval=ieval+1 write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1) else mask_r=.false. nhpb= nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 return endif endif call flush(iout) ! !dte time0=MPI_WTIME() !de change=reduce(var) !de if (check_var(var,info)) then !de write(iout,*) 'error before mask dist' !de call var_to_geom(nvar,var) !de call chainbuild !de call write_pdb(10000+in_pdb,'before mask dist',etot) !de endif !dte call minimize(etot,var,iretcode,nfun) !dte write(iout,*)'SUMSL MASK DIST return code is',iretcode, !dte & ' eval ',nfun !dte ieval=ieval+nfun !dte !dte time1=MPI_WTIME() !dte write (iout,'(a,f6.2,f8.2,a)') !dte & ' Time for mask dist min.',time1-time0, !dte & nfun/(time1-time0),' eval/s' !dte call flush(iout) !el#ifdef MPI if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(4000+in_pdb,'mask dist',etot) endif ! ! switch off freezing of 2D and ! run full UNRES optimization with constrains ! mask_r=.false. !el#ifdef MPI time0=MPI_WTIME() !el#endif !de change=reduce(var) !de if (check_var(var,info)) then !de write(iout,*) 'error before dist' !de call var_to_geom(nvar,var) !de call chainbuild !de call write_pdb(11000+in_pdb,'before dist',etot) !de endif call minimize(etot,var,iretcode,nfun) !de change=reduce(var) !de if (check_var(var,info)) then !de write(iout,*) 'error after dist',ico !de call var_to_geom(nvar,var) !de call chainbuild !de call write_pdb(12000+in_pdb+ico*1000,'after dist',etot) !de endif write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun ieval=ieval+nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,& nfun/(time1-time0),' eval/s' !de call etotal(energy(0)) !de write(iout,*) 'N7 after dist',energy(0) call flush(iout) if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(in_pdb,linia,etot) endif !el#endif ! ! reset constrains ! nhpb= nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 return end subroutine contact_cp_min !----------------------------------------------------------------------------- subroutine softreg use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! include 'COMMON.INTERACT' ! ! include 'COMMON.DISTFIT' integer :: iff(nres) real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) ! logical :: debug,ltest,fail character(len=50) :: linia integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,& iretcode,nfun real(kind=8) :: wstrain0,wang0,etot ! linia='test' debug=.true. in_pdb=0 !------------------------ ! ! freeze sec.elements ! do i=1,nres mask_phi(i)=1 mask_theta(i)=1 mask_side(i)=1 iff(i)=0 enddo do j=1,nbfrag do i=bfrag(1,j),bfrag(2,j) mask_phi(i)=0 mask_theta(i)=0 iff(i)=1 enddo if (bfrag(3,j).le.bfrag(4,j)) then do i=bfrag(3,j),bfrag(4,j) mask_phi(i)=0 mask_theta(i)=0 iff(i)=1 enddo else do i=bfrag(4,j),bfrag(3,j) mask_phi(i)=0 mask_theta(i)=0 iff(i)=1 enddo endif enddo do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) mask_phi(i)=0 mask_theta(i)=0 iff(i)=1 enddo enddo mask_r=.true. nhpb0=nhpb ! ! store dist. constrains ! do i=1,nres-3 do j=i+3,nres if ( iff(i).eq.1.and.iff(j).eq.1 ) then nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=0.1 dhpb(nhpb)=DIST(i,j) endif enddo enddo call hpb_partition if (debug) then call chainbuild call write_pdb(100+in_pdb,'input reg. structure',0d0) endif ipot0=ipot maxmin0=maxmin maxfun0=maxfun wstrain0=wstrain wang0=wang ! ! run soft pot. optimization ! ipot=6 wang=3.0 maxmin=2000 maxfun=4000 call geom_to_var(nvar,var) !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,& nfun/(time1-time0),' SOFT eval/s' if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(300+in_pdb,'soft structure',etot) endif ! ! run full UNRES optimization with constrains and frozen 2D ! the same variables as soft pot. optimizatio ! ipot=ipot0 wang=wang0 maxmin=maxmin0 maxfun=maxfun0 !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL MASK DIST return code is',iretcode,& ' eval ',nfun ieval=nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)') & ' Time for mask dist min.',time1-time0,& nfun/(time1-time0),' eval/s' if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(400+in_pdb,'mask & dist',etot) endif ! ! switch off constrains and ! run full UNRES optimization with frozen 2D ! ! ! reset constrains ! nhpb_c=nhpb nhpb=nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun ieval=ieval+nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for mask min.',time1-time0,& nfun/(time1-time0),' eval/s' if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(500+in_pdb,'mask 2d frozen',etot) endif mask_r=.false. ! ! run full UNRES optimization with constrains and NO frozen 2D ! nhpb=nhpb_c link_start=1 link_end=nhpb maxfun=maxfun0/5 do ico=1,5 wstrain=wstrain0/ico !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,'(a10,f6.3,a14,i3,a6,i5)') & ' SUMSL DIST',wstrain,' return code is',iretcode,& ' eval ',nfun ieval=nfun !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)') & ' Time for dist min.',time1-time0,& nfun/(time1-time0),' eval/s' if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(600+in_pdb+ico,'dist cons',etot) endif enddo ! nhpb=nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 maxfun=maxfun0 ! if (minim) then !el#ifdef MPI time0=MPI_WTIME() !el#endif call minimize(etot,var,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,& '+ DIST eval',ieval !el#ifdef MPI time1=MPI_WTIME() !el#endif write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0,& nfun/(time1-time0),' eval/s' call var_to_geom(nvar,var) call chainbuild call write_pdb(999,'full min',etot) endif !el#endif return end subroutine softreg !----------------------------------------------------------------------------- subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij) use geometry, only:dist !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! include 'COMMON.CHAIN' real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) integer :: jdata(5),isec(nres) ! !el local variables integer :: i1,i2,i3,i4,i5,ieval,ij integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico real(kind=8) :: etot,wscloc0,wstrain0 jdata(1)=i1 jdata(2)=i2 jdata(3)=i3 jdata(4)=i4 jdata(5)=i5 call secondary2(.false.) do i=1,nres isec(i)=0 enddo do j=1,nbfrag do i=bfrag(1,j),bfrag(2,j) isec(i)=1 enddo do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j)) isec(i)=1 enddo enddo do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) isec(i)=2 enddo enddo ! ! cut strands at the ends ! if (jdata(2)-jdata(1).gt.3) then jdata(1)=jdata(1)+1 jdata(2)=jdata(2)-1 if (jdata(3).lt.jdata(4)) then jdata(3)=jdata(3)+1 jdata(4)=jdata(4)-1 else jdata(3)=jdata(3)-1 jdata(4)=jdata(4)+1 endif endif !v call chainbuild !v call etotal(energy(0)) !v etot=energy(0) !v write(iout,*) nnt,nct,etot !v call write_pdb(ij*100,'first structure',etot) !v write(iout,*) 'N16 test',(jdata(i),i=1,5) !------------------------ ! generate constrains ! ishift=jdata(5)-2 if(ishift.eq.0) ishift=-2 nhpb0=nhpb call chainbuild do i=jdata(1),jdata(2) isec(i)=-1 if(jdata(4).gt.jdata(3))then do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2 isec(j)=-1 !d print *,i,j,j+ishift nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=1000.0 dhpb(nhpb)=DIST(i,j+ishift) enddo else do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1 isec(j)=-1 !d print *,i,j,j+ishift nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=1000.0 dhpb(nhpb)=DIST(i,j+ishift) enddo endif enddo do i=nnt,nct-2 do j=i+2,nct if(isec(i).gt.0.or.isec(j).gt.0) then !d print *,i,j nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=0.1 dhpb(nhpb)=DIST(i,j) endif enddo enddo call hpb_partition call geom_to_var(nvar,var) maxfun0=maxfun wstrain0=wstrain maxfun=4000/5 do ico=1,5 wstrain=wstrain0/ico !v time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) write(iout,'(a10,f6.3,a14,i3,a6,i5)') & ' SUMSL DIST',wstrain,' return code is',iretcode,& ' eval ',nfun ieval=ieval+nfun !v time1=MPI_WTIME() !v write (iout,'(a,f6.2,f8.2,a)') !v & ' Time for dist min.',time1-time0, !v & nfun/(time1-time0),' eval/s' !v call var_to_geom(nvar,var) !v call chainbuild !v call write_pdb(ij*100+ico,'dist cons',etot) enddo ! nhpb=nhpb0 call hpb_partition wstrain=wstrain0 maxfun=maxfun0 ! !d print *,etot wscloc0=wscloc wscloc=10.0 call sc_move(nnt,nct,100,100d0,nft_sc,etot) wscloc=wscloc0 !v call chainbuild !v call etotal(energy(0)) !v etot=energy(0) !v call write_pdb(ij*100+10,'sc_move',etot) !d call intout !d print *,nft_sc,etot return end subroutine beta_slide !----------------------------------------------------------------------------- subroutine beta_zip(i1,i2,ieval,ij) !el use minim ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' ! include 'COMMON.SBRIDGE' ! include 'COMMON.CONTROL' ! include 'COMMON.FFIELD' ! include 'COMMON.MINIM' ! include 'COMMON.CHAIN' real(kind=8) :: time0,time1 real(kind=8) :: energy(0:n_ene),ee real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) character(len=10) :: test !el local variables integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0 real(kind=8) :: etot,wstrain0 !v call chainbuild !v call etotal(energy(0)) !v etot=energy(0) !v write(test,'(2i5)') i1,i2 !v call write_pdb(ij*100,test,etot) !v write(iout,*) 'N17 test',i1,i2,etot,ij ! ! generate constrains ! nhpb0=nhpb nhpb=nhpb+1 ihpb(nhpb)=i1 jhpb(nhpb)=i2 forcon(nhpb)=1000.0 dhpb(nhpb)=4.0 call hpb_partition call geom_to_var(nvar,var) maxfun0=maxfun wstrain0=wstrain maxfun=1000/5 do ico=1,5 wstrain=wstrain0/ico !v time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) write(iout,'(a10,f6.3,a14,i3,a6,i5)') & ' SUMSL DIST',wstrain,' return code is',iretcode,& ' eval ',nfun ieval=ieval+nfun !v time1=MPI_WTIME() !v write (iout,'(a,f6.2,f8.2,a)') !v & ' Time for dist min.',time1-time0, !v & nfun/(time1-time0),' eval/s' ! do not comment the next line call var_to_geom(nvar,var) !v call chainbuild !v call write_pdb(ij*100+ico,'dist cons',etot) enddo nhpb=nhpb0 call hpb_partition wstrain=wstrain0 maxfun=maxfun0 !v call etotal(energy(0)) !v etot=energy(0) !v write(iout,*) 'N17 test end',i1,i2,etot,ij return end subroutine beta_zip !----------------------------------------------------------------------------- ! thread.F !----------------------------------------------------------------------------- subroutine thread_seq use geometry, only:dist use random, only:iran_num use control, only:tcpu use regularize_, only:regularize use mcm_data, only: nsave_part,nacc_tot ! Thread the sequence through a database of known structures ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use MPI_data !include 'COMMON.INFO' use MPI_ #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.CHAIN' ! include 'COMMON.DBASE' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.THREAD' ! include 'COMMON.FFIELD' ! include 'COMMON.SBRIDGE' ! include 'COMMON.HEADER' ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' ! include 'COMMON.CONTACTS' ! include 'COMMON.MCM' ! include 'COMMON.NAMES' #ifdef MPI integer :: ThreadId,ThreadType,Kwita #endif real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) real(kind=8) :: przes(3),obr(3,3) real(kind=8) :: time_for_thread logical :: found_pattern,non_conv character(len=32) :: head_pdb real(kind=8) :: energia(0:n_ene) integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,& link_end0,iproc real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1 n_ene_comp=nprint_ene ! ! Body ! #ifdef MPI if (me.eq.king) then do i=1,nctasks nsave_part(i)=0 enddo endif nacc_tot=0 Kwita=0 #endif close(igeom) close(ipdb) close(istat) do i=1,maxthread do j=1,14 ener0(j,i)=0.0D0 ener(j,i)=0.0D0 enddo enddo nres0=nct-nnt+1 ave_time_for_thread=0.0D0 max_time_for_thread=0.0D0 !d print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0 nthread=nexcl+nthread do ithread=1,nthread found_pattern=.false. itrial=0 do while (.not.found_pattern) itrial=itrial+1 if (itrial.gt.1000) then write (iout,'(/a/)') 'Too many attempts to find pattern.' nthread=ithread-1 #ifdef MPI call recv_stop_sig(Kwita) call send_stop_sig(-3) #endif goto 777 endif ! Find long enough chain in the database ii=iran_num(1,nseq) nres_t=nres_base(1,ii) ! Select the starting position to thread. print *,'nseq',nseq,' ii=',ii,' nres_t=',& nres_t,' nres0=',nres0 if (nres_t.ge.nres0) then ist=iran_num(0,nres_t-nres0) #ifdef MPI if (Kwita.eq.0) call recv_stop_sig(Kwita) if (Kwita.lt.0) then write (iout,*) 'Stop signal received. Terminating.' write (*,*) 'Stop signal received. Terminating.' nthread=ithread-1 write (*,*) 'ithread=',ithread,' nthread=',nthread goto 777 endif call pattern_receive #endif do i=1,nexcl if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10 enddo found_pattern=.true. endif ! If this point is reached, the pattern has not yet been examined. 10 continue ! print *,'found_pattern:',found_pattern enddo nexcl=nexcl+1 iexam(1,nexcl)=ii iexam(2,nexcl)=ist #ifdef MPI if (Kwita.eq.0) call recv_stop_sig(Kwita) if (Kwita.lt.0) then write (iout,*) 'Stop signal received. Terminating.' nthread=ithread-1 write (*,*) 'ithread=',ithread,' nthread=',nthread goto 777 endif call pattern_send #endif ipatt(1,ithread)=ii ipatt(2,ithread)=ist #ifdef MPI write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') & 'Processor:',me,' Attempt:',ithread,& ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),& ' start at res.',ist+1 write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,& ' Attempt:',ithread,& ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),& ' start at res.',ist+1 #else write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') & 'Attempt:',ithread,& ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),& ' start at res.',ist+1 write (*,'(a,i5,2a,i3,a,i3,a,i3)') & 'Attempt:',ithread,& ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),& ' start at res.',ist+1 #endif ipattern=ii ! Copy coordinates from the database. ist=ist-(nnt-1) do i=nnt,nct do j=1,3 c(j,i)=cart_base(j,i+ist,ii) ! cref(j,i)=c(j,i) enddo !d write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3) enddo !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr, !d non_conv) !d write (iout,'(a,f10.5)') !d & 'Initial RMS deviation from reference structure:',rms if (itype(nres).eq.ntyp1) then do j=1,3 dcj=c(j,nres-2)-c(j,nres-3) c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo endif if (itype(1).eq.ntyp1) then do j=1,3 dcj=c(j,4)-c(j,3) c(j,1)=c(j,2)-dcj c(j,nres+1)=c(j,1) enddo endif call int_from_cart(.false.,.false.) !d print *,'Exit INT_FROM_CART.' !d print *,'nhpb=',nhpb do i=nss+1,nhpb ii=ihpb(i) jj=jhpb(i) dhpb(i)=dist(ii,jj) ! write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i) enddo ! stop 'End generate' ! Generate SC conformations. call sc_conf ! call intout #ifdef MPI !d print *,'Processor:',me,': exit GEN_SIDE.' #else !d print *,'Exit GEN_SIDE.' #endif ! Calculate initial energy. call chainbuild call etotal(energia) etot=energia(0) do i=1,n_ene_comp ener0(i,ithread)=energia(i) enddo ener0(n_ene_comp+1,ithread)=energia(0) if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.true.) ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,& icont,icont_ref) ener0(n_ene_comp+2,ithread)=rms ener0(n_ene_comp+4,ithread)=frac ener0(n_ene_comp+5,ithread)=frac_nn endif ener0(n_ene_comp+3,ithread)=0.0d0 ! Minimize energy. #ifdef MPI print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.' #else print*,'ithread=',ithread,' Start REGULARIZE.' #endif curr_tim=tcpu() call regularize(nct-nnt+1,etot,rms,& cart_base(1,ist+nnt,ipattern),iretcode) curr_tim1=tcpu() time_for_thread=curr_tim1-curr_tim ave_time_for_thread= & ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread if (time_for_thread.gt.max_time_for_thread) & max_time_for_thread=time_for_thread #ifdef MPI print *,'Processor',me,': Exit REGULARIZE.' if (WhatsUp.eq.2) then write (iout,*) & 'Sufficient number of confs. collected. Terminating.' nthread=ithread-1 goto 777 else if (WhatsUp.eq.-1) then nthread=ithread-1 write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.' if (Kwita.eq.0) call recv_stop_sig(Kwita) call send_stop_sig(-2) goto 777 else if (WhatsUp.eq.-2) then nthread=ithread-1 write (iout,*) 'Timeup signal received. Terminating.' goto 777 else if (WhatsUp.eq.-3) then nthread=ithread-1 write (iout,*) 'Error stop signal received. Terminating.' goto 777 endif #else print *,'Exit REGULARIZE.' if (iretcode.eq.11) then write (iout,'(/a/)') & '******* Allocated time exceeded in SUMSL. The program will stop.' nthread=ithread-1 goto 777 endif #endif head_pdb=titel(:24)//':'//str_nam(ipattern) if (outpdb) call pdbout(etot,head_pdb,ipdb) if (outmol2) call mol2out(etot,head_pdb) ! call intout call briefout(ithread,etot) link_end0=link_end link_end=min0(link_end,nss) write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,& ' nss=',nss call etotal(energia) ! call enerprint(energia(0)) link_end=link_end0 !d call chainbuild !d call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv) !d write (iout,'(a,f10.5)') !d & 'RMS deviation from reference structure:',dsqrt(rms) do i=1,n_ene_comp ener(i,ithread)=energia(i) enddo ener(n_ene_comp+1,ithread)=energia(0) ener(n_ene_comp+3,ithread)=rms if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.true.) ener(n_ene_comp+2,ithread)=rms ener(n_ene_comp+4,ithread)=frac ener(n_ene_comp+5,ithread)=frac_nn endif call write_stat_thread(ithread,ipattern,ist) ! write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)') ! & ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11), ! & (ener(k,ithread),k=12,14) #ifdef MPI if (me.eq.king) then nacc_tot=nacc_tot+1 call pattern_receive call receive_MCM_info if (nacc_tot.ge.nthread) then write (iout,*) & 'Sufficient number of conformations collected nacc_tot=',& nacc_tot,'. Stopping other processors and terminating.' write (*,*) & 'Sufficient number of conformations collected nacc_tot=',& nacc_tot,'. Stopping other processors and terminating.' call recv_stop_sig(Kwita) if (Kwita.eq.0) call send_stop_sig(-1) nthread=ithread goto 777 endif else call send_MCM_info(2) endif #endif if (timlim-curr_tim1-safety .lt. max_time_for_thread) then write (iout,'(/2a)') & '********** There would be not enough time for another thread. ',& 'The program will stop.' write (*,'(/2a)') & '********** There would be not enough time for another thread. ',& 'The program will stop.' write (iout,'(a,1pe14.4/)') & 'Elapsed time for last threading step: ',time_for_thread nthread=ithread #ifdef MPI call recv_stop_sig(Kwita) call send_stop_sig(-2) #endif goto 777 else curr_tim=curr_tim1 write (iout,'(a,1pe14.4)') & 'Elapsed time for this threading step: ',time_for_thread endif #ifdef MPI if (Kwita.eq.0) call recv_stop_sig(Kwita) if (Kwita.lt.0) then write (iout,*) 'Stop signal received. Terminating.' write (*,*) 'Stop signal received. Terminating.' nthread=ithread write (*,*) 'nthread=',nthread,' ithread=',ithread goto 777 endif #endif enddo #ifdef MPI call send_stop_sig(-1) #endif 777 continue #ifdef MPI ! Any messages left for me? call pattern_receive if (Kwita.eq.0) call recv_stop_sig(Kwita) #endif call write_thread_summary #ifdef MPI if (king.eq.king) then Kwita=1 do while (Kwita.ne.0 .or. nacc_tot.ne.0) Kwita=0 nacc_tot=0 call recv_stop_sig(Kwita) call receive_MCM_info enddo do iproc=1,nprocs-1 call receive_thread_results(iproc) enddo call write_thread_summary else call send_thread_results endif #endif return end subroutine thread_seq !----------------------------------------------------------------------------- subroutine sc_conf ! Sample (hopefully) optimal SC orientations given backcone conformation. !el use comm_srutu use random, only:iran_num ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.DBASE' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.THREAD' ! include 'COMMON.FFIELD' ! include 'COMMON.SBRIDGE' ! include 'COMMON.HEADER' ! include 'COMMON.GEO' ! include 'COMMON.IOUNITS' real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) !el integer :: icall !el common /srutu/ icall real(kind=8) :: energia(0:n_ene) logical :: glycine,fail integer :: i,maxsample,link_end0,ind_sc,isample real(kind=8) :: alph0,omeg0,e1,e0 maxsample=10 link_end0=link_end link_end=min0(link_end,nss) do i=nnt,nct if (itype(i).ne.10) then !d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1) call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail) endif enddo call chainbuild call etotal(energia) e0 = energia(0) do isample=1,maxsample ! Choose a non-glycine side chain. glycine=.true. do while(glycine) ind_sc=iran_num(nnt,nct) glycine=(itype(ind_sc).eq.10) enddo alph0=alph(ind_sc) omeg0=omeg(ind_sc) call gen_side(itype(ind_sc),theta(ind_sc+1),alph(ind_sc),& omeg(ind_sc),fail) call chainbuild call etotal(energia) !d write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))') !d & 'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg, !d & ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1 e1=energia(0) if (e0.le.e1) then alph(ind_sc)=alph0 omeg(ind_sc)=omeg0 else e0=e1 endif enddo link_end=link_end0 return end subroutine sc_conf !----------------------------------------------------------------------------- ! minim_jlee.F !----------------------------------------------------------------------------- logical function check_var(var,info) use MPI_data use geometry_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.GEO' ! include 'COMMON.SETUP' real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) integer,dimension(3) :: info integer :: i,j ! AL ------- check_var=.false. do i=nphi+ntheta+1,nphi+ntheta+nside ! Check the side chain "valence" angles alpha if (var(i).lt.1.0d-7) then write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!' write (iout,*) 'Processor',me,'received bad variables!!!!' write (iout,*) 'Variables' write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar) write (iout,*) 'Continuing calculations at this point',& ' could destroy the results obtained so far... ABORTING!!!!!!' write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') & 'valence angle alpha',i-nphi-ntheta,var(i),& 'n it',info(1),info(2),'mv ',info(3) write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!' write (*,*) 'Processor',me,'received bad variables!!!!' write (*,*) 'Variables' write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar) write (*,*) 'Continuing calculations at this point',& ' could destroy the results obtained so far... ABORTING!!!!!!' write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') & 'valence angle alpha',i-nphi-ntheta,var(i),& 'n it',info(1),info(2),'mv ',info(3) check_var=.true. return endif enddo ! Check the backbone "valence" angles theta do i=nphi+1,nphi+ntheta if (var(i).lt.1.0d-7) then write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!' write (iout,*) 'Processor',me,'received bad variables!!!!' write (iout,*) 'Variables' write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar) write (iout,*) 'Continuing calculations at this point',& ' could destroy the results obtained so far... ABORTING!!!!!!' write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') & 'valence angle theta',i-nphi,var(i),& 'n it',info(1),info(2),'mv ',info(3) write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!' write (*,*) 'Processor',me,'received bad variables!!!!' write (*,*) 'Variables' write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar) write (*,*) 'Continuing calculations at this point',& ' could destroy the results obtained so far... ABORTING!!!!!!' write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') & 'valence angle theta',i-nphi,var(i),& 'n it',info(1),info(2),'mv ',info(3) check_var=.true. return endif enddo return end function check_var !----------------------------------------------------------------------------- ! distfit.f !----------------------------------------------------------------------------- subroutine distfit(debug,maxit) use geometry_data, only: phi use compare_data use md_calc ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.DISTFIT' integer :: i,maxit,MAXMAR,IT,IMAR real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold !(maxres) logical :: debug,sing real(kind=8) :: TOL,RL,F0,AIN,F1 !input------------------------------------ ! NX=NRES-3 ! NY=((NRES-4)*(NRES-5))/2 !input------------------------------------ !test MAXIT=20 TOL=0.5 MAXMAR=10 RL=100.0 CALL TRANSFER(NRES,phi,phiold) F0=RDIF() !d WRITE (IOUT,*) 'DISTFIT: F0=',F0 DO IT=1,MAXIT CALL RDERIV CALL HEVAL DO I=1,NX DIAGH(I)=H(I,I) ENDDO RL=RL*0.1 DO IMAR=1,MAXMAR DO I=1,NX H(I,I)=DIAGH(I)+RL ENDDO CALL TRANSFER(NX,XX,X) CALL BANACH(NX,NRES,H,X,sing) AIN=0.0 DO I=1,NX AIN=AIN+DABS(X(I)) ENDDO IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN if (debug) then WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED' WRITE (IOUT,*) 'IT=',it,'F=',F0 endif RETURN ENDIF DO I=4,NRES phi(I)=phiold(I)+mask(i)*X(I-3) ! print *,X(I-3) ENDDO F1=RDIF() !d WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1 IF (F1.LT.F0) THEN CALL TRANSFER(NRES,phi,phiold) F0=F1 GOTO 1 ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN if (debug) then WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT' WRITE (IOUT,*) 'IT=',it,'F=',F1 endif RETURN ENDIF RL=RL*10.0 ENDDO WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED' WRITE (IOUT,*) 'IT=',it,'F=',F0 CALL TRANSFER(NRES,phiold,phi) RETURN 1 continue !d write (iout,*) "it",it," imar",imar," f0",f0 enddo WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit return end subroutine distfit !----------------------------------------------------------------------------- real(kind=8) function RDIF() use compare_data use geometry, only: dist ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.DISTFIT' integer :: i,j,ind real(kind=8) :: suma,DIJ ! print *,'in rdif' suma=0.0 ind=0 call chainbuild do i=1,nres-3 do j=i+3,nres ind=ind+1 if (w(ind).ne.0.0) then DIJ=DIST(i,j) suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind)) DDD(ind)=DIJ ! print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma endif enddo enddo RDIF=suma return end function RDIF !----------------------------------------------------------------------------- subroutine RDERIV use compare_data use geometry_data use geometry, only:dist ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.DISTFIT' ! include 'COMMON.GEO' integer :: i,j,k,l,I1,I2,IND real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU DO I=1,NY DO J=1,NX DRDG(I,J)=0.0 ENDDO ENDDO DO I=1,NX I1=I+1 I2=I+2 CALL VEC(I1,I2,E12) DO J=1,I DO K=1,3 R13(K)=C(K,J)-C(K,I1) ENDDO DO K=I2+1,NRES DO L=1,3 R24(L)=C(L,K)-C(L,I2) ENDDO IND=((J-1)*(2*NRES-J-6))/2+K-3 PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2) PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3) PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1) DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K) ENDDO ENDDO ENDDO return end subroutine RDERIV !----------------------------------------------------------------------------- subroutine HEVAL use compare_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.DISTFIT' integer :: i,k,j real(kind=8) :: XI,HII,BKI,BKIWK,HIJ DO I=1,NX XI=0.0 HII=0.0 DO K=1,NY BKI=DRDG(K,I) BKIWK=w(K)*BKI XI=XI+BKIWK*(D0(K)-DDD(K)) HII=HII+BKI*BKIWK ENDDO H(I,I)=HII XX(I)=XI DO J=I+1,NX HIJ=0.0 DO K=1,NY HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K) ENDDO H(I,J)=HIJ H(J,I)=HIJ ENDDO ENDDO return end subroutine HEVAL !----------------------------------------------------------------------------- subroutine VEC(I,J,U) ! use geometry_data, only: C ! Find the unit vector from atom (I) to atom (J). Store in U. ! ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' integer :: I,J,K real(kind=8),DIMENSION(3) :: U real(kind=8) :: ANORM,UK ANORM=0.0 DO K=1,3 UK=C(K,J)-C(K,I) ANORM=ANORM+UK*UK U(K)=UK ENDDO ANORM=SQRT(ANORM) DO K=1,3 U(K)=U(K)/ANORM ENDDO return end subroutine VEC !----------------------------------------------------------------------------- subroutine TRANSFER(N,X1,X2) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' integer :: N,I real(kind=8),DIMENSION(N) :: X1,X2 DO 1 I=1,N 1 X2(I)=X1(I) return end subroutine TRANSFER !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- subroutine alloc_compare_arrays maxres22=nres*(nres+1)/2 ! common.dbase ! common /struct/ in io_common: read_threadbase ! allocate(cart_base !(3,maxres_base,maxseq) ! allocate(nres_base !(3,maxseq) ! allocate(str_nam !(maxseq) ! common.distfit ! COMMON /c_frag/ in io_conf: readpdb if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3) if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3) ! COMMON /WAGI/ allocate(w(maxres22),d0(maxres22)) !(maxres22) ! COMMON /POCHODNE/ !el allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES) allocate(DDD(maxres22)) !(maxres22) allocate(H(nres,nres)) !(MAXRES,MAXRES) allocate(XX(nres)) !(MAXRES) ! COMMON /frozen/ allocate(mask(nres)) !(maxres) ! common.thread ! common /thread/ allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread) ! common /thread1/ allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread) return end subroutine alloc_compare_arrays !----------------------------------------------------------------------------- #endif !----------------------------------------------------------------------------- end module compare