Merge branch 'devel' into AFM
[unres.git] / source / unres / src_MD-M-newcorr / together.F
diff --git a/source/unres/src_MD-M-newcorr/together.F b/source/unres/src_MD-M-newcorr/together.F
new file mode 100644 (file)
index 0000000..b0e0997
--- /dev/null
@@ -0,0 +1,1222 @@
+      Subroutine together
+c  feeds tasks for parallel processing
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'mpif.h'
+      real ran1,ran2
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.TIME1'
+      include 'COMMON.SETUP'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
+      real tcpu
+      double precision time_start,time_start_c,time0f,time0i
+      logical ovrtim,sync_iter,timeout,flag,timeout1
+      dimension muster(mpi_status_size)
+      dimension t100(0:100),indx(mxio)
+      dimension xout(maxvar),eout(mxch*(mxch+1)/2+1),ind(9)
+      dimension cout(2)
+      parameter (rad=1.745329252d-2)
+
+cccccccccccccccccccccccccccccccccccccccccccccccc
+      IF (ME.EQ.KING) THEN
+
+       time0f=MPI_WTIME()
+       ilastnstep=1
+       sync_iter=.false.
+       numch=1
+       nrmsdb=0
+       nrmsdb1=0
+       rmsdbc1c=rmsdbc1
+       nstep=0
+       call csa_read
+       call make_array
+
+       if(iref.ne.0) call from_int(1,0,idum)
+
+c To minimize input conformation (bank conformation)
+c Output to $mol.reminimized
+       if (irestart.lt.0) then
+        call read_bank(0,nft,cutdifr)
+        if (irestart.lt.-10) then
+         p_cut=nres*4.d0
+         call prune_bank(p_cut)
+         return
+        endif
+        call reminimize(jlee)
+        return
+       endif
+
+       if (irestart.eq.0) then
+        call initial_write
+        nbank=nconf
+        ntbank=nconf
+        if (ntbankm.eq.0) ntbank=0
+        nstep=0
+        nft=0
+        do i=1,mxio
+         ibank(i)=0
+         jbank(i)=0
+        enddo
+       else
+        call restart_write
+c!bankt        call read_bankt(jlee,nft,cutdifr)
+        call read_bank(jlee,nft,cutdifr)
+        call read_rbank(jlee,adif)
+        if(iref.ne.0) call from_int(1,0,idum)
+       endif
+
+       nstmax=nstmax+nstep
+       ntrial=n1+n2+n3+n4+n5+n6+n7+n8
+       ntry=ntrial+1
+       ntry=ntry*nseed
+
+c ntrial : number of trial conformations per seed.
+c ntry : total number of trial conformations including seed conformations.
+
+       idum2=-123
+c       imax=2**31-1
+       imax=huge(0)
+       ENDIF
+
+       call mpi_bcast(jend,1,mpi_integer,0,CG_COMM,ierr)
+cccccccccccccccccccccccccccccccccccccccc
+       do 300 jlee=1,jend
+cccccccccccccccccccccccccccccccccccccccc
+  331   continue 
+        IF (ME.EQ.KING) THEN
+        if(sync_iter) goto 333
+        idum=-  ran2(idum2)*imax
+        if(jlee.lt.jstart) goto 300
+
+C Restart the random number generator for conformation generation
+
+        if(irestart.gt.0) then
+         idum2=idum2+nstep
+         if(idum2.le.0) idum2=-idum2+1
+         idum=-  ran2(idum2)*imax
+        endif
+
+        idumm=idum
+        call vrndst(idumm)
+
+        open(icsa_seed,file=csa_seed,status="old")
+        write(icsa_seed,*) "jlee : ",jlee
+        close(icsa_seed)
+
+      call history_append
+      write(icsa_history,*) "number of procs is ",nodes
+      write(icsa_history,*) jlee,idum,idum2
+      close(icsa_history)
+
+cccccccccccccccccccccccccccccccccccccccccccccccc
+  333 icycle=0
+
+       call history_append
+        write(icsa_history,*) "nbank is ",nbank
+       close(icsa_history)
+
+      if(irestart.eq.1) goto 111
+      if(irestart.eq.2) then
+       icycle=0
+       do i=1,nbank
+        ibank(i)=1
+       enddo
+       do i=nbank+1,nbank+nconf
+        ibank(i)=0
+       enddo
+      endif
+
+c  start energy minimization
+      nconfr=max0(nconf+nadd,nodes-1)
+      if (sync_iter) nconf_in=0
+c  king-emperor - feed input and sort output
+       write (iout,*) "NCONF_IN",nconf_in
+       m=0
+       if (nconf_in.gt.0) then
+c al 7/2/00 - added possibility to read in some of the initial conformations
+        do m=1,nconf_in
+          read (intin,'(i5)',end=11,err=12) iconf
+   12     continue
+          write (iout,*) "write READ_ANGLES",iconf,m
+          call read_angles(intin,*11)
+          if (iref.eq.0) then
+            mm=m
+          else
+            mm=m+1
+          endif
+          do j=2,nres-1
+            dihang_in(1,j,1,mm)=theta(j+1)
+            dihang_in(2,j,1,mm)=phi(j+2)
+            dihang_in(3,j,1,mm)=alph(j)
+            dihang_in(4,j,1,mm)=omeg(j)
+          enddo
+        enddo ! m
+        goto 13
+   11   write (iout,*) nconf_in," conformations requested, but only",
+     &   m-1," found in the angle file."
+        nconf_in=m-1
+   13   continue
+        m=nconf_in
+        write (iout,*) nconf_in,
+     &    " initial conformations have been read in."
+       endif
+       if (iref.eq.0) then
+        if (nconfr.gt.nconf_in) then
+          call make_ranvar(nconfr,m,idum)
+          write (iout,*) nconfr-nconf_in,
+     &     " conformations have been generated randomly."
+        endif
+       else
+        nconfr=nconfr*2
+        call from_int(nconfr,m,idum)
+c       call from_pdb(nconfr,idum)
+       endif
+       write (iout,*) 'Exitted from make_ranvar nconfr=',nconfr
+       write (*,*) 'Exitted from make_ranvar nconfr=',nconfr
+       do m=1,nconfr
+          write (iout,*) 'Initial conformation',m
+          write(iout,'(8f10.4)') (rad2deg*dihang_in(1,j,1,m),j=2,nres-1)
+          write(iout,'(8f10.4)') (rad2deg*dihang_in(2,j,1,m),j=2,nres-1)
+          write(iout,'(8f10.4)') (rad2deg*dihang_in(3,j,1,m),j=2,nres-1)
+          write(iout,'(8f10.4)') (rad2deg*dihang_in(4,j,1,m),j=2,nres-1)
+       enddo 
+       write(iout,*)'Calling FEEDIN NCONF',nconfr
+       time1i=MPI_WTIME()
+       call feedin(nconfr,nft)
+       write (iout,*) ' Time for first bank min.',MPI_WTIME()-time1i
+       call  history_append
+        write(icsa_history,*) jlee,nft,nbank
+        write(icsa_history,851) (etot(i),i=1,nconfr)
+        write(icsa_history,850) (rmsn(i),i=1,nconfr)
+        write(icsa_history,850) (pncn(i),i=1,nconfr)
+        write(icsa_history,*)
+       close(icsa_history)
+      ELSE
+c To minimize input conformation (bank conformation)
+c Output to $mol.reminimized   
+       if (irestart.lt.0) then 
+        call reminimize(jlee)
+        return
+       endif
+       if (irestart.eq.1) goto 111
+c  soldier - perform energy minimization
+ 334   call minim_jlee
+      ENDIF
+
+ccccccccccccccccccccccccccccccccccc
+c need to syncronize all procs
+      call mpi_barrier(CG_COMM,ierr)
+      if (ierr.ne.0) then
+       print *, ' cannot synchronize MPI'
+       stop
+      endif
+ccccccccccccccccccccccccccccccccccc
+
+      IF (ME.EQ.KING) THEN
+
+c      print *,"ok after minim"
+      nstep=nstep+nconf
+      if(irestart.eq.2) then
+       nbank=nbank+nconf
+c      ntbank=ntbank+nconf
+       if(ntbank.gt.ntbankm) ntbank=ntbankm
+      endif
+c      print *,"ok before indexx"
+      if(iref.eq.0) then
+       call indexx(nconfr,etot,indx)
+      else
+c cc/al 7/6/00
+       do k=1,nconfr
+         indx(k)=k
+       enddo
+       call indexx(nconfr-nconf_in,rmsn(nconf_in+1),indx(nconf_in+1))
+       do k=nconf_in+1,nconfr
+         indx(k)=indx(k)+nconf_in
+       enddo
+c cc/al
+c       call indexx(nconfr,rmsn,indx)
+      endif
+c      print *,"ok after indexx"
+      do im=1,nconf
+       m=indx(im)
+       if (m.gt.mxio .or. m.lt.1) then
+         write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,' M',m
+         call mpi_abort(mpi_comm_world,ierror,ierrcode)
+       endif
+       jbank(im+nbank-nconf)=0
+       bene(im+nbank-nconf)=etot(m)
+       rene(im+nbank-nconf)=etot(m)
+c!bankt       btene(im)=etot(m)
+c
+       brmsn(im+nbank-nconf)=rmsn(m)
+       bpncn(im+nbank-nconf)=pncn(m)
+       rrmsn(im+nbank-nconf)=rmsn(m)
+       rpncn(im+nbank-nconf)=pncn(m)
+       if (im+nbank-nconf.gt.mxio .or. im+nbank-nconf.lt.1) then
+         write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,
+     &   ' NBANK',nbank,' NCONF',nconf,' IM+NBANK-NCONF',im+nbank-nconf
+         call mpi_abort(mpi_comm_world,ierror,ierrcode)
+       endif
+       do k=1,numch
+        do j=2,nres-1
+         do i=1,4
+          bvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
+          rvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
+c!bankt          btvar(i,j,k,im)=dihang(i,j,k,m)
+c
+         enddo
+        enddo
+       enddo
+       if(iref.eq.1) then
+        if(brmsn(im+nbank-nconf).gt.rmscut.or.
+     &     bpncn(im+nbank-nconf).lt.pnccut) ibank(im+nbank-nconf)=9
+       endif
+       if(vdisulf) then
+           bvar_ns(im+nbank-nconf)=ns-2*nss
+           k=0
+           do i=1,ns
+             j=1
+             do while( iss(i).ne.ihpb(j)-nres .and. 
+     &                 iss(i).ne.jhpb(j)-nres .and. j.le.nss)
+              j=j+1 
+             enddo
+             if (j.gt.nss) then            
+               k=k+1
+               bvar_s(k,im+nbank-nconf)=iss(i)
+             endif
+           enddo
+       endif
+       bvar_nss(im+nbank-nconf)=nss
+       do i=1,nss
+           bvar_ss(1,i,im+nbank-nconf)=ihpb(i)
+           bvar_ss(2,i,im+nbank-nconf)=jhpb(i)
+       enddo
+      enddo
+      ENDIF
+
+  111 continue
+
+      IF (ME.EQ.KING) THEN
+
+      call find_max
+      call find_min
+
+      call get_diff
+      if(nbank.eq.nconf.and.irestart.eq.0) then
+       adif=avedif
+      endif
+
+      cutdif=adif/cut1
+      ctdif1=adif/cut2
+
+cd      print *,"adif,xctdif,cutdifr"
+cd      print *,adif,xctdif,cutdifr
+       nst=ntotal/ntrial/nseed
+       xctdif=(cutdif/ctdif1)**(-1.0/nst)
+       if(irestart.ge.1) call estimate_cutdif(adif,xctdif,cutdifr)
+c       print *,"ok after estimate"
+
+      irestart=0
+
+       call write_rbank(jlee,adif,nft)
+       call write_bank(jlee,nft)
+c!bankt       call write_bankt(jlee,nft)
+c       call write_bank1(jlee)
+       call  history_append
+        write(icsa_history,*) "xctdif: ", xctdif,nst,adif/cut1,ctdif1
+        write(icsa_history,851) (bene(i),i=1,nbank)
+        write(icsa_history,850) (brmsn(i),i=1,nbank)
+        write(icsa_history,850) (bpncn(i),i=1,nbank)
+        close(icsa_history)
+  850 format(10f8.3)
+  851 format(5e15.6)
+
+      ifar=nseed/4*3+1
+      ifar=nseed+1
+      ENDIF
+    
+
+      finished=.false.
+      iter = 0
+      irecv = 0
+      isent =0
+      ifrom= 0
+      time0i=MPI_WTIME()
+      time1i=time0i
+      time_start_c=time0i
+      if (.not.sync_iter) then 
+        time_start=time0i
+        nft00=nft
+      else
+        sync_iter=.false.
+      endif
+      nft00_c=nft
+      nft0i=nft
+ccccccccccccccccccccccccccccccccccccccc
+      do while (.not. finished)
+ccccccccccccccccccccccccccccccccccccccc
+crc      print *,"iter ", iter,' isent=',isent
+
+      IF (ME.EQ.KING) THEN
+c  start energy minimization
+
+       if (isent.eq.0) then
+c  king-emperor - select seeds & make var & feed input
+cd        print *,'generating new conf',ntrial,MPI_WTIME()
+        call select_is(nseed,ifar,idum)
+
+        open(icsa_seed,file=csa_seed,status="old")
+        write(icsa_seed,39) 
+     &    jlee,icycle,nstep,(is(i),bene(is(i)),i=1,nseed)
+        close(icsa_seed)
+        call  history_append
+        write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,
+     *   ebmin,ebmax,nft,iuse,nbank,ntbank
+        close(icsa_history)
+
+         
+
+        call make_var(ntry,idum,iter)
+cd        print *,'new trial generated',ntrial,MPI_WTIME()
+           time2i=MPI_WTIME()
+           write (iout,'(a20,i4,f12.2)') 
+     &       'Time for make trial',iter+1,time2i-time1i
+       endif
+
+crc        write(iout,*)'1:Calling FEEDIN NTRY',NTRY,' ntrial',ntrial
+crc        call feedin(ntry,nft)
+
+       isent=isent+1
+       if (isent.ge.nodes.or.iter.gt.0)  then
+ct            print *,'waiting ',MPI_WTIME()
+            irecv=irecv+1
+            call recv(0,ifrom,xout,eout,ind,timeout)
+ct            print *,'   ',irecv,' received from',ifrom,MPI_WTIME()
+       else
+            ifrom=ifrom+1
+       endif
+
+ct            print *,'sending to',ifrom,MPI_WTIME()
+       call send(isent,ifrom,iter)
+ct            print *,isent,' sent ',MPI_WTIME()
+
+c store results -----------------------------------------------
+       if (isent.ge.nodes.or.iter.gt.0)  then
+         nft=nft+ind(3)
+         movernx(irecv)=iabs(ind(5))
+         call getx(ind,xout,eout,cout,rad,iw_pdb,irecv)
+         if(vdisulf) then
+             nss_out(irecv)=nss
+             do i=1,nss
+               iss_out(i,irecv)=ihpb(i)
+               jss_out(i,irecv)=jhpb(i)  
+             enddo
+         endif
+         if(iw_pdb.gt.0) 
+     &          call write_csa_pdb(xout,eout,nft,irecv,iw_pdb)
+       endif
+c--------------------------------------------------------------
+       if (isent.eq.ntry) then
+           time1i=MPI_WTIME()
+           write (iout,'(a18,f12.2,a14,f10.2)') 
+     &       'Nonsetup time     ',time1i-time_start_c,
+     &       ' sec, Eval/s =',(nft-nft00_c)/(time1i-time_start_c)
+           write (iout,'(a14,i4,f12.2,a14,f10.2)') 
+     &       'Time for iter ',iter+1,time1i-time0i,
+     &       ' sec, Eval/s =',(nft-nft0i)/(time1i-time0i)
+           time0i=time1i
+           nft0i=nft
+           cutdif=cutdif*xctdif
+           if(cutdif.lt.ctdif1) cutdif=ctdif1
+           if (iter.eq.0) then
+              print *,'UPDATING ',ntry-nodes+1,irecv
+              write(iout,*) 'UPDATING ',ntry-nodes+1
+              iter=iter+1
+c----------------- call update(ntry-nodes+1) -------------------
+              nstep=nstep+ntry-nseed-(nodes-1)
+              call refresh_bank(ntry-nodes+1)
+c!bankt              call refresh_bankt(ntry-nodes+1)
+           else
+c----------------- call update(ntry) ---------------------------
+              iter=iter+1
+              print *,'UPDATING ',ntry,irecv
+              write(iout,*) 'UPDATING ',ntry
+              nstep=nstep+ntry-nseed
+              call refresh_bank(ntry)
+c!bankt              call refresh_bankt(ntry)
+           endif         
+c----------------------------------------------------------------- 
+
+           call write_bank(jlee,nft)
+c!bankt           call write_bankt(jlee,nft)
+           call find_min
+
+           time1i=MPI_WTIME()
+           write (iout,'(a20,i4,f12.2)') 
+     &       'Time for refresh ',iter,time1i-time0i
+
+           if(ebmin.lt.estop) finished=.true.
+           if(icycle.gt.icmax) then
+               call write_bank1(jlee)
+               do i=1,nbank
+                 ibank(i)=2
+                 ibank(i)=1
+               enddo
+               nbank=nbank+nconf
+               if(nbank.gt.1000) then 
+                   finished=.true.
+               else
+crc                   goto 333
+                   sync_iter=.true.
+               endif
+           endif
+           if(nstep.gt.nstmax) finished=.true.
+
+           if(finished.or.sync_iter) then
+            do ij=1,nodes-1
+              call recv(1,ifrom,xout,eout,ind,timeout)
+              if (timeout) then
+                nstep=nstep+ij-1
+                print *,'ERROR worker is not responding'
+                write(iout,*) 'ERROR worker is not responding'
+                time1i=MPI_WTIME()-time_start_c
+                print *,'End of cycle, master time for ',iter,' iters ',
+     &             time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
+                write (iout,*) 'End of cycle, master time for ',iter,
+     &             ' iters ',time1i,' sec'
+                write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
+                print *,'UPDATING ',ij-1
+                write(iout,*) 'UPDATING ',ij-1
+                call flush(iout)
+                call refresh_bank(ij-1)
+c!bankt                call refresh_bankt(ij-1)
+                goto 1002
+              endif
+c              print *,'node ',ifrom,' finished ',ij,nft
+              write(iout,*) 'node ',ifrom,' finished ',ij,nft
+              call flush(iout)
+              nft=nft+ind(3)
+              movernx(ij)=iabs(ind(5))
+              call getx(ind,xout,eout,cout,rad,iw_pdb,ij)
+              if(vdisulf) then
+               nss_out(ij)=nss
+               do i=1,nss
+                 iss_out(i,ij)=ihpb(i)
+                 jss_out(i,ij)=jhpb(i)  
+               enddo
+              endif
+              if(iw_pdb.gt.0) 
+     &          call write_csa_pdb(xout,eout,nft,ij,iw_pdb)
+            enddo
+            nstep=nstep+nodes-1
+crc            print *,'---------bcast finished--------',finished
+            time1i=MPI_WTIME()-time_start_c
+            print *,'End of cycle, master time for ',iter,' iters ',
+     &             time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
+            write (iout,*) 'End of cycle, master time for ',iter,
+     &             ' iters ',time1i,' sec'
+            write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
+
+ctimeout            call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
+ctimeout            call mpi_bcast(sync_iter,1,mpi_logical,0,
+ctimeout     &                                              CG_COMM,ierr)
+            do ij=1,nodes-1 
+               tstart=MPI_WTIME()
+               call mpi_issend(finished,1,mpi_logical,ij,idchar,
+     &             CG_COMM,ireq,ierr)
+               call mpi_issend(sync_iter,1,mpi_logical,ij,idchar,
+     &             CG_COMM,ireq2,ierr)
+               flag=.false.  
+               timeout1=.false.
+               do while(.not. (flag .or. timeout1))
+                 call MPI_TEST(ireq2,flag,muster,ierr)
+                 tend1=MPI_WTIME()
+                 if(tend1-tstart.gt.60) then 
+                  print *,'ERROR worker ',ij,' is not responding'
+                  write(iout,*) 'ERROR worker ',ij,' is not responding'
+                  timeout1=.true.
+                 endif
+               enddo
+               if(timeout1) then
+                write(iout,*) 'worker ',ij,' NOT OK ',tend1-tstart
+                timeout=.true.
+               else
+                write(iout,*) 'worker ',ij,' OK ',tend1-tstart
+               endif
+            enddo
+            print *,'UPDATING ',nodes-1,ij
+            write(iout,*) 'UPDATING ',nodes-1
+            call refresh_bank(nodes-1)
+c!bankt            call refresh_bankt(nodes-1)
+ 1002       continue
+            call write_bank(jlee,nft)
+c!bankt            call write_bankt(jlee,nft)
+            call find_min
+
+            do i=0,mxmv  
+              do j=1,3
+                nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
+                nstatnx(i,j)=0
+              enddo
+            enddo
+
+            write(iout,*)'### Total stats:'
+            do i=0,mxmv
+             if(nstatnx_tot(i,1).ne.0) then
+              if (i.le.9) then
+              write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') 
+     &        '### N',i,' total=',nstatnx_tot(i,1),
+     &      ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',
+     &       (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
+              else
+              write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') 
+     &        '###N',i,' total=',nstatnx_tot(i,1),
+     &      ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',
+     &       (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
+              endif
+             else
+              if (i.le.9) then
+              write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') 
+     &          '### N',i,' total=',nstatnx_tot(i,1),
+     &          ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),
+     &          ' %acc',0.0
+              else
+              write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') 
+     &          '###N',i,' total=',nstatnx_tot(i,1),
+     &          ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),
+     &          ' %acc',0.0
+              endif
+             endif
+            enddo
+
+           endif
+           if(sync_iter) goto 331
+
+   39 format(2i3,i7,5(i4,f8.3,1x),19(/,13x,5(i4,f8.3,1x)))
+   40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
+   43 format(10i8)
+   44 format('jlee =',i3,':',4f10.1,' E =',f8.3,i7,i10)
+
+           isent=0
+           irecv=0
+       endif
+      ELSE
+c  soldier - perform energy minimization
+        call minim_jlee
+        print *,'End of minim, proc',me,'time ',MPI_WTIME()-time_start
+        write (iout,*) 'End of minim, proc',me,'time ',
+     &                  MPI_WTIME()-time_start
+        call flush(iout)
+ctimeout        call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
+ctimeout        call mpi_bcast(sync_iter,1,mpi_logical,0,CG_COMM,ierr)
+         call mpi_recv(finished,1,mpi_logical,0,idchar,
+     *                 CG_COMM,muster,ierr)             
+         call mpi_recv(sync_iter,1,mpi_logical,0,idchar,
+     *                 CG_COMM,muster,ierr)             
+        if(sync_iter) goto 331
+      ENDIF
+
+ccccccccccccccccccccccccccccccccccccccc
+      enddo
+ccccccccccccccccccccccccccccccccccccccc
+
+      IF (ME.EQ.KING) THEN
+        call  history_append
+        write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,
+     *  ebmin,ebmax,nft,iuse,nbank,ntbank
+
+        write(icsa_history,44) jlee,0.0,0.0,0.0,
+     &   0.0,ebmin,nstep,nft
+        write(icsa_history,*)
+       close(icsa_history)
+
+       time1i=MPI_WTIME()-time_start
+       print *,'End of RUN, master time ',
+     &             time1i,'sec, Eval/s ',(nft-nft00)/time1i
+       write (iout,*) 'End of RUN, master time  ',
+     &             time1i,' sec'
+       write (iout,*) 'Total eval/s ',(nft-nft00)/time1i
+
+       if(timeout) then 
+        write(iout,*) '!!!! ERROR worker was not responding'
+        write(iout,*) '!!!! cannot finish work normally'
+        write(iout,*) 'Processor0 is calling MPI_ABORT'
+        print *,'!!!! ERROR worker was not responding'
+        print *,'!!!! cannot finish work normally'
+        print *,'Processor0 is calling MPI_ABORT'
+        call flush(iout)
+        call mpi_abort(mpi_comm_world, 111, ierr)
+       endif
+      ENDIF
+
+cccccccccccccccccccccccccccccc
+  300 continue
+cccccccccccccccccccccccccccccc
+
+      return
+      end
+c-------------------------------------------------
+      subroutine feedin(nconf,nft)
+c  sends out starting conformations and receives results of energy minimization
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      include 'mpif.h'
+      dimension xin(maxvar),xout(maxvar),eout(mxch*(mxch+1)/2+1),
+     *          cout(2),ind(9),info(12)
+      dimension muster(mpi_status_size)
+      include 'COMMON.SETUP'
+      parameter (rad=1.745329252d-2)
+
+      print *,'FEEDIN: NCONF=',nconf
+      mm=0
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      if (nconf .lt. nodes-1) then
+        write (*,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',
+     &   nconf,nodes-1 
+        write (iout,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',
+     &   nconf,nodes-1 
+        call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      endif
+      do n=1,nconf
+c  pull out external and internal variables for next start
+        call putx(xin,n,rad)
+!        write (iout,*) 'XIN from FEEDIN N=',n
+!        write(iout,'(8f10.4)') (xin(j),j=1,nvar)
+        mm=mm+1
+        if (mm.lt.nodes) then
+c  feed task to soldier
+!       print *, ' sending input for start # ',n
+         info(1)=n
+         info(2)=-1
+         info(3)=0
+         info(4)=0
+         info(5)=0
+         info(6)=0
+         call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,
+     *                  ierr)
+         call mpi_send(xin,nvar,mpi_double_precision,mm,
+     *                  idreal,CG_COMM,ierr)
+        else
+c  find an available soldier
+         call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
+     *                 CG_COMM,muster,ierr)
+!        print *, ' receiving output from start # ',ind(1)
+         man=muster(mpi_source)
+c  receive final energies and variables
+         nft=nft+ind(3)
+         call mpi_recv(eout,1,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+!         print *,eout 
+#ifdef CO_BIAS
+         call mpi_recv(co,1,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+         write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
+#endif
+         call mpi_recv(xout,nvar,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+!         print *,nvar , ierr
+c  feed next task to soldier
+!        print *, ' sending input for start # ',n
+         info(1)=n
+         info(2)=-1
+         info(3)=0
+         info(4)=0
+         info(5)=0
+         info(6)=0
+         info(7)=0
+         info(8)=0
+         info(9)=0
+         call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,
+     *                  ierr)
+         call mpi_send(xin,nvar,mpi_double_precision,man,
+     *                  idreal,CG_COMM,ierr)
+c  retrieve latest results
+         call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
+         if(iw_pdb.gt.0) 
+     &        call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
+        endif
+      enddo
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c  no more input
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      do j=1,nodes-1
+c  wait for a soldier
+       call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
+     *               CG_COMM,muster,ierr)
+crc       if (ierr.ne.0) go to 30
+!      print *, ' receiving output from start # ',ind(1)
+       man=muster(mpi_source)
+c  receive final energies and variables
+       nft=nft+ind(3)
+       call mpi_recv(eout,1,
+     *               mpi_double_precision,man,idreal,
+     *               CG_COMM,muster,ierr)
+!       print *,eout
+#ifdef CO_BIAS
+         call mpi_recv(co,1,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+         write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
+#endif
+crc       if (ierr.ne.0) go to 30
+       call mpi_recv(xout,nvar,mpi_double_precision,
+     *               man,idreal,CG_COMM,muster,ierr)
+!       print *,nvar , ierr
+crc       if (ierr.ne.0) go to 30
+c  halt soldier
+       info(1)=0
+       info(2)=-1
+       info(3)=0 
+       info(4)=0
+       info(5)=0
+       info(6)=0
+       call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,
+     *                ierr)
+c  retrieve results
+       call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
+       if(iw_pdb.gt.0) 
+     &          call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
+      enddo
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      return
+   10 print *, ' dispatching error'
+      call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      return
+   20 print *, ' communication error'
+      call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      return
+   30 print *, ' receiving error'
+      call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      return
+      end
+cccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine getx(ind,xout,eout,cout,rad,iw_pdb,k)
+c  receives and stores data from soldiers
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.CONTACTS'
+      dimension ind(9),xout(maxvar),eout(mxch*(mxch+1)/2+1)
+cjlee
+      double precision przes(3),obr(3,3),cout(2)
+      logical non_conv
+cjlee
+      iw_pdb=2
+      if (k.gt.mxio .or. k.lt.1) then 
+        write (iout,*) 
+     &   'ERROR - dimensions of ANGMIN have been exceeded K=',k
+        call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      endif
+c  store ind()
+      do j=1,9
+       indb(k,j)=ind(j)
+      enddo
+c  store energies
+      etot(k)=eout(1)
+c  retrieve dihedral angles etc
+      call var_to_geom(nvar,xout)
+      do j=2,nres-1
+        dihang(1,j,1,k)=theta(j+1)
+        dihang(2,j,1,k)=phi(j+2)
+        dihang(3,j,1,k)=alph(j)
+        dihang(4,j,1,k)=omeg(j)
+      enddo
+      dihang(2,nres-1,1,k)=0.0d0
+cjlee
+      if(iref.eq.0) then 
+       iw_pdb=1
+cd       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4)') 
+cd     &      ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' mv ',
+cd     &      ind(5),ind(4)
+       return
+      endif
+      call chainbuild
+c     call dihang_to_c(dihang(1,1,1,k))
+c     call fitsq(rms,c(1,1),crefjlee(1,1),nres,przes,obr,non_conv)
+c     call fitsq(rms,c(1,2),crefjlee(1,2),nres-1,przes,obr,non_conv)
+c           call fitsq(rms,c(1,nstart_seq),crefjlee(1,nstart_sup),
+c    &                 nsup,przes,obr,non_conv)
+c      rmsn(k)=dsqrt(rms)
+
+       call rmsd_csa(rmsn(k))
+       call contact(.false.,ncont,icont,co)
+       pncn(k)=contact_fract(ncont,ncont_ref,icont,icont_ref)     
+
+cd       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a5
+cd     &      ,0pf5.2,a5,f5.1,a,f6.3,a4,i3,i4)') 
+cd     &    ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' rms ',
+cd     &    rmsn(k),' %NC ',pncn(k)*100,' cont.order',co,' mv ',
+cd     &    ind(5),ind(4)
+
+      if (rmsn(k).gt.rmscut.or.pncn(k).lt.pnccut) iw_pdb=0
+      return
+      end
+cccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine putx(xin,n,rad)
+c  gets starting variables
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      dimension xin(maxvar)
+
+c  pull out starting values for variables
+!       write (iout,*)'PUTX: N=',n
+      do m=1,numch
+!        write (iout,'(8f10.4)') (dihang_in(1,j,m,n),j=2,nres-1)
+!        write (iout,'(8f10.4)') (dihang_in(2,j,m,n),j=2,nres-1)
+!        write (iout,'(8f10.4)') (dihang_in(3,j,m,n),j=2,nres-1)
+!        write (iout,'(8f10.4)') (dihang_in(4,j,m,n),j=2,nres-1)
+       do j=2,nres-1
+        theta(j+1)=dihang_in(1,j,m,n)
+        phi(j+2)=dihang_in(2,j,m,n)
+        alph(j)=dihang_in(3,j,m,n)
+        omeg(j)=dihang_in(4,j,m,n)
+       enddo
+      enddo
+c  set up array of variables
+      call geom_to_var(nvar,xin)
+!       write (iout,*) 'xin in PUTX N=',n 
+!       call intout
+!       write (iout,'(8f10.4)') (xin(i),i=1,nvar) 
+      return
+      end
+c--------------------------------------------------------
+      subroutine putx2(xin,iff,n)
+c  gets starting variables
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      dimension xin(maxvar),iff(maxres)
+
+c  pull out starting values for variables
+      do m=1,numch
+       do j=2,nres-1
+        theta(j+1)=dihang_in2(1,j,m,n)
+        phi(j+2)=dihang_in2(2,j,m,n)
+        alph(j)=dihang_in2(3,j,m,n)
+        omeg(j)=dihang_in2(4,j,m,n)
+       enddo
+      enddo
+c  set up array of variables
+      call geom_to_var(nvar,xin)
+       
+      do i=1,nres
+        iff(i)=iff_in(i,n)
+      enddo
+      return
+      end
+
+c-------------------------------------------------------
+      subroutine prune_bank(p_cut)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'mpif.h'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.TIME1'
+      include 'COMMON.SETUP'
+c---------------------------
+c This subroutine prunes bank conformations using p_cut
+c---------------------------
+
+      nprune=0
+      nprune=nprune+1
+      m=1
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         dihang(i,j,k,nprune)=bvar(i,j,k,m)
+        enddo
+       enddo
+      enddo
+      bene(nprune)=bene(m)
+      brmsn(nprune)=brmsn(m)
+      bpncn(nprune)=bpncn(m) 
+
+      do m=2,nbank
+       ddmin=9.d190
+       do ip=1,nprune
+        call get_diff12(dihang(1,1,1,ip),bvar(1,1,1,m),diff) 
+        if(diff.lt.p_cut) goto 100
+        if(diff.lt.ddmin) ddmin=diff
+       enddo
+       nprune=nprune+1
+       do k=1,numch
+        do j=2,nres-1
+         do i=1,4
+          dihang(i,j,k,nprune)=bvar(i,j,k,m)
+         enddo
+        enddo
+       enddo
+       bene(nprune)=bene(m)
+       brmsn(nprune)=brmsn(m)
+       bpncn(nprune)=bpncn(m)
+  100  continue
+       write (iout,*) 'Pruning :',m,nprune,p_cut,ddmin
+      enddo
+      nbank=nprune
+      print *, 'Pruning :',m,nprune,p_cut
+      call write_bank(0,0)
+
+      return
+      end
+c-------------------------------------------------------
+
+      subroutine reminimize(jlee)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'mpif.h'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.TIME1'
+      include 'COMMON.SETUP'
+c---------------------------
+c This subroutine re-minimizes bank conformations:
+c---------------------------
+
+       ntry=nbank
+
+       call find_max
+       call find_min
+
+       if (me.eq.king) then
+        open(icsa_history,file=csa_history,status="old")
+         write(icsa_history,*) "Re-minimization",nodes,"nodes"
+         write(icsa_history,851) (bene(i),i=1,nbank)
+         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,
+     *   ebmin,ebmax,nft,iuse,nbank,ntbank
+        close(icsa_history)
+        do index=1,ntry
+         do k=1,numch
+          do j=2,nres-1
+           do i=1,4
+            dihang_in(i,j,k,index)=bvar(i,j,k,index)
+           enddo
+          enddo
+         enddo
+        enddo
+        nft=0
+        call feedin(ntry,nft)
+       else
+        call minim_jlee
+       endif
+
+       call find_max
+       call find_min
+
+       if (me.eq.king) then
+        do i=1,ntry
+         call replace_bvar(i,i)
+        enddo
+        open(icsa_history,file=csa_history,status="old")
+         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,
+     *   ebmin,ebmax,nft,iuse,nbank,ntbank
+         write(icsa_history,851) (bene(i),i=1,nbank)
+        close(icsa_history)
+        call write_bank_reminimized(jlee,nft)
+       endif
+
+   40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
+  851 format(5e15.6)
+  850 format(5e15.10)
+c  850 format(10f8.3)
+
+      return
+      end
+c-------------------------------------------------------
+      subroutine send(n,mm,it)
+c  sends out starting conformation for minimization
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'mpif.h'
+      dimension xin(maxvar),xout(maxvar),eout(mxch*(mxch+1)/2+1),
+     *          cout(2),ind(8),xin2(maxvar),iff(maxres),info(12)
+      dimension muster(mpi_status_size)
+      include 'COMMON.SETUP'
+      parameter (rad=1.745329252d-2)
+
+      if (isend2(n).eq.0) then
+c  pull out external and internal variables for next start
+        call putx(xin,n,rad)
+        info(1)=n
+        info(2)=it
+        info(3)=movenx(n)
+        info(4)=nss_in(n)
+        info(5)=parent(1,n)
+        info(6)=parent(2,n)
+
+        if (movenx(n).eq.14.or.movenx(n).eq.17) then
+          info(7)=idata(1,n)
+          info(8)=idata(2,n)
+        else if (movenx(n).eq.16) then
+          info(7)=idata(1,n)
+          info(8)=idata(2,n)
+          info(10)=idata(3,n)
+          info(11)=idata(4,n)
+          info(12)=idata(5,n)
+        else
+         info(7)=0
+         info(8)=0
+         info(10)=0
+         info(11)=0
+         info(12)=0
+        endif
+
+        if (movenx(n).eq.15) then
+         info(9)=parent(3,n)
+        else
+         info(9)=0
+        endif
+        call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,
+     *                  ierr)
+        call mpi_send(xin,nvar,mpi_double_precision,mm,
+     *                  idreal,CG_COMM,ierr)
+      else
+c  distfit & minimization for n7 move
+        info(1)=-n
+        info(2)=it
+        info(3)=movenx(n)
+        info(4)=nss_in(n)
+        info(5)=parent(1,n)
+        info(6)=parent(2,n)
+        info(7)=0
+        info(8)=0
+        info(9)=0
+        call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,
+     *                  ierr)
+        call putx2(xin,iff,isend2(n))
+        call mpi_send(xin,nvar,mpi_double_precision,mm,
+     *                  idreal,CG_COMM,ierr)
+        call mpi_send(iff,nres,mpi_integer,mm,
+     *                  idint,CG_COMM,ierr)
+        call putx(xin2,n,rad)
+        call mpi_send(xin2,nvar,mpi_double_precision,mm,
+     *                  idreal,CG_COMM,ierr)
+      endif 
+      if (vdisulf.and.nss_in(n).ne.0) then
+        call mpi_send(iss_in(1,n),nss_in(n),mpi_integer,mm,
+     *                  idint,CG_COMM,ierr)
+        call mpi_send(jss_in(1,n),nss_in(n),mpi_integer,mm,
+     *                  idint,CG_COMM,ierr)
+      endif
+      return
+      end
+c-------------------------------------------------
+
+      subroutine recv(ihalt,man,xout,eout,ind,tout)
+c  receives results of energy minimization
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'mpif.h'
+      dimension xin(maxvar),xout(maxvar),eout(mxch*(mxch+1)/2+1),
+     *          cout(2),ind(9),info(12)
+      dimension muster(mpi_status_size)
+      include 'COMMON.SETUP'
+      logical tout,flag
+      double precision twait,tstart,tend1
+      parameter(twait=600.0d0)
+
+c  find an available soldier
+       tout=.false.
+       flag=.false.
+       tstart=MPI_WTIME()
+       do while(.not. (flag .or. tout))
+         call MPI_IPROBE(mpi_any_source,idint,CG_COMM,flag, 
+     *            muster,ierr)
+         tend1=MPI_WTIME()
+         if(tend1-tstart.gt.twait .and. ihalt.eq.1) tout=.true.
+c_error         if(tend1-tstart.gt.twait) tout=.true.
+       enddo
+       if (tout) then 
+         write(iout,*) 'ERROR = timeout for recv ',tend1-tstart
+         call flush(iout)
+         return
+       endif
+       man=muster(mpi_source)        
+
+ctimeout         call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
+ctimeout     *                 CG_COMM,muster,ierr)
+!        print *, ' receiving output from start # ',ind(1)
+ct         print *,'receiving ',MPI_WTIME()
+ctimeout         man=muster(mpi_source)
+         call mpi_recv(ind,9,mpi_integer,man,idint,
+     *                 CG_COMM,muster,ierr)
+ctimeout
+c  receive final energies and variables
+         call mpi_recv(eout,1,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+!         print *,eout 
+#ifdef CO_BIAS
+         call mpi_recv(co,1,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+         write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
+#endif
+         call mpi_recv(xout,nvar,mpi_double_precision,
+     *                 man,idreal,CG_COMM,muster,ierr)
+!         print *,nvar , ierr
+         if(vdisulf) nss=ind(6)
+         if(vdisulf.and.nss.ne.0) then
+          call mpi_recv(ihpb,nss,mpi_integer,
+     *                 man,idint,CG_COMM,muster,ierr)         
+          call mpi_recv(jhpb,nss,mpi_integer,
+     *                 man,idint,CG_COMM,muster,ierr)
+         endif
+c  halt soldier
+       if(ihalt.eq.1) then 
+c        print *,'sending halt to ',man
+        write(iout,*) 'sending halt to ',man
+        info(1)=0
+        call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,ierr)
+       endif
+      return
+      end
+
+c---------------------------------------------------------- 
+      subroutine history_append 
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+
+#if defined(AIX) || defined(PGI)
+       open(icsa_history,file=csa_history,position="append")
+#else
+       open(icsa_history,file=csa_history,access="append")
+#endif
+       return
+       end