X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=inline;f=source%2Funres%2Fcompare.f90;fp=source%2Funres%2Fcompare.f90;h=0000000000000000000000000000000000000000;hb=df0d15d1ab81e01e177d3d39354e72364b294e1c;hp=b65e57cb7a73c5ab657bfdb6c02de28a79f4d6db;hpb=513f5fb9e31afc4a12d8ae2d113f7677ef057892;p=unres4.git diff --git a/source/unres/compare.f90 b/source/unres/compare.f90 deleted file mode 100644 index b65e57c..0000000 --- a/source/unres/compare.f90 +++ /dev/null @@ -1,4552 +0,0 @@ - module compare -!----------------------------------------------------------------------------- - use io_units - use names - use geometry_data - use energy_data - use control_data -#if .not. defined WHAM_RUN && .not. 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 .not. defined WHAM_RUN && .not. 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