c***********************************************************
c multihydro_x.for begins here
c***********************************************************

c Program MULTIHYDRO 
c Version 1

c XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX	
c X														   X			
c X   J. Garcia de la Torre, A. Ortega, H.E. Perez Sanchez   X
c X              and J.G Hernandez Cifre                     X
c X "MULTIHYDRO and MONTEHYDRO: Conformational search and    X
c X  Monte Carlo calculation of solution properties of rigid X
c X                flexible bead models                      X
c X          Biophys. Chem. 116, 121-128 (2005)              X
c X                                                          X
c XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX


c***********************************************************
c General section #1

	PROGRAM MULTIHYDRO

	PARAMETER(nmax=20000)
	DIMENSION x(nmax,3),e(nmax),u(3)
	CHARACTER title*20, filename*5, filename2*13

c seed for random numbers to be used by subroutine RAN1
c must be smaller than 50000
	iseed=50000
        
	iconf=0

c***********************************************************
c Specific section #2 (compulsory data)

	temp=20.    
	icase=12
	eta0=0.010   
	rm=110000.   
	vbar=0.710
	solden=1.0   
	title='The bent dumb-bell'
	filename='confi'
	nconf=36

c***********************************************************
c Specific section #3 (Model data - dimension is Angs.)

c This is for the broken dumbbell, with two arms, one 
c with n1=2+1 (1 is the hinge), and n2 beads, , respectively.
c with radius r1 Angs
c At the ends of the arms there are two larger beads of 
c radius r2 and r3 

	n1=3 
	n2=2 
	ntotal=n1+n2+2
	r1=5. 
	r2=25. 
	r3=25.


c***********************************************************
c General section #4

c ntotconf will be the total number of conformations
c generated (eventually discounting from nconf the ones
c that could be discarded for some reason)

	ntotconf= 0

	OPEN(1,FILE='hydro.dat')

	DO ic=1,nconf


c***********************************************************
c Specific section #5 Construction of each conformation

	i=0
c First arm, aligned with (negative) axis Z
	DO i1=1,n1 
	  i=i+1 
	  e(i)=r1 
	  x(i,3)=-(i1-1)*2.*r1 
	  x(i1,1)=0. 
	  x(i1,2)=0. ; 
	ENDDO

c One of the big spheres
	i=i+1
	e(i)=r2  
	x(i,3)=x(i-1,3)-r2-r1 
	x(i,1)=0. 
	x(i,2)=0. ;

c Second arm, making angle theta

	theta=(3.1416/nconf)*(ic-1)
	u(1)=sin(theta) 
	u(2)=0. 
	u(3)=cos(theta) 
	DO i2=1,n2
	  i=i+1 
	  e(i)=r1 
	  rad=i2*2.*r1
	  x(i,1)=rad*u(1) 
	  x(i,2)=rad*u(2) 
	  x(i,3)=rad*u(3) 
	ENDDO

c Second big sphere 
	i=i+1 
	e(i)=r3 
	rad=rad+r3+r1
	x(i,1)=rad*u(1) 
	x(i,2)=rad*u(2) 
	x(i,3)=rad*u(3) 

	print *,'ic,theta',ic,theta

	n=i


c Check overlap between spheres

	do i=1,n
	  call CHECK_OVERLAP(x,e,n,nmax,i,iflag)
	  if(iflag.ne.0) print *,'Found overlap between beads:',i,iflag,  
     &  ' for theta=',theta
	  if(iflag.ne.0) GOTO 222
	enddo


C***********************************************************
c General section #6 (down to the end of the code)

c NOTE: At this stage, the following values must be given 
c n : number of beads
c x(i,k), k=1,3, i=1,n  bead coordinates
c e(i) bead radii

c Write files

	  iconf=iconf+1

	  i1=iconf/100
	  i2=(iconf-i1*100)/10
	  i3=mod(iconf,10)
	  WRITE(1,101) title
101	  FORMAT(a20)
	  WRITE(1,100) filename,i1,i2,i3
	  WRITE(1,100) filename,i1,i2,i3,'.txt'
