Skip to content

Commit

Permalink
Merge pull request #1603 from jmmshn/master
Browse files Browse the repository at this point in the history 
path_finder midpoint
  • Loading branch information
@shyuep
shyuep committed Oct 2, 2019
2 parents b6b6950 + bf9732b commit 32c1b82
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 52 deletions.
155 changes: 104 additions & 51 deletions pymatgen/analysis/path_finder.py
View file
Expand Up @@ -9,15 +9,15 @@
Ceder, The Journal of Chemical Physics 145 (7), 074112
"""

import numpy as np
import numpy.linalg as la
import scipy.signal
import scipy.stats
from scipy.interpolate import interp1d
from abc import ABCMeta
import math
import logging
from scipy.interpolate import interp1d
import scipy.signal
import scipy.stats
import numpy as np
import numpy.linalg as la

from abc import ABCMeta

from pymatgen.core.structure import Structure
from pymatgen.core.sites import PeriodicSite
Expand All @@ -34,7 +34,13 @@


class NEBPathfinder:
def __init__(self, start_struct, end_struct, relax_sites, v, n_images=20):
def __init__(self,
start_struct,
end_struct,
relax_sites,
v,
n_images=20,
mid_struct=None):
"""
General pathfinder for interpolating between two structures, where the
interpolating path is calculated with the elastic band method with
Expand All @@ -47,9 +53,11 @@ def __init__(self, start_struct, end_struct, relax_sites, v, n_images=20):
be relaxed
v: Static potential field to use for the elastic band relaxation
n_images: Number of interpolation images to generate
mid_struct: (optional) additional structure between the start and end structures to help
"""
self.__s1 = start_struct
self.__s2 = end_struct
self.__mid = mid_struct
self.__relax_sites = relax_sites
self.__v = v
self.__n_images = n_images
Expand All @@ -63,26 +71,46 @@ def interpolate(self):
pymatgen.core.structure.interpolate), and for the site indices given
in relax_sites, the path is relaxed by the elastic band method within
the static potential V.
If a mid point is defined we will interpolate from s1--> mid -->s2
The final number of images will still be n_images.
"""
images = self.__s1.interpolate(self.__s2, nimages=self.__n_images,
interpolate_lattices=False)
if self.__mid is not None:
# to make arithmatic easier we will do the interpolation in two parts with n images each
# then just take every other image at the end, this results in exactly n images
images_0 = self.__s1.interpolate(self.__mid,
nimages=self.__n_images,
interpolate_lattices=False)[:-1]
images_1 = self.__mid.interpolate(self.__s2,
nimages=self.__n_images,
interpolate_lattices=False)
images = images_0 + images_1
images = images[::2]
else:
images = self.__s1.interpolate(self.__s2,
nimages=self.__n_images,
interpolate_lattices=False)
for site_i in self.__relax_sites:
start_f = images[0].sites[site_i].frac_coords
end_f = images[-1].sites[site_i].frac_coords

path = NEBPathfinder.string_relax(
NEBPathfinder.__f2d(start_f, self.__v),
NEBPathfinder.__f2d(end_f, self.__v),
self.__v, n_images=(self.__n_images + 1),
dr=[self.__s1.lattice.a / self.__v.shape[0],
self.__v,
n_images=(self.__n_images + 1),
dr=[
self.__s1.lattice.a / self.__v.shape[0],
self.__s1.lattice.b / self.__v.shape[1],
self.__s1.lattice.c / self.__v.shape[2]])
self.__s1.lattice.c / self.__v.shape[2]
])
for image_i, image in enumerate(images):
image.translate_sites(site_i,
NEBPathfinder.__d2f(path[image_i],
self.__v) -
image.sites[site_i].frac_coords,
frac_coords=True, to_unit_cell=True)
image.translate_sites(
site_i,
NEBPathfinder.__d2f(path[image_i], self.__v) -
image.sites[site_i].frac_coords,
frac_coords=True,
to_unit_cell=True)
self.__images = images

