##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# pseudopotential.py #
# Classes for reading pseudopotential data and converting #
# pseudopotentials between file formats. #
# #
# Content summary: #
# Pseudopotentials #
# Class contains a list of all pseudopotentials available if #
# user provides pseudo_dir in settings. #
# #
# PseudoFile #
# Represents a single pseudopotential file. No functionality. #
# #
# gamessPPFile #
# Represents a pseudopotential file for GAMESS. #
# Used during Nexus execution to separate basis and channel #
# information to fill in the input file. #
# #
# Classes below are user-facing, but are not used by Nexus itself.#
# #
# Pseudopotential #
# Represents a generic pseudopotential. #
# #
# SemilocalPP #
# Represents a semi-local pseudopotential. #
# Contains data for each non-local channel. #
# Supports plotting, gaussian fitting, and writing to QMCPACK's #
# grid-based file format. #
# GaussianPP and QmcpackPP are derived classes to read/write #
# GAMESS/Gaussian and QMCPACK pseudopotential files. #
# #
#====================================================================#
import os
from pathlib import Path
import numpy as np
from .execute import execute
from .fileio import TextFile
from .xmlreader import readxml
from .periodic_table import Elements
from .unit_converter import convert
from .developer import DevBase, obj, unavailable, error
from .basisset import process_gaussian_text, GaussianBasisSet
from .physical_system import PhysicalSystem
from .testing import object_eq
from .utilities import path_string, is_valid_filename
try:
import matplotlib.pyplot as plt
except:
plt = unavailable('matplotlib','pyplot')
#end try
[docs]
def pp_elem_label(filename,guard=False):
if guard and not is_valid_filename(filename):
error(f"Pseudopotential file name {filename} is invalid!")
el = ''
for c in filename:
if c=='.' or c=='_' or c=='-':
break
#end if
el+=c
#end for
elem_label = el
is_elem, element = Elements.is_element(el, return_element=True)
if guard:
if not is_elem:
error(
'cannot determine element for pseudopotential file: {0}\n'
'pseudopotential file names must be prefixed by an atomic symbol or label\n'
'(e.g. Si, Si1, etc)'.format(filename)
)
#end if
return elem_label, element.symbol
else:
if isinstance(element, Elements):
return elem_label, element.symbol, is_elem
else:
return elem_label, element, is_elem
#end if
#end def pp_elem_label
# basic interface for nexus, only gamess really needs this for now
[docs]
class PseudoFile(DevBase):
def __init__(self,filepath=None):
self.element = None
self.element_label = None
self.filename = None
self.location = None
if filepath is not None:
self.filename = os.path.basename(filepath)
self.location = os.path.abspath(filepath)
elem_label,symbol,is_elem = pp_elem_label(self.filename)
if not is_elem:
self.error('cannot determine element for pseudopotential file: {0}\npseudopotential file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(filepath))
#end if
self.element = symbol
self.element_label = elem_label
self.read(filepath)
#end if
#end def __init__
[docs]
def read(self,filepath):
None
#end def read
#end class PseudoFile
[docs]
class gamessPPFile(PseudoFile):
def __init__(self,filepath=None):
self.pp_text = None
self.pp_name = None
self.basis_text = None
PseudoFile.__init__(self,filepath)
#end def __init__
[docs]
def read(self,filepath):
with open(filepath, "r") as f:
lines = f.read().splitlines()
new_block = True
tokens = []
block = ''
nline = 0
for line in lines:
nline+=1
ls = line.strip()
if len(ls)>0 and ls[0]!='!' and ls[0]!='#':
if new_block:
tokens = ls.split()
new_block = False
if len(tokens)!=5:
block+=line+'\n'
#end if
else:
block+=line+'\n'
#end if
#end if
if (len(ls)==0 or nline==len(lines)) and len(block)>0:
block = block.rstrip()
if len(tokens)==4:
self.pp_text = block
self.pp_name = tokens[0]
elif len(tokens)==5:
self.basis_text = block
else:
self.error('could not identify text block in {0} as pseudopotential or basis text\btext block:\n{1}'.format(self.filename,block))
#end if
new_block = True
tokens = []
block = ''
#end if
#end for
if self.pp_text is None:
self.error('could not find pseudopotential text in '+self.filename)
#end if
if self.basis_text is None:
self.error('could not find basis text in '+self.filename)
#end if
#end def read
#end class gamessPPFile
[docs]
class Pseudopotentials(DevBase):
def __init__(self,*pseudopotentials):
if len(pseudopotentials)==1 and isinstance(pseudopotentials[0],list):
pseudopotentials = pseudopotentials[0]
#end if
ppfiles = []
pps = []
errors = False
for pp in pseudopotentials:
if isinstance(pp,PseudoFile):
pps.append(pp)
elif isinstance(pp, str | Path):
ppfiles.append(path_string(pp))
else:
self.error('expected PseudoFile type or filepath, got '+str(type(pp)),exit=False)
errors = True
#end if
#end for
if errors:
self.error('cannot create Pseudopotentials object')
#end if
if len(pps)>0:
self.addpp(pps)
#end if
if len(ppfiles)>0:
self.readpp(ppfiles)
#end if
#end def __init__
[docs]
def addpp(self,*pseudopotentials):
if len(pseudopotentials)==1 and isinstance(pseudopotentials[0],list):
pseudopotentials = pseudopotentials[0]
#end if
for pp in pseudopotentials:
self[pp.filename] = pp
#end for
#end def addpp
[docs]
def readpp(self,*ppfiles):
if len(ppfiles)==1 and isinstance(ppfiles[0],list):
ppfiles = ppfiles[0]
#end if
pps = []
self.log('\n Pseudopotentials')
for filepath in ppfiles:
filename = os.path.basename(filepath)
elem_label,symbol,is_elem = pp_elem_label(filename)
is_file = os.path.isfile(filepath)
if is_elem and is_file:
self.log(' reading pp: ',filepath)
ext = filepath.split('.')[-1].lower()
if ext=='gms':
pp = gamessPPFile(filepath)
else:
pp = PseudoFile(filepath)
#end if
pps.append(pp)
elif not is_file:
self.log(' ignoring directory: ',filepath)
elif not is_elem:
self.log(' ignoring file w/o atomic symbol: ',filepath)
#end if
#end for
self.log(' ')
self.addpp(pps)
#end def readpp
[docs]
def pseudos_by_atom(self,*ppfiles):
pps = obj()
for ppfile in ppfiles:
if ppfile in self:
pp = self[ppfile]
pps[pp.element_label] = pp
else:
self.error('pseudopotential file not found\nmissing file: {0}'.format(ppfile))
#end if
#end for
return pps
#end def pseudos_by_atom
#end class Pseudopotentials
# user interface to group sets of pseudopotentials together and refer to them by labels
# labeling should eliminate the need to provide lists of pseudopotential files to each
# simulation object (e.g. via a generate_* call) separately
[docs]
class PPset(DevBase):
instance_counter = 0
known_codes = set('pwscf gamess vasp qmcpack'.split())
default_extensions = obj(
pwscf = ['ncpp','upf'],
gamess = ['gms'],
vasp = ['potcar'],
qmcpack = ['xml'],
)
def __init__(self):
if PPset.instance_counter!=0:
self.error('cannot instantiate more than one PPset object\nintended use follows a singleton pattern')
#end if
PPset.instance_counter+=1
self.pseudos = obj()
#end def __init__
[docs]
def supports_code(self,code):
return code in PPset.known_codes
#end def supports_code
def __call__(self,label,**code_pps):
if not isinstance(label,str):
self.error('incorrect use of ppset\nlabel provided must be a string\nreceived type instead: {0}\nwith value: {1}'.format(label.__class__.__name__,label))
#end if
if label in self.pseudos:
self.error('incorrect use of ppset\npseudopotentials with label "{0}" have already been added to ppset'.format(label))
#end if
pseudos = obj()
self.pseudos[label]=pseudos
for code,pps in code_pps.items():
clow = code.lower()
if clow not in self.known_codes:
self.error('incorrect use of ppset\ninvalid simulation code "{0}" provided with set labeled "{1}"\nknown simulation codes are: {2}'.format(code,label,sorted(self.known_codes)))
#end if
if not isinstance(pps,(list,tuple)):
self.error('incorrect use of ppset\nmust provide a list of pseudopotentials for code "{0}" in set labeled "{1}"\ntype provided instead of list: {2}'.format(code,label,pps.__class__.__name__))
#end if
ppcoll = obj()
for pp in pps:
if not isinstance(pp, (str, Path)):
self.error('incorrect use of ppset\nnon-filename provided with set labeled "{0}" for simulation code "{1}"\neach pseudopential file name must be a string\nreceived type: {2}\nwith value: {3}'.format(label,code,pp.__class__.__name__,pp))
else:
pp = path_string(pp)
elem_label, symbol, is_elem = pp_elem_label(pp)
if not is_elem:
self.error('invalid filename provided to ppset\ncannot determine element for pseudopotential file: {0}\npseudopotential file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(pp))
elif symbol in ppcoll:
self.error('incorrect use of ppset\nmore than one pseudopotential file provided for element "{0}" for code "{1}" in set labeled "{2}"\nfirst file: {3}\nsecond file: {4}'.format(symbol,code,label,ppcoll[symbol],pp))
#end if
ppcoll[symbol] = path_string(pp)
#end for
pseudos[clow] = ppcoll
#end for
#end def __call__
[docs]
def has_set(self,label):
return label in self.pseudos
#end def has_set
# test needed
[docs]
def get(self,label,code,system):
if system is None or not system.pseudized:
return []
#end if
if not isinstance(system,PhysicalSystem):
self.error('system object must be of type PhysicalSystem')
#end if
species_labels,species = system.structure.species(symbol=True)
if not isinstance(label,str):
self.error('incorrect use of ppset\nlabel provided must be a string\nreceived type instead: {0}\nwith value: {1}'.format(label.__class__.__name__,label))
#end if
if not self.has_set(label):
self.error('incorrect use of ppset\npseudopotential set labeled "{0}" is not present in ppset\nset labels present: {1}\nplease either provide pseudopotentials with label "{0}" or correct the provided label'.format(label,sorted(self.pseudos.keys())))
#end if
pseudos = self.pseudos[label]
clow = code.lower()
if clow not in self.known_codes:
self.error('simulation code "{0}" is not known to ppset\nknown codes are: {1}'.format(code,sorted(self.known_codes)))
elif clow not in pseudos:
self.error('incorrect use of ppset\npseudopotentials were not provided for simulation code "{0}" in set labeled "{1}"\npseudopotentials are required for physical system with pseudo-elements: {2}\nplease add these pseudopotentials for code "{0}" in set "{1}"'.format(code,label,sorted(species)))
#end if
ppcoll = pseudos[clow]
pps = []
for symbol in species:
if symbol not in ppcoll:
self.error('incorrect use of ppset\npseudopotentials were not provided for element "{0}" code "{1}" in set labeled "{2}"\nphysical system encountered with pseudo-elements: {3}\nplease ensure that pseudopotentials are provided for these elements in set "{2}" for code "{1}"'.format(symbol,code,label,sorted(species)))
#end if
pps.append(ppcoll[symbol])
#end for
return pps
#end def get
#end class PPset
ppset = PPset()
# real pseudopotentials
[docs]
class Pseudopotential(DevBase):
requires_format = False
formats = None
def __init__(self,filepath=None,format=None):
self.element = None
self.core = None
self.Zval = None
self.Zcore = None
if filepath is not None:
self.read(filepath,format)
#end if
#end def __init__
[docs]
def transfer_core_from(self,other):
self.element = other.element
self.core = other.core
self.Zval = other.Zval
self.Zcore = other.Zcore
#end def transfer_core_from
[docs]
def read(self,filepath,format=None):
filepath = path_string(filepath)
if self.requires_format:
if format is None:
self.error('format keyword must be specified to read file {0}\nvalid options are: {1}'.format(filepath,self.formats))
elif format not in self.formats:
self.error('incorrect format requested: {0}\nvalid options are: {1}'.format(format,self.formats))
#end if
#end if
if not os.path.exists(filepath):
self.error('cannot read {0}, file does not exist'.format(filepath))
#end if
self.element = pp_elem_label(os.path.split(filepath)[1])[0]
with open(filepath, "r") as f:
text = f.read()
self.read_text(text,format,filepath=filepath)
#end def read
[docs]
def write(self,filepath=None,format=None):
if self.requires_format:
if format is None:
self.error('format keyword must be specified to write file {0}\nvalid options are: {1}'.format(filepath,self.formats))
elif format not in self.formats:
self.error('incorrect format requested: {0}\nvalid options are: {1}'.format(format,self.formats))
#end if
#end if
text = self.write_text(format)
if filepath is not None:
with open(filepath, "w") as f:
f.write(text)
#end if
return text
#end def write
[docs]
def read_text(self,text,format=None,filepath=None):
self.not_implemented()
#end def read_text
[docs]
def write_text(self,format=None):
self.not_implemented()
#end def write_text
[docs]
def convert(self,format):
self.not_implemented()
#end def convert
[docs]
def plot(self,r=None,show=True):
self.not_implemented()
#end def plot
#end class Pseudopotential
[docs]
class SemilocalPP(Pseudopotential):
l_channels = tuple('spdfghiklmnoqrtuvwxyz')
channel_colors = obj(s='g',p='r',d='b',f='m',g='c',h='k',i='g',L2='k')
numeric = False
interpolatable = True
formats = ['qmcpack','casino']
channel_indices = obj()
for i,c in enumerate(l_channels):
channel_indices[c] = i
#end for
channel_indices.L2 = -1
def __init__(self,filepath=None,format=None,name=None,src=None):
self.name = name
self.rcut = None
self.lmax = None
self.local = None
if self.numeric:
self.r = None
#end if
self.rcut_L2 = None
self.components = obj()
Pseudopotential.__init__(self,filepath,format)
if src is not None:
self.transfer_core_from(src)
#end if
#end def __init__
# test needed
[docs]
def transfer_core_from(self,other):
self.name = other.name
self.rcut = other.rcut
self.lmax = other.lmax
self.local = other.local
if self.numeric and other.numeric:
self.r = other.r
#end if
Pseudopotential.transfer_core_from(self,other)
#end def transfer_core_from
[docs]
def read(self,filepath,format=None):
Pseudopotential.read(self,filepath,format)
if self.rcut is None:
self.update_rcut()
#end if
#end def read
[docs]
def has_component(self,l):
return l in self.components
#end def has_component
# test needed
[docs]
def set_component(self,l,v,guard=False):
if guard and l in self.components:
self.error('cannot set requested component potential\nrequested potential is already present\nrequested potential: {0}'.format(l))
#end if
self.components[l] = v
#end def set_component
[docs]
def get_component(self,l,guard=False):
v = None
if l in self.components:
v = self.components[l]
elif guard:
self.error('cannot get requested component potential\nrequested potential is not present\nrequested potential: {0}\npotentials present: {1}'.format(l,list(self.components.keys())))
#end if
return v
#end def get_component
# test needed
[docs]
def remove_component(self,l,guard=False):
if l in self.components:
del self.components[l]
elif guard:
self.error('cannot remove requested component potential\nrequested potential is not present\nrequested potential: {0}\npotentials present: {1}'.format(l,list(self.components.keys())))
#end if
#end def remove_component
[docs]
def has_local(self):
return self.local in self.components
#end def has_local
[docs]
def has_nonlocal(self,l=None):
vnl = self.get_nonlocal()
if l is None:
return len(vnl)>0
else:
return l in vnl
#end if
#end def has_nonlocal
[docs]
def has_L2(self):
return 'L2' in self.components
#end def has_L2
[docs]
def get_local(self):
return self.get_component(self.local,guard=True)
#end def get_local
[docs]
def get_nonlocal(self,l=None):
vnl = obj()
for lc,vc in self.components.items():
if lc!=self.local and lc!='L2':
vnl[lc] = vc
#end if
#end for
if l is None:
return vnl
elif l in vnl:
return vnl[l]
else:
self.error('cannot get nonlocal potential\nrequested potential is not nonlocal\nrequested potential: {0}\nnonlocal potentials present: {1}'.format(l,list(vnl.keys())))
#end if
#end def get_nonlocal
# test needed
[docs]
def get_L2(self):
return self.get_component('L2',guard=True)
#end def get_L2
# test needed
[docs]
def add_local(self,l,v):
self.set_component(l,v,guard=True)
self.local = l
#end def add_local
# test needed
[docs]
def add_nonlocal(self,l,v):
if l==self.local:
self.promote_local()
#end if
self.set_component(l,v,guard=True)
#end def add_nonlocal
# test needed
[docs]
def add_L2(self,v):
self.set_component('L2',v,guard=True)
#end def add_L2
# test needed
[docs]
def remove_local(self):
self.remove_component(self.local,guard=True)
self.local = None
#end def remove_local
# test needed
[docs]
def remove_nonlocal(self,l=None):
vnl = self.get_nonlocal()
if l is None:
for l in vnl.keys():
self.remove_component(l,guard=True)
#end for
elif l in vnl:
self.remove_component(l,guard=True)
else:
self.error('cannot remove nonlocal potential\nrequested potential is not present\nrequested potential: {0}\nnonlocal potentials present: {1}'.format(l,list(vnl.keys())))
#end if
#end def remove_nonlocal
# test needed
[docs]
def remove_L2(self):
self.remove_component('L2',guard=True)
#end def remove_L2
[docs]
def assert_numeric(self,loc):
if not self.numeric:
self.error('failing at {0}\n{0} is only supported for numerical pseudopotential formats'.format(loc))
#end if
#end def assert_numeric
# test needed
[docs]
def change_local(self,local):
self.assert_numeric('change_local')
if local==self.local:
return
#end if
if self.has_component('L2'):
self.error('cannot change local channel\nL2 potential is present')
elif not self.has_component(self.local):
self.error('cannot change local potential\ncurrent local potential is not present\ncurrent local potential: {0}\npotentials present: {1}'.format(self.local,list(self.components.keys())))
elif not self.has_component(local):
self.error('cannot change local potential\nrequested local potential is not present\nrequested local potential: {0}\npotentials present: {1}'.format(local,list(self.components.keys())))
#end if
vcs = self.components
vloc = vcs[self.local]
# get the l channels
vls = obj()
for l,vc in self.components.items():
if l==self.local:
vls[l] = vloc
else:
vls[l] = vc+vloc
#end if
#end for
# from l channels, reconstruct nonlocal components
vcs.clear()
vloc = vls[local]
for l,vl in vls.items():
if l==local:
vcs[l] = vloc
else:
vcs[l] = vls[l]-vloc
#end if
#end for
self.local = local
#end def change_local
# test needed
#end def promote_local
# test needed
# ensure that v_l=<l|vpp|l>==v, while v_l' remains unchanged
[docs]
def set_channel(self,l,v):
self.assert_numeric('set_channel')
if not self.has_local():
self.error('cannot enforce channel matching via set_channel\nthe local potential is missing and must be present\nrequested channel: {0}'.format(l))
#end if
if l==self.local:
self.promote_local()
#end if
vloc = self.get_local()
vnl = self.get_nonlocal()
if self.has_L2():
vL2 = self.get_L2()
else:
vL2 = 0*vloc
#end if
li = self.channel_indices[l]
vl = vloc + li*(li+1)*vL2
self.components[l] = v-vl
#end def set_channel
# test needed
[docs]
def expand_L2(self,lmax):
self.assert_numeric('expand_L2')
if lmax not in self.channel_indices:
self.error('cannot expand L2 up to angular momentum "{0}"\nvalid options for lmax are: {1}'.format(lmax,self.l_channels))
#end if
limax = self.channel_indices[lmax]
vps = obj()
for li in range(limax+1):
l = self.l_channels[li]
vps[l] = self.evaluate_channel(l=l,rpow=1)
#end for
if self.has_L2():
self.remove_L2()
#end if
self.components[self.local] = vps[self.local]
for l,v in vps.items():
if l!=self.local:
self.set_channel(l,v)
#end if
#end for
#end def expand_L2
# test needed
[docs]
def angular_channels(self):
channels = []
for l in self.l_channels:
if l in self.components:
channels.append(l)
#end if
#end for
return channels
#end def angular_channels
# evaluate r*potential based on a potential component object
# component representation is specific to each derived class
[docs]
def evaluate_comp_rV(self,r,l,vcomp):
self.not_implemented()
#end def evaluate_comp_rV
[docs]
def find_r_rng(self,r,rmin):
if r is None and self.numeric and not self.interpolatable:
r = self.r
#end if
rng = None
if rmin>1e-12:
rng = r>rmin
r = r[rng]
#end if
return r,rng
#end def find_r_rng
# evaluate potential based on a potential component object
# local, nonlocal, and L2 are all represented by separate component objects
[docs]
def evaluate_comp(self,r,l,vcomp,rpow=0,rmin=0,rret=True):
v = self.evaluate_comp_rV(r,l,vcomp)
r,rng = self.find_r_rng(r,rmin)
if rng is not None:
v = v[rng]
#end if
if rpow!=1:
v = r**(rpow-1)*v
#end if
if rret:
return r,v
else:
return v
#end if
#end def evaluate_comp
# evaluate the local component potential only
[docs]
def evaluate_local(self,r=None,rpow=0,rmin=0,rret=False):
l = self.local
if not self.has_component(l):
self.error('cannot evaluate local potential\nlocal potential is not present')
#end if
vcomp = self.get_component(l)
ret = self.evaluate_comp(r,l,vcomp,rpow,rmin,rret)
return ret
#end def evaluate_local
# evaluate a nonlocal component potential
[docs]
def evaluate_nonlocal(self,r=None,l=None,rpow=0,rmin=0,rret=False):
if l==self.local:
self.error('called evaluate_nonlocal requesting local potential\nthe l index of the local potential is: {0}'.format(self.local))
elif l=='L2':
self.error('called evaluate_nonlocal requesting L2 potential')
elif l is None:
self.error('called evaluate_nonlocal without specifying the angular channel')
elif not self.has_component(l):
self.error('cannot evaluate non-local potential\nlocal potential is not present\nrequested potential: {0}'.format(l))
#end if
vcomp = self.get_component(l)
ret = self.evaluate_comp(r,l,vcomp,rpow,rmin,rret)
return ret
#end def evaluate_nonlocal
# test needed
# evaluate the L2 component potential
[docs]
def evaluate_L2(self,r=None,rpow=0,rmin=0,rret=False):
l = 'L2'
if not self.has_component(l):
self.error('cannot evaluate L2 potential\nL2 potential is not present')
#end if
vcomp = self.get_component(l)
ret = self.evaluate_comp(r,l,vcomp,rpow,rmin,rret)
return ret
#end def evaluate_L2
# evaluate semilocal potential components in isolation
[docs]
def evaluate_component(self,r=None,l=None,rpow=0,rmin=0,rret=False,optional=False):
vcomp = self.get_component(l)
if vcomp is not None:
return self.evaluate_comp(r,l,vcomp,rpow,rmin,rret)
elif optional:
z = 0*self.evaluate_local(r,rpow,rmin,rret=False)
if rret:
return z,z
else:
return z
#end if
else:
self.error('requested evaluation of non-existent component\ncomponent requested: {0}'.format(l))
#end if
#end def evaluate_component
# evaluate angular momentum channel of full potential
[docs]
def evaluate_channel(self,r=None,l=None,rpow=0,rmin=0,rret=False,with_local=True,with_L2=True):
if l not in self.l_channels:
self.error('evaluate_channel must be called with a valid angular momentum label\nvalid options are l=s,p,d,f,...\nyou provided: l={0}'.format(l))
#end if
eval_any = False
loc_present = self.has_component(self.local)
l_present = self.has_component(l)
L2_present = self.has_component('L2')
if (l==self.local or with_local) and loc_present:
vloc = self.evaluate_local(r,rpow,rmin)
eval_any = True
else:
vloc = 0
#end if
if with_L2 and L2_present:
l_int = self.channel_indices[l]
vL2 = self.evaluate_L2(r,rpow,rmin)*l_int*(l_int+1)
eval_any = True
else:
vL2 = 0
#end if
if l!=self.local and l_present:
vnonloc = self.evaluate_nonlocal(r,l,rpow,rmin)
eval_any = True
else:
vnonloc = 0
#end if
if rret or not eval_any:
r,rng = self.find_r_rng(r,rmin)
#end if
if eval_any:
v = vloc+vL2+vnonloc
else:
v = 0*r
#end if
if rret:
return r,v
else:
return v
#end if
#end def evaluate_channel
# similar to evaluate_channel but with defaults appropriate for QMCPACK
[docs]
def numeric_channel(self,l=None,rmin=0.,rmax=10.,npts=10001,rpow=0,with_local=True,with_L2=True):
if self.numeric and not self.interpolatable:
r = None
else:
r = np.linspace(rmin,rmax,npts)
#end if
if l=='L2':
r,v = self.evaluate_L2(r,rpow,rmin,rret=True)
else:
r,v = self.evaluate_channel(r,l,rpow,rmin,rret=True,with_local=with_local,with_L2=with_L2)
#end if
return r,v
#end def numeric_channel
[docs]
def update_rcut(self,tol=1e-5,optional=False):
if not optional or self.rcut is None:
self.rcut = self.find_rcut(tol=tol,with_L2=False)
#end if
has_vL2 = self.has_component('L2')
if has_vL2 and (not optional or self.rcut_L2 is None):
self.rcut_L2 = self.find_rcut(tol=tol,with_L2=True)
#end if
return self.rcut
#end def update_rcut
[docs]
def find_rcut(self,tol=1e-5,with_L2=False):
vnl = self.get_nonlocal()
if with_L2 and self.has_L2():
vnl.L2 = self.get_L2()
#end if
if len(vnl)==0:
rmax = 10.
return rmax
#end if
channels = list(vnl.keys())
rv = obj()
# add a zero potential
l = channels[0]
rvl = self.numeric_channel(l,rmin=0.01,with_local=False,with_L2=l=='L2')
rv[l] = rvl
rv['0'] = rvl[0],0*rvl[1]
for l in channels[1:]:
rv[l] = self.numeric_channel(l,rmin=0.01,with_local=False,with_L2=l=='L2')
#end for
r = None
vmin = None
vmax = None
for l,(rc,vc) in rv.items():
if r is None:
r = rc
vmin = np.array(vc)
vmax = np.array(vc)
elif len(rc)!=len(r):
self.error('numeric representation of channels do not match in length')
else:
vmin = np.minimum(vmin,vc)
vmax = np.maximum(vmax,vc)
#end if
#end for
vspread = vmax-vmin
rcut = r[-1]
nr = len(r)
for i in range(nr):
n = nr-1-i
if vspread[n]>tol:
rcut = r[n]
break
#end if
#end for
#figure()
#plot(r,vspread,'k-')
#plot([rcut,rcut],[0,vspread[1:].max()],'k--')
#title('rcut = {0}'.format(rcut))
#show()
return rcut
#end def find_rcut
[docs]
def plot(self,r=None,show=True,fig=True,linestyle='-',channels=None,with_local=False,rmin=0.01,rmax=5.0,title=None,metric=None,color=None):
if channels is None:
channels = self.l_channels
#end if
if fig:
plt.figure()
#end if
if r is None and self.numeric:
r = self.r
elif r is None:
r = np.linspace(rmin,rmax,1000)
#end if
rin = r
color_in = color
for c in channels:
r = rin
color = color_in
if c in self.components:
if color is None:
color = self.channel_colors[c]
#end if
lab = c
if c==self.local:
lab = self.local+' loc'
#end if
if self.name is not None:
lab = self.name+' '+lab
#end if
v = self.evaluate_channel(r,c,with_local=with_local,rmin=rmin-1e-12)
rng = r>rmin-1e-12
r = r[rng]
if metric=='r2':
v = r**2*v
if c==self.local:
v += self.Zval*r
#end if
elif metric is not None:
self.error('invalid metric for plotting: {0}\nvalid options are: r2'.format(metric))
#end if
plt.plot(r,v,color+linestyle,label=lab)
#end for
#end for
if fig:
if title is None:
title = 'Semilocal {0} PP ({1} core)'.format(self.element,self.core)
#end if
plt.title(title)
plt.ylabel('channel potentials (Ha)')
plt.xlabel('r (Bohr)')
plt.legend()
#end if
if show:
plt.show()
#end if
#end def plot
[docs]
def plot_components(self,r=None,show=True,fig=True,linestyle='-',rmin=0.01,rmax=5.0,title=None,metric=None,color=None,rpow=0):
channels = list(self.l_channels)+['L2']
if fig:
plt.figure(tight_layout=True)
#end if
if r is None and self.numeric:
r = self.r
#r = None
elif r is None:
r = np.linspace(rmin,rmax,1000)
#end if
rin = r
color_in = color
for c in channels:
r = rin
color = color_in
channel = self.get_component(c)
if channel is not None:
if color is None:
color = self.channel_colors[c]
#end if
lab = c
if c==self.local:
lab += ' loc'
elif c!='L2':
lab += '-'+self.local
#end if
if self.name is not None:
lab = self.name+' '+lab
#end if
v = self.evaluate_component(r,c,rpow,rmin-1e-12)
rng = r>rmin-1e-12
r = r[rng]
if metric=='r2':
v = r**2*v
if c==self.local:
v += self.Zval*r
#end if
elif metric is not None:
self.error('invalid metric for plotting: {0}\nvalid options are: r2'.format(metric))
#end if
plt.plot(r,v,color+linestyle,label=lab)
#end for
#end for
if fig:
if title is None:
title = 'Semilocal {0} PP ({1} core)'.format(self.element,self.core)
#end if
plt.title(title)
plt.ylabel('component potentials (Ha)')
plt.xlabel('r (Bohr)')
plt.legend()
#end if
if show:
plt.show()
#end if
#end def plot_components
[docs]
def plot_channels(self,r=None,channels=None,show=True,fig=True,linestyle='-',rmin=0.01,rmax=5.0,title=None,metric=None,color=None,rpow=0,with_local=True,with_L2=True):
if channels is None:
channels = list(self.l_channels)
#end if
if fig:
plt.figure(tight_layout=True)
#end if
if r is None and self.numeric:
r = self.r
#r = None
elif r is None:
r = np.linspace(rmin,rmax,1000)
#end if
rin = r
color_in = color
if not self.has_component('L2'):
loc_label = self.local
for c in channels:
if c not in self.components:
loc_label+=c
#end if
#end for
for c in channels:
r = rin
color = color_in
if self.has_component(c):
if color is None:
color = self.channel_colors[c]
#end if
lab = c
if c==self.local:
lab = loc_label
#end if
if self.name is not None:
lab = self.name+' '+lab
#end if
v = self.evaluate_channel(r,c,rpow,rmin-1e-12,False,with_local,with_L2)
rng = r>rmin-1e-12
r = r[rng]
if metric=='r2':
v = r**2*v
elif metric is not None:
self.error('invalid metric for plotting: {0}\nvalid options are: r2'.format(metric))
#end if
plt.plot(r,v,color+linestyle,label=lab)
#end for
#end for
else:
for c in channels:
r = rin
color = color_in
if color is None:
color = self.channel_colors[c]
#end if
lab = c
if self.name is not None:
lab = self.name+' '+lab
#end if
v = self.evaluate_channel(r,c,rpow,rmin-1e-12,False,with_local,with_L2)
rng = r>rmin-1e-12
r = r[rng]
if metric=='r2':
v = r**2*v
elif metric is not None:
self.error('invalid metric for plotting: {0}\nvalid options are: r2'.format(metric))
#end if
plt.plot(r,v,color+linestyle,label=lab)
#end for
#end if
if fig:
if title is None:
title = 'Semilocal {0} PP angular channels ({1} core)'.format(self.element,self.core)
#end if
plt.title(title)
plt.ylabel('channels')
plt.xlabel('r')
plt.legend()
#end if
if show:
plt.show()
#end if
#end def plot_channels
[docs]
def plot_positive_definite(self,r=None,show=True,fig=True,linestyle='-',rmin=0.01,rmax=5.0,title=None,color='k'):
if not self.has_L2():
self.error('positive definite condition only applies to L2 potentials')
#end if
if fig:
plt.figure(tight_layout=True)
#end if
if r is None and self.numeric:
r = self.r
elif r is None:
r = np.linspace(rmin,rmax,1000)
#end if
vL2 = self.evaluate_L2(r,0,rmin-1e-12)
rng = r>rmin-1e-12
r = r[rng]
b = vL2*(2*r**2)
plt.plot(r,1+b,color+linestyle,label='1+b')
plt.plot(r,0*r,'r-')
if fig:
if title is None:
title = 'L2 positive definite condition {0} PP ({1} core)'.format(self.element,self.core)
#end if
plt.title(title)
plt.ylabel('1+b > 0')
plt.xlabel('r (Bohr)')
plt.legend()
#end if
if show:
plt.show()
#end if
#end def plot_positive_definite
[docs]
def plot_L2(self,show=True,fig=True,r=None,rmin=0.01,rmax=5.0,linestyle='-',title=None,color=None):
if fig:
plt.figure(tight_layout=True)
#end if
if r is None and self.numeric:
r = self.r
elif r is None:
r = np.linspace(rmin,rmax,1000)
#end if
vs = self.evaluate_channel(r,'s',with_local=True,rmin=rmin-1e-12)
for c in self.l_channels[1:]:
if c in self.components:
if color is None:
color = self.channel_colors[c]
#end if
v = self.evaluate_channel(r,c,with_L2=False,rmin=rmin-1e-12)
rng = r>rmin-1e-12
r = r[rng]
l = self.channel_indices[c]
vL2 = (v-vs)/(l*(l+1))
plt.plot(r,vL2,color+linestyle,label='(v{0}-vs)/(l(l+1))'.format(c))
#end if
#end for
if fig:
plt.xlim([0,rmax])
if title is None:
title = 'Semilocal {0} PP ({1} core)'.format(self.element,self.core)
#end if
plt.title(title)
plt.ylabel('vL2 for channels above s')
plt.xlabel('r')
plt.legend()
if show:
plt.show()
#end if
#end if
#end def plot_L2
[docs]
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
plt.figure()
vmin = 1e99
vmax = -1e99
for l in self.l_channels:
if l in vl:
color = self.channel_colors[l]
v = vl[l]
plt.plot(r,v,color+'-',label='v'+l)
vmin = min(v.min(),vmin)
vmax = max(v.max(),vmax)
#end if
#end for
plt.plot([rc,rc],[vmin,vmax],'k--',lw=2)
plt.xlim([-0.1*rc,1.1*rc])
plt.xlabel('r (Bohr)')
plt.ylabel('V NL (Ha)')
plt.title((tlabel+' NL channels').strip())
plt.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 = np.linspace(vmin,vmax,levels)
fig = plt.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()
plt.xlim(lim)
plt.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)
plt.plot(rc*cos,rc*sin,'k--',lw=2)
fig.colorbar(cs, ax=ax, shrink=0.9)
plt.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:
plt.show()
#end if
#end def plot_nonlocal_polar
[docs]
def write_qmcpack(self,filepath=None):
self.update_rcut(tol=1e-5,optional=True)
channels = self.angular_channels()
symbol = self.element
atomic_number = self.Zcore+self.Zval
zval = self.Zval
creator = 'Nexus'
npots_down = len(channels)
l_local = self.channel_indices[self.local]
if l_local == -1:
self.error('Local channel, {}, not coded.'.format(self.local))
#end if
rmin = 1e99
rmax = -1e99
npts = 0
vps = obj()
for l in channels:
r,v = self.numeric_channel(l,rpow=1,with_local=True,with_L2=False)
rmin = min(rmin,r.min())
rmax = max(rmax,r.max())
npts = len(r)
vps[l] = v
#end for
header = '''<?xml version="1.0" encoding="UTF-8"?>
<pseudo version="0.5">
<header symbol="{0}" atomic-number="{1}" zval="{2}" relativistic="unknown"
polarized="unknown" creator="{3}" flavor="unknown"
core-corrections="unknown" xc-functional-type="unknown"
xc-functional-parametrization="unknown"/>
'''.format(symbol,atomic_number,zval,creator)
grid = ' <grid type="linear" units="bohr" ri="{0}" rf="{1}" npts="{2}"/>\n'.format(rmin,rmax,npts)
L2 = ''
if self.has_component('L2'):
dpad = '\n '
L2 += ' <L2 units="hartree" format="r*V" cutoff="{0}">\n'.format(self.rcut_L2)
L2 += ' <radfunc>\n'
L2 += ' '+grid
L2 += ' <data>'
r,v = self.numeric_channel('L2',rpow=1)
n=0
for d in v:
if n%3==0:
L2 += dpad
#end if
L2 += '{0:22.14e}'.format(d)
n+=1
#end for
L2 = L2.rstrip()+'\n'
L2 += ' </data>\n'
L2 += ' </radfunc>\n'
L2 += ' </L2>\n'
#end if
semilocal = ' <semilocal units="hartree" format="r*V" npots-down="{0}" npots-up="0" l-local="{1}">\n'.format(npots_down,l_local)
dpad = '\n '
for l in self.l_channels:
if l in vps:
semilocal+=' <vps principal-n="0" l="{0}" spin="-1" cutoff="{1}" occupation="unknown">\n'.format(l,self.rcut)
semilocal+=' <radfunc>\n'
semilocal+=' '+grid
semilocal+=' <data>'
v = vps[l]
n=0
for d in v:
if n%3==0:
semilocal+=dpad
#end if
semilocal+='{0:22.14e}'.format(d)
n+=1
#end for
semilocal = semilocal.rstrip()+'\n'
semilocal+=' </data>\n'
semilocal+=' </radfunc>\n'
semilocal+=' </vps>\n'
#end if
#end for
semilocal+=' </semilocal>\n'
footer = '</pseudo>\n'
text = header+grid+L2+semilocal+footer
if filepath is not None:
with open (filepath, "w") as f:
f.write(text)
#end if
return text
#end def write_qmcpack
[docs]
def write_casino(self,filepath=None):
if self.has_component('L2'):
self.error('cannot write potential in CASINO format\nan L2 term is present, but this is not supported by CASINO')
#end if
channels = self.angular_channels()
name = self.name
symbol = self.element
atomic_number = self.Zcore+self.Zval
zval = float(self.Zval)
l_local = 'spdfgi'.find(self.local)
if name is None:
name = '{0} pseudopotential converted by Nexus'.format(symbol)
#end if
rmin = 1e99
rmax = -1e99
npts = 0
vps = obj()
for l in channels:
r,v = self.numeric_channel(l,rpow=1,with_local=True,with_L2=False)
rmin = min(rmin,r.min())
rmax = max(rmax,r.max())
npts = len(r)
vps[l] = v
#end for
header = '''{0}
Atomic number and pseudo-charge
{1} {2}
Energy units (rydberg/hartree/ev):
hartree
Angular momentum of local component (0=s,1=p,2=d..)
{3}
NLRULE override (1) VMC/DMC (2) config gen (0 ==> input/default value)
0 0
Number of grid points
{4}
'''.format(name,atomic_number,zval,l_local,npts)
grid = 'R(i) in atomic units\n'
for d in r:
grid += ' {0:20.14e}\n'.format(d)
#end for
channels = ''
for l in self.l_channels:
if l in vps:
channels += 'r*potential (L={0}) in Ha\n'.format('spdfgi'.find(l))
v = vps[l]
for d in v:
channels += ' {0:20.14e}\n'.format(d)
#end for
#end if
#end for
text = header+grid+channels
if filepath is not None:
open(filepath,'w').write(text)
#end if
return text
#end def write_casino
#end class SemilocalPP
[docs]
class GaussianPP(SemilocalPP):
requires_format = True
formats = SemilocalPP.formats + 'gaussian gamess crystal numhf'.split()
[docs]
@staticmethod
def process_float(s):
return float(s.replace('D','e').replace('d','e'))
#end def process_float
def __init__(self,filepath=None,format=None,name=None,src=None):
self.basis = None
SemilocalPP.__init__(self,filepath,format,name,src)
#end def __init__
# test needed for gaussian and crystal
[docs]
def read_text(self,text,format=None,filepath=None):
lines,basis_lines = process_gaussian_text(text,format)
format=format.lower()
channels = []
basis = None
# need lmax, element, and Zcore
if format=='gamess':
i=0
name,type,Zcore,lmax = lines[i].split(); i+=1
Zcore = int(Zcore)
lmax = int(lmax)
element = pp_elem_label(name)[0]
while i<len(lines):
n = int(lines[i]); i+=1
terms = []
for j in range(n):
coeff,rpow,expon = lines[i].split(); i+=1
terms.append((float(coeff),int(rpow),float(expon)))
#end for
channels.append(terms)
#end while
elif format=='gaussian':
i=0
element,token = lines[i].split(); i+=1
label,lmax,Zcore = lines[i].split(); i+=1
lmax = int(lmax)
Zcore = int(Zcore)
while i<len(lines):
i+=1 # skip comment line
n = int(lines[i]); i+=1
terms = []
for j in range(n):
rpow,expon,coeff = lines[i].split(); i+=1
terms.append((float(coeff),int(rpow),float(expon)))
#end for
channels.append(terms)
#end while
elif format=='crystal':
i = 0
conv_atomic_number,nshells = lines[i].split(); i+=1
if len(conv_atomic_number)==1:
atomic_number = int(conv_atomic_number)
else:
atomic_number = int(conv_atomic_number[-2:])
#end if
element = Elements(atomic_number).symbol
if 'input' not in lines[i].lower():
self.error('INPUT must be present for crystal pseudpotential read')
#end if
i+=1
tokens = lines[i].split()
Zval = int(float(tokens[0])); i+=1
Zcore = atomic_number-Zval
nterms = np.array(tokens[1:],dtype=int)
lmax = 0
if nterms[0]==0:
lmax+=1
channels.append([(0.0,2,1.0)])
nterms = nterms[1:]
#end if
for nt in nterms:
lmax += 1
terms = []
for n in range(nt):
expon,coeff,rpow = lines[i].split(); i+=1
terms.append((float(coeff),int(rpow)+2,float(expon)))
#end for
channels.append(terms)
#end for
lmax-=1
elif format=='atomscf':
#i=0
#self.name = lines[i].strip(); i+=1
i=1 # skip title line
element = 'Rn' # text does not contain element (must be corrected downstream)
lmax = -1
Zcore = int(lines[i].strip()); i+=1
while i<len(lines):
n = int(lines[i]); i+=1
terms = []
for j in range(n):
rpow,expon,coeff = lines[i].split(); i+=1
terms.append((float(coeff),int(rpow),float(expon)))
#end for
channels.append(terms)
lmax+=1
#end while
# PPs in atomscf input are s,p,d,f, etc
# rather than, e.g. f,s-f,p-f,d-f
# so rearrange into form similar to f,s-f,p-f,d-f
loc = channels.pop()
for l in range(lmax):
c = channels[l]
channels[l] = c[0:len(c)-len(loc)]
#end for
channels = [loc] + channels
elif format=='numhf':
name = None
i=0
Zval,lmax = lines[i].split(); i+=1
Zval = int(Zval)
lmax = int(lmax)-1
element = self.element
Zcore = int(Elements(element).atomic_number)-Zval
ns = [int(n) for n in lines[i].split()]; i+=1
while i<len(lines):
for n in ns:
terms = []
for j in range(n):
rpow,expon,coeff = lines[i].split(); i+=1
terms.append((float(coeff),int(rpow),float(expon)))
#end for
channels.append(terms)
#end for
#end while
# Bring local channel to front
channels.insert(0,channels.pop())
else:
self.error('ability to read file format {0} has not been implemented'.format(format))
#end if
if basis_lines is not None:
bs = GaussianBasisSet()
bs.read_lines(basis_lines,format)
basis = bs.basis
#end if
if not Elements.is_element(element):
if not Elements.is_element(self.element):
self.error('cannot identify element for pseudopotential file '+path_string(filepath))
#end if
else:
self.element = element
#end if
Zatom = Elements(element).atomic_number
Zval = Zatom-Zcore
if Zcore==0:
core = None
else:
core = Elements(Zcore).symbol
#end if
self.set(
core = core,
Zval = Zval,
Zcore = Zcore,
lmax = lmax
)
for c in range(len(channels)):
if c==0:
#cname = 'loc'
cname = self.l_channels[lmax]
self.local = cname
else:
cname = self.l_channels[c-1]
#end if
channel = obj()
terms = channels[c]
for t in range(len(terms)):
coeff,rpow,expon = terms[t]
channel[t] = obj(coeff=coeff,rpow=rpow,expon=expon)
#end for
self.components[cname] = channel
#end for
self.basis = basis
if len(self.components)!=self.lmax+1:
self.error('number of channels is not lmax+1!')
#end if
#end def read_text
# test needed for crystal
[docs]
def write_text(self,format=None,occ=None):
text = ''
format = format.lower()
if format=='qmcpack':
return self.write_qmcpack()
elif format=='casino':
return self.write_casino()
#end if
channel_order = [self.local]
for c in self.l_channels:
if c in self.components and c!=self.local:
channel_order.append(c)
#end if
#end for
basis = self.basis
if basis is not None:
bs = GaussianBasisSet()
bs.basis = basis
basis = bs
#end if
if format=='gamess':
if basis is not None:
text += '{0} {1} 0. 0. 0.\n'.format(self.element,self.Zcore+self.Zval)
text += basis.write_text(format)
text += '\n'
#end if
text += '{0}-PP GEN {1} {2}\n'.format(self.element,self.Zcore,self.lmax)
for c in channel_order:
channel = self.components[c]
text += '{0}\n'.format(len(channel))
for i in sorted(channel.keys()):
g = channel[i]
text += '{0:12.8f} {1} {2:12.8f}\n'.format(g.coeff,g.rpow,g.expon)
#end for
#end for
text += '\n'
elif format=='gaussian':
if basis is not None:
text += '{0} 0\n'.format(self.element)
text += basis.write_text(format)
text += '\n'
#end if
text += '{0} 0\n'.format(self.element)
text += '{0}_PP {1} {2}\n'.format(self.element,self.lmax,self.Zcore)
for c in channel_order:
channel = self.components[c]
text += '{0} channel\n'.format(c)
text += '{0}\n'.format(len(channel))
for i in sorted(channel.keys()):
g = channel[i]
text += '{0} {1:12.8f} {2:12.8f}\n'.format(g.rpow,g.expon,g.coeff)
#end for
#end for
text += '\n'
elif format=='crystal':
if basis is not None:
conv_atomic_number = 200 + Elements(self.element).atomic_number
text+='{0} {1}\n'.format(conv_atomic_number,basis.size())
btext = basis.write_text(format,occ=occ)
else:
btext = ''
#end if
text += 'INPUT\n'
tline = '{0}'.format(int(self.Zval))
channels = []
cloc = self.components[channel_order[0]]
if len(cloc)==1 and abs(cloc[0].coeff)<1e-8:
tline += ' 0'
else:
tline += ' {0}'.format(len(cloc))
channels.append(cloc)
#end if
ccount = 1
for c in channel_order[1:]:
channel = self.components[c]
tline += ' {0}'.format(len(channel))
channels.append(channel)
ccount += 1
#end for
for i in range(6-ccount): # crystal14 goes up to g (hence the 6)
tline += ' 0'
#end for
text += tline+'\n'
for channel in channels:
for i in sorted(channel.keys()):
g = channel[i]
text += '{0} {1} {2}\n'.format(g.expon,g.coeff,g.rpow-2)
#end for
#end for
text += btext
elif format=='atomscf':
text += '{0} core potential\n'.format(self.element)
text += '{0}\n'.format(self.Zcore)
local_channel = self.components[self.local]
for c in self.l_channels:
if c in self.components:
channel = self.components[c]
if c!=self.local:
text += '{0}\n'.format(len(channel)+len(local_channel))
else:
text += '{0}\n'.format(len(channel))
#end if
for i in sorted(channel.keys()):
g = channel[i]
text += '{0} {1:12.8f} {2:12.8f}\n'.format(g.rpow,g.expon,g.coeff)
#end for
if c!=self.local:
channel = local_channel
for i in sorted(channel.keys()):
g = channel[i]
text += '{0} {1:12.8f} {2:12.8f}\n'.format(g.rpow,g.expon,g.coeff)
#end for
#end if
#end if
#end for
text += '\n'
elif format=='numhf':
channel_order = self.l_channels[:self.lmax+1]
text += '{} {}\n'.format(self.Zval,len(self.components))
for c in channel_order:
text += '{} '.format(len(self.components[c]))
#end for
text = text[:-1]+'\n'
for c in channel_order:
comp = self.components[c]
for i in sorted(comp.keys()):
g = comp[i]
text += '{0} {1:12.8f} {2:12.8f}\n'.format(g.rpow,g.expon,g.coeff)
#end for
#end for
else:
self.error('ability to write file format {0} has not been implemented'.format(format))
#end if
return text
#end def write_text
# test needed
[docs]
def get_basis(self):
bs = None
if self.basis is not None:
bs = GaussianBasisSet()
bs.basis = self.basis.copy()
#end if
return bs
#end def get_basis
[docs]
def set_basis(self,bs):
self.basis = bs.basis
#end def set_basis
# test needed
[docs]
def uncontract(self):
if self.basis is not None:
bs = GaussianBasisSet()
bs.basis = self.basis.copy()
bs.uncontract()
self.basis = bs.basis
#end if
#end def uncontract
# test needed
[docs]
def write_basis(self,filepath=None,format=None):
basis = self.get_basis()
text = ''
if basis is not None:
if format=='gamess':
text += '{0} {1} 0. 0. 0.\n'.format(self.element,self.Zcore+self.Zval)
text += basis.write_text(format)
text += '\n'
elif format=='gaussian':
text += '{0} 0\n'.format(self.element)
text += basis.write_text(format)
text += '\n'
else:
self.error('ability to write basis for file format {0} has not been implemented'.format(format))
#end if
#end if
if filepath is not None:
fobj = open(filepath,'w')
fobj.write(text)
fobj.close()
#end if
return text
#end def write_basis
[docs]
def evaluate_comp_rV(self,r,l,vcomp):
r = np.array(r)
v = np.zeros(r.shape)
if l==self.local or l is None:
v += -self.Zval
#end if
for g in vcomp:
if g.rpow==1:
v += g.coeff * np.exp(-g.expon*r**2)
else:
v += g.coeff * r**(g.rpow-1) * np.exp(-g.expon*r**2)
#end if
#end for
return v
#end def evaluate_comp_rV
# test needed
[docs]
def ppconvert(self,outfile,ref,extra=None):
of = outfile.lower()
if of.endswith('.xml'):
opts = '--xml'
elif of.endswith('.upf'):
opts = '--log_grid --upf'
else:
self.error('output file format unrecognized for {0}\nvalid extensions are .xml and .upf'.format(outfile))
#end if
tmpfile = 'tmp.gamess'
self.write(tmpfile,'gamess')
if extra is not None:
command = 'ppconvert --gamess_pot {0} --s_ref "{1}" --p_ref "{1}" --d_ref "{1}" {2} {3} {4}'.format(tmpfile,ref,extra,opts,outfile)
else:
command = 'ppconvert --gamess_pot {0} --s_ref "{1}" --p_ref "{1}" --d_ref "{1}" {2} {3}'.format(tmpfile,ref,opts,outfile)
execute(command,verbose=True)
os.system('rm '+tmpfile)
#end def ppconvert
# test needed
[docs]
def append_to_component(self,l,coeff,expon,rpow):
'''
This function is used to append a term to a Gaussian ECP component.
l: the angular ccomponent that the Gaussian term will be appended to
coeff, expon, rpow: the coefficient, exponent, and r-power of the Gaussian term
'''
if l>self.lmax:
self.error('component {} not present in PP.'.format(l))
#end if
chan_labels = ['s','p','d','f','g','h','i','j']
self.components[chan_labels[l]].append(obj(coeff=coeff,expon=expon,rpow=rpow))
#end def append_to_component
# test needed
[docs]
def scale_component(self,l,scale):
'''
This function is used to scale a Gaussian ECP component by a factor.
l: the angular component that is scaled.
scale: the scaling factor
'''
if l>self.lmax:
self.error('component {} not present in PP.'.format(l))
#end if
chan_labels = ['s','p','d','f','g','h','i','j']
for term in self.components[chan_labels[l]]:
term.coeff*=scale
#end for
#end def scale_component
# test needed
[docs]
def simplify(self):
'''This function simplifies the Gaussian ECP.
The simplificactions are as follows:
1. Remove all terms with coefficients that are equal to zero -- unless only one term exists.
2. Within each component, look for terms that have matching exponents and r-powers, if any are
present, then sum their coefficicents to make a single term. If the coefficients sum to
zero, then remove the terms (unless that leaves the component empty)
'''
dec=16
# Remove terms with coefficients equivalent to zero
#chan_labels = ['s','p','d','f','g','h','i','k']
chan_labels = list(self.l_channels)
remove = []
for l in np.arange(self.lmax+1):
for term_idx,term in enumerate(self.components[chan_labels[l]]):
if abs(term.coeff)<1e-12 and len(self.components[chan_labels[l]])>1:
remove.append((chan_labels[l],term_idx))
#end if
#end for
#end for
for r in remove:
self.components[r[0]].delete(r[1])
#end for
comps = self.components.copy()
for l in np.arange(self.lmax+1):
comps[chan_labels[l]] = obj()
for term_idx,term in enumerate(self.components[chan_labels[l]]):
comps[chan_labels[l]].append(term)
#end for
#end for
self.components = comps.copy()
comps = self.components.copy()
for l in np.arange(self.lmax+1):
terms = []
comps[chan_labels[l]] = obj()
for term_idx,term in enumerate(self.components[chan_labels[l]]):
terms.append(term.list())
#end for
terms = np.array(terms)
like_terms = []
if len(terms)>1:
rpows = terms[:,2]
expons = terms[:,1]
coeffs = terms[:,0]
for ex_idx,ex in enumerate(expons.round(decimals=dec)):
if any(ex_idx in subl for subl in like_terms):
continue
#end if
match = np.argwhere(expons.round(decimals=dec)==ex)
if len(match)>1:
match = match.flatten()
unique_pows = np.unique(rpows[match].round(decimals=dec))
if len(unique_pows)==1:
like_terms.append(match.tolist())
else:
for uv in unique_pows:
uv_count = rpows[match].round(decimals=dec).tolist().count(uv)
if uv_count>1:
m = match[rpows[match].round(decimals=dec)==uv][0]
if any(m in subl for subl in like_terms):
continue
else:
like_terms.append(match[rpows[match].round(decimals=dec)==uv].tolist())
#end if
#end if
#end for
#end if
#end if
#end for
#end if
# update comps
comps[chan_labels[l]] = obj()
added = []
for term_idx,term in enumerate(self.components[chan_labels[l]]):
if any(term_idx in subl for subl in like_terms):
if term_idx in added:
continue
else:
for mlist in like_terms:
if term_idx in mlist and term_idx not in added:
coeff = 0.0
mod_term = term.copy()
for ti in mlist:
coeff += self.components[chan_labels[l]][ti].coeff
#end for
if abs(coeff)>1e-12:
mod_term.coeff = coeff
comps[chan_labels[l]].append(mod_term)
#end if
added.extend(mlist)
#end if
#end for
#end if
else:
comps[chan_labels[l]].append(term)
added.append(term_idx)
#end if
#end for
if len(comps[chan_labels[l]])==0:
# All terms cancelled. Add placeholder
plcehldr = self.components['s'][0].copy()
plcehldr.coeff = 0.0
plcehldr.rpow = 2
plcehldr.expon = 1.0
comps[chan_labels[l]].append(plcehldr)
#end for
self.components = comps.copy()
#end def simplify
# test needed
[docs]
def is_truncated_L2(self):
'''
Determine if the Gaussian ECP's channels follow an L2 relationship.
'''
# CHECK IF THIS WORKS FOR lmax=1 !!!!!!
# Only checked for lmax=2 and higher
p1 = self.copy()
p1.simplify()
p2 = self.copy()
p2.transform_to_truncated_L2(keep='s p',lmax=p2.lmax)
p2.simplify()
return object_eq(p2,p1)
#end def is_truncated_L2
# test needed
[docs]
def get_unboundedness(self,db,dbs):
'''
This function quantifies how unbounded a truncated L2 potential is.
This is done by constructing a function that corrects VL2 in the unbounded region.
Then integrating the difference between VL2 and the correcting function.
'''
if not self.is_truncated_L2():
self.error('The PP must be in the truncated L2 form.')
#end if
import math
def poly(x,c):
val=0
for ci,cv in enumerate(c):
val+=cv*x**ci
return val
#end def
def Rs(x,dx,s,c):
if x+1-s<-dx:
return 0-(1-s)
elif x+1-s>dx:
return x
else:
return poly(x+1-s,c)-(1-s)
#end def
A=[]
for i in range(8):
row=[]
for j in range(8):
if i<4:
if j-i<0:
dcoeff=0
dpower=0
else:
dcoeff=math.factorial(j)/math.factorial(j-i)
dpower=j-i
#end if
row.append(dcoeff*db**dpower)
else:
if j-(i-4)<0:
dcoeff=0
dpower=0
else:
dcoeff=math.factorial(j)/math.factorial(j-(i-4))
dpower=j-(i-4)
#end if
row.append(dcoeff*(-db)**dpower)
#end if
#end for
A.append(row)
#end for
A = np.array(A)
b = np.array([db,1]+[0]*6)
c = np.linalg.inv(A).dot(b)
ng=3000
gmin=0.02
gmax=0.85
r = np.linspace( gmin, gmax, ng )
# PP
v = []
for l in ['s','p']:
vtmp = []
for ri in r:
vtmp.append(self.evaluate_component(r=ri,l=l))
#for
v.append(vtmp)
#for
v=np.array(v)
# 2*r^2*VL2
if self.lmax>1:
f = r*r*(v[1]-v[0])
elif self.lmax==1:
f = -r*r*v[0]
else:
self.error('Not sure what to do with fully local potential.')
#end if
# 2*r^2*V'L2
fp = [Rs(fr,db,dbs,c) for fr in f]
unboundedness = 0
for fi,fx in enumerate(f):
unboundedness+=(fp[fi]-fx)*(gmax-gmin)/ng
#end for
return unboundedness
#end def get_unboundedness
# test needed
[docs]
def make_L2_bounded(self,db,dbs,exps0=None,plot=False):
'''
For a truncated L2 potential, this function constructs a correction to VL2 in the unbounded region.
Then the correction is fit to a set of Gaussian primitives that are provided in the array 'exps0'.
The fitted Gaussian primitives are then appended to the ECP, resulting in a bounded truncated L2 potential.
'''
if not self.is_truncated_L2():
self.error('The PP must be in the truncated L2 form.')
#end if
if exps0 is None:
self.error('Please provide a aet of exponents to be used for correction.')
import math
def poly(x,c):
val=0
for ci,cv in enumerate(c):
val+=cv*x**ci
return val
#end def
def Rs(x,dx,s,c):
if x+1-s<-dx:
return 0-(1-s)
elif x+1-s>dx:
return x
else:
return poly(x+1-s,c)-(1-s)
#end def
class fitClass:
def __init__(self):
pass
def gauss_correction(self,x,c1,c2,c3):
val = 0
for ci,c in enumerate([c1,c2,c3]):
val+=x**2.*c*np.exp(-self.exps[ci]*x**2.)
#end for
return val
#end def
def gauss_correction_2_param(self,x,c1,c2):
val = 0
for ci,c in enumerate([c1,c2]):
val+=x**2.*c*np.exp(-self.exps[ci]*x**2.)
#end for
return val
#end def
def gauss_correction_1_param(self,x,c1):
val = 0
for ci,c in enumerate([c1]):
val+=x**2.*c*np.exp(-self.exps[ci]*x**2.)
#end for
return val
#end def
#end class
A=[]
for i in range(8):
row=[]
for j in range(8):
if i<4:
if j-i<0:
dcoeff=0
dpower=0
else:
dcoeff=math.factorial(j)/math.factorial(j-i)
dpower=j-i
#end if
row.append(dcoeff*db**dpower)
else:
if j-(i-4)<0:
dcoeff=0
dpower=0
else:
dcoeff=math.factorial(j)/math.factorial(j-(i-4))
dpower=j-(i-4)
#end if
row.append(dcoeff*(-db)**dpower)
#end if
#end for
A.append(row)
#end for
A = np.array(A)
b = np.array([db,1]+[0]*6)
c = np.linalg.inv(A).dot(b)
ng=3000
gmin=0.02
gmax=0.85
r = np.linspace( gmin, gmax, ng )
# PP
v = []
for l in ['s','p']:
vtmp = []
for ri in r:
vtmp.append(self.evaluate_component(r=ri,l=l))
#for
v.append(vtmp)
#for
v=np.array(v)
# 2*r^2*VL2
if self.lmax>1:
f = r*r*(v[1]-v[0])
elif self.lmax==1:
f = -r*r*v[0]
else:
self.error('Not sure what to do with fully local potential.')
#end if
# 2*r^2*V'L2
fp = [Rs(fr,db,dbs,c) for fr in f]
unboundedness = 0
for fi,fx in enumerate(f):
unboundedness+=(fp[fi]-fx)*(gmax-gmin)/ng
#end for
#print('\npseudopotential undoundedness: ',undoundedness)
from scipy.optimize import curve_fit
fit_instance = fitClass()
fit_instance.exps=exps0
if len(exps0)==1:
popt, pcov = curve_fit(fit_instance.gauss_correction_1_param, r, f-fp)
elif len(exps0)==2:
popt, pcov = curve_fit(fit_instance.gauss_correction_2_param, r, f-fp)
elif len(exps0)==3:
popt, pcov = curve_fit(fit_instance.gauss_correction, r, f-fp)
else:
self.error('Number of correction primitives not coded.')
#end if
if plot:
import matplotlib.pyplot as plt
plt.plot(r, fp, 'g-', label='fp')
plt.xlabel('r (bohr)')
plt.ylabel('$2r^2v_{L^2}$')
if len(exps0)==1:
plt.plot(r, f-fit_instance.gauss_correction_1_param(r, *popt), 'r-',label='f-corr')
elif len(exps0)==2:
plt.plot(r, f-fit_instance.gauss_correction_2_param(r, *popt), 'r-',label='f-corr')
elif len(exps0)==3:
plt.plot(r, f-fit_instance.gauss_correction(r, *popt), 'r-',label='f-corr')
#end if
plt.plot(r, f, 'b-',label='f')
plt.plot(r, [-1]*len(r), 'k-',label=None)
plt.legend()
plt.show()
#end if
for expon_idx,expon in enumerate(exps0):
self.components['s'].append(obj(coeff=1.0*popt[expon_idx],expon=expon,rpow=2))
#end if
self.transform_to_truncated_L2(keep='s p',lmax=self.lmax)
self.simplify()
#end def make_L2_bounded
# test needed
#end def transform_to_truncated_L2
#end class GaussianPP
[docs]
class QmcpackPP(SemilocalPP):
requires_format = False
numeric = True
interpolatable = False
[docs]
def read(self,filepath,format=None):
if not os.path.exists(filepath):
self.error('cannot read {0}, file does not exist'.format(filepath))
#end if
x = readxml(filepath,contract_names=True)
x.convert_numeric()
x.condense()
x.remove_hidden()
pp = x.pseudo
h = pp.header
self.element = h.symbol
self.Zval = h.zval
self.Zcore = h.atomic_number-h.zval
if self.Zcore==0:
self.core = '0'
else:
self.core = Elements(self.Zcore).symbol
#end if
g = pp.grid
if g.type=='linear':
self.rmin = g.ri
self.rmax = g.rf
self.r = np.linspace(g.ri,g.rf,g.npts)
else:
self.error('functionality for '+g.type+' grids has not yet been implemented')
#end if
if 'l2' in pp:
l2 = pp.l2
if l2.format!='r*V':
self.error('unrecognized potential format: {0}\nthe only supported format is r*V'.format(l2.format))
#end if
self.components.L2 = l2.radfunc.data.copy()
#end if
sl = pp.semilocal
if sl.format!='r*V':
self.error('unrecognized potential format: {0}\nthe only supported format is r*V'.format(sl.format))
#end if
lloc = self.l_channels[sl.l_local]
self.local = lloc
vps = sl.vps
if not isinstance(vps,list):
vps = [vps]
#end if
self.lmax = len(vps)-1
for vp in vps:
self.components[vp.l] = vp.radfunc.data.copy()
#end for
for l in self.angular_channels():
if l!=self.local:
self.components[l] -= self.components[self.local]
#end if
#end for
#end def read
[docs]
def evaluate_comp_rV(self,r,l,vcomp):
if r is not None:
if len(r)==len(self.r) and abs( (r[1:]-self.r[1:])/self.r[1:] ).max()<1e-6:
r = self.r
else:
self.error('ability to interpolate at arbitrary r has not been implemented\ncalling evaluate_channel() without specifying r will return the potential on a default grid')
#end if
else:
r = self.r
#end if
v = vcomp.copy()
return v
#end def evaluate_comp_rV
[docs]
def v_at_zero(self,l):
#r = self.r
#v = self.get_component(l)/r
#vz = (v[1]*r[2]**2-v[2]*r[1]**2)/(r[2]**2-r[1]**2)
r = self.r[1:3]
v = self.get_component(l)[1:3]/r
vz = (v[0]*r[1]**2-v[1]*r[0]**2)/(r[1]**2-r[0]**2)
return vz
#end def v_at_zero
#end class QmcpackPP
[docs]
class CasinoPP(SemilocalPP):
requires_format = False
numeric = True
interpolatable = False
unitmap = dict(rydberg='Ry',hartree='Ha',ev='eV')
[docs]
def read(self,filepath,format=None):
filepath = path_string(filepath)
if not os.path.exists(filepath):
self.error('cannot read {0}, file does not exist'.format(filepath))
#end if
# open the file
file = TextFile(filepath)
# read scalar values at the top
Zatom,Z = file.readtokensf('Atomic number and pseudo-charge',int,float)
if Zatom > Elements.num_elements():
self.error('element {0} is not in the periodic table')
#end if
element = Elements(Zatom).symbol
units = file.readtokensf('Energy units',str)
if units not in self.unitmap:
self.error('units {0} unrecognized from casino PP file {1}'.format(units,filepath))
#end if
lloc = file.readtokensf('Angular momentum of local component',int)
lloc = self.l_channels[lloc]
ngrid = file.readtokensf('Number of grid points',int)
# read the radial grid
file.seek('R(i)',1)
file.readline()
r = np.empty((ngrid,),dtype=float)
for ir in range(ngrid):
r[ir] = float(file.readline())
#end for
# read each channel, convert to hartree units
lmax = -1
lvals = []
lpots = []
while(file.seek('pot',1)!=-1):
lmax += 1
potline = file.readline() # read the r*potential line
eqloc = potline.find('=')
if eqloc==-1:
self.error('"=" not found in potential line\nline: {0}'.format(potline))
#end if
l = self.l_channels[int(potline[eqloc+1])] # get the l value
lvals.append(l)
v = np.empty((ngrid,),dtype=float)
for ir in range(ngrid):
v[ir] = float(file.readline())
#end for
lpots.append(convert(v,self.unitmap[units],'Ha'))
#end while
file.close()
# fill in SemilocalPP class data obtained from read
self.element = element
self.Zval = int(Z)
self.Zcore = Zatom-self.Zval
if self.Zcore==0:
self.core = '0'
else:
self.core = Elements(self.Zcore).symbol
#end if
self.lmax = lmax
self.local = lloc
self.r = r
for i in range(len(lpots)):
self.components[lvals[i]] = lpots[i]
#end for
for l in self.angular_channels():
if l!=self.local:
self.components[l] -= self.components[self.local]
#end if
#end for
#end def read_file
[docs]
def evaluate_comp_rV(self,r,l,vcomp):
if r is not None:
if len(r)==len(self.r) and abs( (r[1:]-self.r[1:])/self.r[1:] ).max()<1e-6:
r = self.r
else:
self.error('ability to interpolate at arbitrary r has not been implemented\ncalling evaluate_channel() without specifying r will return the potential on a default grid')
#end if
else:
r = self.r
#end if
v = vcomp.copy()
return v
#end def evaluate_comp_rV
#end class CasinoPP
