program linear

! Program for GCMC simulation for one-component chains of arbitrary sequence
! Compilation for production: gfortran -ffree-line-length-none -O2 linear.f90 -o linear.out  [-Ofast leads to seg faults]
! Compilation for debugging:  gfortran -ffree-line-length-none -g -ffpe-trap=denormal,underflow,invalid,overflow linear.f90 -o linear_dbg.out
! Version 1.63 Apr '24 (AZP) - fixed att_tr issue in long runs
! Version 1.62 Apr '24 (AZP) - fixed nequil, nprod issue to allow for long runs
! Version 1.61 Nov '23 (AZP) - changed cluster_count to include all bead types; no clu output
! Version 1.6 Nov '23 (AZP) - renumbered versions to reflect bigger changes, fixed speed of flip move
! Version 1.05, Oct 23 (AZP) - added flip routine for in-box perturbations
! Version 1.04, Oct 23 (AZP) - added transfer3 routine for in-box displacements, z-profile
! Version 1.03, Feb 23 (AZP) - adapted from gcmc_m2 with the following major changes
!   1. removed all remnants of multicomponent code
!   2. transfer routine now does "thermal" Rosenbluth growth
!   3. removed limitation of integer*4
! Following lines are for gcmc_m2
! Version 1.49, Jan 23 (AZP) - output format to screen and input sequence from file
! Version 1.48, Dec 22 (AZP) - output format to cha.dat
! Version 1.47, Oct 22 (AZP) - output to xyz to ensure ovito compatibility
! Version 1.46, May 22 (AZP) - fixed issue with long runs with low N
! Version 1.45, Jul 18 (AZP) - cluster move only if systems is not a single cluster
! Version 1.44, Jul 18 (AZP) - formats, compile: gfortran -ffree-line-length-none  -finteger-4-integer-8 -Ofast gcmc_m2.f90
! Version 1.43, Dec 06 (AZP) - cluster counting routines fixed
! Version 1.42, March 05 (AZP) - minor changes in I/O 
! Version 1.41, May 01 (AZP) - fixed problem with N_repeat > 1 
! Version 1.4, March 01 (AZP) - partial regrowth move 
! Version 1.32, July 00 (AZP) - ability to do z=6 as well as z=26 
! Version 1.31, May 00 (AZP)
! Based on v. 1.1 of gcmc_m
! Version 1.0, June 99 (AZP) based on Version 1.24 of gcmc_r
! design changes:  sequence can be arbitrary (e.g. HTHTHH)
! Following lines are for gcmc_r
! Version 1.24,  September '98 (AZP) 
! Differences from version 1.1:
! (a) Fixed slight problem with acceptance when many configurations have no particles
! (b) Added possibility for z-profile calculation in lamellar phases
! (c) Input data read from screen: run as "gcmc_r.exe < inp.dat > screen.dat"
! (d) Dynamic allocation of arrays
! Difference of v 1.1 from v. 1.0:  rosenbluth weight calc. correctly for long chains
! Based on version 1.2 of amphiphile
!	Differences from version 1.0
!		1. reordering of arrays (10% efficiency gain)
!		2. reptation move
!	Differences from version 1.1: Cluster move
! Design notes:  last species is always the longest amphiphile in the system,
! first species the shortest - for exchange move purposes
! Max number of steps: 2,147,483,647 (due to integer*4 limitation)

implicit none

integer, parameter ::	dim = 3				! dimensionality of space
integer, parameter ::	nmoves = 6			! types of moves: 1 = transfer (ath), 2 = transfer (th),
                                            ! 3 = in-box transfer (ath) 4 = flip 5 = regrow part (ath), 6 = cluster move
integer                 ngroups             ! if 2, H and T; if 3, HH, HT and TT; "water" is equiv. to lattice vacancies
integer					xdim,ydim,zdim		! box dimensions in x, y, z direction
integer					maxmolec			! maximum number of molecules (of any type) in system
integer					n_energy_hist		! number of points for energy discretization
real*8					width_energy_hist	! width of energy histograms; should take into account interactions
real*8, allocatable ::	energy_hist(:,:)	! element 0 of this matrix contains the
											! starting value of energy for the bins at that amph. number. 
											! element -2 is the min value of energy for non-zero observations
											! element -1 is the max value of energy for non-zero observations
integer					hist_index			! position in histogram
integer					max_ncluster		! maximum number of beads in a cluster
integer					nbead		        ! total number of beads in chain
integer, allocatable::  zprof(:,:)          ! [zdim,ngroups] density profile in z-direction
integer, allocatable::	bead_type(:)		! [nbead],group type of each bead
integer, allocatable::  count_bead(:)       ! beads of each group
integer					nmolec	            ! number of molecules present
integer					namph,iamph,istart,iend	! for backwards compatibility
integer*1,allocatable::	space(:,:,:)		! [xdim,ydim,zdim]
											! contains 0 for vacancy, + 1 for H, + 2 for T etc
											! values less than 0 are used for counting clusters
											! note that integer*1 has a range of -128 to +127
real*8, allocatable ::	interact(:,:)	    ! contains the interactions between groups
											! group 0 is vacuum!
integer, allocatable::	xyz(:,:,:)			! [dim,nbead,maxmolec] positions of beads of each molecule
integer					nvectors			! = 26	or 6, nearest neighbor vectors 
integer, allocatable::	vector(:,:)			!(dim,nvectors)! neighbor vectors 
integer					length(dim)			! length of the space array along respective dimensions
real					ran2				! random number function
integer					seed,seed0			! random number function seed, initial seed
real*8					beta				! inverse temperature
real					prob_moves(nmoves)	! probability of each move type (cummulative)
integer*8				nequil,nproduct,istep		! number of equilibration and production steps, step count
integer					i,j,idim,imolec,ibead,ix,iy,iz,zpos ! local counters
real*8                  znorm,zshift               ! used to normalized z position
integer					ncluster        	! number of beads in cluster for a species
integer					icluster        	! number of molecules in cluster for a species
integer, allocatable :: cluster_bin(:,:)	! bins cluster size observations; from the way cluster_count is called, large
											! clusters are sampled more frequently than small, as the probability of
											! finding a cluster is proportional to the aggregation number.  To convert cluster_bin
											! to frequency of clusters, they must be divided by their size.
integer, allocatable :: chains_bin(:)		! bins cluster size observations at the chain level
integer					n_sample_cluster	! counts every how many steps a cluster sample will be taken
integer					isnap				! snapshot counter
integer					n_snap				! when cluster routine is called, a snapshot is printed to file 34 every n_snap times
logical					success				! was move successful?
real*8					e_total,energy_total	! total energy, function for energy
integer					idim1,idim2,idim3,ivectors	! local counters
real*8					e_check				! check final energy
real*8					e_self          	! energy along backbone of chain
real					step_type			! type of move (transfer1, transfer2, regrowth, cluster)
integer*8				att_tr1,att_tr2,att_tr3,att_fl	    ! attempted transfers
integer*8				succ_tr1,succ_tr2,succ_tr3,succ_fl  ! successful transfers
integer*8				att_regr,succ_regr  ! attempted and successful regrowth moves
integer*8				att_clust,succ_clust! attempted and successful cluster moves
integer, parameter ::	nblocks=8			!  number of blocks for block averages
integer					iblock				!  count at which block you are
integer*8				samples_per_block	! how many samples per block
real*8					eblock(nblocks),e_avg,e_std	! energy block averages and statistics
real*8					nmolec_block(nblocks)
real*8					nmolec_avg,nmolec_std! number block averages and statistics
integer, parameter ::	nprint=4			! how many times xyz and cluster distributions are printed
integer*8				n_print				! how many steps between printing of xyz
integer*8				n_print_all			! how many steps between printing of energy and clusters
logical					is_a_continuation	! is run a continuation of previous run?
real*8					m_equil,m_product	! million of equilibration, production steps
real*8					mu		            ! chemical potential of molecules, with respect to non-self avoiding, reversing random walk
integer					istat				! aux variable
character*20			hostname,auxfile
integer					n_repeat,i_repeat	! how many replications of the run are desired, counter
character*30			cont_file			! name of continuation file
character*30            seq_file            ! name of continuation file
character*5				run_label			! label for output files
integer					iostat				! I/O status
real(4)					itim, ttt, tarr(2)	! timing stuff
integer					auxint1,auxint2		! auxiliary integer variables
real					phi,phi_std			! volume fraction, std. deviation 
integer					nsegments       	! number of segments
integer, allocatable :: nbeadseg(:)		    ! [nsegments] number of beads in each segment
integer					beadcount			! counts beads
character(LEN=30)       FMT1                ! for variable format declarations

! START EXECUTABLE PART OF CODE

itim = dtime(tarr)  ! CPU time

! read input data and output to screen
read (*,*) xdim,ydim,zdim,ngroups,nvectors,maxmolec
if (nvectors/=6.and.nvectors/=26) then
	write (*,*) ' Only Z=6 and Z=26 are supported'
	stop
endif
allocate(count_bead(ngroups))
allocate(interact(0:ngroups,0:ngroups))
allocate(zprof(zdim,ngroups))
zprof = 0
znorm = 0

read (*,*) seq_file
open (3,file=trim(seq_file),status='old')
read (3,*) nbead
nsegments = 1  ! works only with single species
allocate (bead_type(nbead))
count_bead = 0
read(3,*) bead_type(1)
count_bead(bead_type(1)) = count_bead(bead_type(1)) + 1
do j = 2,nbead
    read(3,*) bead_type(j)
    count_bead(bead_type(j)) = count_bead(bead_type(j)) + 1
    if (bead_type(j) .ne. bead_type(j-1)) nsegments = nsegments + 1
enddo
allocate(vector(dim,nvectors))

maxmolec = min(maxmolec,int((xdim*ydim*zdim)/nbead))

allocate (space(xdim,ydim,zdim),xyz(dim,nbead,maxmolec),cluster_bin(maxmolec*nbead,0:maxmolec))
allocate (chains_bin(maxmolec))
do idim=1,dim
	length(idim) = size(space,idim)  ! initialize length
enddo
! interactions with group 0 (vacancies) are always zero.
interact = 0.
read (*,*) ((interact(i,j),i=j,ngroups),j=1,ngroups)
do j=1,ngroups
	do i=j+1,ngroups
		interact(j,i) = interact(i,j)
	enddo
