program SYMMOL implicit double precision (a-h,o-z) c logical esiste PARAMETER (NMA=1000) CHARACTER*2 ATOM(NMA) CHARACTER*80 RIGA c CHARACTER*60 inpfile common/matcel/PC(7),PCR(7),O(3,3),OI(3,3),G(3,3),GI(3,3),CS(12) COMMON/AT0/NA,NMO COMMON/AT1/X(3,NMA),AMAS(NMA),RCOV(NMA),MSP(NMA) COMMON/AT2/SX(3,NMA),SIG(NMA),DXM(NMA),DCM,DCME,indwgh,indtol COMMON/AT3/MOL(NMA),MLG(NMA) CHARACTER*6 NAME,SPEC COMMON/AT4/NAME(NMA) integer out COMMON/output/out c________________________________________________________________________ c if you want address all the output on the standard output, suppress the two c sequent instruction and activate the third c out=9 OPEN(9,file='symmol.out') c out=6 c________________________________________________________________________ if(out.eq.9)then write(*,*)' ============' write(*,*)' SYMMOL' write(*,*)' A PROGRAM FOR THE SYMMETRIZATION OF GROUPS OF ATOMS' write(*,*)' By Tullio Pilati and Alessandra Forni' write(*,*)' Version November 4th 2002' write(*,*)' ===================================================' write(*,*) endif write(out,*)' ============' write(out,*)' SYMMOL' write(out,*)' A PROGRAM FOR THE SYMMETRIZATION OF GROUPS OF ATOMS' write(out,*)' By Tullio Pilati and Alessandra Forni' write(out,*)' Version November 4th 2002' write(out,*)' ===================================================' write(out,*) c 900 write(*,*)'input file name?' c read(*,'(a)')inpfile c inquire(file=inpfile,exist=esiste) c if(.not.esiste)then c write(*,'(a60,a)')inpfile,' does not exist' c stop c else c open(1,file=inpfile) c endif 1000 read(*,'(A80)')RIGA if(RIGA(1:1).eq.'#') go to 1000 READ(RIGA,*)(PC(i),i=1,6) c******************************************** 1010 read(*,'(A80)')RIGA if(RIGA(1:1).eq.'#') go to 1010 READ(RIGA,*)indwgh,indtol,DCM,DCME if(DCME.lt.DCM)DCME=DCM write(out,*) if(indwgh.le.1)then indwgh=1 write(out,*)' INDWGH=1 ===> WEIGHTS AS ATOMIC MASS' endif if(indwgh.eq.2)write(out,*)' INDWGH=2 ===> UNITARY WEIGHTS' if(indwgh.eq.3)write(out,*)' INDWGH=3 ===> WEIGHTS BASED ON S.U.' if(indtol.le.1)then indtol=1 write(out,*)' INDTOL=1 ===> TOLERANCE=CONSTANT' endif if(indtol.eq.2)write(out,*) * ' INDTOL=2 ===> TOLERANCE BASED ON DISTANCES' if(indtol.eq.3) * write(out,*)' INDTOL=3 ===> TOLERANCE BASED ON S.U.' write(out,1)DCM,DCME 1 format(' CONSTANTS OF TOLERANCE=',2f7.3) write(out,*) c******************************************** i=1 MOmax=0 1100 read(*,'(A80)',END=1200)RIGA if(i.gt.NMA)then write(out,*)'ERROR: TOO MANY ATOMS IN INPUT' write(out,*)'Enlarge PARAMETER NMA' write(out,*)'Actual value: PARAMETER (NMA=1000)' stop endif c write(out,*)riga if(RIGA(1:1).eq.'#')go to 1100 READ(RIGA,'(A6,I2,6f9.5)')NAME(i),MOL(i),(X(k,i),k=1,3), *(SX(k,i),k=1,3) if(MOL(i).gt.MOmax)Momax=MOL(i) ATOM(i)=' ' SPEC=NAME(i) call COMPATTA(SPEC,6,K) K=1 if(SPEC(2:2).ge.'A'.and.SPEC(2:2).le.'Z')K=2 if(SPEC(2:2).ge.'a'.and.SPEC(2:2).le.'z')K=2 call MAIUSCOL(SPEC(1:1)) IF(K .GT. 1) then call MAIUSCOL(SPEC(2:2)) ATOM(i)=SPEC(1:2) else ATOM(i)(2:2)=SPEC(1:1) endif i=i+1 go to 1100 1200 NA=i-1 if(indwgh.eq.3.or.indtol.eq.3) then sxmax=0.000001d0 do i=1,NA do k=1,3 if(sx(k,i).gt.sxmax) sxmax=sx(k,i) enddo enddo if(sxmax.lt.0.00001d0) then write(out,*) write(out,*) ' ERROR: S.U. ON COORDINATES ARE MISSING.' if(indtol.eq.3) * write(out,*) ' INTRODUCE THE S.U. OR CHANGE INDTOL' if(indwgh.eq.3) * write(out,*) ' INTRODUCE THE S.U. OR CHANGE INDWGH' stop endif endif C call mass(NA,ATOM) c call costanti call cella write(out,*) write(out,*)' ATOM GROUP INPUT COORDINATES AND THEIR S.U.' write(out,*) do I=1,NA write(out,'(1x,A6,i2,1x,6f9.5,2x,i2)')NAME(I),MOL(I) *,(X(k,I),k=1,3),(SX(k,I),k=1,3) enddo C do 1300 MO=1,MOmax call work(MO) 1300 continue end subroutine arccos(cost,ang) C arco coseno in gradi implicit double precision (a-h,o-z) RAD=57.29577951308232D0 if(cost.gt.1.d0)cost=1.d0 if(cost.lt.-1.d0)cost=-1.d0 ang=RAD*DACOS(cost) return end subroutine asymunit(MK,IASU,N,NMS) PARAMETER (NMA=1000) PARAMETER (NMG=120) DIMENSION MK(NMA,NMG),IASU(NMA) do i =1,N IASU(I)=2 enddo DO 2800 I=1,N if(IASU(I).ne.2)go to 2800 DO 2790 J=1,NMS K=IABS(MK(I,J)) if(K.eq.I)go to 2790 if(K.eq.0)then IASU(I)=0 go to 2800 endif IASU(K)=1 2790 CONTINUE 2800 CONTINUE return end subroutine ax_order(A,i,m,msign,invers) C C m = group order for the matrix SIM(i) C msign= 1 asse di rotazione propria C msign=-1 asse di rotazione impropria C PARAMETER (maxorder=8) implicit double precision (a-h,o-z) dimension A(3,3,*),B(3,3),C(3,3,2) integer out COMMON/output/out call azzera(C,0.d0,18) C(1,1,1)=1.d0 C(2,2,1)=1.d0 C(3,3,1)=1.d0 C(1,1,2)=-1.d0 C(2,2,2)=-1.d0 C(3,3,2)=-1.d0 invers=2 msign=NINT(det(A,i)) call prodmm(A,C,B,i,1,1) do m=1,2*maxorder if(ium(C,B,1.d-2,2,1).eq.1)invers=1 if(ium(C,B,1.d-2,1,1).eq.1)go to 1000 call prodmm(A,B,B,i,1,1) enddo write(out,*)'INPUT PARAMETER DCM probably too HIGH. Reduce it!' stop 1000 if(m.lt.6.or.msign.eq.1)return m1=(m/4)*4 if(m1.ne.m)m=m/invers c write(out,*)'ax_order,i,m,msign,invers' c write(out,'(10i5)')i,m,msign,invers c write(out,'(3f10.6)')((A(ii,kk,i),kk=1,3),ii=1,3) return end subroutine azzera(A,COST,N) implicit double precision (a-h,o-z) C azzeraMENTO DIMENSION A(1) DO 1 I=1,N 1 A(I)=COST RETURN END subroutine bond(XO,IASU,N) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) integer out COMMON/output/out COMMON/AT1/X(3,NMA),AMAS(NMA),RCOV(NMA),MSP(NMA) CHARACTER*6 NAME COMMON/AT4/NAME(NMA) COMMON/AT5/LEGAMI(20,NMA),NLEG(NMA) DIMENSION XO(3,NMA),IASU(NMA) DIMENSION T(3) do 1200 i=1,N r1=RCOV(i) N1=0 do 1100 k=1,N if(IASU(k).eq.0)go to 1100 if(k.eq.i)go to 1100 r2=(r1+RCOV(k))*1.2 call combli(XO,XO,T,1.d0,-1.d0,i,k,1) call prods(t,t,d1,1,1,2) c write(out,'(a20,2i3,2f10.5)')'i,k,d1,r2',i,k,d1,r2 if(d1.gt.r2)go to 1100 N1=N1+1 if(N1.gt.20)then go to 1200 endif LEGAMI(N1,I)=k 1100 continue NLEG(i)=N1 1200 continue c write(out,*)'LEGAMI' c do ii=1,N c N1=NLEG(II) c write(out,'(20i3)')(LEGAMI(kk,ii),kk=1,N1) c enddo return end double precision function brent(ax,bx,cx,func,tol,xmin) implicit double precision (a-h,o-z) INTEGER ITMAX c REAL brent,ax,bx,cx,tol,xmin,func,CGOLD,ZEPS EXTERNAL func PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0e-10) INTEGER iter c REAL a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0.d0 fx=func(x) fv=fx fw=fx do 11 iter=1,ITMAX xm=0.5d0*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.d0*tol1 if(abs(x-xm).le.(tol2-.5d0*(b-a))) goto 3 if(abs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.d0) p=-p q=abs(q) etemp=e e=d if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) *goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(abs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=func(u) write(6,*)'brent,fu,fx',fu,fx if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue pause 'brent exceed maximum iterations' 3 xmin=x brent=fx return END subroutine cella implicit double precision (a-h,o-z) integer out COMMON/output/out common/matcel/PC(7),PCR(7),O(3,3),OI(3,3),G(3,3),GI(3,3),CS(12) RAD=57.29577951308232D0 ARAD=1.D0/RAD PC(7)=1.D0 PCR(7)=1.D0 COM=1.D0 DO 1030 I=1,3 K=I+3 CS(I)=COS(PC(K)*ARAD) CS(K)=SIN(PC(K)*ARAD) IF(PC(K).GT.0.D0)GO TO 1030 PC(K)=90.D0 CS(I)=0.D0 CS(K)=1.D0 1030 COM=COM-CS(I)*CS(I) COM=SQRT(COM+2.D0*CS(1)*CS(2)*CS(3)) O(1,1)=PC(1)*CS(6) O(1,2)=0.D0 O(1,3)=PC(3)*(CS(2)-CS(1)*CS(3))/CS(6) O(2,1)=PC(1)*CS(3) O(2,2)=PC(2) O(2,3)=PC(3)*CS(1) O(3,1)=0.D0 O(3,2)=0.D0 O(3,3)=PC(3)*COM/CS(6) IF(ABS(O(2,1)).LT.1.D-10)O(2,1)=0.D0 IF(ABS(O(2,3)).LT.1.D-10)O(2,3)=0.D0 IF(ABS(O(1,3)).LT.1.D-10)O(1,3)=0.D0 PC(7)=det(O,1) call traspo(O,G,1,1) call prodmm(G,O,G,1,1,1) call inv33(O,OI,1,1) PCR(7)=det(OI,1) call inv33(G,GI,1,1) PCR(1)=SQRT(GI(1,1)) PCR(2)=SQRT(GI(2,2)) PCR(3)=SQRT(GI(3,3)) CS(7)=GI(2,3)/(PCR(2)*PCR(3)) CS(8)=GI(1,3)/(PCR(1)*PCR(3)) CS(9)=GI(1,2)/(PCR(1)*PCR(2)) PCR(4)=RAD*DACOS(CS(7)) PCR(5)=RAD*DACOS(CS(8)) PCR(6)=RAD*DACOS(CS(9)) CS(10)=SIN(ARAD*PCR(4)) CS(11)=SIN(ARAD*PCR(5)) CS(12)=SIN(ARAD*PCR(6)) write(out,*) write(out,*)' CELL' write(out,'(3f10.5,3f10.3,f11.5,/)')PC c write(out,*)' RECIPROCAL CELL' c write(out,'(3f10.5,3f10.3,f11.5,/)')PCR c write(out,*)' METRIC MATRIX' c write(out,'(3g16.6)')((G(i,k),k=1,3),i=1,3) c write(out,*) c write(out,*)' ORTHOGONALIZATION MATRIX' c write(out,'(3g16.6)')((O(i,k),k=1,3),i=1,3) c write(out,*) return END subroutine combli(A,B,C,D,E,I,J,K) COMBINAZIONE LINEARE D*A+B*E=C implicit double precision (a-h,o-z) DIMENSION A(3,1),B(3,1),C(3,1) C(1,K)=A(1,I)*D+B(1,J)*E C(2,K)=A(2,I)*D+B(2,J)*E C(3,K)=A(3,I)*D+B(3,J)*E RETURN END subroutine compatta(WORD,L,K) implicit double precision (a-h,o-z) CHARACTER*(*) WORD K=0 DO 1 I=1,L IF(WORD(I:I).LE.' '.OR.WORD(I:I).GT.'~')GO TO 1 K=K+1 WORD(K:K)=WORD(I:I) 1 CONTINUE N=K+1 DO 2 I=N,L 2 WORD(I:I)=' ' RETURN END subroutine compl_gr(MK,N,*) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) PARAMETER (NMG=120) COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) DIMENSION CO(3,3),MK(NMA,NMG) complete the group 2520 NN=NMS DO 2600 I=1,NN DO 2600 J=1,NN call prodmm(SIM,SIM,CO,I,J,1) DO 2560 JJ=1,NN L=ium(CO,SIM,1.d-2,1,JJ) IF(L.EQ.1)GO TO 2590 2560 CONTINUE GO TO 2610 C group multiplication table 2590 MTG(I,J)=JJ 2600 CONTINUE c Gruppo completato RETURN C HA TROVATO UNA NUOVA MATRICE 2610 NMS=NMS+1 IF(NMS.LE.NMG) GO TO 2630 write(6,2) 2 format(' ERROR: TOO MANY MATRICES FOUND') write(6,'(3(3F10.5,/),/)')(((SIM(I,J,K),J=1,3),I=1,3),K=1,NMS) c ritorno per errore RETURN 1 2630 call trasfm(CO,SIM,1,NMS) do 2640 k=1,N k1=MK(k,J) MK(k,NMS)=MK(k1,I) 2640 continue go to 2520 end double precision function crms(t) implicit double precision (a-h,o-z) C C ppu=vettore XO ripetuto NMS volte C ppo=vettore XS trasformato da MK(NMA,NMS) C t(3)=vettore variazione angolare in radianti C PARAMETER (NMA=1000) PARAMETER (NMG=120) PARAMETER (NMV=NMA*NMG) common/rmsmin/RV(3,3),ppo(3,NMV),ppu(3,NMV),npu dimension t(3),po(3),pu(3),rpu(3) dimension RX(3,3),RY(3,3),RZ(3,3),RR(3,3) c write(6,'(6F9.5)')t call azzera(RX,0.d0,9) call azzera(RY,0.d0,9) call azzera(RZ,0.d0,9) call azzera(RR,0.d0,9) RX(1,1)=1.d0 RX(2,2)=dcos(t(1)) RX(2,3)=-dsin(t(1)) RX(3,2)=-RX(2,3) RX(3,3)=dcos(t(1)) RY(2,2)=1.d0 RY(1,1)=dcos(t(2)) RY(1,3)=-dsin(t(2)) RY(3,1)=-RY(1,3) RY(3,3)=dcos(t(2)) RZ(1,1)=dcos(t(3)) RZ(1,2)=-dsin(t(3)) RZ(2,1)=-RZ(1,2) RZ(2,2)=dcos(t(3)) RZ(3,3)=1.d0 c write(6,'(3g12.5)')RX c write(6,*) c write(6,'(3g12.5)')RY c write(6,*) c write(6,'(3g12.5)')RZ c write(6,*) call prodmm(RY,RX,RR,1,1,1) call prodmm(RZ,RR,RV,1,1,1) C ortonormalizzazione di precisione call prodv(RV,RV,RV,1,2,3) call prodv(RV,RV,RV,3,1,2) call prodv(RV,RV,RV,2,3,1) call norm(RV,1) call norm(RV,2) call norm(RV,3) c write(6,'(3g12.5)')RV C fine ortonormalizzazione di precisione func=0.d0 do i=1,npu po(1)=ppo(1,i) po(2)=ppo(2,i) po(3)=ppo(3,i) pu(1)=ppu(1,i) pu(2)=ppu(2,i) pu(3)=ppu(3,i) call prodmv(RV,pu,pu,1,1,1) c write(6,'(6F9.5)')po,pu func=func+(po(1)-pu(1))**2+(po(2)-pu(2))**2+(po(3)-pu(3))**2 enddo crms=func c write(6,*)func return end FUNCTION det(X,N) implicit double precision (a-h,o-z) C DETERMINANTE DELLA MATRICE X DIMENSION X(3,3,1) DET= +X(1,1,N)*X(2,2,N)*X(3,3,N)-X(1,1,N)*X(2,3,N)*X(3,2,N) DET=DET+X(1,2,N)*X(2,3,N)*X(3,1,N)-X(1,2,N)*X(2,1,N)*X(3,3,N) DET=DET+X(1,3,N)*X(2,1,N)*X(3,2,N)-X(1,3,N)*X(2,2,N)*X(3,1,N) RETURN END subroutine eigen(A,VEC,EIG,W,GAM,BET,BSQ,P,Q,IPO,IORD,IVP,NN) implicit double precision (a-h,o-z) C C -------------- C QCPE VERSION C DECEMBER 1971 C -------------- C C MATRIX DIAGNOLIZATION ROUTINE FOR REAL SYMMETRIC CASE C HOUSEHOLDER METHOD C RHO=UPPER LIMIT FOR OFF-DIAGONAL ELEMENT C NN SIZE OF MATRIX C A=MATRIX (ONLY LOWER TRIANGLE IS USED+THIS IS DESTROYED C EIG=RETURNED eigenVALUES IN ALGEBRAIC DESCENDING ORDER C VEC=RETURNED eigenVECTORS IN COLUMNS C C DIMENSION A(NN,NN),VEC(NN,NN),EIG(NN),W(NN),GAM(NN),BET(NN), *BSQ(NN),P(NN),Q(NN),IPO(NN),IORD(NN),IVP(NN) C DATA RHOSQ/1.0D-12/ ADUE=.50D0 ZERO=0.D0 UNO=1.0D0 DUE=2.0D0 C N=NN IF(N)10,550,10 10 N1=N-1 N2=N-2 GAM(1)=A(1,1) IF(N2) 180,170,20 20 DO 160 NR=1,N2 B=A(NR+1,NR) S=ZERO DO 30 I=NR,N2 30 S=S+A(I+2,NR)**2 C PREPARE FOR POSSIBLE BYPASS OF TRANSFORMATION A(NR+1,NR)=ZERO IF (S) 150,150,40 40 S=S+B*B SGN=UNO IF (B) 50,60,60 50 SGN=-UNO 60 SQRTS=SQRT(S) D=SGN/(SQRTS+SQRTS) TEMP=SQRT(ADUE+B*D) W(NR)=TEMP A(NR+1,NR)=TEMP D=D/TEMP B=-SGN*SQRTS C D IS FACTOR OF PROPORTIONALITY. NOW COMPUTE AND SAVE W VECTOR. C EXTRA SINGLY SUBSCRIPTED W VECTOR USED FOR SPEED. DO 70 I=NR,N2 TEMP=D*A(I+2,NR) W(I+1)=TEMP 70 A(I+2,NR)=TEMP C PREMULTIPLY VECTOR W BY MATRIX A TO OBTAIN P VECTOR. C SIMULTANEOUSLY ACCUMULATE DOT PRODUCT WP,(THE SCALAR K) WTAW=ZERO DO 120 I=NR,N1 I1=I+1 SUM=ZERO DO 80 J=NR,I 80 SUM=SUM+A(I1 ,J+1)*W(J) IF(N1-I1) 110,90,90 90 DO 100 J=I1,N1 100 SUM=SUM+A(J+1,I1 )*W(J) 110 P(I)=SUM WWWI=W(I) 120 WTAW=WTAW+SUM*WWWI C P VECTOR AND SCALAR K NOW STORED. NEXT COMPUTE Q VECTOR DO 130 I=NR,N1 130 Q(I)=P(I)-WTAW*W(I) C NOW FORM PAP MATRIX, REQUIRED PART DO 140 J=NR,N1 QJ=Q(J) WJ=W(J) DO 140 I=J,N1 140 A(I+1,J+1)=A(I+1,J+1)-DUE*(W(I)*QJ+WJ*Q(I)) 150 BET(NR)=B BSQ(NR)=B*B 160 GAM(NR+1)=A(NR+1,NR+1) 170 B=A(N,N1) BET(N1)=B BSQ(N1)=B*B GAM(N)=A(N,N) 180 BSQ(N)=ZERO C ADJOIN AN IDENTIFY MATRIX TO BE POSTMULTIPLIED BY ROTATIONS. DO 200 I=1,N DO 190 J=1,N 190 VEC(I,J)=ZERO 200 VEC(I,I)=UNO M=N SUM=ZERO NPAS=1 GO TO 330 210 SUM=SUM+SHIFT COSA=UNO G=GAM(1)-SHIFT PP=G PPBS=PP*PP+BSQ(1) PPBR=SQRT(PPBS) DO 300 J=1,M COSAP=COSA IF(PPBS)230,220,230 220 SINA=ZERO SINA2=ZERO COSA=UNO GO TO 270 230 SINA=BET(J)/PPBR SINA2=BSQ(J)/PPBS COSA=PP/PPBR C POSTMULTIPLY IDENTITY BY P-TRANSPOSE MATRIX NT=J+NPAS IF(NT-N)250,240,240 240 NT=N 250 DO 260 I=1,NT TEMP=COSA*VEC(I,J)+SINA*VEC(I,J+1) VEC(I,J+1)=-SINA*VEC(I,J)+COSA*VEC(I,J+1) 260 VEC(I,J)=TEMP 270 DIA=GAM(J+1)-SHIFT U=SINA2*(G+DIA) GAM(J)=G+U G=DIA-U PP=DIA*COSA-SINA*COSAP*BET(J) IF(J-M)290,280,290 280 BET(J)=SINA*PP BSQ(J)=SINA2*PP*PP GO TO 310 290 PPBS=PP*PP+BSQ(J+1) PPBR=SQRT(PPBS) BET(J)=SINA*PPBR 300 BSQ(J)=SINA2*PPBS 310 GAM(M+1)=G C TEST FOR CONVERGENCE OF LAST DIAGONAL ELEMENT NPAS=NPAS+1 IF(BSQ(M)-RHOSQ)320,320,350 320 EIG(M+1)=GAM(M+1)+SUM 330 BET(M)=ZERO BSQ(M)=ZERO M=M-1 IF(M)340,380,340 340 IF(BSQ(M)-RHOSQ)320,320,350 C TAKE ROOT OF CORNER 2 BY 2 NEAREST TO LOWER DIAGONAL IN VALUE C AS ESTIMATE OF eigenVALUE TO USE FOR SHIFT 350 A2=GAM(M+1) R2=ADUE*A2 R1=ADUE*GAM(M) R12=R1+R2 DIF=R1-R2 TEMP=SQRT(DIF*DIF+BSQ(M)) R1=R12+TEMP R2=R12-TEMP DIF=ABS(A2-R1)-ABS(A2-R2) IF(DIF)370,360,360 360 SHIFT=R2 GO TO 210 370 SHIFT=R1 GO TO 210 380 EIG(1)=GAM(1)+SUM C INITIALIZE AUXILIARY TABLES REQUIRED FOR REARRANGING THE VECTORS DO 390 J=1,N IPO(J)=J IVP(J)=J 390 IORD(J)=J C USE A TRANSPOSITION SORT TO ORDER THE eigenVALUES M=N GO TO 430 400 DO 420 J=1,M IF(EIG(J)-EIG(J+1))410,420,420 410 TEMP=EIG(J) EIG(J)=EIG(J+1) EIG(J+1)=TEMP ITEMP=IORD(J) IORD(J)=IORD(J+1) IORD(J+1)=ITEMP 420 CONTINUE 430 M=M-1 IF(M)400,440,400 440 IF(N1)450,490,450 450 DO 480 L=1,N1 NV=IORD(L) NP=IPO(NV) IF(NP-L)460,480,460 460 LV=IVP(L) IVP(NP)=LV IPO(LV)=NP DO 470 I=1,N TEMP=VEC(I,L) VEC(I,L)=VEC(I,NP) 470 VEC(I,NP)=TEMP 480 CONTINUE C BACK TRANSFORM THE VECTORS OF THE TRIPLE DIAGONAL MATRIX 490 DO 540 NRR=1,N K=N1 500 K=K-1 IF(K)540,540,510 510 SUM=ZERO DO 520 I=K,N1 520 SUM=SUM+VEC(I+1,NRR)*A(I+1,K) SUM=SUM+SUM DO 530 I=K,N1 530 VEC(I+1,NRR)=VEC(I+1,NRR)-SUM*A(I+1,K) GO TO 500 540 CONTINUE 550 RETURN END subroutine eq_plane(t,u,v,a) implicit double precision (a-h,o-z) dimension t(3),u(3),v(3),AA(3,3),DD(3,3),a(4) C Equazione del piano passante per i punti t,u,v nella forma canonica C a(1).x+a(2).y+a(3).z+a(4)=0 C con a(1),a(2),a(3)=coseni direttori a(4)=-distanza piano-origine C vedi International Tables for Crystallography II, p. 43, equaz. 2,3,4, 6 DD(1,1)=t(1) DD(1,2)=t(2) DD(1,3)=t(3) DD(2,1)=u(1) DD(2,2)=u(2) DD(2,3)=u(3) DD(3,1)=v(1) DD(3,2)=v(2) DD(3,3)=v(3) do i=1,3 call trasfm(DD,AA,1,1) call azzera(AA(1,i),1.d0,3) a(i)=det(AA,1) enddo call prods(a,a,dist,1,1,2) a(1)=a(1)/dist a(2)=a(2)/dist a(3)=a(3)/dist a(4)=-det(DD,1)/dist return end subroutine icosahed(XO,PESO,N,MN,MK,MD,II,MDEG,*) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) PARAMETER (NMG=120) integer out COMMON/output/out CHARACTER*6 NAME COMMON/AT4/NAME(NMA) common/matcel/PC(7),PCR(7),O(3,3),OI(3,3),G(3,3),GI(3,3),CS(12) COMMON/AT0/NA,NMO COMMON/AT2/SX(3,NMA),SIG(NMA),DXM(NMA),DCM,DCME,indwgh,indtol COMMON/ORIE/OR(3,3),OT(3,3),OTI(3,3),BARC(3),BARO(3),RIN(3) COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) DIMENSION XO(3,NMA),PESO(NMA),MD(NMA,2),MK(NMA,NMG),MN(NMA) DIMENSION DA(NMA),meq(NMA) dimension eqp(4),dp(5),vd(3,5),mp(5),io(5),A(3,3),B(3,3),V(3) logical ico RAD=57.29577951308232D0 do 2000 I1=1,N-4 if(MD(I1,1).ne.II)go to 2000 mp(1)=I1 io(1)=I1 in1=I1+1 do 1900 I2=in1,N-3 if(MD(I1,1).ne.II)go to 1900 mp(2)=I2 in2=I2+1 do 1800 I3=in2,N-2 if(MD(I3,1).ne.II)go to 1800 mp(3)=I3 call eq_plane(XO(1,I1),XO(1,I2),XO(1,I3),eqp) c write(out,*)'equazione del piano' c write(out,'(4f10.6,3i5)')eqp,mp(1),mp(2),mp(3) c write(out,'(4f10.6,3i5)')DXM(I1),DXM(I2),DXM(I3),eqp(4) C se il piano e' troppo vicino all'origine viene scartato if(eqp(4).lt.0.5)go to 1800 in3=I3+1 do 1700 I4=in3,N-1 if(MD(I4,1).ne.II)go to 1700 mp(4)=I4 d2=eqp(1)*XO(1,I4)+eqp(2)*XO(2,I4)+eqp(3)*XO(3,I4)+eqp(4) if(DABS(d2).gt.DXM(I4))go to 1700 c write(out,'(A20,i5,f10.5)')'IV atomo del piano',mp(4),d2 1300 in4=I4+1 do 1600 I5=in4,N if(MD(I5,1).ne.II)go to 1600 mp(5)=I5 d2=eqp(1)*XO(1,I5)+eqp(2)*XO(2,I5)+eqp(3)*XO(3,I5)+eqp(4) if(DABS(d2).gt.DXM(I5))go to 1600 c write(out,'(A20,i5,f10.5)')' V atomo del piano',mp(5),d2 C in mp ci sono i possibili 5 atomi equivalenti rispetto all'asse 5 dmin=100. dp(1)=dmin do i=2,5 io(i)=0 dp(i)=0.d0 call combli(XO,XO,vd,1.d0,-1.d0,I1,mp(i),1) call prods(vd,vd,dp(i),1,1,2) if(dp(i).lt.dmin)then dmin=dp(i) jj=i endif enddo io(2)=mp(jj) dp(jj)=100. dmin=100. do i=2,5 if(dp(i).lt.dmin.and.io(2).ne.mp(i))then dmin=dp(i) kk=i endif enddo io(5)=mp(kk) do i=2,5 if(mp(i).ne.io(2).and.mp(i).ne.io(5).and.io(3).eq.0)io(3)=mp(i) if(mp(i).ne.io(2).and.mp(i).ne.io(5).and.mp(i).ne.io(3))jj=mp(i) enddo call combli(XO,XO,vd,1.d0,-1.d0,io(2),io(3),1) call prods(vd,vd,d1,1,1,1) call combli(XO,XO,vd,1.d0,-1.d0,io(2),jj,1) call prods(vd,vd,d2,1,1,1) io(4)=jj if(d2.lt.d1)then io(4)=io(3) io(3)=jj endif C ora i cinque atomi sono ordinati C per controllo preliminare prima di passare alla verify, controllo che siano uguali ( a meno della tolleranza) le distanze 1-2 dmed=0. do i=1,5 k=i+1 if(k.gt.5)k=k-5 call combli(XO,XO,vd,1.d0,-1.d0,io(i),io(k),i) call prods(vd,vd,dp(i),i,i,2) dmed=dmed+dp(i) enddo dmed=dmed*0.2d0 c write(out,'(a5,5i5)')'io',io c write(out,'(a5,5f10.5)')'dp',dp c write(out,'(a5,5f10.5)')'dmed',dmed call AZZERA(A,0.d0,9) do i=1,5 if(DABS(dmed-dp(i)).gt.DXM(i)*.5)go to 1600 c asse C5 come somma dei vertici del pentagono A(1,3)=A(1,3)+XO(1,io(i)) A(2,3)=A(2,3)+XO(2,io(i)) A(3,3)=A(3,3)+XO(3,io(i)) enddo c write(out,'(a5,3f10.6)')'A(3)',(A(i,3),i=1,3) C mette il riferimento in modo che il primo asse 5 coincida con z e che C il secondo sia nel piano xz A(1,1)=XO(1,io(1))+XO(1,io(2))-XO(1,io(3))-XO(1,io(4))+ *XO(1,io(5)) A(2,1)=XO(2,io(1))+XO(2,io(2))-XO(2,io(3))-XO(2,io(4))+ *XO(2,io(5)) A(3,1)=XO(3,io(1))+XO(3,io(2))-XO(3,io(3))-XO(3,io(4))+ *XO(3,io(5)) C se la somma dei cinque atomi da' un vettore nullo, cioe' il pentagono C e' sul cerchio massimo (caso possibile per un gruppo di 30 atomi sugli C assi C2) passo ad altro pentagono 1400 call norm(A,3) call norm(A,1) call prodv(A,A,A,3,1,2) call norm(A,2) call prodv(A,A,A,2,3,1) call traspo(A,A,1,1) c write(out,'(a5,3(/,3f10.6))')'A',((A(i,k),i=1,3),k=1,3) call r_frame(A,XO,N) Costruzione dell'asse C5 call azzera(B,0.d0,9) ROT=72.d0/RAD CA=COS(ROT) CB=SIN(ROT) B(1,1)=CA B(2,2)=CA B(1,2)=-CB B(2,1)=CB B(3,3)=1.d0 c write(out,'(a5,3(/,3f10.6))')'B',((B(i,k),i=1,3),k=1,3) c verifica e ottimizza il primo asse 5 call verify(XO,B,MK,MN,MV,N) c write(out,*)'primo asse C5 MV,NMS',MV,NMS c write(out,'(30i3)')(MK(ll,NMS),ll=1,N) if(MV.eq.0)go to 1600 call opt_axis(XO,PESO,V,MK,N,2) c ottimizza l'sse x con la stessa tecnica precedente call azzera(A,0.d0,9) if(V(3).lt.0)then V(1)=-V(1) V(2)=-V(2) V(3)=-V(3) endif call trasfv(V,A,1,3) A(1,1)=XO(1,io(1))+XO(1,io(2))-XO(1,io(3))-XO(1,io(4))+ *XO(1,io(5)) A(2,1)=XO(2,io(1))+XO(2,io(2))-XO(2,io(3))-XO(2,io(4))+ *XO(2,io(5)) A(3,1)=XO(3,io(1))+XO(3,io(2))-XO(3,io(3))-XO(3,io(4))+ *XO(3,io(5)) 1500 call norm(A,3) call norm(A,1) call prodv(A,A,A,3,1,2) call norm(A,2) call prodv(A,A,A,2,3,1) call traspo(A,A,1,1) c write(out,'(a5,3(/,3f10.6))')'A',((A(i,k),i=1,3),k=1,3) call r_frame(A,XO,N) go to 2010 1600 continue 1700 continue 1800 continue 1900 continue 2000 continue c uscita normale (non ha trovato nessun asse C5) if(NMS.eq.1)return completa il gruppo C5. Indispensabile qui!!! 2010 call compl_gr(MK,N,*2015) C scelgo i due pentagoni piu' distanti dall'asse C5 per trovare C un asse C2. nel caso che gli atomi del sottoset non siano C tutti indipendenti, postrebbero esserci piu' di cinque atomi C sul piano parallelo al primo pentagono e quindi devo C testare tutti i possibili pentagoni 2015 DM=0.d0 do 2020 i=1,N if(MD(I1,1).ne.II)go to 2020 DA(i)=XO(1,I)**2+XO(2,I)**2 if(DM.gt.DA(i))go to 2020 kk=i DM=DA(i) 2020 continue io(1)=kk dp1=XO(3,kk) c write(out,*)'io(1),dp1',io(1),dp1 do k=1,4 k1=k+1 io(k1)=MK(io(k),2) dp1=dp1+XO(3,io(k1)) enddo dp1=-0.2d0*dp1 C in io(i) ci sono i cinque atomi del primo pentagono C in meq(i) vanno gli atomi del piano parallelo me=0 do 2030 i=1,N com=DABS(dp1-XO(3,i)) if(com.gt.DXM(i))go to 2030 C eliminazione (eventuale) degli atomi del primo piano do k=1,5 if(i.eq.io(k))go to 2030 enddo me=me+1 meq(me)=i 2030 continue c write(out,*)'me,meq',me c write(out,'(30i4)')(meq(i),i=1,me) C scelta degli atomi del secondo pentagono 2040 do i=1,me if(meq(i).ne.0)go to 2050 enddo C non ci sono assi C2!!! MDEG=1 return C i e' il primo atomo del nuovo pentagono 2050 mp(1)=meq(i) meq(i)=0 k=1 C cerca gli altri atomi del secondo penatgono e annulla i relativi meq do k=1,4 k1=k+1 mp(k1)=MK(mp(k),2) do l=1,me if(mp(k1).eq.meq(l))meq(l)=0 enddo enddo C somma vettori del primo pentagono opportunamente ruotati C costruzione dell'asse C2 perpendicolare a C5 e parallelo a x 2210 call azzera(A,0.d0,9) call azzera(vd,0.d0,15) c write(out,*)'io,mp' c write(out,'(5i4)')io,mp do 2250 i=1,5 k=io(i) do 2220 l=1,NMS if(MK(k,l).ne.io(1))go to 2220 call prodmv(SIM,XO,vd,l,k,3) call combli(vd,vd,vd,1.d0,1.d0,3,1,1) c write(out,'(a20,i5,3f10.5)')'io(i),vd(1)',io(i),(vd(kk,1),kk=1,3) go to 2230 2220 continue 2230 k=mp(i) do 2240 l=1,NMS if(MK(k,l).ne.mp(1))go to 2240 call prodmv(SIM,XO,vd,l,k,3) call combli(vd,vd,vd,1.d0,1.d0,3,2,2) c write(out,'(a20,i5,3f10.5)')'mp(i),vd(2)',mp(i),(vd(kk,2),kk=1,3) go to 2250 2240 continue 2250 continue vd(3,1)=0.d0 vd(3,2)=0.d0 call norm(vd,1) call norm(vd,2) c write(out,'(a20,3f10.5)')'vd(1)',(vd(kk,1),kk=1,3) c write(out,'(a20,3f10.5)')'vd(2)',(vd(kk,2),kk=1,3) cost=DCOS(36.1/RAD) do i=1,5 call prodmv(SIM,vd,vd,2,2,2) call prods(vd,vd,cosa,1,2,1) c write(out,'(a14,8f9.5)')'vd(1),vd(2)',(vd(kk,1),kk=1,3), c *(vd(kk,2),kk=1,3),cosa,cost if(cosa.ge.cost)go to 2300 enddo 2300 call combli(vd,vd,A,1.d0,1.d0,1,2,1) call norm(A,1) c write(out,'(a20,3f10.5)') 'A(1)',(A(kk,1),kk=1,3) A(3,3)=1.d0 call prodv(A,A,A,3,1,2) call norm(A,2) call prodv(A,A,A,2,3,1) call norm(A,1) call traspo(A,A,1,1) call r_frame(A,XO,N) c write(out,'(i4,3f9.5)')(lll,(XO(kk,lll),kk=1,3),lll=1,N) call azzera(B,0.d0,9) B(1,1)= 1.d0 B(2,2)=-1.d0 B(3,3)=-1.d0 call verify(XO,B,MK,MN,MV,N) c write(out,*)'asse C2 MV,NMS',MV,NMS if(MV.eq.1)go to 2310 call traspo(A,A,1,1) call r_frame(A,XO,N) call azzera(A,0.d0,9) A(1,2)=-1.d0 A(2,1)= 1.d0 A(3,3)= 1.d0 call r_frame(A,XO,N) call verify(XO,B,MK,MN,MV,N) c write(out,*)'asse C2 MV,NMS',MV,NMS if(MV.eq.1)go to 2310 call traspo(A,A,1,1) call r_frame(A,XO,N) go to 2040 c ico=.true. c call compl_gr(MK,N,*5000) c write(out,*)'asse C2,NMS',NMS Controlla il centro 2310 call azzera(B,0.d0,9) B(1,1)=-1.d0 B(2,2)=-1.d0 B(3,3)=-1.d0 call verify(XO,B,MK,MN,MV,N) c write(out,*)'CENTRO MV,NMS',MV,NMS Costruzione dell'asse C5 a |x|=63.43494882 da z e la cui priezione e' a 18 da x call azzera(A,0.d0,9) ROT=-63.43494882d0/RAD 3000 CA=COS(ROT) CB=SIN(ROT) call azzera(A,0.d0,9) A(1,1)=1.d0 A(2,2)= CA A(3,3)= CA A(2,3)=-CB A(3,2)= CB c rotazione 18 su z call azzera(B,0.d0,9) ROT18=-18.d0/RAD 3100 CA=COS(ROT18) CB=SIN(ROT18) call azzera(B,0.d0,9) B(1,1)= CA B(2,2)= CA B(1,2)=-CB B(2,1)= CB B(3,3)=1.d0 call prodmm(B,A,A,1,1,1) c write(out,*)' matrice di rotazione' c write(out,'(a12,2f10.6)')'ROT,ROT18',ROT,ROT18 c write(out,'(3f10.6)')((A(i,k),k=1,3),i=1,3) call prodmm(SIM,A,B,2,1,1) call traspo(A,A,1,1) call prodmm(A,B,A,1,1,1) call verify(XO,A,MK,MN,MV,N) c write(out,*)'secondo asse C5 ,MV,NMS',MV,NMS if(MV.eq.1)go to 4000 if(ROT18.lt.0.d0)then ROT18=-ROT18 go to 3100 endif if(ROT.gt.0.d0)return ROT=-ROT go to 3000 4000 call compl_gr(MK,N,*5000) 5000 return 1 END subroutine intpe(RIGA,B,C,IA,NB,*,*) implicit double precision (a-h,o-z) CHARACTER*(*) RIGA CHARACTER*23 CAR,CARA C INTERPRETA IL VETTORE A COME POSIZIONE EQUIVALENTE O C PREPARA LO STESSO PER LA STAMPA A PARTIRE DA B E C DIMENSION B(3,3,1),C(3,1) CAR='1234567890 /,+-xyz[]XYZ' GO TO (10,500),NB C INTERPRETAZIONE DI A 10 MCM=1 LR=LEN(RIGA) call compatta(RIGA,LR,KM) K=1 COST=1.d0 DO 15 I=1,3 C(I,IA)=0.d0 DO 15 J=1,3 15 B(I,J,IA)=0.d0 DO 200 I=1,KM C IDENTIFICAZIONE DEL CARATTERE CARA=RIGA(I:I) DO 20 J=1,23 IF(CAR(J:J).EQ.CARA)GO TO 30 20 CONTINUE 25 return 2 30 IF(J-11)40,200,100 C CARATTERE NUMERICO 40 IF(J.EQ.10)J=0 GO TO(50,60),MCM C NUMERATORE 50 IF(C(K,IA).NE.0.d0)C(K,IA)=C(K,IA)*10 C(K,IA)=C(K,IA)+COST*J COST=1.d0 GO TO 200 C DENOMINATORE 60 C(K,IA)=C(K,IA)/J MCM=1 GO TO 200 CARATTERI ALFABETICI 100 if(J.gt.20)J=J-5 J=J-11 GO TO(110,120,130,140,150,150,150,200,210),J C BARRA 110 MCM=2 GO TO 200 CAMBIO RIGA MATRICE (VIRGOLA) 120 K=K+1 GO TO 200 C SEGNO + 130 COST=1.d0 GO TO 200 C SEGNO - 140 COST=-1.d0 GO TO 200 C X,Y O Z 150 J=J-4 B(K,J,IA)=COST COST=1.d0 200 CONTINUE 210 RETURN 1 C SCRITTURA DEL VETTORE POSIZIONE EQUIVALENTE ALFANUMERICA 500 KM=LEN(RIGA) DO 510 I=1,KM 510 RIGA(i:i)=car(11:11) RIGA(12:12)=',' RIGA(23:23)=',' DO 1000 I=1,3 MCM=11*I DO 600 J=1,3 L=4-J IF(NINT(B(I,L,IA)).EQ.0)GO TO 600 MA=L+15 RIGA(MCM:MCM)=CAR(MA:MA) K=14 IF(B(I,L,IA).LT.(-0.1d0))K=15 MCM=MCM-1 RIGA(MCM:MCM)=CAR(K:K) MCM=MCM-1 600 CONTINUE MS=11 IF(ABS(C(I,IA)).LT.0.1d0)GO TO 1000 IF(C(I,IA).LT.0.d0)MS=15 DO 700 J=1,6 CA=C(I,IA)*J K=NINT(CA) IF(ABS(FLOAT(K)-CA).LT.0.1d0)GO TO 800 700 CONTINUE 800 IF(J.EQ.1)GO TO 900 RIGA(MCM:MCM)=CAR(J:J) MCM=MCM-1 RIGA(MCM:MCM)=CAR(12:12) MCM=MCM-1 900 J=IABS(MOD(K,10)) RIGA(MCM:MCM)=CAR(J:J) K=K/10 MCM=MCM-1 IF(K.GT.0)GO TO 900 RIGA(MCM:MCM)=CAR(MS:MS) 1000 CONTINUE KM1=KM-1 call compatta(RIGA(2:KM),KM1,K) if(RIGA(2:2).eq.'+')RIGA(2:2)=' ' do 1010 i=1,2 k=index(RIGA,',+') k1=k+1 if(k.ne.0)RIGA(k1:k1)=' ' 1010 continue call compatta(RIGA(2:KM),KM1,K) RETURN 1 END subroutine inv33(X,Y,K,L) implicit double precision (a-h,o-z) DIMENSION X(3,3,1),Y(3,3,1) C INVERSIONE DI UNA MATRICE 3X3. C=1.0/det(X,K) Y(1,1,L)=(X(2,2,K)*X(3,3,K)-X(2,3,K)*X(3,2,K))*C Y(2,1,L)=(X(2,3,K)*X(3,1,K)-X(2,1,K)*X(3,3,K))*C Y(3,1,L)=(X(2,1,K)*X(3,2,K)-X(2,2,K)*X(3,1,K))*C Y(1,2,L)=(X(1,3,K)*X(3,2,K)-X(1,2,K)*X(3,3,K))*C Y(2,2,L)=(X(1,1,K)*X(3,3,K)-X(1,3,K)*X(3,1,K))*C Y(3,2,L)=(X(1,2,K)*X(3,1,K)-X(1,1,K)*X(3,2,K))*C Y(1,3,L)=(X(1,2,K)*X(2,3,K)-X(1,3,K)*X(2,2,K))*C Y(2,3,L)=(X(1,3,K)*X(2,1,K)-X(1,1,K)*X(2,3,K))*C Y(3,3,L)=(X(1,1,K)*X(2,2,K)-X(1,2,K)*X(2,1,K))*C RETURN END FUNCTION ium(A,B,C,M,N) CONTROLLA SE A=B A MENO DI C C VERO ium=1 FALSO ium=0 implicit double precision (a-h,o-z) DIMENSION A(3,3,1),B(3,3,1) ium=0 DO 1 I=1,3 DO 1 J=1,3 IF(ABS(A(I,J,M)-B(I,J,N)).GT.C)RETURN 1 CONTINUE ium=1 RETURN END subroutine maiuscol(c) IMPLICIT DOUBLE PRECISION (A-H,O-Z) character c do 1 i=97,122 if(c.ne.CHAR(i))go to 1 c=CHAR(i-32) return 1 continue return end subroutine mass(NAT,ATOM) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) integer out COMMON/output/out CHARACTER*2 ATOM(NMA) CHARACTER*2 SIMBO(103),ITT REAL*8 MASSA(103),RAGGI(103) COMMON/AT1/X(3,NMA),AMAS(NMA),RCOV(NMA),MSP(NMA) DATA (MASSA(i),i=1,2)/1.00794,4.002602/ DATA (MASSA(i),i=3,10)/6.941,9.012182,10.811,12.011,14.00674, 115.9994,18.9984032,20.1797/ DATA (MASSA(i),i=11,18)/22.989768,24.3050,26.981539,28.0855, 130.97362,32.066,35.4527,39.948/ DATA (MASSA(i),i=19,36)/39.0983,40.078,44.955910,47.88, 1 50.9415,51.9961,54.93085,55.847,58.93320,58.69,63.546,65.39, 2 69.723,72.61,74.92159,78.96,79.904,83.80/ DATA (MASSA(i),i=37,54)/85.4678,87.62,88.90585,91.224, 1 92.90638,95.94,98,101.07,102.90550,106.42,107.8682,112.411, 2 114.82,118.710,121.75,127.60,126.90447,131.29/ DATA (MASSA(i),i=55,86)/132.90543,137.327,138.9055,140.115, 1 140.90765,144.24,145,150.36,151.965,157.25,158.92534,162.50, 2 164.93032,167.26,168.93421,173.04,174.967,178.49,180.9479, 3 183.85,186.207,190.2,192.22,195.08,196.96654,200.59,204.3833, 4 207.2,208.98037,209,210.0,222/ DATA (MASSA(i),i=87,103)/223,226.025,227.028,232.0381, 1 231.03588,238.0289,237.048,244.,243.,247.,247.,251.,252., 2 257.,258.,259.,260.0/ DATA (RAGGI(i),i=1,2)/.32,.93/ DATA (RAGGI(i),i=3,10)/1.23,0.90,.82,.770,.75,.73,.72,.71/ DATA (RAGGI(i),i=11,18)/1.54,1.4,1.18,1.11,1.06,1.02,.99,.98/ DATA (RAGGI(i),i=19,36)/2.03,1.74,1.44,1.32,1.22,1.18,1.17,1.17, 1 1.16,1.15,1.17,1.25,1.26,1.22,1.20,1.16,1.14,1.12/ DATA (RAGGI(i),i=37,54)/2.16,1.91,1.62,1.45,1.34,1.30,1.27,1.25, 1 1.25,1.28,1.34,1.48,1.44,1.41,1.40,1.36,1.33,1.31/ C per i raggi degli elementi dal Ce al Lu ho interpolato DATA (RAGGI(i),i=55,86)/2.35,1.98,1.69,1.67,1.65,1.63,1.61,1.59, 1 1.57,1.55,1.53,1.51,1.49,1.47,1.45,1.45,1.45,1.44,1.34,1.30,1.28, 2 1.26,1.27,1.30,1.34,1.49,1.48,1.47,1.46,1.76,1.45,1.45/ C per i raggi degli elementi dal Ce a Lu ho interpolato DATA (RAGGI(i),i=87,103)/2.35,1.98,1.69,1.67,1.65,1.63,1.61,1.59, 1 1.57,1.55,1.53,1.51,1.49,1.47,1.45,1.45,1.45/ DATA (SIMBO(i),i=1,2)/' H','HE'/ DATA (SIMBO(i),i=3,10)/'LI','BE',' B',' C',' N',' O',' F','NE'/ DATA (SIMBO(i),i=11,18)/'NA','MG','AL','SI',' P',' S','CL','AR'/ DATA (SIMBO(i),i=19,36)/' K','CA','SC','TI',' V','CR','MN','FE', 1'CO','NI','CU','ZN','GA','GE','AS','SE','BR','KR'/ DATA (SIMBO(i),i=37,54)/'RB','SR',' Y','ZR','Nb','MO','TC','RU', 1'RH','PD','AG','CD','IN','SN','SB','TE',' I','XE'/ DATA (SIMBO(i),i=55,86)/'CS','BA','LA','CE','PR','ND','PM', 1'SM','EU','GD','TB','DY','HO','ER','TM','YB','LU','HF','TA',' W', 2'RE','OS','IR','PT','AU','HG','TL','PB','BI','PO','AT','RN'/ DATA (SIMBO(i),i=87,103)/'FR','RA','AC','TH','PA',' U','NP', 1'PU','AM','CM','BK','CF','ES','FM','MD','NO','LR'/ DO J=1,NAT DO I=1,103 IF(ATOM(J).EQ.SIMBO(I)) THEN AMAS(J)=MASSA(I) RCOV(J)=RAGGI(I) c write(out,'(a2,2f7.3)')ATOM(J),AMAS(J),RCOV(J) GOTO 2 ENDIF ENDDO write(out,*) write(out,3) ATOM(J) 3 format(' ERROR: MISSING DATA FOR THE ATOM-',A2,'-') STOP 2 CONTINUE ENDDO DO J=1,NAT MSP(J)=0 ENDDO M=1 ITT=ATOM(1) 10 DO 1 K=1,NAT IF(ATOM(K).EQ.ITT)MSP(K)=M 1 CONTINUE M=M+1 DO J=1,NAT IF(MSP(J).EQ.0)THEN ITT=ATOM(J) MSP(J)=M GOTO 10 ENDIF ENDDO RETURN END subroutine momin(XO,P,RO,N) implicit double precision (a-h,o-z) c da un set di coordinate XO che vengono moltiplicate per la matrice RO C che puo' essere una matrice di ortogonalizzazione calcola i momenti di C inerzia utilizzando i pesi P COMMON/ORIE/OR(3,3),OT(3,3),OTI(3,3),BARC(3),BARO(3),RIN(3) DIMENSION RO(3,3),AA(3,3),AMOM(6),CAM(3),ORI(3,3) DIMENSION C1(3),C2(3),C3(3),C4(3),C5(3),C6(3),IC1(3),IC2(3),IC3(3) DIMENSION XO(3,1),P(1) DO 100 I=1,3 RIN(I)=0.d0 AMOM(I)=0.d0 100 AMOM(I+3)=0.d0 SUM=0.d0 DO 110 I=1,N call prodmv(RO,XO,XO,1,I,I) PESO=P(I) RIN(1)=RIN(1)+PESO*XO(1,I) RIN(2)=RIN(2)+PESO*XO(2,I) RIN(3)=RIN(3)+PESO*XO(3,I) 110 SUM=SUM+PESO DO 115 I=1,3 RIN(I)=RIN(I)/SUM 115 BARO(I)=RIN(I) DO 120 I=1,N PESO=P(I) XO(1,I)=XO(1,I)-BARO(1) XO(2,I)=XO(2,I)-BARO(2) XO(3,I)=XO(3,I)-BARO(3) CAM(1)=XO(1,I)*PESO CAM(2)=XO(2,I)*PESO CAM(3)=XO(3,I)*PESO AMOM(1)=AMOM(1)+CAM(1)*XO(1,I) AMOM(2)=AMOM(2)-CAM(1)*XO(2,I) AMOM(3)=AMOM(3)-CAM(1)*XO(3,I) AMOM(4)=AMOM(4)+CAM(2)*XO(2,I) AMOM(5)=AMOM(5)-CAM(2)*XO(3,I) AMOM(6)=AMOM(6)+CAM(3)*XO(3,I) 120 CONTINUE AA(1,1)=AMOM(4)+AMOM(6) AA(2,2)=AMOM(1)+AMOM(6) AA(3,3)=AMOM(1)+AMOM(4) AA(1,2)=AMOM(2) AA(2,1)=AMOM(2) AA(1,3)=AMOM(3) AA(3,1)=AMOM(3) AA(2,3)=AMOM(5) AA(3,2)=AMOM(5) call eigen(AA,OR,RIN,C1,C2,C3,C4,C5,C6,IC1,IC2,IC3,3) C__________________ C modifica per rendere OR machine-independent am1=-1 am2=-1 do i=1,3 if(DABS(OR(i,1)).gt.am1)then am1=DABS(OR(i,1)) n1=i endif if(DABS(OR(i,2)).gt.am2)then am2=DABS(OR(i,2)) n2=i endif enddo if(OR(n1,1).lt.0.0)then do i=1,3 OR(i,1)=-OR(i,1) enddo endif if(OR(n2,1).lt.0.0)then do i=1,3 OR(i,2)=-OR(i,2) enddo endif CALL prodv(OR,OR,OR,1,2,3) C__________________ call traspo(OR,OR,1,1) call prodmm(OR,RO,OT,1,1,1) call inv33(OT,OTI,1,1) call inv33(RO,ORI,1,1) call prodmv(ORI,BARO,BARC,1,1,1) do 130 i=1,N 130 call prodmv(OR,XO,XO,1,I,I) RETURN END subroutine norm(A,I) implicit double precision (a-h,o-z) C NORMALIZZAZIONE DEL VETTORE A DIMENSION A(3,1) call prods(A,A,B,I,I,2) DO 1 J=1,3 1 A(J,I)=A(J,I)/B RETURN END subroutine opt_axis(XO,PESO,V3,MK,N,IAX) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) PARAMETER (NMG=120) integer out COMMON/output/out COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) DIMENSION XO(3,NMA),MK(NMA,NMG),V3(3),PESO(NMA),VS(3) DIMENSION IV(NMA),CO(3,NMA),W(NMA),nge(NMA) C fa in modo che gli la somma delle distanze dall'elemento di simmetria IAX C sia minimo do i=1,N iv(i)=0 enddo ng=0 cost=det(SIM,IAX) 1000 do I=1,N if(iv(i).eq.0)go to 1100 enddo go to 1300 1100 ng=ng+1 C ng = numero del gruppo di atomi equivalente rispetto all'operazione IAX iv(i)=ng do j=1,3 CO(j,ng)=XO(j,I) enddo W(ng)=PESO(I) nge(ng)=1 k=I 1200 l=MK(k,IAX) if(l.eq.I)go to 1000 iv(l)=ng C somma tutti gli atomi del gruppo con segno opportuno do j=1,3 CO(j,ng)=CO(j,ng)+cost*XO(j,l) enddo W(ng)=W(ng)+PESO(l) nge(ng)=nge(ng)+1 k=l go to 1200 C ha assegnato tutti gli atomi ad un gruppo ed il peso ai gruppi 1300 continue do i=1,ng W(I)=W(i)/nge(I) enddo call azzera(V3,0.d0,3) dism=0.d0 cerca il gruppo piu' lontano dall'elemento di simmetria (im) do i=1,ng call prodmv(SIM,CO,VS,IAX,i,1) call combli(CO,VS,VS,1.d0,cost,i,1,1) call prods(VS,VS,dis,1,1,1) c write(out,'(i3,6f10.5)')i,(CO(j,i),j=1,3),VS if(dis.gt.dism)then dism=dis im=i endif enddo call azzera(V3,0.d0,3) call prods(CO,CO,dim,im,im,2) c write(out,'(6him,dim,i3,3f10.5)')im,dim c somma fra loro tutti i gruppi con un segno determinato dall'angolo c fra il gruppo in esame ed il gruppo im do i=1,ng call prods(CO,CO,abcos,im,i,1) call prods(CO,CO,di,i,i,2) absabcos=abs(abcos) if(absabcos.gt.1.d-5)then cost=W(I)*abcos/absabcos call combli(V3,CO,V3,1.d0,cost,1,i,1) endif enddo c il vettore risultante e' normalizzato e diventera' una riga della c matrice di rotazione definitiva call norm(V3,1) if(V3(1).gt.0.d0)return V3(1)=-V3(1) V3(2)=-V3(2) V3(3)=-V3(3) return end subroutine out_bond(XO,IEQAT,IASU,N) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) integer out COMMON/output/out CHARACTER*6 NAME COMMON/AT4/NAME(NMA) COMMON/AT5/LEGAMI(20,NMA),NLEG(NMA) DIMENSION T(3),V(3),IASU(NMA),tableg(20,NMA),XO(3,NMA),IEQAT(NMA) CHARACTER*80 riga write(out,*) c write(out,'(a6,3f10.6)')(NAME(IEQAT(I)),(XO(k,i),k=1,3),i=1,N) RAD=57.29577951308232D0 riga=' ' nb=0 do 1400 i=1,N if(IASU(i).ne.2)go to 1400 N1=NLEG(i) do 1300 k=1,N1 l=LEGAMI(k,i) call combli(XO,XO,T,1.d0,-1.d0,i,l,1) call prods(t,t,tableg(k,i),1,1,2) nc=nb*26+1 nb=nb+1 write(riga(nc:nc+25),'(a6,a1,a6,f7.4,6x)')NAME(IEQAT(I)),'-', *NAME(IEQAT(l)),tableg(k,i) if(nb.eq.3)then write(out,'(a)')riga nb=0 riga=' ' endif 1300 continue 1400 continue if(nb.ne.0)write(out,'(a)')riga write(out,*) riga=' ' nb=0 do 2000 i=1,N N1=NLEG(i) if(IASU(i).ne.2)go to 2000 do 1900 k=1,N1-1 k1=k+1 i1=LEGAMI(k,i) call combli(XO,XO,T,1.d0,-1.d0,i1,i,1) do 1800 l=k1,N1 i2=LEGAMI(l,i) call combli(XO,XO,V,1.d0,-1.d0,i2,i,1) call prods(T,V,cost,1,1,1) cost=cost/(tableg(k,i)*tableg(l,i)) CALL arccos(cost,ang) nc=nb*35+1 nb=nb+1 write(riga(nc:nc+27),'(a6,2(a1,a6),f8.3,7x)')NAME(IEQAT(I1)), *'-',NAME(IEQAT(I)),'-',NAME(IEQAT(I2)),ang if(nb.eq.2)then write(out,'(a)')riga nb=0 riga=' ' endif 1800 continue 1900 continue 2000 continue if(nb.ne.0)write(out,'(a)')riga write(out,*) return END subroutine pratba(B,A,C,I,J,K) implicit double precision (a-h,o-z) C T C TRASFORMAZIONE C=A BA DIMENSION A(3,3,1),B(3,3,1),C(3,3,1),D(3,3) call PRODMM(B,A,C,I,J,K) call TRASPO(A,D,J,1) call PRODMM(D,C,C,1,K,K) RETURN END subroutine prodmm(A,B,C,I,J,N) implicit double precision (a-h,o-z) DIMENSION A(3,3,1),B(3,3,1),C(3,3,1),X(3,3) DO 1 K=1,3 DO 1 L=1,3 X(K,L)=0.0d0 DO 1 M=1,3 Y=A(K,M,I) Z=B(M,L,J) 1 X(K,L)=X(K,L)+Y*Z DO 2 K=1,3 DO 2 L=1,3 2 C(K,L,N)=X(K,L) RETURN END subroutine prodms(A,B,C,M,N) implicit double precision (a-h,o-z) DIMENSION A(3,3,1), B(3,3,1) X=C DO 1 I=1,3 DO 1 J=1,3 Y=A(I,J,M) 1 B(I,J,N)=X*Y RETURN END subroutine prodmv(A,B,C,I,J,K) implicit double precision (a-h,o-z) DIMENSION A(3,3,1),B(3,1),C(3,1),X(3) DO 1 L=1,3 X(L)=0.0d0 DO 1 M=1,3 Y=A(L,M,I) Z=B(M,J) 1 X(L)=X(L)+Y*Z DO 2 L=1,3 2 C(L,K)=X(L) RETURN END subroutine prods(A,B,C,I,J,K) implicit double precision (a-h,o-z) C PRODOTTO SCALARE A*B=C (K=1) C MODULO VETTORE A (K=2) DIMENSION A(3,1),B(3,1) X=0.0d0 DO 1 L=1,3 Y=A(L,I) Z=B(L,J) 1 X=X+Y*Z C=X IF(K.EQ.2)C=DSQRT(X) RETURN END subroutine prodv(A,B,C,I,J,K) implicit double precision (a-h,o-z) C PRODOTTO VETTORIALE A*B=C DIMENSION A(3,1),B(3,1),C(3,1) X1=A(1,I) X2=A(2,I) X3=A(3,I) Y1=B(1,J) Y2=B(2,J) Y3=B(3,J) C(1,K)=X2*Y3-X3*Y2 C(2,K)=-X1*Y3+X3*Y1 C(3,K)=X1*Y2-X2*Y1 RETURN END subroutine r_frame(CO,XO,N) C ruota il riferimento e le coordinate implicit double precision (a-h,o-z) COMMON/ORIE/OR(3,3),OT(3,3),OTI(3,3),BARC(3),BARO(3),RIN(3) DIMENSION CO(3,3),XO(3,*) c write(6,*)'rotazione del riferimento' c write(6,'(3f10.6)')((CO(i,k),k=1,3),i=1,3) DO 1010 I=1,N 1010 call prodmv(CO,XO,XO,1,I,I) call prodmm(CO,OR,OR,1,1,1) call prodmm(CO,OT,OT,1,1,1) call inv33(OT,OTI,1,1) c write(6,*)'nuove coordinate' c write(6,'(i5,3f10.6)')(i,(XO(k,i),k=1,3),i=1,N) return end subroutine rms_min(tm) implicit double precision (a-h,o-z) C PARAMETER (NMA=1000) PARAMETER (NMG=120) PARAMETER (NMV=NMA*NMG) common/rmsmin/RV(3,3),ppo(3,NMV),ppu(3,NMV),npu dimension tm(3),t(3) passo=0.02 call azzera(RV,0.d0,9) call azzera(tm,0.d0,3) rm=crms(tm) 1 ind1=0 do i1=1,3 t(1)=tm(1)+(2-i1)*passo do i2=1,3 t(2)=tm(2)+(2-i2)*passo do i3=1,3 t(3)=tm(3)+(2-i3)*passo arms=crms(t) C write(6,*)t,arm,i1,i2,i3 if(arms.lt.rm)then ind1=1 rm=arms c write(6,*)'rm,passo',rm,passo call trasfv(t,tm,1,1) endif enddo enddo enddo if(ind1.eq.1)goto 1 if(passo.lt.0.0002)return passo=passo*0.5 go to 1 end subroutine s_coor(SI,XO,XS,CO,D,DEL,DLM,CSM,MK,ID,ISU,nes,N) implicit double precision (a-h,o-z) PARAMETER (NMA=1000) PARAMETER (NMG=120) integer out COMMON/output/out COMMON/SIMME1/RMS(3,NMA),RMST DIMENSION CO(3,NMA),XO(3,NMA),XS(3,NMA),SI(3,3,NMG),D(NMA) DIMENSION DEL(4),DLM(4),V3(3),ID(4),MK(NMA,NMG),ISU(NMG) DIMENSION dc(NMA),cf(NMA) calcola le coordinate simmetrizzate XS c write(out,'(i3,3f10.5,)')(kk,(XO(ll,kk),ll=1,3),kk=1,N) NAZ=3*N call azzera(CO,0.d0,NAZ) call azzera(XS,0.d0,NAZ) call azzera(RMS,0.d0,NAZ) call azzera(DEL,0.d0,4) call azzera(DLM,0.d0,4) call azzera(cf,0.d0,N) C ISU=indice delle matrici del sottogruppo in considerazione C nes=numero di elementi del sottogruppo c write(out,'(a3,12i3)')'ISU',(ISU(J),J=1,nes) DO 2800 I=1,N call prods(XO,XO,D(I),I,I,2) DO 2790 J=1,nes K=MK(I,ISU(j)) call prodmv(SI,XO,V3,ISU(J),I,1) call combli(XS,V3,XS,1.d0,1.d0,K,1,K) 2790 continue c write(out,'(i2,4f9.5)')I,(XO(kk,I),kk=1,3),D(I) 2800 CONTINUE do i=1,N do J=1,3 XS(J,I)=XS(J,I)/dfloat(nes) enddo call prods(XS,XS,dc(i),i,i,2) enddo do I=1,N DO J=1,nes K=MK(I,ISU(j)) cf(I)=cf(I)+D(K) enddo if(dc(i).gt.1.d-01)cf(i)=cf(i)/(dc(i)*nes) do J=1,3 XS(J,I)=XS(J,I)*cf(i) enddo enddo c write(out,'(i2,6f10.6)')1,(XS(kk,1),kk=1,3),dc(1),cf(1) C il fattore di correzione cf=media delle distanze atomi equivalenti baricentro, C corregge le coordinate simmetrizzate facendo si che queste equivalgano alla C media per pura rotazione degli atomi equivalenti calcola le rms e gli scostamenti medi (DEL) e massimi DLM do j=1,4 ID(j)=0 enddo DO 2820 i=1,N do j=1,3 V3(j)=XO(j,i)-XS(j,i) com=dabs(V3(j)) DEL(j)=DEL(J)+com if(com.gt.DLM(j))then DLM(j)=com ID(j)=I endif enddo C massima deviazione dalla simmetria in angstrom call prods(V3,V3,com,1,1,2) DEL(4)=DEL(4)+com if(DLM(4).lt.com)then DLM(4)=com ID(4)=I endif 2820 continue calcola la continuous symmetry measure (CSM) CSM=0.d0 do j=1,nes do i=1,N K=MK(i,ISU(j)) call prodmv(SI,XO,V3,ISU(J),i,1) call combli(XS,V3,V3,1.d0,-1.d0,K,1,1) RMS(1,K)=RMS(1,K)+V3(1)*V3(1) RMS(2,K)=RMS(2,K)+V3(2)*V3(2) RMS(3,K)=RMS(3,K)+V3(3)*V3(3) call prods(V3,V3,com,1,1,1) CSM=CSM+com enddo enddo cost=1.d0/dfloat(nes) RMST=0.d0 do I=1,N do j=1,3 RMST=RMST+RMS(j,i) RMS(j,i)=sqrt(RMS(j,i)*cost) enddo enddo RMST=sqrt(RMST/dfloat(nes*N)) CSM=CSM*1.d2/dfloat(nes*N) return end subroutine schoenfl(MDEG) C From the group matrices gives the Schoenflies symbol implicit double precision (a-h,o-z) PARAMETER (NMA=1000) PARAMETER (NMG=120) integer out COMMON/output/out COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CMS(NMG),CSMT,NMS,MTG(NMG,NMG) COMMON/SIMME1/RMS(3,NMA),RMST character*3 lbls,pt*7 character S1,S2,Sn*2,tot*4 COMMON/SIMME3/IPOGROU,lbls(NMG) character*5 pointgr(32) DATA (pointgr(i),i=1,32) */'C1 ','Ci ','C2 ','Cs ','C2h','D2 ','C2v','D2h','C4 ','S4 ','C4h' *,'D4 ','C4v','D2d','D4h','C3 ','C3i','D3 ','C3v','D3d','C6 ','C3h' *,'C6h','D6 ','C6v','D3h','D6h','T ','Th ','O ','Td ','Oh '/ IPOGROU=0 maxasp=0 maxasi=0 mplane=0 invers=0 nplane=0 lbls(1)='E' C_______________________________________________________ C gruppi lineari C per questi gruppi non stampo le matrici di simmetria ma solo il gruppo C if(MDEG.eq.5)then pt='C(inf)v' do i=1,NMS call ax_order(SIM,i,m,msign,inv) if(inv.eq.1)go to 1000 enddo NMS=1 go to 5200 1000 pt='D(inf)h' NMS=1 go to 5200 endif C_______________________________________________________ if(NMS.eq.1)then pt='C1 ' IPOGROU=1 return endif do i=2,NMS tot=' ' S1=' ' S2=' ' Sn=' ' call ax_order(SIM,i,m,msign,inv) if(m.gt.maxasp.and.msign.eq.1)maxasp=m if(m.gt.maxasi.and.msign.eq.-1)maxasi=m if(m.eq.2.and.mplane.eq.0)mplane=1 c write(out,'(a20,4i3)')'i,msign,m,inv',i,msign,m,inv S1='C' if(m.eq.2)then Centro if(inv.eq.1)then invers=1 Sn='i ' go to 1100 endif C piano Sn='2 ' C asse 2 if(msign.eq.-1)then Sn='s ' nplane=nplane+1 endif go to 1100 endif C asse improprio m>2 write(Sn,'(i2)')m if(msign.eq.-1)S1='S' 1100 tot=S1//Sn//S2 call compatta(tot,4,k) lbls(i)=tot(1:3) enddo c write(out,*)'maxasp maxasi mplane nplane invers MDEG NMS' c write(out,'(7i7)')maxasp,maxasi,mplane,nplane,invers,MDEG,NMS pt=' ' if(NMS.eq.48)pt='Oh ' if(NMS.eq.60)pt='I ' if(NMS.eq.120)pt='Ih ' if(pt.ne.' ')go to 5100 if(MDEG.eq.3)then if(NMS.eq.12)pt='T23' if(invers.eq.1)pt='Th ' if(maxasi.eq.4)pt='Td ' if(pt.eq.' ')pt='O ' go to 5100 endif if(NMS.eq.2)then pt=lbls(2) go to 5000 endif S2=' ' S1='C' C Cn e S2n dove n=maxasp write(Sn,'(I2)')maxasp if(MOD(NMS,2).eq.0.and.maxasi.eq.NMS)then S1='S' write(Sn,'(I2)')maxasi if(NMS.ne.6)go to 5000 pt='C3i' go to 5100 endif if(NMS.eq.maxasp)go to 5000 if(NMS.lt.maxasi)then write(Sn,'(I2)')maxasi S1='S' go to 5000 endif C Dn, Cnv e Cnh S1='D' if(NMS.eq.maxasp*2)then if(nplane.eq.0)go to 5000 S1='C' S2='v' if(nplane.eq.maxasp)go to 5000 S2='h' go to 5000 endif S2='h' if(nplane.eq.maxasp)S2='d' 1 format(/,' Schoenflies symbol = ',a7,' CSM =',f8.4,5x, *'Molecular RMS =',f8.4) 2 format(' CSM,see: Zabrodsky et al. (1993) JACS, 115, 8278-8298') 5000 tot=' ' tot=S1//Sn//S2 call compatta(tot,4,k) pt=tot(1:3) 5100 do m=1,32 if(pt.eq.pointgr(m))IPOGROU=m enddo 5200 write(out,1)pt,CSMT,RMST write(out,2) if(out.eq.9)then write(*,1)pt,CSMT,RMST write(*,2) endif RETURN END subroutine stasy implicit double precision (a-h,o-z) C C STAMPA LE POSIZIONI EQUIVALENTI C c character*120 linea PARAMETER (NMG=120) integer out COMMON/output/out character*80 linea character*40 riga COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) dimension TRL(3,NMG) character*3 lbls COMMON/SIMME3/IPOGROU,lbls(NMG) call azzera(TRL,0.d0,3*NMS) i1=-39 I=1 if(NMS.gt.1)I=2 2 format(2a40) write(out,2)(' Symmetry element its CSM and Max.Diff. ',k=1,I) DO 320 M=1,NMS mm=M call intpe(riga(1:39),SIM,TRL,mm,2,*310,*310) 310 i1=i1+40 i2=i1+39 i3=I1+21 i4=I1+37 write(linea(i1:i2),3)M,lbls(M),riga 3 format(I2,') [',a3,'] ',a30) write(linea(i3:i4),'(2f8.4)')CSM(M),DEV(M) c if(i1.eq.81)then if(i1.eq.41)then write(out,'(1x,a)')linea do 315 kk=1,80 c do 315 kk=1,120 linea(kk:kk)=' ' 315 continue i1=-39 endif 320 CONTINUE if(i1.ge.1)write(out,'(1x,a)')linea(1:i1+39) RETURN END subroutine trasfm(A,B,I,J) implicit double precision (a-h,o-z) C TRASFERIMENTO DI UNA MATRICE DIMENSION A(3,3,1),B(3,3,1) DO 1 K=1,3 DO 1 L=1,3 1 B(K,L,J)=A(K,L,I) RETURN END subroutine trasfv(A,B,I,J) implicit double precision (a-h,o-z) C TRASFERIMENTO DI UN VETTORE DIMENSION A(3,1),B(3,1) B(1,J)=A(1,I) B(2,J)=A(2,I) B(3,J)=A(3,I) RETURN END subroutine traspo(A,B,M,N) implicit double precision (a-h,o-z) DIMENSION A(3,3,1),B(3,3,1),C(3,3) DO 1 I=1,3 DO 1 J=1,3 1 C(I,J)=A(J,I,M) call trasfm(C,B,1,N) RETURN END subroutine verify(XO,A,MK,MN,MV,N) implicit double precision (a-h,o-z) C VERIFICA SE A E' UNA MATRICE DI SIMMETRIA A MENO DI DXM C IN CASO AFFERMATIVO PONE MV=1, NMS=NMS + 1, SIM(I,J,NMS=A(I,J) C MK(I,NMS) CONTIENE IL NUMERO DELL'ATOMO OTTENUTO STASFORMANDO I CON C L'OPERAZIONE NMS PARAMETER (NMA=1000) PARAMETER (NMG=120) integer out COMMON/output/out COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) COMMON/AT2/SX(3,NMA),SIG(NMA),DXM(NMA),DCM,DCME,indwgh,indtol DIMENSION XO(3,1),MN(1),MK(NMA,1),A(3,3),V(3),W(3),iat(NMA) MV=1 Controlla per prima cosa che ls matrice A non sia gia' compresa in SIM do i=1,NMS if(ium(A,SIM,1.d-2,1,i).eq.1)go to 6 enddo NMS=NMS+1 c write(out,'(i3,3f10.6,i4)')(i,(XO(k,i),k=1,3),MN(i),i=1,N) c write(out,'(3f10.6)')((A(k,i),i=1,3),k=1,3) DO 4 I=1,N MK(I,NMS)=0 JJ=0 FMIN=DXM(I) call prodmv(A,XO,V,1,I,1) c write(out,'(3f10.6)')V DO 2 J=1,N IF(MN(I).NE.MN(J))GO TO 2 call combli(XO,V,W,1.d0,-1.d0,J,1,1) call prods(W,W,F,1,1,2) c write(out,*)'I,J,NMS,F',I,J,NMS,F IF(F.GT.FMIN)GO TO 2 FMIN=F JJ=J 2 CONTINUE c write(out,*)'JJ, se JJ=0 rifiuta',JJ IF(JJ.EQ.0)GO TO 5 MK(I,NMS)=JJ 4 CONTINUE c write(out,*)'MK' c write(out,'(29i3)')(MK(k,NMS),k=1,N) C analisi della matrice di moltiplicazione do k=1,N iat(k)=0 enddo do 10 k=1,N if(iat(k).ne.0)go to 10 iat(k)=1 neq=1 l=k 7 m=MK(l,NMS) if(m.eq.k)go to 9 iat(m)=1 neq=neq+1 if(neq.gt.N)go to 5 l=m go to 7 9 continue call ax_order(A,1,nord,msign,inv) c write(out,*)'nord,neq,inv',nord,neq,inv C il numero di atomi equivalenti per l'operazione di ordine m puo' essere C solo 1 o 2 (atomi giacenti sull'elemento di simmetria) m (caso normale) C o 2m ( quando si ha un elemento riducibile 3/m 5/m etc.) if(neq.eq.1.or.neq.eq.2)go to 10 if(neq.eq.nord)go to 10 if(neq.eq.2*nord)go to 10 go to 5 10 continue c write(out,*)'accettata, nord=',nord c write(out,'(3f10.5)')((A(i,j),j=1,3),i=1,3) call trasfm(A,SIM,1,NMS) RETURN 5 NMS=NMS-1 6 MV=0 c write(out,*)'return da verify,NMS=',NMS RETURN END subroutine work(MO) C MO = molecola da simmetrizzare C X = coordinate frazionarie originali C NA = numero totale di atomi C MOL(I) = molecola a cui appartiene l'atomo I C MSP(I) = indicatore della specie dell'atomo I C AMAS(I) = massa dell'atomo I C XO = coordinate soli atomi da simmetrizzare C N = numero di atomi da simmetrizzare C SIM = matrici di simmetria del gruppo molecolare C NMS = ordine del gruppo di simmetria trovato implicit double precision (a-h,o-z) PARAMETER (maxorder=8) PARAMETER (NMA=1000) PARAMETER (NMG=120) PARAMETER (NMV=NMA*NMG) common/rmsmin/RV(3,3),ppo(3,NMV),ppu(3,NMV),npu integer out COMMON/output/out CHARACTER*6 NAME,NNAME CHARACTER*52 comment COMMON/AT4/NAME(NMA) common/matcel/PC(7),PCR(7),O(3,3),OI(3,3),G(3,3),GI(3,3),CS(12) COMMON/AT0/NA,NMO COMMON/AT1/X(3,NMA),AMAS(NMA),RCOV(NMA),MSP(NMA) COMMON/AT2/SX(3,NMA),SIG(NMA),DXM(NMA),DCM,DCME,indwgh,indtol COMMON/AT3/MOL(NMA),MLG(NMA) COMMON/ORIE/OR(3,3),OT(3,3),OTI(3,3),BARC(3),BARO(3),RIN(3) COMMON/SIMME/SIM(3,3,NMG),DEV(NMG),CSM(NMG),CSMT,NMS,MTG(NMG,NMG) COMMON/SIMME1/RMS(3,NMA),RMST COMMON/SIMME2/AL(3,3),ALX(3,3) CHARACTER*3 lbls COMMON/SIMME3/IPOGROU,lbls(NMG) DIMENSION MK(NMA,NMG),MN(NMA),XO(3,NMA),XS(3,NMA),MD(NMA,2) DIMENSION D(NMA),AA(3,3),CO(3,NMA),IEQAT(NMA) DIMENSION PESO(NMA),DEL(4),DELM(4),IDM(4),delta(3,NMA) DIMENSION V(3),W(3),NNAME(4),ISU(NMG),IASU(NMA) CHARACTER type(2),car type(1)=' ' type(2)='*' PIG=3.14159265358979D0 deliner=0.1d0 call azzera(SIM,0.d0,9*NMG) call azzera(AL,0.d0,9) SIM(1,1,1)=1.d0 SIM(2,2,1)=1.d0 SIM(3,3,1)=1.d0 AL(1,1)=1.d0 AL(2,2)=1.d0 AL(3,3)=1.d0 F=0.d0 IES=0 N=0 NMS=1 c selezione gruppo atomi da simmetrizzare write(out,1)MO 1 format(//,38H SYMMETRIZATION OF GROUP NR.,i2,//) if(out.eq.9)write(6,1)MO do 450 i=1,NA call azzera(AA,0.d0,9) AA(1,1)=SX(1,I) AA(2,2)=SX(2,I) AA(3,3)=SX(3,I) call prodmv(O,SX,AA,1,I,1) SIG(I)=0.d0 den=0.d0 do 435 k=1,3 if(AA(K,K).le.0.d0)go to 435 SIG(I)=SIG(I)+AA(k,k) den=den+1.d0 435 continue if(den.ne.0.d0)then SIG(I)=SIG(I)/den else if(indtol.eq.3.or.indwgh.eq.3)then write(out,*)' WARNING: there is at least an atom with s.u.=0' write(out,*) write(out,*)' CHANGE THE WEIGHTING AND TOLERANCE SCHEME!!!' STOP endif endif if(MOL(i).ne.MO)go to 450 N=N+1 IASU(i)=2 IEQAT(N)=i do k=1,3 XO(k,N)=X(k,i) enddo MN(N)=MSP(i) if(indwgh.eq.1) PESO(N)=AMAS(i) if(indwgh.eq.2) PESO(N)=1.0D0 if(indwgh.eq.3)then if(SIG(i).lt..00001d0)then write(out,*)' ERROR: THE S.U. OF AN ATOM IS ZERO!!!' stop endif PESO(N)=1.d0/(SIG(I)*SIG(I)) endif MK(N,1)=N 450 continue porig=PESO(1) if(N.le.2)then write(out,*)' GROUP WITH LESS THAN THREE ATOMS' return endif C C ORTOGONALIZZAZIONE DELLE COORDINATE DEGLI ATOMI DA SIMMETRIZZARE C CALCOLO DELLA DISTANZA MEDIA DAL BARICENTRO DEL GRUPPO C 460 call momin(XO,PESO,O,N) c do i=1,3 c write(OUT,'(3f10.5,5x,3f10.5)')(OR(i,k),k=1,3),(OT(i,k),k=1,3) c enddo C C DM=maximum length of vectors XO(I) C DM=0.d0 DO 480 I=1,N call prods(XO,XO,D(I),I,I,2) c write(OUT,'(a9,i3,4f10.5)')'atom,D(I)',I,(XO(k,i),k=1,3),D(i) if(D(i).gt.DM)then DM=D(i) endif 480 continue C do 500 I=1,N DXM(I)=DCM if(indtol.eq.2)DXM(I)=DCM*D(I)/DM if(indtol.eq.3)DXM(I)=SIG(I)*DCM 500 continue MDEG=0 delr=DABS(RIN(1)-RIN(2))/RIN(1) IF(delr.LE.deliner)then MDEG=1 endif delr=DABS(RIN(2)-RIN(3))/RIN(1) IF(delr.LE.deliner)then MDEG=MDEG+2 endif NDEGA=MDEG+1 nd=(MDEG+1)/2+1 write(out,*)'PRINCIPAL INERTIA MOMENTS and DEGENERATION DEGREE' write(out,'(3f10.2,i10,/)')RIN,nd C_______________________________________________________________________ C gruppi lineari if(MDEG.eq.1)then do i=1,N com=DSQRT(XO(1,i)*XO(1,i)+XO(2,i)*XO(2,i)) if(com.gt.DXM(i))go to 990 enddo call azzera(CO,0.d0,9) CO(1,1)=-1.d0 CO(2,2)=-1.d0 CO(3,3)= 1.d0 call verify(XO,CO,MK,MN,MV,N) CO(3,3)=-1.d0 call verify(XO,CO,MK,MN,MV,N) MDEG=5 go to 2520 endif C_______________________________________________________________________ 990 GO TO (1020,1020,1000,2000),NDEGA C C ASSE UNICO=X RUOTA LE COORDINATE IN MODO DA AVERE Z COME ASSE UNICO C MODIFICA DI CONSEGUENZA ANCHE LA MATRICE DI ORIENTAZIONE C 1000 call azzera(CO,0.d0,9) MDEG=1 CO(3,1)=1.d0 CO(1,3)=1.d0 CO(2,2)=-1.d0 call r_frame(CO,XO,N) COM=RIN(3) RIN(3)=RIN(1) RIN(1)=com C ASSE UNICO Z 1020 IU=3 IB=2 IC=1 C RICERCA DELL'ASSE DI ORDINE MORD 1100 MORD=maxorder+1 IF(MDEG.EQ.0)MORD=3 NASS=1 1110 MORD=MORD-1 IF(MORD.EQ.3.AND.MDEG.EQ.3)MORD=2 IF(MORD.EQ.1)GO TO 1210 C GENERAZIONE DELLA MATRICE DI SIMMETRIA 1120 COST=1.d0 call azzera(CO,0.d0,9) ROT=2.d0*PIG/dfloat(MORD) CA=COS(ROT) CB=SIN(ROT) 1140 CO(IB,IB)=CA*COST CO(IC,IC)=CO(IB,IB) CO(IB,IC)=CB*COST CO(IC,IB)=-CO(IB,IC) CO(IU,IU)=COST call verify(XO,CO,MK,MN,MV,N) c write(out,*)' 1140 tentativo con la matrice:',NINT(cost),MORD c write(out,'(3f9.5)')((CO(kkk,lll),lll=1,3),kkk=1,3) c write(out,*)' NMS MORD NASS MV' c write(out,'(5i5)')NMS,MORD,NASS,MV IF(MV.EQ.1)GO TO 1150 IF(COST.EQ.-1.d0)GO TO 1110 COST=-1.d0 mmez=MORD/2 m2=2*mmez if(m2.ne.MORD)go to 1140 C se MORD=2*mmez ed mmez=dispari (assi -6=3/m -10=5/m) si ha un c elemento riducibile quindi e' inutile testarlo m2=2*mmez/2 if(mmez-m2.eq.1)go to 1110 GO TO 1140 C HA TROVATO UN ELEMENTO DI SIMMETRIA 1150 IF(MORD.EQ.3.OR.MORD.EQ.6)IES=1 c write(OUT,*)'MDEG,MORD',MDEG,MORD if(MDEG.eq.1.and.MORD.gt.2)go to 1300 C RICERCA ELEMENTI DI SIMMETRIA SU NUOVI ASSI 1205 if(MORD.gt.2)go to 1110 1210 NASS=NASS+1 C AL COMPLETAMENTO DEL GRUPPO IF(NASS.EQ.4)GO TO 2500 C RICICLO SUL SECONDO O TERZO ASSE MORD=2 I=IU IU=IB IB=IC IC=I GO TO 1120 C RICERCA DELLA ROTAZIONE PER PORTARE EVENTUALI ELEMENTI DI SIMMETRIA C A COINCIDERE CON GLI ASSI DI RIFERIMENTO 1300 IF(N.le.1)then write(out,*)'Single atoms, not a molecule' return endif C SELEZIONE DELL'ATOMO DA USARE PER DETERMINARE LA ROTAZIONE call azzera(CO,0.d0,18) CO(3,3)=1.d0 C divisione in gruppi equidistanti dall'asse unico do I=1,N D(I)=DSQRT(XO(1,I)*XO(1,I)+XO(2,I)*XO(2,I)) MD(I,1)=0 MD(I,2)=0 enddo IAI=0 C IAI=NUMERO DI ENNUPLE DI ATOMI AVENTI DISTANZE COMPRESE IN UN RANGE C X< D 0 2780 if(MDEG.eq.0)go to 2790 do i=1,NMS ISU(i)=i enddo call s_coor(SIM,XO,XS,CO,D,DEL,DELM,CSMT,MK,IDM,ISU,NMS,N) npb=0 do k=1,NMS do i=1,N npu=npb+i call trasfv(XO,ppu,i,npu) npu=npb+MK(i,k) call prodmv(SIM,XS,ppo,k,i,npu) enddo npb=npb+N enddo npu=npb call rms_min(V) arms=crms(V) c write(6,*)'v,RV' c write(6,'(3f10.5)')V,RV call r_frame(RV,XO,N) Calcolo simmetria per sottogruppi di operazioni 2790 CSM(1)=0.d0 DEV(1)=0.d0 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c write(out,*)'GROUP MULTIPLICATION TABLE' c write(out,*) c nst=NMS c nis=25 c if(NMS.gt.24)nst=24 c do i=1,NMS c write(out,'(24i3)')(MTG(i,j),j=1,nst) c enddo c write(out,*) c if(nst.ne.NMS)then c do i=1,NMS c write(out,'(24i3)')(MTG(i,j),j=nis,NMS) c enddo c endif C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c simmetrizzazione per sottogruppi di una solo elemento (e sue potenze) do 2950 I=2,NMS nes=1 ISU(1)=nes k=i 2800 nes=nes+1 ISU(nes)=k l=MTG(i,k) if(l.eq.1)go to 2900 k=l go to 2800 2900 call s_coor(SIM,XO,XS,CO,D,DEL,DELM,comt,MK,IDM,ISU,nes,N) DEV(I)=DELM(4) CSM(I)=comt 2950 continue CALCOLA LE COORDINATEE SIMMETRIZZATE PER L'INTERO GRUPPO 3000 do i=1,NMS ISU(i)=i enddo call s_coor(SIM,XO,XS,CO,D,DEL,DELM,CSMT,MK,IDM,ISU,NMS,N) c write(out,*) c write(out,*)'FRACTIONAL COORDINATES OF THE CENTRE OF MASS' c write(out,'(3f12.6,/)')BARC c write(out,*)'ORTHOGONAL COORDINATES OF THE CENTRE OF MASS' c write(out,'(3f12.6,/)')BARO c write(out,*)'PRINCIPAL INERTIA MOMENTS' c write(out,'(3f10.2,/)')RIN write(out,*) write(out,*)' ORTHOGONALIZATION MATRIX' write(out,'(3f12.6)')((OT(k,i),i=1,3),k=1,3) write(out,*) write(out,*)'ATOM ORTHOGONAL COORDINATES VECTORS TO MAKE * SYMMETRICAL THE GROUP' c do I=1,N c write(out,'(24I4)')(MK(I,K),K=1,NMS) c enddo C Recupero degli atomi con MOL<0 N1=N do 3100 i=1,NA if(MOL(i).ne.-MO)go to 3100 N1=N1+1 IEQAT(N1)=i MN(N1)=MSP(i) MK(N1,1)=-N1 call combli(X,BARC,XO,1.d0,-1.d0,i,1,N1) call prodmv(OT,XO,XO,1,N1,N1) XS(1,N1)=XO(1,N1) XS(2,N1)=XO(2,N1) XS(3,N1)=XO(3,N1) 3100 continue NSTART=N+1 IF(N1.eq.N)go to 3400 do 3300 I=NSTART,N1 do 3200 L=2,NMS MK(i,l)=0 FMIN=10000. call prodmv(SIM,XO,V,L,I,1) DO 3150 K=NSTART,N1 IF(MN(I).NE.MN(K))GO TO 3150 call combli(XO,V,W,1.d0,-1.d0,K,1,1) call prods(W,W,F,1,1,2) if(F.le.FMIN)then KK=K FMIN=F endif 3150 continue if(FMIN.LT.2.*DCME)then MK(I,L)=KK else MK(I,1)=0 endif 3200 continue 3300 continue c do I=1,N1 c write(out,'(24I4)')(MK(I,K),K=1,NMS) c enddo do 3350 i=NSTART,N1 if(MK(I,1).eq.0)go to 3350 do k=2,NMS call prodmv(SIM,XO,V,k,i,1) l=MK(i,k) call combli(XS,V,XS,1.d0,1.d0,l,1,l) enddo 3350 continue do 3360 I=NSTART,N1 if(MK(I,1).eq.0)go to 3360 do k=1,3 XS(k,I)=XS(k,I)/dfloat(NMS) enddo 3360 continue 3400 write(out,*) call asymunit(MK,IASU,N1,NMS) do i=1,N1 call combli(XS,XO,delta,1.d0,-1.d0,i,i,i) write(out,'(a6,2x,2(3f10.5,5x))')NAME(IEQAT(I)), * (XO(k,I),k=1,3),(delta(k,i),k=1,3) enddo write(out,*) write(out,'(5x,a34,12x,a)')'SYMMETRIZED ORTHOGONAL COORDINATES', *'ATOMIC R.M.S.' ind1=0 ind2=0 do I=1,N1 comment=' ' car=' ' if(MK(i,1))3410,3430,3420 3410 comment='#' ind1=1 3420 car=type(IASU(I)) go to 3440 3430 comment='$' ind2=2 3440 if(MK(i,1).gt.0)then write(out,'(a6,i2,2(3f9.5,2x),a1)')NAME(IEQAT(I)), * MOL(IEQAT(I)),(XS(k,I),k=1,3),(RMS(k,I),k=1,3),car else write(out,'(a6,i2,3f9.5,29x,a1,1x,a1)')NAME(IEQAT(I)), * MOL(IEQAT(I)),(XS(k,I),k=1,3),comment,car endif C IASU(I)=2 se l'atomo rappresenta l'unita' asimmetrica altrimenti IASU(I)=1 enddo ind1=ind1+ind2 write(out,*) write(out,*)'* Atom defining the asymmetric unit for the found sym *metry group.' if(ind1.eq.0)go to 3460 if(ind1.eq.2)go to 3450 write(out,*) write(out,*)'# This atom was symmetrized but NOT used to find the *symmetry group and to calculate CMS, RMS and so on.' 3450 if(ind1.eq.1)go to 3460 write(out,*) write(out,*)'$ It was IMPOSSIBLE to symmetrize this atom accordin *g to the found symmetry group and within the given tolerance.' 3460 write(out,*) do i=1,4 DEL(i)=DEL(I)/dfloat(N) id1=idm(i) if(IDM(i).ne.0) then ie1=IEQAT(id1) NNAME(i)=NAME(ie1) else NNAME(i)=' --- ' endif enddo write(out,*) write(out,'(1x,a29,4F10.5,/)')'AVERAGE DIFFERENCE ON X,Y,Z,D',DEL write(out,'(1x,a29,4F10.5)')'MAXIMUM DIFFERENCE ON X,Y,Z,D',DELM write(out,'(1x,a29,4a10,/)')'DUE TO THE ATOMS', *(NNAME(i),i=1,4) if(out.eq.9)then write(*,*) write(*,'(1x,a29,4F10.5,/)')'AVERAGE DIFFERENCE ON X,Y,Z,D',DEL write(*,'(1x,a29,4F10.5)')'MAXIMUM DIFFERENCE ON X,Y,Z,D',DELM write(*,'(1x,a29,4a10,/)')'DUE TO THE ATOMS', * (NNAME(i),i=1,4) endif C do i=1,N do j=1,NMS K=MK(i,j) call prodmv(SIM,XS,V,j,i,1) do l=1,3 diff=abs(XS(l,k)-V(l)) if(diff.gt.0.001d0)then write(out,*) write(out,*) * ' ========> CAUTION!!! <=======' write(out,*) write(out,*) ' THE SYMMETRY OPERATIONS ARE NOT CONSISTENT WITH THE * EQUIVALENTS POINTS.' write(out,'(a27,f10.5)')' DIFFERENCE FOUND =',diff write(out,*) * ' ========> TRY TO REDUCE THE CONSTANT OF TOLERANCE <=======' c write(out,'(3(3f11.7,2x))')(XS(kk,i),kk=1,3),V,(XS(kk,k),kk=1,3) go to 3500 c RETURN endif enddo enddo enddo C 3500 write(out,*)'Bond lengths and bond angles after symmetrization' call bond(XS,IASU,N1) call out_bond(XS,IEQAT,IASU,N) if(PC(7).gt.2.d0)then write(out,'(7x,a)')' SYMMETRIZED FRACTIONAL COORDINATES' do 3600 i=1,N1 k=IEQAT(i) call prodmv(OTI,XS,X,1,i,k) call combli(X,BARC,X,1.d0,1.d0,k,1,k) write(out,'(1x,a6,i2,6f9.5)')NAME(k),MOL(k),(X(j,k),j=1,3), *(SX(j,k),j=1,3) 3600 continue endif call schoenfl(MDEG) if(IPOGROU.ne.0.and.(IPOGROU.lt.16.or.IPOGROU.gt.27))then call stasy else write(out,6) do 3800 k=1,NMS write(out,8)k,CSM(k),DEV(k),lbls(K) 8 format(/,i3,' CSM =',f7.2,5x,'MAX. DIFF. (Angstrom)=',F6.4, *5x,'TYPE ',a3) write(out,7)((SIM(I,J,K),J=1,3),I=1,3) 3800 continue endif IF(IES.NE.1)RETURN if(IPOGROU.le.0.or.IPOGROU.gt.32)RETURN AL(2,2)=SQRT(3.D0)*0.5D0 AL(1,2)=-0.5D0 AL(1,1)=1.d0 AL(3,3)=1.d0 call inv33(AL,ALX,1,1) write(out,*)' SYMMETRY OPERATIONS IN HEXAGONAL COORDINATES' DO 3900 I=1,NMS call prodmm(ALX,SIM,AA,1,I,1) call prodmm(AA,AL,SIM,1,1,I) 3900 continue call stasy write(out,*) write(out,11) 11 format(' OBLIQUE COORDINATES (HEXAGONAL SYSTEM)') do 4000 I=1,N1 call prodmv(ALX,XS,CO,1,I,I) write(out,'(1x,a8,3f10.5)')NAME(IEQAT(I)),(CO(J,I),J=1,3) 4000 continue 4100 RETURN 3 format(' NO SYMMETRY EXISTS WITHIN GIVEN TOLERANCE') 6 format(/,' SYMMETRY GROUP MATRICES',/) 7 format(3(1X,3F16.10,/)) END