program V2_RIS C C--- C Same version as published in W.L. Mattice and U.W. Suter, "Conformational C Theory of Large Molecules," John Wiley & Sons, 1994, Appendic C, pp 429-436. C References to equations are for the same book. C--- C Computes the exact average of a vector V that is the sum of C contributions V(i), rigidly connected to the skeletal bonds of polymer C chains of regular constitution and configuration. C C The program is capable of handling any number of bonds/unit, but not C more than NU rotational isomeric states/bond (NU is a PARAMETER). The C chain length increases step-wise by a factor of 2; the maximum length is C selected by the number of steps ("loops") LMAX [max. X = 2**(LMAX-1)]. C C Matrices are normalized in each step to prevent numerical overflow or C underflow and there is no strict upper bound to the chain length. C Roundoff error, however, can corrupt the results for very long chains; C adjustment of the quantity EPSILON (a parameter in subroutine NRMLZ) C allows for maximum (machine-dependent!) accuracy without underflow C errors (EPSILON = 10**(-24) is a reasonable setting for most computers). C--- C Attention: Line 9 of the subroutine VMATV, p 436 of the book, contains a C printing error: the last variable, declared as "R(*)", should read "VR(*)". C--- C (C) Copyright 1993 : Wayne L. Mattice & Ulrich W. Suter C--- integer I, II, IN, J, JJ, JN, L, LMAX, M, M1, M5, MAXL, MOLD, N, & N5, NBPU, NOLD, NU, NUN, NUNIT, & LNBLNK parameter (NU=6) character TITLE*80 real*8 PHI(NU), THETA(NU), U(NU,NU), UI(NU,NU), U1(NU), UX(NU), & VEC(3), VECN(3), VEC1(3), & G(5,5), GG(5*NU,5*NU), GGI(5*NU,5*NU), GG1(5*NU), & GGX(5*NU), TMP(5*NU,5*NU), & C, CINF, COLD, ENLSQ, FAC, FG, GMAX, LNZ, SUMVEC, UMAX, & VEC1N, VM2, XINV, Z, ZFAC data U1 / NU*0. / data UX / NU*1. / C--- C Format statements too long to fit into the write statements themselves: C 1 format ('1Characteristic Ratio of 0 for', & ' Chains with Regular Constitution',//, & 1x,' +++ ',a,' +++',///, & ' Number of bonds per unit =',i3,/, & ' Longest chain :',i14,' units') 2 format (//,' Bond',i3,' : ',i3,' x',i3,' U-matrix',10x, & 'associated vector =',3f8.4,/) 3 format (1h1,//,1x,' +++ ',a,' +++',///, & ' Dimension of U(unit) =',i3,' x',i3,/) 4 format (//,10x,'X',5x,'1/X',8x,'ln(Z)',14x,'',6x,'C(X)', & 7x,'est. C(inf)',/) C--- C Initialize matrices to unit-matrices for later convenience C U1(1) = 1. 50 do 70 I = 1, NU do 60 J=1,NU 60 U(I,J) = 0. 70 U(I,I) = 1. do 90 I=1,5*NU do 80 J=1,5*NU 80 GG(I,J) = 0. 90 GG(I,I) = 1. C C Read running TITLE for all pages C read (*,'(a)',end=900) TITLE C C Read # of bonds/unit & # of chain length loops [max.X = 2**(LMAX-1)]. C Defaults are: NBPU : none, LMAX : 10 C SUMVEC = 0. read (*,*,end=900) NBPU, LMAX if (LMAX.le.0) LMAX = 10 MAXL = 2**(LMAX-1) write (*,1) TITLE(1:LNBLNK(TITLE)), NBPU, MAXL C C Read U-matrices and vectors associated with the bonds, then C torison angles and bond angles for each unit. Running sum for V**2. C do 300 NUN=1,NBPU read (*,*,end=900) M, N, (VEC(I),I=1,3) SUMVEC = SUMVEC + VEC(1)**2 + VEC(2)**2 + VEC(3)**2 if (NUN.eq.1) then M1 = M MOLD = M NOLD = M endif write (*,2) NUN, M, N, (VEC(I),I=1,3) if (M.ne.NOLD) go to 800 if ((NUN.eq.NBPU).and.(N.ne.M1)) go to 800 if (NUN.eq.1) then VECN(1) = VEC(1) VECN(2) = VEC(2) VECN(3) = VEC(3) endif if (NUN.eq.NBPU) then VEC1(1) = VEC(1) VEC1(2) = VEC(2) VEC1(3) = VEC(3) endif do 130 I=1,M read (*,*,end=900) (UI(I,J),J=1,N) 130 write (*,'(1x,10f12.6)') (UI(I,J),J=1,N) write (*,'(/,a)') ' Torsion angles' read (*,*,end=900) (PHI(I),I=1,N) write (*,'(1x,10f12.1)') (PHI(I),I=1,N) write (*,'(/,a)') ' Bond angles' read (*,*,end=900) (THETA(I),I=1,N) C C - default for any THETA-angle beyond the first: take the previously C defined one (if THETA is the same for all PHI's of a bond, simply C add zero's to complete the number required by the input statement). C do 140 I=2,N 140 if (THETA(I).lt.1.) THETA(I) = THETA(I-1) write (*,'(1x,10f12.1)') (THETA(I),I=1,N) C C Now compute U(unit) and GG(unit) (=Script-G(unit)) as running product. C First U C call MATMAT (U,UI,U,NU,MOLD,M,N,TMP,5*NU) C C Now Script-G (according to equation (VI-51), with A replaced by G) C M5 = 5*(M-1) + 1 N5 = 5*(N-1) + 1 do 160 J=1,N5,5 JN = 1 + J/5 call GMAT (THETA(JN),PHI(JN),VEC,G) do 160 I=1,M5,5 IN = 1 + I/5 do 160 II=1,5 do 160 JJ=1,5 160 GGI(II+I-1,JJ+J-1) = G(II,JJ)*UI(IN,JN) call MATMAT (GG,GGI,GG,5*NU,5*MOLD,5*M,5*N,TMP,5*NU) C C Now compute Script-G for first and last bond (GG1 and GGX) C according to equations (VI-50) and (VI-52) with A replaced by G C (the first bond of the chain is taken as the last bond of the repeat unit, C occurring once before the repetition starts, the last bond as the first C bond of the repeat unit, taken once after the repetition ends - hence, C the chain-length will be equal to the [number of units]*NBPU + 2). C if (MOLD.lt.M) MOLD = M NOLD = N if (NUN.le.1) then do 170 I=1,M5,5 do 170 II=1,5 170 GGX(I+II-1) = G(II,5) endif if (NUN.ge.NBPU) then do 190 I=1,5 190 GG1(I) = G(1,I) do 200 I=6,5*NU 200 GG1(I) = 0. endif 300 continue C C Ready for looping over chain length. C First write summery of setup, then go ... C write (*,3) TITLE(1:LNBLNK(TITLE)), M1, N do 350 I=1,N 350 write (*,'(1x,10f12.6)') (U(I,J),J=1,N) write (*,'(/,1x,a,f10.4,/)') & 'Sum of mean-square length of vectors / unit =', SUMVEC write (*,4) NUNIT = 1 FAC = 0. ZFAC = 0. VEC1N = VEC1(1)**2 + VEC1(2)**2 + VEC1(3)**2 & + VECN(1)**2 + VECN(2)**2 + VECN(3)**2 C C Now loop C do 500 L=1,LMAX C C - calculate Z, Z, , V[av]**2 (V[av] is the rms average C of V(i) over the entire chain), and C(X)=/(N*V[av]**2). C call VMATV (U1,U,UX,NU,N,Z,TMP) call VMATV (GG1,GG,GGX,5*NU,5*N,FG,TMP) LNZ = ZFAC + log(Z) VM2 = exp(FAC)*FG/Z ENLSQ = real(NUNIT)*SUMVEC + VEC1N C = VM2/ENLSQ XINV = 1./real(NUNIT) C C - extrapolate and write results. C if (L.ge.2) then CINF = 2*C - COLD write (*,'(i11,2e12.4,f17.2,f11.4,f15.4)') & NUNIT, XINV, LNZ, VM2, C, CINF else write (*,'(i11,2e12.4,f17.2,f11.4)') & NUNIT, XINV, LNZ, VM2, C endif if (L.lt.LMAX) then COLD = C C - prevent overflow and underflow by conditioning U and GG. call NRMLZ (U,NU,N,UMAX) call NRMLZ (GG,5*NU,5*N,GMAX) FAC = FAC + log(GMAX) - log(UMAX) ZFAC = ZFAC + log(UMAX) C - square U and GG for double chain length. call MATMAT (U,U,U,NU,N,N,N,TMP,5*NU) call MATMAT (GG,GG,GG,5*NU,5*N,5*N,5*N,TMP,5*NU) NUNIT = 2*NUNIT FAC = 2*FAC ZFAC = 2*ZFAC endif 500 continue C and restart. go to 50 C That's it. 800 write (*,'(/,a)') ' ***** WRONG DIMENSIONS *****' 900 write (*,'(///,1x,a)') '***** End of run *****' end C subroutine GMAT (THETA,PHI,V,G) C--- C Sets up the generator matrix G for V**2, given scalar values of C THETA(i), PHI(i), and the vector V(i), according to equation (VI-21). C--- real*8 G(5,*), V(*), T(3,3), & PHI, THETA integer I, J C--- do 100 I=1,5 do 100 J=1,5 100 G(I,J) = 0. C C Transformation matrix T according to equation (VI-5). C call TMAT (THETA,PHI,T) C do 110 I=2,4 do 110 J=2,4 110 G(I,J) = T(I-1,J-1) G(1,1) = 1. G(1,2) = 2.*(V(1)*T(1,1) + V(2)*T(2,1) + V(3)*T(3,1)) G(1,3) = 2.*(V(1)*T(1,2) + V(2)*T(2,2) + V(3)*T(3,2)) G(1,4) = 2.*(V(1)*T(1,3) + V(2)*T(2,3) + V(3)*T(3,3)) G(1,5) = V(1)**2 + V(2)**2 + V(3)**2 G(2,5) = V(1) G(3,5) = V(2) G(4,5) = V(3) G(5,5) = 1. return end C integer function LNBLNK (STRING) C--- C Returns the location of the last non-blank character in string STRING. C--- character STRING*(*) integer I C--- do 100 I=len(STRING),1,-1 100 if (STRING(I:I).ne.' ') go to 200 I = 0 200 LNBLNK = I return end C subroutine MATMAT (A,B,AB,MM,MA,MB,NB,TMP,MTMP) C--- C Multiplies matrix A (MAxMB) with matrix B (MBxNB) and stores C the result in matrix AB (MAxNB). AB can be one of A or B. C The first dimension of the arrays containing A, B, and AB in the C calling program's defining statement is MM. C The first dimension of the array containing TMP in the C calling program's defining statement is MTMP. C--- real*8 A(MM,*), B(MM,*), AB(MM,*), TMP(MTMP,*), & X integer I, J, K, MA, MB, MM, NB C--- do 100 I=1,MA do 100 J=1,NB X = 0. do 90 K=1,MB 90 X = X + A(I,K)*B(K,J) 100 TMP(I,J) = X C do 200 I=1,MA do 200 J=1,NB 200 AB(I,J) = TMP(I,J) return end C subroutine NRMLZ (A,MM,N,FAC) C--- C Divides every element of square matrix A (NxN) by that with the largest C magnitude, FAC. Set elements that are smaller than EPSILON in magnitude C to true zero (EPSILON is a PARAMETER). C The first dimension of the array containing A in the C calling program's defining statement is MM. C--- real*8 A(MM,*), & EPSILON, FAC integer I, J, MM, N parameter (EPSILON=1.0d-24) C--- FAC = 0. do 100 I=1,N do 100 J=1,N 100 FAC = max(FAC,abs(A(I,J))) C do 200 I=1,N do 200 J=1,N A(I,J) = A(I,J)/FAC 200 if (abs(A(I,J)).lt.EPSILON) A(I,J) = 0. return end C subroutine TMAT (THETA,PHI,T) C--- C Sets up the transformation matrix T according to equation (VI-4). C--- real*8 T(3,*), & COSTHT, COSPHI, PHI, SINTHT, SINPHI, THETA integer I, J C--- COSTHT = cos(THETA*3.1415926535/180.) SINTHT = sin(THETA*3.1415926535/180.) COSPHI = cos(PHI*3.1415926535/180.) SINPHI = sin(PHI*3.1415926535/180.) T(1,1) = - COSTHT T(1,2) = SINTHT T(1,3) = 0. T(2,1) = - SINTHT * COSPHI T(2,2) = - COSTHT * COSPHI T(2,3) = - SINPHI T(3,1) = - SINTHT * SINPHI T(3,2) = - COSTHT * SINPHI T(3,3) = COSPHI C return end C subroutine VMATV (VL,A,VR,MM,M,RESULT,TMP) C--- C Forms the bylinear product from the vector VL, C the square matrix A (MxM), and the vector VR: C RESULT = VL(transpose) A VR C The first dimension of the array containing A in the C calling program's defining statement is MM. C--- real*8 A(MM,*), TMP(*), VL(*), VR(*), & RESULT, X integer I, K, MM, M C--- do 100 I=1,M X = 0. do 90 K=1,M 90 X = X + VL(K)*A(K,I) 100 TMP(I) = X X = 0. C do 200 I=1,M 200 X = X + TMP(I)*VR(I) RESULT = X return end