@property
Expand All @@ -101,18 +129,27 @@ def plot_images(self, outfile):
sum_struct = self.__images[0].sites
for image in self.__images:
for site_i in self.__relax_sites:
sum_struct.append(PeriodicSite(image.sites[site_i].specie,
image.sites[site_i].frac_coords,
self.__images[0].lattice,
to_unit_cell=True,
coords_are_cartesian=False))
sum_struct.append(
PeriodicSite(image.sites[site_i].specie,
image.sites[site_i].frac_coords,
self.__images[0].lattice,
to_unit_cell=True,
coords_are_cartesian=False))
sum_struct = Structure.from_sites(sum_struct, validate_proximity=False)
p = Poscar(sum_struct)
p.write_file(outfile)

@staticmethod
def string_relax(start, end, V, n_images=25, dr=None, h=3.0, k=0.17,
min_iter=100, max_iter=10000, max_tol=5e-6):
def string_relax(start,
end,
V,
n_images=25,
dr=None,
h=3.0,
k=0.17,
min_iter=100,
max_iter=10000,
max_tol=5e-6):
"""
Implements path relaxation via the elastic band method. In general, the
method is to define a path by a set of points (images) connected with
Expand Down Expand Up @@ -191,12 +228,14 @@ def string_relax(start, end, V, n_images=25, dr=None, h=3.0, k=0.17,
# Calculate forces acting on string
d = V.shape
s0 = s
edV = np.array([[dV[0][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
int(pt[2]) % d[2]] / dr[0],
dV[1][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
int(pt[2]) % d[2]] / dr[0],
dV[2][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
int(pt[2]) % d[2]] / dr[0]] for pt in s])
edV = np.array([[
dV[0][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]]
/ dr[0],
dV[1][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]]
/ dr[0],
dV[2][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]]
/ dr[0]
] for pt in s])
# if(step % 100 == 0):
# logger.debug(edV)

Expand All @@ -206,10 +245,10 @@ def string_relax(start, end, V, n_images=25, dr=None, h=3.0, k=0.17,
ds_plus[0] = (ds_plus[0] - ds_plus[0])
ds_minus[-1] = (ds_minus[-1] - ds_minus[-1])
Fpot = edV
Fel = keff * (la.norm(ds_plus) - la.norm(ds0_plus)) * (
ds_plus / la.norm(ds_plus))
Fel += keff * (la.norm(ds_minus) - la.norm(ds0_minus)) * (
ds_minus / la.norm(ds_minus))
Fel = keff * (la.norm(ds_plus) -
la.norm(ds0_plus)) * (ds_plus / la.norm(ds_plus))
Fel += keff * (la.norm(ds_minus) -
la.norm(ds0_minus)) * (ds_minus / la.norm(ds_minus))
s -= h * (Fpot + Fel)

# Fix endpoints
Expand Down Expand Up @@ -246,19 +285,22 @@ def __f2d(frac_coords, v):
the grid size of v
"""
# frac_coords = frac_coords % 1
return np.array([int(frac_coords[0] * v.shape[0]),
int(frac_coords[1] * v.shape[1]),
int(frac_coords[2] * v.shape[2])])
return np.array([
int(frac_coords[0] * v.shape[0]),
int(frac_coords[1] * v.shape[1]),
int(frac_coords[2] * v.shape[2])
])

@staticmethod
def __d2f(disc_coords, v):
"""
Converts a point given in discrete coordinates withe respect to the
grid in v to fractional coordinates.
"""
return np.array([disc_coords[0] / v.shape[0],
disc_coords[1] / v.shape[1],
disc_coords[2] / v.shape[2]])
return np.array([
disc_coords[0] / v.shape[0], disc_coords[1] / v.shape[1],
disc_coords[2] / v.shape[2]
])


class StaticPotential(metaclass=ABCMeta):
Expand All @@ -269,6 +311,10 @@ class StaticPotential(metaclass=ABCMeta):
"""

def __init__(self, struct, pot):
"""
:param struct: atomic structure of the potential
:param pot: volumentric data to be used as a potential
"""
self.__v = pot
self.__s = struct

Expand Down Expand Up @@ -298,17 +344,23 @@ def rescale_field(self, new_dim):
"""
v_dim = self.__v.shape
padded_v = np.lib.pad(self.__v, ((0, 1), (0, 1), (0, 1)), mode='wrap')
ogrid_list = np.array([list(c) for c in list(
np.ndindex(v_dim[0] + 1, v_dim[1] + 1, v_dim[2] + 1))])
ogrid_list = np.array([
list(c)
for c in list(np.ndindex(v_dim[0] + 1, v_dim[1] + 1, v_dim[2] + 1))
])
v_ogrid = padded_v.reshape(
((v_dim[0] + 1) * (v_dim[1] + 1) * (v_dim[2] + 1), -1))
ngrid_a, ngrid_b, ngrid_c = (np.mgrid[0: v_dim[0]: v_dim[0] / new_dim[0], 0: v_dim[1]: v_dim[1] / new_dim[1],
0: v_dim[2]: v_dim[2] / new_dim[2]])
ngrid_a, ngrid_b, ngrid_c = (np.mgrid[0:v_dim[0]:v_dim[0] /
new_dim[0], 0:v_dim[1]:v_dim[1] /
new_dim[1], 0:v_dim[2]:v_dim[2] /
new_dim[2]])

