PROGRAM CONTACTP 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,an,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, 1imol,ix,jx,jj,nmol,jmol,istart,IERR,ico,ind,itry,iie, 1nat1,iu,iof,ice,natx,nave,iccub allocatable ::ik(:),r(:,:),iout(:),alist(:,:),natm(:), 1cm(:,:),pm(:,:),qm(:,:,:),nave(:), 2amvm(:,:,:),ri(:,:),rj(:,:),ico(:),ind(:) 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,lmol character*10 s10 c pair variables: real*4 dmin integer*4 np integer*4,allocatable::m1(:),m2(:),mc(:),clist(:,:),natc(:) c WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X, selects',/, 1 ' pairs 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,lmol) 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)) 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) if(lmol)then write(6,608) 608 format(' Reading molecules from MOL.LST') open(9,file='MOL.LST',status='old') read(9,*)nmol allocate(alist(nmol,nat),natm(nmol)) do 14 i=1,nmol 14 read(9,*)natm(i),(alist(i,j),j=1,natm(i)) close(9) else 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 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 endif write(6,606)nmol 606 format(i5,' molecules') if(lwr)then write(6,*)'Atom list:' 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) write(6,*) endif c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule write(6,*)' Looking for close molecular pairs' natx=0 nc=0 np=0 allocate(nave(nmol)) do 45 i=1,nmol 45 nave(i)=0 c first cycle:determine number, second:record do 4000 ic=1,2 if(ic.eq.2)then allocate(m1(np),m2(np),mc(np),natc(np),cm(np,3),pm(np,5), 1 qm(np,3,3),amvm(np,3,3),ico(np),clist(np,natx)) endif np=0 c jmol in ice, imol in all cubes do 4000 iu=1,27 ii=nat1*(iu-1) do 4000 imol=1,nmol do 4000 jmol=1,nmol if(iu.ne.ice.or.jmol.gt.imol)then dmin=1.0e9 do 12 i=1,natm(imol) ia=alist(imol,i)+ii xi=r(ia,1) yi=r(ia,2) zi=r(ia,3) do 12 j=1,natm(jmol) ja=alist(jmol,j)+iof 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 np=np+1 if(ic.eq.1)then natx=max(natx,natm(imol)+natm(jmol)) else m1(np)=jmol m2(np)=imol mc(np)=iu nave(jmol)=nave(jmol)+1 if(iu.eq.ice)nc=nc+1 endif endif endif 4000 continue c 400000000000000000000000000000000000000000000000000000000000 write(6,*)np,' pairs,',nc,' within cell',ice write(6,605)(i,m1(i),ice,m2(i),mc(i),i=1,np) 605 format(4(i5,i3,'(',i2,')',i3,'(',i2,');')) write(6,*) an=0 do 46 i=1,nmol 46 an=an+real(nave(i)) an=an/real(nmol) write(6,607)an 607 format(f12.1,' molecules sourounding each on overage') c clist ... atom list in pair, absolute atom numbers: do 38 i=1,np natc(i)=0 do 381 ia=1,natm(m1(i)) natc(i)=natc(i)+1 381 clist(i,natc(i))=alist(m1(i),ia)+iof do 38 ia=1,natm(m2(i)) natc(i)=natc(i)+1 38 clist(i,natc(i))=alist(m2(i),ia)+nat1*(mc(i)-1) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c calculate moments c 444444444444444444444444444444444444444444444444444444444444 do 4 imol=1,np write(6,*)'Pair:',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') write(19,901)np 901 format('LNAMES',/,'T',/,'POLYMER',/,I9) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 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,np 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 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,lmol) implicit none real*4 d1,d2,tol,tol2,px,py,pz,btol logical lwr,lex,lmol 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. px=0. py=0. pz=0. lmol=.false. 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='CONTACTS.OPT',exist=lex) if(lex)then open(2,file='CONTACTS.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.'LMOL')read(2,*)lmol c LMOL ... determine molecule from MOL.LST 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