source/unres/CSA.F90

   1       module csa
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use math
   6       use random
   7       use geometry_data, only: nres,rad2deg
   8       use geometry
   9       use energy
  10       use MPI_
  11       use csa_data
  12       use compare
  13       use io_config
  14
  15       implicit none
  16 !-----------------------------------------------------------------------------
  17 ! commom.bank
  18 !      common/varin/
  19 !      real(kind=8),dimension(:,:,:,:),allocatable :: dihang_in !(mxang,maxres,mxch,mxio)
  20       integer,dimension(:),allocatable :: nss_in !(mxio)
  21       integer,dimension(:,:),allocatable :: iss_in,jss_in !(maxss,mxio)
  22 !      common/minvar/
  23       real(kind=8),dimension(:,:,:,:),allocatable :: dihang !(mxang,maxres,mxch,mxio)
  24       real(kind=8),dimension(:),allocatable :: etot!,rmsn,pncn !(mxio)
  25       integer,dimension(:),allocatable :: nss_out !(mxio)
  26       integer,dimension(:,:),allocatable ::iss_out,jss_out !(maxss,mxio)
  27 !      common/bank/
  28 !      real(kind=8),dimension(:,:,:,:),allocatable :: bvar !(mxang,maxres,mxch,mxio)
  29 !      real(kind=8),dimension(:),allocatable :: bene,rene,&
  30 !       brmsn,rrmsn,bpncn,rpncn !(mxio)
  31       integer,dimension(:),allocatable :: is,jbank !(mxio)
  32       real(kind=8) :: avedif,difmin,ebmin,ebmax,ebmaxt!,&
  33 !       dele,difcut,cutdif,rmscut,pnccut
  34       real(kind=8),dimension(:,:),allocatable :: dij !(mxio,mxio)
  35 !      common/bank_disulfid/
  36 !      common/mvstat/
  37       integer,dimension(:),allocatable :: movenx,movernx !(mxio)
  38       integer,dimension(:,:),allocatable :: nstatnx,nstatnx_tot !(0:mxmv,3)
  39       integer,dimension(:,:),allocatable :: indb !(mxio,9)
  40       integer,dimension(:,:),allocatable :: parent !(3,mxio)
  41 !      common/send2/
  42       integer,dimension(:),allocatable :: isend2 !(mxio)
  43       integer,dimension(:,:),allocatable :: iff_in !(maxres,mxio2)
  44       integer,dimension(:,:,:,:),allocatable :: dihang_in2 !(mxang,maxres,mxch,mxio2)
  45       integer,dimension(:,:),allocatable :: idata !(5,mxio)
  46 !-----------------------------------------------------------------------------
  47 ! common.csa
  48 !      integer :: irestart,ndiff
  49 !      common/alphaa/
  50       integer,dimension(:),allocatable :: ngroup !(mxgr)
  51       integer,dimension(:,:,:),allocatable :: igroup !(3,mxang,mxgr)
  52       integer :: ntotgr!,numch
  53 !      common/csa_input/
  54 !      common/dih_control/
  55 !      real(kind=8) :: rdih_bias
  56       logical :: ldih_bias
  57 !-----------------------------------------------------------------------------
  58 ! common.distfit
  59 !      COMMON /frag/
  60       integer,dimension(:,:),allocatable :: bvar_frag !(mxio,6)
  61       integer,dimension(:,:),allocatable :: hvar_frag,lvar_frag,svar_frag !(mxio,3)
  62       integer,dimension(:,:),allocatable :: avar_frag !(mxio,5)
  63 !-----------------------------------------------------------------------------
  64 ! commom.hairpin
  65 !      common /spinka/
  66       integer :: nharp_tot
  67       integer,dimension(:),allocatable :: nharp_seed,nharp_use  !(max_seed)
  68       integer,dimension(:,:,:),allocatable :: iharp_seed        !(4,maxres/3,max_seed)
  69       integer,dimension(:,:,:),allocatable :: iharp_use         !(0:4,maxres/3,max_seed)
  70 !-----------------------------------------------------------------------------
  71 ! Maximum number of moves (n1-n8)
  72       integer,parameter :: mxmv=18
  73 !-----------------------------------------------------------------------------
  74 !
  75 !
  76 !-----------------------------------------------------------------------------
  77       contains
  78 !-----------------------------------------------------------------------------
  79 ! bank.F
  80 !-----------------------------------------------------------------------------
  81       subroutine refresh_bank(ntrial)
  82
  83 !      implicit real*8 (a-h,o-z)
  84 !      include 'DIMENSIONS'
  85       include 'mpif.h'
  86 !      include 'COMMON.CSA'
  87 !      include 'COMMON.BANK'
  88 !      include 'COMMON.IOUNITS'
  89 !      include 'COMMON.CHAIN'
  90 !      include 'COMMON.VAR'
  91 !      include 'COMMON.CONTROL'
  92       character :: chacc
  93       integer :: iaccn,ntrial
  94       real(kind=8) :: l_diff(mxio),denep
  95       integer :: i,j,n,m,i1,idmin
  96       real(kind=8) :: del_ene
  97
  98       do i=0,mxmv
  99         do j=1,3
 100           nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
 101           nstatnx(i,j)=0
 102         enddo
 103       enddo
 104
 105 ! loop over all newly obtained conformations
 106       do n=1,ntrial
 107        chacc=' '
 108        iaccn=0
 109        nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1
 110 !ccccccccccccccccccccccccccccccccccccccccccc
 111 !jlee
 112        if(iref.ne.0) then
 113         if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100
 114        endif
 115 !jlee
 116        if(etot(n).gt.ebmax) goto 100
 117 ! Find the conformation closest to the conformation n in the bank
 118        difmin=9.d9
 119        do m=1,nbank
 120         call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m))
 121         if(l_diff(m).lt.difmin) then
 122          difmin=l_diff(m)
 123          idmin=m
 124         endif
 125        enddo
 126
 127        if(difmin.lt.cutdif) then
 128 ! n is redundant to idmin
 129         if(etot(n).lt.bene(idmin)) then
 130          if(etot(n).lt.bene(idmin)-0.01d0) then
 131           ibank(idmin)=0
 132           jbank(idmin)=0
 133          endif
 134          denep=bene(idmin)-etot(n)
 135          call replace_bvar(idmin,n)
 136 !rc Update dij
 137          do i1=1,nbank
 138            if (i1.ne.idmin) then
 139             dij(i1,idmin)=l_diff(i1)
 140             dij(idmin,i1)=l_diff(i1)
 141            endif
 142          enddo
 143          chacc='c'
 144          iaccn=idmin
 145          nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1
 146          if(idmin.eq.ibmax) call find_max
 147         endif
 148        else
 149 ! got new conformation
 150         del_ene=0.0d0
 151         if(ebmax-ebmin.gt.del_ene) then
 152          denep=ebmax-etot(n)
 153          call replace_bvar(ibmax,n)
 154 !rc Update dij
 155          do i1=1,nbank
 156            if (i1.ne.ibmax) then
 157             dij(i1,ibmax)=l_diff(i1)
 158             dij(ibmax,i1)=l_diff(i1)
 159            endif
 160          enddo
 161          chacc='f'
 162          iaccn=ibmax
 163          nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1
 164          ibank(ibmax)=0
 165          jbank(ibmax)=0
 166          call find_max
 167         else
 168          if(del_ene.lt.0.0001) then
 169           write (iout,*) 'ERROR in refresh_bank: '
 170           write (iout,*) 'ebmax: ',ebmax
 171           write (iout,*) 'ebmin: ',ebmin
 172           write (iout,*) 'del_ene: ',del_ene
 173 !rc          call mpi_abort(mpi_comm_world,ierror,ierrcode)
 174          endif
 175 !jp nbmax is never defined so condition below is always false
 176 !         if(nbank.lt.nbmax) then
 177 !          nbank=nbank+1
 178 !          call replace_bvar(nbank,n)
 179 !          ibank(nbank)=0
 180 !          jbank(nbank)=0
 181 !         else
 182           call replace_bvar(ibmax,n)
 183           ibank(ibmax)=0
 184           jbank(ibmax)=0
 185           call find_max
 186 !         endif
 187         endif
 188        endif
 189 !ccccccccccccccccccccccccccccccccccccccccccc
 190   100 continue
 191        if (iaccn.eq.0) then
 192         if (iref.eq.0) then
 193          write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)') &
 194           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
 195           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9)
 196         else
 197       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0)') &
 198           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
 199           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
 200           ' rms ',rmsn(n),' %NC ',pncn(n)*100
 201         endif
 202        else
 203         if (iref.eq.0) then
 204         write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,1x,a1,i4,0pf8.1,0pf8.1)') &
 205           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
 206           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
 207           chacc,iaccn,difmin,denep
 208         else
 209          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)') &
 210           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
 211           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
 212           ' rms ',rmsn(n),' %NC ',pncn(n)*100,&
 213            chacc,iaccn,difmin,denep
 214         endif
 215        endif
 216       enddo
 217 ! end of loop over all newly obtained conformations
 218       do i=0,mxmv
 219         if(nstatnx(i,1).ne.0) then
 220          if (i.le.9) then
 221          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
 222                 '## N',i,' total=',nstatnx(i,1),&
 223                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
 224                 ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
 225          else
 226          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
 227                 '##N',i,' total=',nstatnx(i,1),&
 228                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
 229                 ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
 230          endif
 231         else
 232          if (i.le.9) then
 233          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
 234                 '## N',i,' total=',nstatnx(i,1),&
 235                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
 236                 ' %acc',0.0
 237          else
 238          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
 239                 '##N',i,' total=',nstatnx(i,1),&
 240                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
 241                 ' %acc',0.0
 242          endif
 243         endif
 244       enddo
 245       call flush(iout)
 246 !rc Update dij
 247 !rc moved up, saves some get_diff12 calls
 248 !rc
 249 !rc      do i1=1,nbank-1
 250 !rc       do i2=i1+1,nbank
 251 !rc        if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
 252 !rc         call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
 253 !rc         dij(i1,i2)=diff
 254 !rc         dij(i2,i1)=diff
 255 !rc        endif
 256 !rc       enddo
 257 !rc      enddo
 258
 259       do i=1,nbank
 260        jbank(i)=1
 261       enddo
 262
 263       return
 264       end subroutine refresh_bank
 265 !-----------------------------------------------------------------------------
 266       subroutine replace_bvar(iold,inew)
 267
 268       use control_data, only: vdisulf
 269       use energy_data, only: ns,iss
 270 !      implicit real*8 (a-h,o-z)
 271 !      include 'DIMENSIONS'
 272       include 'mpif.h'
 273 !      include 'COMMON.IOUNITS'
 274 !      include 'COMMON.CSA'
 275 !      include 'COMMON.BANK'
 276 !      include 'COMMON.CHAIN'
 277 !      include 'COMMON.CONTROL'
 278 !      include 'COMMON.SBRIDGE'
 279       integer :: iold,inew,ierror,ierrcode,i,j,k
 280
 281       if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1) &
 282         then
 283         write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold,&
 284         ' INEW',inew
 285         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 286       endif
 287       do k=1,numch
 288        do j=2,nres-1
 289         do i=1,4
 290          bvar(i,j,k,iold)=dihang(i,j,k,inew)
 291         enddo
 292        enddo
 293       enddo
 294       bene(iold)=etot(inew)
 295       brmsn(iold)=rmsn(inew)
 296       bpncn(iold)=pncn(inew)
 297
 298       if(bene(iold).lt.ebmin) then
 299         ebmin=bene(iold)
 300         ibmin=iold
 301       endif
 302
 303       if(vdisulf) then
 304         bvar_nss(iold)=nss_out(inew)
 305 !d        write(iout,*) 'SS BANK',iold,bvar_nss(iold)
 306         do i=1,bvar_nss(iold)
 307           bvar_ss(1,i,iold)=iss_out(i,inew)
 308           bvar_ss(2,i,iold)=jss_out(i,inew)
 309 !d          write(iout,*) 'SS',bvar_ss(1,i,iold)-nres,
 310 !d     &          bvar_ss(2,i,iold)-nres
 311         enddo
 312
 313         bvar_ns(iold)=ns-2*bvar_nss(iold)
 314 !d        write(iout,*) 'CYS #free ', bvar_ns(iold)
 315            k=0
 316            do i=1,ns
 317              j=1
 318              do while( iss(i).ne.iss_out(j,inew)-nres .and. &
 319                        iss(i).ne.jss_out(j,inew)-nres .and. &
 320                        j.le.nss_out(inew))
 321                 j=j+1
 322              enddo
 323              if (j.gt.nss_out(inew)) then
 324                k=k+1
 325                bvar_s(k,iold)=iss(i)
 326              endif
 327            enddo
 328 !d         write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold))
 329       endif
 330
 331       return
 332       end subroutine replace_bvar
 333 !-----------------------------------------------------------------------------
 334       subroutine save_is(ind)
 335
 336 !      implicit real*8 (a-h,o-z)
 337 !      include 'DIMENSIONS'
 338       include 'mpif.h'
 339 !      include 'COMMON.IOUNITS'
 340 !      include 'COMMON.CSA'
 341 !      include 'COMMON.BANK'
 342 !      include 'COMMON.CHAIN'
 343       integer :: ind,i,j,k,index,ierror,ierrcode
 344
 345       index=nbank+ind
 346 !     print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind)
 347       if (index.gt.mxio .or. index.lt.1 .or. &
 348         is(ind).gt.mxio .or. is(ind).lt.1) then
 349         write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index,&
 350         ' IND',ind,' IS',is(ind)
 351         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 352       endif
 353       do k=1,numch
 354        do j=2,nres-1
 355         do i=1,4
 356          bvar(i,j,k,index)=bvar(i,j,k,is(ind))
 357         enddo
 358        enddo
 359       enddo
 360       bene(index)=bene(is(ind))
 361       ibank(is(ind))=1
 362
 363       return
 364       end subroutine save_is
 365 !-----------------------------------------------------------------------------
 366       subroutine select_is(n,ifar,idum)
 367
 368 !      implicit real*8 (a-h,o-z)
 369 !      include 'DIMENSIONS'
 370 !      include 'COMMON.CSA'
 371 !      include 'COMMON.BANK'
 372       integer,dimension(mxio) :: itag
 373       real(kind=8),dimension(mxio) :: adiff
 374       integer :: n,ifar,idum,i,iusesv,imade
 375
 376       iuse=0
 377       do i=1,nbank
 378        if(ibank(i).eq.0) then
 379         iuse=iuse+1
 380         itag(iuse)=i
 381        endif
 382       enddo
 383       iusesv=iuse
 384
 385       if(iuse.eq.0) then
 386        icycle=icycle+1
 387        do i=1,nbank
 388         if(ibank(i).eq.2) then
 389          ibank(i)=1
 390         else
 391          ibank(i)=0
 392         endif
 393        enddo
 394        imade=0
 395        call get_is(idum,ifar,n,imade,0)
 396 !test3       call get_is_max(idum,ifar,n,imade,0)
 397       else if(iuse.eq.n) then
 398        do i=1,iuse
 399         is(i)=itag(i)
 400         call save_is(i)
 401        enddo
 402       else if(iuse.lt.n) then
 403 !      if(icycle.eq.0) then
 404 !       do i=1,n
 405 !        ind=mod(i-1,iuse)+1
 406 !        is(i)=itag(ind)
 407 !        call save_is(i)
 408 !       enddo
 409 !      else
 410 !      endif
 411        do i=1,iuse
 412         is(i)=itag(i)
 413         call save_is(i)
 414        enddo
 415        imade=iuse
 416 !      call get_is_ran(idum,n,imade,1)
 417        call get_is(idum,ifar,n,imade,1)
 418 !test3       call get_is_max(idum,ifar,n,imade,1)
 419 !      if(iusesv.le.n/10) then
 420        if(iusesv.le.0) then
 421         icycle=icycle+1
 422         do i=1,nbank
 423 !        if(ibank(i).eq.2) then
 424 !         ibank(i)=1
 425          if(ibank(i).ge.2) then
 426           ibank(i)=ibank(i)-1
 427          else
 428           ibank(i)=0
 429          endif
 430         enddo
 431        endif
 432       else
 433        imade=0
 434        call get_is(idum,ifar,n,imade,0)
 435 !test3       call get_is_max(idum,ifar,n,imade,0)
 436       endif
 437       iuse=iusesv
 438
 439       return
 440       end subroutine select_is
 441 !-----------------------------------------------------------------------------
 442       subroutine get_is_ran(idum,n,imade,k)
 443
 444 !      implicit real*8 (a-h,o-z)
 445 !      include 'DIMENSIONS'
 446 !      include 'COMMON.CSA'
 447 !      include 'COMMON.BANK'
 448 !      real(kind=4) :: ran1,ran2
 449       integer,dimension(mxio) :: itag
 450       real(kind=8),dimension(mxio) :: adiff
 451       integer :: idum,n,imade,k,j,i,iran
 452
 453       do j=imade+1,n
 454        iuse=0
 455        do i=1,nbank
 456         if(ibank(i).eq.k) then
 457          iuse=iuse+1
 458          itag(iuse)=i
 459         endif
 460        enddo
 461        iran=iuse*  ran1(idum)+1
 462        is(j)=itag(iran)
 463        call save_is(j)
 464       enddo
 465
 466       return
 467       end subroutine get_is_ran
 468 !-----------------------------------------------------------------------------
 469       subroutine get_is(idum,ifar,n,imade,k)
 470
 471 !      implicit real*8 (a-h,o-z)
 472 !      include 'DIMENSIONS'
 473 !      include 'COMMON.CSA'
 474 !      include 'COMMON.BANK'
 475 !      real(kind=4) :: ran1,ran2
 476       integer,dimension(mxio) :: itag
 477       real(kind=8),dimension(mxio) :: adiff
 478       integer :: idum,ifar,n,imade,k,i,iran
 479
 480       iuse=0
 481       do i=1,nbank
 482        if(ibank(i).eq.k) then
 483         iuse=iuse+1
 484         itag(iuse)=i
 485        endif
 486       enddo
 487       iran=iuse*  ran1(idum)+1
 488       imade=imade+1
 489       is(imade)=itag(iran)
 490       call save_is(imade)
 491
 492       do i=imade+1,ifar-1
 493        if(icycle.eq.-1) then
 494         call select_iseed_max(i,k)
 495        else
 496         call select_iseed_min(i,k)
 497 !test4  call select_iseed_max(i,k)
 498        endif
 499        call save_is(i)
 500       enddo
 501
 502       do i=ifar,n
 503        call select_iseed_far(i,k)
 504        call save_is(i)
 505       enddo
 506
 507       return
 508       end subroutine get_is
 509 !-----------------------------------------------------------------------------
 510       subroutine select_iseed_max(imade1,ik)
 511
 512 !      implicit real*8 (a-h,o-z)
 513 !      include 'DIMENSIONS'
 514 !      include 'COMMON.CSA'
 515 !      include 'COMMON.BANK'
 516       integer,dimension(mxio) :: itag
 517       real(kind=8),dimension(mxio) :: adiff
 518       integer :: imade1,ik,i,n,imade,m,itagi
 519       real(kind=8) :: difmax,diff,emax,benei,diffmn
 520
 521       iuse=0
 522       avedif=0.d0
 523       difmax=0.d0
 524       do n=1,nbank
 525        if(ibank(n).eq.ik) then
 526         iuse=iuse+1
 527         diffmn=9.d190
 528         do imade=1,imade1-1
 529 !        m=nbank+imade
 530 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 531          m=is(imade)
 532          diff=dij(n,m)
 533          if(diff.lt.diffmn) diffmn=diff
 534         enddo
 535         if(diffmn.gt.difmax) difmax=diffmn
 536         adiff(iuse)=diffmn
 537         itag(iuse)=n
 538         avedif=avedif+diffmn
 539        endif
 540       enddo
 541
 542       avedif=avedif/iuse
 543 !     avedif=(avedif+difmax)/2
 544       emax=-9.d190
 545       do i=1,iuse
 546        if(adiff(i).ge.avedif) then
 547         itagi=itag(i)
 548         benei=bene(itagi)
 549         if(benei.gt.emax) then
 550          emax=benei
 551          is(imade1)=itagi
 552         endif
 553        endif
 554       enddo
 555
 556       if(ik.eq.0) iuse=iuse-1
 557
 558       return
 559       end subroutine select_iseed_max
 560 !-----------------------------------------------------------------------------
 561       subroutine select_iseed_min(imade1,ik)
 562
 563 !      implicit real*8 (a-h,o-z)
 564 !      include 'DIMENSIONS'
 565 !      include 'COMMON.CSA'
 566 !      include 'COMMON.BANK'
 567       integer,dimension(mxio) :: itag
 568       real(kind=8),dimension(mxio) :: adiff
 569       integer :: imade1,ik,n,imade,m,i,itagi
 570       real(kind=8) :: difmax,diff,diffmn,emin,benei
 571
 572       iuse=0
 573       avedif=0.d0
 574       difmax=0.d0
 575       do n=1,nbank
 576        if(ibank(n).eq.ik) then
 577         iuse=iuse+1
 578         diffmn=9.d190
 579         do imade=1,imade1-1
 580 !        m=nbank+imade
 581 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 582          m=is(imade)
 583          diff=dij(n,m)
 584          if(diff.lt.diffmn) diffmn=diff
 585         enddo
 586         if(diffmn.gt.difmax) difmax=diffmn
 587         adiff(iuse)=diffmn
 588         itag(iuse)=n
 589         avedif=avedif+diffmn
 590        endif
 591       enddo
 592
 593       avedif=avedif/iuse
 594 !     avedif=(avedif+difmax)/2
 595       emin=9.d190
 596       do i=1,iuse
 597 !      print *,"i, adiff(i),avedif : ",i,adiff(i),avedif
 598        if(adiff(i).ge.avedif) then
 599         itagi=itag(i)
 600         benei=bene(itagi)
 601 !       print *,"i, benei,emin : ",i,benei,emin
 602         if(benei.lt.emin) then
 603          emin=benei
 604          is(imade1)=itagi
 605         endif
 606        endif
 607       enddo
 608
 609       if(ik.eq.0) iuse=iuse-1
 610
 611 !     print *, "exiting select_iseed_min",is(imade1)
 612
 613       return
 614       end subroutine select_iseed_min
 615 !-----------------------------------------------------------------------------
 616       subroutine select_iseed_far(imade1,ik)
 617
 618 !      implicit real*8 (a-h,o-z)
 619 !      include 'DIMENSIONS'
 620 !      include 'COMMON.CSA'
 621 !      include 'COMMON.BANK'
 622       integer :: imade1,ik,n,imade,m
 623       real(kind=8) :: dmax,diffmn,diff
 624
 625       dmax=-9.d190
 626       do n=1,nbank
 627        if(ibank(n).eq.ik) then
 628         diffmn=9.d190
 629         do imade=1,imade1-1
 630 !        m=nbank+imade
 631 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 632          m=is(imade)
 633          diff=dij(n,m)
 634          if(diff.lt.diffmn) diffmn=diff
 635         enddo
 636        endif
 637        if(diffmn.gt.dmax) then
 638         dmax=diffmn
 639         is(imade1)=n
 640        endif
 641       enddo
 642
 643       return
 644       end subroutine select_iseed_far
 645 !-----------------------------------------------------------------------------
 646       subroutine find_min
 647
 648 !      implicit real*8 (a-h,o-z)
 649 !      include 'DIMENSIONS'
 650 !      include 'COMMON.CSA'
 651 !      include 'COMMON.BANK'
 652       integer :: i
 653       real(kind=8) :: benei
 654
 655       ebmin=9.d190
 656
 657       do i=1,nbank
 658        benei=bene(i)
 659        if(benei.lt.ebmin) then
 660         ebmin=benei
 661         ibmin=i
 662        endif
 663       enddo
 664
 665       return
 666       end subroutine find_min
 667 !-----------------------------------------------------------------------------
 668       subroutine find_max
 669
 670 !      implicit real*8 (a-h,o-z)
 671 !      include 'DIMENSIONS'
 672 !      include 'COMMON.CSA'
 673 !      include 'COMMON.BANK'
 674       integer :: i
 675       real(kind=8) :: benei
 676
 677       ebmax=-9.d190
 678
 679       do i=1,nbank
 680        benei=bene(i)
 681        if(benei.gt.ebmax) then
 682         ebmax=benei
 683         ibmax=i
 684        endif
 685       enddo
 686
 687       return
 688       end subroutine find_max
 689 !-----------------------------------------------------------------------------
 690       subroutine get_diff
 691
 692 !      implicit real*8 (a-h,o-z)
 693 !      include 'DIMENSIONS'
 694 !      include 'COMMON.CSA'
 695 !      include 'COMMON.BANK'
 696       integer :: i,i1,i2
 697       real(kind=8) :: tdiff,difmin,diff
 698
 699       tdiff=0.d0
 700       difmin=9.d190
 701       do i1=1,nbank-1
 702        do i2=i1+1,nbank
 703         if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
 704          call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
 705          dij(i1,i2)=diff
 706          dij(i2,i1)=diff
 707         else
 708          diff=dij(i1,i2)
 709         endif
 710         tdiff=tdiff+diff
 711         if(diff.lt.difmin) difmin=diff
 712        enddo
 713        dij(i1,i1)=0.0
 714       enddo
 715
 716       do i=1,nbank
 717        jbank(i)=1
 718       enddo
 719
 720       avedif=tdiff/nbank/(nbank-1)*2
 721
 722       return
 723       end subroutine get_diff
 724 !-----------------------------------------------------------------------------
 725       subroutine estimate_cutdif(adif,xct,cutdifr)
 726
 727 !      implicit real*8 (a-h,o-z)
 728 !      include 'DIMENSIONS'
 729 !      include 'COMMON.CSA'
 730 !      include 'COMMON.BANK'
 731       integer :: nexp
 732       real(kind=8) :: adif,xct,cutdifr,ctdif1,exponent
 733
 734       ctdif1=adif/cut2
 735
 736       exponent = cutdifr*cut1/adif
 737       exponent = dlog(exponent)/dlog(xct)
 738
 739       nexp=exponent+0.25
 740       cutdif= adif/cut1*xct**nexp
 741       if(cutdif.lt.ctdif1) cutdif=ctdif1
 742
 743       return
 744       end subroutine estimate_cutdif
 745 !-----------------------------------------------------------------------------
 746       subroutine get_is_max(idum,ifar,n,imade,k)
 747
 748 !      implicit real*8 (a-h,o-z)
 749 !      include 'DIMENSIONS'
 750 !      include 'COMMON.CSA'
 751 !      include 'COMMON.BANK'
 752       integer :: idum,ifar,n,imade,k,i,j
 753       real(kind=8) :: emax
 754
 755       do i=imade+1,n
 756        emax=-9.d190
 757        do j=1,nbank
 758         if(ibank(j).eq.k .and. bene(j).gt.emax) then
 759            emax=bene(j)
 760            is(i)=j
 761         endif
 762        enddo
 763        call save_is(i)
 764       enddo
 765
 766       return
 767       end subroutine get_is_max
 768 !-----------------------------------------------------------------------------
 769 ! csa.f
 770 !-----------------------------------------------------------------------------
 771       subroutine make_array
 772
 773       use energy_data, only: itype
 774 !      implicit real*8 (a-h,o-z)
 775 !      include 'DIMENSIONS'
 776 !      include 'COMMON.IOUNITS'
 777 !      include 'COMMON.CHAIN'
 778 !      include 'COMMON.INTERACT'
 779 !      include 'COMMON.CSA'
 780        integer :: k,j,i,indg
 781 !cccccccccccccccccccccccc
 782 ! Level-2: group
 783 !cccccccccccccccccccccccc
 784
 785       indg=0
 786       do k=1,numch
 787 !cccccccccccccccccccccccccccccccccccccccc
 788 ! Groups the THETAs and the GAMMAs
 789        do j=2,nres-1
 790          indg=indg+1
 791          if (j.lt.nres-1) then
 792            ngroup(indg)=2
 793          else
 794            ngroup(indg)=1
 795          endif
 796          do i=1,ngroup(indg)
 797           igroup(1,i,indg)=i
 798           igroup(2,i,indg)=j
 799           igroup(3,i,indg)=k
 800          enddo
 801        enddo
 802 !cccccccccccccccccccccccccccccccccccccccc
 803       enddo
 804 ! Groups the ALPHAs and the BETAs
 805       do k=1,numch
 806        do j=2,nres-1
 807         if(itype(j).ne.10) then
 808          indg=indg+1
 809          ngroup(indg)=2
 810          do i=1,ngroup(indg)
 811           igroup(1,i,indg)=i+2
 812           igroup(2,i,indg)=j
 813           igroup(3,i,indg)=k
 814          enddo
 815         endif
 816        enddo
 817       enddo
 818
 819       ntotgr=indg
 820       write(iout,*)
 821       write(iout,*) "# of groups: ",ntotgr
 822       do i=1,ntotgr
 823        write(iout,41) i,ngroup(i),((igroup(k,j,i),k=1,3),j=1,ngroup(i))
 824       enddo
 825 !      close(iout)
 826
 827    40 format(i3,3x,3i3)
 828    41 format(2i3,3x,6(3i3,2x))
 829
 830       return
 831       end subroutine make_array
 832 !-----------------------------------------------------------------------------
 833       subroutine make_ranvar(n,m,idum)
 834
 835       use geometry_data
 836 !      implicit real*8 (a-h,o-z)
 837 !      include 'DIMENSIONS'
 838 !      include 'COMMON.IOUNITS'
 839 !      include 'COMMON.CHAIN'
 840 !      include 'COMMON.VAR'
 841 !      include 'COMMON.BANK'
 842       integer :: n,m,j,idum,itrial,jeden
 843
 844 ! al      m=0
 845       print *,'HOHOHOHO Make_RanVar!!!!!',n,m
 846       itrial=0
 847       do while(m.lt.n .and. itrial.le.10000)
 848         itrial=itrial+1
 849         jeden=1
 850         call gen_rand_conf(jeden,*10)
 851 !        call intout
 852         m=m+1
 853         do j=2,nres-1
 854           dihang_in(1,j,1,m)=theta(j+1)
 855           dihang_in(2,j,1,m)=phi(j+2)
 856           dihang_in(3,j,1,m)=alph(j)
 857           dihang_in(4,j,1,m)=omeg(j)
 858         enddo
 859         dihang_in(2,nres-1,1,m)=0.0d0
 860         goto 20
 861   10    write (iout,*) 'Failed to generate conformation #',m+1,&
 862          ' itrial=',itrial
 863   20    continue
 864       enddo
 865       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
 866
 867       return
 868       end subroutine make_ranvar
 869 !-----------------------------------------------------------------------------
 870       subroutine make_ranvar_reg(n,idum)
 871
 872       use geometry_data
 873 !      implicit real*8 (a-h,o-z)
 874 !      include 'DIMENSIONS'
 875 !      include 'COMMON.IOUNITS'
 876 !      include 'COMMON.CHAIN'
 877 !      include 'COMMON.VAR'
 878 !      include 'COMMON.BANK'
 879 !      include 'COMMON.GEO'
 880       integer :: n,idum,j,m,itrial,jeden
 881
 882       m=0
 883       print *,'HOHOHOHO Make_RanVar!!!!!'
 884       itrial=0
 885       do while(m.lt.n .and. itrial.le.10000)
 886         itrial=itrial+1
 887         jeden=1
 888         call gen_rand_conf(jeden,*10)
 889 !        call intout
 890         m=m+1
 891         do j=2,nres-1
 892           dihang_in(1,j,1,m)=theta(j+1)
 893           dihang_in(2,j,1,m)=phi(j+2)
 894           dihang_in(3,j,1,m)=alph(j)
 895           dihang_in(4,j,1,m)=omeg(j)
 896           if(m.le.n*0.1) then
 897            dihang_in(1,j,1,m)=90.0*deg2rad
 898            dihang_in(2,j,1,m)=50.0*deg2rad
 899           endif
 900         enddo
 901         dihang_in(2,nres-1,1,m)=0.0d0
 902         goto 20
 903   10    write (iout,*) 'Failed to generate conformation #',m+1,&
 904          ' itrial=',itrial
 905   20    continue
 906       enddo
 907       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
 908
 909       return
 910       end subroutine make_ranvar_reg
 911 !-----------------------------------------------------------------------------
 912 ! diff12.f
 913 !-----------------------------------------------------------------------------
 914       subroutine get_diff12(aarray,barray,diff)
 915
 916 !      implicit real*8 (a-h,o-z)
 917 !      include 'DIMENSIONS'
 918 !      include 'COMMON.CSA'
 919 !      include 'COMMON.BANK'
 920 !      include 'COMMON.CHAIN'
 921 !      include 'COMMON.GEO'
 922       integer :: k,j,i
 923       real(kind=8),dimension(mxang,nres,mxch) :: aarray,barray  !(mxang,maxres,mxch)
 924       real(kind=8) :: diff,dif
 925
 926       diff=0.d0
 927       do k=1,numch
 928        do j=2,nres-1
 929 !       do i=1,4
 930 !        do i=1,2
 931         do i=1,ndiff
 932          dif=rad2deg*dabs(aarray(i,j,k)-barray(i,j,k))
 933          if(dif.gt.180.) dif=360.-dif
 934          if (dif.gt.diffcut) diff=diff+dif
 935         enddo
 936        enddo
 937       enddo
 938
 939       return
 940       end subroutine get_diff12
 941 !-----------------------------------------------------------------------------
 942 ! indexx.f
 943 !-----------------------------------------------------------------------------
 944       subroutine indexx(n,arr,indx)
 945
 946 !      implicit real*8 (a-h,o-z)
 947       INTEGER :: n,indx(n)
 948       REAL(kind=8) :: arr(n)
 949 !     PARAMETER (M=7,NSTACK=50)
 950       integer,PARAMETER :: M=7,NSTACK=500
 951       INTEGER :: i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
 952       REAL(kind=8) :: a
 953
 954       do 11 j=1,n
 955         indx(j)=j
 956 11    continue
 957       jstack=0
 958       l=1
 959       ir=n
 960 1     if(ir-l.lt.M)then
 961         do 13 j=l+1,ir
 962           indxt=indx(j)
 963           a=arr(indxt)
 964           do 12 i=j-1,1,-1
 965             if(arr(indx(i)).le.a)goto 2
 966             indx(i+1)=indx(i)
 967 12        continue
 968           i=0
 969 2         indx(i+1)=indxt
 970 13      continue
 971         if(jstack.eq.0)return
 972         ir=istack(jstack)
 973         l=istack(jstack-1)
 974         jstack=jstack-2
 975       else
 976         k=(l+ir)/2
 977         itemp=indx(k)
 978         indx(k)=indx(l+1)
 979         indx(l+1)=itemp
 980         if(arr(indx(l+1)).gt.arr(indx(ir)))then
 981           itemp=indx(l+1)
 982           indx(l+1)=indx(ir)
 983           indx(ir)=itemp
 984         endif
 985         if(arr(indx(l)).gt.arr(indx(ir)))then
 986           itemp=indx(l)
 987           indx(l)=indx(ir)
 988           indx(ir)=itemp
 989         endif
 990         if(arr(indx(l+1)).gt.arr(indx(l)))then
 991           itemp=indx(l+1)
 992           indx(l+1)=indx(l)
 993           indx(l)=itemp
 994         endif
 995         i=l+1
 996         j=ir
 997         indxt=indx(l)
 998         a=arr(indxt)
 999 3       continue