v_ngrid = scipy.interpolate.griddata(ogrid_list, v_ogrid,
v_ngrid = scipy.interpolate.griddata(ogrid_list,
v_ogrid,
(ngrid_a, ngrid_b, ngrid_c),
method='linear').reshape(
(new_dim[0], new_dim[1], new_dim[2]))
(new_dim[0], new_dim[1],
new_dim[2]))
self.__v = v_ngrid

def gaussian_smear(self, r):
Expand Down Expand Up @@ -346,13 +398,14 @@ def gaussian_smear(self, r):
1.0):
g = np.array(
[g_a / v_dim[0], g_b / v_dim[1], g_c / v_dim[2]]).T
gauss_dist[int(g_a + r_disc[0])][int(g_b + r_disc[1])][
int(g_c + r_disc[2])] = la.norm(
np.dot(self.__s.lattice.matrix, g)) / r
gauss_dist[int(g_a + r_disc[0])][int(g_b + r_disc[1])][int(
g_c + r_disc[2])] = la.norm(
np.dot(self.__s.lattice.matrix, g)) / r
gauss = scipy.stats.norm.pdf(gauss_dist)
gauss = gauss / np.sum(gauss, dtype=float)
padded_v = np.pad(self.__v, (
(r_disc[0], r_disc[0]), (r_disc[1], r_disc[1]), (r_disc[2], r_disc[2])),
padded_v = np.pad(self.__v,
((r_disc[0], r_disc[0]), (r_disc[1], r_disc[1]),
(r_disc[2], r_disc[2])),
mode='wrap')
smeared_v = scipy.signal.convolve(padded_v, gauss, mode='valid')
self.__v = smeared_v
Expand Down
11 changes: 10 additions & 1 deletion pymatgen/analysis/tests/test_path_finder.py
View file
Expand Up @@ -7,7 +7,7 @@
from pymatgen.analysis.path_finder import NEBPathfinder, ChgcarPotential
from pymatgen.io.vasp import Poscar, Chgcar
from pymatgen.core.periodic_table import Element

from numpy import mean
__author__ = 'Ziqin (Shaun) Rong'
__version__ = '0.1'
__maintainer__ = 'Ziqin (Shaun) Rong'
Expand All @@ -24,6 +24,8 @@ def test_image_num(self):
"path_finder")
start_s = Poscar.from_file(os.path.join(test_file_dir, 'LFP_POSCAR_s')).structure
end_s = Poscar.from_file(os.path.join(test_file_dir, 'LFP_POSCAR_e')).structure
mid_s = start_s.interpolate(end_s, nimages=2,
interpolate_lattices=False)[1]
chg = Chgcar.from_file(os.path.join(test_file_dir, 'LFP_CHGCAR.gz'))
moving_cation_specie = Element('Li')
relax_sites = []
Expand All @@ -38,6 +40,13 @@ def test_image_num(self):
images.append(image)
self.assertEqual(len(images), 9)

pf_mid = NEBPathfinder(start_s, end_s, relax_sites=relax_sites,
v=ChgcarPotential(chg).get_v(), n_images=10, mid_struct=mid_s)
moving_site = relax_sites[0]
dists = [s1.sites[moving_site].distance(s2.sites[moving_site])
for s1, s2 in zip(pf.images[:-1], pf.images[1:])]
# check that all the small distances are about equal
self.assertTrue(abs(min(dists)-max(dists))/mean(dists) < 0.02)

if __name__ == '__main__':
unittest.main()

0 comments on commit 32c1b82

Please sign in to comment.