' ********************************************************* ' *** Free BASIC *** ' *** Program : MD-MSD-PLOT (Mean square displacement *** ' *** of atoms) *** ' *** 1991-12-14 by Kats *** ' *** 1992-05-30 (for MXDNEW) by Kats *** ' *** 1993-11-27 by Kats *** ' *** 1994-02-04 (for LIPS3) by kats *** ' *** 2000-04-29 500,000 steps by Kats *** ' *** 2005-08-15 FreeBASIC version *** ' *** 2011-09-18 60000 atoms *** ' ********************************************************* ' deflng I-N, X-X : defdbl A-h,o-w,y-z ' Screen 19 : xscreen= 800 : yscreen= 600 : lcx=100 : lcy=37 '(VGA) XCL=100 MAXLINE = 37 ' CLS 0 LCOMPO = 8: LPOSI = 2005 DIM atom$(lcompo), ions(2, lcompo), dmax(1, lcompo) DIM fy(lposi), nion(lcompo), dms(lposi, lcompo) DIM tm(lposi), iypos(Maxline), spres(6),al(6), diff(lcompo) ' dim dr0$ dim np,ncompo,njob1,njob2,jx,jy,icol,iy0 dim dt ' hajime: dr$ = "" '-------------------------------------------- Data directory PRINT "Type a directory name (...\), or hit [return] for the current directory." INPUT " (q for end)"; dr$ IF dr$ = "Q" OR dr0$ = "q" THEN END IF LEN(dr$) > 0 AND RIGHT$(dr$, 1) <> ":" THEN i = LEN(dr$) IF RIGHT$(dr$, 1) <> "\" THEN dr$ = dr$ + "\" END IF ' ' L1: GOSUB DREADf07 '------------------------------------- Read file07.dat gosub DREADF06MSD '------------------------------------- READ F06MSD.DAT ' PRINT : PRINT "Max(time span):"; tm(np); " picoseconds" PRINT : PRINT "Min(m.s.d.) "; FOR i = 1 TO ncompo: PRINT USING " #####.###"; dmax(0, i); : NEXT PRINT : PRINT "Max(m.s.d.) "; FOR i = 1 TO ncompo: PRINT USING " #####.###"; dmax(1, i); : NEXT ' PRINT : PRINT INPUT "Time span (default = 10 picoseconds)"; ntime IF ntime < .1 THEN ntime = 10 INPUT "Y(m.s.d.) span (default = 2 A^2)"; Yscl IF Yscl < .1 THEN Yscl = 2 ' OSCREEN = 1 P: cls 0 '-------------------------------------- Plot and Calculate approx. D NEND = NP * 3 / 4 NSTART = 2: if NEND >= 10 then NSTART = 4 if NEND >= 20 then NSTART = 8 if NEND >= 50 then NSTART = 22 if NEND >= 100 then NSTART = 40 if NEND >= 200 then NSTART = 60 if NEND >= 500 then NSTART =100 if NEND >=1000 then NSTART =200 if nend >=2000 then nstart =400 nstart=nend*1/4 ' P1: GOSUB sflame for I = 1 to MAXLINE: IYPOS(I) = 0: next I FOR ion = 1 TO ncompo : ' print "Ncompo=";ncompo;" ion=";ion if NION(ION) > 0 then ICOL = ION gosub SPLOT : ' print "splot" gosub DCALC : ' print "dcalc" diff(ion)=ddd end if NEXT ion OPEN "DIFF.TXT" for output as #1 print #1, "Temp/K"; for io=1 to ncompo : print #1," ";atom$(io);" "; : next io print #1," " print #1, using "####.# ";temp; for io=1 to ncompo : print #1, using " ######.#######E-9";diff(io)*1.0E9; : next io print #1, close #1 ' MENU: color 4,15 if OSCREEN = 1 then locate 3,52: print "D-calc" locate 4,52: print "Quit" end if LOCATE 5,52: print "?" aaa: a$= inkey$: if a$="" or a$= " " then goto aaa IF a$ = "d" or a$ = "D" THEN locate 18,22: input "Tmin and Tmax (ps)";Tmin,Tmax amin = 99999: amax = 99999 for i = 1 to np if abs(Tmin-tm(i)) < amin then nstart = i amin = abs(Tmin - tm(i)) end if if abs(Tmax-tm(i)) < amax then nend = i amax = abs(Tmax - tm(i)) end if next i goto P1 end if ' if A$ = "q" or A$ = "Q" then end if A$ = "o" or A$ = "O" then if OSCREEN=1 then OSCREEN = 0 else OSCREEN = 1 end if goto P end if goto MENU ' ' DREADf07: '-------------------------------------------------------- dread-f07 OPEN dr$ + "FILE07.DAT" FOR INPUT AS #2: ' READ FROM FILE07.DAT LINE INPUT #2, a0$: Title$ = LEFT$(a0$, 60) NHIST = VAL(MID$(a0$, 66, 5)) LINE INPUT #2, a1$: PRINT a1$ ntion = VAL(LEFT$(a1$, 7)) ncompo = VAL(MID$(a1$, 8, 3)) anstep = VAL(MID$(a1$, 11, 10)) npstep = VAL(MID$(a1$, 41, 5)) LINE INPUT #2, a2$: PRINT a2$ LINE INPUT #2, a3$: PRINT a3$ LINE INPUT #2, a4$: PRINT a4$ LINE INPUT #2, a5$: PRINT a5$ FOR i = 1 TO ncompo: i6 = 6 * (i - 1) + 1 atom$(i) = MID$(a2$, i6 + 2, 6) 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 #2, 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)): pres = spres(1) spres(2) = VAL(MID$(a6$, 41, 10)) spres(3) = VAL(MID$(a6$, 51, 10)) LINE INPUT #2, a7$: PRINT a7$ DTIME = 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)) LINE INPUT #2, a8$: PRINT a8$ density = VAL(MID$(a8$, 1, 10)) Htensor$ = MID$(a8$, 11, 10) IF Htensor$ = "H-TENSOR " THEN LINE INPUT #2, a81$ LINE INPUT #2, a82$ LINE INPUT #2, a83$ END IF PRINT Title$ PRINT "NTION="; ntion, "NSTEP="; anstep, PRINT "JOB NO.="; njob1; njob2, "NHIST="; NHIST FOR i = 1 TO ncompo: i10 = 10 * (i - 1) PRINT " "; i; atom$(i), nion(i), ions(1, i); ions(2, i) NEXT i PRINT "T="; temp, "DT="; dt, "AL(1-3)="; al(1); al(2); al(3) PRINT "DELTA="; delta, "TGET ="; tget PRINT "Density="; density, "SPRES(1-3)="; spres(1); spres(2); spres(3) PRINT ' ' PRINT "READING ATOMIC COORDINATES AT "; : ly = CSRLIN: lx = POS(0) ' FOR i = 1 TO ntion: LINE INPUT #1, a$ ' FOR j = 1 TO 3: p(j, i) = VAL(MID$(a$, j * 9 - 8)): NEXT j ' FOR j = 1 TO 3: v(j, i) = VAL(MID$(a$, 21 + j * 8)) * .0625: NEXT j ' FOR j = 1 TO 3: p0(j, i) = VAL(MID$(a$, 45 + j * 9)): NEXT j ' LOCATE ly,lx: PRINT USING "#### "; i; ' NEXT i ' ' mhist = INT(NHIST / 4): PRINT ' 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 CLOSE #2: PRINT RETURN ' ' DREADf06msd: '----------------------------------------------- data read-f06msd print print "============= Read from f06msd.dat ==============" OPEN dr$ + "F06MSD.DAT" FOR INPUT AS #1 FOR i = 1 TO 4: LINE INPUT #1, l$: PRINT l$: NEXT FOR i = 5 TO 6: LINE INPUT #1, l$: PRINT LEFT$(l$, 75): NEXT FOR i = 7 TO 10: LINE INPUT #1, l$: NEXT FOR i = 1 TO ncompo: dmax(0, i) = 99999 dmax(1, i) = 0: NEXT np = 0 D1: np = np + 1 LINE INPUT #1, l$: tm(np) = VAL(LEFT$(l$, 10)) FOR i = 1 TO ncompo dms(np, i) = VAL(MID$(l$, 9*i+2, 9)) IF dmax(0, i) > dms(np, i) THEN dmax(0, i) = dms(np, i) IF dmax(1, i) < dms(np, i) THEN dmax(1, i) = dms(np, i) NEXT IF np MOD 10 = 0 THEN PRINT USING "###### #####.##"; np; tm(np); for k = 1 to ncompo PRINT USING " ####.####"; dms(np, k); next k print END IF IF tm(np) > .001 THEN GOTO D1 ' np = np - 1 CLOSE #1: ' STOP RETURN ' ' sflame: ' -------------------------------------------- Screen title and flame line (0,0)-(xscreen,yscreen),15,bf XO = 43: XL = 700 YO = maxline*16-36: YL = YO-20 '------------------------------------------------------ horizontal axis (time) color 0,15 MST = 1 IF ntime > 10 THEN mst = 2 IF ntime > 20 THEN mst = 5 if NTIME > 50 then MST = 10 if ntime > 140 then mst = 20 if ntime > 299 then mst = 50 if ntime > 599 then mst = 100 if ntime >1199 then mst = 200 FOR i = 0 TO ntime ix = xl / ntime * i + xo IF ix / 8 <= lcx THEN IF i MOD mst = 0 THEN LINE (ix, yo-yl)-(ix, yo), 7 line (ix+1,yo-yl)-(ix+1,yo), 7 if i<1000 then locate maxline-1,IX/8-1 PRINT USING "###"; i; else locate maxline-1,IX/8-2 PRINT USING "####"; i; end if END IF END IF NEXT i '---------------------------------------------------- Vertical axis (m.s.d.) ystep = 0.2: if yscl > 1.1 then ystep = 0.5 if yscl > 2.9 then ystep = 1 IF Yscl > 9.9 THEN ystep = 2 IF Yscl > 19.9 THEN ystep = 5 IF Yscl > 49.9 THEN ystep = 10 IF Yscl > 99.9 THEN ystep = 20 if Yscl > 249.9 then ystep = 50 if Yscl > 999.9 then ystep = 100 if Yscl >2999.9 then ystep = 200 if Yscl >5999.9 then ystep = 500 if Yscl >9999.9 then ystep = 1000 FOR y = 0 TO Yscl STEP ystep '---------------- Vertical axis (m.s.d.) iy = yo - yl / Yscl * y LINE (xo, iy )-(xo + xl, iy ), 1, B LINE (xo, iy-1)-(xo + xl, iy-1), 1, B LIY = (IY+1)/16-0.4: if LIY <= 0 then LIY = 0 if liy>= Maxline then LIY=maxline locate LIY+1,1 if Yscl<1000 then print using "###.#"; y; else print using "#####"; y; end if NEXT y line (XO, YO-YL )-(XO+XL, YO ), 0, b line (XO+1, YO-YL-1)-(XO+XL+1,YO-1), 0, b ' line (68,88)-(247,215),15,bf ' -------- blank for D-coeff locate MAXLINE,40 : print "Time / ps"; locate 10,62 : print "M.s.d. / A2" ' ' line (XO+15,1)-(280,77),15,bf ' ------- blank for title etc. color 0,15 locate 1,6: print TITLE$: print LOCATE 3,8: PRINT USING " T=#### K P=##.#### GPa"; temp; pres; LOCATE 4,8: PRINT USING " D=#.#### g/cm3"; density print RETURN ' SPLOT: '--------------------------------------------------------- screen plot FACT = yl / Yscl: MULT = 1: FACT0 = FACT DM = DMAX(1,ION)*NTIME/TM(NP) if DM > YSCL * 1 then FACT = FACT0 / 5: MULT = 5 if DM > YSCL * 5 then FACT = FACT0 / 10: MULT = 10 if DM > YSCL * 10 then FACT = FACT0 / 20: MULT = 20 if DM > YSCL * 20 then FACT = FACT0 / 25: MULT = 25 if DM > YSCL * 25 then FACT = FACT0 / 50: MULT = 50 if DM > YSCL * 50 then FACT = FACT0 / 100: MULT = 100 if DM > YSCL * 100 then FACT = FACT0 / 200: MULT = 200 if DM > YSCL * 200 then FACT = FACT0 / 250: MULT = 250 if DM > YSCL * 250 then FACT = FACT0 / 500: MULT = 500 if DM > YSCL * 500 then FACT = FACT0 /1000: MULT =1000 : print "NP=";np for J= 0 to 1 for I = 1 to NP ' --------------------------------------- m.s.d. plot IX = XO + TM(I) * XL / NTIME: if IX>XO+XL then goto SP0 iy = yo - dms(i,ion) * FACT - j PSET (ix, iy), 8-ion if I > 1 then line (JX, JY)-(IX, IY),ICOL jx = ix: jy = iy NEXT i SP0: next j color ICOL,15: IY = IY/16: if IY <= 0 then IY=1 IF ion = 1 THEN GOTO SP1 IY0 = iy: IF iy > maxline-3 THEN iy = maxline-3 IF iypos(iy) THEN iy = iy + 1: IF iy > maxline-3 THEN iy = maxline-3 IF iypos(iy) THEN iy = IY0 - 1: IF iy = 0 THEN iy = IY0 + 2 SP1: if iy