enddo
read (*,*) n_energy_hist,width_energy_hist
allocate (energy_hist(-2:n_energy_hist,0:maxmolec))
read (*,*) beta,mu
beta = 1/beta ! since input file says "temperature"
read (*,*) (prob_moves(i),i=1,nmoves-1),n_sample_cluster,n_snap,n_repeat
if (nbead<5 .and. prob_moves(nmoves-1)-prob_moves(nmoves-2)>0.) then
	write (*,*) ' Regrowth moves should not be used for short species'
	stop
endif
read (*,*) m_equil,m_product,seed0
seed0 = -abs(seed0)	! to initialize ran2 properly
nequil = nint(m_equil*1000000,8)
nproduct = nint(m_product*1000000,8)
read (*,*) cont_file
read (*,*) run_label

! ----- end reading data -----
! Echo to screen
write (*,*)	'------------------------------------------------------------------'
write (*,*) ' GCMC of linear chains: Version 1.63 - Apr 2024 (AZP)'
write (*,*)	'------------------------------------------------------------------'
auxint2 = int(log10(10.*maxval(length)))+1
write(FMT1,'("(A,3I", I0, ",A,i2,3A,i5,i5,f6.1)")') auxint2
write (*,FMT1) ' Lattice:',length,' Z=',nvectors,'  Run=',trim(run_label),', Max/pts/w=',maxmolec, &
		 n_energy_hist,width_energy_hist
write (*,'(3A,i4,A,3f6.1)') ' Architecture: ',trim(seq_file), ' nseg=',nsegments, &
            ' groups % ',(100./nbead*count_bead(i),i=1,ngroups)
write (*,'(A,f7.3,f9.3,A,6f7.3)') ' T,μ =',1./beta,mu,' I:',((interact(i,j),i=j,ngroups),j=1,ngroups)
write (*,*) 'P_tr1 tr2 tr3  flip regr  10^6x Eq.     Prod. Clu_sample Snap Seed'
write (*,'(4f5.2,f6.3,f10.2,f11.2,2i8,i4)') (prob_moves(i),i=1,nmoves-1),m_equil,m_product,n_sample_cluster,n_snap,seed0
!------ end echo to screen ------
open (24,file='scr'//trim(run_label)//'.dat') !,status='new')
write (24,*)	'------------------------------------------------------------------'
write (24,*) ' GCMC of linear chains: Version 1.63 - Nov 2024 (AZP)'
write (24,*)	'------------------------------------------------------------------'
write (24,FMT1) ' Lattice:',length,' Z=',nvectors,'  Run=',trim(run_label),', Max/pts/w=',maxmolec, &
		 n_energy_hist,width_energy_hist
write (24,'(3A,i4,A,3f6.1)') ' Architecture: ',trim(seq_file), ' nseg=',nsegments, &
            ' groups % ',(100./nbead*count_bead(i),i=1,ngroups)
write (24,'(A,f7.3,f9.3,A,6f7.3)') ' T,μ =',1./beta,mu,' I:',((interact(i,j),i=j,ngroups),j=1,ngroups)
write (24,*) 'P_tr1 tr2 tr3  flip regr  10^6x Eq.     Prod. Clu_sample Snap Seed'
write (24,'(4f5.2,f6.3,f10.2,f11.2,2i8,i4)') (prob_moves(i),i=1,nmoves-1),m_equil,m_product,n_sample_cluster,n_snap,seed0
!end echo to screen file

samples_per_block = int(nproduct/nblocks)
n_print = int((nequil+nproduct)/nprint)
if (n_print == 0) n_print = 1


! self energy is due to interactions between neighboring beads on the same chain;
! it is subtracted from the total energy when the energy is printed or averaged

e_self = 0
do ibead = 1,nbead-1
    e_self = e_self + interact(bead_type(ibead),bead_type(ibead+1))
enddo

! initialize vectors
if (nvectors==26) then
	ivectors = 0
	do idim1 = -1,1
		do idim2 = -1,1
			do idim3 = -1,1	!#3D
				if (abs(idim1)+abs(idim2)+abs(idim3)>0) then	!#3D
	!			if (abs(idim1)+abs(idim2)>0) then	!#2D
					ivectors = ivectors + 1
					vector(1:dim,ivectors)  = (/ idim1, idim2, idim3/)	!#3D
	!				vector(1:dim,ivectors)  = (/ idim1, idim2/)	!#2D
				endif
			enddo	!#3D
		enddo
	enddo
else ! Z=6
	vector(1:dim,1)  = (/ 1, 0, 0/)	
	vector(1:dim,2)  = (/ 0, 1, 0/)	
	vector(1:dim,3)  = (/ 0, 0, 1/)	
	vector(1:dim,4)  = (/ -1, 0, 0/)	
	vector(1:dim,5)  = (/ 0, -1, 0/)	
	vector(1:dim,6)  = (/ 0, 0, -1/)	
endif


!--------------------- NO INDENT OVER N_REPEAT LOOP -----------
do i_repeat = 1,n_repeat

nmolec = 0 ! always start with zero chains, unless continuation from previous run
eblock = 0
nmolec_block = 0.
cluster_bin = 0
seed = seed0 + i_repeat - 1
energy_hist = 0.
xyz = 0
space = 0 ! initialize space by filling it with water
iostat = 0
open (3,file=trim(cont_file),status='old',iostat=iostat)
if (iostat == 0) then
	is_a_continuation = .true.
else
	is_a_continuation = .false.
endif
if (is_a_continuation) then
	read (3,*) nmolec
	nmolec=nmolec/nbead
    read (3,*)
	write (*,'(3a,i4)') '  Continuation from file ',trim(cont_file),', nmolec=',nmolec
	write (24,'(3a,i4)') '  Continuation from file ',trim(cont_file),', nmolec=',nmolec
    do imolec = 1,nmolec
        do ibead = 1,nbead
            read (3,'(2x,3i5)') xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)
            if (space(xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)) /= 0) then
                write (*,*) '  Overlaps in file: molec/bead: ',trim(cont_file),imolec, ibead
                write (24,*) '  Overlaps in file: molec/bead: ',trim(cont_file),imolec, ibead
                stop
            endif
            space(xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)) = bead_type(ibead)
        enddo
    enddo
	close(3)
endif
! subtract self energy
e_total = energy_total(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,nmolec,nbead,interact,xyz) - nmolec*e_self

att_tr1=0; succ_tr1=0; att_tr2=0; succ_tr2=0; att_tr3=0; succ_tr3=0; att_fl=0; succ_fl=0; att_clust=0; succ_clust=0
att_regr=0; succ_regr=0
max_ncluster = nbead*nmolec
isnap = 0
eblock = eblock/samples_per_block
nmolec_block = nmolec_block/samples_per_block
do istep = 1,nequil+nproduct
	if (istep == nequil + 1) then
		if (succ_tr1+succ_tr2/= 0) then
			! the choice of n_print_all is somewhat arbitrary
			! More molecules, lower acceptance ratio -> higher n_print_all
			n_print_all = (att_tr1+att_tr2)/(succ_tr1+succ_tr2)
			n_print_all = min(n_print_all,95000)	! nothing magical about 95,000
		else
			n_print_all = 95000	! if there were no sucessful transfers in equil.
		endif
		! reset counters
        att_tr1=0; succ_tr1=0; att_tr2=0; succ_tr2=0; att_tr3=0; succ_tr3=0; att_fl=0; succ_fl=0; att_clust=0; succ_clust=0
        att_regr=0; succ_regr=0
		! open files
		if (n_sample_cluster/=0) then
			open (34,file='sna'//trim(run_label)//achar(96+i_repeat)//'.dat')	
			write (34,'(f6.2,f9.3,3i4,i3,i6,A)') 1./beta,mu,xdim,ydim,zdim,nbead,n_snap,' ! T, mu, x,y,zdim, nbead, n_snap'
			write (34,'(A)') '  short lines give: # of molec. in cluster, total, sites in cluster'
		endif
	endif

	step_type = ran2(seed)
	if (step_type<=prob_moves(1)) then	! transfer in / out (athermal)
		call transfer1(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
				nmolec,nbead,bead_type,interact,e_self,xyz,success,seed,beta,e_total,mu,istep)
		att_tr1 = att_tr1 + 1
		if (success) succ_tr1 = succ_tr1 + 1
! DDD
!        write (*,'(i7,A,L2,i5,f11.2,4i6)') istep, ' tr1',success, nmolec, e_total
!          e_check = energy_total(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,nmolec,nbead,interact,xyz)
!          e_check = e_check - nmolec*e_self ! subtract self energy
!          if (e_check/=e_total) then
!               write (*,*) ' **** Energy does not check: nmolec,e_check, e_total, e_self **** '
!               write (*,'(2i7,3f12.2)') istep, nmolec,e_check,e_total,e_self    !debug
!               stop
!         endif
! end DDD

    else if (step_type<=prob_moves(2)) then    ! transfer in / out (thermal)
        call transfer2(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
            nmolec,nbead,bead_type,interact,e_self,xyz,success,seed,beta,e_total,mu)
        att_tr2 = att_tr2 + 1
        if (success) succ_tr2 = succ_tr2 + 1

    else if (step_type<=prob_moves(3).and.nmolec>0) then    ! transfer in box
        call transfer3(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
                nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)
        att_tr3 = att_tr3 + 1
        if (success) succ_tr3 = succ_tr3 + 1

    else if (step_type<=prob_moves(4).and.nmolec>0) then    ! flip
        call flip(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,&
                nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)
        att_fl = att_fl + 1
        if (success) succ_fl = succ_fl + 1

    else if (step_type<=prob_moves(5).and.nmolec>0) then	! partial regrowth
		call regrowth(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
                nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)
		att_regr = att_regr + 1
		if (success) succ_regr = succ_regr + 1

	else if (nmolec>0) then ! cluster moves
		max_ncluster = nbead*nmolec
		call cluster_move(dim,xdim,ydim,zdim,maxmolec,length,nvectors,vector,nmolec,nbead,xyz,space,max_ncluster,success,seed,istep)
		att_clust = att_clust + 1
		if (success) succ_clust = succ_clust + 1
	endif



	if (mod(istep,n_print)==0) then		! print xyz
		auxfile = trim(run_label)//achar(96+i_repeat)//'.xyz'
		open (21,file=auxfile)	
		max_ncluster = nbead*nmolec
		write (21,'(i8)') max_ncluster
		write (21,'(A,9i4,A)') 'Lattice="',xdim,0,0,0,ydim,0,0,0,zdim,'" Properties=S:1:pos:R:3"'
        do imolec = 1,nmolec
            do ibead = 1,nbead
                if (bead_type(ibead) == 1) then
                    write (21,'(A,3i5)') 'N ',xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)
                else if (bead_type(ibead) == ngroups) then
                    write (21,'(A,3i5)') 'O ',xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)
                else
                    write (21,'(A,3i5)') 'C ',xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)
                endif
            enddo
        enddo
		close (21)
	endif

	if (istep>nequil) then	! sample properties

		if (mod(istep,n_print_all)==0) then		! sample histograms
			namph = nmolec
			if (energy_hist(-2,namph) == energy_hist(0,namph)) then	! first energy observed for this namph - set
																	! at mid-point
				energy_hist(0,namph) = width_energy_hist*int(e_total/width_energy_hist - n_energy_hist/2)
				energy_hist(-2,namph) = energy_hist(0,namph) - width_energy_hist*(int((energy_hist(0,namph)-e_total)/width_energy_hist)+1)
				energy_hist(-1,namph) = energy_hist(0,namph) + width_energy_hist*(int((-energy_hist(0,namph)+e_total)/width_energy_hist)+1)
			endif
			if (e_total < energy_hist(-2,namph)) then	! found lower value than min
				energy_hist(-2,namph) = energy_hist(0,namph) - width_energy_hist*(int((energy_hist(0,namph)-e_total)/width_energy_hist)+1)
				if (energy_hist(-2,namph) <= energy_hist(0,namph) + width_energy_hist ) then
					write (*,*) 'Maximum width of energy histograms exceeded (low)'
					write (*,*) 'istep,namph, energy(-2),energy(-1),energy(0),e_total', istep,namph,energy_hist(-2:0,namph),e_total
					write (24,*) 'Maximum width of energy histograms exceeded (low)'
					write (24,*) 'istep,namph, energy(-2),energy(-1),energy(0),e_total', istep,namph,energy_hist(-2:0,namph),e_total
					stop
				endif
			else if (e_total > energy_hist(-1,namph)) then	! found greater value than max
				energy_hist(-1,namph) = energy_hist(0,namph) + width_energy_hist*(int((-energy_hist(0,namph)+e_total)/width_energy_hist)+1)
				if (energy_hist(-1,namph) > energy_hist(0,namph)+width_energy_hist*(n_energy_hist-2)) then
					write (*,*) 'Maximum width of energy histograms exceeded (high)'
					write (*,*) 'istep,namph, energy(-2),energy(-1),energy(0),e_total', istep,namph,energy_hist(-2:0,namph),e_total
					write (24,*) 'Maximum width of energy histograms exceeded (high)'
					write (24,*) 'istep,namph, energy(-2),energy(-1),energy(0),e_total', istep,namph,energy_hist(-2:0,namph),e_total
					stop
				endif
			endif
			! value ok
			hist_index = int((e_total-energy_hist(0,namph))/width_energy_hist + 0.5) 
			energy_hist(hist_index,namph) = energy_hist(hist_index,namph) + 1
		endif

		iblock = int((istep-nequil)/samples_per_block) + 1
		if (iblock<=nblocks) then
			if(nmolec>0) eblock(iblock) = eblock(iblock) + e_total/nmolec
			nmolec_block(iblock) = nmolec_block(iblock) + nmolec
		endif
		if (n_sample_cluster/=0) then
			if (mod(istep-nequil,n_sample_cluster)==0.and.nmolec>0) then
				max_ncluster = nbead*nmolec
                call cluster_count(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,nmolec,&
						nbead,xyz,ncluster,icluster,max_ncluster,isnap,n_snap,seed)
				cluster_bin(ncluster,0) = cluster_bin(ncluster,0) + 1
				cluster_bin(ncluster,icluster) = cluster_bin(ncluster,icluster) + 1
                !  collect density profiles
                znorm = znorm+length(1)*length(2)
                zshift = 0.
