c*********************************************************** c multihydro_x.for begins here c*********************************************************** c Program MULTIHYDRO c Version 2 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=49999 iconf=0 c*********************************************************** c Specific section #2 (compulsory data) icase=.... temp= ... eta0= ... rm= ... vbar=... solden=... title=' ... ' filename='...' nconf= ... c*********************************************************** c Specific section #3 (Model data - dimension is Angs.) c ... INSERT YOUR CODE HERE 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 C ... INSERT YOUR CODE HERE 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***********************************************************