mirror of https://github.com/QMCPACK/qmcpack.git
nexus: add plotting capability for NL PP operator
This commit is contained in:
Jaron Krogel
parent
6c7c3ff755
commit
0905b49367
|
|
@ -2358,16 +2358,20 @@ class SpheroidGrid(StructuredGridWithAxes):
|
|||
elif cells is not None:
|
||||
grid_dim = len(cells)
|
||||
#end if
|
||||
|
||||
# force bconds to match sphere constraints
|
||||
if grid_dim==3:
|
||||
bconds = tuple('oop')
|
||||
elif grid_dim==2:
|
||||
bconds = tuple('op')
|
||||
else:
|
||||
if grid_dim not in (2,3):
|
||||
self.error('only 2 and 3 dimensional spheroids grids are supported\nrequested dimension: {}'.format(grid_dim))
|
||||
#end if
|
||||
|
||||
# bconds match sphere constraints by default
|
||||
bconds = kwargs.get('bconds',None)
|
||||
if bconds is None:
|
||||
if grid_dim==3:
|
||||
bconds = tuple('oop')
|
||||
elif grid_dim==2:
|
||||
bconds = tuple('op')
|
||||
#end if
|
||||
#end if
|
||||
|
||||
endpoint = self.has_endpoints(bconds=bconds,grid_dim=grid_dim)
|
||||
|
||||
points,shape = spheroid_grid_points(axes,shape=shape,cells=cells,centered=centered,endpoint=endpoint,return_shape=True)
|
||||
|
|
@ -3975,8 +3979,8 @@ if __name__=='__main__':
|
|||
|
||||
if demos.plot_surface:
|
||||
|
||||
gp = ParallelotopeGrid(
|
||||
axes = [[1,0],
|
||||
gp = ParallelotopeGrid(
|
||||
axes = [[1,0],
|
||||
[1,1]],
|
||||
bconds = 'pp',
|
||||
cells = (80,80),
|
||||
|
|
|
|||
|
|
@ -39,7 +39,9 @@
|
|||
import os
|
||||
from subprocess import Popen
|
||||
from execute import execute
|
||||
import numpy as np
|
||||
from numpy import linspace,array,zeros,append,mgrid,empty,exp,minimum,maximum,sqrt,arange
|
||||
import matplotlib.pyplot as plt
|
||||
from fileio import TextFile
|
||||
from xmlreader import readxml
|
||||
from superstring import string2val,split_delims
|
||||
|
|
@ -1241,6 +1243,171 @@ class SemilocalPP(Pseudopotential):
|
|||
#end def plot_L2
|
||||
|
||||
|
||||
def plot_nonlocal_polar(self,show=True,lmax=10,rmin=0.01,rmax=2.0,nr=100,nt=100,levels=100,label=''):
|
||||
from scipy.special import eval_legendre as legendre
|
||||
|
||||
tlabel = label
|
||||
rpow = 1
|
||||
|
||||
# update rcut, if needed
|
||||
if self.rcut is None:
|
||||
self.update_rcut()
|
||||
#end if
|
||||
rc = self.rcut
|
||||
|
||||
# get legendre polynomials
|
||||
th = np.linspace(0.0,2*np.pi,nt)
|
||||
cos = np.cos(th)
|
||||
sin = np.sin(th)
|
||||
Pl = obj()
|
||||
Pli = obj()
|
||||
for li in range(lmax):
|
||||
pl_i = (2*li+1)/(4*np.pi)*legendre(li,cos)
|
||||
Pli[li] = pl_i
|
||||
if li<len(self.l_channels):
|
||||
l = self.l_channels[li]
|
||||
Pl[l] = pl_i
|
||||
#end if
|
||||
#end for
|
||||
|
||||
# set the colormap and centre the colorbar
|
||||
import matplotlib.colors as colors
|
||||
class MidNorm(colors.Normalize):
|
||||
def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
|
||||
self.midpoint = midpoint
|
||||
colors.Normalize.__init__(self, vmin, vmax, clip)
|
||||
#end def __init__
|
||||
def __call__(self, value, clip=None):
|
||||
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
|
||||
return np.ma.masked_array(np.interp(value, x, y), np.isnan(value))
|
||||
#end def __call__
|
||||
#end class MidNorm
|
||||
|
||||
# select a subset of the radial points
|
||||
r = None
|
||||
vl = obj()
|
||||
vnl = self.get_nonlocal()
|
||||
if self.numeric and not self.interpolatable:
|
||||
for l in vnl.keys():
|
||||
rf,vf = self.evaluate_nonlocal(l=l,rpow=rpow,rret=True)
|
||||
if r is None:
|
||||
drf = rf[1]-rf[0]
|
||||
dr = rmax/nr
|
||||
ndrf = int(np.round(np.ceil(dr/drf)))
|
||||
dr = drf*ndrf
|
||||
rmax = nr*dr + 0.1*drf
|
||||
rng = rf<rmax
|
||||
r = rf[rng][::ndrf]
|
||||
#end if
|
||||
vl[l] = vf[rng][::ndrf]
|
||||
#end for
|
||||
else:
|
||||
self.error('plot_polar does not yet support non-numeric potentials')
|
||||
#end if
|
||||
|
||||
# plot the radial potentials
|
||||
figure()
|
||||
vmin = 1e99
|
||||
vmax = -1e99
|
||||
for l in self.l_channels:
|
||||
if l in vl:
|
||||
color = self.channel_colors[l]
|
||||
v = vl[l]
|
||||
plot(r,v,color+'-',label='v'+l)
|
||||
vmin = min(v.min(),vmin)
|
||||
vmax = max(v.max(),vmax)
|
||||
#end if
|
||||
#end for
|
||||
plot([rc,rc],[vmin,vmax],'k--',lw=2)
|
||||
xlim([-0.1*rc,1.1*rc])
|
||||
xlabel('r (Bohr)')
|
||||
ylabel('V NL (Ha)')
|
||||
title((tlabel+' NL channels').strip())
|
||||
legend()
|
||||
|
||||
# function for a single polar plot
|
||||
def plot_V(V,label):
|
||||
vmin = V.min()
|
||||
vmax = V.max()
|
||||
vm = max(np.abs(vmin),np.abs(vmax))
|
||||
vmin = -vm
|
||||
vmax = vm
|
||||
lev = linspace(vmin,vmax,levels)
|
||||
|
||||
fig = figure(tight_layout=True)
|
||||
ax = fig.add_subplot(111)
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('z')
|
||||
ax.set_aspect('equal','box')
|
||||
fig.tight_layout()
|
||||
|
||||
xlim(lim)
|
||||
ylim(lim)
|
||||
|
||||
cmap = plt.cm.get_cmap('seismic')
|
||||
|
||||
mid_norm = MidNorm(vmin,vmax,0.0)
|
||||
|
||||
cs = ax.contourf(X,Z,V,levels=lev,cmap=cmap,clim=(vmin,vmax),norm=mid_norm)
|
||||
plot(rc*cos,rc*sin,'k--',lw=2)
|
||||
|
||||
fig.colorbar(cs, ax=ax, shrink=0.9)
|
||||
|
||||
title((tlabel+' V {}'.format(label)).strip())
|
||||
#end def plot_V
|
||||
|
||||
# make a polar plot of each non-local component
|
||||
R,COS = np.array(np.meshgrid(r,cos,indexing='ij'))
|
||||
Z = R*COS
|
||||
R,SIN = np.array(np.meshgrid(r,sin,indexing='ij'))
|
||||
X = R*SIN
|
||||
lim = [-1.1*rc,1.1*rc]
|
||||
VSUM = None
|
||||
for l in self.l_channels:
|
||||
if l in vl:
|
||||
VL,PL = np.array(np.meshgrid(vl[l],Pl[l],indexing='ij'))
|
||||
V = VL*PL
|
||||
|
||||
if VSUM is None:
|
||||
VSUM = V.copy()
|
||||
else:
|
||||
VSUM += V
|
||||
#end if
|
||||
|
||||
plot_V(V,'NL '+l)
|
||||
#end if
|
||||
#end for
|
||||
|
||||
# plot the total non-local operator
|
||||
plot_V(VSUM,'NL sum')
|
||||
|
||||
# plot approximation to L2 operator, if present
|
||||
if self.has_L2():
|
||||
vL2 = self.evaluate_L2(rpow=rpow)
|
||||
vL2 = vL2[rng][::ndrf]
|
||||
VL2SUM = None
|
||||
for li in range(1,lmax):
|
||||
print('V L2 {} of {}'.format(li,lmax))
|
||||
v = li*(li+1)*vL2
|
||||
VL,PL = np.array(np.meshgrid(v,Pli[li],indexing='ij'))
|
||||
V = VL*PL
|
||||
if VL2SUM is None:
|
||||
VL2SUM = V.copy()
|
||||
else:
|
||||
VL2SUM += V
|
||||
#end if
|
||||
#plot_V(V,'L2 '+str(li))
|
||||
#end for
|
||||
plot_V(VL2SUM,'L2 sum (Lmax={})'.format(lmax))
|
||||
#end if
|
||||
|
||||
if show:
|
||||
show_plots()
|
||||
#end if
|
||||
|
||||
#end def plot_nonlocal_polar
|
||||
|
||||
|
||||
def write_qmcpack(self,filepath=None):
|
||||
self.update_rcut(tol=1e-5,optional=True)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue