PROGRAM OPROPAG implicit none integer*4 nat,i,iargc,ia,ic,IERR,ii,ix,j,ja,jj,jx,n,n1,nn,iab real*8 mc(3,3),v1(3),v2(3),v(3),t1,t2,t,mi(3,3),e(3),pi,ri(3), 1k(3,3),vn,ra(3),rr(3),rl(3),sp,a(3,3),a2(3,3),cc(3),twist,flim, 1ft,s2,s3,ct,st,ff(3,3) character*80 name character*8 s8 integer*4 ,allocatable::ik(:),j33(:) real*8,allocatable::r(:,:),f(:,:),fbig(:,:),a33(:),PS(:,:,:), 1AS(:,:,:),PB(:,:,:),ABT(:,:,:),rb(:,:) c pi=4.0d0*datan(1.0d0) WRITE(6,3000) 3000 FORMAT(/, 1' Takes a AAA trimer, propagates it to (A)n n-mer',/, 1' conserving A-A twist and interactions',/,/, 1' Input: CELL.X(.FC,,.TEN,.TTT)',/,/ 3' Output: FILE.X(.FC,.TTT,.TEN,ACRYSTAL.FC)',/,/) if(iargc().ne.1)call report('Usage: opropag ') call getarg(1,s8) read(s8,*)n open(2,file='CELL.X') read(2,2000)name 2000 format(a80) read(2,*)nat allocate(ik(nat),r(nat,3)) do 1 i=1,nat 1 read(2,*)ik(i),(r(i,j),j=1,3) close(2) if(mod(nat,3).ne.0)call report('Nat not divisible by three') n1=nat/3 c centers of each monomer mc=0.0d0 do 2 i=1,3 do 2 ix=1,3 do 3 ia=n1*(i-1)+1,n1*(i-1)+n1 3 mc(i,ix)=mc(i,ix)+r(ia,ix) 2 mc(i,ix)=mc(i,ix)/dble(n1) c vectors A ---- v1 ----> A ---- v2 ----> A do 4 ix=1,3 v1(ix)=mc(2,ix)-mc(1,ix) v2(ix)=mc(3,ix)-mc(2,ix) 4 v(ix)=0.5d0*(v1(ix)+v2(ix)) vn=dsqrt(v(1)**2+v(2)**2+v(3)**2) write(6,601)v1,v2,v 601 format(' V1:',3f12.6,/,' v2:',3f12.6,/,' ave:',3f12.6,/) c 2d moments of inertia of each monomer do 5 i=1,3 write(6,602)i 602 format(' Monomer',i2,':') mi=0.0d0 do 6 ia=n1*(i-1)+1,n1*(i-1)+n1 do 61 ii=1,3 61 ra(ii)=r(ia,ii)-mc(i,ii) sp=ra(1)*v(1)+ra(2)*v(2)+ra(3)*v(3) do 62 ii=1,3 62 ra(ii)=ra(ii)-sp*v(ii)/vn**2 do 6 ix=1,3 do 6 jx=1,3 6 mi(ix,jx)=mi(ix,jx)+ra(ix)*ra(jx) write(6,603)mi 603 format(' moment:',3f10.4,/,2(16x,3f10.4,/)) call TRED12(3,mi,a,e,2,IERR) c does a3 go along v? s2=a(1,2)*v(1)+a(2,2)*v(2)+a(3,2)*v(3) s3=a(1,3)*v(1)+a(2,3)*v(2)+a(3,3)*v(3) if(dabs(s3).lt.dabs(s2))then call sw(s2,s3) call sw(e(2),e(3)) do 67 ix=1,3 67 call sw(a(ix,2),a(ix,3)) s3=a(1,3)*v(1)+a(2,3)*v(2)+a(3,3)*v(3) endif if(s3.lt.0.0d0)then do 68 ix=1,3 68 a(ix,3)=-a(ix,3) endif write(6,604)e,a 604 format(' eigenvalues:',3f10.4,/, 1 ' eigenvectors:',3f10.4,/,2(16x,3f10.4,/)) if(i.eq.2)a2=a do 5 ii=1,3 5 k(i,ii)=a(ii,1) c twist: t1=acos(k(1,1)*k(2,1)+k(1,2)*k(2,2)+k(1,3)*k(2,3)) t2=acos(k(2,1)*k(3,1)+k(2,2)*k(3,2)+k(2,3)*k(3,3)) t=(t1+t2)/2.0d0 write(6,605)t1*180.0d0/pi,t2*180.0d0/pi,t*180.0d0/pi 605 format(' t1 t2 ave:',3f10.4,' deg',/) c propagate medium geometry: allocate(rb(n1*n,3)) open(2,file='FILE.X') write(2,2000)name write(2,*)n1*n twist=-2.0d0*t do 23 ix=1,3 23 cc(ix)=mc(2,ix)-2.0d0*v(ix) iab=0 do 7 ii=1,n twist=twist+t ct=cos(twist) st=sin(twist) c mass center of this unit: do 22 ix=1,3 22 cc(ix)=cc(ix)+v(ix) do 7 ia=n1+1,2*n1 iab=iab+1 c coordinates to the mass center: do 63 ix=1,3 63 ra(ix)=r(ia,ix)-mc(2,ix) c coordinates in moment of inertia system: do 64 ix=1,3 64 ri(ix)=ra(1)*a2(1,ix)+ra(2)*a2(2,ix)+ra(3)*a2(3,ix) c rotate by twist: rr(1)=ri(1)*ct-ri(2)*st rr(2)=ri(1)*st+ri(2)*ct rr(3)=ri(3) c coordinates in lab system: do 65 ix=1,3 65 rl(ix)=rr(1)*a2(ix,1)+rr(2)*a2(ix,2)+rr(3)*a2(ix,3) c add mass center do 66 ix=1,3 rl(ix)=rl(ix)+cc(ix) 66 rb(iab,ix)=rl(ix) 7 write(2,20)ik(ia),(rl(ix),ix=1,3) 20 format(i2,3f15.8,' 0 0 0 0 0 0 0 0.0') close(2) write(6,*)' FILE.X was written' allocate(PS(nat,3,3),AS(nat,3,3),PB(n*n1,3,3),ABT(n*n1,3,3)) CALL READTEN(PS,AS,.false.,r,NAT) ia=0 twist=-2.0d0*t do 12 ic=0,n-1 twist=twist+t ct=cos(twist) st=sin(twist) do 12 i=n1+1,2*n1 ia=ia+1 do 121 ix=1,3 do 121 jx=1,3 121 ff(ix,jx)=PS(i,ix,jx) c ff in moment of inertia system: call t22(ff,a2) c rotate by twist: call r2(ct,st,ff) c ff to lab system: call ti(ff,a2) do 123 ix=1,3 do 123 jx=1,3 123 PB(ia,ix,jx)=ff(ix,jx) do 122 ix=1,3 do 122 jx=1,3 122 ff(ix,jx)=AS(i,ix,jx) call t22(ff,a2) call r2(ct,st,ff) call ti(ff,a2) do 12 ix=1,3 do 12 jx=1,3 12 ABT(ia,ix,jx)=ff(ix,jx) call WRITETEN(PB,ABT,rb,n*n1) allocate(f( 3*nat, 3*nat)) open(20,file='CELL.FC') call READFF(3*nat,3*nat,f) close(20) allocate(FBIG(n*3*n1 ,n*3*n1 )) fbig=0.0d0 twist=-2.0d0*t do 8 ic=0,n-1 twist=twist+t ct=cos(twist) st=sin(twist) c atoms in AAA trimer: do 8 i=1,nat c atoms in polymer: ia=(ic-1)*n1+i if(ia.ge.1.and.ia.le.n*n1)then do 9 j=1,nat ja=(ic-1)*n1+j if(ja.ge.1.and.ja.le.n*n1)then do 74 ix=1,3 do 74 jx=1,3 74 ff(ix,jx)=f(ix+3*(i-1),jx+3*(j-1)) c ff in moment of inertia system: call t22(ff,a2) c rotate by twist: call r2(ct,st,ff) c ff to lab system: call ti(ff,a2) do 10 ix=1,3 ii=ix+3*(ia-1) do 10 jx=1,3 jj=jx+3*(ja-1) 10 fbig(ii,jj)=ff(ix,jx) endif 9 continue endif 8 continue open(20,file='FILE.FC') call WRITEFF(3*n1*n,3*n*n1,fbig) close(20) write(6,*)' FILE.FC written' allocate(a33(3*n1*n),j33(3*n1*n)) open(33,file='ACRYSTAL.FC',form='unformatted') flim=2.0d-7 do 11 i=1,3*n1*n nn=0 do 81 j=i,3*n1*n ft=fbig(i,j) if(abs(ft).gt.flim)then nn=nn+1 j33(nn)=j if(i.eq.j)ft=ft/2.0d0 a33(nn)=ft endif 81 continue 11 write(33)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) close(33) write(6,*)' ACRYSTAL.FC written' end SUBROUTINE WRITETEN(P,A,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3),X(NAT,3) real*8,allocatable::R(:,:) allocate(R(3,NAT)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(L,I)/BOHR DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted back to laboratory system' do L=1,NAT do J=1,3 do I=1,3 if(dabs(P(L,J,I)).gt.1.0d10)then write(6,*)'p',L,J,I,P(L,J,I) endif if(dabs(A(L,J,I)).gt.1.0d10)then write(6,*)'p',L,J,I,A(L,J,I) endif enddo enddo enddo OPEN(15,FILE='FILE.TEN') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(15,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F20.6,I5) DO 220 L=1,NAT DO 220 J=1,3 220 WRITE(15,1501) (A(L,J,I),I=1,3),L Z0=0.0d0 DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(15,1501) (P(L,J,I),I=1,3),L CLOSE(15) WRITE(6,*)' FILE.TEN written ...' RETURN END SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' RETURN END SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END subroutine report(s) character*(*) s write(6,*)s stop end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP MACHEP = 1.0d-12 C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END subroutine sw(a,b) real*8 a,b,t t=a a=b b=t return end subroutine t22(ff,a2) implicit none integer*4 ix,jx real*8 ff(3,3),a2(3,3),ft(3,3) do 68 ix=1,3 do 68 jx=1,3 68 ft(ix,jx)=a2(1,ix)*ff(1,jx)+a2(2,ix)*ff(2,jx)+a2(3,ix)*ff(3,jx) do 69 ix=1,3 do 69 jx=1,3 69 ff(ix,jx)=a2(1,jx)*ft(ix,1)+a2(2,jx)*ft(ix,2)+a2(3,jx)*ft(ix,3) return end subroutine r2(ct,st,ff) implicit none integer*4 ix,jx real*8 st,ct,ff(3,3),ft(3,3) do 72 jx=1,3 ft(1,jx)=ff(1,jx)*ct-ff(2,jx)*st ft(2,jx)=ff(1,jx)*st+ff(2,jx)*ct 72 ft(3,jx)=ff(3,jx) do 73 ix=1,3 ff(ix,1)=ft(ix,1)*ct-ft(ix,2)*st ff(ix,2)=ft(ix,1)*st+ft(ix,2)*ct 73 ff(ix,3)=ft(ix,3) return end subroutine ti(ff,a2) implicit none integer*4 ix,jx real*8 ff(3,3),a2(3,3),ft(3,3) do 70 ix=1,3 do 70 jx=1,3 70 ft(ix,jx)=a2(ix,1)*ff(1,jx)+a2(ix,2)*ff(2,jx)+a2(ix,3)*ff(3,jx) do 71 ix=1,3 do 71 jx=1,3 71 ff(ix,jx)=a2(jx,1)*ft(ix,1)+a2(jx,2)*ft(ix,2)+a2(jx,3)*ft(ix,3) return end SUBROUTINE READTEN(P,A,LAPT,XX,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3),XX(NAT,3) LOGICAL LAPT real*8,allocatable::R(:,:) allocate(R(3,NAT)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=XX(L,I)/BOHR open(15,file='CELL.TEN') READ (15,*) C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2203 A(L,J,I)=SUM WRITE(6,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) 1 close(15) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted to atomic origin' RETURN END FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END