PROGRAM CONTACTC IMPLICIT none logical lex character*4 key real*4 bonding(88),r,d1,d2,xi,xj,yi,yj,zi,zj,rt,rb,pm, 1qm,tol,cm,c1,q1,p1,tol2,ri,rj,amvm,rfac real*8 am(3,3),amv(3,3),rc(3),etem(3) integer*4 ik,iout,i,j,ic,nat,natc,nc,alist,natm,ia,ii,ja, 1clist,imol,ix,jx,jj,nmol,jmol,istart,IERR,ico,ind,itry,iie logical li allocatable ::ik(:),r(:,:),iout(:),alist(:,:),natm(:), 1clist(:,:),natc(:),cm(:,:),pm(:,:),qm(:,:,:),c1(:,:),p1(:,:), 2q1(:,:,:),amvm(:,:,:),ri(:,:),rj(:,:),ico(:),ind(:),li(:) 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 incluster,lwr character*10 s10 c WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X, selects',/, 1 ' molecules with neighboring molecules',/ 1 ' Input: CRYSTAL.X',/, 3 ' Output: cluster geometries, CCT.INP',/,/) d1=4.0 tol=0.5 tol2=0.1 lwr=.true. rfac=1.1 inquire(file='CONTACTC.OPT',exist=lex) if(lex)then open(2,file='CONTACTC.OPT') 100 read(2,2000,end=2200,err=2200)key 2000 format(a4) if(key.eq.'DIST')read(2,*)d1 if(key.eq.'TOLE')read(2,*)tol if(key.eq.'TOL2')read(2,*)tol2 if(key.eq.'WRIT')read(2,*)lwr if(key.eq.'RFAC')read(2,*)lwr goto 100 2200 close(2) endif d2=d1**2 write(6,606)d1,tol 606 format(' Distance limit: ',F10.1,' A',/, 1 ' RMS tolerance : ',F10.1) open(2,file='CRYSTAL.X') read(2,*) read(2,*)nat allocate(ik(nat),r(nat,3),iout(nat),ri(nat,3),rj(nat,3),ind(nat), 1li(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 CRYSTAL.X' 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),clist(nmol,nat),natc(nmol)) do 8 imol=1,nmol natm(imol)=0 8 natc(imol)=0 do 1 i=1,nat natm(iout(i))=natm(iout(i))+1 natc(iout(i))=natc(iout(i))+1 alist(iout(i),natm(iout(i)))=i 1 clist(iout(i),natc(iout(i)))=i if(lwr)then do 9 imol=1,nmol write(6,*)' molecule ',imol,': ',natm(imol),' atoms' 9 write(6,603)(alist(imol,jj),jj=1,natm(imol)) 603 format(15i5) endif c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule c clist .. list of atoms in a cluster write(6,*) allocate(cm(nmol,3),pm(nmol,5),qm(nmol,3,3),amvm(nmol,3,3), 1c1(nmol,3),p1(nmol,5),q1(nmol,3,3),ico(nmol)) c make bigger clusters: c 400000000000000000000000000000000000000000000000000000000000 do 4000 imol=1,nmol write(6,*)' Molecule ',imol,':' do 12 j=1,nat if(iout(j).ne.imol)then xi=r(j,1) yi=r(j,2) zi=r(j,3) do 121 i=1,nat if(iout(i).eq.imol)then xj=(r(i,1)-xi)**2 if(xj.lt.d2)then yj=xj+(r(i,2)-yi)**2 if(yj.lt.d2)then if(yj+(r(i,3)-zi)**2.lt.d2)then c atom j is not in imol, but is close, add it to cluster: natc(imol)=natc(imol)+1 clist(imol,natc(imol))=j goto 12 endif endif endif endif 121 continue endif 12 continue write(6,*)natc(imol)-natm(imol),'close atoms' if(lwr)then write(6,603)(clist(imol,jj),jj=natm(imol)+1,natc(imol)) endif c keep atoms bound in a molecule close to central molecule: do 40 i=1,nat 40 li(i)=incluster(i,imol,natc,clist,nmol,nat) nc=0 992 ic=0 do 14 i=1,nat if(.not.li(i))then c i is not in cluster imol: xi=r(i,1) yi=r(i,2) zi=r(i,3) rb=bonding(ik(i)+1) do 7 j=1,nat if(li(j).and.j.ne.i)then c j is in cluster imol: rt=(rb+bonding(ik(j)+1))**2 xj=(r(j,1)-xi)**2 if(xj.lt.rt)then yj=xj+(r(j,2)-yi)**2 if(yj.lt.rt)then if(yj+(r(j,3)-zi)**2.lt.rt)then c atom i is bound to j, which is in=take i, too: natc(imol)=natc(imol)+1 clist(imol,natc(imol))=i li(i)=.true. ic=ic+1 nc=nc+1 goto 14 endif endif endif endif 7 continue endif 14 continue if(ic.ne.0)goto 992 4000 write(6,*)nc,' additionally bound atoms' c 400000000000000000000000000000000000000000000000000000000000 c imol c order according to the size c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 5 imol=1,nmol do 5 jmol=imol+1,nmol if(natc(imol).lt.natc(jmol))then call isw(natc(imol),natc(jmol)) do 23 jj=1,max(natc(imol),natc(jmol)) call isw(alist(imol,jj),alist(jmol,jj)) 23 call isw(clist(imol,jj),clist(jmol,jj)) endif 5 continue c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write(6,*)'clusters reordered' do 36 imol=1,nmol write(6,*)' Cluster ',imol 36 write(6,6009)(clist(imol,jj),jj=1,natc(imol)) 6009 format(16i5) c calculate moments c 444444444444444444444444444444444444444444444444444444444444 do 4 imol=1,nmol call sets(s10,imol,istart) call wrx('c.'//s10(istart:len(s10))//'.x',imol,nat,r,natc, 1clist,nmol,ik) write(6,*)' Molecule ',imol,':' c cluster center and dipoles: write(6,*)'cluster:' call cd(cm,pm,qm,natc,imol,clist,nmol,nat,r,ik) write(6,*)'molecule:' c molecular center: call cd(c1,p1,q1,natm,imol,alist,nmol,nat,r,ik) c c diagonalize moment of inertia ans save transformation matrix: do 27 ix=1,3 do 27 jx=1,3 27 am(ix,jx)=dble(q1(imol,ix,jx)) IERR=0 CALL TRED12(3,am,amv,rc,2,IERR,etem) c uniform sign(phase) of eigenvectors: c do 28 ix=1,3 c if(amv(1,ix).lt.0.0d0)then c do 281 jx=1,3 c81 amv(jx,ix)=-amv(jx,ix) c endif c8 continue c check that we have a right-handed system 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,/)) do 4 ix=1,3 do 4 jx=1,3 4 amvm(imol,ix,jx)=real(amv(ix,jx)) c 444444444444444444444444444444444444444444444444444444444444 c imol c c write final atomic map: open(19,file='CCT.TEM') write(19,901)nmol 901 format('LNAMES',/,'T',/,'POLYMER',/,I9) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 26 imol=1,nmol ico(imol)=imol c shift cluster to the central molecule and transform to inertia axes: call shm(nat,nmol,imol,amvm,natc,clist,r,ri,c1) c write this to disk: call sets(s10,imol,istart) call wrx('t.'//s10(istart:len(s10))//'.x', 1imol,nat,ri,natc,clist,nmol,ik) c indicate enantiomer: iie=0 c c look if imol is contained in some jmol: do 29 jmol=1,imol-1 do 291 itry=0,7 c transform jmol to inertia axes: call shm(nat,nmol,jmol,amvm,natc,clist,r,rj,c1) c rotate x y z inversion invx invy invz: call rr(itry,nat,nmol,jmol,natc,clist,rj) c find overlapped atoms: call iniind(ind,nat) do 33 ii=1,natc(imol) ia=clist(imol,ii) xi=ri(ia,1) yi=ri(ia,2) zi=ri(ia,3) do 34 jj=1,natc(jmol) ja=clist(jmol,jj) xj=rj(ja,1)-xi yj=rj(ja,2)-yi zj=rj(ja,3)-zi if(ik(ia).eq.ik(ja).and.xj**2+yj**2+zj**2.lt.tol2)then ind(ia)=jj goto 33 endif 34 continue c atom i does not have analogue in jmol, try next orientation orjmol: goto 291 33 continue c all atoms had analogues, fragment imol is part of fragment jmol: if(itry.ge.4)then c for enantiomers - inform only for the first time: if(iie.eq.0)then iie=iie+1 write(6,706)imol,jmol 706 format(' INFO: ',i4,' IS ENANTIOMER OF ',i4) endif goto 29 else c for true isomers - take it: ico(imol)=jmol goto 333 endif 291 continue 29 continue c jmol c 333 if(ico(imol).eq.imol)then c atom i does not have analogue in any jmol, try x y z C2 rotations write(6,701)imol 701 format(i5,' is not contained anywhere') jmol=ico(imol) call sets(s10,jmol,istart) call wrx('i.'//s10(istart:len(s10))//'.x',imol,nat,r,natc, 1 clist,nmol,ik) call iniind(ind,nat) do 35 ii=1,natc(imol) 35 ind(clist(imol,ii))=ii else write(6,702)imol,ico(imol),itry 702 format(i5,' is contained in ',i5,' trial ',i1) 704 if(ico(ico(imol)).ne.ico(imol))then write(6,703)ico(ico(imol)) 703 format(' ... which is contained in ',i5) call report(' .. which is error in program logic!') endif jmol=ico(imol) endif write(19,*)imol call sets(s10,jmol,istart) write(19,'(a)')'i.'//s10(istart:len(s10)) write(19,902)(i-(i/1000)*1000,i=1,nat) 902 format(20i3) do 261 i=1,nat write(19,9021)ind(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 26 continue c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 37 imol=1,nmol 37 write(6,609)imol,natc(imol),ico(imol),pm(imol,4),pm(imol,5) 609 format(3i10,2f10.1) c write(6,*)'CCT.TEM written' close(19) stop end subroutine report(s) character*(*) s write(6,*)s stop end function incluster(j,imol,natc,clist,nmol,nat) c is j in the cluster imol?: implicit none logical incluster integer*4 j,jj,imol,natc(*),nmol,nat,clist(nmol,nat) do 21 jj=1,natc(imol) if(clist(imol,jj).eq.j)then incluster=.true. return endif 21 continue incluster=.false. return end subroutine wrx(s,imol,nat,r,natc,clist,nmol,ik) implicit none integer*4 imol,nat,i,natc(*),ix,nmol,ik(*),clist(nmol,nat) real*4 r(nat,3) character*(*)s open(9,file=s) write(9,*)' cluster ',imol write(9,*)natc(imol) do 1 i=1,natc(imol) 1 write(9,90)ik(clist(imol,i)),(r(clist(imol,i),ix),ix=1,3) 90 format(i5,3f12.6,' 0 0 0 0 0 0 0 0.0') close(9) write(6,*)s//' written' end subroutine cd(cm,pm,qm,natc,imol,clist,nmol,nat,r,ik) 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)+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 isw(i,j) integer*4 i,j,k k=i i=j j=k return end subroutine sw(i,j) real*4 i,j,k k=i i=j j=k 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, 1c1(nmol,3) do 32 ix=1,3 do 32 jx=1,3 32 tr(ix,jx)=0.0d0 c xx=0: identity if(xx.eq.0)return c xx=1: C2x-rotation if(xx.eq.1)then tr(1,1)= 1.0d0 tr(2,2)=-1.0d0 tr(3,3)=-1.0d0 endif c xx=2: C2y-rotation if(xx.eq.2)then tr(1,1)=-1.0d0 tr(2,2)= 1.0d0 tr(3,3)=-1.0d0 endif c xx=3: C2z-rotation if(xx.eq.3)then tr(1,1)=-1.0d0 tr(2,2)=-1.0d0 tr(3,3)= 1.0d0 endif c xx=4: inversion if(xx.eq.4)then tr(1,1)=-1.0d0 tr(2,2)=-1.0d0 tr(3,3)=-1.0d0 endif c xx=5: inversion and C2x rotation if(xx.eq.5)then tr(1,1)=-1.0d0 tr(2,2)= 1.0d0 tr(3,3)= 1.0d0 endif c xx=6: inversion and C2y rotation if(xx.eq.6)then tr(1,1)= 1.0d0 tr(2,2)=-1.0d0 tr(3,3)= 1.0d0 endif c xx=7: inversion and C2z rotation if(xx.eq.7)then tr(1,1)= 1.0d0 tr(2,2)= 1.0d0 tr(3,3)=-1.0d0 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 iniind(ind,nat) implicit none integer*4 nat,ind(nat),i do 1 i=1,nat 1 ind(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