PROGRAM IMPLANT IMPLICIT none logical lex character*4 key real*4 bonding(88),r,xi,xj,yi,yj,zi,zj,rt,rb,ris,rs, 1c1,q1,p1,tol2,rj,amvm,rfac,cm1(3),pm1(5),qm1(3,3), 1d,dmin,told,tolq,rms,rmsmin real*8 am(3,3),amv(3,3),rc(3),etem(3),amv1(3,3),rc1(3),am1(3,3) integer*4 ik,iout,i,j,ic,nat,nc,alist,natm,ia,ja,iks,iop,no, 1ix,jx,jj,nmol,jmol,istart,IERR,ind,itry,iie,nats,imin,jjmin allocatable ::ik(:),r(:,:),iout(:),alist(:,:),natm(:), 1c1(:,:),p1(:,:),iks(:), 2q1(:,:,:),amvm(:,:,:),rj(:,:),ind(:,:), 1rs(:,:),ris(:,:) data bonding/ 10.50,0.32,0.98,1.28,0.95,0.87,0.82,0.80,0.78,0.77, 10.76,1.59,1.41,1.23,1.16,1.11,1.20,1.04,1.03,2.08,1.79,1.49,1.37, 11.27,1.23,1.22,1.22,1.21,1.20,1.22,1.30,1.31,1.27,1.25,1.21,1.19, 11.17,2.21,1.96,1.67,1.50,1.39,1.35,1.32,1.30,1.30,1.33,1.39,1.53, 11.49,1.46,1.45,1.41,1.38,1.36,2.40,2.03,1.74,1.70,1.70,1.69,1.68, 11.67,1.90,1.66,1.64,1.64,1.63,1.62,1.61,1.75,1.61,1.49,1.39,1.35, 11.33,1.31,1.32,1.35,1.39,1.54,1.53,1.52,1.51,1.51,1.50,1.49,0.10/ logical lwr character*10 s10 character*80 fn,fns c WRITE(6,3000) 3000 FORMAT(' Takes geometry in BIG.X, identifies',/, 1 ' molecules, and makes atomic maps with a molecule',/ 1 ' in SMALL.X' 1 ' Input: BIG.X',/, 1 ' SMALL.X',/,/, 3 ' Output: CCT.TEM',/,/) c tolerance dipole: told=2. c tolerance quadrupole: tolq=20.0 c tolerance geometry: tol2=1.0 lwr=.true. rfac=1.1 fn='BIG.X' fns='SMALL.X' inquire(file='IMPLANT.OPT',exist=lex) if(lex)then open(2,file='IMPLANT.OPT') 100 read(2,2000,end=2200,err=2200)key 2000 format(a4) if(key.eq.'TOLD')read(2,*)told if(key.eq.'TOLQ')read(2,*)tolq if(key.eq.'TOL2')read(2,*)tol2 if(key.eq.'WRIT')read(2,*)lwr if(key.eq.'RFAC')read(2,*)lwr if(key.eq.'FILE')read(2,*)fn if(key.eq.'SMAL')read(2,*)fns goto 100 2200 close(2) endif write(6,606)tol2 606 format(' tolerance : ',F10.1,' A') open(2,file=fn) read(2,*) read(2,*)nat allocate(ik(nat),r(nat,3),iout(nat),rj(nat,3),ind(8,nat)) do 3 i=1,nat iout(i)=0 read(2,*)ik(i),(r(i,j),j=1,3) 3 if(ik(i)+1.lt.1.or.ik(i)+1.gt.88)call report('type out of range') close(2) write(6,*)nat,' atoms in '//fn open(2,file=fns) read(2,*) read(2,*)nats allocate(iks(nats),rs(nats,3),ris(nats,3)) do 37 i=1,nats read(2,*)iks(i),(rs(i,j),j=1,3) 37 if(iks(i)+1.lt.1.or.iks(i)+1.gt.88)call report('type off range') close(2) write(6,*)nats,' atoms in '//fns c find free atom, and connected atoms nmol=1 iout(1)=1 991 nc=1 99 ic=0 do 13 i=1,nat c is atom i close to some in the molecule?: if(iout(i).eq.0)then xi=r(i,1) yi=r(i,2) zi=r(i,3) rb=bonding(ik(i)+1)*rfac do 2 j=1,nat if(iout(j).eq.nmol)then rt=(rb+bonding(ik(j)+1)*rfac)**2 xj=r(j,1) yj=r(j,2) zj=r(j,3) if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.rt)then c atom j is bound to i, which is close=take i, too: iout(i)=nmol ic=ic+1 nc=nc+1 goto 13 endif endif 2 continue endif 13 continue if(ic.ne.0)goto 99 do 131 i=1,nat if(iout(i).eq.0)then c some atoms is left, start new molecule: nmol=nmol+1 iout(i)=nmol goto 991 endif 131 continue write(6,*)nmol,' molecules' allocate(alist(nmol,nat),natm(nmol)) do 8 i=1,nmol 8 natm(i)=0 do 1 i=1,nat natm(iout(i))=natm(iout(i))+1 1 alist(iout(i),natm(iout(i)))=i if(lwr)then do 9 i=1,nmol write(6,*)' molecule ',i,': ',natm(i),' atoms' 9 write(6,603)(alist(i,jj),jj=1,natm(i)) 603 format(15i5) endif c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule write(6,*) write(6,*)'Source molecule:' call cd1(cm1,pm1,qm1,nats,rs,iks) c diagonalize moment of inertia ans save transformation matrix: do 27 ix=1,3 do 27 jx=1,3 27 am1(ix,jx)=dble(qm1(ix,jx)) IERR=0 CALL TRED12(3,am1,amv1,rc1,2,IERR,etem) c check that we have a right-handed system: CALL chkxx(amv1,rc1) allocate(amvm(nmol,3,3),c1(nmol,3),p1(nmol,5),q1(nmol,3,3)) c calculate moments c 444444444444444444444444444444444444444444444444444444444444 do 4 i=1,nmol write(6,*)' Molecule ',i,':' call sets(s10,i,istart) c center and dipoles: call cd(c1,p1,q1,natm,i,alist,nmol,nat,r,ik) c diagonalize moment of inertia ans save transformation matrix: do 271 ix=1,3 do 271 jx=1,3 271 am(ix,jx)=dble(q1(i,ix,jx)) IERR=0 CALL TRED12(3,am,amv,rc,2,IERR,etem) c check that we have a right-handed system CALL chkxx(amv,rc) do 4 ix=1,3 do 4 jx=1,3 4 amvm(i,ix,jx)=real(amv(ix,jx)) c 444444444444444444444444444444444444444444444444444444444444 c i c c write final atomic map: open(19,file='CCT.TEM') write(19,901)nmol 901 format('LNAMES',/,'T',/,'POLYMER',/,I9) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c shift molecule to the centre and transform to inertia axes: call shm1(nats,amv1,rs,ris,cm1) c indicate enantiomer: iie=0 iop=0 c c look if contained in some jmol: do 29 jmol=1,nmol write(6,*)' molecule ',jmol if(natm(jmol).eq.nats.and. 1 abs(pm1(4)-p1(jmol,4)).lt.told.and. 1 abs(pm1(5)-p1(jmol,5)).lt.tolq)then write(6,*)'similar dipole and second moment' imin=0 rmsmin=1.0e5 call iniind(ind,nat) do 291 itry=0,7 write(6,*)' itry ',itry c transform jmol to inertia axes: c and rotate x y z inversion invx invy invz: call shm(nat,nmol,jmol,amvm,natm,alist,r,rj,c1) call rr(itry,nat,nmol,jmol,natm,alist,rj) c find overlapped atoms: rms=0.0 do 33 ia=1,nats xi=ris(ia,1) yi=ris(ia,2) zi=ris(ia,3) c minimum distance to atom of the same kind: dmin=1.0e5 jjmin=1 do 34 jj=1,natm(jmol) ja=alist(jmol,jj) if(iks(ia).eq.ik(ja))then d=(rj(ja,1)-xi)**2+(rj(ja,2)-yi)**2+(rj(ja,3)-zi)**2 if(d.lt.dmin)then dmin=d jjmin=ja endif endif 34 continue if(dmin.lt.tol2)ind(itry+1,jjmin)=ia 33 rms=rms+dmin rms=rms/real(nats) write(6,6001)rms 6001 format(' rms:',f12.4) if (rms.lt.rmsmin)then rmsmin=rms imin=itry endif 291 continue c itry write(6,6101)rmsmin,imin 6101 format(' best rms:',f12.4,', itry:',i2) call shm(nat,nmol,jmol,amvm,natm,alist,r,rj,c1) call rr(imin,nat,nmol,jmol,natm,alist,rj) c all atoms had analogues, fragment imol is part of fragment jmol: if(imin.ge.4)then c for enantiomers - inform only for the first time: if(iie.eq.0)then iie=iie+1 write(6,706)jmol 706 format(' INFO: ',i4,' IS ENANTIOMER') endif else c for true isomers - take it,write overlap: iop=iop+1 write(19,*)iop write(19,902)(i-(i/1000)*1000,i=1,nat) 902 format(20i3) no=0 do 261 i=1,nat if(ind(imin+1,i).ne.0)no=no+1 write(19,9021)ind(imin+1,i) 9021 format(i3,$) 261 if(20*(i/20).eq.i)write(19,9022)i 9022 format(' -- ',i6) if(20*(nat/20).ne.nat)write(19,9022)nat write(6,*)no,' atoms overlap' endif else write(6,*)'different dipole or second moment' endif 29 continue c jmol c write(6,*)'CCT.TEM written' close(19) write(6,*)iop,' overlaps' if(iop.ne.nmol)write(6,*)' IOP <> NMOL !!!!' end subroutine report(s) character*(*) s write(6,*)s stop end subroutine cd(cm,pm,qm,natc,imol,clist,nmol,nat,r,ik) c for molecule imol c calculate center (cm) c dipole moment = pm1_3 c size of dipole moment = pm4 c second moment = qm c trace of second moment = pm5 implicit none integer*4 nmol,nat,natc(*),imol,clist(nmol,nat),ix,jj,jx,ik(*) real*4 cm(nmol,3),pm(nmol,5),qm(nmol,3,3),r(nat,3) do 51 ix=1,3 cm(imol,ix)=0.0 do 52 jj=1,natc(imol) 52 cm(imol,ix)=cm(imol,ix)+r(clist(imol,jj),ix) 51 cm(imol,ix)=cm(imol,ix)/real(natc(imol)) do 53 ix=1,3 pm(imol,ix)=0.0 do 53 jj=1,natc(imol) 53 pm(imol,ix)=pm(imol,ix)+real(ik(clist(imol,jj)))* 1(r(clist(imol,jj),ix)-cm(imol,ix)) pm(imol,4)=sqrt(pm(imol,1)**2+pm(imol,2)**2+pm(imol,3)**2) do 54 ix=1,3 do 54 jx=1,3 qm(imol,ix,jx)=0.0 do 54 jj=1,natc(imol) 54 qm(imol,ix,jx)=qm(imol,ix,jx)+(r(clist(imol,jj),ix)-cm(imol,ix)) 1* (r(clist(imol,jj),jx)-cm(imol,jx)) pm(imol,5)=qm(imol,1,1)+qm(imol,2,2)+qm(imol,3,3) write(6,602)(cm(imol,ix),ix=1,3),natc(imol), 1pm(imol,4),pm(imol,5) 602 format(' Center : ',3f12.3,/, 1 ' Mass : ',i12,/, 1 ' Dipole : ',f12.3,/, 1 ' Quadrupole: ',f12.3) return end subroutine cd1(cm,pm,qm,nat,r,ik) c for a molecule c calculate center (cm) c dipole moment = pm1_3 c size of dipole moment = pm4 c second moment = qm c trace of second moment = pm5 implicit none integer*4 nat,ix,jj,jx,ik(*) real*4 cm(3),pm(5),qm(3,3),r(nat,3) c geometry center: do 51 ix=1,3 cm(ix)=0.0 do 52 jj=1,nat 52 cm(ix)=cm(ix)+r(jj,ix) 51 cm(ix)=cm(ix)/real(nat) c first moment: do 53 ix=1,3 pm(ix)=0.0 do 53 jj=1,nat 53 pm(ix)=pm(ix)+real(ik(jj))*(r(jj,ix)-cm(ix)) c mass moment length: pm(4)=sqrt(pm(1)**2+pm(2)**2+pm(3)**2) c second geometry moment: do 54 ix=1,3 do 54 jx=1,3 qm(ix,jx)=0.0 do 54 jj=1,nat 54 qm(ix,jx)=qm(ix,jx)+(r(jj,ix)-cm(ix))*(r(jj,jx)-cm(jx)) c second geometry moment trace: pm(5)=qm(1,1)+qm(2,2)+qm(3,3) write(6,602)(cm(ix),ix=1,3),nat,pm(4),pm(5) 602 format(' Center : ',3f12.3,/, 1 ' Nat : ',i12,/, 1 ' Dipole : ',f12.3,/, 1 ' Quadrupole: ',f12.3) return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION A(N,N),Z(N,N),D(*),E(*) 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 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 EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) 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 rr(xx,nat,nmol,imol,natc,clist,ri) c shift cluster to the central molecule and transform to inertia axes: implicit none integer*4 ix,jx,imol,natc(*),nmol,nat,clist(nmol,nat),ia,ii,xx real*4 tr(3,3),ri(nat,3),xi,yi,zi do 32 ix=1,3 do 32 jx=1,3 32 tr(ix,jx)=0.0 c xx=0: identity if(xx.eq.0)return c xx=1: C2x-rotation if(xx.eq.1)then tr(1,1)= 1.0 tr(2,2)=-1.0 tr(3,3)=-1.0 endif c xx=2: C2y-rotation if(xx.eq.2)then tr(1,1)=-1.0 tr(2,2)= 1.0 tr(3,3)=-1.0 endif c xx=3: C2z-rotation if(xx.eq.3)then tr(1,1)=-1.0 tr(2,2)=-1.0 tr(3,3)= 1.0 endif c xx=4: inversion if(xx.eq.4)then tr(1,1)=-1.0 tr(2,2)=-1.0 tr(3,3)=-1.0 endif c xx=5: inversion and C2x rotation if(xx.eq.5)then tr(1,1)=-1.0 tr(2,2)= 1.0 tr(3,3)= 1.0 endif c xx=6: inversion and C2y rotation if(xx.eq.6)then tr(1,1)= 1.0 tr(2,2)=-1.0 tr(3,3)= 1.0 endif c xx=7: inversion and C2z rotation if(xx.eq.7)then tr(1,1)= 1.0 tr(2,2)= 1.0 tr(3,3)=-1.0 endif do 28 ii=1,natc(imol) ia=clist(imol,ii) xi=ri(ia,1) yi=ri(ia,2) zi=ri(ia,3) ri(ia,1)=tr(1,1)*xi+tr(2,1)*yi+tr(3,1)*zi ri(ia,2)=tr(1,2)*xi+tr(2,2)*yi+tr(3,2)*zi 28 ri(ia,3)=tr(1,3)*xi+tr(2,3)*yi+tr(3,3)*zi return end subroutine shm(nat,nmol,imol,amvm,natc,clist,r,ri,c1) c shift cluster to the central molecule and transform to inertia axes: implicit none integer*4 ix,jx,imol,natc(*),nmol,nat,clist(nmol,nat),ia,ii real*4 amvm(nmol,3,3),tr(3,3),r(nat,3),ri(nat,3),xi,yi,zi, 1c1(nmol,3) do 32 ix=1,3 do 32 jx=1,3 32 tr(ix,jx)=amvm(imol,ix,jx) do 28 ii=1,natc(imol) ia=clist(imol,ii) xi=r(ia,1)-c1(imol,1) yi=r(ia,2)-c1(imol,2) zi=r(ia,3)-c1(imol,3) ri(ia,1)=tr(1,1)*xi+tr(2,1)*yi+tr(3,1)*zi ri(ia,2)=tr(1,2)*xi+tr(2,2)*yi+tr(3,2)*zi 28 ri(ia,3)=tr(1,3)*xi+tr(2,3)*yi+tr(3,3)*zi return end subroutine shm1(nat,amvd,r,ri,c1) c shift molecule to the centre and transform to inertia axes: implicit none integer*4 nat,ia,i,j real*4 amvm(3,3),r(nat,3),ri(nat,3),xi,yi,zi,c1(3) real*8 amvd(3,3) do 1 i=1,3 do 1 j=1,3 1 amvm(i,j)=real(amvd(i,j)) do 28 ia=1,nat xi=r(ia,1)-c1(1) yi=r(ia,2)-c1(2) zi=r(ia,3)-c1(3) ri(ia,1)=amvm(1,1)*xi+amvm(2,1)*yi+amvm(3,1)*zi ri(ia,2)=amvm(1,2)*xi+amvm(2,2)*yi+amvm(3,2)*zi 28 ri(ia,3)=amvm(1,3)*xi+amvm(2,3)*yi+amvm(3,3)*zi return end subroutine iniind(ind,nat) implicit none integer*4 nat,ind(8,nat),i,o do 1 o=1,8 do 1 i=1,nat 1 ind(o,i)=0 return end subroutine sets(s10,i,istart) implicit none character*10 s10 integer*4 i,istart write(s10,5002)i 5002 format(i10) do 504 istart=1,len(s10) 504 if(s10(istart:istart).ne.' ')goto 505 505 return end subroutine chkxx(amv,rc) implicit none integer ix,jx real*8 amv(3,3),rc(3) if(amv(3,3)*(amv(1,1)*amv(2,2)-amv(1,2)*amv(2,1)).lt.0.0d0)then write(6,*)'3 3:',amv(3,3) write(6,*)'1x2:',amv(1,1)*amv(2,2)-amv(1,2)*amv(2,1) amv(1,3)=-amv(1,3) amv(2,3)=-amv(2,3) amv(3,3)=-amv(3,3) write(6,*)'chirality fixed' endif write(6,705)(rc(ix),ix=1,3),((amv(jx,ix),ix=1,3),jx=1,3) 705 format(/,' Moments of inertia and the matrix:',/,/, 13f12.3,/,/,3(3f12.3,/)) return end