--- /dev/null
+ module csa
+!-----------------------------------------------------------------------------
+ use io_units
+ use names
+ use math
+ use random
+ use geometry_data, only: nres,rad2deg
+ use geometry
+ use energy
+ use MPI_
+ use csa_data
+ use compare
+ use io_config
+
+ implicit none
+!-----------------------------------------------------------------------------
+! commom.bank
+! common/varin/
+! real(kind=8),dimension(:,:,:,:),allocatable :: dihang_in !(mxang,maxres,mxch,mxio)
+ integer,dimension(:),allocatable :: nss_in !(mxio)
+ integer,dimension(:,:),allocatable :: iss_in,jss_in !(maxss,mxio)
+! common/minvar/
+ real(kind=8),dimension(:,:,:,:),allocatable :: dihang !(mxang,maxres,mxch,mxio)
+ real(kind=8),dimension(:),allocatable :: etot!,rmsn,pncn !(mxio)
+ integer,dimension(:),allocatable :: nss_out !(mxio)
+ integer,dimension(:,:),allocatable ::iss_out,jss_out !(maxss,mxio)
+! common/bank/
+! real(kind=8),dimension(:,:,:,:),allocatable :: bvar !(mxang,maxres,mxch,mxio)
+! real(kind=8),dimension(:),allocatable :: bene,rene,&
+! brmsn,rrmsn,bpncn,rpncn !(mxio)
+ integer,dimension(:),allocatable :: is,jbank !(mxio)
+ real(kind=8) :: avedif,difmin,ebmin,ebmax,ebmaxt!,&
+! dele,difcut,cutdif,rmscut,pnccut
+ real(kind=8),dimension(:,:),allocatable :: dij !(mxio,mxio)
+! common/bank_disulfid/
+! common/mvstat/
+ integer,dimension(:),allocatable :: movenx,movernx !(mxio)
+ integer,dimension(:,:),allocatable :: nstatnx,nstatnx_tot !(0:mxmv,3)
+ integer,dimension(:,:),allocatable :: indb !(mxio,9)
+ integer,dimension(:,:),allocatable :: parent !(3,mxio)
+! common/send2/
+ integer,dimension(:),allocatable :: isend2 !(mxio)
+ integer,dimension(:,:),allocatable :: iff_in !(maxres,mxio2)
+ integer,dimension(:,:,:,:),allocatable :: dihang_in2 !(mxang,maxres,mxch,mxio2)
+ integer,dimension(:,:),allocatable :: idata !(5,mxio)
+!-----------------------------------------------------------------------------
+! common.csa
+! integer :: irestart,ndiff
+! common/alphaa/
+ integer,dimension(:),allocatable :: ngroup !(mxgr)
+ integer,dimension(:,:,:),allocatable :: igroup !(3,mxang,mxgr)
+ integer :: ntotgr!,numch
+! common/csa_input/
+! common/dih_control/
+! real(kind=8) :: rdih_bias
+ logical :: ldih_bias
+!-----------------------------------------------------------------------------
+! common.distfit
+! COMMON /frag/
+ integer,dimension(:,:),allocatable :: bvar_frag !(mxio,6)
+ integer,dimension(:,:),allocatable :: hvar_frag,lvar_frag,svar_frag !(mxio,3)
+ integer,dimension(:,:),allocatable :: avar_frag !(mxio,5)
+!-----------------------------------------------------------------------------
+! commom.hairpin
+! common /spinka/
+ integer :: nharp_tot
+ integer,dimension(:),allocatable :: nharp_seed,nharp_use !(max_seed)
+ integer,dimension(:,:,:),allocatable :: iharp_seed !(4,maxres/3,max_seed)
+ integer,dimension(:,:,:),allocatable :: iharp_use !(0:4,maxres/3,max_seed)
+!-----------------------------------------------------------------------------
+! Maximum number of moves (n1-n8)
+ integer,parameter :: mxmv=18
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! bank.F
+!-----------------------------------------------------------------------------
+ subroutine refresh_bank(ntrial)
+
+! 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.VAR'
+! include 'COMMON.CONTROL'
+ character :: chacc
+ integer :: iaccn,ntrial
+ real(kind=8) :: l_diff(mxio),denep
+ integer :: i,j,n,m,i1,idmin
+ real(kind=8) :: del_ene
+
+ 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
+
+! loop over all newly obtained conformations
+ do n=1,ntrial
+ chacc=' '
+ iaccn=0
+ nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1
+!ccccccccccccccccccccccccccccccccccccccccccc
+!jlee
+ if(iref.ne.0) then
+ if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100
+ endif
+!jlee
+ if(etot(n).gt.ebmax) goto 100
+! Find the conformation closest to the conformation n in the bank
+ difmin=9.d9
+ do m=1,nbank
+ call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m))
+ if(l_diff(m).lt.difmin) then
+ difmin=l_diff(m)
+ idmin=m
+ endif
+ enddo
+
+ if(difmin.lt.cutdif) then
+! n is redundant to idmin
+ if(etot(n).lt.bene(idmin)) then
+ if(etot(n).lt.bene(idmin)-0.01d0) then
+ ibank(idmin)=0
+ jbank(idmin)=0
+ endif
+ denep=bene(idmin)-etot(n)
+ call replace_bvar(idmin,n)
+!rc Update dij
+ do i1=1,nbank
+ if (i1.ne.idmin) then
+ dij(i1,idmin)=l_diff(i1)
+ dij(idmin,i1)=l_diff(i1)
+ endif
+ enddo
+ chacc='c'
+ iaccn=idmin
+ nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1
+ if(idmin.eq.ibmax) call find_max
+ endif
+ else
+! got new conformation
+ del_ene=0.0d0
+ if(ebmax-ebmin.gt.del_ene) then
+ denep=ebmax-etot(n)
+ call replace_bvar(ibmax,n)
+!rc Update dij
+ do i1=1,nbank
+ if (i1.ne.ibmax) then
+ dij(i1,ibmax)=l_diff(i1)
+ dij(ibmax,i1)=l_diff(i1)
+ endif
+ enddo
+ chacc='f'
+ iaccn=ibmax
+ nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1
+ ibank(ibmax)=0
+ jbank(ibmax)=0
+ call find_max
+ else
+ if(del_ene.lt.0.0001) then
+ write (iout,*) 'ERROR in refresh_bank: '
+ write (iout,*) 'ebmax: ',ebmax
+ write (iout,*) 'ebmin: ',ebmin
+ write (iout,*) 'del_ene: ',del_ene
+!rc call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+!jp nbmax is never defined so condition below is always false
+! if(nbank.lt.nbmax) then
+! nbank=nbank+1
+! call replace_bvar(nbank,n)
+! ibank(nbank)=0
+! jbank(nbank)=0
+! else
+ call replace_bvar(ibmax,n)
+ ibank(ibmax)=0
+ jbank(ibmax)=0
+ call find_max
+! endif
+ endif
+ endif
+!ccccccccccccccccccccccccccccccccccccccccccc
+ 100 continue
+ if (iaccn.eq.0) then
+ if (iref.eq.0) then
+ write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)') &
+ indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
+ indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9)
+ else
+ write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0)') &
+ indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
+ indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
+ ' rms ',rmsn(n),' %NC ',pncn(n)*100
+ endif
+ else
+ if (iref.eq.0) then
+ write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,1x,a1,i4,0pf8.1,0pf8.1)') &
+ indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
+ indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
+ chacc,iaccn,difmin,denep
+ else
+ write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0,1x,a1,i4,0pf8.1,0pf8.1)') &
+ indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
+ indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
+ ' rms ',rmsn(n),' %NC ',pncn(n)*100,&
+ chacc,iaccn,difmin,denep
+ endif
+ endif
+ enddo
+! end of loop over all newly obtained conformations
+ do i=0,mxmv
+ if(nstatnx(i,1).ne.0) then
+ if (i.le.9) then
+ write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
+ '## N',i,' total=',nstatnx(i,1),&
+ ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
+ ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
+ else
+ write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
+ '##N',i,' total=',nstatnx(i,1),&
+ ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
+ ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
+ endif
+ else
+ if (i.le.9) then
+ write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
+ '## N',i,' total=',nstatnx(i,1),&
+ ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
+ ' %acc',0.0
+ else
+ write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
+ '##N',i,' total=',nstatnx(i,1),&
+ ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
+ ' %acc',0.0
+ endif
+ endif
+ enddo
+ call flush(iout)
+!rc Update dij
+!rc moved up, saves some get_diff12 calls
+!rc
+!rc do i1=1,nbank-1
+!rc do i2=i1+1,nbank
+!rc if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
+!rc call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
+!rc dij(i1,i2)=diff
+!rc dij(i2,i1)=diff
+!rc endif
+!rc enddo
+!rc enddo
+
+ do i=1,nbank
+ jbank(i)=1
+ enddo
+
+ return
+ end subroutine refresh_bank
+!-----------------------------------------------------------------------------
+ subroutine replace_bvar(iold,inew)
+
+ use control_data, only: vdisulf
+ use energy_data, only: ns,iss
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ integer :: iold,inew,ierror,ierrcode,i,j,k
+
+ if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1) &
+ then
+ write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold,&
+ ' INEW',inew
+ 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,iold)=dihang(i,j,k,inew)
+ enddo
+ enddo
+ enddo
+ bene(iold)=etot(inew)
+ brmsn(iold)=rmsn(inew)
+ bpncn(iold)=pncn(inew)
+
+ if(bene(iold).lt.ebmin) then
+ ebmin=bene(iold)
+ ibmin=iold
+ endif
+
+ if(vdisulf) then
+ bvar_nss(iold)=nss_out(inew)
+!d write(iout,*) 'SS BANK',iold,bvar_nss(iold)
+ do i=1,bvar_nss(iold)
+ bvar_ss(1,i,iold)=iss_out(i,inew)
+ bvar_ss(2,i,iold)=jss_out(i,inew)
+!d write(iout,*) 'SS',bvar_ss(1,i,iold)-nres,
+!d & bvar_ss(2,i,iold)-nres
+ enddo
+
+ bvar_ns(iold)=ns-2*bvar_nss(iold)
+!d write(iout,*) 'CYS #free ', bvar_ns(iold)
+ k=0
+ do i=1,ns
+ j=1
+ do while( iss(i).ne.iss_out(j,inew)-nres .and. &
+ iss(i).ne.jss_out(j,inew)-nres .and. &
+ j.le.nss_out(inew))
+ j=j+1
+ enddo
+ if (j.gt.nss_out(inew)) then
+ k=k+1
+ bvar_s(k,iold)=iss(i)
+ endif
+ enddo
+!d write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold))
+ endif
+
+ return
+ end subroutine replace_bvar
+!-----------------------------------------------------------------------------
+ subroutine save_is(ind)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+ integer :: ind,i,j,k,index,ierror,ierrcode
+
+ index=nbank+ind
+! print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind)
+ if (index.gt.mxio .or. index.lt.1 .or. &
+ is(ind).gt.mxio .or. is(ind).lt.1) then
+ write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index,&
+ ' IND',ind,' IS',is(ind)
+ 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,index)=bvar(i,j,k,is(ind))
+ enddo
+ enddo
+ enddo
+ bene(index)=bene(is(ind))
+ ibank(is(ind))=1
+
+ return
+ end subroutine save_is
+!-----------------------------------------------------------------------------
+ subroutine select_is(n,ifar,idum)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer,dimension(mxio) :: itag
+ real(kind=8),dimension(mxio) :: adiff
+ integer :: n,ifar,idum,i,iusesv,imade
+
+ iuse=0
+ do i=1,nbank
+ if(ibank(i).eq.0) then
+ iuse=iuse+1
+ itag(iuse)=i
+ endif
+ enddo
+ iusesv=iuse
+
+ if(iuse.eq.0) then
+ icycle=icycle+1
+ do i=1,nbank
+ if(ibank(i).eq.2) then
+ ibank(i)=1
+ else
+ ibank(i)=0
+ endif
+ enddo
+ imade=0
+ call get_is(idum,ifar,n,imade,0)
+!test3 call get_is_max(idum,ifar,n,imade,0)
+ else if(iuse.eq.n) then
+ do i=1,iuse
+ is(i)=itag(i)
+ call save_is(i)
+ enddo
+ else if(iuse.lt.n) then
+! if(icycle.eq.0) then
+! do i=1,n
+! ind=mod(i-1,iuse)+1
+! is(i)=itag(ind)
+! call save_is(i)
+! enddo
+! else
+! endif
+ do i=1,iuse
+ is(i)=itag(i)
+ call save_is(i)
+ enddo
+ imade=iuse
+! call get_is_ran(idum,n,imade,1)
+ call get_is(idum,ifar,n,imade,1)
+!test3 call get_is_max(idum,ifar,n,imade,1)
+! if(iusesv.le.n/10) then
+ if(iusesv.le.0) then
+ icycle=icycle+1
+ do i=1,nbank
+! if(ibank(i).eq.2) then
+! ibank(i)=1
+ if(ibank(i).ge.2) then
+ ibank(i)=ibank(i)-1
+ else
+ ibank(i)=0
+ endif
+ enddo
+ endif
+ else
+ imade=0
+ call get_is(idum,ifar,n,imade,0)
+!test3 call get_is_max(idum,ifar,n,imade,0)
+ endif
+ iuse=iusesv
+
+ return
+ end subroutine select_is
+!-----------------------------------------------------------------------------
+ subroutine get_is_ran(idum,n,imade,k)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! real(kind=4) :: ran1,ran2
+ integer,dimension(mxio) :: itag
+ real(kind=8),dimension(mxio) :: adiff
+ integer :: idum,n,imade,k,j,i,iran
+
+ do j=imade+1,n
+ iuse=0
+ do i=1,nbank
+ if(ibank(i).eq.k) then
+ iuse=iuse+1
+ itag(iuse)=i
+ endif
+ enddo
+ iran=iuse* ran1(idum)+1
+ is(j)=itag(iran)
+ call save_is(j)
+ enddo
+
+ return
+ end subroutine get_is_ran
+!-----------------------------------------------------------------------------
+ subroutine get_is(idum,ifar,n,imade,k)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! real(kind=4) :: ran1,ran2
+ integer,dimension(mxio) :: itag
+ real(kind=8),dimension(mxio) :: adiff
+ integer :: idum,ifar,n,imade,k,i,iran
+
+ iuse=0
+ do i=1,nbank
+ if(ibank(i).eq.k) then
+ iuse=iuse+1
+ itag(iuse)=i
+ endif
+ enddo
+ iran=iuse* ran1(idum)+1
+ imade=imade+1
+ is(imade)=itag(iran)
+ call save_is(imade)
+
+ do i=imade+1,ifar-1
+ if(icycle.eq.-1) then
+ call select_iseed_max(i,k)
+ else
+ call select_iseed_min(i,k)
+!test4 call select_iseed_max(i,k)
+ endif
+ call save_is(i)
+ enddo
+
+ do i=ifar,n
+ call select_iseed_far(i,k)
+ call save_is(i)
+ enddo
+
+ return
+ end subroutine get_is
+!-----------------------------------------------------------------------------
+ subroutine select_iseed_max(imade1,ik)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer,dimension(mxio) :: itag
+ real(kind=8),dimension(mxio) :: adiff
+ integer :: imade1,ik,i,n,imade,m,itagi
+ real(kind=8) :: difmax,diff,emax,benei,diffmn
+
+ iuse=0
+ avedif=0.d0
+ difmax=0.d0
+ do n=1,nbank
+ if(ibank(n).eq.ik) then
+ iuse=iuse+1
+ diffmn=9.d190
+ do imade=1,imade1-1
+! m=nbank+imade
+! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
+ m=is(imade)
+ diff=dij(n,m)
+ if(diff.lt.diffmn) diffmn=diff
+ enddo
+ if(diffmn.gt.difmax) difmax=diffmn
+ adiff(iuse)=diffmn
+ itag(iuse)=n
+ avedif=avedif+diffmn
+ endif
+ enddo
+
+ avedif=avedif/iuse
+! avedif=(avedif+difmax)/2
+ emax=-9.d190
+ do i=1,iuse
+ if(adiff(i).ge.avedif) then
+ itagi=itag(i)
+ benei=bene(itagi)
+ if(benei.gt.emax) then
+ emax=benei
+ is(imade1)=itagi
+ endif
+ endif
+ enddo
+
+ if(ik.eq.0) iuse=iuse-1
+
+ return
+ end subroutine select_iseed_max
+!-----------------------------------------------------------------------------
+ subroutine select_iseed_min(imade1,ik)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer,dimension(mxio) :: itag
+ real(kind=8),dimension(mxio) :: adiff
+ integer :: imade1,ik,n,imade,m,i,itagi
+ real(kind=8) :: difmax,diff,diffmn,emin,benei
+
+ iuse=0
+ avedif=0.d0
+ difmax=0.d0
+ do n=1,nbank
+ if(ibank(n).eq.ik) then
+ iuse=iuse+1
+ diffmn=9.d190
+ do imade=1,imade1-1
+! m=nbank+imade
+! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
+ m=is(imade)
+ diff=dij(n,m)
+ if(diff.lt.diffmn) diffmn=diff
+ enddo
+ if(diffmn.gt.difmax) difmax=diffmn
+ adiff(iuse)=diffmn
+ itag(iuse)=n
+ avedif=avedif+diffmn
+ endif
+ enddo
+
+ avedif=avedif/iuse
+! avedif=(avedif+difmax)/2
+ emin=9.d190
+ do i=1,iuse
+! print *,"i, adiff(i),avedif : ",i,adiff(i),avedif
+ if(adiff(i).ge.avedif) then
+ itagi=itag(i)
+ benei=bene(itagi)
+! print *,"i, benei,emin : ",i,benei,emin
+ if(benei.lt.emin) then
+ emin=benei
+ is(imade1)=itagi
+ endif
+ endif
+ enddo
+
+ if(ik.eq.0) iuse=iuse-1
+
+! print *, "exiting select_iseed_min",is(imade1)
+
+ return
+ end subroutine select_iseed_min
+!-----------------------------------------------------------------------------
+ subroutine select_iseed_far(imade1,ik)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: imade1,ik,n,imade,m
+ real(kind=8) :: dmax,diffmn,diff
+
+ dmax=-9.d190
+ do n=1,nbank
+ if(ibank(n).eq.ik) then
+ diffmn=9.d190
+ do imade=1,imade1-1
+! m=nbank+imade
+! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
+ m=is(imade)
+ diff=dij(n,m)
+ if(diff.lt.diffmn) diffmn=diff
+ enddo
+ endif
+ if(diffmn.gt.dmax) then
+ dmax=diffmn
+ is(imade1)=n
+ endif
+ enddo
+
+ return
+ end subroutine select_iseed_far
+!-----------------------------------------------------------------------------
+ subroutine find_min
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: i
+ real(kind=8) :: benei
+
+ ebmin=9.d190
+
+ do i=1,nbank
+ benei=bene(i)
+ if(benei.lt.ebmin) then
+ ebmin=benei
+ ibmin=i
+ endif
+ enddo
+
+ return
+ end subroutine find_min
+!-----------------------------------------------------------------------------
+ subroutine find_max
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: i
+ real(kind=8) :: benei
+
+ ebmax=-9.d190
+
+ do i=1,nbank
+ benei=bene(i)
+ if(benei.gt.ebmax) then
+ ebmax=benei
+ ibmax=i
+ endif
+ enddo
+
+ return
+ end subroutine find_max
+!-----------------------------------------------------------------------------
+ subroutine get_diff
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: i,i1,i2
+ real(kind=8) :: tdiff,difmin,diff
+
+ tdiff=0.d0
+ difmin=9.d190
+ do i1=1,nbank-1
+ do i2=i1+1,nbank
+ if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
+ call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
+ dij(i1,i2)=diff
+ dij(i2,i1)=diff
+ else
+ diff=dij(i1,i2)
+ endif
+ tdiff=tdiff+diff
+ if(diff.lt.difmin) difmin=diff
+ enddo
+ dij(i1,i1)=0.0
+ enddo
+
+ do i=1,nbank
+ jbank(i)=1
+ enddo
+
+ avedif=tdiff/nbank/(nbank-1)*2
+
+ return
+ end subroutine get_diff
+!-----------------------------------------------------------------------------
+ subroutine estimate_cutdif(adif,xct,cutdifr)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: nexp
+ real(kind=8) :: adif,xct,cutdifr,ctdif1,exponent
+
+ ctdif1=adif/cut2
+
+ exponent = cutdifr*cut1/adif
+ exponent = dlog(exponent)/dlog(xct)
+
+ nexp=exponent+0.25
+ cutdif= adif/cut1*xct**nexp
+ if(cutdif.lt.ctdif1) cutdif=ctdif1
+
+ return
+ end subroutine estimate_cutdif
+!-----------------------------------------------------------------------------
+ subroutine get_is_max(idum,ifar,n,imade,k)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+ integer :: idum,ifar,n,imade,k,i,j
+ real(kind=8) :: emax
+
+ do i=imade+1,n
+ emax=-9.d190
+ do j=1,nbank
+ if(ibank(j).eq.k .and. bene(j).gt.emax) then
+ emax=bene(j)
+ is(i)=j
+ endif
+ enddo
+ call save_is(i)
+ enddo
+
+ return
+ end subroutine get_is_max
+!-----------------------------------------------------------------------------
+! csa.f
+!-----------------------------------------------------------------------------
+ subroutine make_array
+
+ use energy_data, only: itype
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CSA'
+ integer :: k,j,i,indg
+!cccccccccccccccccccccccc
+! Level-2: group
+!cccccccccccccccccccccccc
+
+ indg=0
+ do k=1,numch
+!cccccccccccccccccccccccccccccccccccccccc
+! Groups the THETAs and the GAMMAs
+ do j=2,nres-1
+ indg=indg+1
+ if (j.lt.nres-1) then
+ ngroup(indg)=2
+ else
+ ngroup(indg)=1
+ endif
+ do i=1,ngroup(indg)
+ igroup(1,i,indg)=i
+ igroup(2,i,indg)=j
+ igroup(3,i,indg)=k
+ enddo
+ enddo
+!cccccccccccccccccccccccccccccccccccccccc
+ enddo
+! Groups the ALPHAs and the BETAs
+ do k=1,numch
+ do j=2,nres-1
+ if(itype(j).ne.10) then
+ indg=indg+1
+ ngroup(indg)=2
+ do i=1,ngroup(indg)
+ igroup(1,i,indg)=i+2
+ igroup(2,i,indg)=j
+ igroup(3,i,indg)=k
+ enddo
+ endif
+ enddo
+ enddo
+
+ ntotgr=indg
+ write(iout,*)
+ write(iout,*) "# of groups: ",ntotgr
+ do i=1,ntotgr
+ write(iout,41) i,ngroup(i),((igroup(k,j,i),k=1,3),j=1,ngroup(i))
+ enddo
+! close(iout)
+
+ 40 format(i3,3x,3i3)
+ 41 format(2i3,3x,6(3i3,2x))
+
+ return
+ end subroutine make_array
+!-----------------------------------------------------------------------------
+ subroutine make_ranvar(n,m,idum)
+
+ use geometry_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.BANK'
+ integer :: n,m,j,idum,itrial,jeden
+
+! al m=0
+ print *,'HOHOHOHO Make_RanVar!!!!!',n,m
+ itrial=0
+ do while(m.lt.n .and. itrial.le.10000)
+ itrial=itrial+1
+ jeden=1
+ call gen_rand_conf(jeden,*10)
+! call intout
+ m=m+1
+ do j=2,nres-1
+ dihang_in(1,j,1,m)=theta(j+1)
+ dihang_in(2,j,1,m)=phi(j+2)
+ dihang_in(3,j,1,m)=alph(j)
+ dihang_in(4,j,1,m)=omeg(j)
+ enddo
+ dihang_in(2,nres-1,1,m)=0.0d0
+ goto 20
+ 10 write (iout,*) 'Failed to generate conformation #',m+1,&
+ ' itrial=',itrial
+ 20 continue
+ enddo
+ print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
+
+ return
+ end subroutine make_ranvar
+!-----------------------------------------------------------------------------
+ subroutine make_ranvar_reg(n,idum)
+
+ use geometry_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.BANK'
+! include 'COMMON.GEO'
+ integer :: n,idum,j,m,itrial,jeden
+
+ m=0
+ print *,'HOHOHOHO Make_RanVar!!!!!'
+ itrial=0
+ do while(m.lt.n .and. itrial.le.10000)
+ itrial=itrial+1
+ jeden=1
+ call gen_rand_conf(jeden,*10)
+! call intout
+ m=m+1
+ do j=2,nres-1
+ dihang_in(1,j,1,m)=theta(j+1)
+ dihang_in(2,j,1,m)=phi(j+2)
+ dihang_in(3,j,1,m)=alph(j)
+ dihang_in(4,j,1,m)=omeg(j)
+ if(m.le.n*0.1) then
+ dihang_in(1,j,1,m)=90.0*deg2rad
+ dihang_in(2,j,1,m)=50.0*deg2rad
+ endif
+ enddo
+ dihang_in(2,nres-1,1,m)=0.0d0
+ goto 20
+ 10 write (iout,*) 'Failed to generate conformation #',m+1,&
+ ' itrial=',itrial
+ 20 continue
+ enddo
+ print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
+
+ return
+ end subroutine make_ranvar_reg
+!-----------------------------------------------------------------------------
+! diff12.f
+!-----------------------------------------------------------------------------
+ subroutine get_diff12(aarray,barray,diff)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+ integer :: k,j,i
+ real(kind=8),dimension(mxang,nres,mxch) :: aarray,barray !(mxang,maxres,mxch)
+ real(kind=8) :: diff,dif
+
+ diff=0.d0
+ do k=1,numch
+ do j=2,nres-1
+! do i=1,4
+! do i=1,2
+ do i=1,ndiff
+ dif=rad2deg*dabs(aarray(i,j,k)-barray(i,j,k))
+ if(dif.gt.180.) dif=360.-dif
+ if (dif.gt.diffcut) diff=diff+dif
+ enddo
+ enddo
+ enddo
+
+ return
+ end subroutine get_diff12
+!-----------------------------------------------------------------------------
+! indexx.f
+!-----------------------------------------------------------------------------
+ subroutine indexx(n,arr,indx)
+
+! implicit real*8 (a-h,o-z)
+ INTEGER :: n,indx(n)
+ REAL(kind=8) :: arr(n)
+! PARAMETER (M=7,NSTACK=50)
+ integer,PARAMETER :: M=7,NSTACK=500
+ INTEGER :: i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
+ REAL(kind=8) :: a
+
+ do 11 j=1,n
+ indx(j)=j
+11 continue
+ jstack=0
+ l=1
+ ir=n
+1 if(ir-l.lt.M)then
+ do 13 j=l+1,ir
+ indxt=indx(j)
+ a=arr(indxt)
+ do 12 i=j-1,1,-1
+ if(arr(indx(i)).le.a)goto 2
+ indx(i+1)=indx(i)
+12 continue
+ i=0
+2 indx(i+1)=indxt
+13 continue
+ if(jstack.eq.0)return
+ ir=istack(jstack)
+ l=istack(jstack-1)
+ jstack=jstack-2
+ else
+ k=(l+ir)/2
+ itemp=indx(k)
+ indx(k)=indx(l+1)
+ indx(l+1)=itemp
+ if(arr(indx(l+1)).gt.arr(indx(ir)))then
+ itemp=indx(l+1)
+ indx(l+1)=indx(ir)
+ indx(ir)=itemp
+ endif
+ if(arr(indx(l)).gt.arr(indx(ir)))then
+ itemp=indx(l)
+ indx(l)=indx(ir)
+ indx(ir)=itemp
+ endif
+ if(arr(indx(l+1)).gt.arr(indx(l)))then
+ itemp=indx(l+1)
+ indx(l+1)=indx(l)
+ indx(l)=itemp
+ endif
+ i=l+1
+ j=ir
+ indxt=indx(l)
+ a=arr(indxt)
+3 continue
+ i=i+1
+ if(arr(indx(i)).lt.a)goto 3
+4 continue
+ j=j-1
+ if(arr(indx(j)).gt.a)goto 4
+ if(j.lt.i)goto 5
+ itemp=indx(i)
+ indx(i)=indx(j)
+ indx(j)=itemp
+ goto 3
+5 indx(l)=indx(j)
+ indx(j)=indxt
+ jstack=jstack+2
+ if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
+ if(ir-i+1.ge.j-l)then
+ istack(jstack)=ir
+ istack(jstack-1)=i
+ ir=j-1
+ else
+ istack(jstack)=j-1
+ istack(jstack-1)=l
+ l=i
+ endif
+ endif
+ goto 1
+ end subroutine indexx
+! (C) Copr. 1986-92 Numerical Recipes Software *11915aZ%.
+!-----------------------------------------------------------------------------
+! minim_jlee.F
+!-----------------------------------------------------------------------------
+ subroutine minim_jlee
+
+ use minim_data
+ use MPI_data
+ use energy_data
+ use compare_data
+ use control_data
+ use geometry_data, only: nvar,nphi
+ use geometry, only:dist
+ use energy, only:fdum
+ use control, only:init_int_table
+ use minimm, only:sumsl,deflt
+! controls minimization and sorting routines
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MINIM'
+! include 'COMMON.CONTROL'
+ include 'mpif.h'
+ integer,parameter :: liv=60
+ integer :: lv
+! external func,gradient!,fdum !use minim & energy
+! real(kind=4) :: ran1,ran2,ran3
+! include 'COMMON.SETUP'
+! include 'COMMON.GEO'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.DISTFIT'
+! include 'COMMON.CHAIN'
+ integer,dimension(mpi_status_size) :: muster
+ real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
+ real(kind=8),dimension(6*nres) :: var2 !(maxvar) (maxvar=6*maxres)
+ integer,dimension(nres) :: iffr !(maxres)
+ integer,dimension((nres-1)*(nres-2)/2) :: ihpbt,jhpbt !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2)
+ real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres)
+!el real(kind=8),dimension(1:lv+1) :: v
+ real(kind=8) :: energia(0:n_ene),time0s,time1s
+ integer,dimension(9) :: indx
+ integer,dimension(12) :: info
+ integer,dimension(liv) :: iv
+ integer :: idum(1)
+ real(kind=8) :: rdum(1)
+ integer,dimension(2,12*nres) :: icont_ !(2,maxcont)(maxcont=12*maxres)
+ logical :: fail !check_var,
+ integer :: iloop(2)
+!el common /przechowalnia/ v
+ integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode
+ real(kind=8) :: rad,eee,etot !,fdum
+!el from subroutine parmread
+! Define the constants of the disulfide bridge
+! Old arbitrary potential
+ real(kind=8),parameter :: dbr=4.20D0
+ real(kind=8),parameter :: fbr=3.30D0
+!-----------------
+ lv=77+(6*nres)*(6*nres+17)/2 !77+maxvar*(maxvar+17)/2 (maxvar=6*maxres)
+ data rad /1.745329252d-2/
+! receive # of start
+! print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
+! & ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
+ if (.not. allocated(v)) allocate(v(1:lv))
+ nhpb0=nhpb
+ 10 continue
+ time0s=MPI_WTIME()
+! print *, 'MINIM_JLEE: ',me,' is waiting'
+ call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,&
+ muster,ierr)
+ time1s=MPI_WTIME()
+ write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
+ call flush(iout)
+ n=info(1)
+! print *, 'MINIM_JLEE: ',me,' received: ',n
+
+!rc if (ierr.ne.0) go to 100
+! if # = 0, return
+ if (n.eq.0) then
+ write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
+ call flush(iout)
+ return
+ endif
+
+ nfun=0
+ IF (n.lt.0) THEN
+ call mpi_recv(var,nvar,mpi_double_precision,&
+ king,idreal,CG_COMM,muster,ierr)
+ call mpi_recv(iffr,nres,mpi_integer,&
+ king,idint,CG_COMM,muster,ierr)
+ call mpi_recv(var2,nvar,mpi_double_precision,&
+ king,idreal,CG_COMM,muster,ierr)
+ ELSE
+! receive initial values of variables
+ call mpi_recv(var,nvar,mpi_double_precision,&
+ king,idreal,CG_COMM,muster,ierr)
+!rc if (ierr.ne.0) go to 100
+ ENDIF
+
+ if(vdisulf.and.info(2).ne.-1) then
+ if(info(4).ne.0)then
+ call mpi_recv(ihpbt,info(4),mpi_integer,&
+ king,idint,CG_COMM,muster,ierr)
+ call mpi_recv(jhpbt,info(4),mpi_integer,&
+ king,idint,CG_COMM,muster,ierr)
+ endif
+ endif
+
+ IF (n.lt.0) THEN
+ n=-n
+ nhpb=nhpb0
+ link_start=1
+ link_end=nhpb
+ call init_int_table
+ call contact_cp(var,var2,iffr,nfun,n)
+ ENDIF
+
+ if(vdisulf.and.info(2).ne.-1) then
+ nss=0
+ if(info(4).ne.0)then
+!d write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
+ call var_to_geom(nvar,var)
+ call chainbuild
+ do i=1,info(4)
+ if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
+ nss=nss+1
+ ihpb(nss)=ihpbt(i)
+ jhpb(nss)=jhpbt(i)
+!d write(iout,*) 'SS mv=',info(3),
+!d & ihpb(nss)-nres,jhpb(nss)-nres,
+!d & dist(ihpb(nss),jhpb(nss))
+ dhpb(nss)=dbr
+ forcon(nss)=fbr
+ else
+!d write(iout,*) 'rm SS mv=',info(3),
+!d & ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
+ endif
+ enddo
+ endif
+ nhpb=nss
+ link_start=1
+ link_end=nhpb
+ call init_int_table
+ endif
+
+ if (info(3).eq.14) then
+ write(iout,*) 'calling local_move',info(7),info(8)
+ call local_move_init(.false.)
+ call var_to_geom(nvar,var)
+ call local_move(info(7),info(8),20d0,50d0)
+ call geom_to_var(nvar,var)
+ endif
+
+
+ if (info(3).eq.16) then
+ write(iout,*) 'calling beta_slide',info(7),info(8),&
+ info(10), info(11), info(12)
+ call var_to_geom(nvar,var)
+ call beta_slide(info(7),info(8),info(10),info(11),info(12), &
+ nfun,n)
+ call geom_to_var(nvar,var)
+ endif
+
+
+ if (info(3).eq.17) then
+ write(iout,*) 'calling beta_zip',info(7),info(8)
+ call var_to_geom(nvar,var)
+ call beta_zip(info(7),info(8),nfun,n)
+ call geom_to_var(nvar,var)
+ endif
+
+
+!rc overlap test
+
+ if (overlapsc) then
+
+ call var_to_geom(nvar,var)
+ call chainbuild
+ call etotal(energia)
+ nfun=nfun+1
+ if (energia(1).eq.1.0d20) then
+ info(3)=-info(3)
+ write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
+ call overlap_sc(fail)
+ if(.not.fail) then
+ call geom_to_var(nvar,var)
+ call etotal(energia)
+ nfun=nfun+1
+ write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
+ else
+ v(10)=1.0d20
+ iv(1)=-1
+ goto 201
+ endif
+ endif
+ endif
+
+ if (searchsc) then
+ call var_to_geom(nvar,var)
+ call sc_move(2,nres-1,1,10d0,nft_sc,etot)
+ call geom_to_var(nvar,var)
+!d write(iout,*) 'sc_move',nft_sc,etot
+ endif
+
+ if (check_var(var,info)) then
+ v(10)=1.0d21
+ iv(1)=6
+ goto 201
+ endif
+
+
+!rc
+
+! write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
+! write (iout,'(8f10.4)') (var(i),i=1,nvar)
+! write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
+! write (*,'(8f10.4)') (var(i),i=1,nvar)
+
+ do i=1,nvar
+ garbage(i)=var(i)
+ enddo
+
+ call deflt(2,iv,liv,lv,v)
+! 12 means fresh start, dont call deflt
+ iv(1)=12
+! max num of fun calls
+ if (maxfun.eq.0) maxfun=500
+ iv(17)=maxfun
+! max num of iterations
+ if (maxmin.eq.0) maxmin=1000
+ iv(18)=maxmin
+! controls output
+ iv(19)=2
+! selects output unit
+!d iv(21)=iout
+ iv(21)=0
+! 1 means to print out result
+ iv(22)=0
+!d iv(22)=1
+! 1 means to print out summary stats
+ iv(23)=0
+! 1 means to print initial x and d
+ iv(24)=0
+
+! if(me.eq.3.and.n.eq.255) then
+! print *,' CHUJ: stoi'
+! iv(21)=6
+! iv(22)=1
+! iv(23)=1
+! iv(24)=1
+! endif
+
+! min val for v(radfac) default is 0.1
+ v(24)=0.1D0
+! max val for v(radfac) default is 4.0
+ v(25)=2.0D0
+! v(25)=4.0D0
+! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+! the sumsl default is 0.1
+ v(26)=0.1D0
+! false conv if (act fnctn decrease) .lt. v(34)
+! the sumsl default is 100*machep
+ v(34)=v(34)/100.0D0
+! absolute convergence
+ if (tolf.eq.0.0D0) tolf=1.0D-4
+ v(31)=tolf
+! relative convergence
+ if (rtolf.eq.0.0D0) rtolf=1.0D-4
+ v(32)=rtolf
+! controls initial step size
+ v(35)=1.0D-1
+! large vals of d correspond to small components of step
+ do i=1,nphi
+ d(i)=1.0D-1
+ enddo
+ do i=nphi+1,nvar
+ d(i)=1.0D-1
+ enddo
+! minimize energy
+! write (iout,*) 'Processor',me,' nvar',nvar
+! write (iout,*) 'Variables BEFORE minimization:'
+! write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+
+! print *, 'MINIM_JLEE: ',me,' before SUMSL '
+
+ call func(nvar,var,nf,eee,idum,rdum,fdum)
+ nfun=nfun+1
+ if(eee.ge.1.0d20) then
+! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
+! print *,' energy before SUMSL =',eee
+! print *,' aborting local minimization'
+ iv(1)=-1
+ v(10)=eee
+ go to 201
+ endif
+
+!t time0s=MPI_WTIME()
+ call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
+!t write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
+! print *, 'MINIM_JLEE: ',me,' after SUMSL '
+
+! find which conformation was returned from sumsl
+ nfun=nfun+iv(7)
+! print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
+! & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
+! if (iv(1).ne.4 .or. nf.le.1) then
+! write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
+! write (*,*) 'Initial Variables'
+! write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
+! write (*,*) 'Variables'
+! write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+! write (*,*) 'Vector d'
+! write (*,'(8f10.4)') (d(i),i=1,nvar)
+! write (iout,*) 'Processor',me,' something bad in SUMSL',
+! & iv(1),nf
+! write (iout,*) 'Initial Variables'
+! write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
+! write (iout,*) 'Variables'
+! write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
+! write (iout,*) 'Vector d'
+! write (iout,'(8f10.4)') (d(i),i=1,nvar)
+! endif
+! if (nf.lt.iv(6)-1) then
+! recalculate intra- and interchain energies
+! call func(nvar,var,nf,v(10),iv,v,fdum)
+! else if (nf.eq.iv(6)-1) then
+! regenerate conformation
+! call var_to_geom(nvar,var)
+! call chainbuild
+! endif
+! change origin and axes to standard ECEPP format
+! call var_to_geom(nvar,var)
+! write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
+! write (iout,'(8f10.4)') (var(i),i=1,nvar)
+! write (iout,*) 'Energy:',v(10)
+! send back output
+! print *, 'MINIM_JLEE: ',me,' minimized: ',n
+ 201 continue
+ indx(1)=n
+! return code: 6-gradient 9-number of ftn evaluation, etc
+ indx(2)=iv(1)
+! total # of ftn evaluations (for iwf=0, it includes all minimizations).
+ indx(3)=nfun
+ indx(4)=info(2)
+ indx(5)=info(3)
+ indx(6)=nss
+ indx(7)=info(5)
+ indx(8)=info(6)
+ indx(9)=info(9)
+ call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,&
+ ierr)
+! send back energies
+! al & cc
+! calculate contact order
+#ifdef CO_BIAS
+ call contact(.false.,ncont,icont_,co)
+ erg(1)=v(10)-1.0d2*co
+#else
+ erg(1)=v(10)
+#endif
+ j=1
+ call mpi_send(erg,j,mpi_double_precision,king,idreal,&
+ CG_COMM,ierr)
+#ifdef CO_BIAS
+ call mpi_send(co,j,mpi_double_precision,king,idreal,&
+ CG_COMM,ierr)
+#endif
+! send back values of variables
+ call mpi_send(var,nvar,mpi_double_precision,&
+ king,idreal,CG_COMM,ierr)
+! print * , 'MINIM_JLEE: Processor',me,' send erg and var '
+
+ if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
+!d call intout
+!d call chainbuild
+!d call etotal(energia(0))
+!d etot=energia(0)
+!d call enerprint(energia(0))
+ call mpi_send(ihpb,nss,mpi_integer,&
+ king,idint,CG_COMM,ierr)
+ call mpi_send(jhpb,nss,mpi_integer,&
+ king,idint,CG_COMM,ierr)
+ endif
+
+ go to 10
+ 100 print *, ' error in receiving message from emperor', me
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 200 print *, ' error in sending message to emperor'
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 300 print *, ' error in communicating with emperor'
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ return
+ 956 format (' initial energy could not be calculated',41x)
+ 957 format (80x)
+ 965 format (' convergence code ',i2,' # of function calls ',&
+ i4,' # of gradient calls ',i4,10x)
+ 975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
+ end subroutine minim_jlee
+!-----------------------------------------------------------------------------
+! newconf.f
+!-----------------------------------------------------------------------------
+ subroutine make_var(n,idum,iter_csa)
+
+ use MD_data
+ use energy_data
+ use compare_data
+ use control_data, only: vdisulf
+ use geometry_data
+ use geometry, only: dist
+ include 'mpif.h'
+
+! use random, only: iran_num,ran_number
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.HAIRPIN'
+! include 'COMMON.VAR'
+! include 'COMMON.DISTFIT'
+! include 'COMMON.GEO'
+! include 'COMMON.CONTROL'
+ logical :: nicht_getan,nicht_getan1,fail,lfound
+ integer :: nharp,iharp(4,nres/3),nconf_harp
+ integer :: iisucc(mxio)
+ logical :: ifused(mxio)
+ integer :: nhx_seed(nseed),ihx_seed(4,nres/3,nseed) !max_seed
+ integer :: nhx_use(nseed),ihx_use(0:4,nres/3,nseed)
+ integer :: nlx_seed(nseed),ilx_seed(2,nres/3,nseed),&
+ nlx_use(nseed),ilx_use(nres/3,nseed)
+! real(kind=4) :: ran1,ran2
+
+ integer :: i,j,k,n,idum,iter_csa,iran,index,n7frag,n8frag,n14frag,&
+ n15frag,nbefrag,nlx_tot,iters,i1,i2,i3,ntot_gen,ngen,iih,&
+ ij,jr,iim,nhx_tot,idummy,iter,iif,iig,icheck,ishift,iang,&
+ n8c,ih_start,ih_end,n7c,index2,isize,nsucc,nacc,j1,nran,&
+ ierror,ierrcode
+ real(kind=8) :: d
+
+ write (iout,*) 'make_var : nseed=',nseed,'ntry=',n
+ index=0
+
+!-----------------------------------------
+ if (n7.gt.0.or.n8.gt.0.or.n9.gt.0.or.n14.gt.0.or.n15.gt.0 &
+ .or.n16.gt.0.or.n17.gt.0.or.n18.gt.0) &
+ call select_frag(n7frag,n8frag,n14frag,&
+ n15frag,nbefrag,iter_csa)
+
+!---------------------------------------------------
+! N18 - random perturbation of one phi(=gamma) angle in a loop
+!
+ IF (n18.gt.0) THEN
+ nlx_tot=0
+ do iters=1,nseed
+ i1=is(iters)
+ nlx_seed(iters)=0
+ do i2=1,n14frag
+ if (lvar_frag(i2,1).eq.i1) then
+ nlx_seed(iters)=nlx_seed(iters)+5
+ ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
+ ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
+ ilx_use(nlx_seed(iters),iters)=5
+ endif
+ enddo
+ nlx_use(iters)=nlx_seed(iters)
+ nlx_tot=nlx_tot+nlx_seed(iters)
+ enddo
+
+ if (nlx_tot .ge. n18*nseed) then
+ ntot_gen=n18*nseed
+ else
+ ntot_gen=(nlx_tot/nseed)*nseed
+ endif
+
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+ if (nlx_use(iters).gt.0) then
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nlx_seed(iters))
+ if (ilx_use(iih,iters).gt.0) then
+ nicht_getan=.false.
+ ilx_use(iih,iters)=ilx_use(iih,iters)-1
+ nlx_use(iters)=nlx_use(iters)-1
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=18
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ jr=iran_num(ilx_seed(1,iih,iters),ilx_seed(2,iih,iters))
+ d=ran_number(-pi,pi)
+ dihang_in(2,jr-2,1,index)=pinorm(dihang_in(2,jr-2,1,index)+d)
+
+
+ if (ngen.eq.ntot_gen) goto 145
+ endif
+ enddo
+ enddo
+ 145 continue
+
+ ENDIF
+
+
+!-----------------------------------------
+! N17 : zip a beta in a seed by forcing one additional p-p contact
+!
+ IF (n17.gt.0) THEN
+ nhx_tot=0
+ do iters=1,nseed
+ i1=is(iters)
+ nhx_seed(iters)=0
+ nhx_use(iters)=0
+ do i2=1,nbefrag
+ if (avar_frag(i2,1).eq.i1) then
+ nhx_seed(iters)=nhx_seed(iters)+1
+ ihx_use(2,nhx_seed(iters),iters)=1
+ if (avar_frag(i2,5)-avar_frag(i2,3).le.3.and. &
+ avar_frag(i2,2).gt.1.and.avar_frag(i2,4).lt.nres) then
+ ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+ ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
+ ihx_use(0,nhx_seed(iters),iters)=1
+ ihx_use(1,nhx_seed(iters),iters)=0
+ nhx_use(iters)=nhx_use(iters)+1
+ else
+ if (avar_frag(i2,4).gt.avar_frag(i2,5)) then
+ if (avar_frag(i2,2).gt.1.and. &
+ avar_frag(i2,4).lt.nres) then
+ ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+ ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
+ ihx_use(0,nhx_seed(iters),iters)=1
+ ihx_use(1,nhx_seed(iters),iters)=0
+ nhx_use(iters)=nhx_use(iters)+1
+ endif
+ if (avar_frag(i2,3).lt.nres.and. &
+ avar_frag(i2,5).gt.1) then
+ ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
+ ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)-1
+ ihx_use(0,nhx_seed(iters),iters)= &
+ ihx_use(0,nhx_seed(iters),iters)+1
+ ihx_use(2,nhx_seed(iters),iters)=0
+ nhx_use(iters)=nhx_use(iters)+1
+ endif
+ else
+ if (avar_frag(i2,2).gt.1.and. &
+ avar_frag(i2,4).gt.1) then
+ ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+ ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)-1
+ ihx_use(0,nhx_seed(iters),iters)=1
+ ihx_use(1,nhx_seed(iters),iters)=0
+ nhx_use(iters)=nhx_use(iters)+1
+ endif
+ if (avar_frag(i2,3).lt.nres.and. &
+ avar_frag(i2,5).lt.nres) then
+ ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
+ ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)+1
+ ihx_use(0,nhx_seed(iters),iters)= &
+ ihx_use(0,nhx_seed(iters),iters)+1
+ ihx_use(2,nhx_seed(iters),iters)=0
+ nhx_use(iters)=nhx_use(iters)+1
+ endif
+ endif
+ endif
+ endif
+ enddo
+
+ nhx_tot=nhx_tot+nhx_use(iters)
+!d write (iout,*) "debug N17",iters,nhx_seed(iters),
+!d & nhx_use(iters),nhx_tot
+ enddo
+
+ if (nhx_tot .ge. n17*nseed) then
+ ntot_gen=n17*nseed
+ else if (nhx_tot .ge. nseed) then
+ ntot_gen=(nhx_tot/nseed)*nseed
+ else
+ ntot_gen=nhx_tot
+ endif
+!d write (iout,*) "debug N17==",ntot_gen,nhx_tot,nseed
+
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+ if (nhx_use(iters).gt.0) then
+!d write (iout,*) "debug N17",nhx_use(iters),ngen,ntot_gen
+!d write (iout,*) "debugN17^",
+!d & (ihx_use(0,k,iters),k=1,nhx_use(iters))
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nhx_seed(iters))
+!d write (iout,*) "debugN17^",iih
+ if (ihx_use(0,iih,iters).gt.0) then
+ iim=iran_num(1,2)
+!d write (iout,*) "debugN17=",iih,nhx_seed(iters)
+!d write (iout,*) "debugN17-",iim,'##',
+!d & (ihx_use(k,iih,iters),k=0,2)
+!d call flush(iout)
+ do while (ihx_use(iim,iih,iters).eq.1)
+ iim=iran_num(1,2)
+!d write (iout,*) "debugN17-",iim,'##',
+!d & (ihx_use(k,iih,iters),k=0,2)
+!d call flush(iout)
+ enddo
+ nicht_getan=.false.
+ ihx_use(iim,iih,iters)=1
+ ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
+ nhx_use(iters)=nhx_use(iters)-1
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=17
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ if (iim.eq.1) then
+ idata(1,index)=ihx_seed(1,iih,iters)
+ idata(2,index)=ihx_seed(2,iih,iters)
+ else
+ idata(1,index)=ihx_seed(3,iih,iters)
+ idata(2,index)=ihx_seed(4,iih,iters)
+ endif
+
+ if (ngen.eq.ntot_gen) goto 115
+ endif
+ enddo
+ enddo
+ 115 continue
+ write (iout,*) "N17",n17," ngen/nseed",ngen/nseed,&
+ ngen,nseed
+
+
+ ENDIF
+!-----------------------------------------
+! N16 : slide non local beta in a seed by +/- 1 or +/- 2
+!
+ IF (n16.gt.0) THEN
+ nhx_tot=0
+ do iters=1,nseed
+ i1=is(iters)
+ nhx_seed(iters)=0
+ do i2=1,n7frag
+ if (bvar_frag(i2,1).eq.i1) then
+ nhx_seed(iters)=nhx_seed(iters)+1
+ ihx_seed(1,nhx_seed(iters),iters)=bvar_frag(i2,3)
+ ihx_seed(2,nhx_seed(iters),iters)=bvar_frag(i2,4)
+ ihx_seed(3,nhx_seed(iters),iters)=bvar_frag(i2,5)
+ ihx_seed(4,nhx_seed(iters),iters)=bvar_frag(i2,6)
+ ihx_use(0,nhx_seed(iters),iters)=4
+ do i3=1,4
+ ihx_use(i3,nhx_seed(iters),iters)=0
+ enddo
+ endif
+ enddo
+ nhx_use(iters)=4*nhx_seed(iters)
+ nhx_tot=nhx_tot+nhx_seed(iters)
+!d write (iout,*) "debug N16",iters,nhx_seed(iters)
+ enddo
+
+ if (4*nhx_tot .ge. n16*nseed) then
+ ntot_gen=n16*nseed
+ else if (4*nhx_tot .ge. nseed) then
+ ntot_gen=(4*nhx_tot/nseed)*nseed
+ else
+ ntot_gen=4*nhx_tot
+ endif
+ write (iout,*) "debug N16",ntot_gen,4*nhx_tot,nseed
+
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+ if (nhx_use(iters).gt.0) then
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nhx_seed(iters))
+ if (ihx_use(0,iih,iters).gt.0) then
+ iim=iran_num(1,4)
+ do while (ihx_use(iim,iih,iters).eq.1)
+!d write (iout,*) iim,
+!d & ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
+ iim=iran_num(1,4)
+ enddo
+ nicht_getan=.false.
+ ihx_use(iim,iih,iters)=1
+ ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
+ nhx_use(iters)=nhx_use(iters)-1
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=16
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do i=1,4
+ idata(i,index)=ihx_seed(i,iih,iters)
+ enddo
+ idata(5,index)=iim
+
+ if (ngen.eq.ntot_gen) goto 116
+ endif
+ enddo
+ enddo
+ 116 continue
+ write (iout,*) "N16",n16," ngen/nseed",ngen/nseed,&
+ ngen,nseed
+ ENDIF
+!-----------------------------------------
+! N15 : copy two 2nd structure elements from 1 or 2 conf. in bank to a seed
+!
+ IF (n15.gt.0) THEN
+
+ do iters=1,nseed
+ iseed=is(iters)
+ do i=1,mxio
+ ifused(i)=.false.
+ enddo
+
+ do idummy=1,n15
+ iter=0
+ 84 continue
+
+ iran=0
+ iif=iran_num(1,n15frag)
+ do while( (ifused(iif) .or. svar_frag(iif,1).eq.iseed) .and. &
+ iran.le.mxio )
+ iif=iran_num(1,n15frag)
+ iran=iran+1
+ enddo
+ if(iran.ge.mxio) goto 811
+
+ iran=0
+ iig=iran_num(1,n15frag)
+ do while( (ifused(iig) .or. svar_frag(iig,1).eq.iseed .or. &
+ .not.(svar_frag(iif,3).lt.svar_frag(iig,2).or. &
+ svar_frag(iig,3).lt.svar_frag(iif,2)) ) .and. &
+ iran.le.mxio )
+ iig=iran_num(1,n15frag)
+ iran=iran+1
+ enddo
+ if(iran.ge.mxio) goto 811
+
+ index=index+1
+ movenx(index)=15
+ parent(1,index)=iseed
+ parent(2,index)=svar_frag(iif,1)
+ parent(3,index)=svar_frag(iig,1)
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ ifused(iif)=.true.
+ ifused(iig)=.true.
+ call newconf_copy(idum,dihang_in(1,1,1,index),&
+ svar_frag(iif,1),svar_frag(iif,2),svar_frag(iif,3))
+
+ do j=svar_frag(iig,2),svar_frag(iig,3)
+ do i=1,4
+ dihang_in(i,j,1,index)=bvar(i,j,1,svar_frag(iig,1))
+ enddo
+ enddo
+
+
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) then
+ index=index-1
+ ifused(iif)=.false.
+ goto 84
+ endif
+ endif
+
+ 811 continue
+ enddo
+ enddo
+ ENDIF
+
+!-----------------------------------------
+! N14 local_move (Maurizio) for loops in a seed
+!
+ IF (n14.gt.0) THEN
+ nlx_tot=0
+ do iters=1,nseed
+ i1=is(iters)
+ nlx_seed(iters)=0
+ do i2=1,n14frag
+ if (lvar_frag(i2,1).eq.i1) then
+ nlx_seed(iters)=nlx_seed(iters)+3
+ ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
+ ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
+ ilx_use(nlx_seed(iters),iters)=3
+ endif
+ enddo
+ nlx_use(iters)=nlx_seed(iters)
+ nlx_tot=nlx_tot+nlx_seed(iters)
+!d write (iout,*) "debug N14",iters,nlx_seed(iters)
+ enddo
+
+ if (nlx_tot .ge. n14*nseed) then
+ ntot_gen=n14*nseed
+ else
+ ntot_gen=(nlx_tot/nseed)*nseed
+ endif
+!d write (iout,*) "debug N14",ntot_gen,n14frag,nseed
+
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+ if (nlx_use(iters).gt.0) then
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nlx_seed(iters))
+ if (ilx_use(iih,iters).gt.0) then
+ nicht_getan=.false.
+ ilx_use(iih,iters)=ilx_use(iih,iters)-1
+ nlx_use(iters)=nlx_use(iters)-1
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=14
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+ idata(1,index)=ilx_seed(1,iih,iters)
+ idata(2,index)=ilx_seed(2,iih,iters)
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ if (ngen.eq.ntot_gen) goto 131
+ endif
+ enddo
+ enddo
+ 131 continue
+!d write (iout,*) "N14",n14," ngen/nseed",ngen/nseed,
+!d & ngen,nseed
+
+ ENDIF
+!-----------------------------------------
+! N9 : shift a helix in a seed
+!
+ IF (n9.gt.0) THEN
+ nhx_tot=0
+ do iters=1,nseed
+ i1=is(iters)
+ nhx_seed(iters)=0
+ do i2=1,n8frag
+ if (hvar_frag(i2,1).eq.i1) then
+ nhx_seed(iters)=nhx_seed(iters)+1
+ ihx_seed(1,nhx_seed(iters),iters)=hvar_frag(i2,2)
+ ihx_seed(2,nhx_seed(iters),iters)=hvar_frag(i2,3)
+ ihx_use(0,nhx_seed(iters),iters)=4
+ do i3=1,4
+ ihx_use(i3,nhx_seed(iters),iters)=0
+ enddo
+ endif
+ enddo
+ nhx_use(iters)=4*nhx_seed(iters)
+ nhx_tot=nhx_tot+nhx_seed(iters)
+!d write (iout,*) "debug N9",iters,nhx_seed(iters)
+ enddo
+
+ if (4*nhx_tot .ge. n9*nseed) then
+ ntot_gen=n9*nseed
+ else
+ ntot_gen=(4*nhx_tot/nseed)*nseed
+ endif
+!d write (iout,*) "debug N9",ntot_gen,n8frag,nseed
+
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+ if (nhx_use(iters).gt.0) then
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nhx_seed(iters))
+ if (ihx_use(0,iih,iters).gt.0) then
+ iim=iran_num(1,4)
+ do while (ihx_use(iim,iih,iters).eq.1)
+!d write (iout,*) iim,
+!d & ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
+ iim=iran_num(1,4)
+ enddo
+ nicht_getan=.false.
+ ihx_use(iim,iih,iters)=1
+ ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
+ nhx_use(iters)=nhx_use(iters)-1
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=9
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ jstart=max(nnt,ihx_seed(1,iih,iters)+1)
+ jend=min(nct,ihx_seed(2,iih,iters))
+!d write (iout,*) "debug N9",iters,iih,jstart,jend
+ if (iim.eq.1) then
+ ishift=-2
+ else if (iim.eq.2) then
+ ishift=-1
+ else if (iim.eq.3) then
+ ishift=1
+ else if (iim.eq.4) then
+ ishift=2
+ else
+ write (iout,*) 'CHUJ NASTAPIL: iim=',iim
+#ifdef MPI !el
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+#endif
+ endif
+ do j=jstart,jend
+ if (itype(j).eq.10) then
+ iang=2
+ else
+ iang=4
+ endif
+ do i=1,iang
+ if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
+ dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
+ enddo
+ enddo
+ if (ishift.gt.0) then
+ do j=0,ishift-1
+ if (itype(jend+j).eq.10) then
+ iang=2
+ else
+ iang=4
+ endif
+ do i=1,iang
+ if (jend+j.ge.nnt.and.jend+j.le.nct) &
+ dihang_in(i,jstart+j,1,index)=bvar(i,jend+j,1,iseed)
+ enddo
+ enddo
+ else
+ do j=0,-ishift-1
+ if (itype(jstart+j).eq.10) then
+ iang=2
+ else
+ iang=4
+ endif
+ do i=1,iang
+ if (jend+j.ge.nnt.and.jend+j.le.nct) &
+ dihang_in(i,jend+j,1,index)=bvar(i,jstart+j,1,iseed)
+ enddo
+ enddo
+ endif
+ if (ngen.eq.ntot_gen) goto 133
+ endif
+ enddo
+ enddo
+ 133 continue
+!d write (iout,*) "N9",n9," ngen/nseed",ngen/nseed,
+!d & ngen,nseed
+
+ ENDIF
+!-----------------------------------------
+! N8 : copy a helix from bank to seed
+!
+ if (n8.gt.0) then
+ if (n8frag.lt.n8) then
+ write (iout,*) "N8: only ",n8frag,'helices'
+ n8c=n8frag
+ else
+ n8c=n8
+ endif
+
+ do iters=1,nseed
+ iseed=is(iters)
+ do i=1,mxio
+ ifused(i)=.false.
+ enddo
+
+
+ do idummy=1,n8c
+ iter=0
+ 94 continue
+ iran=0
+ iif=iran_num(1,n8frag)
+ do while( (ifused(iif) .or. hvar_frag(iif,1).eq.iseed) .and. &
+ iran.le.mxio )
+ iif=iran_num(1,n8frag)
+ iran=iran+1
+ enddo
+
+ if(iran.ge.mxio) goto 911
+
+ index=index+1
+ movenx(index)=8
+ parent(1,index)=iseed
+ parent(2,index)=hvar_frag(iif,1)
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ ifused(iif)=.true.
+ if (hvar_frag(iif,3)-hvar_frag(iif,2).le.6) then
+ call newconf_copy(idum,dihang_in(1,1,1,index),&
+ hvar_frag(iif,1),hvar_frag(iif,2),hvar_frag(iif,3))
+ else
+ ih_start=iran_num(hvar_frag(iif,2),hvar_frag(iif,3)-6)
+ ih_end=iran_num(ih_start,hvar_frag(iif,3))
+ call newconf_copy(idum,dihang_in(1,1,1,index),&
+ hvar_frag(iif,1),ih_start,ih_end)
+ endif
+ iter=iter+1
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) then
+ index=index-1
+ ifused(iif)=.false.
+ goto 94
+ endif
+ endif
+
+
+ 911 continue
+
+ enddo
+ enddo
+
+ endif
+
+!-----------------------------------------
+! N7 : copy nonlocal beta fragment from bank to seed
+!
+ if (n7.gt.0) then
+ if (n7frag.lt.n7) then
+ write (iout,*) "N7: only ",n7frag,'nonlocal fragments'
+ n7c=n7frag
+ else
+ n7c=n7
+ endif
+
+ do i=1,nres
+ do j=1,mxio2
+ iff_in(i,j)=0
+ enddo
+ enddo
+ index2=0
+ do i=1,mxio
+ isend2(i)=0
+ enddo
+
+ do iters=1,nseed
+ iseed=is(iters)
+ do i=1,mxio
+ ifused(i)=.false.
+ enddo
+
+ do idummy=1,n7c
+ iran=0
+ iif=iran_num(1,n7frag)
+ do while( (ifused(iif) .or. bvar_frag(iif,1).eq.iseed) .and. &
+ iran.le.mxio )
+ iif=iran_num(1,n7frag)
+ iran=iran+1
+ enddo
+
+!d write (*,'(3i5,l,4i5)'),iters,idummy,iif,ifused(iif),
+!d & bvar_frag(iif,1),iseed,iran,index2
+
+ if(iran.ge.mxio) goto 999
+ if(index2.ge.mxio2) goto 999
+
+ index=index+1
+ movenx(index)=7
+ parent(1,index)=iseed
+ parent(2,index)=bvar_frag(iif,1)
+ index2=index2+1
+ isend2(index)=index2
+ ifused(iif)=.true.
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in2(i,j,k,index2)=bvar(i,j,k,bvar_frag(iif,1))
+ enddo
+ enddo
+ enddo
+
+ if (bvar_frag(iif,2).eq.4) then
+ do i=bvar_frag(iif,3),bvar_frag(iif,4)
+ iff_in(i,index2)=1
+ enddo
+ if (bvar_frag(iif,5).lt.bvar_frag(iif,6)) then
+!d print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
+!d & bvar_frag(iif,5),bvar_frag(iif,6)
+ do i=bvar_frag(iif,5),bvar_frag(iif,6)
+ iff_in(i,index2)=1
+ enddo
+ else
+!d print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
+!d & bvar_frag(iif,6),bvar_frag(iif,5)
+ do i=bvar_frag(iif,6),bvar_frag(iif,5)
+ iff_in(i,index2)=1
+ enddo
+ endif
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+
+ 999 continue
+
+ enddo
+ enddo
+
+ endif
+!-----------------------------------------------
+! N6 : copy random continues fragment from bank to seed
+!
+ do iters=1,nseed
+ iseed=is(iters)
+ do idummy=1,n6
+ isize=(is2-is1+1)*ran1(idum)+is1
+ index=index+1
+ movenx(index)=6
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ iter=0
+ 104 continue
+ if(icycle.le.0) then
+ i1=nconf* ran1(idum)+1
+ i1=nbank-nconf+i1
+ else
+ i1=nbank* ran1(idum)+1
+ endif
+ if(i1.eq.iseed) goto 104
+ iter=iter+1
+ call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+ parent(1,index)=iseed
+ parent(2,index)=i1
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 104
+ endif
+ enddo
+ enddo
+!-----------------------------------------
+ if (n3.gt.0.or.n4.gt.0) call gen_hairpin
+ nconf_harp=0
+ do iters=1,nseed
+ if (nharp_seed(iters).gt.0) nconf_harp=nconf_harp+1
+ enddo
+!-----------------------------------------
+! N3 : copy hairpin from bank to seed
+!
+ do iters=1,nseed
+ iseed=is(iters)
+ nsucc=0
+ nacc=0
+ do idummy=1,n3
+ index=index+1
+ iter=0
+ 124 continue
+ if(icycle.le.0) then
+ i1=nconf* ran1(idum)+1
+ i1=nbank-nconf+i1
+ else
+ i1=nbank* ran1(idum)+1
+ endif
+ if(i1.eq.iseed) goto 124
+ do k=1,nsucc
+ if (i1.eq.iisucc(k).and.nsucc.lt.nconf_harp-1) goto 124
+ enddo
+ nsucc=nsucc+1
+ iisucc(nsucc)=i1
+ iter=iter+1
+ call newconf_residue_hairpin(idum,dihang_in(1,1,1,index),&
+ i1,fail)
+ if (fail) then
+ if (icycle.le.0 .and. nsucc.eq.nconf .or. &
+ icycle.gt.0 .and. nsucc.eq.nbank) then
+ index=index-1
+ goto 125
+ else
+ goto 124
+ endif
+ endif
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 124
+ endif
+ movenx(index)=3
+ parent(1,index)=iseed
+ parent(2,index)=i1
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ nacc=nacc+1
+ enddo
+! if not enough hairpins, supplement with windows
+ 125 continue
+!dd if (n3.ne.0) write (iout,*) "N3",n3," nsucc",nsucc," nacc",nacc
+ do idummy=nacc+1,n3
+ isize=(is2-is1+1)*ran1(idum)+is1
+ index=index+1
+ movenx(index)=6
+ parent(1,index)=iseed
+ parent(2,index)=i1
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ iter=0
+ 114 continue
+ if(icycle.le.0) then
+ i1=nconf* ran1(idum)+1
+ i1=nbank-nconf+i1
+ else
+ i1=nbank* ran1(idum)+1
+ endif
+ if(i1.eq.iseed) goto 114
+ iter=iter+1
+ call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 114
+ endif
+ enddo
+ enddo
+!-----------------------------------------
+! N4 : shift a turn in hairpin in seed
+!
+ IF (N4.GT.0) THEN
+ if (4*nharp_tot .ge. n4*nseed) then
+ ntot_gen=n4*nseed
+ else
+ ntot_gen=(4*nharp_tot/nseed)*nseed
+ endif
+ ngen=0
+ do while (ngen.lt.ntot_gen)
+ do iters=1,nseed
+ iseed=is(iters)
+! write (iout,*) 'iters',iters,' iseed',iseed,' nharp_seed',
+! & nharp_seed(iters),' nharp_use',nharp_use(iters),
+! & ' ntot_gen',ntot_gen
+! write (iout,*) 'iharp_use(0)',
+! & (iharp_use(0,k,iters),k=1,nharp_seed(iters))
+ if (nharp_use(iters).gt.0) then
+ nicht_getan=.true.
+ do while (nicht_getan)
+ iih=iran_num(1,nharp_seed(iters))
+! write (iout,*) 'iih',iih,' iharp_use',
+! & (iharp_use(k,iih,iters),k=1,4)
+ if (iharp_use(0,iih,iters).gt.0) then
+ nicht_getan1=.true.
+ do while (nicht_getan1)
+ iim=iran_num(1,4)
+ nicht_getan1=iharp_use(iim,iih,iters).eq.1
+ enddo
+ nicht_getan=.false.
+ iharp_use(iim,iih,iters)=1
+ iharp_use(0,iih,iters)=iharp_use(0,iih,iters)-1
+ nharp_use(iters)=nharp_use(iters)-1
+!dd write (iout,'(a16,i3,a5,i2,a10,2i4)')
+!dd & 'N4 selected hairpin',iih,' move',iim,' iharp_seed',
+!dd & iharp_seed(1,iih,iters),iharp_seed(2,iih,iters)
+ endif
+ enddo
+ ngen=ngen+1
+ index=index+1
+ movenx(index)=4
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+ jstart=iharp_seed(1,iih,iters)+1
+ jend=iharp_seed(2,iih,iters)
+ if (iim.eq.1) then
+ ishift=-2
+ else if (iim.eq.2) then
+ ishift=-1
+ else if (iim.eq.3) then
+ ishift=1
+ else if (iim.eq.4) then
+ ishift=2
+ else
+ write (iout,*) 'CHUJ NASTAPIL: iim=',iim
+#ifdef MPI !el
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+#endif !el
+ endif
+! write (iout,*) 'jstart',jstart,' jend',jend,' ishift',ishift
+! write (iout,*) 'Before turn shift'
+! do j=2,nres-1
+! theta(j+1)=dihang_in(1,j,1,index)
+! phi(j+2)=dihang_in(2,j,1,index)
+! alph(j)=dihang_in(3,j,1,index)
+! omeg(j)=dihang_in(4,j,1,index)
+! enddo
+! call intout
+ do j=jstart,jend
+ if (itype(j).eq.10) then
+ iang=2
+ else
+ iang=4
+ endif
+ do i=1,iang
+ if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
+ dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
+ enddo
+ enddo
+! write (iout,*) 'After turn shift'
+! do j=2,nres-1
+! theta(j+1)=dihang_in(1,j,1,index)
+! phi(j+2)=dihang_in(2,j,1,index)
+! alph(j)=dihang_in(3,j,1,index)
+! omeg(j)=dihang_in(4,j,1,index)
+! enddo
+! call intout
+ if (ngen.eq.ntot_gen) goto 135
+ endif
+ enddo
+ enddo
+! if not enough hairpins, supplement with windows
+! write (iout,*) 'end of enddo'
+ 135 continue
+!dd write (iout,*) "N4",n4," ngen/nseed",ngen/nseed,
+!dd & ngen,nseed
+ do iters=1,nseed
+ iseed=is(iters)
+ do idummy=ngen/nseed+1,n4
+ isize=(is2-is1+1)*ran1(idum)+is1
+ index=index+1
+ movenx(index)=6
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+
+ iter=0
+ 134 continue
+ if(icycle.le.0) then
+ i1=nconf* ran1(idum)+1
+ i1=nbank-nconf+i1
+ else
+ i1=nbank* ran1(idum)+1
+ endif
+ if(i1.eq.iseed) goto 134
+ iter=iter+1
+ call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+ parent(1,index)=iseed
+ parent(2,index)=i1
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 134
+ endif
+ enddo
+ enddo
+ ENDIF
+!-----------------------------------------
+! N5 : copy one residue from bank to seed (normally switched off - use N1)
+!
+ do iters=1,nseed
+ iseed=is(iters)
+ isize=1
+ do i=1,n5
+ index=index+1
+ movenx(index)=5
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+
+ iter=0
+ 105 continue
+ if(icycle.le.0) then
+ i1=nconf* ran1(idum)+1
+ i1=nbank-nconf+i1
+ else
+ i1=nbank* ran1(idum)+1
+ endif
+ if(i1.eq.iseed) goto 105
+ iter=iter+1
+ call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+ parent(1,index)=iseed
+ parent(2,index)=i1
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 105
+ endif
+ enddo
+ enddo
+!-----------------------------------------
+! N2 : copy backbone of one residue from bank or first bank to seed
+! (normally switched off - use N1)
+!
+ do iters=1,nseed
+ iseed=is(iters)
+ do i=n2,1,-1
+ if(icycle.le.0.and.iuse.gt.nconf-irr) then
+ iseed=ran1(idum)*nconf+1
+ iseed=nbank-nconf+iseed
+ endif
+ index=index+1
+ movenx(index)=2
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ iter=0
+ 102 i1= ran1(idum)*nbank+1
+ if(i1.eq.iseed) goto 102
+ iter=iter+1
+ if(icycle.le.0.and.iuse.gt.nconf-irr) then
+ nran=mod(i-1,nran0)+3
+ call newconf1arr(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=-iseed
+ parent(2,index)=-i1
+ else if(icycle.le.0.and.iters.le.iuse) then
+ nran=mod(i-1,nran0)+1
+ call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=-i1
+ else
+ nran=mod(i-1,nran1)+1
+ if(ran1(idum).lt.0.5) then
+ call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=-i1
+ else
+ call newconf1abb(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=i1
+ endif
+ endif
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 102
+ endif
+ enddo
+ enddo
+!-----------------------------------------
+! N1 : copy backbone or sidechain of one residue from bank or
+! first bank to seed
+!
+ do iters=1,nseed
+ iseed=is(iters)
+ do i=n1,1,-1
+ if(icycle.le.0.and.iuse.gt.nconf-irr) then
+ iseed=ran1(idum)*nconf+1
+ iseed=nbank-nconf+iseed
+ endif
+ index=index+1
+ movenx(index)=1
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ iter=0
+ 101 i1= ran1(idum)*nbank+1
+
+ if(i1.eq.iseed) goto 101
+ iter=iter+1
+ if(icycle.le.0.and.iuse.gt.nconf-irr) then
+ nran=mod(i-1,nran0)+3
+ call newconf1rr(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=-iseed
+ parent(2,index)=-i1
+ else if(icycle.le.0.and.iters.le.iuse) then
+ nran=mod(i-1,nran0)+1
+ call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=-i1
+ else
+ nran=mod(i-1,nran1)+1
+ if(ran1(idum).lt.0.5) then
+ call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=-i1
+ else
+ call newconf1bb(idum,dihang_in(1,1,1,index),nran,i1)
+ parent(1,index)=iseed
+ parent(2,index)=i1
+ endif
+ endif
+ if(iter.lt.10) then
+ call check_old(icheck,index)
+ if(icheck.eq.1) goto 101
+ endif
+ enddo
+ enddo
+!-----------------------------------------
+! N0 just all seeds
+!
+ IF (n0.gt.0) THEN
+ do iters=1,nseed
+ iseed=is(iters)
+ index=index+1
+ movenx(index)=0
+ parent(1,index)=iseed
+ parent(2,index)=0
+
+ if (vdisulf) then
+ nss_in(index)=bvar_nss(iseed)
+ do ij=1,nss_in(index)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ endif
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+ enddo
+ ENDIF
+!-----------------------------------------
+ if (vdisulf) then
+ do iters=1,nseed
+ iseed=is(iters)
+
+ do k=1,numch
+ do j=2,nres-1
+ theta(j+1)=bvar(1,j,k,iseed)
+ phi(j+2)=bvar(2,j,k,iseed)
+ alph(j)=bvar(3,j,k,iseed)
+ omeg(j)=bvar(4,j,k,iseed)
+ enddo
+ enddo
+ call chainbuild
+
+!d write(iout,*) 'makevar DYNSS',iseed,'#',bvar_ns(iseed),
+!d & (bvar_s(k,iseed),k=1,bvar_ns(iseed)),
+!d & bvar_nss(iseed),
+!d & (bvar_ss(1,k,iseed)-nres,'-',
+!d & bvar_ss(2,k,iseed)-nres,k=1,bvar_nss(iseed))
+
+ do i1=1,bvar_ns(iseed)
+!
+! N10 fussion of free halfcysteines in seed
+! first select CYS with distance < 7A
+!
+ do j1=i1+1,bvar_ns(iseed)
+ if (dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres) &
+ .lt.7.0.and. &
+ iabs(bvar_s(i1,iseed)-bvar_s(j1,iseed)).gt.3) then
+
+ index=index+1
+ movenx(index)=10
+ parent(1,index)=iseed
+ parent(2,index)=0
+ do ij=1,bvar_nss(iseed)
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ enddo
+ ij=bvar_nss(iseed)+1
+ nss_in(index)=ij
+ iss_in(ij,index)=bvar_s(i1,iseed)+nres
+ jss_in(ij,index)=bvar_s(j1,iseed)+nres
+
+!d write(iout,*) 'makevar NSS0',index,
+!d & dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres),
+!d & nss_in(index),iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ endif
+ enddo
+!
+! N11 type I transdisulfidation
+!
+ do j1=1,bvar_nss(iseed)
+ if (dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)) &
+ .lt.7.0.and. &
+ iabs(bvar_s(i1,iseed)-(bvar_ss(1,j1,iseed)-nres)) &
+ .gt.3) then
+
+ index=index+1
+ movenx(index)=11
+ parent(1,index)=iseed
+ parent(2,index)=0
+ do ij=1,bvar_nss(iseed)
+ if (ij.ne.j1) then
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ endif
+ enddo
+ nss_in(index)=bvar_nss(iseed)
+ iss_in(j1,index)=bvar_s(i1,iseed)+nres
+ jss_in(j1,index)=bvar_ss(1,j1,iseed)
+ if (iss_in(j1,index).gt.jss_in(j1,index)) then
+ iss_in(j1,index)=bvar_ss(1,j1,iseed)
+ jss_in(j1,index)=bvar_s(i1,iseed)+nres
+ endif
+
+!d write(iout,*) 'makevar NSS1 #1',index,
+!d & bvar_s(i1,iseed),bvar_ss(1,j1,iseed)-nres,
+!d & dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)),
+!d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+!d & ij=1,nss_in(index))
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+ endif
+ if (dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)) &
+ .lt.7.0.and. &
+ iabs(bvar_s(i1,iseed)-(bvar_ss(2,j1,iseed)-nres)) &
+ .gt.3) then
+
+ index=index+1
+ movenx(index)=11
+ parent(1,index)=iseed
+ parent(2,index)=0
+ do ij=1,bvar_nss(iseed)
+ if (ij.ne.j1) then
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ endif
+ enddo
+ nss_in(index)=bvar_nss(iseed)
+ iss_in(j1,index)=bvar_s(i1,iseed)+nres
+ jss_in(j1,index)=bvar_ss(2,j1,iseed)
+ if (iss_in(j1,index).gt.jss_in(j1,index)) then
+ iss_in(j1,index)=bvar_ss(2,j1,iseed)
+ jss_in(j1,index)=bvar_s(i1,iseed)+nres
+ endif
+
+
+!d write(iout,*) 'makevar NSS1 #2',index,
+!d & bvar_s(i1,iseed),bvar_ss(2,j1,iseed)-nres,
+!d & dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)),
+!d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+!d & ij=1,nss_in(index))
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ endif
+ enddo
+ enddo
+
+!
+! N12 type II transdisulfidation
+!
+ do i1=1,bvar_nss(iseed)
+ do j1=i1+1,bvar_nss(iseed)
+ if (dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)) &
+ .lt.7.0.and. &
+ dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)) &
+ .lt.7.0.and. &
+ iabs(bvar_ss(1,i1,iseed)-bvar_ss(1,j1,iseed)) &
+ .gt.3.and. &
+ iabs(bvar_ss(2,i1,iseed)-bvar_ss(2,j1,iseed)) &
+ .gt.3) then
+ index=index+1
+ movenx(index)=12
+ parent(1,index)=iseed
+ parent(2,index)=0
+ do ij=1,bvar_nss(iseed)
+ if (ij.ne.i1 .and. ij.ne.j1) then
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ endif
+ enddo
+ nss_in(index)=bvar_nss(iseed)
+ iss_in(i1,index)=bvar_ss(1,i1,iseed)
+ jss_in(i1,index)=bvar_ss(1,j1,iseed)
+ if (iss_in(i1,index).gt.jss_in(i1,index)) then
+ iss_in(i1,index)=bvar_ss(1,j1,iseed)
+ jss_in(i1,index)=bvar_ss(1,i1,iseed)
+ endif
+ iss_in(j1,index)=bvar_ss(2,i1,iseed)
+ jss_in(j1,index)=bvar_ss(2,j1,iseed)
+ if (iss_in(j1,index).gt.jss_in(j1,index)) then
+ iss_in(j1,index)=bvar_ss(2,j1,iseed)
+ jss_in(j1,index)=bvar_ss(2,i1,iseed)
+ endif
+
+
+!d write(iout,*) 'makevar NSS2 #1',index,
+!d & bvar_ss(1,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
+!d & dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)),
+!d & bvar_ss(2,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
+!d & dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)),
+!d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+!d & ij=1,nss_in(index))
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ endif
+
+ if (dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)) &
+ .lt.7.0.and. &
+ dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)) &
+ .lt.7.0.and. &
+ iabs(bvar_ss(1,i1,iseed)-bvar_ss(2,j1,iseed)) &
+ .gt.3.and. &
+ iabs(bvar_ss(2,i1,iseed)-bvar_ss(1,j1,iseed)) &
+ .gt.3) then
+ index=index+1
+ movenx(index)=12
+ parent(1,index)=iseed
+ parent(2,index)=0
+ do ij=1,bvar_nss(iseed)
+ if (ij.ne.i1 .and. ij.ne.j1) then
+ iss_in(ij,index)=bvar_ss(1,ij,iseed)
+ jss_in(ij,index)=bvar_ss(2,ij,iseed)
+ endif
+ enddo
+ nss_in(index)=bvar_nss(iseed)
+ iss_in(i1,index)=bvar_ss(1,i1,iseed)
+ jss_in(i1,index)=bvar_ss(2,j1,iseed)
+ if (iss_in(i1,index).gt.jss_in(i1,index)) then
+ iss_in(i1,index)=bvar_ss(2,j1,iseed)
+ jss_in(i1,index)=bvar_ss(1,i1,iseed)
+ endif
+ iss_in(j1,index)=bvar_ss(2,i1,iseed)
+ jss_in(j1,index)=bvar_ss(1,j1,iseed)
+ if (iss_in(j1,index).gt.jss_in(j1,index)) then
+ iss_in(j1,index)=bvar_ss(1,j1,iseed)
+ jss_in(j1,index)=bvar_ss(2,i1,iseed)
+ endif
+
+
+!d write(iout,*) 'makevar NSS2 #2',index,
+!d & bvar_ss(1,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
+!d & dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)),
+!d & bvar_ss(2,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
+!d & dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)),
+!d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+!d & ij=1,nss_in(index))
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ endif
+
+
+ enddo
+ enddo
+!
+! N13 removal of disulfide bond
+!
+ if (bvar_nss(iseed).gt.0) then
+ i1=bvar_nss(iseed)*ran1(idum)+1
+
+ index=index+1
+ movenx(index)=13
+ parent(1,index)=iseed
+ parent(2,index)=0
+ ij=0
+ do j1=1,bvar_nss(iseed)
+ if (j1.ne.i1) then
+ ij=ij+1
+ iss_in(ij,index)=bvar_ss(1,j1,iseed)
+ jss_in(ij,index)=bvar_ss(2,j1,iseed)
+ endif
+ enddo
+ nss_in(index)=bvar_nss(iseed)-1
+
+!d write(iout,*) 'NSS3',index,i1,
+!d & bvar_ss(1,i1,iseed)-nres,'=',bvar_ss(2,i1,iseed)-nres,'#',
+!d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+!d & ij=1,nss_in(index))
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ endif
+
+ enddo
+ endif
+!-----------------------------------------
+
+
+
+ if(index.ne.n) write(iout,*)'make_var : ntry=',index
+
+ n=index
+!d do ii=1,n
+!d write (istat,*) "======== ii=",ii," the dihang array"
+!d do i=1,nres
+!d write (istat,'(i5,4f15.5)') i,(dihang_in(k,i,1,ii)*rad2deg,k=1,4)
+!d enddo
+!d enddo
+ return
+ end subroutine make_var
+!-----------------------------------------------------------------------------
+ subroutine check_old(icheck,n)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+ integer :: icheck,n,i1,i2,m,j,i
+ real(kind=8) :: ctdif,ctdiff,diff,dif
+
+ data ctdif /10./
+ data ctdiff /60./
+
+ i1=n
+ do i2=1,n-1
+ diff=0.d0
+ do m=1,numch
+ do j=2,nres-1
+ do i=1,4
+ dif=rad2deg*dabs(dihang_in(i,j,m,i1)-dihang_in(i,j,m,i2))
+ if(dif.gt.180.0) dif=360.0-dif
+ if(dif.gt.ctdif) goto 100
+ diff=diff+dif
+ if(diff.gt.ctdiff) goto 100
+ enddo
+ enddo
+ enddo
+ icheck=1
+ return
+ 100 continue
+ enddo
+
+ icheck=0
+
+ return
+ end subroutine check_old
+!-----------------------------------------------------------------------------
+ subroutine newconf1rr(idum,vvar,nran,i1)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
+ real(kind=8) :: ctdif,dif
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=rvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=ntotgr
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=rvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1rr
+!-----------------------------------------------------------------------------
+ subroutine newconf1br(idum,vvar,nran,i1)
+
+ use energy_data, only: ndih_nconstr,idih_nconstr
+ use control_data, only: i2ndstr
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,juhc,ind
+ real(kind=8) :: ctdif,dif,rtmp
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=ntotgr
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(i2ndstr.gt.0) then
+ rtmp=ran1(idum)
+ if(rtmp.le.rdih_bias) then
+ i=0
+ do j=1,ndih_nconstr
+ if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
+ enddo
+ if(i.eq.0) then
+ juhc=0
+4321 juhc=juhc+1
+ iran= ran1(idum)*number+1
+ i=0
+ do j=1,ndih_nconstr
+ if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
+ enddo
+ if(i.eq.0.or.juhc.lt.1000)goto 4321
+ if(juhc.eq.1000) then
+ print *, 'move 6 : failed to find unconstrained group'
+ write(iout,*) 'move 6 : failed to find unconstrained group'
+ endif
+ endif
+ endif
+ endif
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=rvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1br
+!-----------------------------------------------------------------------------
+ subroutine newconf1bb(idum,vvar,nran,i1)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
+ real(kind=8) :: ctdif,dif
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=ntotgr
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=bvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1bb
+!-----------------------------------------------------------------------------
+ subroutine newconf1arr(idum,vvar,nran,i1)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
+ real(kind=8) :: ctdif,dif
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=rvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=nres-2
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=rvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1arr
+!-----------------------------------------------------------------------------
+ subroutine newconf1abr(idum,vvar,nran,i1)
+
+ use energy_data, only: ndih_nconstr,idih_nconstr
+ use control_data, only: i2ndstr
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
+ real(kind=8) :: ctdif,dif,rtmp
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=nres-2
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(i2ndstr.gt.0) then
+ rtmp=ran1(idum)
+ if(rtmp.le.rdih_bias) then
+ iran=ran1(idum)*ndih_nconstr+1
+ iran=idih_nconstr(iran)
+ endif
+ endif
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=rvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1abr
+!-----------------------------------------------------------------------------
+ subroutine newconf1abb(idum,vvar,nran,i1)
+
+ use energy_data, only: ndih_nconstr,idih_nconstr
+ use control_data, only: i2ndstr
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
+ real(kind=8) :: ctdif,dif,rtmp
+
+ ctdif=10.
+
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+ do index=1,nran
+ iold(index) = 0
+ enddo
+
+ number=nres-2
+
+ iter=0
+ do index=1,nran
+ 10 iran= ran1(idum)*number+1
+ if(i2ndstr.gt.0) then
+ rtmp=ran1(idum)
+ if(rtmp.le.rdih_bias) then
+ iran=ran1(idum)*ndih_nconstr+1
+ iran=idih_nconstr(iran)
+ endif
+ endif
+ if(iter.gt.number) return
+ iter=iter+1
+ if(iter.eq.1) goto 11
+ do ind=1,index-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ if(iter.gt.number) goto 20
+ goto 10
+ 20 continue
+ do ind=1,ngroup(iran)
+ i=igroup(1,ind,iran)
+ j=igroup(2,ind,iran)
+ k=igroup(3,ind,iran)
+ vvar(i,j,k)=bvar(i,j,k,i1)
+ enddo
+ iold(index)=iran
+ enddo
+
+ return
+ end subroutine newconf1abb
+!-----------------------------------------------------------------------------
+ subroutine newconf_residue(idum,vvar,i1,isize)
+
+ use energy_data, only: ndih_nconstr,idih_nconstr
+ use control_data, only: i2ndstr
+ use MPI_data
+ include 'mpif.h'
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,i1,isize,iran,number,iter,ind,iend,istart,&
+ ierror,ierrcode
+ real(kind=8) :: ctdif,dif,rtmp
+
+ ctdif=10.
+
+ if (iseed.gt.mxio .or. iseed.lt.1) then
+ write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+
+ k=1
+ number=nres+isize-2
+ iter=1
+ 10 iran= ran1(idum)*number+1
+ if(i2ndstr.gt.0) then
+ rtmp=ran1(idum)
+ if(rtmp.le.rdih_bias) then
+ iran=ran1(idum)*ndih_nconstr+1
+ iran=idih_nconstr(iran)
+ endif
+ endif
+ istart=iran-isize+1
+ iend=iran
+ if(istart.lt.2) istart=2
+ if(iend.gt.nres-1) iend=nres-1
+
+ if(iter.eq.1) goto 11
+ do ind=1,iter-1
+ if(iran.eq.iold(ind)) goto 10
+ enddo
+ 11 continue
+
+ do j=istart,iend
+ do i=1,4
+ dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ enddo
+ iold(iter)=iran
+ iter=iter+1
+ if(iter.gt.number) goto 20
+ goto 10
+
+ 20 continue
+ do j=istart,iend
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,i1)
+ enddo
+ enddo
+
+ return
+ end subroutine newconf_residue
+!-----------------------------------------------------------------------------
+ subroutine newconf_copy(idum,vvar,i1,istart,iend)
+
+ use MPI_data
+ include 'mpif.h'
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: i,j,k,idum,i1,istart,iend,ierror,ierrcode
+ real(kind=8) :: ctdif,dif
+
+ ctdif=10.
+
+ if (iseed.gt.mxio .or. iseed.lt.1) then
+ write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+
+
+ do j=istart,iend
+ do i=1,4
+ vvar(i,j,1)=bvar(i,j,1,i1)
+ enddo
+ enddo
+
+ return
+ end subroutine newconf_copy
+!-----------------------------------------------------------------------------
+ subroutine newconf_residue_hairpin(idum,vvar,i1,fail)
+
+ use geometry_data
+! use random, only: iran_num
+ use MPI_data
+ use compare, only:hairpin
+
+ include 'mpif.h'
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! real(kind=4) :: ran1,ran2
+ real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch)
+ integer,dimension(ntotal) :: iold
+ integer :: nharp,iharp(4,nres/3),icipa(nres/3)
+ logical :: fail,not_done
+ integer :: idum,i,j,k,i1,iend,istart,iii,n_used,icount,iih,&
+ ierror,ierrcode
+ real(kind=8) :: ctdif,dif
+
+ ctdif=10.
+
+ fail=.false.
+ if (iseed.gt.mxio .or. iseed.lt.1) then
+ write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+ do k=1,numch
+ do j=2,nres-1
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,iseed)
+ enddo
+ enddo
+ enddo
+ do k=1,numch
+ do j=2,nres-1
+ theta(j+1)=bvar(1,j,k,i1)
+ phi(j+2)=bvar(2,j,k,i1)
+ alph(j)=bvar(3,j,k,i1)
+ omeg(j)=bvar(4,j,k,i1)
+ enddo
+ enddo
+! call intout
+ call chainbuild
+ call hairpin(.false.,nharp,iharp)
+
+ if (nharp.eq.0) then
+ fail=.true.
+ return
+ endif
+
+ n_used=0
+
+ DO III=1,NHARP
+
+ not_done = .true.
+ icount=0
+ do while (not_done)
+ icount=icount+1
+ iih=iran_num(1,nharp)
+ do k=1,n_used
+ if (iih.eq.icipa(k)) then
+ iih=0
+ goto 22
+ endif
+ enddo
+ not_done=.false.
+ n_used=n_used+1
+ icipa(n_used)=iih
+ 22 continue
+ not_done = not_done .and. icount.le.nharp
+ enddo
+
+ if (iih.eq.0) then
+ write (iout,*) "CHUJ NASTAPIL W NEWCONF_RESIDUE_HAIRPIN!!!!"
+ fail=.true.
+ return
+ endif
+
+ istart=iharp(1,iih)+1
+ iend=iharp(2,iih)
+
+!dd write (iout,*) "newconf_residue_hairpin: iih",iih,
+!dd & " istart",istart," iend",iend
+
+ do k=1,numch
+ do j=istart,iend
+ do i=1,4
+ dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+ if(dif.gt.180.) dif=360.-dif
+ if(dif.gt.ctdif) goto 20
+ enddo
+ enddo
+ enddo
+ goto 10
+ 20 continue
+ do k=1,numch
+ do j=istart,iend
+ do i=1,4
+ vvar(i,j,k)=bvar(i,j,k,i1)
+ enddo
+ enddo
+ enddo
+! do j=1,numch
+! do l=2,nres-1
+! write (iout,'(4f8.3)') (rad2deg*vvar(i,l,j),i=1,4)
+! enddo
+! enddo
+ return
+ 10 continue
+ ENDDO
+
+ fail=.true.
+
+ return
+ end subroutine newconf_residue_hairpin
+!-----------------------------------------------------------------------------
+ subroutine gen_hairpin
+
+ use geometry_data
+ use MD_data
+ use compare, only:hairpin
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.HAIRPIN'
+ integer :: i1,j,k,iters
+
+! write (iout,*) 'Entering GEN_HAIRPIN'
+ do iters=1,nseed
+ i1=is(iters)
+ do k=1,numch
+ do j=2,nres-1
+ theta(j+1)=bvar(1,j,k,i1)
+ phi(j+2)=bvar(2,j,k,i1)
+ alph(j)=bvar(3,j,k,i1)
+ omeg(j)=bvar(4,j,k,i1)
+ enddo
+ enddo
+ call chainbuild
+ call hairpin(.false.,nharp_seed(iters),iharp_seed(1,1,iters))
+ enddo
+
+ nharp_tot=0
+ do iters=1,nseed
+ nharp_tot=nharp_tot+nharp_seed(iters)
+ nharp_use(iters)=4*nharp_seed(iters)
+ do j=1,nharp_seed(iters)
+ iharp_use(0,j,iters)=4
+ do k=1,4
+ iharp_use(k,j,iters)=0
+ enddo
+ enddo
+ enddo
+
+ write (iout,*) 'GEN_HAIRPIN: nharp_tot',nharp_tot
+!dd do i=1,nseed
+!dd write (iout,*) 'seed',i
+!dd write (iout,*) 'nharp_seed',nharp_seed(i),
+!dd & ' nharp_use',nharp_use(i)
+!d write (iout,*) 'iharp_seed, iharp_use'
+!d do j=1,nharp_seed(i)
+!d write (iout,'(7i3)') iharp_seed(1,j,i),iharp_seed(2,j,i),
+!d & (iharp_use(k,j,i),k=0,4)
+!d enddo
+!dd enddo
+ return
+ end subroutine gen_hairpin
+!-----------------------------------------------------------------------------
+ subroutine select_frag(nn,nh,nl,ns,nb,i_csa)
+
+ use geometry_data
+ use MD_data
+ use compare_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.HAIRPIN'
+! include 'COMMON.DISTFIT'
+ character(len=50) :: linia
+ integer :: isec(nres)
+ integer :: i,j,i1,k,nn,nh,nl,ns,nb,i_csa,nl1,ns1
+
+ nn=0
+ nh=0
+ nl=0
+ ns=0
+ nb=0
+!d write (iout,*) 'Entering select_frag'
+ do i1=1,nbank
+ do i=1,nres
+ isec(i)=0
+ enddo
+ do k=1,numch
+ do j=2,nres-1
+ theta(j+1)=bvar(1,j,k,i1)
+ phi(j+2)=bvar(2,j,k,i1)
+ alph(j)=bvar(3,j,k,i1)
+ omeg(j)=bvar(4,j,k,i1)
+ enddo
+ enddo
+ call chainbuild
+!d write (iout,*) ' -- ',i1,' -- '
+ call secondary2(.false.)
+!
+! bvar_frag nn==pair of nonlocal strands in beta sheet (loop>4)
+! strands > 4 residues; used by N7 and N16
+!
+ do j=1,nbfrag
+!
+!test 09/12/02 bfrag(2,j)-bfrag(1,j).gt.3
+!
+ 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
+
+ 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.4 ) then
+ nn=nn+1
+
+
+ if (bfrag(3,j).lt.bfrag(4,j)) then
+ write(linia,'(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(linia,'(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
+!d call write_pdb(i_csa*1000+nn+nh,linia,0d0)
+
+ bvar_frag(nn,1)=i1
+ bvar_frag(nn,2)=4
+ do i=1,4
+ bvar_frag(nn,i+2)=bfrag(i,j)
+ enddo
+ endif
+ enddo
+
+!
+! hvar_frag nh==helices; used by N8 and N9
+!
+ do j=1,nhfrag
+
+ do i=hfrag(1,j),hfrag(2,j)
+ isec(i)=2
+ enddo
+
+ if ( hfrag(2,j)-hfrag(1,j).gt.4 ) then
+ nh=nh+1
+
+!d write(linia,'(a6,i3,a1,i3)')
+!d & "select",hfrag(1,j)-1,"-",hfrag(2,j)-1
+!d call write_pdb(i_csa*1000+nn+nh,linia,0d0)
+
+ hvar_frag(nh,1)=i1
+ hvar_frag(nh,2)=hfrag(1,j)
+ hvar_frag(nh,3)=hfrag(2,j)
+ endif
+ enddo
+
+
+!v write(iout,'(i4,1pe12.4,1x,1000i1)')
+!v & i1,bene(i1),(isec(i),i=1,nres)
+!v write(linia,'(i4,1x,1000i1)')
+!v & i1,(isec(i),i=1,nres)
+!v call write_pdb(i_csa*1000+i1,linia,bene(i1))
+!
+! lvar_frag nl==loops; used by N14
+!
+ i=1
+ nl1=nl
+ do while (i.lt.nres)
+ if (isec(i).eq.0) then
+ nl=nl+1
+ lvar_frag(nl,1)=i1
+ lvar_frag(nl,2)=i
+ i=i+1
+ do while (isec(i).eq.0.and.i.le.nres)
+ i=i+1
+ enddo
+ lvar_frag(nl,3)=i-1
+ if (lvar_frag(nl,3)-lvar_frag(nl,2).lt.1) nl=nl-1
+ endif
+ i=i+1
+ enddo
+!d write(iout,'(4i5)') (i,(lvar_frag(i,ii),ii=1,3),i=nl1+1,nl)
+
+!
+! svar_frag ns==an secondary structure element; used by N15
+!
+ i=1
+ ns1=ns
+ do while (i.lt.nres)
+ if (isec(i).gt.0) then
+ ns=ns+1
+ svar_frag(ns,1)=i1
+ svar_frag(ns,2)=i
+ i=i+1
+ do while (isec(i).gt.0.and.isec(i-1).eq.isec(i) &
+ .and.i.le.nres)
+ i=i+1
+ enddo
+ svar_frag(ns,3)=i-1
+ if (svar_frag(ns,3)-svar_frag(ns,2).lt.1) ns=ns-1
+ endif
+ if (isec(i).eq.0) i=i+1
+ enddo
+!d write(iout,'(4i5)') (i,(svar_frag(i,ii),ii=1,3),i=ns1+1,ns)
+
+!
+! avar_frag nb==any pair of beta strands; used by N17
+!
+ do j=1,nbfrag
+ nb=nb+1
+ avar_frag(nb,1)=i1
+ do i=1,4
+ avar_frag(nb,i+1)=bfrag(i,j)
+ enddo
+ enddo
+
+ enddo
+
+ return
+ end subroutine select_frag
+!-----------------------------------------------------------------------------
+! together.F
+!-----------------------------------------------------------------------------
+ subroutine together
+
+! feeds tasks for parallel processing
+ use MPI_data
+ use geometry_data
+ use control_data, only: vdisulf
+ use energy_data
+ use io, only:from_int,write_csa_pdb
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! real(kind=4) :: 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(kind=4) :: tcpu
+ real(kind=8) :: time_start,time_start_c,time0f,time0i
+ logical :: ovrtim,sync_iter,timeout,flag,timeout1
+ integer,dimension(mpi_status_size) :: muster
+ real(kind=8),dimension(0:100) :: t100
+ integer,dimension(mxio) :: indx
+ real(kind=8),dimension(6*nres) :: xout !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
+ integer,dimension(9) :: ind
+ real(kind=8),dimension(2) :: cout
+ real(kind=8),parameter :: rad=1.745329252d-2
+
+ integer :: i,m,j,jlee,nft,idum,nrmsdb,nrmsdb1,ierr,ierror,ierrcode,&
+ ntrial,ntry,idum2,imax,idumm,nconfr,iconf,mm,k,im,nst,ifar,&
+ iter,irecv,isent,iw_pdb,nft0i,nft00_c,nft00,ifrom,ij,&
+ ireq,ireq2
+ real(kind=8) :: adif,p_cut,cutdifr,rmsdbc1c,time1i,ctdif1,xctdif,&
+ time2i,tstart,tend1
+!ccccccccccccccccccccccccccccccccccccccccccccccc
+ sync_iter=.true. !el
+ nft=0 !el
+ time_start=0.0d0
+ 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)
+
+! To minimize input conformation (bank conformation)
+! 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
+!!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
+
+! ntrial : number of trial conformations per seed.
+! ntry : total number of trial conformations including seed conformations.
+
+ idum2=-123
+! imax=2**31-1
+ imax=huge(0)
+ ENDIF
+
+ call mpi_bcast(jend,1,mpi_integer,0,CG_COMM,ierr)
+!ccccccccccccccccccccccccccccccccccccccc
+ do 300 jlee=1,jend
+!ccccccccccccccccccccccccccccccccccccccc
+ 331 continue
+ IF (ME.EQ.KING) THEN
+ if(sync_iter) goto 333
+ idum=- ran2(idum2)*imax
+ if(jlee.lt.jstart) goto 300
+
+! 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)
+
+!ccccccccccccccccccccccccccccccccccccccccccccccc
+ 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
+
+! start energy minimization
+ nconfr=max0(nconf+nadd,nodes-1)
+ if (sync_iter) nconf_in=0
+! king-emperor - feed input and sort output
+ write (iout,*) "NCONF_IN",nconf_in
+ m=0
+ if (nconf_in.gt.0) then
+! 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)
+! 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
+! To minimize input conformation (bank conformation)
+! Output to $mol.reminimized
+ if (irestart.lt.0) then
+ call reminimize(jlee)
+ return
+ endif
+ if (irestart.eq.1) goto 111
+! soldier - perform energy minimization
+ 334 call minim_jlee
+ ENDIF
+
+!cccccccccccccccccccccccccccccccccc
+! need to syncronize all procs
+ call mpi_barrier(CG_COMM,ierr)
+ if (ierr.ne.0) then
+ print *, ' cannot synchronize MPI'
+ stop
+ endif
+!cccccccccccccccccccccccccccccccccc
+
+ IF (ME.EQ.KING) THEN
+
+! print *,"ok after minim"
+ nstep=nstep+nconf
+ if(irestart.eq.2) then
+ nbank=nbank+nconf
+! ntbank=ntbank+nconf
+ if(ntbank.gt.ntbankm) ntbank=ntbankm
+ endif
+! print *,"ok before indexx"
+ if(iref.eq.0) then
+ call indexx(nconfr,etot,indx)
+ else
+! 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
+! cc/al
+! call indexx(nconfr,rmsn,indx)
+ endif
+! 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)
+!!bankt btene(im)=etot(m)
+!
+ 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)
+!!bankt btvar(i,j,k,im)=dihang(i,j,k,m)
+!
+ 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
+
+!d print *,"adif,xctdif,cutdifr"
+!d print *,adif,xctdif,cutdifr
+ nst=ntotal/ntrial/nseed
+ xctdif=(cutdif/ctdif1)**(-1.0/nst)
+ if(irestart.ge.1) call estimate_cutdif(adif,xctdif,cutdifr)
+! print *,"ok after estimate"
+
+ irestart=0
+
+ call write_rbank(jlee,adif,nft)
+ call write_bank(jlee,nft)
+!!bankt call write_bankt(jlee,nft)
+! 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
+
+!cccccccccccccccccccccccccccccccccccccc
+ do while (.not. finished)
+!cccccccccccccccccccccccccccccccccccccc
+!rc print *,"iter ", iter,' isent=',isent
+
+ IF (ME.EQ.KING) THEN
+! start energy minimization
+
+ if (isent.eq.0) then
+! king-emperor - select seeds & make var & feed input
+!d 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)
+!d 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
+
+!rc write(iout,*)'1:Calling FEEDIN NTRY',NTRY,' ntrial',ntrial
+!rc call feedin(ntry,nft)
+
+ isent=isent+1
+ if (isent.ge.nodes.or.iter.gt.0) then
+!t print *,'waiting ',MPI_WTIME()
+ irecv=irecv+1
+ call recv(0,ifrom,xout,eout,ind,timeout)
+!t print *,' ',irecv,' received from',ifrom,MPI_WTIME()
+ else
+ ifrom=ifrom+1
+ endif
+
+!t print *,'sending to',ifrom,MPI_WTIME()
+ call send(isent,ifrom,iter)
+!t print *,isent,' sent ',MPI_WTIME()
+
+! 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
+!--------------------------------------------------------------
+ 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
+!----------------- call update(ntry-nodes+1) -------------------
+ nstep=nstep+ntry-nseed-(nodes-1)
+ call refresh_bank(ntry-nodes+1)
+!!bankt call refresh_bankt(ntry-nodes+1)
+ else
+!----------------- call update(ntry) ---------------------------
+ iter=iter+1
+ print *,'UPDATING ',ntry,irecv
+ write(iout,*) 'UPDATING ',ntry
+ nstep=nstep+ntry-nseed
+ call refresh_bank(ntry)
+!!bankt call refresh_bankt(ntry)
+ endif
+!-----------------------------------------------------------------
+
+ call write_bank(jlee,nft)
+!!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
+!rc 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)
+!!bankt call refresh_bankt(ij-1)
+ goto 1002
+ endif
+! 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
+!rc 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
+
+!timeout call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
+!timeout call mpi_bcast(sync_iter,1,mpi_logical,0,
+!timeout & 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)
+!!bankt call refresh_bankt(nodes-1)
+ 1002 continue
+ call write_bank(jlee,nft)
+!!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
+! 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)
+!timeout call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
+!timeout 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
+
+!cccccccccccccccccccccccccccccccccccccc
+ enddo
+!cccccccccccccccccccccccccccccccccccccc
+
+ 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
+
+!ccccccccccccccccccccccccccccc
+ 300 continue
+!ccccccccccccccccccccccccccccc
+
+ return
+ end subroutine together
+!-----------------------------------------------------------------------------
+ subroutine feedin(nconf,nft)
+
+ use MPI_data
+ use geometry_data, only:nvar
+ use io, only:write_csa_pdb
+! 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'
+ real(kind=8),dimension(6*nres) :: xin,xout !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
+ real(kind=8),dimension(2) :: cout
+ integer,dimension(9) :: ind
+ integer,dimension(12) :: info
+ integer,dimension(mpi_status_size) :: muster
+! include 'COMMON.SETUP'
+ real(kind=8),parameter :: rad=1.745329252d-2
+ integer :: j,nconf,nft,mm,n,ierror,ierrcode,ierr,iw_pdb,&
+ man
+
+ print *,'FEEDIN: NCONF=',nconf
+ mm=0
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 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
+! 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
+! 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
+! 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)
+! 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
+! 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)
+! 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
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! no more input
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ do j=1,nodes-1
+! wait for a soldier
+ call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,&
+ CG_COMM,muster,ierr)
+!rc if (ierr.ne.0) go to 30
+! print *, ' receiving output from start # ',ind(1)
+ man=muster(mpi_source)
+! 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
+!rc if (ierr.ne.0) go to 30
+ call mpi_recv(xout,nvar,mpi_double_precision,&
+ man,idreal,CG_COMM,muster,ierr)
+! print *,nvar , ierr
+!rc if (ierr.ne.0) go to 30
+! 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)
+! 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
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 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 subroutine feedin
+!-----------------------------------------------------------------------------
+ subroutine getx(ind,xout,eout,cout,rad,iw_pdb,k)
+
+ use geometry_data
+ use energy_data
+ use compare, only: contact_fract
+ use MPI_data
+ include 'mpif.h'
+! 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'
+ integer,dimension(9) :: ind
+ real(kind=8),dimension(6*nres) :: xout !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
+!jlee
+ real(kind=8) :: przes(3),obr(3,3),cout(2)
+ logical :: non_conv
+ integer :: iw_pdb,k,j,ierror,ierrcode
+ real(kind=8) :: rad,co
+!jlee
+ 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
+! store ind()
+ do j=1,9
+ indb(k,j)=ind(j)
+ enddo
+! store energies
+ etot(k)=eout(1)
+! 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
+!jlee
+ if(iref.eq.0) then
+ iw_pdb=1
+!d write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4)')
+!d & ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' mv ',
+!d & ind(5),ind(4)
+ return
+ endif
+ call chainbuild
+! call dihang_to_c(dihang(1,1,1,k))
+! call fitsq(rms,c(1,1),crefjlee(1,1),nres,przes,obr,non_conv)
+! call fitsq(rms,c(1,2),crefjlee(1,2),nres-1,przes,obr,non_conv)
+! call fitsq(rms,c(1,nstart_seq),crefjlee(1,nstart_sup),
+! & nsup,przes,obr,non_conv)
+! 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)
+
+!d write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a5
+!d & ,0pf5.2,a5,f5.1,a,f6.3,a4,i3,i4)')
+!d & ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' rms ',
+!d & rmsn(k),' %NC ',pncn(k)*100,' cont.order',co,' mv ',
+!d & ind(5),ind(4)
+
+
+ if (rmsn(k).gt.rmscut.or.pncn(k).lt.pnccut) iw_pdb=0
+ return
+ end subroutine getx
+!-----------------------------------------------------------------------------
+ subroutine putx(xin,n,rad)
+
+ use geometry_data
+! 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'
+ integer :: n,m,j
+ real(kind=8),dimension(6*nres) :: xin !(maxvar) (maxvar=6*maxres)
+ real(kind=8) :: rad
+
+! 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
+! 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 subroutine putx
+!-----------------------------------------------------------------------------
+ subroutine putx2(xin,iff,n)
+
+ use geometry_data
+! 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'
+ integer :: n,m,j,i
+ real(kind=8),dimension(6*nres) :: xin !(maxvar) (maxvar=6*maxres)
+ integer,dimension(nres) :: iff !(maxres)
+
+! 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
+! set up array of variables
+ call geom_to_var(nvar,xin)
+
+ do i=1,nres
+ iff(i)=iff_in(i,n)
+ enddo
+ return
+ end subroutine putx2
+!-----------------------------------------------------------------------------
+ subroutine prune_bank(p_cut)
+
+ use MPI_data
+! 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'
+ integer :: k,j,i,m,ip,nprune
+ real(kind=8) :: p_cut,diff,ddmin
+!---------------------------
+! This subroutine prunes bank conformations using p_cut
+!---------------------------
+
+ 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 subroutine prune_bank
+!-----------------------------------------------------------------------------
+ subroutine reminimize(jlee)
+
+ use MPI_data
+! 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'
+ integer :: i,j,k,jlee,index,nft,ntry
+!---------------------------
+! This subroutine re-minimizes bank conformations:
+!---------------------------
+
+ 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)
+! 850 format(10f8.3)
+
+ return
+ end subroutine reminimize
+!-----------------------------------------------------------------------------
+ subroutine send(n,mm,it)
+
+ use MPI_data
+ use geometry_data, only: nvar
+ use control_data, only: vdisulf
+! 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'
+ real(kind=8),dimension(6*nres) :: xin,xout,xin2 !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
+ real(kind=8),dimension(2) :: cout
+ integer,dimension(9) :: ind
+ integer,dimension(nres) :: iff !(maxres)
+ integer,dimension(12) :: info
+ integer,dimension(mpi_status_size) :: muster
+! include 'COMMON.SETUP'
+ real(kind=8),parameter :: rad=1.745329252d-2
+ integer :: n,mm,it,ierr
+
+ if (isend2(n).eq.0) then
+! 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
+! 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 subroutine send
+!-----------------------------------------------------------------------------
+ subroutine recv(ihalt,man,xout,eout,ind,tout)
+
+ use MPI_data
+ use energy_data
+ use geometry_data, only: nvar
+ use control_data, only: vdisulf
+! 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'
+ real(kind=8),dimension(6*nres) :: xin,xout !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
+ real(kind=8),dimension(2) :: cout
+ integer,dimension(9) :: ind
+ integer,dimension(12) :: info
+ integer,dimension(mpi_status_size) :: muster
+! include 'COMMON.SETUP'
+ logical :: tout,flag
+ real(kind=8) :: tstart,tend1
+ real(kind=8),parameter :: twait=600.0d0
+ integer :: ihalt,man,ierr
+
+! 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.
+!_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)
+
+!timeout call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
+!timeout * CG_COMM,muster,ierr)
+! print *, ' receiving output from start # ',ind(1)
+!t print *,'receiving ',MPI_WTIME()
+!timeout man=muster(mpi_source)
+ call mpi_recv(ind,9,mpi_integer,man,idint,&
+ CG_COMM,muster,ierr)
+!timeout
+! 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
+! halt soldier
+ if(ihalt.eq.1) then
+! 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 subroutine recv
+!-----------------------------------------------------------------------------
+ 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 subroutine history_append
+!-----------------------------------------------------------------------------
+ subroutine alloc_CSA_arrays
+
+ use energy_data, only: ns
+
+ mxgr=2*nres
+
+ if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3))
+! commom.bank
+! common/varin/
+!el allocate(dihang_in(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
+ allocate(dihang_in(mxang,nres,mxch,5000)) !(mxang,maxres,mxch,mxio)
+ allocate(nss_in(mxio)) !(mxio)
+ allocate(iss_in(ns,mxio),jss_in(ns,mxio)) !(maxss,mxio)
+! common/minvar/
+ allocate(dihang(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
+ allocate(rmsn(mxio),pncn(mxio)) !(mxio)
+ allocate(etot(mxio)) !(mxio)
+ allocate(nss_out(mxio)) !(mxio)
+ allocate(iss_out(ns,mxio),jss_out(ns,mxio)) !(maxss,mxio)
+! common/bank/
+ allocate(rvar(mxang,nres,mxch,mxio),bvar(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
+ allocate(bene(mxio),rene(mxio),brmsn(mxio),rrmsn(mxio))
+ allocate(bpncn(mxio),rpncn(mxio)) !(mxio)
+ allocate(ibank(mxio),is(mxio),jbank(mxio)) !(mxio)
+ allocate(dij(mxio,mxio)) !(mxio,mxio)
+! common/bank_disulfid/
+ allocate(bvar_nss(mxio),bvar_ns(mxio)) !(mxio)
+ allocate(bvar_s(ns,mxio)) !(maxss,mxio)
+ allocate(bvar_ss(2,ns,mxio)) !(2,maxss,mxio)
+! common/mvstat/
+ allocate(movenx(mxio),movernx(mxio)) !(mxio)
+ allocate(nstatnx(0:mxmv,3),nstatnx_tot(0:mxmv,3)) !(0:mxmv,3)
+ allocate(indb(mxio,9)) !(mxio,9)
+ allocate(parent(3,mxio)) !(3,mxio)
+! common/send2/
+ allocate(isend2(mxio)) !(mxio)
+ allocate(iff_in(nres,mxio2)) !(maxres,mxio2)
+ allocate(dihang_in2(mxang,nres,mxch,mxio2)) !(mxang,maxres,mxch,mxio2)
+ allocate(idata(5,mxio)) !(5,mxio)
+! common.csa
+! common/alphaa/
+ allocate(ngroup(mxgr)) !(mxgr)
+ allocate(igroup(3,mxang,mxgr)) !(3,mxang,mxgr)
+! common.distfit
+! COMMON /frag/
+ allocate(bvar_frag(mxio,6)) !(mxio,6)
+ allocate(hvar_frag(mxio,3),lvar_frag(mxio,3),svar_frag(mxio,3)) !(mxio,3)
+ allocate(avar_frag(mxio,5)) !(mxio,5)
+! commom.hairpin
+! common /spinka/
+ allocate(nharp_seed(nseed),nharp_use(nseed)) !(max_seed)
+ allocate(iharp_seed(4,nres/3,nseed)) !(4,maxres/3,max_seed)
+ allocate(iharp_use(0:4,nres/3,nseed)) !(0:4,maxres/3,max_seed)
+
+ return
+ end subroutine alloc_CSA_arrays
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module csa
projects / unres4.git / blobdiff