program subttt implicit none integer*4 nat1,nat2,n,iargc,n2 character*80 f1,f2,f3,s character*1 pm character*28 sd real*8 si real*8,allocatable::L1(:,:,:),A1(:,:,:,:),G1(:,:,:),V1(:,:,:), 1 L2(:,:,:),A2(:,:,:,:),G2(:,:,:),V2(:,:,:), 1 L3(:,:,:),A3(:,:,:,:),G3(:,:,:),V3(:,:,:) logical lcompl1,lcompl2 if(iargc().lt.3)then write(6,222) 222 format(' Usage: subttt <+->') stop endif if(iargc().lt.3)stop call getarg(1,f1) call getarg(2,f2) call getarg(3,f3) si=-1.0d0 if(iargc().gt.3)then call getarg(4,pm) if(pm.eq.'+')si=1.0d0 endif if(si.lt.0.0d0)then sd=' Difference of two FILE.TTTs' else sd=' Sum of two FILE.TTTs ' endif write(6,200)sd open(11,file=f1,status='old') read(11,*) read(11,*)nat1 open(12,file=f2,status='old') read(12,*) read(12,*)nat2 if(nat1.ne.nat2)then write(6,*)nat1,nat2 call report(' Different atom numbers') endif n=3*nat1 n2=2*n allocate(L1(n2,3,3),A1(n2,3,3,3),G1(n2,3,3),V1(n2,3,3)) allocate(L2(n2,3,3),A2(n2,3,3,3),G2(n2,3,3),V2(n2,3,3)) allocate(L3(n2,3,3),A3(n2,3,3,3),G3(n2,3,3),V3(n2,3,3)) rewind 11 rewind 12 CALL READTEN(11,n,L1,G1,A1,V1,lcompl1,f1) CALL READTEN(12,n,L2,G2,A2,V2,lcompl2,f2) close(12) L3=L1+si*L2 call subtrg(n2,lcompl1,lcompl2,si,G1,G2,G3) A3=A1+si*A2 V3=V1+si*V2 open(22,file=f3) write(22,200)sd 200 format(a28) rewind 11 read(11,*) read(11,80)s write(22,80)s 80 format(a80) read(11,80)s write(22,80)s read(11,80)s write(22,80)s close(11) CALL WRITETTTrest(N,L3,A3,G3,V3,lcompl1.or.lcompl2,f3) end subroutine subtrg(n2,lcompl1,lcompl2,si,G1,G2,G3) implicit none integer*4 n2,n,l,i,j real*8 G1(n2,3,3),G2(n2,3,3),G3(n2,3,3),si logical lcompl1,lcompl2 n=n2/2 if(lcompl1)then if(lcompl2)then G3=G1+si*G2 else DO 1 l=1,n DO 1 i=1,3 DO 1 j=1,3 G3(l ,i,j)=G1(l ,i,j) c second contains only the imaginary part: 1 G3(l+n,i,j)=G1(l+n,i,j)+si*G2(l,i,j) endif else if(lcompl2)then DO 2 l=1,n DO 2 i=1,3 DO 2 j=1,3 G3(l ,i,j)= si*G2(l ,i,j) c first contains only the imaginary part: 2 G3(l+n,i,j)=G1(l,i,j)+si*G2(l+n,i,j) else G3=G1+si*G2 endif endif return end C ========================================== SUBROUTINE WRITETTTrest(N,ALPHA,A,G,V,lcompl,f) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(2*N,3,3),G(2*N,3,3),A(2*N,3,3,3),V(2*N,3,3) logical lcompl character*80 f NAT=N/3 do 5 j=1,len(f) 5 if(f(j:j).ne.' ')write(6,55)f(j:j) 55 format(a1,$) if(lcompl)then write(6,56)'complex.' 56 format(' written as ',a8) ic=3 else ic=0 write(6,56)'real. ' endif c DO 1 I=1,3 WRITE(22,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 II=3*(L-1)+IX 1 WRITE(22,2001)L,IX,(ALPHA(II,I,J),J=1,3),(ALPHA(II+N,I,J),J=1,ic) 2001 FORMAT(I5,1H ,I1,6g15.7) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(22,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 II=3*(L-1)+IX c last inner index magnetic 2 WRITE(22,2001)L,IX,(G(II,I,J),J=1,3),(G(II+N,I,J),J=1,ic) WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(22,2006)I,J 2006 FORMAT(' A(',I1,'(u),',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 II=3*(L-1)+IX 3 WRITE(22,2001)L,IX,(A(II,I,J,K),K=1,3),(A(II+N,I,J,K),K=1,ic) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' Atom/x vx vy vz') DO 11 I=1,3 WRITE(22,2012)I 2012 FORMAT(' AlphaV(',I1,',J):') DO 11 L=1,NAT DO 11 IX=1,3 II=3*(L-1)+IX 11 WRITE(22,2001)L,IX,(V(II,I,J),J=1,3),(V(II+N,I,J),J=1,ic) close(22) RETURN END C ================================================================ SUBROUTINE READTEN(io,N,P,G,A,V,lcompl,f) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(2*N,3,3),G(2*N,3,3),A(2*N,3,3,3),V(2*N,3,3) logical lcompl character*80 f character*160 s160 READ(io,*) READ(io,*) READ(io,*) READ(io,*) lcompl=.false. ic=0 DO 1 I=1,3 READ(io,*) c Test if complex: if(I.eq.1)then READ(io,160)s160 160 format(a160) nd=0 do 6 j=1,len(s160) 6 if(s160(j:j).eq.'.')nd=nd+1 lcompl=nd.gt.3 do 5 j=1,len(f) 5 if(f(j:j).ne.' ')write(6,55)f(j:j) 55 format(a1,$) if(lcompl)then write(6,56)'Complex' 56 format(': ',a7) ic=3 else write(6,56)'Real ' endif backspace io endif DO 1 L=1,N/3 DO 1 IX=1,3 M=3*(L-1)+IX 1 READ(io,*)(P(M,I,J),J=1,2),(P(M,I,J),J=1,3),(P(M+N,I,J),J=1,ic) READ(io,*) READ(io,*) DO 2 I=1,3 READ(io,*) DO 2 L=1,N/3 DO 2 IX=1,3 M=3*(L-1)+IX c last index magnetic, inner 2 READ(io,*)(G(M,I,J),J=1,2),(G(M,I,J),J=1,3),(G(M+N,I,J),J=1,ic) READ(io,*) READ(io,*) DO 3 I=1,3 DO 3 J=1,3 READ(io,*) DO 3 L=1,N/3 DO 3 IX=1,3 M=3*(L-1)+IX 3 READ(io,*)(A(M,I,J,K),K=1,2),(A(M,I,J,K),K=1,3), 1(A(M+N,I,J,K),K=1,ic) READ(io,*) READ(io,*) DO 4 I=1,3 READ(io,*) DO 4 L=1,N/3 DO 4 IX=1,3 M=3*(L-1)+IX 4 READ(io,*)(V(M,I,J),J=1,2),(V(M,I,J),J=1,3),(V(M+N,I,J),J=1,ic) RETURN END subroutine report(s) character*(*) s write(6,*)s stop end