!                do imolec = 1,nmolec   !  deactivated zshift / may reactivate at some point
!                    zshift = zshift+xyz(3,int(nbead/2)+1,imolec)
!                enddo
!                zshift = zshift/nmolec
!                if (zshift > length(3)/2.) zshift = length(3)/2. + 0.5 - zshift
                do imolec = 1,nmolec
                    do ibead = 1,nbead
                        zpos = xyz(3,ibead,imolec) + zshift
!                        if (zpos <=0) zpos = zpos + length(3)
!                        if (zpos >length(3)) zpos = zpos - length(3)
                        zprof(zpos,bead_type(ibead)) = zprof(zpos,bead_type(ibead)) + 1
                    enddo
                enddo
			endif
		endif
		if (mod(istep-nequil,samples_per_block)==samples_per_block-1) then
			if (iblock == 1) then
				if (n_repeat>1) then
					write (*,'(A,i3)') '  Repetition number ',i_repeat
					write (24,'(A,i3)') '  Repetition number ',i_repeat
				endif
				write (*,*) 'Block       <U/N>                <N>             N/Phi'
				write (24,*) 'Block       <U/N>                <N>              N/Phi'
			endif
			write (*,'(i3,5x,f10.2,10x,f10.2,8x,i7)') iblock,eblock(iblock)/samples_per_block &
							,nmolec_block(iblock)/samples_per_block,nmolec
			write (24,'(i3,5x,f10.2,10x,f10.2,8x,i7)') iblock,eblock(iblock)/samples_per_block &
							,nmolec_block(iblock)/samples_per_block,nmolec
		endif
	endif
enddo

!  Checking for various errors, once per run
e_check = energy_total(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,nmolec,nbead,interact,xyz)
e_check = e_check - nmolec*e_self ! subtract self energy
if(abs(e_check-e_total)>1.e-6*abs(e_total)) then
	write (*,*) ' **** PROBLEM WITH ENERGY CONSERVATION: e_check, e_total **** ',e_check,e_total
	write (24,*) ' **** PROBLEM WITH ENERGY CONSERVATION: e_check, e_total **** ',e_check,e_total
endif
! check for broken chains
do imolec = 1,nmolec
    do ibead = 1,nbead-1
        do idim = 1,dim
            i = abs (xyz(idim,ibead,imolec) - xyz(idim,ibead+1,imolec))
            if (i > length(idim)/2) i = length(idim) - i
            if (i > 1) then
                write (*,'(A,3i5)') ' **** Broken chain: imolec, ibead, dist **** ',imolec,ibead,i
            endif
        enddo
    enddo
enddo
! check for inconsistent space
do imolec = 1,nmolec
    do ibead = 1,nbead-1
        if (space(xyz(1,ibead,imolec),xyz(2,ibead,imolec),xyz(3,ibead,imolec)) /= bead_type(ibead)) &
        write (*,'(A,5i5)') ' **** Space is inconsistent: imolec, ibead, xyz **** ',imolec,ibead,(xyz(i,imolec,ibead),i=1,dim)
    enddo
enddo