100	  format(a5,3i1,a4)
	  BACKSPACE 1
	  READ(1,*) filename2
        WRITE(1,250) icase
	  WRITE(1,150) temp
	  WRITE(1,160) eta0
	  WRITE(1,150) rm
	  WRITE(1,160) vbar
	  WRITE(1,150) solden
	  WRITE(1,250) 0
	  WRITE(1,250) 0
	  WRITE(1,250) 0
	  WRITE(1,250) 1
150	  FORMAT(F10.2,' ,')
160	  FORMAT(F6.3)
250     FORMAT(I4,' ,')

	  PRINT *,filename2
	  OPEN(2,FILE=filename2)

	  WRITE(2,*) '1.E-8'
	  WRITE(2,300) n
300	  FORMAT(i5)
	  DO i=1,n
	    WRITE(2,*) (x(i,k),k=1,3),e(i)
	  ENDDO
	  CLOSE(2)

222	CONTINUE
	ENDDO

	print *,'Number of valid conformations generated:',iconf
	WRITE(1,301) '*' ; CLOSE(1)
301	FORMAT(a1)
	END

c Utility subroutines
c --------------------


	SUBROUTINE CHECK_OVERLAP(x,e,n,nmax,i,iflag)

c      Check overlap of bead i with all the other beads
c      iflag=0 if there is no overlap
c      Otherwise, iflag equals the index of the first
c         bead, found by the program, with which bead i overlaps

	DIMENSION x(nmax,3),e(nmax)
	iflag=0
	DO ii=1,n
	  IF(ii.NE.i) THEN
	    dist2=(x(i,1)-x(ii,1))**2+(x(i,2)-x(ii,2))**2
     &      +(x(i,3)-x(ii,3))**2
	    dist=sqrt(dist2)
		distmin=e(i)+e(ii)
		IF(dist.LT.distmin*0.9999) THEN
		  iflag=1
		  RETURN
		ENDIF
	  ENDIF
	ENDDO


	END

c ----------------------------------------------

c This is a random numbers generation subroutine
c taken from the "Numerical Recipes" Book
c Instead of this, you may prefer to use your compiler's RAN 
c subroutine, or your favourite one, 
 
C Subroutine for random number generation. 
C IDUM is a dummy integer number (initial seed
c    smaller than 50000)

      FUNCTION RAN1(IDUM)
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
      DATA IFF /0/
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
          IX1=MOD(IA1*IX1+IC1,M1)
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11      CONTINUE
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
      IF(J.GT.97.OR.J.LT.1)PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END

c------------------
	SUBROUTINE GENER_FJEV(INDEV,ISEED,NMAX,N,A,B,X,E)		

c Generate freely jointed chain with excluded volume
 
c Data:
c ------
c ISEED seed for random numbers
c INDEV =0 no excluded volume
c      not 0 excluded volume hard spheres
c N: number of beads (N-1 segments)
c NMAX : dimension of X and E in calling subroutine
c B: length of segments
c A: diameter of beads
c Results
c ------
c X(i,k) Coordinate k of bead i
c E(i) Radius of bead i

c Conformations in which there is a pair of beads
c with centers separated less than A are rejected.
 
	IMPLICIT NONE

	INTEGER NMAX,N,ISEED,INDEV
	REAL X(NMAX,3),E(NMAX)
	REAL A,B
	REAL SINTHE,COSTHE,PHI,RAN1
	INTEGER I,IFLAG

	I=1
	X(I,1)=0.
	X(I,2)=0.
	X(I,3)=0.
	E(I)=A/2

100	CONTINUE
	DO I=2,N
	  COSTHE=-1.+2.*RAN1(ISEED)
        SINTHE=SQRT(1.-COSTHE**2)
	  PHI=2.*3.1416*RAN1(ISEED)
	  X(I,1)=X(I-1,1) + B*COS(PHI)*SINTHE
	  X(I,2)=X(I-1,2) + B*SIN(PHI)*SINTHE
	  X(I,3)=X(I-1,3) + B*COSTHE
	  E(I)=A/2
	  IF(INDEV.NE.0) THEN
	    CALL CHECK_OVERLAP(X,E,I,NMAX,I,IFLAG)
	    IF(IFLAG.NE.0) THEN
CX	      PRINT *,'OVERLAPPING BEADS:',I,IFLAG
	      GO TO 100
	    ENDIF
	  ENDIF
	ENDDO

	RETURN
	END
	  
c***********************************************************
c multyhydro_x.for ends here
c***********************************************************
