Merge branch 'devel' into feature-ga
[unres.git] / source / unres / src_MD / src / thread.F
diff --git a/source/unres/src_MD/src/thread.F b/source/unres/src_MD/src/thread.F
deleted file mode 100644 (file)
index 9f169a0..0000000
+++ /dev/null
@@ -1,549 +0,0 @@
-      subroutine thread_seq
-C Thread the sequence through a database of known structures
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-#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
-      include 'COMMON.INFO'
-      integer ThreadId,ThreadType,Kwita
-#endif
-      double precision varia(maxvar)
-      double precision przes(3),obr(3,3)
-      double precision time_for_thread
-      logical found_pattern,non_conv
-      character*32 head_pdb
-      double precision energia(0:n_ene)
-      n_ene_comp=nprint_ene
-C   
-C Body
-C
-#ifdef MPI
-      if (me.eq.king) then
-        do i=1,nctasks
-          nsave_part(i)=0
-        enddo
-      endif
-      nacc_tot=0
-#endif
-      Kwita=0
-      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
-cd    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
-C Find long enough chain in the database
-          ii=iran_num(1,nseq)
-          nres_t=nres_base(1,ii)
-C 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
-C If this point is reached, the pattern has not yet been examined.
-   10     continue
-c         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
-C 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)
-c           cref(j,i)=c(j,i)
-          enddo
-cd        write (iout,'(a,i4,3f10.5)') restyp(itype(i)),i,(c(j,i),j=1,3)
-        enddo
-cd      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
-cd             non_conv) 
-cd      write (iout,'(a,f10.5)') 
-cd   &  'Initial RMS deviation from reference structure:',rms
-        if (itype(nres).eq.21) 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.21) 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.)
-cd      print *,'Exit INT_FROM_CART.'
-cd      print *,'nhpb=',nhpb
-        do i=nss+1,nhpb
-          ii=ihpb(i)
-          jj=jhpb(i)
-          dhpb(i)=dist(ii,jj)
-c         write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
-        enddo
-c       stop 'End generate'
-C Generate SC conformations.
-        call sc_conf
-c       call intout
-#ifdef MPI
-cd      print *,'Processor:',me,': exit GEN_SIDE.'
-#else
-cd      print *,'Exit GEN_SIDE.'
-#endif
-C Calculate initial energy.
-        call chainbuild
-        call etotal(energia(0))
-        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
-C 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)
-c       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(0))
-c       call enerprint(energia(0))
-        link_end=link_end0
-cd      call chainbuild
-cd      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv) 
-cd      write (iout,'(a,f10.5)') 
-cd   &  '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)
-c        write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)') 
-c     &  ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
-c     &  (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
-C 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
-c--------------------------------------------------------------------------
-      subroutine write_thread_summary
-C Thread the sequence through a database of known structures
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-#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.NAMES'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.TIME1'
-#ifdef MPI
-      include 'COMMON.INFO'
-#endif
-      dimension ip(maxthread)
-      double precision energia(0:n_ene)
-      write (iout,'(30x,a/)') 
-     & '  *********** Summary threading statistics ************'
-      write (iout,'(a)') 'Initial energies:'
-      write (iout,'(a4,2x,a12,14a14,3a8)') 
-     & 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',
-     & 'RMSnat','NatCONT','NNCONT','RMS'
-C Energy sort patterns
-      do i=1,nthread
-        ip(i)=i
-      enddo
-      do i=1,nthread-1
-        enet=ener(n_ene-1,ip(i))
-        jj=i
-        do j=i+1,nthread
-          if (ener(n_ene-1,ip(j)).lt.enet) then
-            jj=j
-            enet=ener(n_ene-1,ip(j))
-          endif
-        enddo
-        if (jj.ne.i) then
-          ipj=ip(jj)
-          ip(jj)=ip(i)
-          ip(i)=ipj
-        endif
-      enddo
-      do ik=1,nthread
-        i=ip(ik)
-        ii=ipatt(1,i)
-        ist=nres_base(2,ii)+ipatt(2,i)
-        do kk=1,n_ene_comp
-          energia(i)=ener0(kk,i)
-        enddo
-        etot=ener0(n_ene_comp+1,i)
-        rmsnat=ener0(n_ene_comp+2,i)
-        rms=ener0(n_ene_comp+3,i)
-        frac=ener0(n_ene_comp+4,i)
-        frac_nn=ener0(n_ene_comp+5,i)
-
-        if (refstr) then 
-        write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') 
-     &  i,str_nam(ii),ist+1,
-     &  (energia(print_order(kk)),kk=1,nprint_ene),
-     &  etot,rmsnat,frac,frac_nn,rms
-        else
-        write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') 
-     &  i,str_nam(ii),ist+1,
-     &  (energia(print_order(kk)),kk=1,nprint_ene),etot
-        endif
-      enddo
-      write (iout,'(//a)') 'Final energies:'
-      write (iout,'(a4,2x,a12,17a14,3a8)') 
-     & 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',
-     & 'RMSnat','NatCONT','NNCONT','RMS'
-      do ik=1,nthread
-        i=ip(ik)
-        ii=ipatt(1,i)
-        ist=nres_base(2,ii)+ipatt(2,i)
-        do kk=1,n_ene_comp
-          energia(kk)=ener(kk,ik)
-        enddo
-        etot=ener(n_ene_comp+1,i)
-        rmsnat=ener(n_ene_comp+2,i)
-        rms=ener(n_ene_comp+3,i)
-        frac=ener(n_ene_comp+4,i)
-        frac_nn=ener(n_ene_comp+5,i)
-        write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') 
-     &  i,str_nam(ii),ist+1,
-     &  (energia(print_order(kk)),kk=1,nprint_ene),
-     &  etot,rmsnat,frac,frac_nn,rms
-      enddo
-      write (iout,'(/a/)') 'IEXAM array:'
-      write (iout,'(i5)') nexcl
-      do i=1,nexcl
-        write (iout,'(2i5)') iexam(1,i),iexam(2,i)
-      enddo
-      write (iout,'(/a,1pe14.4/a,1pe14.4/)') 
-     & 'Max. time for threading step ',max_time_for_thread,
-     & 'Average time for threading step: ',ave_time_for_thread
-      return
-      end
-c----------------------------------------------------------------------------
-      subroutine sc_conf
-C Sample (hopefully) optimal SC orientations given backcone conformation.
-      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'
-      double precision varia(maxvar)
-      common /srutu/ icall
-      double precision energia(0:n_ene)
-      logical glycine,fail
-      maxsample=10
-      link_end0=link_end
-      link_end=min0(link_end,nss)
-      do i=nnt,nct
-        if (itype(i).ne.10) then
-cd        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(0))
-      do isample=1,maxsample
-C 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(0))
-cd      write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))') 
-cd   &   'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
-cd   &   ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
-        e1=0.0d0
-        if (e0.le.e1) then
-          alph(ind_sc)=alph0
-          omeg(ind_sc)=omeg0 
-        else
-          e0=e1
-        endif
-      enddo
-      link_end=link_end0
-      return
-      end
-c---------------------------------------------------------------------------
-      subroutine write_stat_thread(ithread,ipattern,ist)
-      implicit real*8 (a-h,o-z)
-      include "DIMENSIONS"
-      include "COMMON.CONTROL"
-      include "COMMON.IOUNITS"
-      include "COMMON.THREAD"
-      include "COMMON.FFIELD"
-      include "COMMON.DBASE"
-      include "COMMON.NAMES"
-      double precision energia(0:n_ene)
-
-#if defined(AIX) || defined(PGI)
-      open(istat,file=statname,position='append')
-#else
-      open(istat,file=statname,access='append')
-#endif
-      do i=1,n_ene_comp
-        energia(i)=ener(i,ithread)
-      enddo
-      etot=ener(n_ene_comp+1,ithread)
-      rmsnat=ener(n_ene_comp+2,ithread)
-      rms=ener(n_ene_comp+3,ithread)
-      frac=ener(n_ene_comp+4,ithread)
-      frac_nn=ener(n_ene_comp+5,ithread)
-      write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') 
-     &  ithread,str_nam(ipattern),ist+1,
-     &  (energia(print_order(i)),i=1,nprint_ene),
-     &  etot,rmsnat,frac,frac_nn,rms
-      close (istat)
-      return
-      end