open (23,file='his'//trim(run_label)//achar(96+i_repeat)//'.dat')

write (23,'(a,i7,a)') ' nprint=',n_print_all,'    T       mu          width     x- y- zdim  '
write (23,*) 1./beta,mu,width_energy_hist,xdim,ydim,zdim
do iamph = 0,maxmolec
	istart = (energy_hist(-2,iamph)-energy_hist(0,iamph))/width_energy_hist
	iend = (energy_hist(-1,iamph)-energy_hist(0,iamph))/width_energy_hist
	if (istart/=iend) then
		write (23,'(2i7,f15.4)') iamph,iend-istart+1,energy_hist(-2,iamph)
		write (23,'(8i10)') (nint(energy_hist(i,iamph)),i=istart,iend)
	endif
enddo
close (23)

if (n_sample_cluster/=0) then
!	open (22,file='clu'//trim(run_label)//achar(96+i_repeat)//'.dat')
!	write (22,'(A,f6.2,f9.3,5(1x,i0),A)') '/*', 1./beta,mu,xdim,ydim,zdim,nbead,n_sample_cluster,' ! T, mu, x,y,zdim, nbead, n_sample_cluster  */'
!	write (22,'(A,f7.2,A,f6.2,A,f10.6,A8,f8.3,A)') '/* Istep = ',real(istep)/1000000, &
!		 '  (million); T = ',1./beta,'; Rho @ step = ', real(nbead*nmolec)/size(space),' ; Mu = ',mu,' */'
 !   write (22,'(A)') '/* columns: # of beads in cluster; # observed; start(chains); counts  */'

!	do i = 1,size(cluster_bin,1)
!		if (cluster_bin(i,0)/=0) then
!			istart = 1
!			iend = maxmolec
!			do while (cluster_bin(i,istart)==0 .and. istart <= maxmolec)
!				istart = istart + 1
!			enddo
!			do while (cluster_bin(i,iend)==0 .and. iend >= 1)
!				iend = iend - 1
!			enddo
!			write (22,'(i0,1x,i0,1x,i0,60(1x,i0))') i,cluster_bin(i,0),istart,cluster_bin(i,istart:iend)
!		endif
!	enddo
!	close (22)

	open (22,file='cha'//trim(run_label)//achar(96+i_repeat)//'.dat')
	write (22,'(A,f6.2,f9.3,3i4,i3,i5,A)') '/*', 1./beta,mu,xdim,ydim,zdim,nbead,n_sample_cluster,' ! T, mu, x,y,zdim, nbead, n_sample_cluster */'
	write (22,'(A,f7.2,A,f6.2,A,f10.6,A8,f8.3,A)') '/* Istep = ',real(istep)/1000000, &
		 '  (million); T = ',1./beta,'; Rho @ step = ',&
		 real(nbead*nmolec)/size(space),' ; Mu = ',mu,' */'
	write (22,'(A)') '/* columns are: number of chains in cluster; # observations  */'
	chains_bin = 0
	do i = 1,maxmolec
		do j = 1,size(cluster_bin,1)
			chains_bin(i) = chains_bin(i) + cluster_bin(j,i)
		enddo
		if (chains_bin(i) /= 0) write (22,'(i4,i8)') i,chains_bin(i)
	enddo
	close (22)

endif

eblock = eblock/samples_per_block 
nmolec_block = nmolec_block/samples_per_block
call stats(nblocks,eblock,e_avg,e_std)
call stats(nblocks,nmolec_block(1:nblocks),nmolec_avg,nmolec_std)
phi = nmolec_avg*nbead/size(space)
phi_std = nmolec_std*nbead/size(space)
write (*,'(A,5x,f10.2,A,f5.2,1x,f10.2,A,f5.2,f10.4,A,f5.4)') ' <>',e_avg,' +- ',e_std,nmolec_avg,' +- ',nmolec_std,phi,' +- ',phi_std
write (24,'(A,5x,f10.2,A,f5.2,1x,f10.2,A,f5.2,f10.4,A,f5.4)') ' <>',e_avg,' +- ',e_std,nmolec_avg,' +- ',nmolec_std,phi,' +- ',phi_std

if (att_tr1 > 2E6 .or. att_tr2 > 2E6 .or. att_tr3 > 2E6 .or. att_fl > 2E6 .or. att_clust > 10000) then
    if (att_tr1 > 0) then
    write (*,'(A,f8.3,A,f8.3,A)') ' Transfers(ath) (10^6):',(dble(att_tr1)/1000000.),' att., ',(dble(succ_tr1)/1000000.),'  succ.'
	write (24,'(A,f8.3,A,f8.3,A)') ' Transfers(ath) (10^6):',(dble(att_tr1)/1000000.),' att., ',(dble(succ_tr1)/1000000.),'  succ.'
    endif
    if (att_tr2 > 0) then
    write (*,'(A,f8.3,A,f8.3,A)') ' Transfers(th) (10^6) :',(real(att_tr2)/1000000),' att., ',(real(succ_tr2)/1000000),'  succ.'
    write (24,'(A,f8.3,A,f8.3,A)') ' Transfers(th) (10^6) :',(real(att_tr2)/1000000),' att., ',(real(succ_tr2)/1000000),'  succ.'
    endif
    if (att_tr3 > 0) then
    write (*,'(A,f8.3,A,f8.3,A)') ' Int. Tr. (ath) (10^6):',(real(att_tr3)/1000000),' att., ',(real(succ_tr3)/1000000),'  succ.'
    write (24,'(A,f8.3,A,f8.3,A)') ' Int. Tr. (ath) (10^6):',(real(att_tr3)/1000000),' att., ',(real(succ_tr3)/1000000),'  succ.'
    endif
    if (att_fl > 0) then
    write (*,'(A,f8.3,A,f8.3,A)') ' Seq. Flips     (10^6):',(real(att_fl)/1000000),' att., ',(real(succ_fl)/1000000),'  succ.'
    write (24,'(A,f8.3,A,f8.3,A)') ' Seq. Flips     (10^6):',(real(att_fl)/1000000),' att., ',(real(succ_fl)/1000000),'  succ.'
    endif
	if (att_regr>0) then
		write (*,'(A,f8.3,A,f8.3,A)') ' Regrowth moves (10^6):',real(att_regr)/1000000,' att., ',real(succ_regr)/1000000,'  succ.'
		write (24,'(A,f8.3,A,f8.3,A)') ' Regrowth moves (10^6):',real(att_regr)/1000000,' att., ',real(succ_regr)/1000000,'  succ.'
	endif
	if (att_clust>0) then
		write (*,'(A,f8.3,A,f8.3,A)') ' Clust. moves (10^4):  ',real(att_clust)/10000,' att., ',real(succ_clust)/10000,'  succ.'
		write (24,'(A,f8.3,A,f8.3,A)') ' Clust. moves (10^4):  ',real(att_clust)/10000,' att., ',real(succ_clust)/10000,'  succ.'
	endif

else
    if (att_tr1>0)  then
        write (*,'(A,i8,A,i8,A)')     ' Transfers(ath) : ',att_tr1,' attempted, ',succ_tr1,'  successful'
        write (24,'(A,i8,A,i8,A)')     ' Transfers(ath)   : ',att_tr1,' attempted, ',succ_tr1,'  successful'
    endif
    if (att_tr2>0)  then
        write (*,'(A,i8,A,i8,A)')     ' Transfers(th)  : ',att_tr2,' attempted, ',succ_tr2,'  successful'
        write (24,'(A,i8,A,i8,A)')     ' Transfers(th)   : ',att_tr2,' attempted, ',succ_tr2,'  successful'
    endif
    if (att_tr3>0)  then
        write (*,'(A,i8,A,i8,A)')     ' Int. Tr. (ath) : ',att_tr3,' attempted, ',succ_tr3,'  successful'
        write (24,'(A,i8,A,i8,A)')    ' Int. Tr. (ath) : ',att_tr3,' attempted, ',succ_tr3,'  successful'
    endif
    if (att_fl>0)  then
        write (*,'(A,i8,A,i8,A)')     ' Seq. Flips     : ',att_fl,' attempted, ',succ_fl,'  successful'
        write (24,'(A,i8,A,i8,A)')    ' Seq. Flips     : ',att_fl,' attempted, ',succ_fl,'  successful'
    endif
    if (att_regr>0)  then
		write (*,'(A,i8,A,i8,A)') ' Regrowth Moves : ',att_regr,' attempted, ',succ_regr,' successful'
		write (24,'(A,i8,A,i8,A)') ' Regrowth Moves : ',att_regr,' attempted, ',succ_regr,' successful'
	endif
	if (att_clust>0)  then
		write (*,'(A,i8,A,i8,A)') ' Cluster Moves  : ',att_clust,' attempted, ',succ_clust,' successful'
		write (24,'(A,i8,A,i8,A)') ' Cluster Moves  : ',att_clust,' attempted, ',succ_clust,' successful'
	endif
endif	
if (znorm>1) then
    open(22,file='zpr'//trim(run_label)//achar(96+i_repeat)//'.dat')
    write (22,'(A,f6.2,f9.3,5(1x,i0),A)') '/*', 1./beta,mu,xdim,ydim,zdim,nbead,n_sample_cluster,' ! T, mu, x,y,zdim, nbead, n_sample_cluster  */'
    write (22,'(A,f7.2,A,f6.2,A,f10.6,A8,f8.3,A)') '/* Istep = ',real(istep)/1000000, &
         '  (million); T = ',1./beta,'; Rho @ step = ', real(nbead*nmolec)/size(space),' ; Mu = ',mu,' */'
    do iz = 1,length(3)
        write(22,'(5f10.5)') zprof(iz,1:ngroups)/znorm
    enddo
endif
enddo !--------- enddo on i_repeat w/o indent -------------------------



itim = dtime(tarr)
istat = hostnm(hostname)

if (itim<300) then
	write (*,'(A,f6.1,2A)') ' Program run in', itim, ' CPU s, on ',trim(hostname)
	write (24,'(A,f6.1,2A)') ' Program run in', itim, ' CPU s, on ',trim(hostname)
else if (itim<3600) then
	write (*,'(A,f6.1,2A)') ' Program run in', itim/60, ' CPU min on ',trim(hostname)
	write (24,'(A,f6.1,2A)') ' Program run in', itim/60, ' CPU min on ',trim(hostname)
else
	write (*,'(A,f6.2,2A)') ' Program run in', itim/3600, ' CPU hr on ',trim(hostname)
	write (24,'(A,f6.2,2A)') ' Program run in', itim/3600, ' CPU hr on ',trim(hostname)
    close (24)
endif

stop
end


function energy_total(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,nmolec,nbead,interact,xyz)

! routine picks sums up the energy in the initial configuration
implicit none
integer,	intent(in) ::	dim,xdim,ydim,zdim,maxmolec
integer*1,	intent(in) ::	space(xdim,ydim,zdim)
integer,	intent(in) ::	length(dim)
integer,	intent(in) ::	nvectors
integer,	intent(in) ::	ngroups
integer,	intent(in) ::	nmolec,nbead
real*8,		intent(in) ::	interact(0:ngroups,0:ngroups)
integer,	intent(inout):: xyz(dim,nbead,maxmolec)
real*8 energy_total

integer pos(dim),pos_neighb(dim,nvectors) ! local auxiliary variable
integer i,iprev,ibead,ipos,imolec ! local counters

iprev = 0
energy_total = 0.
do imolec = 1,nmolec
    do ibead = 1,nbead
        pos(1:dim) = xyz(1:dim,ibead,imolec)
        call find_neighb(dim,length,nvectors,pos,pos_neighb)
        do ipos = 1,nvectors
            energy_total = energy_total + interact(space(pos(1),pos(2),pos(3)),&
                space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
        enddo
    enddo
enddo

energy_total = energy_total / 2.

return 
end


subroutine transfer1(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
				nmolec,nbead,bead_type,interact,e_self,xyz,success,seed,beta,e_total,mu,istep)

! routine inserts a molecule in new random position
! or removes a randomly selected one with equal probabilities

implicit none
integer,	intent(in) ::		dim,xdim,ydim,zdim,maxmolec,nbead
integer*1,	intent(inout) ::	space(xdim,ydim,zdim)
integer,	intent(in) ::		length(dim)
integer,	intent(in) ::		nvectors
integer,	intent(in) ::		vector(dim,nvectors)
integer,	intent(in) ::		ngroups
integer,	intent(inout) ::	nmolec
integer,	intent(in) ::		bead_type(nbead)
real*8,		intent(in) ::		interact(0:ngroups,0:ngroups)
real*8,		intent(in) ::		e_self
integer,	intent(inout) ::	xyz(dim,nbead,maxmolec)
logical,	intent(out) ::		success
integer,	intent(inout) ::	seed
real*8,		intent(in) ::		beta
real*8,		intent(inout) ::	e_total
real*8,		intent(in) ::		mu
integer*8,    intent(in) ::       istep

! local variables
real			ran2		! function
integer			pos(dim)								! local auxiliary variable
integer			i,ibead,idim,ipos,ivector	! local counters
integer			removed							! removed molecule
integer			pos_neighb(dim,nvectors)
integer			pos_history(dim,nbead)
real*8			e_old, e_new	! energy of old and new configurations
real*8			rosen	! Logarithm of Rosenbluth factor; set to special value of -99 / beta if zero
integer			n_unocc	! how many sites are unoccupied
integer			v_unocc(nvectors)	! unoccupied vectors
integer			bead_count	! counts how many beads were successfully inserted

success = .false.
if (ran2(seed)<0.5) then	! remove chain

	if (nmolec>=1) then

		! pick random chain
		removed = int(nmolec*ran2(seed)) + 1
		rosen = 0.
		e_old = - e_self	! since e_self is included later
		do ibead = nbead,1,-1
			pos(1:dim) = xyz(1:dim,ibead,removed)
			n_unocc = 0
			call find_neighb(dim,length,nvectors,pos,pos_neighb)
			do ipos = 1,nvectors
				if (ibead/=nbead.and.space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) n_unocc = n_unocc + 1
				e_old = e_old + interact(space(pos(1),pos(2),pos(3)),space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
			enddo
			space(pos(1),pos(2),pos(3)) = 0 ! delete entity	
			if (ibead/=nbead) rosen = rosen + log(real(n_unocc))-log(real(nvectors))
		enddo

		success = .true.

		! accept/reject deletion move
        if (-rosen+beta*(e_old-mu)<-99.) rosen = 99.+beta*(e_old-mu) ! avoids underflow
		if (nmolec/real(size(space))*exp(-rosen+beta*(e_old-mu)) > ran2(seed)) then
			!accept deletion move
			e_total = e_total - e_old
			! exchange last molecule w/ removed one
			xyz(1:dim,1:nbead,removed) = xyz(1:dim,1:nbead,nmolec)
			nmolec = nmolec - 1
		else  ! reject deletion move
			do ibead = 1,nbead  ! restore chain in its old position
				space(xyz(1,ibead,removed),xyz(2,ibead,removed),xyz(3,ibead,removed)) = bead_type(ibead)
			enddo
			success = .false.
		endif
	endif

else	! insert molecule
	
	bead_count = 0
	! pick random position
	do idim = 1,dim
		pos(idim) = int(ran2(seed)*length(idim)) + 1
	enddo
	success = .true.
    e_new = - e_self    ! since e_self is included later
    rosen = 0.
	if (space(pos(1),pos(2),pos(3))/=0) then
		success = .false.
	else
		space(pos(1),pos(2),pos(3))= bead_type(1)  ! place first bead
		bead_count = 1 ! fisrt bead placed
		pos_history(1:dim,bead_count) = pos(1:dim)
		call find_neighb(dim,length,nvectors,pos,pos_neighb)
		do ipos = 1,nvectors
			e_new = e_new + interact(space(pos(1),pos(2),pos(3)),space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
		enddo
        ibead = 2
        do while (ibead <= nbead .and. success)
            n_unocc = 0
            v_unocc(1) = 1
            call find_neighb(dim,length,nvectors,pos,pos_neighb)
            do ipos = 1,nvectors
                if (space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) then
                    n_unocc = n_unocc + 1
                    v_unocc(n_unocc) = ipos
                endif
            enddo
            if (n_unocc /= 0) then
                rosen = rosen + log(real(n_unocc)) - log(real(nvectors))
            else
                rosen = -99. / beta  ! small value for the exponential, but avoids underflow
            endif
            ivector = v_unocc(int(ran2(seed)*n_unocc) + 1) ! select random unocc. position
            do idim = 1,dim ! periodic boundary conditions
                pos(idim) = pos(idim) + vector(idim,ivector)
                if (pos(idim)>length(idim)) pos(idim) = 1
                if (pos(idim)==0) pos(idim) = length(idim)
            enddo
            if (space(pos(1),pos(2),pos(3))==0) then ! found empty spot
                space(pos(1),pos(2),pos(3))=bead_type(ibead)
                bead_count = ibead
                pos_history(1:dim,bead_count) = pos(1:dim)
                call find_neighb(dim,length,nvectors,pos,pos_neighb)
                do ipos = 1,nvectors
                    e_new = e_new + interact(space(pos(1),pos(2),pos(3)),&
                    space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
                enddo
            else
                success = .false.  ! overlapped
            endif
            ibead = ibead + 1
		enddo
	endif
	! accept/reject insertion move
    if (rosen+beta*(-e_new+mu)<-99.) rosen = -99.-beta*(-e_new+mu) ! avoids underflow
	if (real(size(space))/(nmolec+1)*exp(rosen+beta*(-e_new+mu)) > ran2(seed).and.success) then	!accept move
		if (nmolec == maxmolec) then ! cannot add another molecule
			write (*,*) 'Maximum number of molecules exceeded in transfer routine'
			write (24,*) 'Maximum number of molecules exceeded in transfer routine'
			stop
		endif
		e_total = e_total + e_new
        ! put new chain into last position
		do ibead = 1,nbead
			xyz(1:dim,ibead,nmolec+1) = pos_history(1:dim,ibead)
		enddo
		nmolec = nmolec + 1
	else  ! reject move
		do ibead = 1,bead_count ! remove partially or fully grown new chain
			space(pos_history(1,ibead),pos_history(2,ibead),pos_history(3,ibead)) = 0
		enddo
		success = .false.
	endif
endif

end

subroutine transfer2(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
                nmolec,nbead,bead_type,interact,e_self,xyz,success,seed,beta,e_total,mu)

! routine inserts a molecule in new random position or removes a randomly selected one with equal probabilities
! Thermal Rosenbluth factors!

implicit none
integer,    intent(in) ::        dim,xdim,ydim,zdim,maxmolec,nbead
integer*1,    intent(inout) ::    space(xdim,ydim,zdim)
integer,    intent(in) ::        length(dim)
integer,    intent(in) ::        nvectors
integer,    intent(in) ::        vector(dim,nvectors)
integer,    intent(in) ::        ngroups
integer,    intent(inout) ::    nmolec
integer,    intent(in) ::        bead_type(nbead)
real*8,        intent(in) ::        interact(0:ngroups,0:ngroups)
real*8,        intent(in) ::        e_self
integer,    intent(inout) ::    xyz(dim,nbead,maxmolec)
logical,    intent(out) ::        success
integer,    intent(inout) ::    seed
real*8,        intent(in) ::        beta
real*8,        intent(inout) ::    e_total
real*8,        intent(in) ::        mu

! local variables
real               ran2                 ! function
real*8               aux                  ! for picking direction to go to
integer            pos(dim),npos(dim)   ! local auxiliary variables
integer            i,ibead,idim,ipos,jpos    ! local counters
integer            removed              ! removed molecule
integer            pos_neighb(dim,nvectors) ! neighbor positions
integer            pos_nneigh(dim,nvectors) ! next-neighbor positions
integer            inext                ! counter for next-neighbor positions
integer            pos_history(dim,nbead)
real*8             ene(nvectors),ene_b(nbead)        ! energy of insertion of a bead, also along chosen vector
real*8             rosen                ! Actual Rosenbluth factor (unlike transfer1 where it was log)
real*8             w_vect(nvectors)   ! cummulative weight up to nvector for picking direction
real*8             w_max
integer            bead_count    ! counts how many beads were successfully inserted

success = .false.
if (ran2(seed)<0.5) then    ! remove chain

    if (nmolec>=1) then

        ! pick random chain
        removed = int(nmolec*ran2(seed)) + 1
        rosen = 1.
        pos(1:dim) = xyz(1:dim,nbead,removed)
        space(pos(1),pos(2),pos(3)) = 0 ! delete last bead
        call find_neighb(dim,length,nvectors,pos,pos_neighb)
        ene_b(nbead) = 0.
        do ipos = 1,nvectors
            ene_b(nbead) = ene_b(nbead) + interact(bead_type(nbead),&
                    space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
        enddo
        do ibead = nbead-1,1,-1
            pos(1:dim) = xyz(1:dim,ibead,removed)
            call find_neighb(dim,length,nvectors,pos,pos_neighb)
            ene_b(ibead) = 0.
            do ipos = 1,nvectors
                ene_b(ibead) = ene_b(ibead) + interact(bead_type(ibead),&
                        space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
            enddo
            do ipos = 1,nvectors
                npos = (/ pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos) /)
                ene(ipos) = 99./beta ! large energy, avoids overflowing the exp
                if (space(npos(1),npos(2),npos(3))==0) then
                    ene(ipos) = 0.
                    call find_neighb(dim,length,nvectors,npos,pos_nneigh)
                    do jpos = 1,nvectors
                        ene(ipos) = ene(ipos) + interact(bead_type(ibead+1),&
                                space(pos_nneigh(1,jpos),pos_nneigh(2,jpos),pos_nneigh(3,jpos)))
                    enddo
                endif
                w_vect(ipos) = exp(-beta*ene(ipos))
            enddo
            do ipos = 2,nvectors
                w_vect(ipos) = w_vect(ipos-1) + w_vect(ipos)
            enddo
            space(pos(1),pos(2),pos(3)) = 0 ! delete bead
            rosen = rosen*w_vect(nvectors)/nvectors
        enddo
        success = .true.

        ! accept/reject deletion move
        if (nmolec/real(size(space))/rosen*exp(beta*(ene_b(1) - e_self - mu)) > ran2(seed)) then
            !accept deletion move
            e_total = e_total - ene_b(1) + e_self
            do ibead = 2,nbead
                e_total = e_total - ene_b(ibead)
            enddo
            ! exchange last molecule w/ removed one
            xyz(1:dim,1:nbead,removed) = xyz(1:dim,1:nbead,nmolec)
            nmolec = nmolec - 1
        else  ! reject deletion move
            do ibead = 1,nbead  ! restore chain in its old position
                space(xyz(1,ibead,removed),xyz(2,ibead,removed),xyz(3,ibead,removed)) = bead_type(ibead)
            enddo
            success = .false.
        endif
    endif

else    ! insert chain
    
    bead_count = 0
    ! pick random position
    do idim = 1,dim
        pos(idim) = int(ran2(seed)*length(idim)) + 1
    enddo
    success = .true.
    rosen = 1.
    if (space(pos(1),pos(2),pos(3))/=0) then
        success = .false.
    else
        space(pos(1),pos(2),pos(3))= bead_type(1)  ! place first bead
        bead_count = 1 ! fisrt bead placed
        pos_history(1:dim,bead_count) = pos(1:dim)
        ene_b(1) = 0.
        call find_neighb(dim,length,nvectors,pos,pos_neighb)
        do ipos = 1,nvectors
            ene_b(1) = ene_b(1) + interact(bead_type(1),space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
        enddo
        ibead = 2
        do while (ibead <= nbead .and. success)
            call find_neighb(dim,length,nvectors,pos,pos_neighb)
            do ipos = 1,nvectors
                npos = (/ pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos) /)
                ene(ipos) = 99./beta ! large energy, avoids overflowing the exp
                if (space(npos(1),npos(2),npos(3))==0) then
                    ene(ipos) = 0.
                    call find_neighb(dim,length,nvectors,npos,pos_nneigh)
                    do jpos = 1,nvectors
                        ene(ipos) =  ene(ipos) + interact(bead_type(ibead),&
                                space(pos_nneigh(1,jpos),pos_nneigh(2,jpos),pos_nneigh(3,jpos)))
                    enddo
                endif
                w_vect(ipos) = exp(-beta*ene(ipos))
            enddo
            do ipos = 2,nvectors
                w_vect(ipos) = w_vect(ipos-1) + w_vect(ipos)
            enddo
            rosen = rosen*w_vect(nvectors)/nvectors
            aux=ran2(seed)*w_vect(nvectors)
            ipos = 1
            do while ( w_vect(ipos) < aux .and. ipos < nvectors)
                ipos = ipos + 1
            enddo
            do idim = 1,dim ! periodic boundary conditions
                pos(idim) = pos(idim) + vector(idim,ipos)
                if (pos(idim)>length(idim)) pos(idim) = 1
                if (pos(idim)==0) pos(idim) = length(idim)
            enddo
            if (space(pos(1),pos(2),pos(3))==0) then ! found empty spot
                space(pos(1),pos(2),pos(3))=bead_type(ibead)
                ene_b(ibead) = ene(ipos)
                bead_count = ibead
                pos_history(1:dim,bead_count) = pos(1:dim)
            else
                success = .false.  ! overlapped
            endif
            ibead = ibead + 1
        enddo
    endif
    ! accept/reject insertion move
    if (real(size(space))/(nmolec+1) * rosen * exp(beta*(-ene_b(1) + e_self + mu)) > ran2(seed).and.success) then    !accept move
        if (nmolec == maxmolec) then ! cannot add another molecule
            write (*,*) 'Maximum number of molecules exceeded in transfer routine'
            write (24,*) 'Maximum number of molecules exceeded in transfer routine'
            stop
        endif
        e_total = e_total + ene_b(1) - e_self             ! energy for first bead + e_self
        do ibead = 2, nbead
            e_total = e_total + ene_b(ibead)  ! adds energies taken into account in Rosenbluth factor
        enddo
        ! put new chain into last position
        do ibead = 1,nbead
            xyz(1:dim,ibead,nmolec+1) = pos_history(1:dim,ibead)
        enddo
        nmolec = nmolec + 1
    else  ! reject move
        do ibead = 1,bead_count ! remove partially or fully grown new chain
            space(pos_history(1,ibead),pos_history(2,ibead),pos_history(3,ibead)) = 0
        enddo
        success = .false.
    endif
endif

end

subroutine transfer3(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
				nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)

! routine removes a randomly selected molecule and inserts it in new position using athermal Rosenbluth growth

implicit none
integer,	intent(in) ::		dim,xdim,ydim,zdim,maxmolec
integer*1,	intent(inout) ::	space(xdim,ydim,zdim)
integer,	intent(in) ::		length(dim)
integer,	intent(in) ::		nvectors
integer,	intent(in) ::		vector(dim,nvectors)
integer,	intent(in) ::		ngroups
integer,	intent(in) ::	    nmolec,nbead
integer,	intent(in) ::		bead_type(nbead)
real*8,		intent(in) ::		interact(0:ngroups,0:ngroups)
integer,	intent(inout) ::	xyz(dim,nbead,maxmolec)
logical,	intent(out) ::		success
integer,	intent(inout) ::	seed
real*8,		intent(in) ::		beta
real*8,		intent(inout) ::	e_total
integer*8,    intent(in) ::       istep

! local variables
real			ran2		! function
integer			pos(dim)								! local auxiliary variable
integer			ibead,idim,ipos,ivector	! local counters
integer			removed							! removed molecule
integer			pos_neighb(dim,nvectors)
integer			pos_history(dim,nbead)
real*8			e_old, e_new	! energy of old and new configurations
real*8			rosen	! Logarithm of Rosenbluth factor; set to special value of -99/beta if zero
integer			n_unocc	! how many sites are unoccupied
integer			v_unocc(nvectors)	! unoccupied vectors
integer			bead_count	! counts how many beads were successfully inserted
integer			start_bead,end_bead,increment	! loop controls


if (ran2(seed)<0.5) then	! start from last bead
	start_bead = nbead
	end_bead = 1
	increment = -1
else	! start from first bead
	start_bead = 1
	end_bead = nbead
	increment = 1
endif

! pick random chain
removed = int(nmolec*ran2(seed)) + 1
rosen = 0.
e_old = 0.
e_new = 0.
success = .true.

do ibead = start_bead,end_bead,increment
	pos(1:dim) = xyz(1:dim,ibead,removed)
	n_unocc = 0
	call find_neighb(dim,length,nvectors,pos,pos_neighb)
	do ipos = 1,nvectors
		if (ibead/=start_bead .and. space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) &
			n_unocc = n_unocc + 1
		e_old = e_old + interact(space(pos(1),pos(2),pos(3)),&	
			space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
	enddo
	space(pos(1),pos(2),pos(3)) = 0 ! delete entity	
	if (ibead/=start_bead) rosen = rosen - log(real(n_unocc))
enddo

! pick random position
do idim = 1,dim
    pos(idim) = int(ran2(seed)*length(idim)) + 1
enddo

if (space(pos(1),pos(2),pos(3))/=0) then
    success = .false.
    bead_count = 0
else
    space(pos(1),pos(2),pos(3)) = bead_type(end_bead)
    pos_history(1:dim,end_bead) = pos(1:dim)
    bead_count = end_bead
    call find_neighb(dim,length,nvectors,pos,pos_neighb)
    do ipos = 1,nvectors
        e_new = e_new + interact(space(pos(1),pos(2),pos(3)),&
        space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
    enddo
    do ibead = end_bead-increment,start_bead,-increment
        if (success) then
            n_unocc = 0
            v_unocc(1) = 1
            call find_neighb(dim,length,nvectors,pos,pos_neighb)
            do ipos = 1,nvectors
                if (space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) then
                    n_unocc = n_unocc + 1
                    v_unocc(n_unocc) = ipos
                endif
            enddo
            bead_count = ibead
            if (n_unocc /= 0) then
                rosen = rosen + log(real(n_unocc))
            else
                success = .false.
                bead_count = ibead + increment
            endif
            ivector = v_unocc(int(ran2(seed)*n_unocc) + 1) ! select random unocc. position
            do idim = 1,dim
                pos(idim) = pos(idim) + vector(idim,ivector)
                if (pos(idim)>length(idim)) pos(idim) = 1
                if (pos(idim)==0) pos(idim) = length(idim)
            enddo
            if (space(pos(1),pos(2),pos(3))==0) then
                space(pos(1),pos(2),pos(3))=bead_type(ibead)
                pos_history(1:dim,bead_count) = pos(1:dim)
                call find_neighb(dim,length,nvectors,pos,pos_neighb)
                do ipos = 1,nvectors
                    e_new = e_new + interact(space(pos(1),pos(2),pos(3)),&
                    space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
                enddo
            endif
        endif
    enddo
endif

! accept/reject regrowth move
if (rosen-beta*(e_new-e_old)<-99.) rosen = -99.+beta*(e_new-e_old) ! avoids underflow
if (exp(+rosen-beta*(e_new-e_old)) > ran2(seed) .and. success) then
    !accept move
    e_total = e_total - e_old + e_new
    do ibead = start_bead,end_bead,increment
        xyz(1:dim,ibead,removed) = pos_history(1:dim,ibead)
    enddo
else  ! reject move
    success = .false.
endif
if (.not.success) then ! restore
    if (bead_count > 0) then
        do ibead = end_bead,bead_count,-increment ! remove partially or fully grown new chain
            space(pos_history(1,ibead),pos_history(2,ibead),pos_history(3,ibead)) = 0
        enddo
    endif
    do ibead = start_bead,end_bead,increment  ! restore chain in its old position
		space(xyz(1,ibead,removed),xyz(2,ibead,removed),xyz(3,ibead,removed)) = &
			bead_type(ibead)
	enddo
endif

end

subroutine flip(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,ngroups,&
                nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)

! routine flips sequence of a randomly selected molecule

implicit none
integer,    intent(in) ::       dim,xdim,ydim,zdim
integer*1,  intent(inout) ::    space(xdim,ydim,zdim)
integer,    intent(in) ::       length(dim)
integer,    intent(in) ::       nvectors
integer,    intent(in) ::       ngroups
integer,    intent(in) ::       maxmolec
integer,    intent(in) ::       nmolec,nbead
integer,    intent(in) ::       bead_type(nbead)
real*8,     intent(in) ::       interact(0:ngroups,0:ngroups)
integer,    intent(inout) ::    xyz(dim,nbead,maxmolec)
logical,    intent(out) ::      success
integer,    intent(inout) ::    seed
real*8,     intent(in) ::       beta
real*8,     intent(inout) ::    e_total
integer*8,    intent(in) ::       istep   !DDD for debugging

! local variables
real            ran2        ! function
integer         pos(dim)    ! local auxiliary variable
integer         ivector, ibead, jbead, removed, ix, iy, iz    ! local counters
integer         pos_neighb(dim,nvectors)
integer         pos_history(dim,nbead)
integer         pos_his_n(dim,nvectors,nbead)  ! to avoid calling find_neighb twice w/ same positions
real*8          e_old, e_new    ! energy of old and new configurations
logical         nn ! two beads interact


! pick random chain
removed = int(nmolec*ran2(seed)) + 1
e_old = 0.
e_new = 0.
do ibead = 1,nbead ! old interactions
    pos(1:dim) = xyz(1:dim,ibead,removed)
    pos_history(1:dim,ibead) = pos(1:dim)
    call find_neighb(dim,length,nvectors,pos,pos_neighb)
    do ivector = 1,nvectors
        pos_his_n(1:dim,ivector,nbead-ibead+1) = pos_neighb(1:dim,ivector)
        e_old = e_old + interact(space(pos(1),pos(2),pos(3)),&
            space(pos_neighb(1,ivector),pos_neighb(2,ivector),pos_neighb(3,ivector)))
    enddo
    space(pos(1),pos(2),pos(3)) = 0
enddo

do ibead = 1,nbead ! flip chain and compute new interactions
    pos(1:dim) = xyz(1:dim,ibead,removed)
    space(pos(1),pos(2),pos(3)) = bead_type(nbead-ibead+1)
    do ivector = 1,nvectors
        e_new = e_new + interact(space(pos(1),pos(2),pos(3)),&
            space(pos_his_n(1,ivector,nbead-ibead+1),pos_his_n(2,ivector,nbead-ibead+1),pos_his_n(3,ivector,nbead-ibead+1)))
    enddo
enddo

! accept/reject flip move
if (exp(beta*(e_old-e_new)) > ran2(seed)) then !accept move
    success = .true.
    e_total = e_total + e_new - e_old
    do ibead = 1,nbead
        xyz(1:dim,ibead,removed) = pos_history(1:dim,nbead-ibead+1)
    enddo
else  ! reject move, restoring space
    success = .false.
    do ibead = 1,nbead
        space(xyz(1,ibead,removed),xyz(2,ibead,removed),xyz(3,ibead,removed)) = bead_type(ibead)
    enddo
endif

end

subroutine regrowth(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,vector,ngroups,&
                nmolec,nbead,bead_type,interact,xyz,success,seed,beta,e_total,istep)

! routine removes half of a randomly selected molecule and inserts it using athermal Rosenbluth growth

implicit none
integer,    intent(in) ::        dim,xdim,ydim,zdim,maxmolec
integer*1,    intent(inout) ::    space(xdim,ydim,zdim)
integer,    intent(in) ::        length(dim)
integer,    intent(in) ::        nvectors
integer,    intent(in) ::        vector(dim,nvectors)
integer,    intent(in) ::        ngroups
integer,    intent(in) ::        nmolec,nbead
integer,    intent(in) ::        bead_type(nbead)
real*8,     intent(in) ::        interact(0:ngroups,0:ngroups)
integer,    intent(inout) ::    xyz(dim,nbead,maxmolec)
logical,    intent(out) ::        success
integer,    intent(inout) ::    seed
real*8,     intent(in) ::        beta
real*8,     intent(inout) ::    e_total
integer*8,    intent(in) ::       istep   !DDD for debugging

! local variables
real            ran2        ! function
integer            pos(dim)                                ! local auxiliary variable
integer            ibead,idim,ipos,ivector,ix,iy,iz    ! local counters
integer            removed                            ! removed molecule
integer            pos_neighb(dim,nvectors)
integer            pos_history(dim,nbead)
real*8            e_old, e_new    ! energy of old and new configurations
real*8            rosen    ! Logarithm of Rosenbluth factor; set to special value of -99/beta if zero
integer            n_unocc    ! how many sites are unoccupied
integer            v_unocc(nvectors)    ! unoccupied vectors
integer            bead_count    ! counts how many beads were successfully inserted
integer            start_bead,end_bead,increment    ! loop controls

end_bead = ran2(seed)*(nbead-3)+2 ! change introduced in v. 1.05, was fixed at nbead/2 previously
if (ran2(seed)<0.5) then    ! forward
    start_bead = nbead
    increment = -1
else    ! first half
    start_bead = 1
    increment = 1
endif
! pick random chain
removed = int(nmolec*ran2(seed)) + 1
rosen = 0.
e_old = 0.
e_new = 0.
success = .true.
do ibead = start_bead,end_bead,increment
    pos(1:dim) = xyz(1:dim,ibead,removed)
    n_unocc = 0
    call find_neighb(dim,length,nvectors,pos,pos_neighb)
    do ipos = 1,nvectors
        if (ibead/=nbead .and. ibead/=1 .and. space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) &
            n_unocc = n_unocc + 1
        e_old = e_old + interact(space(pos(1),pos(2),pos(3)),&
            space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
    enddo
    space(pos(1),pos(2),pos(3)) = 0 ! delete entity
    if (ibead/=nbead.and.ibead/=1) rosen = rosen - log(real(n_unocc))
enddo


! locate last non-removed bead
pos(1:dim) = xyz(1:dim,end_bead+increment,removed)
n_unocc = 0
call find_neighb(dim,length,nvectors,pos,pos_neighb)
do ipos = 1,nvectors
    if (space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) n_unocc = n_unocc + 1
enddo
rosen = rosen - log(real(n_unocc))

do ibead = end_bead,start_bead,-increment
    if (success) then
        n_unocc = 0
        v_unocc(1) = 1
        call find_neighb(dim,length,nvectors,pos,pos_neighb)
        do ipos = 1,nvectors
            if (space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos))==0) then
                n_unocc = n_unocc + 1
                v_unocc(n_unocc) = ipos
            endif
        enddo
        if (n_unocc /= 0) then
            rosen = rosen + log(real(n_unocc))
        else
            rosen = -99. / beta  ! small value for the exponential, but avoids underflow
        endif
        ivector = v_unocc(int(ran2(seed)*n_unocc) + 1) ! select random unocc. position
        do idim = 1,dim
            pos(idim) = pos(idim) + vector(idim,ivector)
            if (pos(idim)>length(idim)) pos(idim) = 1
            if (pos(idim)==0) pos(idim) = length(idim)
        enddo
        if (space(pos(1),pos(2),pos(3))==0) then
            space(pos(1),pos(2),pos(3))=bead_type(ibead)
            bead_count = ibead
            pos_history(1:dim,bead_count) = pos(1:dim)
            call find_neighb(dim,length,nvectors,pos,pos_neighb)
            do ipos = 1,nvectors
                e_new = e_new + interact(space(pos(1),pos(2),pos(3)),&
                space(pos_neighb(1,ipos),pos_neighb(2,ipos),pos_neighb(3,ipos)))
            enddo
        else
            success = .false.  ! overlapped
        endif
    endif
enddo

! accept/reject regrowth move
if (rosen-beta*(e_new-e_old)<-99.) rosen = -99.+beta*(e_new-e_old) ! avoids underflow
if (exp(+rosen-beta*(e_new-e_old)) > ran2(seed).and.success) then
    !accept move
    e_total = e_total - e_old + e_new
    do ibead = end_bead,start_bead,-increment
        xyz(1:dim,ibead,removed) = pos_history(1:dim,ibead)
    enddo
else  ! reject move
    do ibead = end_bead,bead_count,-increment ! remove partially or fully grown new chain
        space(pos_history(1,ibead),pos_history(2,ibead),pos_history(3,ibead)) = 0
    enddo
    do ibead = start_bead,end_bead,increment  ! restore chain in its old position
        space(xyz(1,ibead,removed),xyz(2,ibead,removed),xyz(3,ibead,removed)) = &
            bead_type(ibead)
    enddo
    success = .false.
endif

end



subroutine cluster_move(dim,xdim,ydim,zdim,maxmolec,length,nvectors,vector,nmolec,nbead,xyz,space,max_ncluster,success,seed,istep)

! routine picks up a random amphiphile molecule and moves the cluster to which it belongs by a fixed amount.

implicit none

integer,	intent(in) ::		dim,xdim,ydim,zdim,maxmolec
integer*1,	intent(inout) ::	space(xdim,ydim,zdim)
integer,	intent(in) ::		length(dim)
integer,	intent(in) ::		nvectors
integer,	intent(in) ::		vector(dim,nvectors)
integer,	intent(in) ::	    nmolec
integer,	intent(in) ::		nbead
integer,	intent(inout) ::	xyz(dim,nbead,maxmolec)
integer,	intent(in) ::		max_ncluster ! (BEADS!!!)
logical,	intent(out) ::		success
integer,	intent(inout) ::	seed
integer*8,    intent(in) ::       istep   ! used for debugging

real		ran2				! function
integer		ncluster 			! number of cluster sites (beads)
integer		imolec,ipos,ivector,idim,iprev,i,ibead	! local counters
integer		pos(dim)			! local auxiliary position variable
integer		pos_cluster(dim,max_ncluster)	! positions of cluster sites found
integer*1	save_cluster(max_ncluster)	! keeps type of site in cluster from space for recovery
integer		new_pos_cluster(dim,max_ncluster)	! new positions of cluster sites found
integer		newpos(dim,nvectors)	! local auxiliary variables
integer		displ_vector	! vector along which cluster is displaced
integer		rmolec			! randomly selected molecule

!pick random molecule
rmolec = int(ran2(seed)*nmolec) + 1

! find cluster in which it is located: connectivity defined by ALL sites
ncluster = 1
save_cluster(1) = space(xyz(1,1,rmolec),xyz(2,1,rmolec),xyz(3,1,rmolec))
pos_cluster(1:dim,1) = xyz(1:dim,1,rmolec)
! the next construct sets a position to vacuum so that it is not counted twice
space(pos_cluster(1,1),pos_cluster(2,1),pos_cluster(3,1)) = 0
ipos = 1
do while (ipos<=ncluster)
	! look around all newly found cluster sites
	pos(1:dim) = pos_cluster(1:dim,ipos)
	call find_neighb(dim,length,nvectors,pos,newpos)
	do ivector = 1, nvectors
		if ( space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector)) /= 0) then
			! found a new position belonging to the same cluster
			ncluster = ncluster + 1
			if (ncluster>max_ncluster) then
				write (*,'(A,I6)') "Max_ncluster exceeded in cluster_move routine, ncluster=", ncluster
				write (24,'(A,I6)') "Max_ncluster exceeded in cluster_move routine, ncluster=", ncluster
				stop
			endif
			save_cluster(ncluster) = &
				space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector))
			pos_cluster(1:dim,ncluster) = newpos(1:dim,ivector)
			space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector)) = 0
		endif
	enddo
	ipos = ipos + 1
enddo
success=.false.

if(ncluster<nmolec*nbead) then ! only do cluster move if system as a whole is not one cluster
    ! select random vector along which to displace cluster
    displ_vector = int(ran2(seed)*nvectors)+1
    success = .true.
    ipos = 1
    do while(ipos <= ncluster .and. success)
        new_pos_cluster(1:dim,ipos) = pos_cluster(1:dim,ipos) + vector(1:dim,displ_vector)
        do idim = 1,dim ! periodic boundary conditions
            if (new_pos_cluster(idim,ipos)>length(idim)) new_pos_cluster(idim,ipos) = 1
            if (new_pos_cluster(idim,ipos)==0) new_pos_cluster(idim,ipos) = length(idim)
        enddo
        pos(1:dim) = new_pos_cluster(1:dim,ipos)
        call find_neighb(dim,length,nvectors,pos,newpos)
        do ivector = 1,nvectors
            if ( space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector)) /= 0) success = .false.	!#3D
        enddo
        ipos = ipos + 1
    enddo
endif
if (success) then	! accept move
    do imolec = 1,nmolec
        if (space(xyz(1,1,imolec),xyz(2,1,imolec),xyz(3,1,imolec))==0) then
            ! found a chain that was part of the cluster, need to displace
            do ibead = 1,nbead
                xyz(1:dim,ibead,imolec) = xyz(1:dim,ibead,imolec) + vector(1:dim,displ_vector)
                do idim = 1,dim ! periodic boundary conditions
                    if (xyz(idim,ibead,imolec)>length(idim)) xyz(idim,ibead,imolec) = 1
                    if (xyz(idim,ibead,imolec)==0) xyz(idim,ibead,imolec) = length(idim)
                enddo
            enddo
        endif
    enddo
    do ipos = 1,ncluster
        space(new_pos_cluster(1,ipos),new_pos_cluster(2,ipos),new_pos_cluster(3,ipos)) = &
            save_cluster(ipos)
    enddo
else ! reject move
	do ipos = 1,ncluster
		space(pos_cluster(1,ipos),pos_cluster(2,ipos),pos_cluster(3,ipos)) = &
			save_cluster(ipos)
	enddo
endif

return 
end


subroutine cluster_count(dim,xdim,ydim,zdim,maxmolec,space,length,nvectors,nmolec,&
						nbead,xyz,ncluster,icluster,max_ncluster,isnap,n_snap,seed)

! This routine picks a random chain and returns the cluster size, ncluster,
! in which it belongs.  Cluster connectivity is defined as any two beads within interaction range.

implicit none

integer,	intent(in) ::		dim,xdim,ydim,zdim,maxmolec,nbead
integer*1,	intent(inout) ::	space(xdim,ydim,zdim)
integer,	intent(in) ::		length(dim)
integer,	intent(in) ::		nvectors
integer,	intent(in) ::	    nmolec
integer,	intent(inout) ::	xyz(dim,nbead,maxmolec)
integer,	intent(out) ::		ncluster    ! beads in cluster
integer,	intent(out) ::		icluster    ! counts how many different chains in the cluster
integer,	intent(in) ::		max_ncluster ! (BEADS!!!)
integer,	intent(inout) ::	isnap ! counts snapshots taken
integer,	intent(in) ::		n_snap ! determines how frequently snapshots are taken
integer,	intent(inout) ::	seed

real		ran2				    ! function
integer		i,imolec,ipos	    ! local counters
integer		pos(dim)			    ! local auxiliary position variable
integer		pos_cluster(dim,max_ncluster)	! positions of cluster sites found
integer*1	save_cluster(max_ncluster)	! keeps type of site in cluster from space for recovery
integer		newpos(dim,nvectors)	! local auxiliary variables
integer		ivector                 ! local variable
integer		ibead		            ! counter
logical*1	is_in_cluster(maxmolec)	
integer		tot_cluster				! total sites in cluster
integer		rmolec					! randomly picked molecule

!pick random molecule
rmolec = int(ran2(seed)*nmolec) + 1

! find cluster in which it is located: connectivity defined by ALL sites
tot_cluster = 1
save_cluster(1) = space(xyz(1,1,rmolec),xyz(2,1,rmolec),xyz(3,1,rmolec))
pos_cluster(1:dim,1) = xyz(1:dim,1,rmolec)
! Set position to vacuum so that it is not counted twice
space(pos_cluster(1,1),pos_cluster(2,1),pos_cluster(3,1)) = 0
ncluster = 1
ipos = 1
do while (ipos<=tot_cluster)
	! look around all newly found cluster sites
	pos(1:dim) = pos_cluster(1:dim,ipos)
	call find_neighb(dim,length,nvectors,pos,newpos)
	do ivector = 1, nvectors
		if ( space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector)) /= 0) then
			! found a new position belonging to the same cluster
			tot_cluster = tot_cluster + 1
			if (tot_cluster > max_ncluster) then
				write (*,*) "Max_ncluster exceeded in cluster_count routine",tot_cluster
				write (24,*) "Max_ncluster exceeded in cluster_count routine",tot_cluster
				stop
			endif
			save_cluster(tot_cluster) = space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector))
			pos_cluster(1:dim,tot_cluster) = newpos(1:dim,ivector)
			space(newpos(1,ivector),newpos(2,ivector),newpos(3,ivector)) = 0
			ncluster = ncluster + 1
		endif
	enddo
	ipos = ipos + 1
