A series of Python3 script to lower the barrier of computing and simulating molecular and material systems.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

69 lines
2.2 KiB

2 years ago
import warnings
warnings.filterwarnings("ignore")
import MDAnalysis
import MDAnalysis.analysis.rdf
import numpy as np
from pymatgen.core import lattice
import builtins
import os
def rdf(traject,latt,pair,endrange,resolution):
cell = []
with open(latt, 'r') as f:
for line in f:
if 'TV' in line:
cell.append(list(map(float, line.split()[1:4])))
box = lattice.Lattice(cell)
(lx, ly, lz), (a, b, c) = box.lengths_and_angles
u = MDAnalysis.Universe(traject, format='XYZ')
u.dimensions = np.array([lx, ly, lz, a, b, c], dtype=np.float32)
builtins.u = u
ag1 = None
ag2 = None
elems = pair.split('-')
# Definisi file keluaran
outfile = open('rdf_{}-{}.dat'.format(elems[0],elems[1]), 'w')
ag1 = u.select_atoms(f'name {elems[0]}')
ag2 = u.select_atoms(f'name {elems[1]}')
nbins = int((endrange - 0.05)/resolution)
rdf = MDAnalysis.analysis.rdf.InterRDF(ag1, ag2, nbins=nbins, range=(0.05, endrange), verbose=True)
rdf.run()
vol = np.power(rdf.edges[1:], 3) - np.power(rdf.edges[:-1], 3)
vol *= 4/3.0 * np.pi
box_vol = box.volume
pair_den = box_vol / (len(ag1.atoms) * len(ag2.atoms))
integral = rdf.rdf / (pair_den / vol)
total = 0.0
step = 0
for r, v, d in zip(rdf.bins, rdf.rdf, integral):
step += 1
total += d/len(ag1.atoms)
if (step == 1):
continue
print('{:8.4f} {:8.4f} {:8.4f}'.format(r, v, total), file=outfile)
outfile.close()
with open('rdf_{}-{}.gp'.format(elems[0],elems[1]),'w') as fout:
print("""set terminal pdf
set output "rdf_{}-{}.pdf"
set key
set yrange [0:10]
set ytics (0,1,2,3,4,5,6,7,8,9,10)
set xlabel "Distance [Angstroms]"
set ylabel "RDF [arb.]"
set style line 1 lc rgb '#00008b' lt 1 lw 2 pt 7
set style line 2 lc rgb '#FF0000' lt 1 lw 2 pt 7
plot 'rdf_{}-{}.dat' u 1:2 with line title "RDF" linestyle 1,\
'' u 1:3 with line title "Coordination Number" linestyle 2""".format(elems[0],elems[1],elems[0],elems[1]),file=fout)
os.system('gnuplot rdf_{}-{}.gp'.format(elems[0],elems[1]))