PROGRAM CONTACTN IMPLICIT none real*4 bonding(88),r,d1,d2,xi,xj,yi,yj,zi,zj,rt,rb,pm, 1qm,tol,cm,tol2,ri,rj,amvm,px,py,pz,btol real*8 am(3,3),amv(3,3),rc(3),etem(3) integer*4 ik,iout,i,j,ic,nat,nc,alist,natm,ia,ii,ja,np0, 1imol,ix,jx,jj,nmol,jmol,istart,IERR,ico,ind,itry,iie,lo,lot, 1nat1,iu,iof,ice,natx,iccub,indn,nt,nn,NM,STATE,si,lm,lmt,icen, 1icube,iex,ige,jge,im,natt,nexc,nv allocatable ::ik(:),r(:,:),iout(:),alist(:,:),natm(:), 1cm(:,:),pm(:,:),qm(:,:,:),STATE(:),si(:),lm(:),lmt(:),lot(:,:), 2amvm(:,:,:),ri(:,:),rj(:,:),ico(:),ind(:),indn(:),lo(:,:) 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,lall,notisit character*10 s10 c pair variables: real*4 dmin integer*4 np integer*4,allocatable::clist(:,:),natc(:) c WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X, selects',/, 1 ' clusters of neighboring molecules',/ 1 ' Input: CRYSTAL.X',/, 1 ' CELL.TXT',/, 3 ' Output: cluster geometries, CCT.TEM',/,/) call readopt(d1,d2,tol,tol2,px,py,pz,lwr,btol,lall,NM) 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), 1indn(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' ice=iccub(nat,r) nat1=nat/27 iof=nat1*(ice-1) write(6,602) 602 format(' Finding molecules') c find free atom, and connected atoms, in the first cube actually: nmol=1 iout(1)=1 991 nc=1 99 ic=0 do 13 i=1,nat1 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)*btol do 2 j=1,nat1 if(iout(j).eq.nmol)then rt=(rb+btol*bonding(ik(j)+1))**2 xj=r(j,1)-xi yj=r(j,2)-yi zj=r(j,3)-zi if(xj**2+yj**2+zj**2.lt.rt)then c atom i is bound to j, which is in, take i, too: iout(i)=nmol ic=ic+1 nc=nc+1 goto 13 endif endif 2 continue endif 13 continue c if still adding to nmol, continue: if(ic.ne.0)goto 99 do 131 i=1,nat1 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,606)nmol 606 format(i5,' molecules found, atom lists:') c alist .. atom list related to the first cube!!! allocate(alist(nmol,nat),natm(nmol)) do 8 imol=1,nmol 8 natm(imol)=0 do 1 i=1,nat1 natm( iout(i))=natm(iout(i))+1 1 alist(iout(i),natm(iout(i)))=i if(lwr)then do 9 imol=1,nmol write(6,604)imol 604 Format(i4,':',$) 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 write(6,*) write(6,608)NM 608 format(/,' Looking for close molecular clusters, NM =',i2,/) allocate(STATE(nat),si(NM),lm(nat),lmt(nat),lo(nat,NM)) np0=nat lo=0 natx=0 nc=0 np=0 c imol in ice do 4000 imol=1,nmol c ige... generalized number ige=imol+nmol*(ice-1) c molecules in the vicinity nv=0 do 4001 iu=1,27 ii=nat1*(iu-1) do 4001 jmol=1,nmol c jge... generalized number jge=jmol+nmol*(iu-1) dmin=1.0e9 do 12 i=1,natm(imol) ia=alist(imol,i)+iof xi=r(ia,1) yi=r(ia,2) zi=r(ia,3) do 12 j=1,natm(jmol) ja=alist(jmol,j)+ii xj=r(ja,1)-xi yj=r(ja,2)-yi zj=r(ja,3)-zi 12 if(xj**2+yj**2+zj**2.lt.dmin)dmin=xj**2+yj**2+zj**2 if(dmin.lt.d2)then nv=nv+1 lm(nv)=jge endif 4001 continue if(nv+1.lt.NM)call report('Not enough molecules around') if(lwr)then write(6,609)imol,ige,nv 609 format(2i6,' molecule,',i6,' molecules around (total inices):') write(6,612)(lm(j),j=1,nv) 612 format(10i5) endif c make MN-molecule clusters from nv molecules c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx Nexc=NM do 4100 iex=1,Nexc 4100 si(iex)=1 5 DO 1966 j=1,nv 1966 STATE(j)=0 do 4009 iex=1,Nexc 4009 STATE(si(iex))=STATE(si(iex))+1 c check that the state is allowed: c # of centers=# of excitations: icen=0 do 11 j=1,nv 11 if(STATE(j).gt.0)icen=icen+1 if(icen.ne.Nexc)goto 1212 c must contain imol~ige ic=0 do 4010 j=1,NM 4010 if(lm(si(j)).eq.ige)ic=1 if(ic.eq.0)goto 1212 c suggested cluster: do 4002 j=1,NM 4002 lmt(j)=lm(si(j)) call sort(lmt,NM) c is it in the list, if not, record: if(notisit(lo,lmt,NM,nat,np))then np=np+1 if(np.gt.np0)then allocate(lot(np0,NM)) lot=lo deallocate(lo) allocate(lo(np,NM)) np0=np do 4006 ii=1,np-1 do 4006 j=1,NM 4006 lo(ii,j)=lot(ii,j) deallocate(lot) endif do 4003 j=1,NM 4003 lo(np,j)=lmt(j) endif c c find index to be changed 1212 do 4007 ic=Nexc,1,-1 4007 if(si(ic).lt.nv)goto 129 goto 3333 129 do 10 j=ic+1,Nexc 10 si(j)=si(ic)+1 si(ic)=si(ic)+1 c goto 5 3333 continue c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx do 4008 i=1,np natt=0 do 4005 j=1,NM icube=(lo(i,j)-1)/nmol+1 im=lo(i,j)-(icube-1)*nmol 4005 natt=natt+natm(im) 4008 if(natx.lt.natt)natx=natt 4000 continue c 400000000000000000000000000000000000000000000000000000000000 write(6,*)np,' clusters, natx = ',natx if(lwr)then do 4004 i=1,np write(6,605)i 605 format(i6,':',$) 4004 write(6,611)(lo(i,j),j=1,NM) 611 format(10i5) write(6,*) endif allocate(natc(np),cm(np,3),pm(np,5), 1qm(np,3,3),amvm(np,3,3),ico(np),clist(np,natx)) write(6,*)'clist allocated' c clist ... atom list in pair, absolute atom numbers: do 38 i=1,np natc(i)=0 do 38 j=1,NM icube=(lo(i,j)-1)/nmol+1 imol=lo(i,j)-(icube-1)*nmol do 38 ia=1,natm(imol) natc(i)=natc(i)+1 38 clist(i,natc(i))=alist(imol,ia)+nat1*(icube-1) c open(9,file='29t.x') c write(9,*)'test' c write(9,*)natc(29) c do 9009 ia=1,natc(29) c ii=clist(29,ia) c9009 write(9,9019)ik(ii),(r(ii,j),j=1,3) c9019 format(i6,3f12.6) c close(9) c stop c c calculate moments c 444444444444444444444444444444444444444444444444444444444444 do 4 imol=1,np write(6,*)'Cluster:',imol call sets(s10,imol,istart) call wrx('p.'//s10(istart:len(s10))//'.x',imol,nat,r,natc, 1clist,np,natx,ik) c cluster center and dipoles: call cd(cm,pm,qm,natc,imol,clist,np,natx,nat,r,ik) c diagonalize moment of inertia ans save transformation matrix: do 27 ix=1,3 do 27 jx=1,3 27 am(ix,jx)=dble(qm(imol,ix,jx)) IERR=0 CALL TRED12(3,am,amv,rc,2,IERR,etem) 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') if(lall)then write(19,903) write(6,903) 903 format('LNAMES',/,'T',/,'POLYMER',/,'set manually') else write(19,901)np 901 format('LNAMES',/,'T',/,'POLYMER',/,I9) endif c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc nn=0 do 26 imol=1,np ico(imol)=imol c shift cluster to the centre and transform to inertia axes: call shm(nat,np,natx,imol,amvm,natc,clist,r,ri,cm) c write this to disk: 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,np,natx,jmol,amvm,natc,clist,r,rj,cm) c rotate x y z inversion invx invy invz: call rr(itry,nat,np,natx,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) if(ik(ia).eq.ik(ja))then xj=rj(ja,1)-xi yj=rj(ja,2)-yi zj=rj(ja,3)-zi if(xj**2+yj**2+zj**2.lt.tol2)then ind(ia)=jj goto 33 endif 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)=ico(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,np,natx,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) 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 c c center cube: nn=nn+1 write(19,*)nn 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 c c other cubes: if(lall)then do 39 ic=1,27 if(ic.ne.ice)then call iniind(indn,nat) nt=0 do 41 i=1,nat ii=i+nat1*(ic-ice) if(ind(i).ne.0.and.ii.gt.0.and.ii.le.nat)then indn(ii)=ind(i) nt=nt+1 endif 41 continue if(nt.eq.natc(imol))then nn=nn+1 write(19,*)nn write(19,'(a)')'i.'//s10(istart:len(s10)) write(19,902)(i-(i/1000)*1000,i=1,nat) do 262 i=1,nat write(19,9021)indn(i) 262 if(20*(i/20).eq.i)write(19,9022)i if(20*(nat/20).ne.nat)write(19,9022)nat endif endif 39 continue endif 26 continue c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(nn.ne.np)write(6,*)nn,' overlaps' close(19) write(6,*)'CCT.TEM written' end subroutine report(s) character*(*) s write(6,*)s stop end subroutine wrx(s,imol,nat,r,natc,clist,nmol,natx,ik) implicit none integer*4 imol,nat,i,natc(*),ix,nmol,natx,ik(*),clist(nmol,natx) 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,natx,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,natx,nat,natc(*),imol,clist(nmol,natx), 1ix,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 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,natx,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,natx,clist(nmol,natx), 1ia,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,natx,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,natx,clist(nmol,natx),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 subroutine readopt(d1,d2,tol,tol2,px,py,pz,lwr,btol,lall,NM) implicit none real*4 d1,d2,tol,tol2,px,py,pz,btol logical lwr,lex,lall integer*4 NM character*4 key c bonding cutoff multiplied by btol, to accomodate for MD distorted c bonds btol=1.1 d1=4.0 tol=0.5 tol2=0.1 lwr=.true. c pairs to all atoms in the CRYSTAL, not only central cell: lall=.false. px=0. py=0. pz=0. c number of molecules in the clusters: NM=3 inquire(file='CELL.TXT',exist=lex) if(lex)then open(2,file='CELL.TXT') read(2,*)px read(2,*)py,py read(2,*)pz,pz,pz close(2) endif inquire(file='CONTACTN.OPT',exist=lex) if(lex)then open(2,file='CONTACTN.OPT') 100 read(2,2000,end=2200,err=2200)key 2000 format(a4) if(key.eq.'CELL')read(2,*)px,py,pz 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.'BTOL')read(2,*)btol if(key.eq.'LALL')read(2,*)lall if(key.eq.'MOLE')read(2,*)NM goto 100 2200 close(2) endif if(px.lt.0.1)call report('PX must be set') if(py.lt.0.1)call report('PY must be set') if(pz.lt.0.1)call report('PZ must be set') d2=d1**2 write(6,606)d1,tol 606 format(' Distance limit: ',F10.1,' A',/, 1 ' RMS tolerance : ',F10.1) return end function iccub(nat,r) implicit none integer*4 nat,nat1,ix,iu,ii,iccub,ice,i real*4 cet(3),cei(27,3),r(nat,3),xi,yi,zi if(mod(nat,27).ne.0)call report('mod nat 27 <> 0') nat1=nat/27 c determine central cube do 40 ix=1,3 cet(ix)=0.0 do 42 iu=1,27 cei(iu,ix)=0.0 ii=nat1*(iu-1) do 41 i=ii+1,ii+nat1 cet( ix)=cet( ix)+r(i,ix) 41 cei(iu,ix)=cei(iu,ix)+r(i,ix) 42 cei(iu,ix)=cei(iu,ix)/real(nat1) 40 cet(ix)=cet(ix)/real(nat) write(6,600)cet 600 format(/,' CRYSTAL.X center (A):',7x,3f7.2,/,' Cube centers (A):') write(6,601)(iu,(cei(iu,ix),ix=1,3),iu=1,27) 601 format(3(i3,':',3f7.2)) ice=0 do 43 iu=1,27 xi=cei(iu,1)-cet(1) yi=cei(iu,2)-cet(2) zi=cei(iu,3)-cet(3) 43 if(sqrt(xi**2+yi**2+zi**2).lt.0.1)ice=iu if(ice.eq.0)then call report('Central cube not found') else write(6,602)ice 602 format(' Central cube:',i3,', finding molecules') iccub=ice endif return end subroutine sort(lmt,NM) implicit none integer*4 lmt(*),NM,i,j,t do 1 i=1,NM do 1 j=i+1,NM if(lmt(j).lt.lmt(i))then t=lmt(i) lmt(i)=lmt(j) lmt(j)=t endif 1 continue return end function notisit(lo,lmt,NM,nat,np) implicit none logical notisit integer*4 lmt(*),NM,nat,lo(nat,NM),np,i,j do 1 i=1,np do 11 j=1,NM 11 if(lo(i,j).ne.lmt(j))goto 1 c lmt is equal to i_th element of lo: notisit=.false. return 1 continue c lmt has not been found anywhere: notisit=.true. return end