enddo	

icluster = ncluster/nbead
isnap = isnap+1

if (mod(isnap,n_snap)==0) then	! save snapshot of cluster to disk
	write (34,'((3i5))') icluster,nmolec,ncluster
	do imolec = 1,nmolec
		if (is_in_cluster(imolec)) then
			write (34,'(70i4)') (xyz(1:dim,ibead,imolec),ibead=1,nbead)
		endif
	enddo
endif

do ipos = 1,tot_cluster
		space(pos_cluster(1,ipos),pos_cluster(2,ipos),pos_cluster(3,ipos)) = save_cluster(ipos) ! reset space
enddo

return 
end


subroutine find_neighb(dim,length,nvectors,pos,pos_neighb)
! finds neighbors around pos, taking into account periodic boundary
! conditions; ugly, fast version; ONLY 3D

implicit none
integer, intent(in) :: dim,pos(dim),length(dim)
integer, intent(in) :: nvectors
integer, dimension(dim,nvectors), intent(out) :: pos_neighb
integer ivector ! local counter

if (nvectors==26) then
	do ivector = 1,9
		pos_neighb(1,ivector) = pos(1) - 1
		if (pos_neighb(1,ivector) == 0) pos_neighb(1,ivector) = length(1)
	enddo
	do ivector = 10,17
		pos_neighb(1,ivector) = pos(1) 
	enddo
	do ivector = 18,26
		pos_neighb(1,ivector) = pos(1) + 1
		if (pos_neighb(1,ivector) > length(1)) pos_neighb(1,ivector) = 1
	enddo
	do ivector = 1,3
		pos_neighb(2,ivector) = pos(2) - 1
		if (pos_neighb(2,ivector) == 0) pos_neighb(2,ivector) = length(2)
	enddo
	do ivector = 10,12
		pos_neighb(2,ivector) = pos(2) - 1
		if (pos_neighb(2,ivector) == 0) pos_neighb(2,ivector) = length(2)
	enddo
	do ivector = 18,20
		pos_neighb(2,ivector) = pos(2) - 1
		if (pos_neighb(2,ivector) == 0) pos_neighb(2,ivector) = length(2)
	enddo
	do ivector = 4,6
		pos_neighb(2,ivector) = pos(2) 
	enddo
	do ivector = 13,14
		pos_neighb(2,ivector) = pos(2) 
	enddo
	do ivector = 21,23
		pos_neighb(2,ivector) = pos(2) 
	enddo
	do ivector = 7,9
		pos_neighb(2,ivector) = pos(2) + 1
		if (pos_neighb(2,ivector) > length(2)) pos_neighb(2,ivector) = 1
	enddo
	do ivector = 15,17
		pos_neighb(2,ivector) = pos(2) + 1
		if (pos_neighb(2,ivector) > length(2)) pos_neighb(2,ivector) = 1
	enddo
	do ivector = 24,26
		pos_neighb(2,ivector) = pos(2) + 1
		if (pos_neighb(2,ivector) > length(2)) pos_neighb(2,ivector) = 1
	enddo
	do ivector = 1,13,3
		pos_neighb(3,ivector) = pos(3) - 1 
		if (pos_neighb(3,ivector) == 0) pos_neighb(3,ivector) = length(3)
	enddo
	do ivector = 15,24,3
		pos_neighb(3,ivector) = pos(3) - 1 
		if (pos_neighb(3,ivector) == 0) pos_neighb(3,ivector) = length(3)
	enddo
	do ivector = 2,11,3
		pos_neighb(3,ivector) = pos(3) 
	enddo
	do ivector = 16,25,3
		pos_neighb(3,ivector) = pos(3)
	enddo
	do ivector = 3,12,3
		pos_neighb(3,ivector) = pos(3) + 1 
		if (pos_neighb(3,ivector) > length(3)) pos_neighb(3,ivector) = 1
	enddo
	do ivector = 14,26,3
		pos_neighb(3,ivector) = pos(3) + 1 
		if (pos_neighb(3,ivector) > length(3)) pos_neighb(3,ivector) = 1
	enddo
