Get MSD (mean square displacement) from MD trajectory =========== The first script convert XDATCAR to XYZ file The second script extract MSD info from XYZ file :: import numpy as np franum=0 count=0 latt=[] axyz=[] ele=[] etnum=[] #fixflag=input(" Seletive Dynamics is T or F: ") fixflag=False tlnum=0 pf=open("XDATCAR","r") fileend=0 tmpstr=pf.readline() count=count+1 tmpstr=pf.readline() count=count+1 #### Begin Read Lattice Vector ##### latt.append([]) for j in range(0,3): tl=pf.readline() count=count+1 tlst=tl.split() tlf=[float(x) for x in tlst[0:3]] latt[franum].append(tlf) print latt #### Begin Read Element Num ##### tmpstr=pf.readline() count=count+1 [ele.append(x) for x in tmpstr.split()] print ele tnumstr=pf.readline() print "read ele num" print tnumstr count=count+1 etnum=[int(x) for x in tnumstr.split()] print etnum totnum=sum(etnum) if fixflag: tmpstr=pf.readline() count=count+1 #### Begin Read Atom Fraction Coordination ##### while not fileend: tmpstr=pf.readline() if tmpstr == '': print "file read done" fileend=1 break count=count+1 axyz.append([]) for j in range(0,totnum): txyz=pf.readline() if txyz == '': print "Error:Missing Atom In Frame %d"%franum fileend=1 break count=count+1 taxyz=txyz.split() xyz=[float(x) for x in taxyz[0:3]] axyz[franum].append(xyz) franum=franum+1 LATTICE=np.array(latt) ATXYZ=np.array(axyz) print ATXYZ[0] CXYZ=[] outfile="LiXY" of=open(outfile,"w") for i in range(0,franum,10): # ATXYZ=np.array(axyz[i]) # LATT=LATTICE[i] # tCXYZ=np.dot(ATXYZ,LATT) # CXYZ.append(tCXYZ) # print >> of," %d "%(enum[i][2]) # print >> of," %s %d %s "%("### ",i+1," ####") for j in range(0,etnum[0]): print >> of," %f %f "%(ATXYZ[i][j][0],ATXYZ[i][j][1]) of.close() :: # The MIT License (MIT) # # Copyright (c) 2014 Muratahan Aykol # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE import numpy as np from copy import deepcopy # This function reads an XYZ file and a list of lattice vectors L = [x,y,z] and gives MSD + unwrapped coordinates def MSD(xyz_file,L): a = []; l = []; a.append(L[0]); a.append(L[1]); a.append(L[2]); #basis vectors in cartesian coords l.append(np.sqrt(np.dot(a[0],a[0]))); l.append(np.sqrt(np.dot(a[1],a[1]))); l.append(np.sqrt(np.dot(a[2],a[2]))); #basis vector lengths file = open(xyz_file, 'r') recorder = open("msd.out", 'w') coord_rec = open("unwrapped.xyz", 'w') origin_list = [] # Stores the origin as [element,[coords]] prev_list = [] # Stores the wrapped previous step unwrapped_list = [] # Stores the instantenous unwrapped msd = [] #Stores atom-wise MSD Stores msd as [msd] msd_dict ={} #Stores element-wise MSD msd_lattice = [] msd_dict_lattice ={} element_list = [] # element list element_dict = {} # number of elements stored content = file.readline() N = int(content) for i in range(N): msd.append(np.float64('0.0')) msd_lattice.append([0.0, 0.0, 0.0 ]) file.readline() step = 0 while True: step += 1 # Get and store the origin coordinates in origin_dict at first step if step == 1: for i in xrange(N): t = file.readline().rstrip('\n').split() element = t[0] if element not in element_list: element_list.append(element) if element not in element_dict: element_dict[element] = 1.0 else: element_dict[element] += 1.0 coords = np.array( [ float(s) for s in t[1:] ] ) origin_list.append([element,coords]) # Copy the first set of coordinates as prev_dict and unwrapped unwrapped_list = deepcopy(origin_list) prev_list = deepcopy(origin_list) recorder.write("step ") for element in element_list: recorder.write(element+" ") recorder.write("\n") # Read wrapped coordinates into wrapped_dict content = file.readline() if len(content) == 0: print "\n---End of file---\n" break N = int(content) file.readline() wrapped_list = [] # Erease the previous set of coordinates for i in xrange(N): t = file.readline().rstrip('\n').split() element = t[0] coords = np.array( [ float(s) for s in t[1:] ] ) wrapped_list.append([element,coords]) coord_rec.write(str(N)+ "\ncomment\n") # Unwrap coodinates and get MSD for atom in range(N): msd[atom] = 0.0 coord_rec.write(wrapped_list[atom][0]) # decompose wrapped atom coordinates to onto lattice vectors: w1 = wrapped_list[atom][1][0] w2 = wrapped_list[atom][1][1] w3 = wrapped_list[atom][1][2] # decompose prev atom coordinates to onto lattice vectors: p1 = prev_list[atom][1][0] p2 = prev_list[atom][1][1] p3 = prev_list[atom][1][2] #get distance between periodic images and use the smallest one if np.fabs(w1 - p1) > 0.5: u1 = w1 - p1 - np.sign(w1 - p1) else: u1 = w1 - p1 if np.fabs(w2 - p2) > 0.5: u2 = w2 - p2 - np.sign(w2 - p2) else: u2 = w2 - p2 if np.fabs(w3 - p3) > 0.5: u3 = w3 - p3 - np.sign(w3 - p3) else: u3 = w3 - p3 #add unwrapped displacements to unwrapped coords unwrapped_list[atom][1][0] += u1 unwrapped_list[atom][1][1] += u2 unwrapped_list[atom][1][2] += u3 uw = unwrapped_list[atom][1][0]*a[0] + unwrapped_list[atom][1][1]*a[1] +unwrapped_list[atom][1][2]*a[2] ol = origin_list[atom][1][0]*a[0] + origin_list[atom][1][1]*a[1] + origin_list[atom][1][2]*a[2] msd[atom] = np.linalg.norm(uw-ol)**2 msd_lattice[atom] = [np.linalg.norm(uw[0]-ol[0])**2,np.linalg.norm(uw[1]-ol[1])**2,np.linalg.norm(uw[2]-ol[2])**2] coord_rec.write(" " + np.array_str(uw).replace("[","").replace("]","")) coord_rec.write("\n") prev_list = [] # Store current wrapped coordinates for the next step prev_list = deepcopy(wrapped_list) # record msd recorder.write(str(step) + " ") for el in element_list: msd_dict[el] = 0.0 msd_dict_lattice[el]=[0.,0.,0.] for atom in range(len(msd)): msd_dict[wrapped_list[atom][0]] += msd[atom]/element_dict[wrapped_list[atom][0]] for i in range(3): msd_dict_lattice[wrapped_list[atom][0]][i] += msd_lattice[atom][i]/element_dict[wrapped_list[atom][0]] for el in element_list: recorder.write(str(msd_dict[el])+ " " + str(msd_dict_lattice[el][0])+ " " + str(msd_dict_lattice[el][1])+ " " + str(msd_dict_lattice[el][2])+ " ") recorder.write("\n") if step % 10 == 0: print step recorder.close() file.close() coord_rec.close() def read_lat_vec(): lat_file = open('lattice.vectors','r') line = [] for i in xrange(3): line.append([float(x) for x in lat_file.readline().rstrip('\n').split()]) print line[i] lattice = np.array([line[0],line[1],line[2]]) return lattice #You can read the lattice vectors from lattice.vector file lattice = read_lat_vec() #Or define the lattice vector manually as in #lattice =np.array([[-12.181156,-4.306689,7.459404],[0.000000,-12.920067,7.459404],[0.000000,0.000000,14.918808]]) #Run the MSD calculator with XDATCAR_fract.xyz and the lattice vector defined above MSD("XDATCAR_fract.xyz",lattice)