subroutine test implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif 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' double precision time0,time1 double precision energy(0:n_ene),ee double precision var(maxvar),var1(maxvar) integer j1,j2 logical debug,accepted debug=.true. call geom_to_var(nvar,var1) call chainbuild call etotal(energy(0)) 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) c phi(jr)=pinorm(phi(jr)+d) call chainbuild call etotal(energy(0)) 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(0)) 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. c 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) c minimize c 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) c -------- else phi(jr)=phiold endif enddo c minimize c call sc_move(2,nres-1,1,10d0,nft_sc,etot) c call geom_to_var(nvar,var) c c call chainbuild c call write_pdb(998 ,'sc min',etot) c c call minimize(etot,var,iretcode,nfun) c write(iout,*)'------------------------------------------------' c write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun c c call var_to_geom(nvar,var) c call chainbuild c call write_pdb(999,'full min',etot) return end subroutine test_n16 implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif 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' double precision time0,time1 double precision energy(0:n_ene),ee double precision var(maxvar),var1(maxvar) integer jdata(5) logical debug debug=.true. c call geom_to_var(nvar,var1) call chainbuild call etotal(energy(0)) 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 #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, & '+ DIST eval',ieval #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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_local implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' double precision time0,time1 double precision energy(0:n_ene),ee double precision varia(maxvar) c call chainbuild c call geom_to_var(nvar,varia) call write_pdb(1,'first structure',0d0) call etotal(energy(0)) 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_sc implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' double precision time0,time1 double precision energy(0:n_ene),ee double precision varia(maxvar) c call chainbuild c call geom_to_var(nvar,varia) call write_pdb(1,'first structure',0d0) call etotal(energy(0)) 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 c-------------------------------------------------------- subroutine bgrow(bstrand,nbstrand,in,ind,new) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' integer bstrand(maxres/3,6) 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 c------------------------------------------ subroutine test11 implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.GEO' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.FFIELD' include 'COMMON.MINIM' c include 'COMMON.DISTFIT' integer if(20,maxres),nif,ifa(20) integer ibc(0:maxres,0:maxres),istrand(20) integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0 integer itmp(20,maxres) double precision time0,time1 double precision energy(0:n_ene),ee double precision varia(maxvar),vorg(maxvar) c logical debug,ltest,usedbfrag(maxres/3) character*50 linia c integer betasheet(maxres),ibetasheet(maxres),nbetasheet integer bstrand(maxres/3,6),nbstrand c------------------------ debug=.true. c------------------------ 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 c 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 c 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 c------------------------ 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 c read (inp,*) (ifa(i),i=1,4) c do i=1,nres c read (inp,*,err=20,end=20) (if(j,i),j=1,4) c enddo c 20 nif=i-1 stop c------------------------ isa=4 is=2*isa-1 iconf=0 cccccccccccccccccccccccccccccccccc DO ig=1,is**isa-1 cccccccccccccccccccccccccccccccccc 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 c ifa(1),ifa(2),ifa(3),ifa(4) c 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 c 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 c------------------------ c c freeze sec.elements c 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. c------------------------ c generate constrains c 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 dd(ind)=d0(ind) enddo enddo call hpb_partition cd-------------------------- 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) cd nhpb=0 cd goto 901 c c call contact_cp_min(varia,ifun,iconf,linia,debug) if (minim) then #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, & '+ DIST eval',ifun #ifdef MPI time1=MPI_WTIME() #else time0=tcpu() #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 call etotal(energy(0)) etot=energy(0) call enerprint(energy(0)) cd call intout cd call briefout(0,etot) cd call secondary2(.true.) 901 CONTINUE ctest return cccccccccccccccccccccccccccccccccccc ENDIF ENDDO cccccccccccccccccccccccccccccccccccc return 10 write (iout,'(a)') 'Error reading test structure.' return end c-------------------------------------------------------- subroutine test3 implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.GEO' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.FFIELD' include 'COMMON.MINIM' c include 'COMMON.DISTFIT' integer if(3,maxres),nif integer ibc(maxres,maxres),istrand(20) integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0 double precision time0,time1 double precision energy(0:n_ene),ee double precision varia(maxvar) c logical debug,ltest character*50 linia c 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) c------------------------ call secondary2(debug) c------------------------ do i=1,nres do j=1,nres ibc(i,j)=0 enddo enddo c c freeze sec.elements and store indexes for beta constrains c 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. c ---------------- 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 c------------------------ 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 cd-------------------------- 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 c call contact_cp_min(varia,ieval,in_pdb,linia,debug) if (minim) then #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, & '+ DIST eval',ieval #ifdef MPI time1=MPI_WTIME() #else time0=tcpu() #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 call etotal(energy(0)) etot=energy(0) call enerprint(energy(0)) call intout call briefout(0,etot) call secondary2(.true.) return 10 write (iout,'(a)') 'Error reading test structure.' return end subroutine test__ implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.GEO' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.FFIELD' include 'COMMON.MINIM' c include 'COMMON.DISTFIT' integer if(2,2),ind integer iff(maxres) double precision time0,time1 double precision energy(0:n_ene),ee double precision theta2(maxres),phi2(maxres),alph2(maxres), & omeg2(maxres), & theta1(maxres),phi1(maxres),alph1(maxres), & omeg1(maxres) double precision varia(maxvar),varia2(maxvar) c 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 c------------------------ 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) c------------------------------------------------------- ifun=-1 call contact_cp(varia,varia2,iff,ifun,7) if (minim) then #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,varia,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, & '+ DIST eval',ifun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 call etotal(energy(0)) etot=energy(0) call enerprint(energy(0)) call intout call briefout(0,etot) return 10 write (iout,'(a)') 'Error reading test structure.' return end c------------------------------------------------- c------------------------------------------------- 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,maxres*maxres/2),isec(maxres,3) logical lprint,not_done real dcont(maxres*maxres/2),d real rcomp /7.0/ real rbeta /5.2/ real ralfa /5.2/ real r310 /6.6/ double precision xpi(3),xpj(3) call chainbuild cd 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 cd d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) + cd & (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) + cd & (c(3,i)-c(3,j))*(c(3,i)-c(3,j)) cd 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 c finding parallel beta cd 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 cd 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 cd 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 c finding antiparallel beta cd 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 cd 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 cd 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 c 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 cd if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i) cd 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 cd write (iout,*) i1,j1,not_done enddo j1=j1-1 if (j1-ii1.gt.4) then nhelix=nhelix+1 cd 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 c---------------------------------------------------------------------------- subroutine write_pdb(npdb,titelloc,ee) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' character*50 titelloc1 character*(*) titelloc character*3 zahl character*5 liczba5 double precision ee integer npdb,ilen external ilen titelloc1=titelloc lenpre=ilen(prefix) if (npdb.lt.1000) then call numstr(npdb,zahl) open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb') else if (npdb.lt.10000) then write(liczba5,'(i1,i4)') 0,npdb else write(liczba5,'(i5)') npdb endif open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb') endif call pdbout(ee,titelloc1,ipdb) close(ipdb) return end c----------------------------------------------------------- subroutine contact_cp2(var,var2,iff,ieval,in_pdb) implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SBRIDGE' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.DISTFIT' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.MINIM' character*50 linia integer nf,ij(4) double precision var(maxvar),var2(maxvar) double precision time0,time1 integer iff(maxres),ieval double precision theta1(maxres),phi1(maxres),alph1(maxres), & omeg1(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 cd call write_pdb(3,'combined structure',0d0) cd time0=MPI_WTIME() NX=NRES-3 NY=((NRES-4)*(NRES-5))/2 call distfit(.true.,200) cd time1=MPI_WTIME() cd 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) cd time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun cd time1=MPI_WTIME() cd write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0, cd & 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 cd write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') cd & "select",ij(1),"-",ij(2), cd & ",",ij(3),"-",ij(4) cd call write_pdb(in_pdb,linia,etot) ipot=ipot0 maxmin=maxmin0 maxfun=maxfun0 cd time0=MPI_WTIME() call minimize(etot,var,iretcode,nfun) cd write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun ieval=nfun cd time1=MPI_WTIME() cd write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0, cd & nfun/(time1-time0),' eval/s' cd call var_to_geom(nvar,var) cd call chainbuild cd call write_pdb(6,'dist structure',etot) nhpb= nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 return end c----------------------------------------------------------- subroutine contact_cp(var,var2,iff,ieval,in_pdb) 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*50 linia integer nf,ij(4) double precision energy(0:n_ene) double precision var(maxvar),var2(maxvar) double precision time0,time1 integer iff(maxres),ieval double precision theta1(maxres),phi1(maxres),alph1(maxres), & omeg1(maxres) logical debug debug=.false. c debug=.true. if (ieval.eq.-1) debug=.true. c c store selected dist. constrains from 1st structure c #ifdef OSF c Intercept NaNs in the coordinates c 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 c c freeze sec.elements from 2nd structure c 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. c c copy selected res from 1st to 2nd structure c 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 c c prepare description in linia variable c 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 c c run optimization c call contact_cp_min(var,ieval,in_pdb,linia,debug) return end subroutine contact_cp_min(var,ieval,in_pdb,linia,debug) c c input : theta,phi,alph,omeg,in_pdb,linia,debug c output : var,ieval c implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SBRIDGE' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' include 'COMMON.DISTFIT' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.MINIM' character*50 linia integer nf,ij(4) double precision energy(0:n_ene) double precision var(maxvar) double precision time0,time1 integer ieval,info(3) logical debug,fail,check_var,reduce,change write(iout,'(a20,i6,a20)') & '------------------',in_pdb,'-------------------' if (debug) then call chainbuild call write_pdb(1000+in_pdb,'combined structure',0d0) #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif endif c c run optimization of distances c c uses d0(),w() and mask() for frozen 2D c ctest--------------------------------------------- ctest NX=NRES-3 ctest NY=((NRES-4)*(NRES-5))/2 ctest call distfit(debug,5000) do i=1,nres mask_side(i)=0 enddo ipot01=ipot maxmin01=maxmin maxfun01=maxfun c wstrain01=wstrain wsc01=wsc wscp01=wscp welec01=welec wvdwpp01=wvdwpp c wang01=wang wscloc01=wscloc wtor01=wtor wtor_d01=wtor_d ipot=6 maxmin=2000 maxfun=4000 c wstrain=1.0 wsc=0.0 wscp=0.0 welec=0.0 wvdwpp=0.0 c wang=0.0 wscloc=0.0 wtor=0.0 wtor_d=0.0 call geom_to_var(nvar,var) cde change=reduce(var) if (check_var(var,info)) then write(iout,*) 'cp_min error in input' print *,'cp_min error in input' return endif cd call etotal(energy(0)) cd call enerprint(energy(0)) cd call check_eint #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif cdtest call minimize(etot,var,iretcode,nfun) cdtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #endif cd call etotal(energy(0)) cd call enerprint(energy(0)) cd call check_eint do i=1,nres mask_side(i)=1 enddo ipot=ipot01 maxmin=maxmin01 maxfun=maxfun01 c wstrain=wstrain01 wsc=wsc01 wscp=wscp01 welec=welec01 wvdwpp=wvdwpp01 c wang=wang01 wscloc=wscloc01 wtor=wtor01 wtor_d=wtor_d01 ctest-------------------------------------------------- if(debug) then #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 c c run soft pot. optimization c with constrains: c nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition c and frozen 2D: c mask_phi(),mask_theta(),mask_side(),mask_r c ipot=6 maxmin=2000 maxfun=4000 cde change=reduce(var) cde if (check_var(var,info)) write(iout,*) 'error before soft' #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 c c run full UNRES optimization with constrains and frozen 2D c the same variables as soft pot. optimizatio c ipot=ipot0 maxmin=maxmin0 maxfun=maxfun0 c c check overlaps before calling full UNRES minim c call var_to_geom(nvar,var) call chainbuild call etotal(energy(0)) #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(0)) 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) c cdte time0=MPI_WTIME() cde change=reduce(var) cde if (check_var(var,info)) then cde write(iout,*) 'error before mask dist' cde call var_to_geom(nvar,var) cde call chainbuild cde call write_pdb(10000+in_pdb,'before mask dist',etot) cde endif cdte call minimize(etot,var,iretcode,nfun) cdte write(iout,*)'SUMSL MASK DIST return code is',iretcode, cdte & ' eval ',nfun cdte ieval=ieval+nfun cdte cdte time1=MPI_WTIME() cdte write (iout,'(a,f6.2,f8.2,a)') cdte & ' Time for mask dist min.',time1-time0, cdte & nfun/(time1-time0),' eval/s' cdte call flush(iout) if (debug) then call var_to_geom(nvar,var) call chainbuild call write_pdb(4000+in_pdb,'mask dist',etot) endif c c switch off freezing of 2D and c run full UNRES optimization with constrains c mask_r=.false. #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif cde change=reduce(var) cde if (check_var(var,info)) then cde write(iout,*) 'error before dist' cde call var_to_geom(nvar,var) cde call chainbuild cde call write_pdb(11000+in_pdb,'before dist',etot) cde endif call minimize(etot,var,iretcode,nfun) cde change=reduce(var) cde if (check_var(var,info)) then cde write(iout,*) 'error after dist',ico cde call var_to_geom(nvar,var) cde call chainbuild cde call write_pdb(12000+in_pdb+ico*1000,'after dist',etot) cde endif write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun ieval=ieval+nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #endif write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0, & nfun/(time1-time0),' eval/s' cde call etotal(energy(0)) cde 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 c c reset constrains c nhpb= nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 return end c-------------------------------------------------------- subroutine softreg implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif 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' c include 'COMMON.DISTFIT' integer iff(maxres) double precision time0,time1 double precision energy(0:n_ene),ee double precision var(maxvar) integer ieval c logical debug,ltest,fail character*50 linia c linia='test' debug=.true. in_pdb=0 c------------------------ c c freeze sec.elements c 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 c c store dist. constrains c 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 c c run soft pot. optimization c ipot=6 wang=3.0 maxmin=2000 maxfun=4000 call geom_to_var(nvar,var) #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 c c run full UNRES optimization with constrains and frozen 2D c the same variables as soft pot. optimizatio c ipot=ipot0 wang=wang0 maxmin=maxmin0 maxfun=maxfun0 #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL MASK DIST return code is',iretcode, & ' eval ',nfun ieval=nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 c c switch off constrains and c run full UNRES optimization with frozen 2D c c c reset constrains c nhpb_c=nhpb nhpb=nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun ieval=ieval+nfun #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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. c c run full UNRES optimization with constrains and NO frozen 2D c nhpb=nhpb_c link_start=1 link_end=nhpb maxfun=maxfun0/5 do ico=1,5 wstrain=wstrain0/ico #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #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 #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 c nhpb=nhpb0 link_start=1 link_end=nhpb wstrain=wstrain0 maxfun=maxfun0 c if (minim) then #ifdef MPI time0=MPI_WTIME() #else time0=tcpu() #endif call minimize(etot,var,iretcode,nfun) write(iout,*)'------------------------------------------------' write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, & '+ DIST eval',ieval #ifdef MPI time1=MPI_WTIME() #else time1=tcpu() #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 return end subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij) implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif 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' double precision time0,time1 double precision energy(0:n_ene),ee double precision var(maxvar) integer jdata(5),isec(maxres) c 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 c c cut strands at the ends c 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 cv call chainbuild cv call etotal(energy(0)) cv etot=energy(0) cv write(iout,*) nnt,nct,etot cv call write_pdb(ij*100,'first structure',etot) cv write(iout,*) 'N16 test',(jdata(i),i=1,5) c------------------------ c generate constrains c 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 cd 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 cd 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 cd 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 cv 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 cv time1=MPI_WTIME() cv write (iout,'(a,f6.2,f8.2,a)') cv & ' Time for dist min.',time1-time0, cv & nfun/(time1-time0),' eval/s' cv call var_to_geom(nvar,var) cv call chainbuild cv call write_pdb(ij*100+ico,'dist cons',etot) enddo c nhpb=nhpb0 call hpb_partition wstrain=wstrain0 maxfun=maxfun0 c cd print *,etot wscloc0=wscloc wscloc=10.0 call sc_move(nnt,nct,100,100d0,nft_sc,etot) wscloc=wscloc0 cv call chainbuild cv call etotal(energy(0)) cv etot=energy(0) cv call write_pdb(ij*100+10,'sc_move',etot) cd call intout cd print *,nft_sc,etot return end subroutine beta_zip(i1,i2,ieval,ij) implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif 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' double precision time0,time1 double precision energy(0:n_ene),ee double precision var(maxvar) character*10 test cv call chainbuild cv call etotal(energy(0)) cv etot=energy(0) cv write(test,'(2i5)') i1,i2 cv call write_pdb(ij*100,test,etot) cv write(iout,*) 'N17 test',i1,i2,etot,ij c c generate constrains c 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 cv 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 cv time1=MPI_WTIME() cv write (iout,'(a,f6.2,f8.2,a)') cv & ' Time for dist min.',time1-time0, cv & nfun/(time1-time0),' eval/s' c do not comment the next line call var_to_geom(nvar,var) cv call chainbuild cv call write_pdb(ij*100+ico,'dist cons',etot) enddo nhpb=nhpb0 call hpb_partition wstrain=wstrain0 maxfun=maxfun0 cv call etotal(energy(0)) cv etot=energy(0) cv write(iout,*) 'N17 test end',i1,i2,etot,ij return end