Program HCOOR C************************************************************************ C* January 1998 * C* Program HCOOR * C* * C* by Mario NARDELLI * C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica, * della Universita` degli Studi di C* Chimica Fisica della Universita' degli Studi di Parma * C* Centro di Studio per la Strutturistica Diffrattometrica del CNR, * C* Viale delle Scienze 78, I-43100 Parma (Italy) * C* * C* This routine calculates the fractional coordinates of atom X lying * C* on the line joining atom A to atom B at a requested distance D from * C* atom A and nearest to atom B. * C* * C* INPUT * C* - Title (80A1) * C* - a,b,c,alpha,beta,gamma (free format) * C* - 'ATOM(A)',x/a,y/b,z/c (free format; ATOM(A):A4) * C* - 'ATOM(B)',x/a,y/b,z/c (free format; ATOM(B):A4) * C* - 'ATOM(X)',D (free format; ATOM(X):A4) * C* * C* OUTPUT * C* Title * C* Table of fractional coordinates * C* Distance ATOM(A)...ATOM(X) * C* * C************************************************************************ CHARACTER*1 TIT(80) CHARACTER*4 ATM(3) CHARACTER*14 FILIN,FILOUT DOUBLE PRECISION O(3,3),OI(3,3) DIMENSION X(3,3),C(3),A(3),XO(3,3),SN(3),CS(3),DDC(3),X1(3),X2(3) RD=57.295779513 WRITE(*,100) 100 FORMAT(1X,'Key in the name of the input file'/) READ(*,'(A14)')FILIN WRITE(*,101) 101 FORMAT(1X,'Key in the name of the output file'/) READ(*,'(A14)')FILOUT OPEN(UNIT=13,FILE=FILIN,FORM='FORMATTED') REWIND(UNIT=13) OPEN(UNIT=15,FILE=FILOUT,FORM='FORMATTED') REWIND(UNIT=15) READ(13,102)TIT 102 FORMAT(80A1) WRITE(15,103)TIT 103 FORMAT(1X,80A1/) READ(13,*)C,A READ(13,*)(ATM(I),(X(I,K),K=1,3),I=1,2) READ(13,*)ATM(3),D13 DO 1 I=1,3 AA=A(I)/RD SN(I)=SIN(AA) 1 CS(I)=COS(AA) CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3)) SNAR=SQRT(1.0-CSAR**2) O(1,1)=C(1) O(1,2)=C(2)*CS(3) O(1,3)=C(3)*CS(2) O(2,1)=0.0 O(2,2)=C(2)*SN(3) O(2,3)=-C(3)*SN(2)*CSAR O(3,1)=0.0 O(3,2)=0.0 O(3,3)=C(3)*SN(2)*SNAR DO 2 I=1,2 DO 2 J=1,3 XO(I,J)=0.0 DO 2 K=1,3 2 XO(I,J)=XO(I,J)+O(J,K)*X(I,K) DO 3 K=1,3 DDC(K)=XO(2,K)-XO(1,K) 3 D12=D12+DDC(K)**2 D12=SQRT(D12) DO 4 K=1,3 DDC(K)=DDC(K)/D12 X1(K)=XO(1,K)+DDC(K)*D13 4 X2(K)=XO(1,K)-DDC(K)*D13 DH1=0.0 DH2=0.0 DO 5 K=1,3 DH1=DH1+(XO(2,K)-X1(K))**2 5 DH2=DH2+(XO(2,K)-X2(K))**2 DH1=SQRT(DH1) DH2=SQRT(DH2) DO 7 K=1,3 XO(3,K)=X1(K) 7 IF(DH1.GT.DH2)XO(3,K)=X2(K) DH=0.0 DO 8 K=1,3 8 DH=DH+(XO(1,K)-XO(3,K))**2 DH=SQRT(DH) CALL INVMAT3(O,OI) DO 6 J=1,3 X(3,J)=0.0 DO 6 K=1,3 6 X(3,J)=X(3,J)+OI(J,K)*XO(3,K) WRITE(15,104)(ATM(I),(X(I,K),K=1,3),I=1,3),ATM(1),ATM(3),DH 104 FORMAT(1X,'Atom',9X,'x',9X,'y',9X,'z'/3(1X,A4,3X,3F10.5/)/ 1 ' Distance ',A4,'-',A4,' = ',F5.3) STOP END SUBROUTINE INVMAT3(A,B) C-----Inverte una matrice 3x3 con il metodo di partizione. DOUBLE PRECISION A(3,3),B(3,3),D,AM1,AM2,AM3,AM4,C,X D=A(1,1)*A(2,2)-A(1,2)*A(2,1) IF(DABS(D).LT.1D-10) GO TO 1 AM1=A(1,2)*A(2,3)-A(1,3)*A(2,2) AM2=A(1,3)*A(2,1)-A(1,1)*A(2,3) AM3=A(2,1)*A(3,2)-A(2,2)*A(3,1) AM4=A(1,2)*A(3,1)-A(1,1)*A(3,2) X=A(3,3)+(AM1*A(3,1)+AM2*A(3,2))/D B(3,3)=1D0/X X=X*D B(1,3)=AM1/X B(1,1)=(B(1,3)*AM3+A(2,2))/D B(1,2)=(B(1,3)*AM4-A(1,2))/D B(2,3)=AM2/X B(2,1)=(B(2,3)*AM3-A(2,1))/D B(2,2)=(B(2,3)*AM4+A(1,1))/D B(3,1)=AM3/X B(3,2)=AM4/X GO TO 2 1 WRITE(15,100) 100 FORMAT(' The matrix cannot be inverted') STOP 2 RETURN END