Merge branch 'devel' into AFM
[unres.git] / source / unres / src_MD-M-newcorr / newconf.f
diff --git a/source/unres/src_MD-M-newcorr/newconf.f b/source/unres/src_MD-M-newcorr/newconf.f
new file mode 100644 (file)
index 0000000..5f93b95
--- /dev/null
@@ -0,0 +1,2454 @@
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine make_var(n,idum,iter_csa)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.HAIRPIN'
+      include 'COMMON.VAR'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      logical nicht_getan,nicht_getan1,fail,lfound
+      integer nharp,iharp(4,maxres/3),nconf_harp
+      integer iisucc(mxio)
+      logical ifused(mxio)
+      integer nhx_seed(max_seed),ihx_seed(4,maxres/3,max_seed)
+      integer nhx_use(max_seed),ihx_use(0:4,maxres/3,max_seed)
+      integer nlx_seed(max_seed),ilx_seed(2,maxres/3,max_seed),
+     &         nlx_use(max_seed),ilx_use(maxres/3,max_seed)
+      real ran1,ran2
+
+      write (iout,*) 'make_var : nseed=',nseed,'ntry=',n
+      index=0
+
+c-----------------------------------------
+      if (n7.gt.0.or.n8.gt.0.or.n9.gt.0.or.n14.gt.0.or.n15.gt.0
+     &    .or.n16.gt.0.or.n17.gt.0.or.n18.gt.0) 
+     &           call select_frag(n7frag,n8frag,n14frag,
+     &           n15frag,nbefrag,iter_csa) 
+
+c---------------------------------------------------
+c N18 - random perturbation of one phi(=gamma) angle in a loop
+c
+      IF (n18.gt.0) THEN
+      nlx_tot=0
+      do iters=1,nseed
+        i1=is(iters)
+        nlx_seed(iters)=0
+        do i2=1,n14frag
+          if (lvar_frag(i2,1).eq.i1) then
+            nlx_seed(iters)=nlx_seed(iters)+5
+            ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
+            ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
+            ilx_use(nlx_seed(iters),iters)=5
+          endif
+        enddo
+        nlx_use(iters)=nlx_seed(iters)
+        nlx_tot=nlx_tot+nlx_seed(iters)
+      enddo
+
+      if (nlx_tot .ge. n18*nseed) then
+        ntot_gen=n18*nseed
+      else
+        ntot_gen=(nlx_tot/nseed)*nseed
+      endif
+
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+        if (nlx_use(iters).gt.0) then
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nlx_seed(iters))
+            if (ilx_use(iih,iters).gt.0) then
+              nicht_getan=.false.
+              ilx_use(iih,iters)=ilx_use(iih,iters)-1
+              nlx_use(iters)=nlx_use(iters)-1
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=18
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+           
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+
+          jr=iran_num(ilx_seed(1,iih,iters),ilx_seed(2,iih,iters))
+          d=ran_number(-pi,pi)
+          dihang_in(2,jr-2,1,index)=pinorm(dihang_in(2,jr-2,1,index)+d)
+
+
+          if (ngen.eq.ntot_gen) goto 145
+        endif
+      enddo
+      enddo
+  145 continue
+
+      ENDIF
+
+
+c-----------------------------------------
+c N17 : zip a beta in a seed by forcing one additional p-p contact
+c
+      IF (n17.gt.0) THEN
+      nhx_tot=0
+      do iters=1,nseed
+        i1=is(iters)
+        nhx_seed(iters)=0
+        nhx_use(iters)=0
+        do i2=1,nbefrag
+          if (avar_frag(i2,1).eq.i1) then
+            nhx_seed(iters)=nhx_seed(iters)+1
+            ihx_use(2,nhx_seed(iters),iters)=1
+            if (avar_frag(i2,5)-avar_frag(i2,3).le.3.and.
+     &           avar_frag(i2,2).gt.1.and.avar_frag(i2,4).lt.nres) then
+             ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+             ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
+             ihx_use(0,nhx_seed(iters),iters)=1
+             ihx_use(1,nhx_seed(iters),iters)=0
+             nhx_use(iters)=nhx_use(iters)+1
+            else
+              if (avar_frag(i2,4).gt.avar_frag(i2,5)) then
+               if (avar_frag(i2,2).gt.1.and.
+     &                     avar_frag(i2,4).lt.nres) then
+                ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+                ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
+                ihx_use(0,nhx_seed(iters),iters)=1
+                ihx_use(1,nhx_seed(iters),iters)=0
+                nhx_use(iters)=nhx_use(iters)+1
+               endif
+               if (avar_frag(i2,3).lt.nres.and.
+     &                     avar_frag(i2,5).gt.1) then
+                ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
+                ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)-1
+                ihx_use(0,nhx_seed(iters),iters)=
+     &                    ihx_use(0,nhx_seed(iters),iters)+1
+                ihx_use(2,nhx_seed(iters),iters)=0
+                nhx_use(iters)=nhx_use(iters)+1
+               endif
+              else
+               if (avar_frag(i2,2).gt.1.and.
+     &                     avar_frag(i2,4).gt.1) then
+                ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
+                ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)-1
+                ihx_use(0,nhx_seed(iters),iters)=1
+                ihx_use(1,nhx_seed(iters),iters)=0
+                nhx_use(iters)=nhx_use(iters)+1
+               endif
+               if (avar_frag(i2,3).lt.nres.and.
+     &                     avar_frag(i2,5).lt.nres) then
+                ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
+                ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)+1
+                ihx_use(0,nhx_seed(iters),iters)=
+     &                    ihx_use(0,nhx_seed(iters),iters)+1
+                ihx_use(2,nhx_seed(iters),iters)=0
+                nhx_use(iters)=nhx_use(iters)+1
+               endif
+              endif
+            endif
+          endif
+        enddo
+
+        nhx_tot=nhx_tot+nhx_use(iters)
+cd        write (iout,*) "debug N17",iters,nhx_seed(iters),
+cd     &                     nhx_use(iters),nhx_tot
+      enddo
+
+      if (nhx_tot .ge. n17*nseed) then
+        ntot_gen=n17*nseed
+      else if (nhx_tot .ge. nseed) then
+        ntot_gen=(nhx_tot/nseed)*nseed
+      else
+        ntot_gen=nhx_tot
+      endif
+cd      write (iout,*) "debug N17==",ntot_gen,nhx_tot,nseed
+
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+        if (nhx_use(iters).gt.0) then
+cd          write (iout,*) "debug N17",nhx_use(iters),ngen,ntot_gen
+cd          write (iout,*) "debugN17^",
+cd     &                    (ihx_use(0,k,iters),k=1,nhx_use(iters))
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nhx_seed(iters))
+cd            write (iout,*) "debugN17^",iih
+            if (ihx_use(0,iih,iters).gt.0) then
+              iim=iran_num(1,2)
+cd                write (iout,*) "debugN17=",iih,nhx_seed(iters)
+cd                write (iout,*) "debugN17-",iim,'##',
+cd     &                           (ihx_use(k,iih,iters),k=0,2)
+cd                call flush(iout)
+              do while (ihx_use(iim,iih,iters).eq.1)
+                iim=iran_num(1,2)
+cd                write (iout,*) "debugN17-",iim,'##',
+cd     &                           (ihx_use(k,iih,iters),k=0,2)
+cd                call flush(iout)
+              enddo
+              nicht_getan=.false.
+              ihx_use(iim,iih,iters)=1
+              ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
+              nhx_use(iters)=nhx_use(iters)-1
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=17
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+
+          if (iim.eq.1) then
+            idata(1,index)=ihx_seed(1,iih,iters)
+            idata(2,index)=ihx_seed(2,iih,iters)
+          else
+            idata(1,index)=ihx_seed(3,iih,iters)
+            idata(2,index)=ihx_seed(4,iih,iters)
+          endif
+
+          if (ngen.eq.ntot_gen) goto 115
+        endif
+      enddo
+      enddo
+  115 continue
+      write (iout,*) "N17",n17," ngen/nseed",ngen/nseed,
+     &                           ngen,nseed
+
+
+      ENDIF
+c-----------------------------------------
+c N16 : slide non local beta in a seed by +/- 1 or +/- 2
+c
+      IF (n16.gt.0) THEN
+      nhx_tot=0
+      do iters=1,nseed
+        i1=is(iters)
+        nhx_seed(iters)=0
+        do i2=1,n7frag
+          if (bvar_frag(i2,1).eq.i1) then
+            nhx_seed(iters)=nhx_seed(iters)+1
+            ihx_seed(1,nhx_seed(iters),iters)=bvar_frag(i2,3)
+            ihx_seed(2,nhx_seed(iters),iters)=bvar_frag(i2,4)
+            ihx_seed(3,nhx_seed(iters),iters)=bvar_frag(i2,5)
+            ihx_seed(4,nhx_seed(iters),iters)=bvar_frag(i2,6)
+            ihx_use(0,nhx_seed(iters),iters)=4
+            do i3=1,4
+              ihx_use(i3,nhx_seed(iters),iters)=0
+            enddo
+          endif
+        enddo
+        nhx_use(iters)=4*nhx_seed(iters)
+        nhx_tot=nhx_tot+nhx_seed(iters)
+cd        write (iout,*) "debug N16",iters,nhx_seed(iters)
+      enddo
+
+      if (4*nhx_tot .ge. n16*nseed) then
+        ntot_gen=n16*nseed
+      else if (4*nhx_tot .ge. nseed) then
+        ntot_gen=(4*nhx_tot/nseed)*nseed
+      else
+        ntot_gen=4*nhx_tot
+      endif
+      write (iout,*) "debug N16",ntot_gen,4*nhx_tot,nseed
+
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+        if (nhx_use(iters).gt.0) then
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nhx_seed(iters))
+            if (ihx_use(0,iih,iters).gt.0) then
+              iim=iran_num(1,4)
+              do while (ihx_use(iim,iih,iters).eq.1)
+cd                write (iout,*) iim,
+cd     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
+                iim=iran_num(1,4)
+              enddo
+              nicht_getan=.false.
+              ihx_use(iim,iih,iters)=1
+              ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
+              nhx_use(iters)=nhx_use(iters)-1
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=16
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+
+          do i=1,4
+           idata(i,index)=ihx_seed(i,iih,iters)
+          enddo
+          idata(5,index)=iim           
+
+          if (ngen.eq.ntot_gen) goto 116
+        endif
+      enddo
+      enddo
+  116 continue
+      write (iout,*) "N16",n16," ngen/nseed",ngen/nseed,
+     &                           ngen,nseed
+      ENDIF
+c-----------------------------------------
+c N15 : copy two 2nd structure elements from 1 or 2 conf. in bank to a seed
+c
+      IF (n15.gt.0) THEN
+
+       do iters=1,nseed
+         iseed=is(iters)
+         do i=1,mxio
+           ifused(i)=.false.
+         enddo
+
+         do idummy=1,n15
+           iter=0
+   84      continue
+
+           iran=0
+           iif=iran_num(1,n15frag)
+           do while( (ifused(iif) .or. svar_frag(iif,1).eq.iseed) .and.
+     &                      iran.le.mxio )
+            iif=iran_num(1,n15frag)
+            iran=iran+1
+           enddo
+           if(iran.ge.mxio) goto 811
+
+           iran=0
+           iig=iran_num(1,n15frag)
+           do while( (ifused(iig) .or. svar_frag(iig,1).eq.iseed  .or. 
+     &             .not.(svar_frag(iif,3).lt.svar_frag(iig,2).or.
+     &                   svar_frag(iig,3).lt.svar_frag(iif,2)) ) .and.
+     &                      iran.le.mxio )
+            iig=iran_num(1,n15frag)
+            iran=iran+1
+           enddo
+           if(iran.ge.mxio) goto 811
+
+           index=index+1
+           movenx(index)=15
+           parent(1,index)=iseed
+           parent(2,index)=svar_frag(iif,1)
+           parent(3,index)=svar_frag(iig,1)
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+           ifused(iif)=.true.
+           ifused(iig)=.true.
+           call newconf_copy(idum,dihang_in(1,1,1,index),
+     &             svar_frag(iif,1),svar_frag(iif,2),svar_frag(iif,3))
+
+           do j=svar_frag(iig,2),svar_frag(iig,3)
+             do i=1,4
+              dihang_in(i,j,1,index)=bvar(i,j,1,svar_frag(iig,1))
+             enddo
+           enddo
+
+
+           if(iter.lt.10) then
+            call check_old(icheck,index)
+            if(icheck.eq.1) then 
+              index=index-1
+              ifused(iif)=.false. 
+              goto 84
+            endif
+           endif
+
+  811     continue
+         enddo
+       enddo
+      ENDIF
+
+c-----------------------------------------
+c N14 local_move (Maurizio) for loops in a seed
+c
+      IF (n14.gt.0) THEN
+      nlx_tot=0
+      do iters=1,nseed
+        i1=is(iters)
+        nlx_seed(iters)=0
+        do i2=1,n14frag
+          if (lvar_frag(i2,1).eq.i1) then
+            nlx_seed(iters)=nlx_seed(iters)+3
+            ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
+            ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
+            ilx_use(nlx_seed(iters),iters)=3
+          endif
+        enddo
+        nlx_use(iters)=nlx_seed(iters)
+        nlx_tot=nlx_tot+nlx_seed(iters)
+cd        write (iout,*) "debug N14",iters,nlx_seed(iters)
+      enddo
+
+      if (nlx_tot .ge. n14*nseed) then
+        ntot_gen=n14*nseed
+      else
+        ntot_gen=(nlx_tot/nseed)*nseed
+      endif
+cd      write (iout,*) "debug N14",ntot_gen,n14frag,nseed
+
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+        if (nlx_use(iters).gt.0) then
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nlx_seed(iters))
+            if (ilx_use(iih,iters).gt.0) then
+              nicht_getan=.false.
+              ilx_use(iih,iters)=ilx_use(iih,iters)-1
+              nlx_use(iters)=nlx_use(iters)-1
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=14
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+          idata(1,index)=ilx_seed(1,iih,iters)
+          idata(2,index)=ilx_seed(2,iih,iters)
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+           
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+
+          if (ngen.eq.ntot_gen) goto 131
+        endif
+      enddo
+      enddo
+  131 continue
+cd      write (iout,*) "N14",n14," ngen/nseed",ngen/nseed,
+cd     &                           ngen,nseed
+
+      ENDIF
+c-----------------------------------------
+c N9 : shift a helix in a seed
+c
+      IF (n9.gt.0) THEN
+      nhx_tot=0
+      do iters=1,nseed
+        i1=is(iters)
+        nhx_seed(iters)=0
+        do i2=1,n8frag
+          if (hvar_frag(i2,1).eq.i1) then
+            nhx_seed(iters)=nhx_seed(iters)+1
+            ihx_seed(1,nhx_seed(iters),iters)=hvar_frag(i2,2)
+            ihx_seed(2,nhx_seed(iters),iters)=hvar_frag(i2,3)
+            ihx_use(0,nhx_seed(iters),iters)=4
+            do i3=1,4
+              ihx_use(i3,nhx_seed(iters),iters)=0
+            enddo
+          endif
+        enddo
+        nhx_use(iters)=4*nhx_seed(iters)
+        nhx_tot=nhx_tot+nhx_seed(iters)
+cd        write (iout,*) "debug N9",iters,nhx_seed(iters)
+      enddo
+
+      if (4*nhx_tot .ge. n9*nseed) then
+        ntot_gen=n9*nseed
+      else
+        ntot_gen=(4*nhx_tot/nseed)*nseed
+      endif
+cd      write (iout,*) "debug N9",ntot_gen,n8frag,nseed
+
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+        if (nhx_use(iters).gt.0) then
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nhx_seed(iters))
+            if (ihx_use(0,iih,iters).gt.0) then
+              iim=iran_num(1,4)
+              do while (ihx_use(iim,iih,iters).eq.1)
+cd                write (iout,*) iim,
+cd     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
+                iim=iran_num(1,4)
+              enddo
+              nicht_getan=.false.
+              ihx_use(iim,iih,iters)=1
+              ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
+              nhx_use(iters)=nhx_use(iters)-1
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=9
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+
+          jstart=max(nnt,ihx_seed(1,iih,iters)+1)
+          jend=min(nct,ihx_seed(2,iih,iters))
+cd          write (iout,*) "debug N9",iters,iih,jstart,jend
+          if (iim.eq.1) then
+            ishift=-2
+          else if (iim.eq.2) then
+            ishift=-1
+          else if (iim.eq.3) then
+            ishift=1
+          else if (iim.eq.4) then
+            ishift=2
+          else
+            write (iout,*) 'CHUJ NASTAPIL: iim=',iim
+            call mpi_abort(mpi_comm_world,ierror,ierrcode)
+          endif
+          do j=jstart,jend
+            if (itype(j).eq.10) then
+              iang=2
+            else
+              iang=4
+            endif
+            do i=1,iang
+              if (j+ishift.ge.nnt.and.j+ishift.le.nct)
+     &          dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
+            enddo
+          enddo
+          if (ishift.gt.0) then
+           do j=0,ishift-1
+            if (itype(jend+j).eq.10) then
+              iang=2
+            else
+              iang=4
+            endif
+            do i=1,iang
+              if (jend+j.ge.nnt.and.jend+j.le.nct)
+     &          dihang_in(i,jstart+j,1,index)=bvar(i,jend+j,1,iseed)
+            enddo
+           enddo
+          else
+           do j=0,-ishift-1
+            if (itype(jstart+j).eq.10) then
+              iang=2
+            else
+              iang=4
+            endif
+            do i=1,iang
+              if (jend+j.ge.nnt.and.jend+j.le.nct)
+     &          dihang_in(i,jend+j,1,index)=bvar(i,jstart+j,1,iseed)
+            enddo
+           enddo
+          endif
+          if (ngen.eq.ntot_gen) goto 133
+        endif
+      enddo
+      enddo
+  133 continue
+cd      write (iout,*) "N9",n9," ngen/nseed",ngen/nseed,
+cd     &                           ngen,nseed
+
+      ENDIF
+c-----------------------------------------
+c N8 : copy a helix from bank to seed
+c
+      if (n8.gt.0) then 
+       if (n8frag.lt.n8) then 
+          write (iout,*) "N8: only ",n8frag,'helices'
+          n8c=n8frag
+       else
+          n8c=n8
+       endif
+
+       do iters=1,nseed
+         iseed=is(iters)
+         do i=1,mxio
+           ifused(i)=.false.
+         enddo
+
+
+         do idummy=1,n8c
+           iter=0
+   94      continue
+           iran=0
+           iif=iran_num(1,n8frag)
+           do while( (ifused(iif) .or. hvar_frag(iif,1).eq.iseed) .and.
+     &                      iran.le.mxio )
+            iif=iran_num(1,n8frag)
+            iran=iran+1
+           enddo
+            
+           if(iran.ge.mxio) goto 911
+
+           index=index+1
+           movenx(index)=8
+           parent(1,index)=iseed
+           parent(2,index)=hvar_frag(iif,1)
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+           ifused(iif)=.true.
+           if (hvar_frag(iif,3)-hvar_frag(iif,2).le.6) then
+             call newconf_copy(idum,dihang_in(1,1,1,index),
+     &             hvar_frag(iif,1),hvar_frag(iif,2),hvar_frag(iif,3))
+           else
+              ih_start=iran_num(hvar_frag(iif,2),hvar_frag(iif,3)-6)
+              ih_end=iran_num(ih_start,hvar_frag(iif,3))
+             call newconf_copy(idum,dihang_in(1,1,1,index),
+     &             hvar_frag(iif,1),ih_start,ih_end)
+           endif
+           iter=iter+1
+           if(iter.lt.10) then
+            call check_old(icheck,index)
+            if(icheck.eq.1) then 
+              index=index-1
+              ifused(iif)=.false. 
+              goto 94
+            endif
+           endif
+
+
+  911     continue
+
+         enddo
+       enddo
+
+      endif
+
+c-----------------------------------------
+c N7 : copy nonlocal beta fragment from bank to seed
+c
+      if (n7.gt.0) then 
+       if (n7frag.lt.n7) then 
+          write (iout,*) "N7: only ",n7frag,'nonlocal fragments'
+          n7c=n7frag
+       else
+          n7c=n7
+       endif
+
+       do i=1,maxres
+         do j=1,mxio2
+            iff_in(i,j)=0
+         enddo
+       enddo
+       index2=0
+       do i=1,mxio
+         isend2(i)=0
+       enddo
+
+       do iters=1,nseed
+         iseed=is(iters)
+         do i=1,mxio
+           ifused(i)=.false.
+         enddo
+
+         do idummy=1,n7c
+           iran=0
+           iif=iran_num(1,n7frag)
+           do while( (ifused(iif) .or. bvar_frag(iif,1).eq.iseed) .and.
+     &                      iran.le.mxio )
+            iif=iran_num(1,n7frag)
+            iran=iran+1
+           enddo
+            
+cd       write (*,'(3i5,l,4i5)'),iters,idummy,iif,ifused(iif),
+cd     &                    bvar_frag(iif,1),iseed,iran,index2
+         
+           if(iran.ge.mxio) goto 999
+           if(index2.ge.mxio2) goto 999 
+
+           index=index+1
+           movenx(index)=7
+           parent(1,index)=iseed
+           parent(2,index)=bvar_frag(iif,1)
+           index2=index2+1
+           isend2(index)=index2
+           ifused(iif)=.true.
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+           do k=1,numch
+            do j=2,nres-1
+             do i=1,4
+              dihang_in2(i,j,k,index2)=bvar(i,j,k,bvar_frag(iif,1))
+             enddo
+            enddo
+           enddo
+
+           if (bvar_frag(iif,2).eq.4) then                                                                  
+             do i=bvar_frag(iif,3),bvar_frag(iif,4)
+               iff_in(i,index2)=1                                                              
+             enddo                                                                   
+             if (bvar_frag(iif,5).lt.bvar_frag(iif,6)) then
+cd               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
+cd     &                 bvar_frag(iif,5),bvar_frag(iif,6)
+               do i=bvar_frag(iif,5),bvar_frag(iif,6)
+                 iff_in(i,index2)=1                                                              
+               enddo                                                                   
+             else
+cd               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
+cd     &                 bvar_frag(iif,6),bvar_frag(iif,5)
+               do i=bvar_frag(iif,6),bvar_frag(iif,5)
+                 iff_in(i,index2)=1                                                              
+               enddo                                                                   
+             endif
+           endif
+
+           do k=1,numch
+            do j=2,nres-1
+             do i=1,4
+              dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+             enddo
+            enddo
+           enddo
+
+
+  999     continue
+
+         enddo
+       enddo
+
+      endif
+c-----------------------------------------------
+c N6 : copy random continues fragment from bank to seed
+c
+      do iters=1,nseed
+       iseed=is(iters)
+       do idummy=1,n6
+        isize=(is2-is1+1)*ran1(idum)+is1
+        index=index+1
+        movenx(index)=6
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+        iter=0
+  104   continue
+        if(icycle.le.0) then
+         i1=nconf*  ran1(idum)+1
+         i1=nbank-nconf+i1
+        else
+         i1=nbank*  ran1(idum)+1
+        endif
+        if(i1.eq.iseed) goto 104
+        iter=iter+1
+        call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+        parent(1,index)=iseed
+        parent(2,index)=i1
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 104
+        endif
+       enddo
+      enddo
+c-----------------------------------------
+      if (n3.gt.0.or.n4.gt.0) call gen_hairpin
+      nconf_harp=0
+      do iters=1,nseed
+        if (nharp_seed(iters).gt.0) nconf_harp=nconf_harp+1
+      enddo      
+c-----------------------------------------
+c N3 : copy hairpin from bank to seed
+c
+      do iters=1,nseed
+       iseed=is(iters)
+       nsucc=0
+       nacc=0
+       do idummy=1,n3
+        index=index+1
+        iter=0
+  124   continue
+        if(icycle.le.0) then
+         i1=nconf*  ran1(idum)+1
+         i1=nbank-nconf+i1
+        else
+         i1=nbank*  ran1(idum)+1
+        endif
+        if(i1.eq.iseed) goto 124
+        do k=1,nsucc
+          if (i1.eq.iisucc(k).and.nsucc.lt.nconf_harp-1) goto 124
+        enddo
+        nsucc=nsucc+1
+        iisucc(nsucc)=i1
+        iter=iter+1
+        call newconf_residue_hairpin(idum,dihang_in(1,1,1,index),
+     &     i1,fail)
+        if (fail) then
+          if (icycle.le.0 .and. nsucc.eq.nconf .or. 
+     &        icycle.gt.0 .and. nsucc.eq.nbank)  then
+            index=index-1
+            goto 125
+          else
+            goto 124
+          endif
+        endif
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 124
+        endif
+        movenx(index)=3
+        parent(1,index)=iseed
+        parent(2,index)=i1
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+        nacc=nacc+1
+       enddo
+c if not enough hairpins, supplement with windows
+  125  continue
+cdd       if (n3.ne.0) write (iout,*) "N3",n3," nsucc",nsucc," nacc",nacc 
+       do idummy=nacc+1,n3
+        isize=(is2-is1+1)*ran1(idum)+is1
+        index=index+1
+        movenx(index)=6
+        parent(1,index)=iseed
+        parent(2,index)=i1
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+        iter=0
+  114   continue
+        if(icycle.le.0) then
+         i1=nconf*  ran1(idum)+1
+         i1=nbank-nconf+i1
+        else
+         i1=nbank*  ran1(idum)+1
+        endif
+        if(i1.eq.iseed) goto 114
+        iter=iter+1
+        call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 114
+        endif
+       enddo
+      enddo
+c-----------------------------------------
+c N4 : shift a turn in hairpin in seed
+c
+      IF (N4.GT.0) THEN
+      if (4*nharp_tot .ge. n4*nseed) then
+        ntot_gen=n4*nseed
+      else
+        ntot_gen=(4*nharp_tot/nseed)*nseed
+      endif
+      ngen=0
+      do while (ngen.lt.ntot_gen)
+      do iters=1,nseed
+        iseed=is(iters)
+c        write (iout,*) 'iters',iters,' iseed',iseed,' nharp_seed',
+c     &     nharp_seed(iters),' nharp_use',nharp_use(iters),
+c     &     ' ntot_gen',ntot_gen
+c        write (iout,*) 'iharp_use(0)',
+c     &     (iharp_use(0,k,iters),k=1,nharp_seed(iters))
+        if (nharp_use(iters).gt.0) then
+          nicht_getan=.true.
+          do while (nicht_getan)
+            iih=iran_num(1,nharp_seed(iters))
+c            write (iout,*) 'iih',iih,' iharp_use',
+c     &            (iharp_use(k,iih,iters),k=1,4)
+            if (iharp_use(0,iih,iters).gt.0) then
+              nicht_getan1=.true.
+              do while (nicht_getan1)
+                iim=iran_num(1,4)
+                nicht_getan1=iharp_use(iim,iih,iters).eq.1
+              enddo
+              nicht_getan=.false.
+              iharp_use(iim,iih,iters)=1
+              iharp_use(0,iih,iters)=iharp_use(0,iih,iters)-1 
+              nharp_use(iters)=nharp_use(iters)-1
+cdd              write (iout,'(a16,i3,a5,i2,a10,2i4)') 
+cdd     &           'N4 selected hairpin',iih,' move',iim,' iharp_seed',
+cdd     &            iharp_seed(1,iih,iters),iharp_seed(2,iih,iters)
+            endif
+          enddo
+          ngen=ngen+1
+          index=index+1
+          movenx(index)=4
+          parent(1,index)=iseed
+          parent(2,index)=0
+
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+          do k=1,numch
+            do j=2,nres-1
+              do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+              enddo
+            enddo
+          enddo
+          jstart=iharp_seed(1,iih,iters)+1
+          jend=iharp_seed(2,iih,iters)
+          if (iim.eq.1) then
+            ishift=-2
+          else if (iim.eq.2) then
+            ishift=-1
+          else if (iim.eq.3) then
+            ishift=1
+          else if (iim.eq.4) then
+            ishift=2
+          else
+            write (iout,*) 'CHUJ NASTAPIL: iim=',iim
+            call mpi_abort(mpi_comm_world,ierror,ierrcode)
+          endif
+c          write (iout,*) 'jstart',jstart,' jend',jend,' ishift',ishift
+c          write (iout,*) 'Before turn shift'
+c          do j=2,nres-1
+c            theta(j+1)=dihang_in(1,j,1,index)
+c            phi(j+2)=dihang_in(2,j,1,index)
+c            alph(j)=dihang_in(3,j,1,index)
+c            omeg(j)=dihang_in(4,j,1,index)
+c          enddo
+c          call intout
+          do j=jstart,jend
+            if (itype(j).eq.10) then
+              iang=2
+            else
+              iang=4
+            endif
+            do i=1,iang
+              if (j+ishift.ge.nnt.and.j+ishift.le.nct)
+     &          dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
+            enddo
+          enddo
+c          write (iout,*) 'After turn shift'
+c          do j=2,nres-1
+c            theta(j+1)=dihang_in(1,j,1,index)
+c            phi(j+2)=dihang_in(2,j,1,index)
+c            alph(j)=dihang_in(3,j,1,index)
+c            omeg(j)=dihang_in(4,j,1,index)
+c          enddo
+c          call intout
+          if (ngen.eq.ntot_gen) goto 135
+        endif
+      enddo
+      enddo
+c if not enough hairpins, supplement with windows
+c      write (iout,*) 'end of enddo'
+  135 continue
+cdd      write (iout,*) "N4",n4," ngen/nseed",ngen/nseed,
+cdd     &                           ngen,nseed
+      do iters=1,nseed
+        iseed=is(iters)
+        do idummy=ngen/nseed+1,n4
+          isize=(is2-is1+1)*ran1(idum)+is1
+          index=index+1
+          movenx(index)=6
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+
+          iter=0
+  134     continue
+          if(icycle.le.0) then
+          i1=nconf*  ran1(idum)+1
+          i1=nbank-nconf+i1
+          else
+            i1=nbank*  ran1(idum)+1
+          endif
+          if(i1.eq.iseed) goto 134
+            iter=iter+1
+            call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+            parent(1,index)=iseed
+            parent(2,index)=i1
+            if(iter.lt.10) then
+              call check_old(icheck,index)
+            if(icheck.eq.1) goto 134
+          endif
+        enddo
+      enddo
+      ENDIF
+c-----------------------------------------
+c N5 : copy one residue from bank to seed (normally switched off - use N1)
+c
+      do iters=1,nseed
+       iseed=is(iters)
+       isize=1
+       do i=1,n5
+        index=index+1
+        movenx(index)=5
+
+           if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+           endif
+
+
+        iter=0
+  105   continue
+        if(icycle.le.0) then
+         i1=nconf*  ran1(idum)+1
+         i1=nbank-nconf+i1
+        else
+         i1=nbank*  ran1(idum)+1
+        endif
+        if(i1.eq.iseed) goto 105
+        iter=iter+1
+        call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
+        parent(1,index)=iseed
+        parent(2,index)=i1
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 105
+        endif
+       enddo
+      enddo
+c-----------------------------------------
+c N2 : copy backbone of one residue from bank or first bank to seed
+c (normally switched off - use N1)
+c
+      do iters=1,nseed
+       iseed=is(iters)
+       do i=n2,1,-1
+        if(icycle.le.0.and.iuse.gt.nconf-irr) then
+         iseed=ran1(idum)*nconf+1
+         iseed=nbank-nconf+iseed
+        endif
+        index=index+1
+        movenx(index)=2
+
+        if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+        endif
+
+        iter=0
+  102   i1=  ran1(idum)*nbank+1
+        if(i1.eq.iseed) goto 102
+        iter=iter+1
+        if(icycle.le.0.and.iuse.gt.nconf-irr) then
+         nran=mod(i-1,nran0)+3
+         call newconf1arr(idum,dihang_in(1,1,1,index),nran,i1)
+         parent(1,index)=-iseed
+         parent(2,index)=-i1
+        else if(icycle.le.0.and.iters.le.iuse) then
+         nran=mod(i-1,nran0)+1
+         call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
+         parent(1,index)=iseed
+         parent(2,index)=-i1
+        else
+         nran=mod(i-1,nran1)+1
+         if(ran1(idum).lt.0.5) then
+          call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
+          parent(1,index)=iseed
+          parent(2,index)=-i1
+         else
+          call newconf1abb(idum,dihang_in(1,1,1,index),nran,i1)
+          parent(1,index)=iseed
+          parent(2,index)=i1
+         endif
+        endif
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 102
+        endif
+       enddo
+      enddo
+c-----------------------------------------
+c N1 : copy backbone or sidechain of one residue from bank or 
+c      first bank to seed
+c
+      do iters=1,nseed
+       iseed=is(iters)
+       do i=n1,1,-1
+        if(icycle.le.0.and.iuse.gt.nconf-irr) then
+         iseed=ran1(idum)*nconf+1
+         iseed=nbank-nconf+iseed
+        endif
+        index=index+1
+        movenx(index)=1
+
+        if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+        endif
+
+        iter=0
+  101   i1=  ran1(idum)*nbank+1
+
+        if(i1.eq.iseed) goto 101
+        iter=iter+1
+        if(icycle.le.0.and.iuse.gt.nconf-irr) then
+         nran=mod(i-1,nran0)+3
+         call newconf1rr(idum,dihang_in(1,1,1,index),nran,i1)
+         parent(1,index)=-iseed
+         parent(2,index)=-i1
+        else if(icycle.le.0.and.iters.le.iuse) then
+         nran=mod(i-1,nran0)+1
+         call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
+         parent(1,index)=iseed
+         parent(2,index)=-i1
+        else
+         nran=mod(i-1,nran1)+1
+         if(ran1(idum).lt.0.5) then
+          call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
+          parent(1,index)=iseed
+          parent(2,index)=-i1
+         else
+          call newconf1bb(idum,dihang_in(1,1,1,index),nran,i1)
+          parent(1,index)=iseed
+          parent(2,index)=i1
+         endif
+        endif
+        if(iter.lt.10) then
+         call check_old(icheck,index)
+         if(icheck.eq.1) goto 101
+        endif
+       enddo
+      enddo
+c-----------------------------------------
+c N0 just all seeds
+c
+      IF (n0.gt.0) THEN
+      do iters=1,nseed
+       iseed=is(iters)
+       index=index+1
+       movenx(index)=0
+       parent(1,index)=iseed
+       parent(2,index)=0
+
+       if (vdisulf) then
+             nss_in(index)=bvar_nss(iseed)
+             do ij=1,nss_in(index)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+       endif
+
+       do k=1,numch
+        do j=2,nres-1
+         do i=1,4
+          dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+         enddo
+        enddo
+       enddo
+      enddo
+      ENDIF
+c-----------------------------------------
+      if (vdisulf) then
+      do iters=1,nseed
+       iseed=is(iters)
+
+        do k=1,numch
+          do j=2,nres-1
+            theta(j+1)=bvar(1,j,k,iseed)
+            phi(j+2)=bvar(2,j,k,iseed)
+            alph(j)=bvar(3,j,k,iseed)
+            omeg(j)=bvar(4,j,k,iseed)
+          enddo
+        enddo
+        call chainbuild
+
+cd       write(iout,*) 'makevar DYNSS',iseed,'#',bvar_ns(iseed),
+cd     &     (bvar_s(k,iseed),k=1,bvar_ns(iseed)),
+cd     &     bvar_nss(iseed),
+cd     &     (bvar_ss(1,k,iseed)-nres,'-',
+cd     &      bvar_ss(2,k,iseed)-nres,k=1,bvar_nss(iseed))
+
+       do i1=1,bvar_ns(iseed)
+c
+c N10 fussion of free halfcysteines in seed 
+c      first select CYS with distance < 7A
+c
+         do j1=i1+1,bvar_ns(iseed)
+           if (dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres)
+     &         .lt.7.0.and.
+     &         iabs(bvar_s(i1,iseed)-bvar_s(j1,iseed)).gt.3) then
+
+             index=index+1
+             movenx(index)=10
+             parent(1,index)=iseed
+             parent(2,index)=0
+             do ij=1,bvar_nss(iseed)
+               iss_in(ij,index)=bvar_ss(1,ij,iseed)
+               jss_in(ij,index)=bvar_ss(2,ij,iseed)
+             enddo
+             ij=bvar_nss(iseed)+1
+             nss_in(index)=ij
+             iss_in(ij,index)=bvar_s(i1,iseed)+nres
+             jss_in(ij,index)=bvar_s(j1,iseed)+nres
+
+cd             write(iout,*) 'makevar NSS0',index,
+cd     &   dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres),
+cd     &   nss_in(index),iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+
+           endif
+         enddo
+c
+c N11 type I transdisulfidation
+c
+         do j1=1,bvar_nss(iseed)
+           if (dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed))
+     &         .lt.7.0.and.
+     &         iabs(bvar_s(i1,iseed)-(bvar_ss(1,j1,iseed)-nres))
+     &         .gt.3) then
+
+             index=index+1
+             movenx(index)=11
+             parent(1,index)=iseed
+             parent(2,index)=0
+             do ij=1,bvar_nss(iseed)
+               if (ij.ne.j1) then
+                iss_in(ij,index)=bvar_ss(1,ij,iseed)
+                jss_in(ij,index)=bvar_ss(2,ij,iseed)
+               endif
+             enddo
+             nss_in(index)=bvar_nss(iseed)
+             iss_in(j1,index)=bvar_s(i1,iseed)+nres
+             jss_in(j1,index)=bvar_ss(1,j1,iseed)
+             if (iss_in(j1,index).gt.jss_in(j1,index)) then
+               iss_in(j1,index)=bvar_ss(1,j1,iseed)
+               jss_in(j1,index)=bvar_s(i1,iseed)+nres
+             endif
+
+cd             write(iout,*) 'makevar NSS1 #1',index,
+cd     &       bvar_s(i1,iseed),bvar_ss(1,j1,iseed)-nres,
+cd     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)),
+cd     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+cd     &       ij=1,nss_in(index))
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+           endif
+           if (dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed))
+     &         .lt.7.0.and.
+     &         iabs(bvar_s(i1,iseed)-(bvar_ss(2,j1,iseed)-nres))
+     &         .gt.3) then
+
+             index=index+1
+             movenx(index)=11
+             parent(1,index)=iseed
+             parent(2,index)=0
+             do ij=1,bvar_nss(iseed)
+               if (ij.ne.j1) then
+                iss_in(ij,index)=bvar_ss(1,ij,iseed)
+                jss_in(ij,index)=bvar_ss(2,ij,iseed)
+               endif
+             enddo
+             nss_in(index)=bvar_nss(iseed)
+             iss_in(j1,index)=bvar_s(i1,iseed)+nres
+             jss_in(j1,index)=bvar_ss(2,j1,iseed)
+             if (iss_in(j1,index).gt.jss_in(j1,index)) then
+               iss_in(j1,index)=bvar_ss(2,j1,iseed)
+               jss_in(j1,index)=bvar_s(i1,iseed)+nres
+             endif
+
+
+cd             write(iout,*) 'makevar NSS1 #2',index,
+cd     &       bvar_s(i1,iseed),bvar_ss(2,j1,iseed)-nres,
+cd     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)),
+cd     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+cd     &       ij=1,nss_in(index))
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+
+           endif
+         enddo
+       enddo
+
+c
+c N12 type II transdisulfidation
+c
+       do i1=1,bvar_nss(iseed)
+         do j1=i1+1,bvar_nss(iseed)
+            if (dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed))
+     &         .lt.7.0.and.
+     &          dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed))
+     &         .lt.7.0.and.
+     &         iabs(bvar_ss(1,i1,iseed)-bvar_ss(1,j1,iseed))
+     &          .gt.3.and.
+     &         iabs(bvar_ss(2,i1,iseed)-bvar_ss(2,j1,iseed))
+     &          .gt.3) then
+             index=index+1
+             movenx(index)=12
+             parent(1,index)=iseed
+             parent(2,index)=0
+             do ij=1,bvar_nss(iseed)
+               if (ij.ne.i1 .and. ij.ne.j1) then
+                iss_in(ij,index)=bvar_ss(1,ij,iseed)
+                jss_in(ij,index)=bvar_ss(2,ij,iseed)
+               endif
+             enddo
+             nss_in(index)=bvar_nss(iseed)
+             iss_in(i1,index)=bvar_ss(1,i1,iseed)
+             jss_in(i1,index)=bvar_ss(1,j1,iseed)
+             if (iss_in(i1,index).gt.jss_in(i1,index)) then
+               iss_in(i1,index)=bvar_ss(1,j1,iseed)
+               jss_in(i1,index)=bvar_ss(1,i1,iseed)
+             endif
+             iss_in(j1,index)=bvar_ss(2,i1,iseed)
+             jss_in(j1,index)=bvar_ss(2,j1,iseed)
+             if (iss_in(j1,index).gt.jss_in(j1,index)) then
+               iss_in(j1,index)=bvar_ss(2,j1,iseed)
+               jss_in(j1,index)=bvar_ss(2,i1,iseed)
+             endif
+
+
+cd             write(iout,*) 'makevar NSS2 #1',index,
+cd     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres, 
+cd     &       dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)),
+cd     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
+cd     &       dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)),
+cd     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+cd     &       ij=1,nss_in(index))
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+
+           endif
+
+            if (dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed))
+     &         .lt.7.0.and.
+     &          dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed))
+     &         .lt.7.0.and.
+     &         iabs(bvar_ss(1,i1,iseed)-bvar_ss(2,j1,iseed))
+     &          .gt.3.and.
+     &         iabs(bvar_ss(2,i1,iseed)-bvar_ss(1,j1,iseed))
+     &          .gt.3) then
+             index=index+1
+             movenx(index)=12
+             parent(1,index)=iseed
+             parent(2,index)=0
+             do ij=1,bvar_nss(iseed)
+               if (ij.ne.i1 .and. ij.ne.j1) then
+                iss_in(ij,index)=bvar_ss(1,ij,iseed)
+                jss_in(ij,index)=bvar_ss(2,ij,iseed)
+               endif
+             enddo
+             nss_in(index)=bvar_nss(iseed)
+             iss_in(i1,index)=bvar_ss(1,i1,iseed)
+             jss_in(i1,index)=bvar_ss(2,j1,iseed)
+             if (iss_in(i1,index).gt.jss_in(i1,index)) then
+               iss_in(i1,index)=bvar_ss(2,j1,iseed)
+               jss_in(i1,index)=bvar_ss(1,i1,iseed)
+             endif
+             iss_in(j1,index)=bvar_ss(2,i1,iseed)
+             jss_in(j1,index)=bvar_ss(1,j1,iseed)
+             if (iss_in(j1,index).gt.jss_in(j1,index)) then
+               iss_in(j1,index)=bvar_ss(1,j1,iseed)
+               jss_in(j1,index)=bvar_ss(2,i1,iseed)
+             endif
+
+
+cd             write(iout,*) 'makevar NSS2 #2',index,
+cd     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres, 
+cd     &       dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)),
+cd     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
+cd     &       dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)),
+cd     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+cd     &       ij=1,nss_in(index))
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+
+           endif
+
+
+         enddo
+       enddo
+c
+c N13 removal of disulfide bond
+c
+      if (bvar_nss(iseed).gt.0) then
+        i1=bvar_nss(iseed)*ran1(idum)+1
+
+             index=index+1
+             movenx(index)=13
+             parent(1,index)=iseed
+             parent(2,index)=0
+             ij=0
+             do j1=1,bvar_nss(iseed)
+              if (j1.ne.i1) then
+               ij=ij+1
+               iss_in(ij,index)=bvar_ss(1,j1,iseed)
+               jss_in(ij,index)=bvar_ss(2,j1,iseed)
+              endif
+             enddo
+             nss_in(index)=bvar_nss(iseed)-1
+
+cd             write(iout,*) 'NSS3',index,i1,
+cd     &   bvar_ss(1,i1,iseed)-nres,'=',bvar_ss(2,i1,iseed)-nres,'#',
+cd     &   (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
+cd     &   ij=1,nss_in(index))
+
+             do k=1,numch
+              do j=2,nres-1
+               do i=1,4
+                dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
+               enddo
+              enddo
+             enddo
+
+      endif
+
+      enddo
+      endif
+c-----------------------------------------
+
+
+
+      if(index.ne.n) write(iout,*)'make_var : ntry=',index
+
+      n=index
+cd      do ii=1,n
+cd        write (istat,*) "======== ii=",ii," the dihang array"
+cd        do i=1,nres
+cd       write (istat,'(i5,4f15.5)') i,(dihang_in(k,i,1,ii)*rad2deg,k=1,4)
+cd        enddo
+cd      enddo
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine check_old(icheck,n)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      data ctdif /10./
+      data ctdiff /60./
+
+      i1=n
+      do i2=1,n-1
+       diff=0.d0
+       do m=1,numch
+        do j=2,nres-1
+         do i=1,4
+          dif=rad2deg*dabs(dihang_in(i,j,m,i1)-dihang_in(i,j,m,i2))
+          if(dif.gt.180.0) dif=360.0-dif
+          if(dif.gt.ctdif) goto 100
+          diff=diff+dif
+          if(diff.gt.ctdiff) goto 100
+         enddo
+        enddo
+       enddo
+       icheck=1
+       return
+  100  continue
+      enddo
+
+      icheck=0
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1rr(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=rvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=ntotgr
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=rvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1br(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=ntotgr
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+       if(i2ndstr.gt.0) then
+        rtmp=ran1(idum)
+        if(rtmp.le.rdih_bias) then
+         i=0
+         do j=1,ndih_nconstr
+          if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
+         enddo
+         if(i.eq.0) then
+          juhc=0
+4321      juhc=juhc+1
+          iran=  ran1(idum)*number+1
+          i=0
+          do j=1,ndih_nconstr
+           if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
+          enddo
+          if(i.eq.0.or.juhc.lt.1000)goto 4321
+          if(juhc.eq.1000) then
+           print *, 'move 6 : failed to find unconstrained group'
+           write(iout,*) 'move 6 : failed to find unconstrained group'
+          endif
+         endif
+        endif
+       endif
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=rvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1bb(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=ntotgr
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=bvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1arr(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=rvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=nres-2
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=rvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1abr(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=nres-2
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+      if(i2ndstr.gt.0) then
+       rtmp=ran1(idum)
+       if(rtmp.le.rdih_bias) then
+        iran=ran1(idum)*ndih_nconstr+1
+        iran=idih_nconstr(iran)
+       endif
+      endif
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=rvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf1abb(idum,vvar,nran,i1)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+      do index=1,nran
+       iold(index) = 0
+      enddo
+
+      number=nres-2
+
+      iter=0
+      do index=1,nran
+   10  iran=  ran1(idum)*number+1
+      if(i2ndstr.gt.0) then
+       rtmp=ran1(idum)
+       if(rtmp.le.rdih_bias) then
+        iran=ran1(idum)*ndih_nconstr+1
+        iran=idih_nconstr(iran)
+       endif
+      endif
+       if(iter.gt.number) return
+       iter=iter+1
+       if(iter.eq.1) goto 11
+       do ind=1,index-1
+        if(iran.eq.iold(ind)) goto 10
+       enddo
+   11  continue
+
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+        if(dif.gt.180.) dif=360.-dif
+        if(dif.gt.ctdif) goto 20
+       enddo
+       if(iter.gt.number) goto 20
+       goto 10
+   20  continue
+       do ind=1,ngroup(iran)
+        i=igroup(1,ind,iran)
+        j=igroup(2,ind,iran)
+        k=igroup(3,ind,iran)
+        vvar(i,j,k)=bvar(i,j,k,i1)
+       enddo
+       iold(index)=iran
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf_residue(idum,vvar,i1,isize)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      if (iseed.gt.mxio .or. iseed.lt.1) then
+        write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+        call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      endif
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+
+      k=1
+      number=nres+isize-2
+      iter=1
+   10 iran=  ran1(idum)*number+1
+      if(i2ndstr.gt.0) then
+       rtmp=ran1(idum)
+       if(rtmp.le.rdih_bias) then
+        iran=ran1(idum)*ndih_nconstr+1
+        iran=idih_nconstr(iran)
+       endif
+      endif
+      istart=iran-isize+1
+      iend=iran
+      if(istart.lt.2) istart=2
+      if(iend.gt.nres-1) iend=nres-1
+
+      if(iter.eq.1) goto 11
+      do ind=1,iter-1
+       if(iran.eq.iold(ind)) goto 10
+      enddo
+   11 continue
+
+      do j=istart,iend
+       do i=1,4
+         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+         if(dif.gt.180.) dif=360.-dif
+         if(dif.gt.ctdif) goto 20
+       enddo
+      enddo
+      iold(iter)=iran
+      iter=iter+1
+      if(iter.gt.number) goto 20
+      goto 10
+
+   20 continue
+      do j=istart,iend
+       do i=1,4
+        vvar(i,j,k)=bvar(i,j,k,i1)
+       enddo
+      enddo
+
+      return
+      end
+
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf_copy(idum,vvar,i1,istart,iend)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      ctdif=10.
+
+      if (iseed.gt.mxio .or. iseed.lt.1) then
+        write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+        call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      endif
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+
+
+      do j=istart,iend
+       do i=1,4
+        vvar(i,j,1)=bvar(i,j,1,i1)
+       enddo
+      enddo
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine newconf_residue_hairpin(idum,vvar,i1,fail)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      real ran1,ran2
+      dimension vvar(mxang,maxres,mxch),iold(ntotal)
+      integer nharp,iharp(4,maxres/3),icipa(maxres/3)
+      logical fail,not_done
+      ctdif=10.
+
+      fail=.false.
+      if (iseed.gt.mxio .or. iseed.lt.1) then
+        write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
+        call mpi_abort(mpi_comm_world,ierror,ierrcode)
+      endif
+      do k=1,numch
+       do j=2,nres-1
+        do i=1,4
+         vvar(i,j,k)=bvar(i,j,k,iseed)
+        enddo
+       enddo
+      enddo
+      do k=1,numch
+      do j=2,nres-1
+        theta(j+1)=bvar(1,j,k,i1)
+        phi(j+2)=bvar(2,j,k,i1)
+        alph(j)=bvar(3,j,k,i1)
+        omeg(j)=bvar(4,j,k,i1)
+      enddo
+      enddo
+c      call intout
+      call chainbuild
+      call hairpin(.false.,nharp,iharp)
+
+      if (nharp.eq.0) then
+        fail=.true.
+        return
+      endif
+
+      n_used=0
+
+      DO III=1,NHARP
+
+      not_done = .true.
+      icount=0
+      do while (not_done)
+        icount=icount+1
+        iih=iran_num(1,nharp)
+        do k=1,n_used
+          if (iih.eq.icipa(k)) then
+            iih=0
+            goto 22
+          endif
+        enddo
+        not_done=.false.
+        n_used=n_used+1
+        icipa(n_used)=iih
+   22   continue
+        not_done = not_done .and. icount.le.nharp
+      enddo
+
+      if (iih.eq.0) then
+        write (iout,*) "CHUJ NASTAPIL W NEWCONF_RESIDUE_HAIRPIN!!!!"
+        fail=.true.
+        return
+      endif
+      
+      istart=iharp(1,iih)+1
+      iend=iharp(2,iih)
+
+cdd      write (iout,*) "newconf_residue_hairpin: iih",iih,
+cdd     &   " istart",istart," iend",iend
+
+      do k=1,numch
+      do j=istart,iend
+       do i=1,4
+         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
+         if(dif.gt.180.) dif=360.-dif
+         if(dif.gt.ctdif) goto 20
+       enddo
+      enddo
+      enddo
+      goto 10
+   20 continue
+      do k=1,numch
+      do j=istart,iend
+       do i=1,4
+        vvar(i,j,k)=bvar(i,j,k,i1)
+       enddo
+      enddo
+      enddo
+c      do j=1,numch
+c       do l=2,nres-1
+c        write (iout,'(4f8.3)') (rad2deg*vvar(i,l,j),i=1,4)
+c       enddo
+c      enddo
+      return
+   10 continue
+      ENDDO
+
+      fail=.true.
+
+      return
+      end
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine gen_hairpin
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.HAIRPIN'
+
+c      write (iout,*) 'Entering GEN_HAIRPIN'
+      do iters=1,nseed
+        i1=is(iters)
+        do k=1,numch
+          do j=2,nres-1
+            theta(j+1)=bvar(1,j,k,i1)
+            phi(j+2)=bvar(2,j,k,i1)
+            alph(j)=bvar(3,j,k,i1)
+            omeg(j)=bvar(4,j,k,i1)
+          enddo
+        enddo
+        call chainbuild
+        call hairpin(.false.,nharp_seed(iters),iharp_seed(1,1,iters))
+      enddo
+
+      nharp_tot=0
+      do iters=1,nseed
+        nharp_tot=nharp_tot+nharp_seed(iters)
+        nharp_use(iters)=4*nharp_seed(iters)
+        do j=1,nharp_seed(iters)
+          iharp_use(0,j,iters)=4
+          do k=1,4
+            iharp_use(k,j,iters)=0
+          enddo
+        enddo
+      enddo
+
+      write (iout,*) 'GEN_HAIRPIN: nharp_tot',nharp_tot
+cdd      do i=1,nseed
+cdd        write (iout,*) 'seed',i 
+cdd        write (iout,*) 'nharp_seed',nharp_seed(i),
+cdd     &     ' nharp_use',nharp_use(i)
+cd        write (iout,*) 'iharp_seed, iharp_use'
+cd        do j=1,nharp_seed(i)
+cd          write (iout,'(7i3)') iharp_seed(1,j,i),iharp_seed(2,j,i),
+cd     &      (iharp_use(k,j,i),k=0,4) 
+cd        enddo
+cdd      enddo
+      return
+      end
+
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+ccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine select_frag(nn,nh,nl,ns,nb,i_csa)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CSA'
+      include 'COMMON.BANK'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.HAIRPIN'
+      include 'COMMON.DISTFIT'
+      character*50 linia
+      integer isec(maxres)
+
+
+      nn=0
+      nh=0
+      nl=0
+      ns=0
+      nb=0
+cd      write (iout,*) 'Entering select_frag'
+      do i1=1,nbank
+        do i=1,nres
+          isec(i)=0
+        enddo
+        do k=1,numch
+          do j=2,nres-1
+            theta(j+1)=bvar(1,j,k,i1)
+            phi(j+2)=bvar(2,j,k,i1)
+            alph(j)=bvar(3,j,k,i1)
+            omeg(j)=bvar(4,j,k,i1)
+          enddo
+        enddo
+        call chainbuild
+cd        write (iout,*) ' -- ',i1,' -- ' 
+        call secondary2(.false.)
+c
+c bvar_frag nn==pair of nonlocal strands in beta sheet (loop>4)
+c               strands > 4 residues; used by N7 and N16
+c
+        do j=1,nbfrag
+c
+Ctest 09/12/02 bfrag(2,j)-bfrag(1,j).gt.3
+c
+              do i=bfrag(1,j),bfrag(2,j)
+                isec(i)=1
+              enddo
+              do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
+                isec(i)=1
+              enddo
+
+           if ( (bfrag(3,j).lt.bfrag(4,j) .or.
+     &          bfrag(4,j)-bfrag(2,j).gt.4) .and. 
+     &          bfrag(2,j)-bfrag(1,j).gt.4 ) then
+             nn=nn+1
+
+
+             if (bfrag(3,j).lt.bfrag(4,j)) then
+               write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') 
+     &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
+     &                         ,",",bfrag(3,j)-1,"-",bfrag(4,j)-1
+             else
+               write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') 
+     &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
+     &                         ,",",bfrag(4,j)-1,"-",bfrag(3,j)-1
+
+             endif
+cd             call write_pdb(i_csa*1000+nn+nh,linia,0d0)
+
+             bvar_frag(nn,1)=i1
+             bvar_frag(nn,2)=4
+             do i=1,4
+               bvar_frag(nn,i+2)=bfrag(i,j)
+             enddo
+           endif
+        enddo
+
+c
+c hvar_frag nh==helices; used by N8 and N9
+c
+        do j=1,nhfrag
+
+              do i=hfrag(1,j),hfrag(2,j)
+                isec(i)=2
+              enddo
+
+          if ( hfrag(2,j)-hfrag(1,j).gt.4 ) then
+            nh=nh+1
+
+cd            write(linia,'(a6,i3,a1,i3)')
+cd     &                "select",hfrag(1,j)-1,"-",hfrag(2,j)-1
+cd            call write_pdb(i_csa*1000+nn+nh,linia,0d0)
+
+            hvar_frag(nh,1)=i1
+            hvar_frag(nh,2)=hfrag(1,j)
+            hvar_frag(nh,3)=hfrag(2,j)
+          endif
+        enddo
+
+        
+cv        write(iout,'(i4,1pe12.4,1x,1000i1)') 
+cv     &         i1,bene(i1),(isec(i),i=1,nres)
+cv            write(linia,'(i4,1x,1000i1)') 
+cv     &         i1,(isec(i),i=1,nres)
+cv            call write_pdb(i_csa*1000+i1,linia,bene(i1))
+c
+c lvar_frag nl==loops; used by N14
+c
+        i=1
+        nl1=nl
+        do while (i.lt.nres)
+          if (isec(i).eq.0) then
+           nl=nl+1
+           lvar_frag(nl,1)=i1
+           lvar_frag(nl,2)=i
+           i=i+1
+           do while (isec(i).eq.0.and.i.le.nres)
+             i=i+1
+           enddo
+           lvar_frag(nl,3)=i-1
+           if (lvar_frag(nl,3)-lvar_frag(nl,2).lt.1) nl=nl-1
+          endif
+          i=i+1
+        enddo
+cd        write(iout,'(4i5)') (i,(lvar_frag(i,ii),ii=1,3),i=nl1+1,nl)
+
+c
+c svar_frag ns==an secondary structure element; used by N15
+c
+        i=1
+        ns1=ns
+        do while (i.lt.nres)
+          if (isec(i).gt.0) then
+           ns=ns+1
+           svar_frag(ns,1)=i1
+           svar_frag(ns,2)=i
+           i=i+1
+           do while (isec(i).gt.0.and.isec(i-1).eq.isec(i)
+     &                         .and.i.le.nres)
+             i=i+1
+           enddo
+           svar_frag(ns,3)=i-1
+           if (svar_frag(ns,3)-svar_frag(ns,2).lt.1) ns=ns-1
+          endif
+          if (isec(i).eq.0) i=i+1
+        enddo
+cd        write(iout,'(4i5)') (i,(svar_frag(i,ii),ii=1,3),i=ns1+1,ns)
+
+c
+c avar_frag nb==any pair of beta strands; used by N17
+c
+        do j=1,nbfrag
+             nb=nb+1
+             avar_frag(nb,1)=i1
+             do i=1,4
+               avar_frag(nb,i+1)=bfrag(i,j)
+             enddo
+        enddo
+
+      enddo
+
+      return
+      end