1000           i=i+1
1001         if(arr(indx(i)).lt.a)goto 3
1002 4       continue
1003           j=j-1
1004         if(arr(indx(j)).gt.a)goto 4
1005         if(j.lt.i)goto 5
1006         itemp=indx(i)
1007         indx(i)=indx(j)
1008         indx(j)=itemp
1009         goto 3
1010 5       indx(l)=indx(j)
1011         indx(j)=indxt
1012         jstack=jstack+2
1013         if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
1014         if(ir-i+1.ge.j-l)then
1015           istack(jstack)=ir
1016           istack(jstack-1)=i
1017           ir=j-1
1018         else
1019           istack(jstack)=j-1
1020           istack(jstack-1)=l
1021           l=i
1022         endif
1023       endif
1024       goto 1
1025       end subroutine indexx
1026 !  (C) Copr. 1986-92 Numerical Recipes Software *11915aZ%.
1027 !-----------------------------------------------------------------------------
1028 ! minim_jlee.F
1029 !-----------------------------------------------------------------------------
1030       subroutine minim_jlee
1031
1032       use minim_data
1033       use MPI_data
1034       use energy_data
1035       use compare_data
1036       use control_data
1037       use geometry_data, only: nvar,nphi
1038       use geometry, only:dist
1039       use energy, only:fdum
1040       use control, only:init_int_table
1041       use minimm, only:sumsl,deflt
1042 !  controls minimization and sorting routines
1043 !      implicit real*8 (a-h,o-z)
1044 !      include 'DIMENSIONS'
1045 !      include 'COMMON.VAR'
1046 !      include 'COMMON.IOUNITS'
1047 !      include 'COMMON.MINIM'
1048 !      include 'COMMON.CONTROL'
1049       include 'mpif.h'
1050       integer,parameter :: liv=60
1051       integer :: lv
1052 !      external func,gradient!,fdum     !use minim & energy
1053 !      real(kind=4) :: ran1,ran2,ran3
1054 !      include 'COMMON.SETUP'
1055 !      include 'COMMON.GEO'
1056 !      include 'COMMON.FFIELD'
1057 !      include 'COMMON.SBRIDGE'
1058 !      include 'COMMON.DISTFIT'
1059 !      include 'COMMON.CHAIN'
1060       integer,dimension(mpi_status_size) :: muster
1061       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
1062       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
1063       real(kind=8),dimension(6*nres) :: var2    !(maxvar) (maxvar=6*maxres)
1064       integer,dimension(nres) :: iffr   !(maxres)
1065       integer,dimension((nres-1)*(nres-2)/2) :: ihpbt,jhpbt     !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2)
1066       real(kind=8),dimension(6*nres) :: d,garbage       !(maxvar) (maxvar=6*maxres)
1067 !el      real(kind=8),dimension(1:lv+1) :: v
1068       real(kind=8) :: energia(0:n_ene),time0s,time1s
1069       integer,dimension(9) :: indx
1070       integer,dimension(12) :: info
1071       integer,dimension(liv) :: iv
1072       integer :: idum(1)
1073       real(kind=8) :: rdum(1)
1074       integer,dimension(2,12*nres) :: icont_    !(2,maxcont)(maxcont=12*maxres)
1075       logical :: fail   !check_var,
1076       integer :: iloop(2)
1077 !el      common /przechowalnia/ v
1078       integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode
1079       real(kind=8) :: rad,eee,etot      !,fdum
1080 !el  from subroutine parmread
1081 ! Define the constants of the disulfide bridge
1082 ! Old arbitrary potential
1083       real(kind=8),parameter :: dbr=4.20D0
1084       real(kind=8),parameter :: fbr=3.30D0
1085 !-----------------
1086       lv=77+(6*nres)*(6*nres+17)/2      !77+maxvar*(maxvar+17)/2 (maxvar=6*maxres)
1087       data rad /1.745329252d-2/
1088 !  receive # of start
1089 !      print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
1090 !     &   ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
1091       if (.not. allocated(v)) allocate(v(1:lv))
1092       nhpb0=nhpb
1093    10 continue
1094       time0s=MPI_WTIME()
1095 !      print *, 'MINIM_JLEE: ',me,' is waiting'
1096       call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,&
1097                     muster,ierr)
1098       time1s=MPI_WTIME()
1099       write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
1100       call flush(iout)
1101        n=info(1)
1102 !      print *, 'MINIM_JLEE: ',me,' received: ',n
1103
1104 !rc      if (ierr.ne.0) go to 100
1105 !  if # = 0, return
1106       if (n.eq.0) then
1107         write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
1108         call flush(iout)
1109         return
1110       endif
1111
1112       nfun=0
1113       IF (n.lt.0) THEN
1114        call mpi_recv(var,nvar,mpi_double_precision,&
1115                     king,idreal,CG_COMM,muster,ierr)
1116        call mpi_recv(iffr,nres,mpi_integer,&
1117                     king,idint,CG_COMM,muster,ierr)
1118        call mpi_recv(var2,nvar,mpi_double_precision,&
1119                     king,idreal,CG_COMM,muster,ierr)
1120       ELSE
1121 !  receive initial values of variables
1122        call mpi_recv(var,nvar,mpi_double_precision,&
1123                     king,idreal,CG_COMM,muster,ierr)
1124 !rc       if (ierr.ne.0) go to 100
1125       ENDIF
1126
1127       if(vdisulf.and.info(2).ne.-1) then
1128         if(info(4).ne.0)then
1129            call mpi_recv(ihpbt,info(4),mpi_integer,&
1130                     king,idint,CG_COMM,muster,ierr)
1131            call mpi_recv(jhpbt,info(4),mpi_integer,&
1132                     king,idint,CG_COMM,muster,ierr)
1133         endif
1134       endif
1135
1136       IF (n.lt.0) THEN
1137         n=-n
1138         nhpb=nhpb0
1139         link_start=1
1140         link_end=nhpb
1141         call init_int_table
1142         call contact_cp(var,var2,iffr,nfun,n)
1143       ENDIF
1144
1145       if(vdisulf.and.info(2).ne.-1) then
1146         nss=0
1147         if(info(4).ne.0)then
1148 !d          write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
1149           call var_to_geom(nvar,var)
1150           call chainbuild
1151           do i=1,info(4)
1152            if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
1153             nss=nss+1
1154             ihpb(nss)=ihpbt(i)
1155             jhpb(nss)=jhpbt(i)
1156 !d            write(iout,*) 'SS  mv=',info(3),
1157 !d     &         ihpb(nss)-nres,jhpb(nss)-nres,
1158 !d     &         dist(ihpb(nss),jhpb(nss))
1159             dhpb(nss)=dbr
1160             forcon(nss)=fbr
1161            else
1162 !d            write(iout,*) 'rm SS  mv=',info(3),
1163 !d     &         ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
1164            endif
1165           enddo
1166         endif
1167         nhpb=nss
1168         link_start=1
1169         link_end=nhpb
1170         call init_int_table
1171       endif
1172
1173       if (info(3).eq.14) then
1174        write(iout,*) 'calling local_move',info(7),info(8)
1175        call local_move_init(.false.)
1176        call var_to_geom(nvar,var)
1177        call local_move(info(7),info(8),20d0,50d0)
1178        call geom_to_var(nvar,var)
1179       endif
1180
1181
1182       if (info(3).eq.16) then
1183        write(iout,*) 'calling beta_slide',info(7),info(8),&
1184                   info(10), info(11), info(12)
1185        call var_to_geom(nvar,var)
1186        call beta_slide(info(7),info(8),info(10),info(11),info(12), &
1187                                                            nfun,n)
1188        call geom_to_var(nvar,var)
1189       endif
1190
1191
1192       if (info(3).eq.17) then
1193        write(iout,*) 'calling beta_zip',info(7),info(8)
1194        call var_to_geom(nvar,var)
1195        call beta_zip(info(7),info(8),nfun,n)
1196        call geom_to_var(nvar,var)
1197       endif
1198
1199
1200 !rc overlap test
1201
1202       if (overlapsc) then
1203
1204           call var_to_geom(nvar,var)
1205           call chainbuild
1206           call etotal(energia)
1207           nfun=nfun+1
1208           if (energia(1).eq.1.0d20) then
1209             info(3)=-info(3)
1210             write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
1211             call overlap_sc(fail)
1212             if(.not.fail) then
1213              call geom_to_var(nvar,var)
1214              call etotal(energia)
1215              nfun=nfun+1
1216              write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
1217             else
1218              v(10)=1.0d20
1219              iv(1)=-1
1220              goto 201
1221             endif
1222           endif
1223       endif
1224
1225       if (searchsc) then
1226           call var_to_geom(nvar,var)
1227           call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1228           call geom_to_var(nvar,var)
1229 !d          write(iout,*) 'sc_move',nft_sc,etot
1230       endif
1231
1232       if (check_var(var,info)) then
1233            v(10)=1.0d21
1234            iv(1)=6
1235            goto 201
1236       endif
1237
1238
1239 !rc
1240
1241 !      write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
1242 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
1243 !      write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
1244 !      write (*,'(8f10.4)') (var(i),i=1,nvar)
1245
1246        do i=1,nvar
1247          garbage(i)=var(i)
1248        enddo
1249
1250       call deflt(2,iv,liv,lv,v)
1251 ! 12 means fresh start, dont call deflt
1252       iv(1)=12
1253 ! max num of fun calls
1254       if (maxfun.eq.0) maxfun=500
1255       iv(17)=maxfun
1256 ! max num of iterations
1257       if (maxmin.eq.0) maxmin=1000
1258       iv(18)=maxmin
1259 ! controls output
1260       iv(19)=2
1261 ! selects output unit
1262 !d      iv(21)=iout
1263       iv(21)=0
1264 ! 1 means to print out result
1265       iv(22)=0
1266 !d        iv(22)=1
1267 ! 1 means to print out summary stats
1268       iv(23)=0
1269 ! 1 means to print initial x and d
1270       iv(24)=0
1271
1272 !      if(me.eq.3.and.n.eq.255) then
1273 !       print *,' CHUJ: stoi'
1274 !       iv(21)=6
1275 !       iv(22)=1
1276 !       iv(23)=1
1277 !       iv(24)=1
1278 !      endif
1279
1280 ! min val for v(radfac) default is 0.1
1281       v(24)=0.1D0
1282 ! max val for v(radfac) default is 4.0
1283       v(25)=2.0D0
1284 !     v(25)=4.0D0
1285 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
1286 ! the sumsl default is 0.1
1287       v(26)=0.1D0
1288 ! false conv if (act fnctn decrease) .lt. v(34)
1289 ! the sumsl default is 100*machep
1290       v(34)=v(34)/100.0D0
1291 ! absolute convergence
1292       if (tolf.eq.0.0D0) tolf=1.0D-4
1293       v(31)=tolf
1294 ! relative convergence
1295       if (rtolf.eq.0.0D0) rtolf=1.0D-4
1296       v(32)=rtolf
1297 ! controls initial step size
1298        v(35)=1.0D-1
1299 ! large vals of d correspond to small components of step
1300       do i=1,nphi
1301         d(i)=1.0D-1
1302       enddo
1303       do i=nphi+1,nvar
1304         d(i)=1.0D-1
1305       enddo
1306 !  minimize energy
1307 !      write (iout,*) 'Processor',me,' nvar',nvar
1308 !      write (iout,*) 'Variables BEFORE minimization:'
1309 !      write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1310
1311 !      print *, 'MINIM_JLEE: ',me,' before SUMSL '
1312
1313       call func(nvar,var,nf,eee,idum,rdum,fdum)
1314       nfun=nfun+1
1315       if(eee.ge.1.0d20) then
1316 !       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
1317 !       print *,' energy before SUMSL =',eee
1318 !       print *,' aborting local minimization'
1319        iv(1)=-1
1320        v(10)=eee
1321        go to 201
1322       endif
1323
1324 !t      time0s=MPI_WTIME()
1325       call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
1326 !t      write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
1327 !      print *, 'MINIM_JLEE: ',me,' after SUMSL '
1328
1329 !  find which conformation was returned from sumsl
1330         nfun=nfun+iv(7)
1331 !      print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
1332 !     & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
1333 !        if (iv(1).ne.4 .or. nf.le.1) then
1334 !          write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
1335 !         write (*,*) 'Initial Variables'
1336 !         write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
1337 !         write (*,*) 'Variables'
1338 !         write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1339 !         write (*,*) 'Vector d'
1340 !         write (*,'(8f10.4)') (d(i),i=1,nvar)
1341 !         write (iout,*) 'Processor',me,' something bad in SUMSL',
1342 !    &       iv(1),nf
1343 !         write (iout,*) 'Initial Variables'
1344 !         write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
1345 !         write (iout,*) 'Variables'
1346 !         write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1347 !         write (iout,*) 'Vector d'
1348 !         write (iout,'(8f10.4)') (d(i),i=1,nvar)
1349 !        endif
1350 !        if (nf.lt.iv(6)-1) then
1351 !  recalculate intra- and interchain energies
1352 !         call func(nvar,var,nf,v(10),iv,v,fdum)
1353 !        else if (nf.eq.iv(6)-1) then
1354 !  regenerate conformation
1355 !         call var_to_geom(nvar,var)
1356 !         call chainbuild
1357 !        endif
1358 !  change origin and axes to standard ECEPP format
1359 !        call var_to_geom(nvar,var)
1360 !      write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
1361 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
1362 !      write (iout,*) 'Energy:',v(10)
1363 !  send back output
1364 !       print *, 'MINIM_JLEE: ',me,' minimized: ',n
1365   201  continue
1366         indx(1)=n
1367 ! return code: 6-gradient 9-number of ftn evaluation, etc
1368         indx(2)=iv(1)
1369 ! total # of ftn evaluations (for iwf=0, it includes all minimizations).
1370         indx(3)=nfun
1371         indx(4)=info(2)
1372         indx(5)=info(3)
1373         indx(6)=nss
1374         indx(7)=info(5)
1375         indx(8)=info(6)
1376         indx(9)=info(9)
1377         call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,&
1378                        ierr)
1379 !  send back energies
1380 ! al & cc
1381 ! calculate contact order
1382 #ifdef CO_BIAS
1383         call contact(.false.,ncont,icont_,co)
1384         erg(1)=v(10)-1.0d2*co
1385 #else
1386         erg(1)=v(10)
1387 #endif
1388         j=1
1389         call mpi_send(erg,j,mpi_double_precision,king,idreal,&
1390                        CG_COMM,ierr)
1391 #ifdef CO_BIAS
1392         call mpi_send(co,j,mpi_double_precision,king,idreal,&
1393                        CG_COMM,ierr)
1394 #endif
1395 !  send back values of variables
1396         call mpi_send(var,nvar,mpi_double_precision,&
1397                      king,idreal,CG_COMM,ierr)
1398 !        print * , 'MINIM_JLEE: Processor',me,' send erg and var '
1399
1400         if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
1401 !d          call intout
1402 !d          call chainbuild
1403 !d          call etotal(energia(0))
1404 !d          etot=energia(0)
1405 !d          call enerprint(energia(0))
1406          call mpi_send(ihpb,nss,mpi_integer,&
1407                      king,idint,CG_COMM,ierr)
1408          call mpi_send(jhpb,nss,mpi_integer,&
1409                      king,idint,CG_COMM,ierr)
1410         endif
1411
1412         go to 10
1413   100 print *, ' error in receiving message from emperor', me
1414       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1415       return
1416   200 print *, ' error in sending message to emperor'
1417       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1418       return
1419   300 print *, ' error in communicating with emperor'
1420       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1421       return
1422   956 format (' initial energy could not be calculated',41x)
1423   957 format (80x)
1424   965 format (' convergence code ',i2,' # of function calls ',&
1425         i4,' # of gradient calls ',i4,10x)
1426   975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
1427       end subroutine minim_jlee
1428 !-----------------------------------------------------------------------------
1429 ! newconf.f
1430 !-----------------------------------------------------------------------------
1431       subroutine make_var(n,idum,iter_csa)
1432
1433       use MD_data
1434       use energy_data
1435       use compare_data
1436       use control_data, only: vdisulf
1437       use geometry_data
1438       use geometry, only: dist
1439       include 'mpif.h'
1440
1441 !      use random, only: iran_num,ran_number
1442 !      implicit real*8 (a-h,o-z)
1443 !      include 'DIMENSIONS'
1444 !      include 'COMMON.IOUNITS'
1445 !      include 'COMMON.CSA'
1446 !      include 'COMMON.BANK'
1447 !      include 'COMMON.CHAIN'
1448 !      include 'COMMON.INTERACT'
1449 !      include 'COMMON.HAIRPIN'
1450 !      include 'COMMON.VAR'
1451 !      include 'COMMON.DISTFIT'
1452 !      include 'COMMON.GEO'
1453 !      include 'COMMON.CONTROL'
1454       logical :: nicht_getan,nicht_getan1,fail,lfound
1455       integer :: nharp,iharp(4,nres/3),nconf_harp
1456       integer :: iisucc(mxio)
1457       logical :: ifused(mxio)
1458       integer :: nhx_seed(nseed),ihx_seed(4,nres/3,nseed) !max_seed
1459       integer :: nhx_use(nseed),ihx_use(0:4,nres/3,nseed)
1460       integer :: nlx_seed(nseed),ilx_seed(2,nres/3,nseed),&
1461                nlx_use(nseed),ilx_use(nres/3,nseed)
1462 !      real(kind=4) :: ran1,ran2
1463
1464       integer :: i,j,k,n,idum,iter_csa,iran,index,n7frag,n8frag,n14frag,&
1465             n15frag,nbefrag,nlx_tot,iters,i1,i2,i3,ntot_gen,ngen,iih,&
1466             ij,jr,iim,nhx_tot,idummy,iter,iif,iig,icheck,ishift,iang,&
1467             n8c,ih_start,ih_end,n7c,index2,isize,nsucc,nacc,j1,nran,&
1468             ierror,ierrcode
1469       real(kind=8) :: d
1470
1471       write (iout,*) 'make_var : nseed=',nseed,'ntry=',n
1472       index=0
1473
1474 !-----------------------------------------
1475       if (n7.gt.0.or.n8.gt.0.or.n9.gt.0.or.n14.gt.0.or.n15.gt.0 &
1476           .or.n16.gt.0.or.n17.gt.0.or.n18.gt.0) &
1477                  call select_frag(n7frag,n8frag,n14frag,&
1478                  n15frag,nbefrag,iter_csa)
1479
1480 !---------------------------------------------------
1481 ! N18 - random perturbation of one phi(=gamma) angle in a loop
1482 !
1483       IF (n18.gt.0) THEN
1484       nlx_tot=0
1485       do iters=1,nseed
1486         i1=is(iters)
1487         nlx_seed(iters)=0
1488         do i2=1,n14frag
1489           if (lvar_frag(i2,1).eq.i1) then
1490             nlx_seed(iters)=nlx_seed(iters)+5
1491             ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
1492             ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
1493             ilx_use(nlx_seed(iters),iters)=5
1494           endif
1495         enddo
1496         nlx_use(iters)=nlx_seed(iters)
1497         nlx_tot=nlx_tot+nlx_seed(iters)
1498       enddo
1499
1500       if (nlx_tot .ge. n18*nseed) then
1501         ntot_gen=n18*nseed
1502       else
1503         ntot_gen=(nlx_tot/nseed)*nseed
1504       endif
1505
1506       ngen=0
1507       do while (ngen.lt.ntot_gen)
1508       do iters=1,nseed
1509         iseed=is(iters)
1510         if (nlx_use(iters).gt.0) then
1511           nicht_getan=.true.
1512           do while (nicht_getan)
1513             iih=iran_num(1,nlx_seed(iters))
1514             if (ilx_use(iih,iters).gt.0) then
1515               nicht_getan=.false.
1516               ilx_use(iih,iters)=ilx_use(iih,iters)-1
1517               nlx_use(iters)=nlx_use(iters)-1
1518             endif
1519           enddo
1520           ngen=ngen+1
1521           index=index+1
1522           movenx(index)=18
1523           parent(1,index)=iseed
1524           parent(2,index)=0
1525
1526
1527            if (vdisulf) then
1528              nss_in(index)=bvar_nss(iseed)
1529              do ij=1,nss_in(index)
1530                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1531                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1532              enddo
1533            endif
1534
1535
1536           do k=1,numch
1537             do j=2,nres-1
1538               do i=1,4
1539                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1540               enddo
1541             enddo
1542           enddo
1543
1544           jr=iran_num(ilx_seed(1,iih,iters),ilx_seed(2,iih,iters))
1545           d=ran_number(-pi,pi)
1546           dihang_in(2,jr-2,1,index)=pinorm(dihang_in(2,jr-2,1,index)+d)
1547
1548
1549           if (ngen.eq.ntot_gen) goto 145
1550         endif
1551       enddo
1552       enddo
1553   145 continue
1554
1555       ENDIF
1556
1557
1558 !-----------------------------------------
1559 ! N17 : zip a beta in a seed by forcing one additional p-p contact
1560 !
1561       IF (n17.gt.0) THEN
1562       nhx_tot=0
1563       do iters=1,nseed
1564         i1=is(iters)
1565         nhx_seed(iters)=0
1566         nhx_use(iters)=0
1567         do i2=1,nbefrag
1568           if (avar_frag(i2,1).eq.i1) then
1569             nhx_seed(iters)=nhx_seed(iters)+1
1570             ihx_use(2,nhx_seed(iters),iters)=1
1571             if (avar_frag(i2,5)-avar_frag(i2,3).le.3.and. &
1572                  avar_frag(i2,2).gt.1.and.avar_frag(i2,4).lt.nres) then
1573              ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1574              ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
1575              ihx_use(0,nhx_seed(iters),iters)=1
1576              ihx_use(1,nhx_seed(iters),iters)=0
1577              nhx_use(iters)=nhx_use(iters)+1
1578             else
1579               if (avar_frag(i2,4).gt.avar_frag(i2,5)) then
1580                if (avar_frag(i2,2).gt.1.and. &
1581                            avar_frag(i2,4).lt.nres) then
1582                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1583                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
1584                 ihx_use(0,nhx_seed(iters),iters)=1
1585                 ihx_use(1,nhx_seed(iters),iters)=0
1586                 nhx_use(iters)=nhx_use(iters)+1
1587                endif
1588                if (avar_frag(i2,3).lt.nres.and. &
1589                            avar_frag(i2,5).gt.1) then
1590                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
1591                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)-1
1592                 ihx_use(0,nhx_seed(iters),iters)= &
1593                           ihx_use(0,nhx_seed(iters),iters)+1
1594                 ihx_use(2,nhx_seed(iters),iters)=0
1595                 nhx_use(iters)=nhx_use(iters)+1
1596                endif
1597               else
1598                if (avar_frag(i2,2).gt.1.and. &
1599                            avar_frag(i2,4).gt.1) then
1600                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1601                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)-1
1602                 ihx_use(0,nhx_seed(iters),iters)=1
1603                 ihx_use(1,nhx_seed(iters),iters)=0
1604                 nhx_use(iters)=nhx_use(iters)+1
1605                endif
1606                if (avar_frag(i2,3).lt.nres.and. &
1607                            avar_frag(i2,5).lt.nres) then
1608                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
1609                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)+1
1610                 ihx_use(0,nhx_seed(iters),iters)= &
1611                           ihx_use(0,nhx_seed(iters),iters)+1
1612                 ihx_use(2,nhx_seed(iters),iters)=0
1613                 nhx_use(iters)=nhx_use(iters)+1
1614                endif
1615               endif
1616             endif
1617           endif
1618         enddo
1619
1620         nhx_tot=nhx_tot+nhx_use(iters)
1621 !d        write (iout,*) "debug N17",iters,nhx_seed(iters),
1622 !d     &                     nhx_use(iters),nhx_tot
1623       enddo
1624
1625       if (nhx_tot .ge. n17*nseed) then
1626         ntot_gen=n17*nseed
1627       else if (nhx_tot .ge. nseed) then
1628         ntot_gen=(nhx_tot/nseed)*nseed
1629       else
1630         ntot_gen=nhx_tot
1631       endif
1632 !d      write (iout,*) "debug N17==",ntot_gen,nhx_tot,nseed
1633
1634       ngen=0
1635       do while (ngen.lt.ntot_gen)
1636       do iters=1,nseed
1637         iseed=is(iters)
1638         if (nhx_use(iters).gt.0) then
1639 !d          write (iout,*) "debug N17",nhx_use(iters),ngen,ntot_gen
1640 !d          write (iout,*) "debugN17^",
1641 !d     &                    (ihx_use(0,k,iters),k=1,nhx_use(iters))
1642           nicht_getan=.true.
1643           do while (nicht_getan)
1644             iih=iran_num(1,nhx_seed(iters))
1645 !d            write (iout,*) "debugN17^",iih
1646             if (ihx_use(0,iih,iters).gt.0) then
1647               iim=iran_num(1,2)
1648 !d                write (iout,*) "debugN17=",iih,nhx_seed(iters)
1649 !d                write (iout,*) "debugN17-",iim,'##',
1650 !d     &                           (ihx_use(k,iih,iters),k=0,2)
1651 !d                call flush(iout)
1652               do while (ihx_use(iim,iih,iters).eq.1)
1653                 iim=iran_num(1,2)
1654 !d                write (iout,*) "debugN17-",iim,'##',
1655 !d     &                           (ihx_use(k,iih,iters),k=0,2)
1656 !d                call flush(iout)
1657               enddo
1658               nicht_getan=.false.
1659               ihx_use(iim,iih,iters)=1
1660               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
1661               nhx_use(iters)=nhx_use(iters)-1
1662             endif
1663           enddo
1664           ngen=ngen+1
1665           index=index+1
1666           movenx(index)=17
1667           parent(1,index)=iseed
1668           parent(2,index)=0
1669
1670            if (vdisulf) then
1671              nss_in(index)=bvar_nss(iseed)
1672              do ij=1,nss_in(index)
1673                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1674                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1675              enddo
1676            endif
1677
1678           do k=1,numch
1679             do j=2,nres-1
1680               do i=1,4
1681                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1682               enddo
1683             enddo
1684           enddo
1685
1686           if (iim.eq.1) then
1687             idata(1,index)=ihx_seed(1,iih,iters)
1688             idata(2,index)=ihx_seed(2,iih,iters)
1689           else
1690             idata(1,index)=ihx_seed(3,iih,iters)
1691             idata(2,index)=ihx_seed(4,iih,iters)
1692           endif
1693
1694           if (ngen.eq.ntot_gen) goto 115
1695         endif
1696       enddo
1697       enddo
1698   115 continue
1699       write (iout,*) "N17",n17," ngen/nseed",ngen/nseed,&
1700                                  ngen,nseed
1701
1702
1703       ENDIF
1704 !-----------------------------------------
1705 ! N16 : slide non local beta in a seed by +/- 1 or +/- 2
1706 !
1707       IF (n16.gt.0) THEN
1708       nhx_tot=0
1709       do iters=1,nseed
1710         i1=is(iters)
1711         nhx_seed(iters)=0
1712         do i2=1,n7frag
1713           if (bvar_frag(i2,1).eq.i1) then
1714             nhx_seed(iters)=nhx_seed(iters)+1
1715             ihx_seed(1,nhx_seed(iters),iters)=bvar_frag(i2,3)
1716             ihx_seed(2,nhx_seed(iters),iters)=bvar_frag(i2,4)
1717             ihx_seed(3,nhx_seed(iters),iters)=bvar_frag(i2,5)
1718             ihx_seed(4,nhx_seed(iters),iters)=bvar_frag(i2,6)
1719             ihx_use(0,nhx_seed(iters),iters)=4
1720             do i3=1,4
1721               ihx_use(i3,nhx_seed(iters),iters)=0
1722             enddo
1723           endif
1724         enddo
1725         nhx_use(iters)=4*nhx_seed(iters)
1726         nhx_tot=nhx_tot+nhx_seed(iters)
1727 !d        write (iout,*) "debug N16",iters,nhx_seed(iters)
1728       enddo
1729
1730       if (4*nhx_tot .ge. n16*nseed) then
1731         ntot_gen=n16*nseed
1732       else if (4*nhx_tot .ge. nseed) then
1733         ntot_gen=(4*nhx_tot/nseed)*nseed
1734       else
1735         ntot_gen=4*nhx_tot
1736       endif
1737       write (iout,*) "debug N16",ntot_gen,4*nhx_tot,nseed
1738
1739       ngen=0
1740       do while (ngen.lt.ntot_gen)
1741       do iters=1,nseed
1742         iseed=is(iters)
1743         if (nhx_use(iters).gt.0) then
1744           nicht_getan=.true.
1745           do while (nicht_getan)
1746             iih=iran_num(1,nhx_seed(iters))
1747             if (ihx_use(0,iih,iters).gt.0) then
1748               iim=iran_num(1,4)
1749               do while (ihx_use(iim,iih,iters).eq.1)
1750 !d                write (iout,*) iim,
1751 !d     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
1752                 iim=iran_num(1,4)
1753               enddo
1754               nicht_getan=.false.
1755               ihx_use(iim,iih,iters)=1
1756               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
1757               nhx_use(iters)=nhx_use(iters)-1
1758             endif
1759           enddo
1760           ngen=ngen+1
1761           index=index+1
1762           movenx(index)=16
1763           parent(1,index)=iseed
1764           parent(2,index)=0
1765
1766            if (vdisulf) then
1767              nss_in(index)=bvar_nss(iseed)
1768              do ij=1,nss_in(index)
1769                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1770                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1771              enddo
1772            endif
1773
1774           do k=1,numch
1775             do j=2,nres-1
1776               do i=1,4
1777                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1778               enddo
1779             enddo
1780           enddo
1781
1782           do i=1,4
1783            idata(i,index)=ihx_seed(i,iih,iters)
1784           enddo
1785           idata(5,index)=iim
1786
1787           if (ngen.eq.ntot_gen) goto 116
1788         endif
1789       enddo
1790       enddo
1791   116 continue
1792       write (iout,*) "N16",n16," ngen/nseed",ngen/nseed,&
1793                                  ngen,nseed
1794       ENDIF
1795 !-----------------------------------------
1796 ! N15 : copy two 2nd structure elements from 1 or 2 conf. in bank to a seed
1797 !
1798       IF (n15.gt.0) THEN
1799
1800        do iters=1,nseed
1801          iseed=is(iters)
1802          do i=1,mxio
1803            ifused(i)=.false.
1804          enddo
1805
1806          do idummy=1,n15
1807            iter=0
1808    84      continue
1809
1810            iran=0
1811            iif=iran_num(1,n15frag)
1812            do while( (ifused(iif) .or. svar_frag(iif,1).eq.iseed) .and. &
1813                             iran.le.mxio )
1814             iif=iran_num(1,n15frag)
1815             iran=iran+1
1816            enddo
1817            if(iran.ge.mxio) goto 811
1818
1819            iran=0
1820            iig=iran_num(1,n15frag)
1821            do while( (ifused(iig) .or. svar_frag(iig,1).eq.iseed  .or. &
1822                    .not.(svar_frag(iif,3).lt.svar_frag(iig,2).or. &
1823                          svar_frag(iig,3).lt.svar_frag(iif,2)) ) .and. &
1824                             iran.le.mxio )
1825             iig=iran_num(1,n15frag)
1826             iran=iran+1
1827            enddo
1828            if(iran.ge.mxio) goto 811
1829
1830            index=index+1
1831            movenx(index)=15
1832            parent(1,index)=iseed
1833            parent(2,index)=svar_frag(iif,1)
1834            parent(3,index)=svar_frag(iig,1)
1835
1836
1837            if (vdisulf) then
1838              nss_in(index)=bvar_nss(iseed)
1839              do ij=1,nss_in(index)
1840                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1841                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1842              enddo
1843            endif
1844
1845            ifused(iif)=.true.
1846            ifused(iig)=.true.
1847            call newconf_copy(idum,dihang_in(1,1,1,index),&
1848                    svar_frag(iif,1),svar_frag(iif,2),svar_frag(iif,3))
1849
1850            do j=svar_frag(iig,2),svar_frag(iig,3)
1851              do i=1,4
1852               dihang_in(i,j,1,index)=bvar(i,j,1,svar_frag(iig,1))
1853              enddo
1854            enddo
1855
1856
1857            if(iter.lt.10) then
1858             call check_old(icheck,index)
1859             if(icheck.eq.1) then
1860               index=index-1
1861               ifused(iif)=.false.
1862               goto 84
1863             endif
1864            endif
1865
1866   811     continue
1867          enddo
1868        enddo
1869       ENDIF
1870
1871 !-----------------------------------------
1872 ! N14 local_move (Maurizio) for loops in a seed
1873 !
1874       IF (n14.gt.0) THEN
1875       nlx_tot=0
1876       do iters=1,nseed
1877         i1=is(iters)
1878         nlx_seed(iters)=0
1879         do i2=1,n14frag
1880           if (lvar_frag(i2,1).eq.i1) then
1881             nlx_seed(iters)=nlx_seed(iters)+3
1882             ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
1883             ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
1884             ilx_use(nlx_seed(iters),iters)=3
1885           endif
1886         enddo
1887         nlx_use(iters)=nlx_seed(iters)
1888         nlx_tot=nlx_tot+nlx_seed(iters)
1889 !d        write (iout,*) "debug N14",iters,nlx_seed(iters)
1890       enddo
1891
1892       if (nlx_tot .ge. n14*nseed) then
1893         ntot_gen=n14*nseed
1894       else
1895         ntot_gen=(nlx_tot/nseed)*nseed
1896       endif
1897 !d      write (iout,*) "debug N14",ntot_gen,n14frag,nseed
1898
1899       ngen=0
1900       do while (ngen.lt.ntot_gen)
1901       do iters=1,nseed
1902         iseed=is(iters)
1903         if (nlx_use(iters).gt.0) then
1904           nicht_getan=.true.
1905           do while (nicht_getan)
1906             iih=iran_num(1,nlx_seed(iters))
1907             if (ilx_use(iih,iters).gt.0) then
1908               nicht_getan=.false.
1909               ilx_use(iih,iters)=ilx_use(iih,iters)-1
1910               nlx_use(iters)=nlx_use(iters)-1
1911             endif
1912           enddo
1913           ngen=ngen+1
1914           index=index+1
1915           movenx(index)=14
1916           parent(1,index)=iseed
1917           parent(2,index)=0
1918
1919           idata(1,index)=ilx_seed(1,iih,iters)
1920           idata(2,index)=ilx_seed(2,iih,iters)
1921
1922
1923            if (vdisulf) then
1924              nss_in(index)=bvar_nss(iseed)
1925              do ij=1,nss_in(index)
1926                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1927                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1928              enddo
1929            endif
1930
1931
1932           do k=1,numch
1933             do j=2,nres-1
1934               do i=1,4
1935                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1936               enddo
1937             enddo
1938           enddo
1939
1940           if (ngen.eq.ntot_gen) goto 131
1941         endif
1942       enddo
1943       enddo
1944   131 continue
1945 !d      write (iout,*) "N14",n14," ngen/nseed",ngen/nseed,
1946 !d     &                           ngen,nseed
1947
1948       ENDIF
1949 !-----------------------------------------
1950 ! N9 : shift a helix in a seed
1951 !
1952       IF (n9.gt.0) THEN
1953       nhx_tot=0
1954       do iters=1,nseed
1955         i1=is(iters)
1956         nhx_seed(iters)=0
1957         do i2=1,n8frag
1958           if (hvar_frag(i2,1).eq.i1) then
1959             nhx_seed(iters)=nhx_seed(iters)+1
1960             ihx_seed(1,nhx_seed(iters),iters)=hvar_frag(i2,2)
1961             ihx_seed(2,nhx_seed(iters),iters)=hvar_frag(i2,3)
1962             ihx_use(0,nhx_seed(iters),iters)=4
1963             do i3=1,4
1964               ihx_use(i3,nhx_seed(iters),iters)=0
1965             enddo
1966           endif
1967         enddo
1968         nhx_use(iters)=4*nhx_seed(iters)
1969         nhx_tot=nhx_tot+nhx_seed(iters)
1970 !d        write (iout,*) "debug N9",iters,nhx_seed(iters)
1971       enddo
1972
1973       if (4*nhx_tot .ge. n9*nseed) then
1974         ntot_gen=n9*nseed
1975       else
1976         ntot_gen=(4*nhx_tot/nseed)*nseed
1977       endif
1978 !d      write (iout,*) "debug N9",ntot_gen,n8frag,nseed
1979
1980       ngen=0
1981       do while (ngen.lt.ntot_gen)
1982       do iters=1,nseed
1983         iseed=is(iters)
1984         if (nhx_use(iters).gt.0) then
1985           nicht_getan=.true.
1986           do while (nicht_getan)
1987             iih=iran_num(1,nhx_seed(iters))
1988             if (ihx_use(0,iih,iters).gt.0) then
1989               iim=iran_num(1,4)
1990               do while (ihx_use(iim,iih,iters).eq.1)
1991 !d                write (iout,*) iim,
1992 !d     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
1993                 iim=iran_num(1,4)
1994               enddo
1995               nicht_getan=.false.
1996               ihx_use(iim,iih,iters)=1
1997               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1
1998               nhx_use(iters)=nhx_use(iters)-1
1999             endif
2000           enddo
2001           ngen=ngen+1
2002           index=index+1
2003           movenx(index)=9
2004           parent(1,index)=iseed
2005           parent(2,index)=0
2006
2007            if (vdisulf) then
2008              nss_in(index)=bvar_nss(iseed)
2009              do ij=1,nss_in(index)
2010                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2011                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2012              enddo
2013            endif
2014
2015           do k=1,numch
2016             do j=2,nres-1
2017               do i=1,4
2018                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2019               enddo
2020             enddo
2021           enddo
2022
2023           jstart=max(nnt,ihx_seed(1,iih,iters)+1)
2024           jend=min(nct,ihx_seed(2,iih,iters))
2025 !d          write (iout,*) "debug N9",iters,iih,jstart,jend
2026           if (iim.eq.1) then
2027             ishift=-2
2028           else if (iim.eq.2) then
2029             ishift=-1
2030           else if (iim.eq.3) then
2031             ishift=1
2032           else if (iim.eq.4) then
2033             ishift=2
2034           else
2035             write (iout,*) 'CHUJ NASTAPIL: iim=',iim
2036 #ifdef MPI !el
2037             call mpi_abort(mpi_comm_world,ierror,ierrcode)
2038 #endif
2039           endif
2040           do j=jstart,jend
2041             if (itype(j).eq.10) then
2042               iang=2
2043             else
2044               iang=4
2045             endif
2046             do i=1,iang
2047               if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
2048                 dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
2049             enddo
2050           enddo
2051           if (ishift.gt.0) then
2052            do j=0,ishift-1
2053             if (itype(jend+j).eq.10) then
2054               iang=2
2055             else
2056               iang=4
2057             endif
2058             do i=1,iang
2059               if (jend+j.ge.nnt.and.jend+j.le.nct) &
2060                 dihang_in(i,jstart+j,1,index)=bvar(i,jend+j,1,iseed)
2061             enddo
2062            enddo
2063           else
2064            do j=0,-ishift-1
2065             if (itype(jstart+j).eq.10) then
2066               iang=2
2067             else
2068               iang=4
2069             endif
2070             do i=1,iang
2071               if (jend+j.ge.nnt.and.jend+j.le.nct) &
2072                 dihang_in(i,jend+j,1,index)=bvar(i,jstart+j,1,iseed)
2073             enddo
2074            enddo
2075           endif
2076           if (ngen.eq.ntot_gen) goto 133
2077         endif
2078       enddo
2079       enddo
2080   133 continue
2081 !d      write (iout,*) "N9",n9," ngen/nseed",ngen/nseed,
2082 !d     &                           ngen,nseed
2083
2084       ENDIF
2085 !-----------------------------------------
2086 ! N8 : copy a helix from bank to seed
2087 !
2088       if (n8.gt.0) then
2089        if (n8frag.lt.n8) then
2090           write (iout,*) "N8: only ",n8frag,'helices'
2091           n8c=n8frag
2092        else
2093           n8c=n8
2094        endif
2095
2096        do iters=1,nseed
2097          iseed=is(iters)
2098          do i=1,mxio
2099            ifused(i)=.false.
2100          enddo
2101
2102
2103          do idummy=1,n8c
2104            iter=0
2105    94      continue
2106            iran=0
2107            iif=iran_num(1,n8frag)
2108            do while( (ifused(iif) .or. hvar_frag(iif,1).eq.iseed) .and. &
2109                             iran.le.mxio )
2110             iif=iran_num(1,n8frag)
2111             iran=iran+1
2112            enddo
2113
2114            if(iran.ge.mxio) goto 911
2115
2116            index=index+1
2117            movenx(index)=8
2118            parent(1,index)=iseed
2119            parent(2,index)=hvar_frag(iif,1)
2120
2121
2122            if (vdisulf) then
2123              nss_in(index)=bvar_nss(iseed)
2124              do ij=1,nss_in(index)
2125                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2126                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2127              enddo
2128            endif
2129
2130            ifused(iif)=.true.
2131            if (hvar_frag(iif,3)-hvar_frag(iif,2).le.6) then
2132              call newconf_copy(idum,dihang_in(1,1,1,index),&
2133                    hvar_frag(iif,1),hvar_frag(iif,2),hvar_frag(iif,3))
2134            else
2135               ih_start=iran_num(hvar_frag(iif,2),hvar_frag(iif,3)-6)
2136               ih_end=iran_num(ih_start,hvar_frag(iif,3))
2137              call newconf_copy(idum,dihang_in(1,1,1,index),&
2138                    hvar_frag(iif,1),ih_start,ih_end)
2139            endif
2140            iter=iter+1
2141            if(iter.lt.10) then
2142             call check_old(icheck,index)
2143             if(icheck.eq.1) then
2144               index=index-1
2145               ifused(iif)=.false.
2146               goto 94
2147             endif
2148            endif
2149
2150
2151   911     continue
2152
2153          enddo
2154        enddo
2155
2156       endif
2157
2158 !-----------------------------------------
2159 ! N7 : copy nonlocal beta fragment from bank to seed
2160 !
2161       if (n7.gt.0) then
2162        if (n7frag.lt.n7) then
2163           write (iout,*) "N7: only ",n7frag,'nonlocal fragments'
2164           n7c=n7frag
2165        else
2166           n7c=n7
2167        endif
2168
2169        do i=1,nres
2170          do j=1,mxio2
2171             iff_in(i,j)=0
2172          enddo
2173        enddo
2174        index2=0
2175        do i=1,mxio
2176          isend2(i)=0
2177        enddo
2178
2179        do iters=1,nseed
2180          iseed=is(iters)
2181          do i=1,mxio
2182            ifused(i)=.false.
2183          enddo
2184
2185          do idummy=1,n7c
2186            iran=0
2187            iif=iran_num(1,n7frag)
2188            do while( (ifused(iif) .or. bvar_frag(iif,1).eq.iseed) .and. &
2189                             iran.le.mxio )
2190             iif=iran_num(1,n7frag)
2191             iran=iran+1
2192            enddo
2193
2194 !d       write (*,'(3i5,l,4i5)'),iters,idummy,iif,ifused(iif),
2195 !d     &                    bvar_frag(iif,1),iseed,iran,index2
2196
2197            if(iran.ge.mxio) goto 999
2198            if(index2.ge.mxio2) goto 999
2199
2200            index=index+1
2201            movenx(index)=7
2202            parent(1,index)=iseed
2203            parent(2,index)=bvar_frag(iif,1)
2204            index2=index2+1
2205            isend2(index)=index2
2206            ifused(iif)=.true.
2207
2208            if (vdisulf) then
2209              nss_in(index)=bvar_nss(iseed)
2210              do ij=1,nss_in(index)
2211                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2212                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2213              enddo
2214            endif
2215
2216            do k=1,numch
2217             do j=2,nres-1
2218              do i=1,4
2219               dihang_in2(i,j,k,index2)=bvar(i,j,k,bvar_frag(iif,1))
2220              enddo
2221             enddo
2222            enddo
2223
2224            if (bvar_frag(iif,2).eq.4) then
2225              do i=bvar_frag(iif,3),bvar_frag(iif,4)
2226                iff_in(i,index2)=1
2227              enddo
2228              if (bvar_frag(iif,5).lt.bvar_frag(iif,6)) then
2229 !d               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
2230 !d     &                 bvar_frag(iif,5),bvar_frag(iif,6)
2231                do i=bvar_frag(iif,5),bvar_frag(iif,6)
2232                  iff_in(i,index2)=1
2233                enddo
2234              else
2235 !d               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
2236 !d     &                 bvar_frag(iif,6),bvar_frag(iif,5)
2237                do i=bvar_frag(iif,6),bvar_frag(iif,5)
2238                  iff_in(i,index2)=1
2239                enddo
2240              endif
2241            endif
2242
2243            do k=1,numch
2244             do j=2,nres-1
2245              do i=1,4
2246               dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2247              enddo
2248             enddo
2249            enddo
2250
2251
2252   999     continue
2253
2254          enddo
2255        enddo
2256
2257       endif
2258 !-----------------------------------------------
2259 ! N6 : copy random continues fragment from bank to seed
2260 !
2261       do iters=1,nseed
2262        iseed=is(iters)
2263        do idummy=1,n6
2264         isize=(is2-is1+1)*ran1(idum)+is1
2265         index=index+1
2266         movenx(index)=6
2267
2268
2269            if (vdisulf) then
2270              nss_in(index)=bvar_nss(iseed)
2271              do ij=1,nss_in(index)
2272                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2273                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2274              enddo
2275            endif
2276
2277         iter=0
2278   104   continue
2279         if(icycle.le.0) then
2280          i1=nconf*  ran1(idum)+1
2281          i1=nbank-nconf+i1
2282         else
2283          i1=nbank*  ran1(idum)+1
2284         endif
2285         if(i1.eq.iseed) goto 104
2286         iter=iter+1
2287         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2288         parent(1,index)=iseed
2289         parent(2,index)=i1
2290         if(iter.lt.10) then
2291          call check_old(icheck,index)
2292          if(icheck.eq.1) goto 104
2293         endif
2294        enddo
2295       enddo
2296 !-----------------------------------------
2297       if (n3.gt.0.or.n4.gt.0) call gen_hairpin
2298       nconf_harp=0
2299       do iters=1,nseed
2300         if (nharp_seed(iters).gt.0) nconf_harp=nconf_harp+1
2301       enddo
2302 !-----------------------------------------
2303 ! N3 : copy hairpin from bank to seed
2304 !
2305       do iters=1,nseed
2306        iseed=is(iters)
2307        nsucc=0
2308        nacc=0
2309        do idummy=1,n3
2310         index=index+1
2311         iter=0
2312   124   continue
2313         if(icycle.le.0) then
2314          i1=nconf*  ran1(idum)+1
2315          i1=nbank-nconf+i1
2316         else
2317          i1=nbank*  ran1(idum)+1
2318         endif
2319         if(i1.eq.iseed) goto 124
2320         do k=1,nsucc
2321           if (i1.eq.iisucc(k).and.nsucc.lt.nconf_harp-1) goto 124
2322         enddo
2323         nsucc=nsucc+1
2324         iisucc(nsucc)=i1
2325         iter=iter+1
2326         call newconf_residue_hairpin(idum,dihang_in(1,1,1,index),&
2327            i1,fail)
2328         if (fail) then
2329           if (icycle.le.0 .and. nsucc.eq.nconf .or. &
2330               icycle.gt.0 .and. nsucc.eq.nbank)  then
2331             index=index-1
2332             goto 125
2333           else
2334             goto 124
2335           endif
2336         endif
2337         if(iter.lt.10) then
2338          call check_old(icheck,index)
2339          if(icheck.eq.1) goto 124
2340         endif
2341         movenx(index)=3
2342         parent(1,index)=iseed
2343         parent(2,index)=i1
2344
2345
2346            if (vdisulf) then
2347              nss_in(index)=bvar_nss(iseed)
2348              do ij=1,nss_in(index)
2349                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2350                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2351              enddo
2352            endif
2353
2354         nacc=nacc+1
2355        enddo
2356 ! if not enough hairpins, supplement with windows
2357   125  continue
2358 !dd       if (n3.ne.0) write (iout,*) "N3",n3," nsucc",nsucc," nacc",nacc
2359        do idummy=nacc+1,n3
2360         isize=(is2-is1+1)*ran1(idum)+is1
2361         index=index+1
2362         movenx(index)=6
2363         parent(1,index)=iseed
2364         parent(2,index)=i1
2365
2366
2367            if (vdisulf) then
2368              nss_in(index)=bvar_nss(iseed)
2369              do ij=1,nss_in(index)
2370                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2371                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2372              enddo
2373            endif
2374
2375         iter=0
2376   114   continue
2377         if(icycle.le.0) then
2378          i1=nconf*  ran1(idum)+1
2379          i1=nbank-nconf+i1
2380         else
2381          i1=nbank*  ran1(idum)+1
2382         endif
2383         if(i1.eq.iseed) goto 114
2384         iter=iter+1
2385         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2386         if(iter.lt.10) then
2387          call check_old(icheck,index)
2388          if(icheck.eq.1) goto 114
2389         endif
2390        enddo
2391       enddo
2392 !-----------------------------------------
2393 ! N4 : shift a turn in hairpin in seed
2394 !
2395       IF (N4.GT.0) THEN
2396       if (4*nharp_tot .ge. n4*nseed) then
2397         ntot_gen=n4*nseed
2398       else
2399         ntot_gen=(4*nharp_tot/nseed)*nseed
2400       endif
2401       ngen=0
2402       do while (ngen.lt.ntot_gen)
2403       do iters=1,nseed
2404         iseed=is(iters)
2405 !        write (iout,*) 'iters',iters,' iseed',iseed,' nharp_seed',
2406 !     &     nharp_seed(iters),' nharp_use',nharp_use(iters),
2407 !     &     ' ntot_gen',ntot_gen
2408 !        write (iout,*) 'iharp_use(0)',
2409 !     &     (iharp_use(0,k,iters),k=1,nharp_seed(iters))
2410         if (nharp_use(iters).gt.0) then
2411           nicht_getan=.true.
2412           do while (nicht_getan)
2413             iih=iran_num(1,nharp_seed(iters))
2414 !            write (iout,*) 'iih',iih,' iharp_use',
2415 !     &            (iharp_use(k,iih,iters),k=1,4)
2416             if (iharp_use(0,iih,iters).gt.0) then
2417               nicht_getan1=.true.
2418               do while (nicht_getan1)
2419                 iim=iran_num(1,4)
2420                 nicht_getan1=iharp_use(iim,iih,iters).eq.1
2421               enddo
2422               nicht_getan=.false.
2423               iharp_use(iim,iih,iters)=1
2424               iharp_use(0,iih,iters)=iharp_use(0,iih,iters)-1
2425               nharp_use(iters)=nharp_use(iters)-1
2426 !dd              write (iout,'(a16,i3,a5,i2,a10,2i4)')
2427 !dd     &           'N4 selected hairpin',iih,' move',iim,' iharp_seed',
2428 !dd     &            iharp_seed(1,iih,iters),iharp_seed(2,iih,iters)
2429             endif
2430           enddo
2431           ngen=ngen+1
2432           index=index+1
2433           movenx(index)=4
2434           parent(1,index)=iseed
2435           parent(2,index)=0
2436
2437
2438            if (vdisulf) then
2439              nss_in(index)=bvar_nss(iseed)
2440              do ij=1,nss_in(index)
2441                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2442                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2443              enddo
2444            endif
2445
2446           do k=1,numch
2447             do j=2,nres-1
2448               do i=1,4
2449                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2450               enddo
2451             enddo
2452           enddo
2453           jstart=iharp_seed(1,iih,iters)+1
2454           jend=iharp_seed(2,iih,iters)
2455           if (iim.eq.1) then
2456             ishift=-2
2457           else if (iim.eq.2) then
2458             ishift=-1
2459           else if (iim.eq.3) then
2460             ishift=1
2461           else if (iim.eq.4) then
2462             ishift=2
2463           else
2464             write (iout,*) 'CHUJ NASTAPIL: iim=',iim
2465 #ifdef MPI !el
2466             call mpi_abort(mpi_comm_world,ierror,ierrcode)
2467 #endif !el
2468           endif
2469 !          write (iout,*) 'jstart',jstart,' jend',jend,' ishift',ishift
2470 !          write (iout,*) 'Before turn shift'
2471 !          do j=2,nres-1
2472 !            theta(j+1)=dihang_in(1,j,1,index)
2473 !            phi(j+2)=dihang_in(2,j,1,index)
2474 !            alph(j)=dihang_in(3,j,1,index)
2475 !            omeg(j)=dihang_in(4,j,1,index)
2476 !          enddo
2477 !          call intout
2478           do j=jstart,jend
2479             if (itype(j).eq.10) then
2480               iang=2
2481             else
2482               iang=4
2483             endif
2484             do i=1,iang
2485               if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
2486                 dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
2487             enddo
2488           enddo
2489 !          write (iout,*) 'After turn shift'
2490 !          do j=2,nres-1
2491 !            theta(j+1)=dihang_in(1,j,1,index)
2492 !            phi(j+2)=dihang_in(2,j,1,index)
2493 !            alph(j)=dihang_in(3,j,1,index)
2494 !            omeg(j)=dihang_in(4,j,1,index)
2495 !          enddo
2496 !          call intout
2497           if (ngen.eq.ntot_gen) goto 135
2498         endif
2499       enddo
2500       enddo
2501 ! if not enough hairpins, supplement with windows
2502 !      write (iout,*) 'end of enddo'
2503   135 continue
2504 !dd      write (iout,*) "N4",n4," ngen/nseed",ngen/nseed,
2505 !dd     &                           ngen,nseed
2506       do iters=1,nseed
2507         iseed=is(iters)
2508         do idummy=ngen/nseed+1,n4
2509           isize=(is2-is1+1)*ran1(idum)+is1
2510           index=index+1
2511           movenx(index)=6
2512
2513            if (vdisulf) then
2514              nss_in(index)=bvar_nss(iseed)
2515              do ij=1,nss_in(index)
2516                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2517                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2518              enddo
2519            endif
2520
2521
2522           iter=0
2523   134     continue
2524           if(icycle.le.0) then
2525           i1=nconf*  ran1(idum)+1
2526           i1=nbank-nconf+i1
2527           else
2528             i1=nbank*  ran1(idum)+1
2529           endif
2530           if(i1.eq.iseed) goto 134
2531             iter=iter+1
2532             call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2533             parent(1,index)=iseed
2534             parent(2,index)=i1
2535             if(iter.lt.10) then
2536               call check_old(icheck,index)
2537             if(icheck.eq.1) goto 134
2538           endif
2539         enddo
2540       enddo
2541       ENDIF
2542 !-----------------------------------------
2543 ! N5 : copy one residue from bank to seed (normally switched off - use N1)
2544 !
2545       do iters=1,nseed
2546        iseed=is(iters)
2547        isize=1
2548        do i=1,n5
2549         index=index+1
2550         movenx(index)=5
2551
2552            if (vdisulf) then
2553              nss_in(index)=bvar_nss(iseed)
2554              do ij=1,nss_in(index)
2555                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2556                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2557              enddo
2558            endif
2559
2560
2561         iter=0
2562   105   continue
2563         if(icycle.le.0) then
2564          i1=nconf*  ran1(idum)+1
2565          i1=nbank-nconf+i1
2566         else
2567          i1=nbank*  ran1(idum)+1
2568         endif
2569         if(i1.eq.iseed) goto 105
2570         iter=iter+1
2571         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2572         parent(1,index)=iseed
2573         parent(2,index)=i1
2574         if(iter.lt.10) then
2575          call check_old(icheck,index)
2576          if(icheck.eq.1) goto 105
2577         endif
2578        enddo
2579       enddo
2580 !-----------------------------------------
2581 ! N2 : copy backbone of one residue from bank or first bank to seed
2582 ! (normally switched off - use N1)
2583 !
2584       do iters=1,nseed
2585        iseed=is(iters)
2586        do i=n2,1,-1
2587         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2588          iseed=ran1(idum)*nconf+1
2589          iseed=nbank-nconf+iseed
2590         endif
2591         index=index+1
2592         movenx(index)=2
2593
2594         if (vdisulf) then
2595              nss_in(index)=bvar_nss(iseed)
2596              do ij=1,nss_in(index)
2597                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2598                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2599              enddo
2600         endif
2601
2602         iter=0
2603   102   i1=  ran1(idum)*nbank+1
2604         if(i1.eq.iseed) goto 102
2605         iter=iter+1
2606         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2607          nran=mod(i-1,nran0)+3
2608          call newconf1arr(idum,dihang_in(1,1,1,index),nran,i1)
2609          parent(1,index)=-iseed
2610          parent(2,index)=-i1
2611         else if(icycle.le.0.and.iters.le.iuse) then
2612          nran=mod(i-1,nran0)+1
2613          call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
2614          parent(1,index)=iseed
2615          parent(2,index)=-i1
2616         else
2617          nran=mod(i-1,nran1)+1
2618          if(ran1(idum).lt.0.5) then
2619           call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
2620           parent(1,index)=iseed
2621           parent(2,index)=-i1
2622          else
2623           call newconf1abb(idum,dihang_in(1,1,1,index),nran,i1)
2624           parent(1,index)=iseed
2625           parent(2,index)=i1
2626          endif
2627         endif
2628         if(iter.lt.10) then
2629          call check_old(icheck,index)
2630          if(icheck.eq.1) goto 102
2631         endif
2632        enddo
2633       enddo
2634 !-----------------------------------------
2635 ! N1 : copy backbone or sidechain of one residue from bank or
2636 !      first bank to seed
2637 !
2638       do iters=1,nseed
2639        iseed=is(iters)
2640        do i=n1,1,-1
2641         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2642          iseed=ran1(idum)*nconf+1
2643          iseed=nbank-nconf+iseed
2644         endif
2645         index=index+1
2646         movenx(index)=1
2647
2648         if (vdisulf) then
2649              nss_in(index)=bvar_nss(iseed)
2650              do ij=1,nss_in(index)
2651                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2652                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2653              enddo
2654         endif
2655
2656         iter=0
2657   101   i1=  ran1(idum)*nbank+1
2658
2659         if(i1.eq.iseed) goto 101
2660         iter=iter+1
2661         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2662          nran=mod(i-1,nran0)+3
2663          call newconf1rr(idum,dihang_in(1,1,1,index),nran,i1)
2664          parent(1,index)=-iseed
2665          parent(2,index)=-i1
2666         else if(icycle.le.0.and.iters.le.iuse) then
2667          nran=mod(i-1,nran0)+1
2668          call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
2669          parent(1,index)=iseed
2670          parent(2,index)=-i1
2671         else
2672          nran=mod(i-1,nran1)+1
2673          if(ran1(idum).lt.0.5) then
2674           call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
2675           parent(1,index)=iseed
2676           parent(2,index)=-i1
2677          else
2678           call newconf1bb(idum,dihang_in(1,1,1,index),nran,i1)
2679           parent(1,index)=iseed
2680           parent(2,index)=i1
2681          endif
2682         endif
2683         if(iter.lt.10) then
2684          call check_old(icheck,index)
2685          if(icheck.eq.1) goto 101
2686         endif
2687        enddo
2688       enddo
2689 !-----------------------------------------
2690 ! N0 just all seeds
2691 !
2692       IF (n0.gt.0) THEN
2693       do iters=1,nseed
2694        iseed=is(iters)
2695        index=index+1
2696        movenx(index)=0
2697        parent(1,index)=iseed
2698        parent(2,index)=0
2699
2700        if (vdisulf) then
2701              nss_in(index)=bvar_nss(iseed)
2702              do ij=1,nss_in(index)
2703                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2704                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2705              enddo
2706        endif
2707
2708        do k=1,numch
2709         do j=2,nres-1
2710          do i=1,4
2711           dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2712          enddo
2713         enddo
2714        enddo
2715       enddo
2716       ENDIF
2717 !-----------------------------------------
2718       if (vdisulf) then
2719       do iters=1,nseed
2720        iseed=is(iters)
2721
2722         do k=1,numch
2723           do j=2,nres-1
2724             theta(j+1)=bvar(1,j,k,iseed)
2725             phi(j+2)=bvar(2,j,k,iseed)
2726             alph(j)=bvar(3,j,k,iseed)
2727             omeg(j)=bvar(4,j,k,iseed)
2728           enddo
2729         enddo
2730         call chainbuild
2731
2732 !d       write(iout,*) 'makevar DYNSS',iseed,'#',bvar_ns(iseed),
2733 !d     &     (bvar_s(k,iseed),k=1,bvar_ns(iseed)),
2734 !d     &     bvar_nss(iseed),
2735 !d     &     (bvar_ss(1,k,iseed)-nres,'-',
2736 !d     &      bvar_ss(2,k,iseed)-nres,k=1,bvar_nss(iseed))
2737
2738        do i1=1,bvar_ns(iseed)
2739 !
2740 ! N10 fussion of free halfcysteines in seed
2741 !      first select CYS with distance < 7A
2742 !
2743          do j1=i1+1,bvar_ns(iseed)
2744            if (dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres) &
2745                .lt.7.0.and. &
2746                iabs(bvar_s(i1,iseed)-bvar_s(j1,iseed)).gt.3) then
2747
2748              index=index+1
2749              movenx(index)=10
2750              parent(1,index)=iseed
2751              parent(2,index)=0
2752              do ij=1,bvar_nss(iseed)
2753                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2754                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2755              enddo
2756              ij=bvar_nss(iseed)+1
2757              nss_in(index)=ij
2758              iss_in(ij,index)=bvar_s(i1,iseed)+nres
2759              jss_in(ij,index)=bvar_s(j1,iseed)+nres
2760
2761 !d             write(iout,*) 'makevar NSS0',index,
2762 !d     &   dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres),
2763 !d     &   nss_in(index),iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres
2764
2765              do k=1,numch
2766               do j=2,nres-1
2767                do i=1,4
2768                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2769                enddo
2770               enddo
2771              enddo
2772
2773            endif
2774          enddo
2775 !
2776 ! N11 type I transdisulfidation
2777 !
2778          do j1=1,bvar_nss(iseed)
2779            if (dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)) &
2780                .lt.7.0.and. &
2781                iabs(bvar_s(i1,iseed)-(bvar_ss(1,j1,iseed)-nres)) &
2782                .gt.3) then
2783
2784              index=index+1
2785              movenx(index)=11
2786              parent(1,index)=iseed
2787              parent(2,index)=0
2788              do ij=1,bvar_nss(iseed)
2789                if (ij.ne.j1) then
2790                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2791                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2792                endif
2793              enddo
2794              nss_in(index)=bvar_nss(iseed)
2795              iss_in(j1,index)=bvar_s(i1,iseed)+nres
2796              jss_in(j1,index)=bvar_ss(1,j1,iseed)
2797              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2798                iss_in(j1,index)=bvar_ss(1,j1,iseed)
2799                jss_in(j1,index)=bvar_s(i1,iseed)+nres
2800              endif
2801
2802 !d             write(iout,*) 'makevar NSS1 #1',index,
2803 !d     &       bvar_s(i1,iseed),bvar_ss(1,j1,iseed)-nres,
2804 !d     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)),
2805 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2806 !d     &       ij=1,nss_in(index))
2807
2808              do k=1,numch
2809               do j=2,nres-1
2810                do i=1,4
2811                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2812                enddo
2813               enddo
2814              enddo
2815            endif
2816            if (dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)) &
2817                .lt.7.0.and. &
2818                iabs(bvar_s(i1,iseed)-(bvar_ss(2,j1,iseed)-nres)) &
2819                .gt.3) then
2820
2821              index=index+1
2822              movenx(index)=11
2823              parent(1,index)=iseed
2824              parent(2,index)=0
2825              do ij=1,bvar_nss(iseed)
2826                if (ij.ne.j1) then
2827                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2828                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2829                endif
2830              enddo
2831              nss_in(index)=bvar_nss(iseed)
2832              iss_in(j1,index)=bvar_s(i1,iseed)+nres
2833              jss_in(j1,index)=bvar_ss(2,j1,iseed)
2834              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2835                iss_in(j1,index)=bvar_ss(2,j1,iseed)
2836                jss_in(j1,index)=bvar_s(i1,iseed)+nres
2837              endif
2838
2839
2840 !d             write(iout,*) 'makevar NSS1 #2',index,
2841 !d     &       bvar_s(i1,iseed),bvar_ss(2,j1,iseed)-nres,
2842 !d     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)),
2843 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2844 !d     &       ij=1,nss_in(index))
2845
2846              do k=1,numch
2847               do j=2,nres-1
2848                do i=1,4
2849                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2850                enddo
2851               enddo
2852              enddo
2853
2854            endif
2855          enddo
2856        enddo
2857
2858 !
2859 ! N12 type II transdisulfidation
2860 !
2861        do i1=1,bvar_nss(iseed)
2862          do j1=i1+1,bvar_nss(iseed)
2863             if (dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)) &
2864                 .lt.7.0.and. &
2865                 dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)) &
2866                 .lt.7.0.and. &
2867                 iabs(bvar_ss(1,i1,iseed)-bvar_ss(1,j1,iseed)) &
2868                 .gt.3.and. &
2869                 iabs(bvar_ss(2,i1,iseed)-bvar_ss(2,j1,iseed)) &
2870                 .gt.3) then
2871              index=index+1
2872              movenx(index)=12
2873              parent(1,index)=iseed
2874              parent(2,index)=0
2875              do ij=1,bvar_nss(iseed)
2876                if (ij.ne.i1 .and. ij.ne.j1) then
2877                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2878                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2879                endif
2880              enddo
2881              nss_in(index)=bvar_nss(iseed)
2882              iss_in(i1,index)=bvar_ss(1,i1,iseed)
2883              jss_in(i1,index)=bvar_ss(1,j1,iseed)
2884              if (iss_in(i1,index).gt.jss_in(i1,index)) then
2885                iss_in(i1,index)=bvar_ss(1,j1,iseed)
2886                jss_in(i1,index)=bvar_ss(1,i1,iseed)
2887              endif
2888              iss_in(j1,index)=bvar_ss(2,i1,iseed)
2889              jss_in(j1,index)=bvar_ss(2,j1,iseed)
2890              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2891                iss_in(j1,index)=bvar_ss(2,j1,iseed)
2892                jss_in(j1,index)=bvar_ss(2,i1,iseed)
2893              endif
2894
2895
2896 !d             write(iout,*) 'makevar NSS2 #1',index,
2897 !d     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
2898 !d     &       dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)),
2899 !d     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
2900 !d     &       dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)),
2901 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2902 !d     &       ij=1,nss_in(index))
2903
2904              do k=1,numch
2905               do j=2,nres-1
2906                do i=1,4
2907                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2908                enddo
2909               enddo
2910              enddo
2911
2912            endif
2913
2914             if (dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)) &
2915                 .lt.7.0.and. &
2916                 dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)) &
2917                 .lt.7.0.and. &
2918                 iabs(bvar_ss(1,i1,iseed)-bvar_ss(2,j1,iseed)) &
2919                 .gt.3.and. &
2920                 iabs(bvar_ss(2,i1,iseed)-bvar_ss(1,j1,iseed)) &
2921                 .gt.3) then
2922              index=index+1
2923              movenx(index)=12
2924              parent(1,index)=iseed
2925              parent(2,index)=0
2926              do ij=1,bvar_nss(iseed)
2927                if (ij.ne.i1 .and. ij.ne.j1) then
2928                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2929                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2930                endif
2931              enddo
2932              nss_in(index)=bvar_nss(iseed)
2933              iss_in(i1,index)=bvar_ss(1,i1,iseed)
2934              jss_in(i1,index)=bvar_ss(2,j1,iseed)
2935              if (iss_in(i1,index).gt.jss_in(i1,index)) then
2936                iss_in(i1,index)=bvar_ss(2,j1,iseed)
2937                jss_in(i1,index)=bvar_ss(1,i1,iseed)
2938              endif
2939              iss_in(j1,index)=bvar_ss(2,i1,iseed)
2940              jss_in(j1,index)=bvar_ss(1,j1,iseed)
2941              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2942                iss_in(j1,index)=bvar_ss(1,j1,iseed)
2943                jss_in(j1,index)=bvar_ss(2,i1,iseed)
2944              endif
2945
2946
2947 !d             write(iout,*) 'makevar NSS2 #2',index,
2948 !d     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
2949 !d     &       dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)),
2950 !d     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
2951 !d     &       dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)),
2952 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2953 !d     &       ij=1,nss_in(index))
2954
2955              do k=1,numch
2956               do j=2,nres-1
2957                do i=1,4
2958                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2959                enddo
2960               enddo
2961              enddo
2962
2963            endif
2964
2965
2966          enddo
2967        enddo
2968 !
2969 ! N13 removal of disulfide bond
2970 !
2971       if (bvar_nss(iseed).gt.0) then
2972         i1=bvar_nss(iseed)*ran1(idum)+1
2973
2974              index=index+1
2975              movenx(index)=13
2976              parent(1,index)=iseed
2977              parent(2,index)=0
2978              ij=0
2979              do j1=1,bvar_nss(iseed)
2980               if (j1.ne.i1) then
2981                ij=ij+1
2982                iss_in(ij,index)=bvar_ss(1,j1,iseed)
2983                jss_in(ij,index)=bvar_ss(2,j1,iseed)
2984               endif
2985              enddo
2986              nss_in(index)=bvar_nss(iseed)-1
2987
2988 !d             write(iout,*) 'NSS3',index,i1,
2989 !d     &   bvar_ss(1,i1,iseed)-nres,'=',bvar_ss(2,i1,iseed)-nres,'#',
2990 !d     &   (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2991 !d     &   ij=1,nss_in(index))
2992
2993              do k=1,numch
2994               do j=2,nres-1
2995                do i=1,4
2996                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2997                enddo
2998               enddo
2999              enddo
3000
3001       endif
3002
3003       enddo
3004       endif
3005 !-----------------------------------------
3006
3007
3008
3009       if(index.ne.n) write(iout,*)'make_var : ntry=',index
3010
3011       n=index
3012 !d      do ii=1,n
3013 !d        write (istat,*) "======== ii=",ii," the dihang array"
3014 !d        do i=1,nres
3015 !d       write (istat,'(i5,4f15.5)') i,(dihang_in(k,i,1,ii)*rad2deg,k=1,4)
3016 !d        enddo
3017 !d      enddo
3018       return
3019       end subroutine make_var
3020 !-----------------------------------------------------------------------------
3021       subroutine check_old(icheck,n)
3022
3023 !      implicit real*8 (a-h,o-z)
3024 !      include 'DIMENSIONS'
3025 !      include 'COMMON.CSA'
3026 !      include 'COMMON.BANK'
3027 !      include 'COMMON.CHAIN'
3028 !      include 'COMMON.GEO'
3029       integer :: icheck,n,i1,i2,m,j,i
3030       real(kind=8) :: ctdif,ctdiff,diff,dif
3031
3032       data ctdif /10./
3033       data ctdiff /60./
3034
3035       i1=n
3036       do i2=1,n-1
3037        diff=0.d0
3038        do m=1,numch
3039         do j=2,nres-1
3040          do i=1,4
3041           dif=rad2deg*dabs(dihang_in(i,j,m,i1)-dihang_in(i,j,m,i2))
3042           if(dif.gt.180.0) dif=360.0-dif
3043           if(dif.gt.ctdif) goto 100
3044           diff=diff+dif
3045           if(diff.gt.ctdiff) goto 100
3046          enddo
3047         enddo
3048        enddo
3049        icheck=1
3050        return
3051   100  continue
3052       enddo
3053
3054       icheck=0
3055
3056       return
3057       end subroutine check_old
3058 !-----------------------------------------------------------------------------
3059       subroutine newconf1rr(idum,vvar,nran,i1)
3060
3061 !      implicit real*8 (a-h,o-z)
3062 !      include 'DIMENSIONS'
3063 !      include 'COMMON.IOUNITS'
3064 !      include 'COMMON.CSA'
3065 !      include 'COMMON.BANK'
3066 !      include 'COMMON.CHAIN'
3067 !      include 'COMMON.GEO'
3068 !      real(kind=4) :: ran1,ran2
3069       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3070       integer,dimension(ntotal) :: iold
3071       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3072       real(kind=8) :: ctdif,dif
3073
3074       ctdif=10.
3075
3076       do k=1,numch
3077        do j=2,nres-1
3078         do i=1,4
3079          vvar(i,j,k)=rvar(i,j,k,iseed)
3080         enddo
3081        enddo
3082       enddo
3083
3084       do index=1,nran
3085        iold(index) = 0
3086       enddo
3087
3088       number=ntotgr
3089
3090       iter=0
3091       do index=1,nran
3092    10  iran=  ran1(idum)*number+1
3093        if(iter.gt.number) return
3094        iter=iter+1
3095        if(iter.eq.1) goto 11
3096        do ind=1,index-1
3097         if(iran.eq.iold(ind)) goto 10
3098        enddo
3099    11  continue
3100
3101        do ind=1,ngroup(iran)
3102         i=igroup(1,ind,iran)
3103         j=igroup(2,ind,iran)
3104         k=igroup(3,ind,iran)
3105         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3106         if(dif.gt.180.) dif=360.-dif
3107         if(dif.gt.ctdif) goto 20
3108        enddo
3109        if(iter.gt.number) goto 20
3110        goto 10
3111    20  continue
3112        do ind=1,ngroup(iran)
3113         i=igroup(1,ind,iran)
3114         j=igroup(2,ind,iran)
3115         k=igroup(3,ind,iran)
3116         vvar(i,j,k)=rvar(i,j,k,i1)
3117        enddo
3118        iold(index)=iran
3119       enddo
3120
3121       return
3122       end subroutine newconf1rr
3123 !-----------------------------------------------------------------------------
3124       subroutine newconf1br(idum,vvar,nran,i1)
3125
3126       use energy_data, only: ndih_nconstr,idih_nconstr
3127       use control_data, only: i2ndstr
3128 !      implicit real*8 (a-h,o-z)
3129 !      include 'DIMENSIONS'
3130 !      include 'COMMON.IOUNITS'
3131 !      include 'COMMON.CSA'
3132 !      include 'COMMON.BANK'
3133 !      include 'COMMON.CHAIN'
3134 !      include 'COMMON.GEO'
3135 !      include 'COMMON.TORCNSTR'
3136 !      include 'COMMON.CONTROL'
3137 !      real(kind=4) :: ran1,ran2
3138       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3139       integer,dimension(ntotal) :: iold
3140       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,juhc,ind
3141       real(kind=8) :: ctdif,dif,rtmp
3142
3143       ctdif=10.
3144
3145       do k=1,numch
3146        do j=2,nres-1
3147         do i=1,4
3148          vvar(i,j,k)=bvar(i,j,k,iseed)
3149         enddo
3150        enddo
3151       enddo
3152
3153       do index=1,nran
3154        iold(index) = 0
3155       enddo
3156
3157       number=ntotgr
3158
3159       iter=0
3160       do index=1,nran
3161    10  iran=  ran1(idum)*number+1
3162        if(i2ndstr.gt.0) then
3163         rtmp=ran1(idum)
3164         if(rtmp.le.rdih_bias) then
3165          i=0
3166          do j=1,ndih_nconstr
3167           if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
3168          enddo
3169          if(i.eq.0) then
3170           juhc=0
3171 4321      juhc=juhc+1
3172           iran=  ran1(idum)*number+1
3173           i=0
3174           do j=1,ndih_nconstr
3175            if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
3176           enddo
3177           if(i.eq.0.or.juhc.lt.1000)goto 4321
3178           if(juhc.eq.1000) then
3179            print *, 'move 6 : failed to find unconstrained group'
3180            write(iout,*) 'move 6 : failed to find unconstrained group'
3181           endif
3182          endif
3183         endif
3184        endif
3185        if(iter.gt.number) return
3186        iter=iter+1
3187        if(iter.eq.1) goto 11
3188        do ind=1,index-1
3189         if(iran.eq.iold(ind)) goto 10
3190        enddo
3191    11  continue
3192
3193        do ind=1,ngroup(iran)
3194         i=igroup(1,ind,iran)
3195         j=igroup(2,ind,iran)
3196         k=igroup(3,ind,iran)
3197         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3198         if(dif.gt.180.) dif=360.-dif
3199         if(dif.gt.ctdif) goto 20
3200        enddo
3201        if(iter.gt.number) goto 20
3202        goto 10
3203    20  continue
3204        do ind=1,ngroup(iran)
3205         i=igroup(1,ind,iran)
3206         j=igroup(2,ind,iran)
3207         k=igroup(3,ind,iran)
3208         vvar(i,j,k)=rvar(i,j,k,i1)
3209        enddo
3210        iold(index)=iran
3211       enddo
3212
3213       return
3214       end subroutine newconf1br
3215 !-----------------------------------------------------------------------------
3216       subroutine newconf1bb(idum,vvar,nran,i1)
3217
3218 !      implicit real*8 (a-h,o-z)
3219 !      include 'DIMENSIONS'
3220 !      include 'COMMON.IOUNITS'
3221 !      include 'COMMON.CSA'
3222 !      include 'COMMON.BANK'
3223 !      include 'COMMON.CHAIN'
3224 !      include 'COMMON.GEO'
3225 !      real(kind=4) :: ran1,ran2
3226       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3227       integer,dimension(ntotal) :: iold
3228       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3229       real(kind=8) :: ctdif,dif
3230
3231       ctdif=10.
3232
3233       do k=1,numch
3234        do j=2,nres-1
3235         do i=1,4
3236          vvar(i,j,k)=bvar(i,j,k,iseed)
3237         enddo
3238        enddo
3239       enddo
3240
3241       do index=1,nran
3242        iold(index) = 0
3243       enddo
3244
3245       number=ntotgr
3246
3247       iter=0
3248       do index=1,nran
3249    10  iran=  ran1(idum)*number+1
3250        if(iter.gt.number) return
3251        iter=iter+1
3252        if(iter.eq.1) goto 11
3253        do ind=1,index-1
3254         if(iran.eq.iold(ind)) goto 10
3255        enddo
3256    11  continue
3257
3258        do ind=1,ngroup(iran)
3259         i=igroup(1,ind,iran)
3260         j=igroup(2,ind,iran)
3261         k=igroup(3,ind,iran)
3262         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3263         if(dif.gt.180.) dif=360.-dif
3264         if(dif.gt.ctdif) goto 20
3265        enddo
3266        if(iter.gt.number) goto 20
3267        goto 10
3268    20  continue
3269        do ind=1,ngroup(iran)
3270         i=igroup(1,ind,iran)
3271         j=igroup(2,ind,iran)
3272         k=igroup(3,ind,iran)
3273         vvar(i,j,k)=bvar(i,j,k,i1)
3274        enddo
3275        iold(index)=iran
3276       enddo
3277
3278       return
3279       end subroutine newconf1bb
3280 !-----------------------------------------------------------------------------
3281       subroutine newconf1arr(idum,vvar,nran,i1)
3282
3283 !      implicit real*8 (a-h,o-z)
3284 !      include 'DIMENSIONS'
3285 !      include 'COMMON.IOUNITS'
3286 !      include 'COMMON.CSA'
3287 !      include 'COMMON.BANK'
3288 !      include 'COMMON.CHAIN'
3289 !      include 'COMMON.GEO'
3290 !      real(kind=4) :: ran1,ran2
3291       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3292       integer,dimension(ntotal) :: iold
3293       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3294       real(kind=8) :: ctdif,dif
3295
3296       ctdif=10.
3297
3298       do k=1,numch
3299        do j=2,nres-1
3300         do i=1,4
3301          vvar(i,j,k)=rvar(i,j,k,iseed)
3302         enddo
3303        enddo
3304       enddo
3305
3306       do index=1,nran
3307        iold(index) = 0
3308       enddo
3309
3310       number=nres-2
3311
3312       iter=0
3313       do index=1,nran
3314    10  iran=  ran1(idum)*number+1
3315        if(iter.gt.number) return
3316        iter=iter+1
3317        if(iter.eq.1) goto 11
3318        do ind=1,index-1
3319         if(iran.eq.iold(ind)) goto 10
3320        enddo
3321    11  continue
3322
3323        do ind=1,ngroup(iran)
3324         i=igroup(1,ind,iran)
3325         j=igroup(2,ind,iran)
3326         k=igroup(3,ind,iran)
3327         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3328         if(dif.gt.180.) dif=360.-dif
3329         if(dif.gt.ctdif) goto 20
3330        enddo
3331        if(iter.gt.number) goto 20
3332        goto 10
3333    20  continue
3334        do ind=1,ngroup(iran)
3335         i=igroup(1,ind,iran)
3336         j=igroup(2,ind,iran)
3337         k=igroup(3,ind,iran)
3338         vvar(i,j,k)=rvar(i,j,k,i1)
3339        enddo
3340        iold(index)=iran
3341       enddo
3342
3343       return
3344       end subroutine newconf1arr
3345 !-----------------------------------------------------------------------------
3346       subroutine newconf1abr(idum,vvar,nran,i1)
3347
3348       use energy_data, only: ndih_nconstr,idih_nconstr
3349       use control_data, only: i2ndstr
3350 !      implicit real*8 (a-h,o-z)
3351 !      include 'DIMENSIONS'
3352 !      include 'COMMON.IOUNITS'
3353 !      include 'COMMON.CSA'
3354 !      include 'COMMON.BANK'
3355 !      include 'COMMON.CHAIN'
3356 !      include 'COMMON.GEO'
3357 !      include 'COMMON.TORCNSTR'
3358 !      include 'COMMON.CONTROL'
3359 !      real(kind=4) :: ran1,ran2
3360       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3361       integer,dimension(ntotal) :: iold
3362       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3363       real(kind=8) :: ctdif,dif,rtmp
3364
3365       ctdif=10.
3366
3367       do k=1,numch
3368        do j=2,nres-1
3369         do i=1,4
3370          vvar(i,j,k)=bvar(i,j,k,iseed)
3371         enddo
3372        enddo
3373       enddo
3374
3375       do index=1,nran
3376        iold(index) = 0
3377       enddo
3378
3379       number=nres-2
3380
3381       iter=0
3382       do index=1,nran
3383    10  iran=  ran1(idum)*number+1
3384       if(i2ndstr.gt.0) then
3385        rtmp=ran1(idum)
3386        if(rtmp.le.rdih_bias) then
3387         iran=ran1(idum)*ndih_nconstr+1
3388         iran=idih_nconstr(iran)
3389        endif
3390       endif
3391        if(iter.gt.number) return
3392        iter=iter+1
3393        if(iter.eq.1) goto 11
3394        do ind=1,index-1
3395         if(iran.eq.iold(ind)) goto 10
3396        enddo
3397    11  continue
3398
3399        do ind=1,ngroup(iran)
3400         i=igroup(1,ind,iran)
3401         j=igroup(2,ind,iran)
3402         k=igroup(3,ind,iran)
3403         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3404         if(dif.gt.180.) dif=360.-dif
3405         if(dif.gt.ctdif) goto 20
3406        enddo
3407        if(iter.gt.number) goto 20
3408        goto 10
3409    20  continue
3410        do ind=1,ngroup(iran)
3411         i=igroup(1,ind,iran)
3412         j=igroup(2,ind,iran)
3413         k=igroup(3,ind,iran)
3414         vvar(i,j,k)=rvar(i,j,k,i1)
3415        enddo
3416        iold(index)=iran
3417       enddo
3418
3419       return
3420       end subroutine newconf1abr
3421 !-----------------------------------------------------------------------------
3422       subroutine newconf1abb(idum,vvar,nran,i1)
3423
3424       use energy_data, only: ndih_nconstr,idih_nconstr
3425       use control_data, only: i2ndstr
3426 !      implicit real*8 (a-h,o-z)
3427 !      include 'DIMENSIONS'
3428 !      include 'COMMON.IOUNITS'
3429 !      include 'COMMON.CSA'
3430 !      include 'COMMON.BANK'
3431 !      include 'COMMON.CHAIN'
3432 !      include 'COMMON.GEO'
3433 !      include 'COMMON.TORCNSTR'
3434 !      include 'COMMON.CONTROL'
3435 !      real(kind=4) :: ran1,ran2
3436       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3437       integer,dimension(ntotal) :: iold
3438       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3439       real(kind=8) :: ctdif,dif,rtmp
3440
3441       ctdif=10.
3442
3443       do k=1,numch
3444        do j=2,nres-1
3445         do i=1,4
3446          vvar(i,j,k)=bvar(i,j,k,iseed)
3447         enddo
3448        enddo
3449       enddo
3450
3451       do index=1,nran
3452        iold(index) = 0
3453       enddo
3454
3455       number=nres-2
3456
3457       iter=0
3458       do index=1,nran
3459    10  iran=  ran1(idum)*number+1
3460       if(i2ndstr.gt.0) then
3461        rtmp=ran1(idum)
3462        if(rtmp.le.rdih_bias) then
3463         iran=ran1(idum)*ndih_nconstr+1
3464         iran=idih_nconstr(iran)
3465        endif
3466       endif
3467        if(iter.gt.number) return
3468        iter=iter+1
3469        if(iter.eq.1) goto 11
3470        do ind=1,index-1
3471         if(iran.eq.iold(ind)) goto 10
3472        enddo
3473    11  continue
3474
3475        do ind=1,ngroup(iran)
3476         i=igroup(1,ind,iran)
3477         j=igroup(2,ind,iran)
3478         k=igroup(3,ind,iran)
3479         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3480         if(dif.gt.180.) dif=360.-dif
3481         if(dif.gt.ctdif) goto 20
3482        enddo
3483        if(iter.gt.number) goto 20
3484        goto 10
3485    20  continue
3486        do ind=1,ngroup(iran)
3487         i=igroup(1,ind,iran)
3488         j=igroup(2,ind,iran)
3489         k=igroup(3,ind,iran)
3490         vvar(i,j,k)=bvar(i,j,k,i1)
3491        enddo
3492        iold(index)=iran
3493       enddo
3494
3495       return
3496       end subroutine newconf1abb
3497 !-----------------------------------------------------------------------------
3498       subroutine newconf_residue(idum,vvar,i1,isize)
3499
3500       use energy_data, only: ndih_nconstr,idih_nconstr
3501       use control_data, only: i2ndstr
3502       use MPI_data
3503       include 'mpif.h'
3504 !      implicit real*8 (a-h,o-z)
3505 !      include 'DIMENSIONS'
3506 !      include 'COMMON.IOUNITS'
3507 !      include 'COMMON.CSA'
3508 !      include 'COMMON.BANK'
3509 !      include 'COMMON.CHAIN'
3510 !      include 'COMMON.GEO'
3511 !      include 'COMMON.TORCNSTR'
3512 !      include 'COMMON.CONTROL'
3513 !      real(kind=4) :: ran1,ran2
3514       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3515       integer,dimension(ntotal) :: iold
3516       integer :: i,j,k,idum,i1,isize,iran,number,iter,ind,iend,istart,&
3517             ierror,ierrcode
3518       real(kind=8) :: ctdif,dif,rtmp
3519
3520       ctdif=10.
3521
3522       if (iseed.gt.mxio .or. iseed.lt.1) then
3523         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3524         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3525       endif
3526       do k=1,numch
3527        do j=2,nres-1
3528         do i=1,4
3529          vvar(i,j,k)=bvar(i,j,k,iseed)
3530         enddo
3531        enddo
3532       enddo
3533
3534
3535       k=1
3536       number=nres+isize-2
3537       iter=1
3538    10 iran=  ran1(idum)*number+1
3539       if(i2ndstr.gt.0) then
3540        rtmp=ran1(idum)
3541        if(rtmp.le.rdih_bias) then
3542         iran=ran1(idum)*ndih_nconstr+1
3543         iran=idih_nconstr(iran)
3544        endif
3545       endif
3546       istart=iran-isize+1
3547       iend=iran
3548       if(istart.lt.2) istart=2
3549       if(iend.gt.nres-1) iend=nres-1
3550
3551       if(iter.eq.1) goto 11
3552       do ind=1,iter-1
3553        if(iran.eq.iold(ind)) goto 10
3554       enddo
3555    11 continue
3556
3557       do j=istart,iend
3558        do i=1,4
3559          dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3560          if(dif.gt.180.) dif=360.-dif
3561          if(dif.gt.ctdif) goto 20
3562        enddo
3563       enddo
3564       iold(iter)=iran
3565       iter=iter+1
3566       if(iter.gt.number) goto 20
3567       goto 10
3568
3569    20 continue
3570       do j=istart,iend
3571        do i=1,4
3572         vvar(i,j,k)=bvar(i,j,k,i1)
3573        enddo
3574       enddo
3575
3576       return
3577       end subroutine newconf_residue
3578 !-----------------------------------------------------------------------------
3579       subroutine newconf_copy(idum,vvar,i1,istart,iend)
3580
3581       use MPI_data
3582       include 'mpif.h'
3583 !      implicit real*8 (a-h,o-z)
3584 !      include 'DIMENSIONS'
3585 !      include 'COMMON.IOUNITS'
3586 !      include 'COMMON.CSA'
3587 !      include 'COMMON.BANK'
3588 !      include 'COMMON.CHAIN'
3589 !      include 'COMMON.GEO'
3590 !      include 'COMMON.TORCNSTR'
3591 !      include 'COMMON.CONTROL'
3592 !      real(kind=4) :: ran1,ran2
3593       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3594       integer,dimension(ntotal) :: iold
3595       integer :: i,j,k,idum,i1,istart,iend,ierror,ierrcode
3596       real(kind=8) :: ctdif,dif
3597
3598       ctdif=10.
3599
3600       if (iseed.gt.mxio .or. iseed.lt.1) then
3601         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3602         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3603       endif
3604       do k=1,numch
3605        do j=2,nres-1
3606         do i=1,4
3607          vvar(i,j,k)=bvar(i,j,k,iseed)
3608         enddo
3609        enddo
3610       enddo
3611
3612
3613       do j=istart,iend
3614        do i=1,4
3615         vvar(i,j,1)=bvar(i,j,1,i1)
3616        enddo
3617       enddo
3618
3619       return
3620       end subroutine newconf_copy
3621 !-----------------------------------------------------------------------------
3622       subroutine newconf_residue_hairpin(idum,vvar,i1,fail)
3623
3624       use geometry_data
3625 !      use random, only: iran_num
3626       use MPI_data
3627       use compare, only:hairpin
3628
3629       include 'mpif.h'
3630 !      implicit real*8 (a-h,o-z)
3631 !      include 'DIMENSIONS'
3632 !      include 'COMMON.IOUNITS'
3633 !      include 'COMMON.CSA'
3634 !      include 'COMMON.BANK'
3635 !      include 'COMMON.CHAIN'
3636 !      include 'COMMON.GEO'
3637 !      include 'COMMON.VAR'
3638 !      real(kind=4) :: ran1,ran2
3639       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3640       integer,dimension(ntotal) :: iold
3641       integer :: nharp,iharp(4,nres/3),icipa(nres/3)
3642       logical :: fail,not_done
3643       integer :: idum,i,j,k,i1,iend,istart,iii,n_used,icount,iih,&
3644             ierror,ierrcode
3645       real(kind=8) :: ctdif,dif
3646
3647       ctdif=10.
3648
3649       fail=.false.
3650       if (iseed.gt.mxio .or. iseed.lt.1) then
3651         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3652         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3653       endif
3654       do k=1,numch
3655        do j=2,nres-1
3656         do i=1,4
3657          vvar(i,j,k)=bvar(i,j,k,iseed)
3658         enddo
3659        enddo
3660       enddo
3661       do k=1,numch
3662       do j=2,nres-1
3663         theta(j+1)=bvar(1,j,k,i1)
3664         phi(j+2)=bvar(2,j,k,i1)
3665         alph(j)=bvar(3,j,k,i1)
3666         omeg(j)=bvar(4,j,k,i1)
3667       enddo
3668       enddo
3669 !      call intout
3670       call chainbuild
3671       call hairpin(.false.,nharp,iharp)
3672
3673       if (nharp.eq.0) then
3674         fail=.true.
3675         return
3676       endif
3677
3678       n_used=0
3679
3680       DO III=1,NHARP
3681
3682       not_done = .true.
3683       icount=0
3684       do while (not_done)
3685         icount=icount+1
3686         iih=iran_num(1,nharp)
3687         do k=1,n_used
3688           if (iih.eq.icipa(k)) then
3689             iih=0
3690             goto 22
3691           endif
3692         enddo
3693         not_done=.false.
3694         n_used=n_used+1
3695         icipa(n_used)=iih
3696    22   continue
3697         not_done = not_done .and. icount.le.nharp
3698       enddo
3699
3700       if (iih.eq.0) then
3701         write (iout,*) "CHUJ NASTAPIL W NEWCONF_RESIDUE_HAIRPIN!!!!"
3702         fail=.true.
3703         return
3704       endif
3705
3706       istart=iharp(1,iih)+1
3707       iend=iharp(2,iih)
3708
3709 !dd      write (iout,*) "newconf_residue_hairpin: iih",iih,
3710 !dd     &   " istart",istart," iend",iend
3711
3712       do k=1,numch
3713       do j=istart,iend
3714        do i=1,4
3715          dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3716          if(dif.gt.180.) dif=360.-dif
3717          if(dif.gt.ctdif) goto 20
3718        enddo
3719       enddo
3720       enddo
3721       goto 10
3722    20 continue
3723       do k=1,numch
3724       do j=istart,iend
3725        do i=1,4
3726         vvar(i,j,k)=bvar(i,j,k,i1)
3727        enddo
3728       enddo
3729       enddo
3730 !      do j=1,numch
3731 !       do l=2,nres-1
3732 !        write (iout,'(4f8.3)') (rad2deg*vvar(i,l,j),i=1,4)
3733 !       enddo
3734 !      enddo
3735       return
3736    10 continue
3737       ENDDO
3738
3739       fail=.true.
3740
3741       return
3742       end subroutine newconf_residue_hairpin
3743 !-----------------------------------------------------------------------------
3744       subroutine gen_hairpin
3745
3746       use geometry_data
3747       use MD_data
3748       use compare, only:hairpin
3749 !      implicit real*8 (a-h,o-z)
3750 !      include 'DIMENSIONS'
3751 !      include 'COMMON.IOUNITS'
3752 !      include 'COMMON.CSA'
3753 !      include 'COMMON.BANK'
3754 !      include 'COMMON.CHAIN'
3755 !      include 'COMMON.GEO'
3756 !      include 'COMMON.VAR'
3757 !      include 'COMMON.HAIRPIN'
3758       integer :: i1,j,k,iters
3759
3760 !      write (iout,*) 'Entering GEN_HAIRPIN'
3761       do iters=1,nseed
3762         i1=is(iters)
3763         do k=1,numch
3764           do j=2,nres-1
3765             theta(j+1)=bvar(1,j,k,i1)
3766             phi(j+2)=bvar(2,j,k,i1)
3767             alph(j)=bvar(3,j,k,i1)
3768             omeg(j)=bvar(4,j,k,i1)
3769           enddo
3770         enddo
3771         call chainbuild
3772         call hairpin(.false.,nharp_seed(iters),iharp_seed(1,1,iters))
3773       enddo
3774
3775       nharp_tot=0
3776       do iters=1,nseed
3777         nharp_tot=nharp_tot+nharp_seed(iters)
3778         nharp_use(iters)=4*nharp_seed(iters)
3779         do j=1,nharp_seed(iters)
3780           iharp_use(0,j,iters)=4
3781           do k=1,4
3782             iharp_use(k,j,iters)=0
3783           enddo
3784         enddo
3785       enddo
3786
3787       write (iout,*) 'GEN_HAIRPIN: nharp_tot',nharp_tot
3788 !dd      do i=1,nseed
3789 !dd        write (iout,*) 'seed',i
3790 !dd        write (iout,*) 'nharp_seed',nharp_seed(i),
3791 !dd     &     ' nharp_use',nharp_use(i)
3792 !d        write (iout,*) 'iharp_seed, iharp_use'
3793 !d        do j=1,nharp_seed(i)
3794 !d          write (iout,'(7i3)') iharp_seed(1,j,i),iharp_seed(2,j,i),
3795 !d     &      (iharp_use(k,j,i),k=0,4)
3796 !d        enddo
3797 !dd      enddo
3798       return
3799       end subroutine gen_hairpin
3800 !-----------------------------------------------------------------------------
3801       subroutine select_frag(nn,nh,nl,ns,nb,i_csa)
3802
3803       use geometry_data
3804       use MD_data
3805       use compare_data
3806 !      implicit real*8 (a-h,o-z)
3807 !      include 'DIMENSIONS'
3808 !      include 'COMMON.IOUNITS'
3809 !      include 'COMMON.CSA'
3810 !      include 'COMMON.BANK'
3811 !      include 'COMMON.CHAIN'
3812 !      include 'COMMON.GEO'
3813 !      include 'COMMON.VAR'
3814 !      include 'COMMON.HAIRPIN'
3815 !      include 'COMMON.DISTFIT'
3816       character(len=50) :: linia
3817       integer :: isec(nres)
3818       integer :: i,j,i1,k,nn,nh,nl,ns,nb,i_csa,nl1,ns1
3819
3820       nn=0
3821       nh=0
3822       nl=0
3823       ns=0
3824       nb=0
3825 !d      write (iout,*) 'Entering select_frag'
3826       do i1=1,nbank
3827         do i=1,nres
3828           isec(i)=0
3829         enddo
3830         do k=1,numch
3831           do j=2,nres-1
3832             theta(j+1)=bvar(1,j,k,i1)
3833             phi(j+2)=bvar(2,j,k,i1)
3834             alph(j)=bvar(3,j,k,i1)
3835             omeg(j)=bvar(4,j,k,i1)
3836           enddo
3837         enddo
3838         call chainbuild
3839 !d        write (iout,*) ' -- ',i1,' -- '
3840         call secondary2(.false.)
3841 !
3842 ! bvar_frag nn==pair of nonlocal strands in beta sheet (loop>4)
3843 !               strands > 4 residues; used by N7 and N16
3844 !
3845         do j=1,nbfrag
3846 !
3847 !test 09/12/02 bfrag(2,j)-bfrag(1,j).gt.3
3848 !
3849               do i=bfrag(1,j),bfrag(2,j)
3850                 isec(i)=1
3851               enddo
3852               do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3853                 isec(i)=1
3854               enddo
3855
3856            if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
3857                 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
3858                 bfrag(2,j)-bfrag(1,j).gt.4 ) then
3859              nn=nn+1
3860
3861
3862              if (bfrag(3,j).lt.bfrag(4,j)) then
3863                write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3864                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
3865                                 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
3866              else
3867                write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3868                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
3869                                 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
3870
3871              endif
3872 !d             call write_pdb(i_csa*1000+nn+nh,linia,0d0)
3873
3874              bvar_frag(nn,1)=i1
3875              bvar_frag(nn,2)=4
3876              do i=1,4
3877                bvar_frag(nn,i+2)=bfrag(i,j)
3878              enddo
3879            endif
3880         enddo
3881
3882 !
3883 ! hvar_frag nh==helices; used by N8 and N9
3884 !
3885         do j=1,nhfrag
3886
3887               do i=hfrag(1,j),hfrag(2,j)
3888                 isec(i)=2
3889               enddo
3890
3891           if ( hfrag(2,j)-hfrag(1,j).gt.4 ) then
3892             nh=nh+1
3893
3894 !d            write(linia,'(a6,i3,a1,i3)')
3895 !d     &                "select",hfrag(1,j)-1,"-",hfrag(2,j)-1
3896 !d            call write_pdb(i_csa*1000+nn+nh,linia,0d0)
3897
3898             hvar_frag(nh,1)=i1
3899             hvar_frag(nh,2)=hfrag(1,j)
3900             hvar_frag(nh,3)=hfrag(2,j)
3901           endif
3902         enddo
3903
3904
3905 !v        write(iout,'(i4,1pe12.4,1x,1000i1)')
3906 !v     &         i1,bene(i1),(isec(i),i=1,nres)
3907 !v            write(linia,'(i4,1x,1000i1)')
3908 !v     &         i1,(isec(i),i=1,nres)
3909 !v            call write_pdb(i_csa*1000+i1,linia,bene(i1))
3910 !
3911 ! lvar_frag nl==loops; used by N14
3912 !
3913         i=1
3914         nl1=nl
3915         do while (i.lt.nres)
3916           if (isec(i).eq.0) then
3917            nl=nl+1
3918            lvar_frag(nl,1)=i1
3919            lvar_frag(nl,2)=i
3920            i=i+1
3921            do while (isec(i).eq.0.and.i.le.nres)
3922              i=i+1
3923            enddo
3924            lvar_frag(nl,3)=i-1
3925            if (lvar_frag(nl,3)-lvar_frag(nl,2).lt.1) nl=nl-1
3926           endif
3927           i=i+1
3928         enddo
3929 !d        write(iout,'(4i5)') (i,(lvar_frag(i,ii),ii=1,3),i=nl1+1,nl)
3930
3931 !
3932 ! svar_frag ns==an secondary structure element; used by N15
3933 !
3934         i=1
3935         ns1=ns
3936         do while (i.lt.nres)
3937           if (isec(i).gt.0) then
3938            ns=ns+1
3939            svar_frag(ns,1)=i1
3940            svar_frag(ns,2)=i
3941            i=i+1
3942            do while (isec(i).gt.0.and.isec(i-1).eq.isec(i) &
3943                                .and.i.le.nres)
3944              i=i+1
3945            enddo
3946            svar_frag(ns,3)=i-1
3947            if (svar_frag(ns,3)-svar_frag(ns,2).lt.1) ns=ns-1
3948           endif
3949           if (isec(i).eq.0) i=i+1
3950         enddo
3951 !d        write(iout,'(4i5)') (i,(svar_frag(i,ii),ii=1,3),i=ns1+1,ns)
3952
3953 !
3954 ! avar_frag nb==any pair of beta strands; used by N17
3955 !
3956         do j=1,nbfrag
3957              nb=nb+1
3958              avar_frag(nb,1)=i1
3959              do i=1,4
3960                avar_frag(nb,i+1)=bfrag(i,j)
3961              enddo
3962         enddo
3963
3964       enddo
3965
3966       return
3967       end subroutine select_frag
3968 !-----------------------------------------------------------------------------
3969 ! together.F
3970 !-----------------------------------------------------------------------------
3971       subroutine together
3972
3973 !  feeds tasks for parallel processing
3974       use MPI_data
3975       use geometry_data
3976       use control_data, only: vdisulf
3977       use energy_data
3978       use io, only:from_int,write_csa_pdb
3979 !      implicit real*8 (a-h,o-z)
3980 !      include 'DIMENSIONS'
3981       include 'mpif.h'
3982 !      real(kind=4) :: ran1,ran2
3983 !      include 'COMMON.CSA'
3984 !      include 'COMMON.BANK'
3985 !      include 'COMMON.IOUNITS'
3986 !      include 'COMMON.CHAIN'
3987 !      include 'COMMON.TIME1'
3988 !      include 'COMMON.SETUP'
3989 !      include 'COMMON.VAR'
3990 !      include 'COMMON.GEO'
3991 !      include 'COMMON.CONTROL'
3992 !      include 'COMMON.SBRIDGE'
3993       real(kind=4) :: tcpu
3994       real(kind=8) :: time_start,time_start_c,time0f,time0i
3995       logical :: ovrtim,sync_iter,timeout,flag,timeout1
3996       integer,dimension(mpi_status_size) :: muster
3997       real(kind=8),dimension(0:100) :: t100
3998       integer,dimension(mxio) :: indx
3999       real(kind=8),dimension(6*nres) :: xout    !(maxvar) (maxvar=6*maxres)
4000       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4001       integer,dimension(9) :: ind
4002       real(kind=8),dimension(2) :: cout
4003       real(kind=8),parameter :: rad=1.745329252d-2
4004
4005       integer :: i,m,j,jlee,nft,idum,nrmsdb,nrmsdb1,ierr,ierror,ierrcode,&
4006             ntrial,ntry,idum2,imax,idumm,nconfr,iconf,mm,k,im,nst,ifar,&
4007             iter,irecv,isent,iw_pdb,nft0i,nft00_c,nft00,ifrom,ij,&
4008             ireq,ireq2
4009       real(kind=8) :: adif,p_cut,cutdifr,rmsdbc1c,time1i,ctdif1,xctdif,&
4010                  time2i,tstart,tend1
4011 !ccccccccccccccccccccccccccccccccccccccccccccccc
4012       sync_iter=.true.   !el
4013       nft=0     !el
4014       time_start=0.0d0
4015       IF (ME.EQ.KING) THEN
4016
4017        time0f=MPI_WTIME()
4018        ilastnstep=1
4019        sync_iter=.false.
4020        numch=1
4021        nrmsdb=0
4022        nrmsdb1=0
4023        rmsdbc1c=rmsdbc1
4024        nstep=0
4025        call csa_read
4026        call make_array
4027
4028        if(iref.ne.0) call from_int(1,0,idum)
4029
4030 ! To minimize input conformation (bank conformation)
4031 ! Output to $mol.reminimized
4032        if (irestart.lt.0) then
4033         call read_bank(0,nft,cutdifr)
4034         if (irestart.lt.-10) then
4035          p_cut=nres*4.d0
4036          call prune_bank(p_cut)
4037          return
4038         endif
4039         call reminimize(jlee)
4040         return
4041        endif
4042
4043        if (irestart.eq.0) then
4044         call initial_write
4045         nbank=nconf
4046         ntbank=nconf
4047         if (ntbankm.eq.0) ntbank=0
4048         nstep=0
4049         nft=0
4050         do i=1,mxio
4051          ibank(i)=0
4052          jbank(i)=0
4053         enddo
4054        else
4055         call restart_write
4056 !!bankt        call read_bankt(jlee,nft,cutdifr)
4057         call read_bank(jlee,nft,cutdifr)
4058         call read_rbank(jlee,adif)
4059         if(iref.ne.0) call from_int(1,0,idum)
4060        endif
4061
4062        nstmax=nstmax+nstep
4063        ntrial=n1+n2+n3+n4+n5+n6+n7+n8
4064        ntry=ntrial+1
4065        ntry=ntry*nseed
4066
4067 ! ntrial : number of trial conformations per seed.
4068 ! ntry : total number of trial conformations including seed conformations.
4069
4070        idum2=-123
4071 !       imax=2**31-1
4072        imax=huge(0)
4073        ENDIF
4074
4075        call mpi_bcast(jend,1,mpi_integer,0,CG_COMM,ierr)
4076 !ccccccccccccccccccccccccccccccccccccccc
4077        do 300 jlee=1,jend
4078 !ccccccccccccccccccccccccccccccccccccccc
4079   331   continue
4080         IF (ME.EQ.KING) THEN
4081         if(sync_iter) goto 333
4082         idum=-  ran2(idum2)*imax
4083         if(jlee.lt.jstart) goto 300
4084
4085 ! Restart the random number generator for conformation generation
4086
4087         if(irestart.gt.0) then
4088          idum2=idum2+nstep
4089          if(idum2.le.0) idum2=-idum2+1
4090          idum=-  ran2(idum2)*imax
4091         endif
4092
4093         idumm=idum
4094         call vrndst(idumm)
4095
4096         open(icsa_seed,file=csa_seed,status="old")
4097         write(icsa_seed,*) "jlee : ",jlee
4098         close(icsa_seed)
4099
4100       call history_append
4101       write(icsa_history,*) "number of procs is ",nodes
4102       write(icsa_history,*) jlee,idum,idum2
4103       close(icsa_history)
4104
4105 !ccccccccccccccccccccccccccccccccccccccccccccccc
4106   333 icycle=0
4107
4108        call history_append
4109         write(icsa_history,*) "nbank is ",nbank
4110        close(icsa_history)
4111
4112       if(irestart.eq.1) goto 111
4113       if(irestart.eq.2) then
4114        icycle=0
4115        do i=1,nbank
4116         ibank(i)=1
4117        enddo
4118        do i=nbank+1,nbank+nconf
4119         ibank(i)=0
4120        enddo
4121       endif
4122
4123 !  start energy minimization
4124       nconfr=max0(nconf+nadd,nodes-1)
4125       if (sync_iter) nconf_in=0
4126 !  king-emperor - feed input and sort output
4127        write (iout,*) "NCONF_IN",nconf_in
4128        m=0
4129        if (nconf_in.gt.0) then
4130 ! al 7/2/00 - added possibility to read in some of the initial conformations
4131         do m=1,nconf_in
4132           read (intin,'(i5)',end=11,err=12) iconf
4133    12     continue
4134           write (iout,*) "write READ_ANGLES",iconf,m
4135           call read_angles(intin,*11)
4136           if (iref.eq.0) then
4137             mm=m
4138           else
4139             mm=m+1
4140           endif
4141           do j=2,nres-1
4142             dihang_in(1,j,1,mm)=theta(j+1)
4143             dihang_in(2,j,1,mm)=phi(j+2)
4144             dihang_in(3,j,1,mm)=alph(j)
4145             dihang_in(4,j,1,mm)=omeg(j)
4146           enddo
4147         enddo ! m
4148         goto 13
4149    11   write (iout,*) nconf_in," conformations requested, but only",&
4150          m-1," found in the angle file."
4151         nconf_in=m-1
4152    13   continue
4153         m=nconf_in
4154         write (iout,*) nconf_in,&
4155           " initial conformations have been read in."
4156        endif
4157        if (iref.eq.0) then
4158         if (nconfr.gt.nconf_in) then
4159           call make_ranvar(nconfr,m,idum)
4160           write (iout,*) nconfr-nconf_in,&
4161            " conformations have been generated randomly."
4162         endif
4163        else
4164         nconfr=nconfr*2
4165         call from_int(nconfr,m,idum)
4166 !       call from_pdb(nconfr,idum)
4167        endif
4168        write (iout,*) 'Exitted from make_ranvar nconfr=',nconfr
4169        write (*,*) 'Exitted from make_ranvar nconfr=',nconfr
4170        do m=1,nconfr
4171           write (iout,*) 'Initial conformation',m
4172           write(iout,'(8f10.4)') (rad2deg*dihang_in(1,j,1,m),j=2,nres-1)
4173           write(iout,'(8f10.4)') (rad2deg*dihang_in(2,j,1,m),j=2,nres-1)
4174           write(iout,'(8f10.4)') (rad2deg*dihang_in(3,j,1,m),j=2,nres-1)
4175           write(iout,'(8f10.4)') (rad2deg*dihang_in(4,j,1,m),j=2,nres-1)
4176        enddo
4177        write(iout,*)'Calling FEEDIN NCONF',nconfr
4178        time1i=MPI_WTIME()
4179        call feedin(nconfr,nft)
4180        write (iout,*) ' Time for first bank min.',MPI_WTIME()-time1i
4181        call  history_append
4182         write(icsa_history,*) jlee,nft,nbank
4183         write(icsa_history,851) (etot(i),i=1,nconfr)
4184         write(icsa_history,850) (rmsn(i),i=1,nconfr)
4185         write(icsa_history,850) (pncn(i),i=1,nconfr)
4186         write(icsa_history,*)
4187        close(icsa_history)
4188       ELSE
4189 ! To minimize input conformation (bank conformation)
4190 ! Output to $mol.reminimized
4191        if (irestart.lt.0) then
4192         call reminimize(jlee)
4193         return
4194        endif
4195        if (irestart.eq.1) goto 111
4196 !  soldier - perform energy minimization
4197  334   call minim_jlee
4198       ENDIF
4199
4200 !cccccccccccccccccccccccccccccccccc
4201 ! need to syncronize all procs
4202       call mpi_barrier(CG_COMM,ierr)
4203       if (ierr.ne.0) then
4204        print *, ' cannot synchronize MPI'
4205        stop
4206       endif
4207 !cccccccccccccccccccccccccccccccccc
4208
4209       IF (ME.EQ.KING) THEN
4210
4211 !      print *,"ok after minim"
4212       nstep=nstep+nconf
4213       if(irestart.eq.2) then
4214        nbank=nbank+nconf
4215 !      ntbank=ntbank+nconf
4216        if(ntbank.gt.ntbankm) ntbank=ntbankm
4217       endif
4218 !      print *,"ok before indexx"
4219       if(iref.eq.0) then
4220        call indexx(nconfr,etot,indx)
4221       else
4222 ! cc/al 7/6/00
4223        do k=1,nconfr
4224          indx(k)=k
4225        enddo
4226        call indexx(nconfr-nconf_in,rmsn(nconf_in+1),indx(nconf_in+1))
4227        do k=nconf_in+1,nconfr
4228          indx(k)=indx(k)+nconf_in
4229        enddo
4230 ! cc/al
4231 !       call indexx(nconfr,rmsn,indx)
4232       endif
4233 !      print *,"ok after indexx"
4234       do im=1,nconf
4235        m=indx(im)
4236        if (m.gt.mxio .or. m.lt.1) then
4237          write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,' M',m
4238          call mpi_abort(mpi_comm_world,ierror,ierrcode)
4239        endif
4240        jbank(im+nbank-nconf)=0
4241        bene(im+nbank-nconf)=etot(m)
4242        rene(im+nbank-nconf)=etot(m)
4243 !!bankt       btene(im)=etot(m)
4244 !
4245        brmsn(im+nbank-nconf)=rmsn(m)
4246        bpncn(im+nbank-nconf)=pncn(m)
4247        rrmsn(im+nbank-nconf)=rmsn(m)
4248        rpncn(im+nbank-nconf)=pncn(m)
4249        if (im+nbank-nconf.gt.mxio .or. im+nbank-nconf.lt.1) then
4250          write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,&
4251          ' NBANK',nbank,' NCONF',nconf,' IM+NBANK-NCONF',im+nbank-nconf
4252          call mpi_abort(mpi_comm_world,ierror,ierrcode)
4253        endif
4254        do k=1,numch
4255         do j=2,nres-1
4256          do i=1,4
4257           bvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
4258           rvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
4259 !!bankt          btvar(i,j,k,im)=dihang(i,j,k,m)
4260 !
4261          enddo
4262         enddo
4263        enddo
4264        if(iref.eq.1) then
4265         if(brmsn(im+nbank-nconf).gt.rmscut.or. &
4266            bpncn(im+nbank-nconf).lt.pnccut) ibank(im+nbank-nconf)=9
4267        endif
4268        if(vdisulf) then
4269            bvar_ns(im+nbank-nconf)=ns-2*nss
4270            k=0
4271            do i=1,ns
4272              j=1
4273              do while( iss(i).ne.ihpb(j)-nres .and. &
4274                        iss(i).ne.jhpb(j)-nres .and. j.le.nss)
4275               j=j+1
4276              enddo
4277              if (j.gt.nss) then
4278                k=k+1
4279                bvar_s(k,im+nbank-nconf)=iss(i)
4280              endif
4281            enddo
4282        endif
4283        bvar_nss(im+nbank-nconf)=nss
4284        do i=1,nss
4285            bvar_ss(1,i,im+nbank-nconf)=ihpb(i)
4286            bvar_ss(2,i,im+nbank-nconf)=jhpb(i)
4287        enddo
4288       enddo
4289       ENDIF
4290
4291   111 continue
4292
4293       IF (ME.EQ.KING) THEN
4294
4295       call find_max
4296       call find_min
4297
4298       call get_diff
4299       if(nbank.eq.nconf.and.irestart.eq.0) then
4300        adif=avedif
4301       endif
4302
4303       cutdif=adif/cut1
4304       ctdif1=adif/cut2
4305
4306 !d      print *,"adif,xctdif,cutdifr"
4307 !d      print *,adif,xctdif,cutdifr
4308        nst=ntotal/ntrial/nseed
4309        xctdif=(cutdif/ctdif1)**(-1.0/nst)
4310        if(irestart.ge.1) call estimate_cutdif(adif,xctdif,cutdifr)
4311 !       print *,"ok after estimate"
4312
4313       irestart=0
4314
4315        call write_rbank(jlee,adif,nft)
4316        call write_bank(jlee,nft)
4317 !!bankt       call write_bankt(jlee,nft)
4318 !       call write_bank1(jlee)
4319        call  history_append
4320         write(icsa_history,*) "xctdif: ", xctdif,nst,adif/cut1,ctdif1
4321         write(icsa_history,851) (bene(i),i=1,nbank)
4322         write(icsa_history,850) (brmsn(i),i=1,nbank)
4323         write(icsa_history,850) (bpncn(i),i=1,nbank)
4324         close(icsa_history)
4325   850 format(10f8.3)
4326   851 format(5e15.6)
4327
4328       ifar=nseed/4*3+1
4329       ifar=nseed+1
4330       ENDIF
4331
4332
4333       finished=.false.
4334       iter = 0
4335       irecv = 0
4336       isent =0
4337       ifrom= 0
4338       time0i=MPI_WTIME()
4339       time1i=time0i
4340       time_start_c=time0i
4341       if (.not.sync_iter) then
4342         time_start=time0i
4343         nft00=nft
4344       else
4345         sync_iter=.false.
4346       endif
4347       nft00_c=nft
4348       nft0i=nft
4349
4350 !cccccccccccccccccccccccccccccccccccccc
4351       do while (.not. finished)
4352 !cccccccccccccccccccccccccccccccccccccc
4353 !rc      print *,"iter ", iter,' isent=',isent
4354
4355       IF (ME.EQ.KING) THEN
4356 !  start energy minimization
4357
4358        if (isent.eq.0) then
4359 !  king-emperor - select seeds & make var & feed input
4360 !d        print *,'generating new conf',ntrial,MPI_WTIME()
4361         call select_is(nseed,ifar,idum)
4362
4363         open(icsa_seed,file=csa_seed,status="old")
4364         write(icsa_seed,39)  &
4365           jlee,icycle,nstep,(is(i),bene(is(i)),i=1,nseed)
4366         close(icsa_seed)
4367         call  history_append
4368         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
4369          ebmin,ebmax,nft,iuse,nbank,ntbank
4370         close(icsa_history)
4371
4372
4373
4374         call make_var(ntry,idum,iter)
4375 !d        print *,'new trial generated',ntrial,MPI_WTIME()
4376            time2i=MPI_WTIME()
4377            write (iout,'(a20,i4,f12.2)') &
4378              'Time for make trial',iter+1,time2i-time1i
4379        endif
4380
4381 !rc        write(iout,*)'1:Calling FEEDIN NTRY',NTRY,' ntrial',ntrial
4382 !rc        call feedin(ntry,nft)
4383
4384        isent=isent+1
4385        if (isent.ge.nodes.or.iter.gt.0)  then
4386 !t            print *,'waiting ',MPI_WTIME()
4387             irecv=irecv+1
4388             call recv(0,ifrom,xout,eout,ind,timeout)
4389 !t            print *,'   ',irecv,' received from',ifrom,MPI_WTIME()
4390        else
4391             ifrom=ifrom+1
4392        endif
4393
4394 !t            print *,'sending to',ifrom,MPI_WTIME()
4395        call send(isent,ifrom,iter)
4396 !t            print *,isent,' sent ',MPI_WTIME()
4397
4398 ! store results -----------------------------------------------
4399        if (isent.ge.nodes.or.iter.gt.0)  then
4400          nft=nft+ind(3)
4401          movernx(irecv)=iabs(ind(5))
4402          call getx(ind,xout,eout,cout,rad,iw_pdb,irecv)
4403          if(vdisulf) then
4404              nss_out(irecv)=nss
4405              do i=1,nss
4406                iss_out(i,irecv)=ihpb(i)
4407                jss_out(i,irecv)=jhpb(i)
4408              enddo
4409          endif
4410          if(iw_pdb.gt.0) &
4411                 call write_csa_pdb(xout,eout,nft,irecv,iw_pdb)
4412        endif
4413 !--------------------------------------------------------------
4414        if (isent.eq.ntry) then
4415            time1i=MPI_WTIME()
4416            write (iout,'(a18,f12.2,a14,f10.2)') &
4417              'Nonsetup time     ',time1i-time_start_c,&
4418              ' sec, Eval/s =',(nft-nft00_c)/(time1i-time_start_c)
4419            write (iout,'(a14,i4,f12.2,a14,f10.2)') &
4420              'Time for iter ',iter+1,time1i-time0i,&
4421              ' sec, Eval/s =',(nft-nft0i)/(time1i-time0i)
4422            time0i=time1i
4423            nft0i=nft
4424            cutdif=cutdif*xctdif
4425            if(cutdif.lt.ctdif1) cutdif=ctdif1
4426            if (iter.eq.0) then
4427               print *,'UPDATING ',ntry-nodes+1,irecv
4428               write(iout,*) 'UPDATING ',ntry-nodes+1
4429               iter=iter+1
4430 !----------------- call update(ntry-nodes+1) -------------------
4431               nstep=nstep+ntry-nseed-(nodes-1)
4432               call refresh_bank(ntry-nodes+1)
4433 !!bankt              call refresh_bankt(ntry-nodes+1)
4434            else
4435 !----------------- call update(ntry) ---------------------------
4436               iter=iter+1
4437               print *,'UPDATING ',ntry,irecv
4438               write(iout,*) 'UPDATING ',ntry
4439               nstep=nstep+ntry-nseed
4440               call refresh_bank(ntry)
4441 !!bankt              call refresh_bankt(ntry)
4442            endif
4443 !-----------------------------------------------------------------
4444
4445            call write_bank(jlee,nft)
4446 !!bankt           call write_bankt(jlee,nft)
4447            call find_min
4448
4449            time1i=MPI_WTIME()
4450            write (iout,'(a20,i4,f12.2)') &
4451              'Time for refresh ',iter,time1i-time0i
4452
4453            if(ebmin.lt.estop) finished=.true.
4454            if(icycle.gt.icmax) then
4455                call write_bank1(jlee)
4456                do i=1,nbank
4457                  ibank(i)=2
4458                  ibank(i)=1
4459                enddo
4460                nbank=nbank+nconf
4461                if(nbank.gt.1000) then
4462                    finished=.true.
4463                else
4464 !rc                   goto 333
4465                    sync_iter=.true.
4466                endif
4467            endif
4468            if(nstep.gt.nstmax) finished=.true.
4469
4470            if(finished.or.sync_iter) then
4471             do ij=1,nodes-1
4472               call recv(1,ifrom,xout,eout,ind,timeout)
4473               if (timeout) then
4474                 nstep=nstep+ij-1
4475                 print *,'ERROR worker is not responding'
4476                 write(iout,*) 'ERROR worker is not responding'
4477                 time1i=MPI_WTIME()-time_start_c
4478                 print *,'End of cycle, master time for ',iter,' iters ',&
4479                    time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
4480                 write (iout,*) 'End of cycle, master time for ',iter,&
4481                    ' iters ',time1i,' sec'
4482                 write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
4483                 print *,'UPDATING ',ij-1
4484                 write(iout,*) 'UPDATING ',ij-1
4485                 call flush(iout)
4486                 call refresh_bank(ij-1)
4487 !!bankt                call refresh_bankt(ij-1)
4488                 goto 1002
4489               endif
4490 !              print *,'node ',ifrom,' finished ',ij,nft
4491               write(iout,*) 'node ',ifrom,' finished ',ij,nft
4492               call flush(iout)
4493               nft=nft+ind(3)
4494               movernx(ij)=iabs(ind(5))
4495               call getx(ind,xout,eout,cout,rad,iw_pdb,ij)
4496               if(vdisulf) then
4497                nss_out(ij)=nss
4498                do i=1,nss
4499                  iss_out(i,ij)=ihpb(i)
4500                  jss_out(i,ij)=jhpb(i)
4501                enddo
4502               endif
4503               if(iw_pdb.gt.0) &
4504                 call write_csa_pdb(xout,eout,nft,ij,iw_pdb)
4505             enddo
4506             nstep=nstep+nodes-1
4507 !rc            print *,'---------bcast finished--------',finished
4508             time1i=MPI_WTIME()-time_start_c
4509             print *,'End of cycle, master time for ',iter,' iters ',&
4510                    time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
4511             write (iout,*) 'End of cycle, master time for ',iter,&
4512                    ' iters ',time1i,' sec'
4513             write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
4514
4515 !timeout            call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
4516 !timeout            call mpi_bcast(sync_iter,1,mpi_logical,0,
4517 !timeout     &                                              CG_COMM,ierr)
4518             do ij=1,nodes-1
4519                tstart=MPI_WTIME()
4520                call mpi_issend(finished,1,mpi_logical,ij,idchar,&
4521                    CG_COMM,ireq,ierr)
4522                call mpi_issend(sync_iter,1,mpi_logical,ij,idchar,&
4523                    CG_COMM,ireq2,ierr)
4524                flag=.false.
4525                timeout1=.false.
4526                do while(.not. (flag .or. timeout1))
4527                  call MPI_TEST(ireq2,flag,muster,ierr)
4528                  tend1=MPI_WTIME()
4529                  if(tend1-tstart.gt.60) then
4530                   print *,'ERROR worker ',ij,' is not responding'
4531                   write(iout,*) 'ERROR worker ',ij,' is not responding'
4532                   timeout1=.true.
4533                  endif
4534                enddo
4535                if(timeout1) then
4536                 write(iout,*) 'worker ',ij,' NOT OK ',tend1-tstart
4537                 timeout=.true.
4538                else
4539                 write(iout,*) 'worker ',ij,' OK ',tend1-tstart
4540                endif
4541             enddo
4542             print *,'UPDATING ',nodes-1,ij
4543             write(iout,*) 'UPDATING ',nodes-1
4544             call refresh_bank(nodes-1)
4545 !!bankt            call refresh_bankt(nodes-1)
4546  1002       continue
4547             call write_bank(jlee,nft)
4548 !!bankt            call write_bankt(jlee,nft)
4549             call find_min
4550
4551             do i=0,mxmv
4552               do j=1,3
4553                 nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
4554                 nstatnx(i,j)=0
4555               enddo
4556             enddo
4557
4558             write(iout,*)'### Total stats:'
4559             do i=0,mxmv
4560              if(nstatnx_tot(i,1).ne.0) then
4561               if (i.le.9) then
4562               write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4563               '### N',i,' total=',nstatnx_tot(i,1),&
4564             ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',&
4565              (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
4566               else
4567               write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4568               '###N',i,' total=',nstatnx_tot(i,1),&
4569             ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',&
4570              (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
4571               endif
4572              else
4573               if (i.le.9) then
4574               write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4575                 '### N',i,' total=',nstatnx_tot(i,1),&
4576                 ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),&
4577                 ' %acc',0.0
4578               else
4579               write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4580                 '###N',i,' total=',nstatnx_tot(i,1),&
4581                 ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),&
4582                 ' %acc',0.0
4583               endif
4584              endif
4585             enddo
4586
4587            endif
4588            if(sync_iter) goto 331
4589
4590    39 format(2i3,i7,5(i4,f8.3,1x),19(/,13x,5(i4,f8.3,1x)))
4591    40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
4592    43 format(10i8)
4593    44 format('jlee =',i3,':',4f10.1,' E =',f8.3,i7,i10)
4594
4595            isent=0
4596            irecv=0
4597        endif
4598       ELSE
4599 !  soldier - perform energy minimization
4600         call minim_jlee
4601         print *,'End of minim, proc',me,'time ',MPI_WTIME()-time_start
4602         write (iout,*) 'End of minim, proc',me,'time ',&
4603                         MPI_WTIME()-time_start
4604         call flush(iout)
4605 !timeout        call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
4606 !timeout        call mpi_bcast(sync_iter,1,mpi_logical,0,CG_COMM,ierr)
4607          call mpi_recv(finished,1,mpi_logical,0,idchar,&
4608                        CG_COMM,muster,ierr)
4609          call mpi_recv(sync_iter,1,mpi_logical,0,idchar,&
4610                        CG_COMM,muster,ierr)
4611         if(sync_iter) goto 331
4612       ENDIF
4613
4614 !cccccccccccccccccccccccccccccccccccccc
4615       enddo
4616 !cccccccccccccccccccccccccccccccccccccc
4617
4618       IF (ME.EQ.KING) THEN
4619         call  history_append
4620         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
4621         ebmin,ebmax,nft,iuse,nbank,ntbank
4622
4623         write(icsa_history,44) jlee,0.0,0.0,0.0,&
4624          0.0,ebmin,nstep,nft
4625         write(icsa_history,*)
4626        close(icsa_history)
4627
4628        time1i=MPI_WTIME()-time_start
4629        print *,'End of RUN, master time ',&
4630                    time1i,'sec, Eval/s ',(nft-nft00)/time1i
4631        write (iout,*) 'End of RUN, master time  ',&
4632                    time1i,' sec'
4633        write (iout,*) 'Total eval/s ',(nft-nft00)/time1i
4634
4635        if(timeout) then
4636         write(iout,*) '!!!! ERROR worker was not responding'
4637         write(iout,*) '!!!! cannot finish work normally'
4638         write(iout,*) 'Processor0 is calling MPI_ABORT'
4639         print *,'!!!! ERROR worker was not responding'
4640         print *,'!!!! cannot finish work normally'
4641         print *,'Processor0 is calling MPI_ABORT'
4642         call flush(iout)
4643         call mpi_abort(mpi_comm_world, 111, ierr)
4644        endif
4645       ENDIF
4646
4647 !ccccccccccccccccccccccccccccc
4648   300 continue
4649 !ccccccccccccccccccccccccccccc
4650
4651       return
4652       end subroutine together
4653 !-----------------------------------------------------------------------------
4654       subroutine feedin(nconf,nft)
4655
4656       use MPI_data
4657       use geometry_data, only:nvar
4658       use io, only:write_csa_pdb
4659 !  sends out starting conformations and receives results of energy minimization
4660 !      implicit real*8 (a-h,o-z)
4661 !      include 'DIMENSIONS'
4662 !      include 'COMMON.VAR'
4663 !      include 'COMMON.IOUNITS'
4664 !      include 'COMMON.CONTROL'
4665       include 'mpif.h'
4666       real(kind=8),dimension(6*nres) :: xin,xout        !(maxvar) (maxvar=6*maxres)
4667       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4668       real(kind=8),dimension(2) :: cout
4669       integer,dimension(9) :: ind
4670       integer,dimension(12) :: info
4671       integer,dimension(mpi_status_size) :: muster
4672 !      include 'COMMON.SETUP'
4673       real(kind=8),parameter :: rad=1.745329252d-2
4674       integer :: j,nconf,nft,mm,n,ierror,ierrcode,ierr,iw_pdb,&
4675             man
4676
4677       print *,'FEEDIN: NCONF=',nconf
4678       mm=0
4679 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4680       if (nconf .lt. nodes-1) then
4681         write (*,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',&
4682          nconf,nodes-1
4683         write (iout,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',&
4684          nconf,nodes-1
4685         call mpi_abort(mpi_comm_world,ierror,ierrcode)
4686       endif
4687       do n=1,nconf
4688 !  pull out external and internal variables for next start
4689         call putx(xin,n,rad)
4690 !        write (iout,*) 'XIN from FEEDIN N=',n
4691 !        write(iout,'(8f10.4)') (xin(j),j=1,nvar)
4692         mm=mm+1
4693         if (mm.lt.nodes) then
4694 !  feed task to soldier
4695 !       print *, ' sending input for start # ',n
4696          info(1)=n
4697          info(2)=-1
4698          info(3)=0
4699          info(4)=0
4700          info(5)=0
4701          info(6)=0
4702          call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
4703                         ierr)
4704          call mpi_send(xin,nvar,mpi_double_precision,mm,&
4705                         idreal,CG_COMM,ierr)
4706         else
4707 !  find an available soldier
4708          call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,&
4709                        CG_COMM,muster,ierr)
4710 !        print *, ' receiving output from start # ',ind(1)
4711          man=muster(mpi_source)
4712 !  receive final energies and variables
4713          nft=nft+ind(3)
4714          call mpi_recv(eout,1,mpi_double_precision,&
4715                        man,idreal,CG_COMM,muster,ierr)
4716 !         print *,eout
4717 #ifdef CO_BIAS
4718          call mpi_recv(co,1,mpi_double_precision,&
4719                        man,idreal,CG_COMM,muster,ierr)
4720          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
4721 #endif
4722          call mpi_recv(xout,nvar,mpi_double_precision,&
4723                        man,idreal,CG_COMM,muster,ierr)
4724 !         print *,nvar , ierr
4725 !  feed next task to soldier
4726 !        print *, ' sending input for start # ',n
4727          info(1)=n
4728          info(2)=-1
4729          info(3)=0
4730          info(4)=0
4731          info(5)=0
4732          info(6)=0
4733          info(7)=0
4734          info(8)=0
4735          info(9)=0
4736          call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,&
4737                         ierr)
4738          call mpi_send(xin,nvar,mpi_double_precision,man,&
4739                         idreal,CG_COMM,ierr)
4740 !  retrieve latest results
4741          call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
4742          if(iw_pdb.gt.0) &
4743               call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
4744         endif
4745       enddo
4746 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4747 !  no more input
4748 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4749       do j=1,nodes-1
4750 !  wait for a soldier
4751        call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,&
4752                      CG_COMM,muster,ierr)
4753 !rc       if (ierr.ne.0) go to 30
4754 !      print *, ' receiving output from start # ',ind(1)
4755        man=muster(mpi_source)
4756 !  receive final energies and variables
4757        nft=nft+ind(3)
4758        call mpi_recv(eout,1,&
4759                      mpi_double_precision,man,idreal,&
4760                      CG_COMM,muster,ierr)
4761 !       print *,eout
4762 #ifdef CO_BIAS
4763          call mpi_recv(co,1,mpi_double_precision,&
4764                        man,idreal,CG_COMM,muster,ierr)
4765          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
4766 #endif
4767 !rc       if (ierr.ne.0) go to 30
4768        call mpi_recv(xout,nvar,mpi_double_precision,&
4769                      man,idreal,CG_COMM,muster,ierr)
4770 !       print *,nvar , ierr
4771 !rc       if (ierr.ne.0) go to 30
4772 !  halt soldier
4773        info(1)=0
4774        info(2)=-1
4775        info(3)=0
4776        info(4)=0
4777        info(5)=0
4778        info(6)=0
4779        call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,&
4780                       ierr)
4781 !  retrieve results
4782        call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
4783        if(iw_pdb.gt.0) &
4784                 call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
4785       enddo
4786 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4787       return
4788    10 print *, ' dispatching error'
4789       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4790       return
4791    20 print *, ' communication error'
4792       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4793       return
4794    30 print *, ' receiving error'
4795       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4796
4797       return
4798       end subroutine feedin
4799 !-----------------------------------------------------------------------------
4800       subroutine getx(ind,xout,eout,cout,rad,iw_pdb,k)
4801
4802       use geometry_data
4803       use energy_data
4804       use compare, only: contact_fract
4805       use MPI_data
4806       include 'mpif.h'
4807 !  receives and stores data from soldiers
4808 !      implicit real*8 (a-h,o-z)
4809 !      include 'DIMENSIONS'
4810 !      include 'COMMON.IOUNITS'
4811 !      include 'COMMON.CSA'
4812 !      include 'COMMON.BANK'
4813 !      include 'COMMON.VAR'
4814 !      include 'COMMON.CHAIN'
4815 !      include 'COMMON.CONTACTS'
4816       integer,dimension(9) :: ind
4817       real(kind=8),dimension(6*nres) :: xout    !(maxvar) (maxvar=6*maxres)
4818       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4819 !jlee
4820       real(kind=8) :: przes(3),obr(3,3),cout(2)
4821       logical :: non_conv
4822       integer :: iw_pdb,k,j,ierror,ierrcode
4823       real(kind=8) :: rad,co
4824 !jlee
4825       iw_pdb=2
4826       if (k.gt.mxio .or. k.lt.1) then
4827         write (iout,*) &
4828          'ERROR - dimensions of ANGMIN have been exceeded K=',k
4829         call mpi_abort(mpi_comm_world,ierror,ierrcode)
4830       endif
4831 !  store ind()
4832       do j=1,9
4833        indb(k,j)=ind(j)
4834       enddo
4835 !  store energies
4836       etot(k)=eout(1)
4837 !  retrieve dihedral angles etc
4838       call var_to_geom(nvar,xout)
4839       do j=2,nres-1
4840         dihang(1,j,1,k)=theta(j+1)
4841         dihang(2,j,1,k)=phi(j+2)
4842         dihang(3,j,1,k)=alph(j)
4843         dihang(4,j,1,k)=omeg(j)
4844       enddo
4845       dihang(2,nres-1,1,k)=0.0d0
4846 !jlee
4847       if(iref.eq.0) then
4848        iw_pdb=1
4849 !d       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4)')
4850 !d     &      ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' mv ',
4851 !d     &      ind(5),ind(4)
4852        return
4853       endif
4854       call chainbuild
4855 !     call dihang_to_c(dihang(1,1,1,k))
4856 !     call fitsq(rms,c(1,1),crefjlee(1,1),nres,przes,obr,non_conv)
4857 !     call fitsq(rms,c(1,2),crefjlee(1,2),nres-1,przes,obr,non_conv)
4858 !           call fitsq(rms,c(1,nstart_seq),crefjlee(1,nstart_sup),
4859 !    &                 nsup,przes,obr,non_conv)
4860 !      rmsn(k)=dsqrt(rms)
4861
4862        call rmsd_csa(rmsn(k))
4863        call contact(.false.,ncont,icont,co)
4864        pncn(k)=contact_fract(ncont,ncont_ref,icont,icont_ref)
4865
4866 !d       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a5
4867 !d     &      ,0pf5.2,a5,f5.1,a,f6.3,a4,i3,i4)')
4868 !d     &    ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' rms ',
4869 !d     &    rmsn(k),' %NC ',pncn(k)*100,' cont.order',co,' mv ',
4870 !d     &    ind(5),ind(4)
4871
4872
4873       if (rmsn(k).gt.rmscut.or.pncn(k).lt.pnccut) iw_pdb=0
4874       return
4875       end subroutine getx
4876 !-----------------------------------------------------------------------------
4877       subroutine putx(xin,n,rad)
4878
4879       use geometry_data
4880 !  gets starting variables
4881 !      implicit real*8 (a-h,o-z)
4882 !      include 'DIMENSIONS'
4883 !      include 'COMMON.CSA'
4884 !      include 'COMMON.BANK'
4885 !      include 'COMMON.VAR'
4886 !      include 'COMMON.CHAIN'
4887 !      include 'COMMON.IOUNITS'
4888       integer :: n,m,j
4889       real(kind=8),dimension(6*nres) :: xin     !(maxvar) (maxvar=6*maxres)
4890       real(kind=8) :: rad
4891
4892 !  pull out starting values for variables
4893 !       write (iout,*)'PUTX: N=',n
4894       do m=1,numch
4895 !        write (iout,'(8f10.4)') (dihang_in(1,j,m,n),j=2,nres-1)
4896 !        write (iout,'(8f10.4)') (dihang_in(2,j,m,n),j=2,nres-1)
4897 !        write (iout,'(8f10.4)') (dihang_in(3,j,m,n),j=2,nres-1)
4898 !        write (iout,'(8f10.4)') (dihang_in(4,j,m,n),j=2,nres-1)
4899        do j=2,nres-1
4900         theta(j+1)=dihang_in(1,j,m,n)
4901         phi(j+2)=dihang_in(2,j,m,n)
4902         alph(j)=dihang_in(3,j,m,n)
4903         omeg(j)=dihang_in(4,j,m,n)
4904        enddo
4905       enddo
4906 !  set up array of variables
4907       call geom_to_var(nvar,xin)
4908 !       write (iout,*) 'xin in PUTX N=',n
4909 !       call intout
4910 !       write (iout,'(8f10.4)') (xin(i),i=1,nvar)
4911       return
4912       end subroutine putx
4913 !-----------------------------------------------------------------------------
4914       subroutine putx2(xin,iff,n)
4915
4916       use geometry_data
4917 !  gets starting variables
4918 !      implicit real*8 (a-h,o-z)
4919 !      include 'DIMENSIONS'
4920 !      include 'COMMON.CSA'
4921 !      include 'COMMON.BANK'
4922 !      include 'COMMON.VAR'
4923 !      include 'COMMON.CHAIN'
4924 !      include 'COMMON.IOUNITS'
4925       integer :: n,m,j,i
4926       real(kind=8),dimension(6*nres) :: xin     !(maxvar) (maxvar=6*maxres)
4927       integer,dimension(nres) :: iff    !(maxres)
4928
4929 !  pull out starting values for variables
4930       do m=1,numch
4931        do j=2,nres-1
4932         theta(j+1)=dihang_in2(1,j,m,n)
4933         phi(j+2)=dihang_in2(2,j,m,n)
4934         alph(j)=dihang_in2(3,j,m,n)
4935         omeg(j)=dihang_in2(4,j,m,n)
4936        enddo
4937       enddo
4938 !  set up array of variables
4939       call geom_to_var(nvar,xin)
4940
4941       do i=1,nres
4942         iff(i)=iff_in(i,n)
4943       enddo
4944       return
4945       end subroutine putx2
4946 !-----------------------------------------------------------------------------
4947       subroutine prune_bank(p_cut)
4948
4949       use MPI_data
4950 !      implicit real*8 (a-h,o-z)
4951 !      include 'DIMENSIONS'
4952       include 'mpif.h'
4953 !      include 'COMMON.CSA'
4954 !      include 'COMMON.BANK'
4955 !      include 'COMMON.IOUNITS'
4956 !      include 'COMMON.CHAIN'
4957 !      include 'COMMON.TIME1'
4958 !      include 'COMMON.SETUP'
4959       integer :: k,j,i,m,ip,nprune
4960       real(kind=8) :: p_cut,diff,ddmin
4961 !---------------------------
4962 ! This subroutine prunes bank conformations using p_cut
4963 !---------------------------
4964
4965       nprune=0
4966       nprune=nprune+1
4967       m=1
4968       do k=1,numch
4969        do j=2,nres-1
4970         do i=1,4
4971          dihang(i,j,k,nprune)=bvar(i,j,k,m)
4972         enddo
4973        enddo
4974       enddo
4975       bene(nprune)=bene(m)
4976       brmsn(nprune)=brmsn(m)
4977       bpncn(nprune)=bpncn(m)
4978
4979       do m=2,nbank
4980        ddmin=9.d190
4981        do ip=1,nprune
4982         call get_diff12(dihang(1,1,1,ip),bvar(1,1,1,m),diff)
4983         if(diff.lt.p_cut) goto 100
4984         if(diff.lt.ddmin) ddmin=diff
4985        enddo
4986        nprune=nprune+1
4987        do k=1,numch
4988         do j=2,nres-1
4989          do i=1,4
4990           dihang(i,j,k,nprune)=bvar(i,j,k,m)
4991          enddo
4992         enddo
4993        enddo
4994        bene(nprune)=bene(m)
4995        brmsn(nprune)=brmsn(m)
4996        bpncn(nprune)=bpncn(m)
4997   100  continue
4998        write (iout,*) 'Pruning :',m,nprune,p_cut,ddmin
4999       enddo
5000       nbank=nprune
5001       print *, 'Pruning :',m,nprune,p_cut
5002       call write_bank(0,0)
5003
5004       return
5005       end subroutine prune_bank
5006 !-----------------------------------------------------------------------------
5007       subroutine reminimize(jlee)
5008
5009       use MPI_data
5010 !      implicit real*8 (a-h,o-z)
5011 !      include 'DIMENSIONS'
5012       include 'mpif.h'
5013 !      include 'COMMON.CSA'
5014 !      include 'COMMON.BANK'
5015 !      include 'COMMON.IOUNITS'
5016 !      include 'COMMON.CHAIN'
5017 !      include 'COMMON.TIME1'
5018 !      include 'COMMON.SETUP'
5019       integer :: i,j,k,jlee,index,nft,ntry
5020 !---------------------------
5021 ! This subroutine re-minimizes bank conformations:
5022 !---------------------------
5023
5024        ntry=nbank
5025
5026        call find_max
5027        call find_min
5028
5029        if (me.eq.king) then
5030         open(icsa_history,file=csa_history,status="old")
5031          write(icsa_history,*) "Re-minimization",nodes,"nodes"
5032          write(icsa_history,851) (bene(i),i=1,nbank)
5033          write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
5034           ebmin,ebmax,nft,iuse,nbank,ntbank
5035         close(icsa_history)
5036         do index=1,ntry
5037          do k=1,numch
5038           do j=2,nres-1
5039            do i=1,4
5040             dihang_in(i,j,k,index)=bvar(i,j,k,index)
5041            enddo
5042           enddo
5043          enddo
5044         enddo
5045         nft=0
5046         call feedin(ntry,nft)
5047        else
5048         call minim_jlee
5049        endif
5050
5051        call find_max
5052        call find_min
5053
5054        if (me.eq.king) then
5055         do i=1,ntry
5056          call replace_bvar(i,i)
5057         enddo
5058         open(icsa_history,file=csa_history,status="old")
5059          write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
5060           ebmin,ebmax,nft,iuse,nbank,ntbank
5061          write(icsa_history,851) (bene(i),i=1,nbank)
5062         close(icsa_history)
5063         call write_bank_reminimized(jlee,nft)
5064        endif
5065
5066    40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
5067   851 format(5e15.6)
5068   850 format(5e15.10)
5069 !  850 format(10f8.3)
5070
5071       return
5072       end subroutine reminimize
5073 !-----------------------------------------------------------------------------
5074       subroutine send(n,mm,it)
5075
5076       use MPI_data
5077       use geometry_data, only: nvar
5078       use control_data, only: vdisulf
5079 !  sends out starting conformation for minimization
5080 !      implicit real*8 (a-h,o-z)
5081 !      include 'DIMENSIONS'
5082 !      include 'COMMON.VAR'
5083 !      include 'COMMON.IOUNITS'
5084 !      include 'COMMON.CONTROL'
5085 !      include 'COMMON.BANK'
5086 !      include 'COMMON.CHAIN'
5087       include 'mpif.h'
5088       real(kind=8),dimension(6*nres) :: xin,xout,xin2   !(maxvar) (maxvar=6*maxres)
5089       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
5090       real(kind=8),dimension(2) :: cout
5091       integer,dimension(9) :: ind
5092       integer,dimension(nres) :: iff    !(maxres)
5093       integer,dimension(12) :: info
5094       integer,dimension(mpi_status_size) :: muster
5095 !      include 'COMMON.SETUP'
5096       real(kind=8),parameter :: rad=1.745329252d-2
5097       integer :: n,mm,it,ierr
5098
5099       if (isend2(n).eq.0) then
5100 !  pull out external and internal variables for next start
5101         call putx(xin,n,rad)
5102         info(1)=n
5103         info(2)=it
5104         info(3)=movenx(n)
5105         info(4)=nss_in(n)
5106         info(5)=parent(1,n)
5107         info(6)=parent(2,n)
5108
5109         if (movenx(n).eq.14.or.movenx(n).eq.17) then
5110           info(7)=idata(1,n)
5111           info(8)=idata(2,n)
5112         else if (movenx(n).eq.16) then
5113           info(7)=idata(1,n)
5114           info(8)=idata(2,n)
5115           info(10)=idata(3,n)
5116           info(11)=idata(4,n)
5117           info(12)=idata(5,n)
5118         else
5119          info(7)=0
5120          info(8)=0
5121          info(10)=0
5122          info(11)=0
5123          info(12)=0
5124         endif
5125
5126         if (movenx(n).eq.15) then
5127          info(9)=parent(3,n)
5128         else
5129          info(9)=0
5130         endif
5131         call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
5132                         ierr)
5133         call mpi_send(xin,nvar,mpi_double_precision,mm,&
5134                         idreal,CG_COMM,ierr)
5135       else
5136 !  distfit & minimization for n7 move
5137         info(1)=-n
5138         info(2)=it
5139         info(3)=movenx(n)
5140         info(4)=nss_in(n)
5141         info(5)=parent(1,n)
5142         info(6)=parent(2,n)
5143         info(7)=0
5144         info(8)=0
5145         info(9)=0
5146         call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
5147                         ierr)
5148         call putx2(xin,iff,isend2(n))
5149         call mpi_send(xin,nvar,mpi_double_precision,mm,&
5150                         idreal,CG_COMM,ierr)
5151         call mpi_send(iff,nres,mpi_integer,mm,&
5152                         idint,CG_COMM,ierr)
5153         call putx(xin2,n,rad)
5154         call mpi_send(xin2,nvar,mpi_double_precision,mm,&
5155                         idreal,CG_COMM,ierr)
5156       endif
5157       if (vdisulf.and.nss_in(n).ne.0) then
5158         call mpi_send(iss_in(1,n),nss_in(n),mpi_integer,mm,&
5159                         idint,CG_COMM,ierr)
5160         call mpi_send(jss_in(1,n),nss_in(n),mpi_integer,mm,&
5161                         idint,CG_COMM,ierr)
5162       endif
5163       return
5164       end subroutine send
5165 !-----------------------------------------------------------------------------
5166       subroutine recv(ihalt,man,xout,eout,ind,tout)
5167
5168       use MPI_data
5169       use energy_data
5170       use geometry_data, only: nvar
5171       use control_data, only: vdisulf
5172 !  receives results of energy minimization
5173 !      implicit real*8 (a-h,o-z)
5174 !      include 'DIMENSIONS'
5175 !      include 'COMMON.VAR'
5176 !      include 'COMMON.IOUNITS'
5177 !      include 'COMMON.CONTROL'
5178 !      include 'COMMON.SBRIDGE'
5179 !      include 'COMMON.BANK'
5180 !      include 'COMMON.CHAIN'
5181       include 'mpif.h'
5182       real(kind=8),dimension(6*nres) :: xin,xout        !(maxvar) (maxvar=6*maxres)
5183       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
5184       real(kind=8),dimension(2) :: cout
5185       integer,dimension(9) :: ind
5186       integer,dimension(12) :: info
5187       integer,dimension(mpi_status_size) :: muster
5188 !      include 'COMMON.SETUP'
5189       logical :: tout,flag
5190       real(kind=8) :: tstart,tend1
5191       real(kind=8),parameter :: twait=600.0d0
5192       integer :: ihalt,man,ierr
5193
5194 !  find an available soldier
5195        tout=.false.
5196        flag=.false.
5197        tstart=MPI_WTIME()
5198        do while(.not. (flag .or. tout))
5199          call MPI_IPROBE(mpi_any_source,idint,CG_COMM,flag, &
5200                   muster,ierr)
5201          tend1=MPI_WTIME()
5202          if(tend1-tstart.gt.twait .and. ihalt.eq.1) tout=.true.
5203 !_error         if(tend1-tstart.gt.twait) tout=.true.
5204        enddo
5205        if (tout) then
5206          write(iout,*) 'ERROR = timeout for recv ',tend1-tstart
5207          call flush(iout)
5208          return
5209        endif
5210        man=muster(mpi_source)
5211
5212 !timeout         call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
5213 !timeout     *                 CG_COMM,muster,ierr)
5214 !        print *, ' receiving output from start # ',ind(1)
5215 !t         print *,'receiving ',MPI_WTIME()
5216 !timeout         man=muster(mpi_source)
5217          call mpi_recv(ind,9,mpi_integer,man,idint,&
5218                        CG_COMM,muster,ierr)
5219 !timeout
5220 !  receive final energies and variables
5221          call mpi_recv(eout,1,mpi_double_precision,&
5222                        man,idreal,CG_COMM,muster,ierr)
5223 !         print *,eout
5224 #ifdef CO_BIAS
5225          call mpi_recv(co,1,mpi_double_precision,&
5226                        man,idreal,CG_COMM,muster,ierr)
5227          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
5228 #endif
5229          call mpi_recv(xout,nvar,mpi_double_precision,&
5230                        man,idreal,CG_COMM,muster,ierr)
5231 !         print *,nvar , ierr
5232          if(vdisulf) nss=ind(6)
5233          if(vdisulf.and.nss.ne.0) then
5234           call mpi_recv(ihpb,nss,mpi_integer,&
5235                        man,idint,CG_COMM,muster,ierr)
5236           call mpi_recv(jhpb,nss,mpi_integer,&
5237                        man,idint,CG_COMM,muster,ierr)
5238          endif
5239 !  halt soldier
5240        if(ihalt.eq.1) then
5241 !        print *,'sending halt to ',man
5242         write(iout,*) 'sending halt to ',man
5243         info(1)=0
5244         call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,ierr)
5245        endif
5246       return
5247       end subroutine recv
5248 !-----------------------------------------------------------------------------
5249       subroutine history_append
5250
5251 !      implicit real*8 (a-h,o-z)
5252 !      include 'DIMENSIONS'
5253 !      include 'COMMON.IOUNITS'
5254
5255 #if defined(AIX) || defined(PGI)
5256        open(icsa_history,file=csa_history,position="append")
5257 #else
5258        open(icsa_history,file=csa_history,access="append")
5259 #endif
5260       return
5261       end subroutine history_append
5262 !-----------------------------------------------------------------------------
5263       subroutine alloc_CSA_arrays
5264
5265       use energy_data, only: ns
5266
5267       mxgr=2*nres
5268
5269       if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3))
5270 ! commom.bank
5271 !      common/varin/
5272 !el      allocate(dihang_in(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5273       allocate(dihang_in(mxang,nres,mxch,5000)) !(mxang,maxres,mxch,mxio)
5274       allocate(nss_in(mxio)) !(mxio)
5275       allocate(iss_in(ns,mxio),jss_in(ns,mxio)) !(maxss,mxio)
5276 !      common/minvar/
5277       allocate(dihang(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5278       allocate(rmsn(mxio),pncn(mxio)) !(mxio)
5279       allocate(etot(mxio)) !(mxio)
5280       allocate(nss_out(mxio)) !(mxio)
5281       allocate(iss_out(ns,mxio),jss_out(ns,mxio)) !(maxss,mxio)
5282 !      common/bank/
5283       allocate(rvar(mxang,nres,mxch,mxio),bvar(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5284       allocate(bene(mxio),rene(mxio),brmsn(mxio),rrmsn(mxio))
5285       allocate(bpncn(mxio),rpncn(mxio)) !(mxio)
5286       allocate(ibank(mxio),is(mxio),jbank(mxio)) !(mxio)
5287       allocate(dij(mxio,mxio)) !(mxio,mxio)
5288 !      common/bank_disulfid/
5289       allocate(bvar_nss(mxio),bvar_ns(mxio)) !(mxio)
5290       allocate(bvar_s(ns,mxio)) !(maxss,mxio)
5291       allocate(bvar_ss(2,ns,mxio)) !(2,maxss,mxio)
5292 !      common/mvstat/
5293       allocate(movenx(mxio),movernx(mxio)) !(mxio)
5294       allocate(nstatnx(0:mxmv,3),nstatnx_tot(0:mxmv,3)) !(0:mxmv,3)
5295       allocate(indb(mxio,9)) !(mxio,9)
5296       allocate(parent(3,mxio)) !(3,mxio)
5297 !      common/send2/
5298       allocate(isend2(mxio)) !(mxio)
5299       allocate(iff_in(nres,mxio2)) !(maxres,mxio2)
5300       allocate(dihang_in2(mxang,nres,mxch,mxio2)) !(mxang,maxres,mxch,mxio2)
5301       allocate(idata(5,mxio)) !(5,mxio)
5302 ! common.csa
5303 !      common/alphaa/
5304       allocate(ngroup(mxgr)) !(mxgr)
5305       allocate(igroup(3,mxang,mxgr)) !(3,mxang,mxgr)
5306 ! common.distfit
5307 !      COMMON /frag/
5308       allocate(bvar_frag(mxio,6)) !(mxio,6)
5309       allocate(hvar_frag(mxio,3),lvar_frag(mxio,3),svar_frag(mxio,3)) !(mxio,3)
5310       allocate(avar_frag(mxio,5)) !(mxio,5)
5311 ! commom.hairpin
5312 !      common /spinka/
5313       allocate(nharp_seed(nseed),nharp_use(nseed))      !(max_seed)
5314       allocate(iharp_seed(4,nres/3,nseed))      !(4,maxres/3,max_seed)
5315       allocate(iharp_use(0:4,nres/3,nseed))     !(0:4,maxres/3,max_seed)
5316
5317       return
5318       end subroutine alloc_CSA_arrays
5319 !-----------------------------------------------------------------------------
5320 !-----------------------------------------------------------------------------
5321       end module csa