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_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 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 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