else	! nvectors = 6
!	vector(1:dim,1)  = (/ 1, 0, 0/)	
!	vector(1:dim,2)  = (/ 0, 1, 0/)	
!	vector(1:dim,3)  = (/ 0, 0, 1/)	
!	vector(1:dim,4)  = (/ -1, 0, 0/)	
!	vector(1:dim,5)  = (/ 0, -1, 0/)	
!	vector(1:dim,6)  = (/ 0, 0, -1/)	
	do ivector = 1,6
		pos_neighb(1:dim,ivector) = pos(1:dim)
	enddo
	pos_neighb(1,1) = pos(1) + 1
	if (pos_neighb(1,1) > length(1)) pos_neighb(1,1) = 1
	pos_neighb(2,2) = pos(2) + 1
	if (pos_neighb(2,2) > length(2)) pos_neighb(2,2) = 1
	pos_neighb(3,3) = pos(3) + 1
	if (pos_neighb(3,3) > length(3)) pos_neighb(3,3) = 1
	pos_neighb(1,4) = pos(1) - 1
	if (pos_neighb(1,4) == 0) pos_neighb(1,4) = length(1)
	pos_neighb(2,5) = pos(2) - 1
	if (pos_neighb(2,5) == 0) pos_neighb(2,5) = length(2)
	pos_neighb(3,6) = pos(3) - 1
	if (pos_neighb(3,6) == 0) pos_neighb(3,6) = length(3)
endif

return
end		


subroutine stats (N,X,Xavg,Std)

! This routine calculates the average and standard deviation of a sample

IMPLICIT NONE
INTEGER				 :: N   ! Number of points in sample (input)
REAL*8, DIMENSION(N)   :: X	! vector containing sample values (input)
REAL*8                 :: Xavg ! mean value	(output)
REAL*8                 :: Std	! standard deviation (output)
INTEGER              :: I	! index (auxiliary variable)

Std = 0 ; Xavg = 0
DO I = 1, N
  Xavg = Xavg + X(I)
END DO

Xavg = Xavg / N

DO I = 1, N
  Std = Std + (X(I) - Xavg) ** 2
END DO

Std = SQRT( Std / (N - 1) )

return
end

FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,&
IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,&
NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
	idum=max(-idum,1)
	idum2=idum
	do j=NTAB+8,1,-1
		k=idum/IQ1
		idum=IA1*(idum-k*IQ1)-k*IR1
		if (idum.lt.0) idum=idum+IM1
		if (j.le.NTAB) iv(j)=idum
	enddo
	iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END

