' ******************************************************* ' *** Program: MXD-STEREO.BAS *** ' *** GRAPHIC REPRESENTATION OF MD-RESULTS *** ' *** USING FILE05.DAT and FILE07.DAT *** ' *** *** ' *** First XY-plotter version (1993-12-27 by Kats) *** ' *** LIPS3 (1994-02-08 by Kats) *** ' *** MXDTRICL (1995-10-17 by Kats) *** ' *** Paint (1996-02-07 by Kats) *** ' *** Cross and Color (1997-08-25 by Kats) *** ' *** Di-atomic molecule (1997-10-20 by Kats) *** ' *** Radius and color (2000-08-10 by Kats) *** ' *** width lcx,lcy (2000-10-17 by kats) *** ' *** Free BASIC version (2005-08-14 by Kats) *** ' ******************************************************* deflng I-N defdbl a-h,o-z ' ' screen 19 : xscreen= 800: yscreen=600 screen 20 : xscreen=1024: yscreen=768 ' LCX=xscreen/8 : LCY=yscreen/16 LPX=LCX*8-1 : LPY=LCY*16-1 ' dim XYORGN(3),XYZRANGE(2,3) LEL=10: dim AL(6), ZI(LEL), NET(LEL), ATOM$(LEL), IONS(2,LEL) dim BL(6), RI1(LEL), ICC(LEL), NION(LEL), SPRES(3) dim RI2(LEL), H(3,3), BOX(3,8),LBOX(3,8), alc(6),nbox(3) LNI=62654: dim XLL(8), KBOX(3,8), IELE(LNI), IBOND(6,LNI) dim XLR(8), ISEQ(LNI), DLIM2(LEL), DLIM2I(LEL) dim YLS(8), JSEQ(LNI), PP(3), P(3,LNI), IP0(LNI) dim IX$(3), IDIATOM(LNI),XYZ(3,LNI),P0(3,LNI) LSQ=60: dim XL(LNI), XR(LNI), YS(LNI), ZS(LNI), IPZ(LNI) L10=7000 : dim JJQ(LSQ),JSQ(L10,LSQ) dim title$,a$,a0$,a1$,a2$,a3$,a4$.a5$,a6$,a7$,a8$ dim temp,dnsty,rcuta,rcuta2,dlimio,dlim2io,dlim2iio,ziio,dx,dy,dz,r,r2,dd,z,x1,x2,y dim pxi,pyi,pzi dim ncompo,ntion,jcol,io,jo,ix,iy,ixl,ixr,ic,i,j ' IX$(1) = "a": XYZRANGE(1,1)=0.0: XYZRANGE(2,1)=1.0 IX$(2) = "b": XYZRANGE(1,2)=0.0: XYZRANGE(2,2)=1.0 IX$(3) = "c": XYZRANGE(1,3)=0.0: XYZRANGE(2,3)=1.0 ' F05F07: cls 0 : line (0,0)-(LPX,LPY),15,bf,15 : color 3,15 print lpx,lpy print print "Program : MXD-STRO is running" print print "Type the directory name of MD-calculated data " print " (ex. b:\, a:\dirname\ ) " print " ([return] for current dir. ) " print " (q to quit from mxd-stro ) "; ' input DR$ : color 0 if DR$ = "q" or DR$ = "Q" then stop: end if DR$ <> "" and right$(DR$,1) <> "\" then DR$ = DR$ + "\" ' gosub FILE05 : DT$ = DATe$ gosub FILE07 : TM$ = time$ ' gosub ATOMS ' IXX=1 : IXY=2 : IXZ=3 if AL(2)>AL(1) then IXX=2 : IXY=1 if AL(3)>AL(IXX) then IXZ=IXX : IXX=3 if AL(IXY)0 and IA<=NCOMPO then locate 12,15 : print ia;":";ATOM$(IA);"is changed" locate 13,15 : input "radius";RI1(IA) locate 14,15 : input "Color ";ICC(IA) end if if A$="Q" or A$="q" then a$="" : return goto ATOMSetting end ' ' STRUCTURE: '====================================================== STRUCUTRE for IO = 1 to LEL if NION(IO) <> 0 then for J = IONS(1,IO) to IONS(2,IO) IELE(J) = IO next J end if next IO ' gosub CALC_POSITIONS ' '------------------------------------------------Search Bonds print print "Searching bonds at "; : LY = csrlin LX = pos() for IO = 1 to LEL if NION(IO) < 1 then goto ST7 if NET(IO) < .1 and RCUTA<0.1 then goto ST7 DLIM2IO = DLIM2(IO): DLIMIO=sqr(DLIM2IO) DLIM2IIO = DLIM2I(IO): ZIIO = ZI(IO) for I = IONS(1, IO) to IONS(2,IO) IDIATOM(I) = 0 PXI = P(1,I) PYI = P(2,I) PZI = P(3,I): NBI = 0 for K = 1 to 6: IBOND(K,I) = 0: next K if net(io)<0.1 then goto st4 for J = 1 to NTION if I = J then goto ST3 JO = IELE(J) 'if ZIIO * ZI(JO) <= 0 then DX = PXI - P(1,J): if abs(DX)>DLIMIO then goto ST3 DY = PYI - P(2,J): if abs(DY)>DLIMIO then goto ST3 DZ = PZI - P(3,J): if abs(DZ)>DLIMIO then goto ST3 R2 = DX * DX + DY * DY + DZ * DZ if R2 <= DLIM2IO then NBI = NBI + 1: if NBI > 6 then goto ST4 IBOND(NBI,I) = J if DLIM2IIO > .1 then if R2 > DLIM2IIO then IBOND(NBI,I)=-J end if end if 'end if ST3: next J ST4: if RCUTA<0.1 goto ST6 'print RCUTA,IONS(1,IO), IONS(1,IO) RCUTA2 = RCUTA for J=IONS(1,IO) to IONS(2,IO) : ' ----- Diatomic molecule if I = J then goto ST5 DX = PXI - P(1,J): if abs(DX)>RCUTA then goto ST5 DY = PYI - P(2,J): if abs(DY)>RCUTA then goto ST5 DZ = PZI - P(3,J): if abs(DZ)>RCUTA then goto ST5 R2 = DX*DX + DY*DY + DZ*DZ if R2 <= RCUTA2 then 'print I,J,R2 IDIATOM(I) = J: goto ST6 end if ST5: next J ST6: locate LX,LY print using "##### (#) "; I; NBI; next I ST7: next IO ' return end ' ' CALC_POSITIONS: '============================================ CALC COORDINATES AL4S = sqr(1-ALC(IXX+3)^2) H(IXX,IXZ) = 0 H(IXY,IXZ) = 0 H(IXZ,IXZ) = AL(IXZ) H(IXX,IXX) = AL(IXX) * AL4S H(IXY,IXX) = 0 H(IXZ,IXX) = AL(IXX) * ALC(IXY+3) H(IXX,IXY) = AL(IXY) * ALC(IXZ+3) * AL4S H(IXY,IXY) = AL(IXY) * sqr(1-ALC(IXX+3)^2 - (ALC(IXZ+3)*AL4S)^2) H(IXZ,IXY) = AL(IXY) * ALC(IXX+3) ' for I = 1 to NTION '------------------------- Calc Atom positions XX = P0(1,I) - .5 YY = P0(2,I) - .5 ZZ = P0(3,I) - .5 P(1,I) = H(1,1)*XX + H(1,2)*YY + H(1,3)*ZZ P(2,I) = H(2,1)*XX + H(2,2)*YY + H(2,3)*ZZ P(3,I) = H(3,1)*XX + H(3,2)*YY + H(3,3)*ZZ next I N=0 '---------------------------------------- Calc Lattice Points for I=0 to 1: for J=0 to 1: for K = 0 to 1 N=N+1 LBOX(1,N)=I: LBOX(2,N)=J: LBOX(3,N)=K XX= I-.5: YY= J-.5: ZZ=K-.5 BOX(1,N) = H(1,1)*XX + H(1,2)*YY + H(1,3) * ZZ BOX(2,N) = H(2,1)*XX + H(2,2)*YY + H(2,3) * ZZ BOX(3,N) = H(3,1)*XX + H(3,2)*YY + H(3,3) * ZZ next K: next J: next I return end ' ' SEQUENCE: '=========================================================== Sequence APMAX = 1 / ((AL(1)^2 + AL(2)^2 + AL(3)^2)*2) for I=1 to LSQ : JJQ(I)=0: next I for I=1 to NTION : J=IP0(I)*APMAX*LSQ+1 : 'print ip0(i); if J<1 then J=1 if J>LSQ then J=LSQ JJQ(J)=JJQ(J)+1 'print I,J,JJQ(J) JSQ(JJQ(J),J)=I next I print for I=1 to LSQ : print using " ##:####";I;JJQ(I); : next print print using"Sorting atom coordinates (1 to ##) at ";LSQ; LY = csrlin : LX = pos() for KKk= 1 to LSQ kk=kkk locate LX,LY: print KK; if JJQ(KK)>1 then for II = 1 to JJQ(KK)-1 IPZMIN = 9999999 IJ = II for JJ = II+1 to JJQ(KK) J = JSQ(JJ,KK): ' print ii,jj,j if IP0(J) < IPZMIN then IJ=JJ: IPZMIN=IP0(J) next JJ if IP0(JSQ(II,KK))>IP0(JSQ(IJ,KK)) then swap JSQ(II,KK), JSQ(IJ,KK) 'print II;IJ end if next II end if 'print KK;JJQ(KK);: for I=1 to JJQ(KK): print JSQ(I,KK); : next : print next KKk ' N=0 for KK=1 to LSQ if JJQ(KK)>0 then for I=1 to JJQ(KK) N=N+1: ISEQ(N)=JSQ(I,KK) 'print KK;N;ISEQ(N) next I end if next KK 'print "SEQ";:for I= 1 to NTION : print ISEQ(I);: next : print 'stop return end ' SDRAW: '=============================================== Draw Single Projection inum=0 LXO = LPX/2 : IMENU=1 :iline=0 LYO = LPY/2+20 : CELL$="Y" SCALEX = LXO*2/AL(IXX)*0.96: SCALE = SCALEX SCALEY = (LYO-72)*2/AL(IXY)*0.96 if SCALEYXYZRANGE(2,1) then goto SD2 if XYZ(2,I)XYZRANGE(2,2) then goto SD2 if XYZ(3,I)XYZRANGE(2,3) then goto SD2 io = iele(i) ic = icc(io): 'print ii,i,io,ic '====== CHECK !!!!! IX = LXO + P(IXX,I) * SCALE IY = LYO - P(IXY,I) * SCALE R = RI1(IO) * (1 + P(IXZ,I) * .01) *SCALE/10 if r<=0 then r=0.1 circle (IX,IY), R, ic,,,,f circle (IX,IY), R, 1 'circle (IX,IY), R, 2 'print II;i;IX;IY;IO;IC;r '====================??????? if NET(IO) > .5 then for JJ = 1 to 6 j = ibond(jj,i) if J <> 0 then IF j < 0 THEN j = -j: jcol = 1 JX = LXO + P(IXX,J) * SCALE JY = LYO - P(IXY,J) * SCALE Z = P(IXZ, J) line (IX,IY)-(JX,JY), IC end if next JJ end if if RCUTA>0.1 and IDIATOM(I)>0 then '---------- DRAW DIATOMIC BONDS J = IDIATOM(I) JX = LXO + P(IXX,J) * SCALE JY = LYO - P(IXY,J) * SCALE Z = P(IXZ, J) for K=-1 to 1 for L= -1 to 1 line (IX+K,IY+L)-(JX+K,JY+L), IC next L: next K end if if inum=1 then locate iy/16,ix/8 if i<10 then print using "#";I; if I>=10 and ii<100 then print using "##";i; if I>=100 and ii<1000 then print using "###";i; if i>=1000 and ii<10000 then print using "####";i; if i>=10000 then print using "#####";i; end if SD2: next II if CELL$="Y" then for L1 = 1 to 7 if LBOX(IXZ,L1)=1 then I1=LBOX(1,L1): J1=LBOX(2,L1): K1=LBOX(3,L1) for L2 = L1+1 to 8 if LBOX(IXZ,L2)=1 then if abs(I1-LBOX(1,L2))+abs(J1-LBOX(2,L2))+abs(K1-LBOX(3,L2))=1 then IX = LXO+BOX(IXX,L1)*SCALE IY = LYO-BOX(IXY,L1)*SCALE JX = LXO+BOX(IXX,L2)*SCALE JY = LYO-BOX(IXY,L2)*SCALE line (IX,IY)-(JX,JY),1 if iline=1 then IX = LXO+BOX(IXX,L1)*SCALE IY = LYO-BOX(IXY,L1)*SCALE JX = LXO+BOX(IXX,L2)*SCALE JY = LYO-BOX(IXY,L2)*SCALE for a=-0.8 to 0.8 step 0.2 kX = LXO+BOX(IXX,L1)*SCALE*a kY = LYO-BOX(IXY,L1)*SCALE*a line (kx,IY)-(kx,JY),27 line (ix,ky)-(jx,ky),27 next a line (lxo,IY)-(lxo,JY),1 line (ix,lyo)-(jx,lyo),1 end if end if end if next L2 end if next L1 end if ' if iMENU=1 then LY=LCY: color 1 locate ly, 30 : print "q:Menu "; print "r:redraw "; print "+-:zoom "; print "abc:proj. "; print "xyz:range "; print "1-9:shift "; print "s:scale "; print "k:set_atoms "; print "lines"; end if SD9: A$ = inkey$: if A$ = "" then goto SD9 if A$ = "+" then SCALE = SCALE * 1.05: goto SD0 IF a$ = "-" THEN scale = scale * .95: GOTO sd0 IF a$ = "q" OR a$ = "Q" THEN RETURN if A$ = "r" or A$ = "R" then goto SD0 if A$ = "a" or A$ = "A" then IXZ = 1 IXX = 2 IXY = 3 if AL(IXX) < AL(IXY) then swap IXX,IXY gosub CALC_POSITIONS 'gosub SEQUENCE goto SD0 end if IF a$ = "b" OR a$ = "B" THEN IXZ = 2 IXX = 1 IXY = 3 if AL(IXX) < AL(IXY) then swap IXX,IXY gosub CALC_POSITIONS 'gosub SEQUENCE goto SD0 end if IF a$ = "c" OR a$ = "C" THEN IXZ = 3 IXX = 1 IXY = 2 if AL(IXX) < AL(IXY) then swap IXX,IXY gosub calc_positions 'gosub SEQUENCE: goto sD0 end if D=5 if A$="x" or A$="X" then input "range of x, 0 to 1 ? ";XYZRANGE(1,1),XYZRANGE(2,1) goto SD0 end if if A$="y" or A$="Y" then input "range of y, 0 to 1 ? ";XYZRANGE(1,2),XYZRANGE(2,2) goto SD0 end if if A$="z" or A$="Z" then input "range of z, 0 to 1 ? ";XYZRANGE(1,3),XYZRANGE(2,3) goto SD0 end if if A$="1" then LXO=LXO-D: LYO=LYO+D: goto SD0 if A$="2" then LYO=LYO+D: goto SD0 if A$="3" then LXO=LXO+D: LYO=LYO+D: goto SD0 if A$="4" then LXO=LXO-D: goto SD0 if A$="5" then goto SDRAW if A$="6" then LXO=LXO+D: goto SD0 if A$="7" then LXO=LXO-D: LYO=LYO-D: goto SD0 if A$="8" then LYO=LYO-D: goto SD0 if A$="9" then LXO=LXO+D: LYO=LYO-D: goto SD0 if A$="O" or A$="o" then if iMENU=1 then iMENU=0 else iMENU=1 end if end if if A$="s" or A$="S" then input "Scale";SCALE goto sD0 end if if A$="k" or A$="K" then gosub ATOMSETTING goto sd0 end if if a$="n" or a$="N" then if inum=0 then inum=1 else inum=0 end if goto sd0 end if if a$="l" or a$="L" then if iline=0 then iline=1 else iline=0 end if goto sd0 end if goto SD9 end ' ' ' STEREO: '========================================= Draw Stereoscopic Pair-view EYES = 8.8: HEYE = EYES / 2: ANGLEROT = 0 DIST = 60: NATOM = 0: AMAG = 1: IBOX = 0 SCALE = 10: XSFT=0 YSFT=0 ZSFT=0 ' XLIM = 14.99 ' width YLIM = 24.1 ' hight ZLIM = 9.45 ' depth LXO = LPX/2+30 LYO = (LPY+20)/2 LYO = (LPY+20)/2 XLO = LXO - 175: XRO = LXO + 175 ' STR0: print "Calculate stereo coordinates at "; : LY = csrlin LX = pos() for I = 1 to NTION : '----------- Calculate stereo coordinates of atoms X = (P(IXX,I)+XSFT)*AMAG Y = (P(IXY,I)+YSFT)*AMAG Z = (P(IXZ,I)+ZSFT)*AMAG '----------------------------------- Rotation around y-axis ANGLR = ANGLEROT * 3.14159 / 180 CANGL = cos(ANGLR) SANGL = sin(ANGLR) XS = X * CANGL - Z * SANGL Z = Z * CANGL + X * SANGL X = XS '----------------------------------------------------------- if abs(X) > XLIM then YS(I) = -9999: goto STR1: ' x horizontal if abs(Y) > YLIM then YS(I) = -9999: goto STR1: ' y vertical if abs(Z) > ZLIM then YS(I) = -9999: goto STR1: ' z depth NATOM = NATOM + 1 DD = DIST / (DIST - Z) if sVIEW$="PARA" then XL(I) = (X + HEYE) * DD - HEYE XR(I) = (X - HEYE) * DD + HEYE end if if sVIEW$="CROS" then XL(I) = (X - HEYE) * DD + HEYE XR(I) = (X + HEYE) * DD - HEYE end if YS(I) = Y * DIST / (DIST - Z) ZS(I) = Z locate ly,LX: print I; STR1: next I print print "Rearrange sequence of atoms at "; : LY = csrlin LX = pos() for I=1 to NTION: JSEQ(I)=I: next for II = 1 to NTION-1 PZMIN = 99999.9: IJ = II for JJ = II+1 to NTION J = JSEQ(JJ) if ZS(J) < PZMIN then IJ=JJ: PZMIN=ZS(J) next JJ if ZS(JSEQ(II))>ZS(JSEQ(IJ)) then swap JSEQ(II), JSEQ(IJ) end if locate ly,LX : print II; next II ' '-------------------------- Calculate stereo coord. of lattice points STR2: NL = 0 KX = NBOX(1) / 2: AKX = KX: 'IF nbox(1) MOD 2 = 0 THEN akx = kx - .5 KY = NBOX(2) / 2: AKY = KY: 'IF nbox(2) MOD 2 = 0 THEN aky = ky - .5 KZ = NBOX(3) / 2: AKZ = KZ: 'IF nbox(3) MOD 2 = 0 THEN akz = kz - .5 for IX = 0 to 1: for IY = 0 to 1: for IZ = 0 to 1 NL = NL + 1 AX = AKX + IX: XX = AX / NBOX(1) - 0.5 AY = AKY + IY: YY = AY / NBOX(2) - 0.5 AZ = AKZ + IZ: ZZ = AZ / NBOX(3) - 0.5 PP(1) = (H(1,1) * XX + H(1,2) * YY + H(1,3) * ZZ) * AMAG PP(2) = (H(2,1) * XX + H(2,2) * YY + H(2,3) * ZZ) * AMAG PP(3) = (H(3,1) * XX + H(3,2) * YY + H(3,3) * ZZ) * AMAG X = PP(IXX): KBOX(1,NL) = IX Y = PP(IXY): KBOX(2,NL) = IY Z = PP(IXZ): KBOX(3,NL) = IZ '----------------------------------- Rotation around y-axis ANGLR = ANGLEROT * 3.14159 / 180 CANGL = cos(ANGLR) SANGL = sin(ANGLR) XS = X * CANGL - Z * SANGL Z = Z * CANGL + X * SANGL X = XS '----------------------------------------------------------- DD = DIST / (DIST - Z) if sVIEW$="PARA" then XLL(NL) = (X + HEYE) * DD - HEYE XLR(NL) = (X - HEYE) * DD + HEYE end if if sVIEW$="CROS" then XLL(NL) = (X - HEYE) * DD + HEYE XLR(NL) = (X + HEYE) * DD - HEYE end if YLS(NL) = Y * DD next IZ: next IY: next IX ' gosub SCREENTITLE NA = 0: color 6 ' for II = 1 to NTION : ' --------------------------------------- Plot atoms I = JSEQ(II) IO = IELE(I): if YS(I) < -9998 then goto STR3 NA = NA + 1 'locate 0,16: print using "####/####"; NA; NATOM; 'locate 0,17: print using " /####"; NTION; IC = ICC(IO) Z = ZS(I): X1 = XL(I) X2 = XR(I): Y = YS(I) r = ri1(io) * (1 + z * .05) * amag iy = lyo - y * scale IXL = X1 * SCALE + XLO: circle (IXL,IY), R, ic,,,,f circle (ixl,iy), R, 1 IXR = X2 * SCALE + XRO: circle (IXR,IY), R, ic,,,,f circle (ixr,iy), R, 1 if NET(IO) > .5 then for JJ = 1 to 5 if IBOND(JJ,I) <> 0 then J = IBOND(JJ,I) JCOL = IC if J < 0 then J = -J: JCOL = 1 if YS(J) > -9998 then JY = LYO - YS(J) * SCALE JX = XL(J)*SCALE+XLO: line (IXL,IY)-(JX,JY),JCOL JX = XR(J)*SCALE+XRO: line (IXR,IY)-(JX,JY),JCOL end if end if next JJ end if if RCUTA>0.1 and IDIATOM(I)>0 then J = IDIATOM(I) if YS(J)> -9998 then JY = LYO - YS(J) * SCALE JX = XL(J)*SCALE+XLO: line (IXL,IY)-(JX,JY),JCOL JX = XR(J)*SCALE+XRO: line (IXR,IY)-(JX,JY),JCOL end if end if STR3: next II ' if IBOX = 1 then '----------------------------------------------- Display unit cell for I = 1 to 7: for J = I to 8 NN = abs(KBOX(1,I)-KBOX(1,J))+abs(KBOX(2,I)-KBOX(2,J))+abs(KBOX(3,I) - KBOX(3,J)) if NN > .9 and NN < 1.1 then IY = LYO - YLS(I) * SCALE JY = LYO - YLS(J) * SCALE IX = XLL(I) * SCALE + XLO JX = XLL(J) * SCALE + XLO: line (IX,IY)-(JX,JY),1 IX = XLR(I) * SCALE + XRO JX = XLR(J) * SCALE + XRO: line (IX,IY)-(JX,JY),1 end if next J: next I end if ' color 0 locate 16,1: print "[q] menu"; locate 17,1: print "[l] cell"; locate 18,1: print "[w] area" locate 19,1: print "[r] angle"; locate 20,1: print "[+-] size"; locate 21,1: print "[d] redraw"; locate 22,1: print "[abc] proj."; STR9: A$ = inkey$: if A$ = "" then goto STR9 IF a$ = "l" OR a$ = "L" THEN IF ibox = 0 THEN ibox = 1 ELSE ibox = 0 END IF GOTO str3 end if if A$="W" or A$="w" then print "Crrent area: Hori=";XLIM;" Ver=";YLIM;" Dep=";ZLIM input "H,V,D";XLIM,YLIM,ZLIM goto str0 end if IF a$ = "r" OR a$ = "R" THEN color 1 input "Angle "; ANGLEROT color 7 goto STR0 END IF IF a$ = "+" THEN scale = scale * 1.05: GOTO str0 IF a$ = "-" THEN scale = scale * .95: GOTO str0 IF a$ = "d" OR a$ = "D" THEN GOTO stereo if A$ = "a" or A$ = "A" then ixz = 1 ixx = 2: ixy = 3 IF al(ixx) > al(ixy) THEN ixx = 3 IXY = 6 - IXX - IXZ gosub CALC_POSITIONS GOSUB sequence: GOTO stereo end if if A$ = "b" or A$ = "B" then ixz = 2 ixx = 1: ixy = 3 IF al(ixx) > al(ixx) THEN ixx = 3 IXY = 6 - IXX - IXZ gosub CALC_POSITIONS GOSUB sequence: GOTO stereo end if if A$ = "c" or A$ = "C" then ixz = 3 ixx = 1: ixy = 2 IF al(ixx) > al(ixx) THEN ixx = 2 IXY = 6 - IXX - IXZ gosub calc_positions GOSUB sequence: GOTO stereo end if if A$="s" or A$="S" then input "Xshift,Yshift,Zshift";XSFT,YSFT,ZSFT goto STR0 end if if A$ = "q" or A$ = "Q" then goto STR99 if A$ = "e" or A$ = "E" then goto STR99 d=5 if A$="1" then LXO=LXO-D: LYO=LYO+D: goto STR0 if A$="2" then LYO=LYO+D: goto STR0 if A$="3" then LXO=LXO+D: LYO=LYO+D: goto STR0 if A$="4" then LXO=LXO-D: goto STR0 if A$="5" then goto STEREO if A$="6" then LXO=LXO+D: goto STR0 if A$="7" then LXO=LXO-D: LYO=LYO-D: goto STR0 if A$="8" then LYO=LYO-D: goto STR0 if A$="9" then LXO=LXO+D: LYO=LYO-D: goto STR0 goto STR9 STR99: RETURN end ' ' SCREENTITLE: '============================================= Plot Screen title cls 0: line (0,0)-(LPX,LPY),15,bf,15 color 0 locate 1,2: print TITLE$ locate 2,2: print using "T=#####K";TEMP locate 3,2 if SPRES(1)>0.99 then if SPRES(1) < 10 then print using "P=##.###GPa"; SPRES(1) if SPRES(1) >= 10 then print using "P=###.##GPa"; SPRES(1) else print using "P=###.##MPa";SPRES(1)*1000 end if locate 2,16: print "D/g/cm3":locate 3,15: print using "##.####"; DNSTY locate 2,27: print "Cell/A": locate 2,31: print using " ###.###"; AL(1) locate 3,27: print " ": locate 2,39: print using " ###.###"; AL(2) locate 2,47: print using " ###.###"; AL(3) locate 3,31: print using " ###.###"; AL(4) locate 3,39: print using " ###.###"; AL(5) locate 3,47: print using " ###.###"; AL(6) locate 2,70: print TM$ locate 3,70: print mid$(DT$,1,8) ' locate 1,72: color 2: if sVIEW$="PARA" then print "[Parallel]"; if sVIEW$="CROS" then print " [Cross]"; ' XX0=20: YY0=LPY-10: color 1 : '-------------------------- Index of axes line (XX0, YY0)-(XX0+30, YY0), 1: locate LCY+0,8 : print IX$(IXX); line (XX0, YY0)-(XX0, YY0-35), 1: locate LCY-3,2 : print IX$(IXY); locate LCY+0,2 : print IX$(IXZ); return end ' ' FILE05: '===================================================== Read File05.dat open DR$ + "FILE05.DAT" for input as #1 line input #1, A0$ : print a0$ for I = 2 to 8: line input #1, A$: next for I = 1 to LEL line input #1, A$: print " "; A$ IO = val(left$(A$,1)): if IO = 0 then goto F51 ATOM$(IO) = mid$(A$,3,3) ZI(IO) = val(mid$(A$,11,10)) ' PRINT io, atom$(io), Zi(io) next I line input #1, A$ ' sleep F51: for II = 1 to LEL: NET(II) = 0: next II MET = 0 IDIA1 = 0 IDIA2 = 0 ' sleep F5A: if eof(1)<0 then goto F52 line input #1, A$: ' PRINT a$ if left$(A$,4) = "STOP" then goto F52 if left$(A$,4) = "END " then goto F52 end if if left$(A$, 7) = "NETWORK" then NET1 = val(mid$(A$,11,10)) if NET1 > 0 then MET = 1: NET(NET1) = 1 Net2 = val(mid$(A$,21,10)) if NET2 > 0 then MET = 2: NET(NET2) = 1 end if print "Network former(s) is (are) "; ATOM$(NET1); if NET2 > 0 then print ATOM$(NET2) end if if left$(A$,8) = "DIATOMIC" then RCUTA = val(mid$(A$,11,20)) IDIA1 = val(mid$(A$,21,30)) IDIA2 = val(mid$(A$,31,40)) end if if left$(A$,9) = "POLYATOMS" then RCUTA = val(mid$(A$,11,20)) : if rcuta<0 then rcuta=-rcuta IDIA1 = val(mid$(A$,21,30)) IDIA2 = val(mid$(A$,31,40)) end if goto F5A F52: close #1: print NBOX(1) = 1: NBOX(2) = 1: NBOX(3) = 1 if left$(A0$,2) = "XD" then open DR$+"FILE10.DAT" for input as #1 line input #1, A$ input #1, NBOX(1), NBOX(2), NBOX(3) print NBOX(1), NBOX(2), NBOX(3) close #1 end if return end ' ' FILE07: '=================================================== Read File07.dat IFORM7=0 open DR$ + "FILE07.DAT" for input as #1 line input #1, A0$: TITLE$ = left$(A0$, 60) NJOB1 = val(mid$(A0$,61,5)) NJOB2 = val(mid$(A0$,66,5)) line input #1, A1$: 'PRINT a1$ NTION = val(left$(A1$,7)) NCOMPO = val(mid$(A1$,8,3)) NSTEP = val(mid$(A1$,11,10)) NHIST = val(mid$(A1$,51,5)) line input #1, A2$: 'PRINT a2$ line input #1, A3$: 'PRINT a3$ line input #1, A4$: 'PRINT a4$ line input #1, A5$: 'PRINT a5$ for I = 1 to NCOMPO: I6 = 6*(I-1)+1 ATOM$(I) = mid$(A2$,I6+2,4) NION(I) = val(mid$(A3$,I6,6)) IONS(1,I) = val(mid$(A4$,I6,6)) IONS(2,I) = val(mid$(A5$,I6,6)) next I ' line input #1, A6$: 'PRINT a6$ TEMP = val(mid$(A6$, 1,10)) DELTA = val(mid$(A6$,11,10)) TGET = val(mid$(A6$,21,10)) SPRES(1) = val(mid$(A6$,31,10)) SPRES(2) = val(mid$(A6$,41,10)) SPRES(3) = val(mid$(A6$,51,10)) line input #1, A7$: print A7$ DT = val(mid$(A7$, 1,10)) AL(1) = val(mid$(A7$,21,10)) AL(2) = val(mid$(A7$,31,10)) AL(3) = val(mid$(A7$,41,10)) V$ = mid$(A7$,51,10): AL(4)=val(V$) V$ = mid$(A7$,61,10): AL(5)=val(V$) V$ = mid$(A7$,71,10): AL(6)=val(V$) if AL(4)<2.0 then for I=4 to 6 ALC(I)=AL(I) S=sqr(1-AL(I)^2) if AL(I)=0 then AL(I) = 90.0 else AL(I)=atn(S/AL(I))*180/3.1415926 if AL(I)<0 then AL(I)=AL(I)+180 end if next I end if line input #1, A8$: 'PRINT a8$ DNSTY = val(mid$(A8$,1,10)) HTENSOR$ = mid$(A8$,11,10) if HTENSOR$ = "H-TENSOR " then line input #1, A81$ line input #1, A82$ line input #1, A83$ for J = 1 to 3 JJ = 11 + (J-1) * 20 V$=mid$(A81$,JJ,20): H(1,J)=val(V$) V$=mid$(A82$,JJ,20): H(2,J)=val(V$) V$=mid$(A83$,JJ,20): H(3,J)=val(V$) next J else AL4S = sqr(1-ALC(4)^2) H(1,3) = 0 H(2,3) = 0 H(3,3) = AL(3) H(1,2) = 0 H(2,2) = AL(2) * AL4S H(3,2) = AL(2) * ALC(4) H(3,1) = AL(1) * ALC(5) H(2,1) = AL(1) * ALC(6) * AL4S H(1,1) = AL(1) * sqr(1-ALC(5)^2 - (ALC(6)*AL4S)^2) end if ' print TITLE$ print "NTION="; NTION, "NSTEP="; NSTEP, print "JOB NO.=("; NJOB1; ","; NJOB2; ") NHIST="; NHIST for I = 1 to 5: I10 = 10 * (I-1) if NION(I) > 0 then print using " ## "; I; : print ATOM$(I); print using " #####"; NION(I); print using " ##### #####"; IONS(1, I); IONS(2, I) end if next I print "DT="; DT; "x10-15sec "; print "Box(1 - 3) = "; AL(1); AL(2); AL(3); "" print "Temp.=";TEMP; " DeltaT=";DELTA;"K/10step Tget =";TGET;"K" print "Density=";DNSTY;"g/cm3 SPRES(1-3)=";SPRES(1);SPRES(2);SPRES(3) print ' print "Reading atomic coordinates at "; : LY = csrlin: LX = pos() for I = 1 to NTION: line input #1, A$ : 'print a$ if len(a$)>86 and len(a$)<105 then iform7=1 for J = 1 to 3 if IFORM7=0 then V = val(mid$(A$, J*9-8,9)) if IFORM7=1 then V = val(mid$(A$, J*10-9,10)) P0(J,I) = V XYZ(J,I) = V next J 'print i;p0(1,i);p0(2,i);p0(3,i) IP0(I) = ((P0(1,I)*AL(1))^2 + (P0(2,I)*AL(2))^2 + (P0(3,I)*AL(3))^2)*2 ' FOR j = 1 TO 3: v(j,i) = (VAL(MID$(a$,21+j*8))-5) * .1: NEXT j ' FOR J = 1 TO 3: p0(J,I) = VAL(MID$(a$,45+j*9)): NEXT j if I mod 10 = 0 then locate LX,LY : print using "##### "; I; FOR J=1 TO 3 : PRINT USING " #.####";xyz(J,I);: NEXT J : PRINT " "; ' FOR J=1 TO 3 : PRINT USING " #.####";v(J,I); : NEXT J : PRINT " "; ' FOR J=1 TO 3 : PRINT USING " #.####";p0(J,I); : NEXT J : PRINT end if next I ' MHIST = int(NHIST / 4): print if NHIST>0 then for I = 1 to MHIST: line input #1, A$: print A$ ' INPUT #1,(7800) ((JHIST(J,I),J=1,4),I=1,NHIST) next I end if close #1: print : ' sleep return end ' ' ATOMS: for I = 1 to LEL DLIM2(I) =0: RI1(I)=0: ICC(I)=0 DLIM2I(I)=0: RI2(I) =0 ATM$ = left$(ATOM$(I),2) ' if ATM$="H " then RI1(I)=2: ICC(I)= 3: DLIM2(I) = 1.3^2 end if if ATM$="D " then RI1(I)=2: ICC(I)= 5: DLIM2(I) = 1.3^2 end if if ATM$="T " then RI1(I)=2: ICC(I)= 0: DLIM2(I) = 1.3^2 end if if ATM$="HE" or ATM$="He" then RI1(I)= 5: ICC(I)= 5 ' if ATM$="LI" or ATM$="Li" then RI1(I)=4: ICC(I)=2 if ATM$="BE" or ATM$="Be" then RI1(I)=3: ICC(I)=4 if ATM$="B " then RI1(I)=3: ICC(I)=4: DLIM2(I)=1.7^2 if ATM$="C " then RI1(I)=3.5: ICC(I)=11: DLIM2(I)=1.8^2 if ZI(I) < 0 then RI1(I)=7: RI2(I)=5: ICC(I)=5 end if if ATM$="N " then RI1(I)=4: ICC(I)=5: DLIM2(I)=1.4^2 if ZI(I) < 0 then RI1(I)=7: RI2(I)=5: ICC(I)=5 end if if ATM$="O " then RI1(I)=7: ICC(I)=44 if ATM$="F " then RI1(I)=7: ICC(I)=43 if ATM$="NE" or ATM$="Ne" then RI1(I)=7: ICC(I)= 88 ' if ATM$="NA" or ATM$="Na" then RI1(I)=5: ICC(I)=42 if ATM$="MG" or ATM$="Mg" then RI1(I)=5.5: ICC(I)= 2 if ATM$="AL" or ATM$="Al" then RI1(I)=4: ICC(I)= 5: DLIM2(I)=2.1^2 if ATM$="SI" or ATM$="Si" then RI1(I)=3.5: ICC(I)= 3: DLIM2(I)=1.9^2 if ATM$="P " then RI1(I)=3.5: ICC(I)= 3: DLIM2(I)=1.8^2 if ATM$="S " then RI1(I)=9: ICC(I)= 44: DLIM2(I)=2.2^2 if ATM$="CL" or ATM$="Cl" then RI1(I)=9: ICC(I)= 92 if ATM$="AR" or ATM$="Ar" then RI1(I)=9: ICC(I)= 85 ' if ATM$="K " then RI1(I)=7: ICC(I)= 40 if ATM$="CA" or ATM$="Ca" then RI1(I)=6.5: ICC(I)= 37 if ATM$="SC" or ATM$="Sc" then RI1(I)=6: ICC(I)= 58 if ATM$="TI" or ATM$="Ti" then RI1(I)=5: ICC(I)= 3: DLIM2(I)=2.2^2 if ATM$="MN" or ATM$="Mn" then RI1(I)=5.5: ICC(I)= 34 if ATM$="FE" or ATM$="Fe" then RI1(I)=5.5: ICC(I)= 34 if ATM$="CU" or ATM$="Cu" then RI1(I)=5.0: ICC(I)= 54 if ATM$="GA" or ATM$="Ga" then RI1(I)=4: ICC(I)= 34 if ATM$="GE" or ATM$="Ge" then RI1(I)=4: ICC(I)= 3 ' if ATM$="RB" or ATM$="Rb" then RI1(I)=7: ICC(I)= 38 if ATM$="SR" or ATM$="Sr" then RI1(I)=7: ICC(I)= 36 if ATM$="Y " then RI1(I)=6: ICC(I)= 9 if ATM$="ZR" or ATM$="Zr" then RI1(I)=6: ICC(I)= 49 if ATM$="AG" or ATM$="Ag" then RI1(I)=6: ICC(I)= 28 if ATM$="I " then RI1(I)=8: icc(i)= 59 ' if ATM$="CS" or ATM$="Cs" then RI1(I)=8: ICC(I)= 35 if ATM$="BA" or ATM$="Ba" then RI1(I)=8: ICC(I)= 161 if ATM$="LA" or ATM$="La" then RI1(I)=5: ICC(I)= 4 if ATM$="CE" or ATM$="Ce" then RI1(I)=4: ICC(I)= 53 if ATM$="GD" or ATM$="Gd" then RI1(I)=5: ICC(I)= 50 if ATM$="PB" or ATM$="Pb" then RI1(I)=6: ICC(I)= 24 if ATM$="VA" or ATM$="Va" then RI1(I)=7: ICC(I)= 99 circle (100+I*20,100),RI1(I),ICC(I) ' print I;ATM$;RI1(I) next I return end '===============================================================================