mirror of
https://github.com/triqs/dft_tools
synced 2024-07-11 05:43:48 +02:00
Removing examples from source code directories.
This commit is contained in:
Alexander Hampel
parent
18abc77e33
commit
c591fe351e
|
|
@ -1,9 +0,0 @@
|
|||
*
|
||||
!.gitignore
|
||||
!INCAR
|
||||
!POSCAR
|
||||
!POTCAR
|
||||
!KPOINTS
|
||||
!plo.cfg
|
||||
!*.sh
|
||||
!*.py
|
||||
|
|
@ -1,49 +0,0 @@
|
|||
ADDGRID = TRUE
|
||||
AMIX = 0.1
|
||||
AMIN = 0.01
|
||||
BMIX = 0.0001
|
||||
#IMIX = 1
|
||||
#MIXFIRST = True
|
||||
ALGO = NORMAL
|
||||
EDIFF = 0.00001
|
||||
EDIFFG = -0.002
|
||||
EMAX = 13.0
|
||||
EMIN = -2.0
|
||||
ENCUT = 600
|
||||
ENAUG = 1100
|
||||
IBRION = -1
|
||||
ICHARG = 5
|
||||
ISIF = 4
|
||||
ISMEAR = 2
|
||||
ISPIN = 2
|
||||
ISYM = -1
|
||||
LASPH = TRUE
|
||||
LCHARG = TRUE
|
||||
LDAU = True
|
||||
LDAUJ = 0
|
||||
LDAUL = 3
|
||||
LDAUTYPE = 2
|
||||
LDAUU = 0.0
|
||||
LDAUPRINT = 2
|
||||
LMAXMIX = 6
|
||||
LORBIT = 13
|
||||
LREAL = False
|
||||
ROPT = 1e-3 1e-3 1e-3
|
||||
LWAVE = TRUE
|
||||
MAGMOM = 1.1
|
||||
#NBANDS = 16
|
||||
NEDOS = 501
|
||||
NELM = 100
|
||||
NELMDL = -7
|
||||
NELMIN = 1
|
||||
NPAR = 4
|
||||
NSIM = 1
|
||||
NSW = 0
|
||||
POTIM = 0.4
|
||||
PREC = Accurate
|
||||
SIGMA = 0.1
|
||||
SMASS = 0.5
|
||||
SYMPREC = 1.0e-6
|
||||
|
||||
# Projectors
|
||||
LOCPROJ = 1 : f : Pr 1
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
Automatic kpoint scheme
|
||||
0
|
||||
Gamma
|
||||
4 4 4
|
||||
|
|
@ -1,9 +0,0 @@
|
|||
gamma-Ce
|
||||
1.0
|
||||
2.5790000000000000 2.5790000000000000 0.0000000000000000
|
||||
0.0000000000000000 2.5790000000000000 2.5790000000000000
|
||||
2.5790000000000000 0.0000000000000000 2.5790000000000000
|
||||
Ce
|
||||
1
|
||||
Direct
|
||||
0.0 0.0 0.0
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,252 +0,0 @@
|
|||
from math import sqrt
|
||||
from scipy.misc import factorial as fact
|
||||
from itertools import product
|
||||
import numpy as np
|
||||
|
||||
# The interaction matrix in desired basis
|
||||
# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
|
||||
def U_matrix(l, radial_integrals=None, U_int=None, J_hund=None, basis="spherical", T=None):
|
||||
"""Calculate U matrix being given either radial_integrals or U_int and J_hund.
|
||||
l = angular momentum of shell being treated (l=2 for d shell, l=3 for f shell)
|
||||
radial_integrals = [F0,F2,F4,..] (default None)
|
||||
U_int, J_hund = values to use to compute radial_integrals (default None),
|
||||
basis = "spherical", "cubic", or "other",
|
||||
T = transformation matrix from spherical to desired basis, if basis='other' (default None)"""
|
||||
|
||||
# Check all necessary information is present and consistent
|
||||
if radial_integrals is None and (U_int is None and J_hund is None):
|
||||
raise ValueError("U_matrix: provide either the radial_integrals or U_int and J_hund.")
|
||||
if radial_integrals is None and (U_int is not None and J_hund is not None):
|
||||
radial_integrals = U_J_to_radial_integrals(l, U_int, J_hund)
|
||||
if radial_integrals is not None and (U_int is not None and J_hund is not None):
|
||||
if len(radial_integrals)-1 != l:
|
||||
raise ValueError("U_matrix: inconsistency in l and number of radial_integrals provided.")
|
||||
if (radial_integrals - U_J_to_radial_integrals(l, U_int, J_hund)).any() != 0.0:
|
||||
print "Warning: U_matrix: radial_integrals provided do not match U_int and J_hund. Using radial_integrals to calculate spherical U_matrix."
|
||||
|
||||
# Full interaction matrix
|
||||
# Basis of spherical harmonics Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}
|
||||
# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
|
||||
U_matrix = np.zeros((2*l+1,2*l+1,2*l+1,2*l+1),dtype=float)
|
||||
|
||||
m_range = range(-l,l+1)
|
||||
for n, F in enumerate(radial_integrals):
|
||||
k = 2*n
|
||||
for m1, m2, m3, m4 in product(m_range,m_range,m_range,m_range):
|
||||
U_matrix[m1+l,m2+l,m3+l,m4+l] += F * angular_matrix_element(l,k,m1,m2,m3,m4)
|
||||
|
||||
# Transform from spherical basis if needed
|
||||
if basis == "cubic": T = spherical_to_cubic(l)
|
||||
if basis == "other" and T is None:
|
||||
raise ValueError("U_matrix: provide T for other bases.")
|
||||
if T is not None: U_matrix = transform_U_matrix(U_matrix, T)
|
||||
|
||||
return U_matrix
|
||||
|
||||
# Convert full 4-index U matrix to 2-index density-density form
|
||||
def reduce_4index_to_2index(U_4index):
|
||||
"""Reduces the four-index matrix to two-index matrices."""
|
||||
|
||||
size = len(U_4index) # 2l+1
|
||||
U = np.zeros((size,size),dtype=float) # matrix for same spin
|
||||
Uprime = np.zeros((size,size),dtype=float) # matrix for opposite spin
|
||||
|
||||
m_range = range(size)
|
||||
for m,mp in product(m_range,m_range):
|
||||
U[m,mp] = U_4index[m,mp,m,mp].real - U_4index[m,mp,mp,m].real
|
||||
Uprime[m,mp] = U_4index[m,mp,m,mp].real
|
||||
|
||||
return U, Uprime
|
||||
|
||||
# Construct the 2-index matrices for the density-density form
|
||||
def U_matrix_kanamori(n_orb, U_int, J_hund):
|
||||
"""Calculate the Kanamori U and Uprime matrices."""
|
||||
|
||||
U = np.zeros((n_orb,n_orb),dtype=float) # matrix for same spin
|
||||
Uprime = np.zeros((n_orb,n_orb),dtype=float) # matrix for opposite spin
|
||||
|
||||
m_range = range(n_orb)
|
||||
for m,mp in product(m_range,m_range):
|
||||
if m == mp:
|
||||
Uprime[m,mp] = U_int
|
||||
else:
|
||||
U[m,mp] = U_int - 3.0*J_hund
|
||||
Uprime[m,mp] = U_int - 2.0*J_hund
|
||||
|
||||
return U, Uprime
|
||||
|
||||
# Get t2g or eg components
|
||||
def t2g_submatrix(U, convention=''):
|
||||
"""Return only the t2g part of the full d-manifold two- or four-index U matrix."""
|
||||
if convention == 'wien2k':
|
||||
return subarray(U, len(U.shape)*[(2,3,4)])
|
||||
else:
|
||||
return subarray(U, len(U.shape)*[(0,1,3)])
|
||||
|
||||
def eg_submatrix(U, convention=''):
|
||||
"""Return only the eg part of the full d-manifold two- or four-index U matrix."""
|
||||
if convention == 'wien2k':
|
||||
return subarray(U, len(U.shape)*[(0,1)])
|
||||
else:
|
||||
return subarray(U, len(U.shape)*[(2,4)])
|
||||
|
||||
# Transform the interaction matrix into another basis
|
||||
def transform_U_matrix(U_matrix, T):
|
||||
"""Transform the interaction matrix into another basis by applying matrix T."""
|
||||
return np.einsum("ij,kl,jlmo,mn,op",np.conj(T),np.conj(T),U_matrix,np.transpose(T),np.transpose(T))
|
||||
|
||||
# Rotation matrices: complex harmonics to cubic harmonics
|
||||
# Complex harmonics basis: ..., Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}, ...
|
||||
def spherical_to_cubic(l, convention=''):
|
||||
"""Returns the spherical harmonics to cubic harmonics rotation matrix."""
|
||||
size = 2*l+1
|
||||
T = np.zeros((size,size),dtype=complex)
|
||||
if convention != 'wien2k' and l == 2:
|
||||
raise ValueError("spherical_to_cubic: standard convention not implemented for l=1")
|
||||
if l == 0:
|
||||
cubic_names = ("s")
|
||||
elif l == 1:
|
||||
cubic_names = ("x","y","z")
|
||||
T[0,0] = 1.0/sqrt(2); T[0,2] = -1.0/sqrt(2)
|
||||
T[1,0] = 1j/sqrt(2); T[1,2] = 1j/sqrt(2)
|
||||
T[2,1] = 1.0
|
||||
elif l == 2:
|
||||
if convention == 'wien2k':
|
||||
cubic_names = ("z^2","x^2-y^2","xy","yz","xz")
|
||||
T[0,2] = 1.0
|
||||
T[1,0] = 1.0/sqrt(2); T[1,4] = 1.0/sqrt(2)
|
||||
T[2,0] = -1j/sqrt(2); T[2,4] = 1j/sqrt(2)
|
||||
T[3,1] = 1j/sqrt(2); T[3,3] = -1j/sqrt(2)
|
||||
T[4,1] = 1.0/sqrt(2); T[4,3] = 1.0/sqrt(2)
|
||||
else:
|
||||
cubic_names = ("xy","yz","z^2","xz","x^2-y^2")
|
||||
T[0,0] = 1j/sqrt(2); T[0,4] = -1j/sqrt(2)
|
||||
T[1,1] = 1j/sqrt(2); T[1,3] = 1j/sqrt(2)
|
||||
T[2,2] = 1.0
|
||||
T[3,1] = 1.0/sqrt(2); T[3,3] = -1.0/sqrt(2)
|
||||
T[4,0] = 1.0/sqrt(2); T[4,4] = 1.0/sqrt(2)
|
||||
elif l == 3:
|
||||
if convention == 'wien2k':
|
||||
cubic_names = ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
|
||||
T[0,0] = 1.0/sqrt(2); T[0,6] = -1.0/sqrt(2)
|
||||
T[1,1] = 1.0/sqrt(2); T[1,5] = 1.0/sqrt(2)
|
||||
T[2,2] = 1.0/sqrt(2); T[2,4] = -1.0/sqrt(2)
|
||||
T[3,5] = 1.0
|
||||
T[4,2] = 1j/sqrt(2); T[4,4] = 1j/sqrt(2)
|
||||
T[5,1] = 1j/sqrt(2); T[5,5] = -1j/sqrt(2)
|
||||
T[6,0] = 1j/sqrt(2); T[6,6] = 1j/sqrt(2)
|
||||
else:
|
||||
cubic_names = ("y(3x^2-y^2)","xyz","yz^2","z^3","xz^2","z(x^2-y^2)","x(x^2-3y^2)")
|
||||
# Verified with the matrix generated in VASP
|
||||
T[0,0] = 1j/sqrt(2); T[0,6] = 1j/sqrt(2)
|
||||
T[1,1] = 1j/sqrt(2); T[1,5] = -1j/sqrt(2)
|
||||
T[2,2] = 1j/sqrt(2); T[2,4] = 1j/sqrt(2)
|
||||
T[3,3] = 1.0
|
||||
T[4,2] = 1.0/sqrt(2); T[4,4] = -1.0/sqrt(2)
|
||||
T[5,1] = 1.0/sqrt(2); T[5,5] = 1.0/sqrt(2)
|
||||
T[6,0] = 1.0/sqrt(2); T[6,6] = -1.0/sqrt(2)
|
||||
else: raise ValueError("spherical_to_cubic: implemented only for l=0,1,2,3")
|
||||
|
||||
return T
|
||||
|
||||
# Names of cubic harmonics
|
||||
def cubic_names(l):
|
||||
if l == 0 or l == 's':
|
||||
return ("s")
|
||||
elif l == 1 or l == 'p':
|
||||
return ("x","y","z")
|
||||
elif l == 2 or l == 'd':
|
||||
return ("xy","yz","z^2","xz","x^2-y^2")
|
||||
elif l == 't2g':
|
||||
return ("xy","yz","xz")
|
||||
elif l == 'eg':
|
||||
return ("z^2","x^2-y^2")
|
||||
elif l == 3 or l == 'f':
|
||||
return ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
|
||||
else: raise ValueError("cubic_names: implemented only for l=0,1,2,3")
|
||||
|
||||
# Convert U,J -> radial integrals F_k
|
||||
def U_J_to_radial_integrals(l, U_int, J_hund):
|
||||
"""Determines the radial integrals F_k from U_int and J_hund."""
|
||||
|
||||
F = np.zeros((l+1),dtype=float)
|
||||
if l == 2:
|
||||
F[0] = U_int
|
||||
F[1] = J_hund * 14.0 / (1.0 + 0.63)
|
||||
F[2] = 0.630 * F[1]
|
||||
elif l == 3:
|
||||
F[0] = U_int
|
||||
F[1] = 6435.0 * J_hund / (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0)
|
||||
F[2] = 451.0 * F[1] / 675.0
|
||||
F[3] = 1001.0 * F[1] / 2025.0
|
||||
else: raise ValueError("U_J_to_radial_integrals: implemented only for l=2,3")
|
||||
|
||||
return F
|
||||
|
||||
# Convert radial integrals F_k -> U,J
|
||||
def radial_integrals_to_U_J(l, F):
|
||||
"""Determines U_int and J_hund from the radial integrals."""
|
||||
|
||||
if l == 2:
|
||||
U_int = F[0]
|
||||
J_hund = F[1] * (1.0 + 0.63) / 14.0
|
||||
elif l == 3:
|
||||
U_int = F[0]
|
||||
J_hund = F[1] * (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0) / 6435.0
|
||||
else: raise ValueError("radial_integrals_to_U_J: implemented only for l=2,3")
|
||||
|
||||
return U_int,J_hund
|
||||
|
||||
# Angular matrix elements of particle-particle interaction
|
||||
# (2l+1)^2 ((l 0) (k 0) (l 0))^2 \sum_{q=-k}^{k} (-1)^{m1+m2+q} ((l -m1) (k q) (l m3)) ((l -m2) (k -q) (l m4))
|
||||
def angular_matrix_element(l, k, m1, m2, m3, m4):
|
||||
result = 0
|
||||
for q in range(-k,k+1):
|
||||
result += three_j_symbol((l,-m1),(k,q),(l,m3))*three_j_symbol((l,-m2),(k,-q),(l,m4))*(-1.0 if (m1+q+m2) % 2 else 1.0)
|
||||
result *= (2*l+1)**2 * (three_j_symbol((l,0),(k,0),(l,0))**2)
|
||||
return result
|
||||
|
||||
# Wigner 3-j symbols
|
||||
# ((j1 m1) (j2 m2) (j3 m3))
|
||||
def three_j_symbol(jm1, jm2, jm3):
|
||||
j1, m1 = jm1
|
||||
j2, m2 = jm2
|
||||
j3, m3 = jm3
|
||||
|
||||
if (m1+m2+m3 != 0 or
|
||||
m1 < -j1 or m1 > j1 or
|
||||
m2 < -j2 or m2 > j2 or
|
||||
m3 < -j3 or m3 > j3 or
|
||||
j3 > j1 + j2 or
|
||||
j3 < abs(j1-j2)):
|
||||
return .0
|
||||
|
||||
result = -1.0 if (j1-j2-m3) % 2 else 1.0
|
||||
result *= sqrt(fact(j1+j2-j3)*fact(j1-j2+j3)*fact(-j1+j2+j3)/fact(j1+j2+j3+1))
|
||||
result *= sqrt(fact(j1-m1)*fact(j1+m1)*fact(j2-m2)*fact(j2+m2)*fact(j3-m3)*fact(j3+m3))
|
||||
|
||||
t_min = max(j2-j3-m1,j1-j3+m2,0)
|
||||
t_max = min(j1-m1,j2+m2,j1+j2-j3)
|
||||
|
||||
t_sum = 0
|
||||
for t in range(t_min,t_max+1):
|
||||
t_sum += (-1.0 if t % 2 else 1.0)/(fact(t)*fact(j3-j2+m1+t)*fact(j3-j1-m2+t)*fact(j1+j2-j3-t)*fact(j1-m1-t)*fact(j2+m2-t))
|
||||
|
||||
result *= t_sum
|
||||
return result
|
||||
|
||||
# Clebsch-Gordan coefficients
|
||||
# < j1 m1 j2 m2 | j3 m3 > = (-1)^{j1-j2+m3} \sqrt{2j3+1} ((j1 m1) (j2 m2) (j3 -m3))
|
||||
def clebsch_gordan(jm1, jm2, jm3):
|
||||
norm = sqrt(2*jm3[0]+1)*(-1 if jm1[0]-jm2[0]+jm3[1] % 2 else 1)
|
||||
return norm*three_j_symbol(jm1,jm2,(jm3[0],-jm3[1]))
|
||||
|
||||
# Create subarray containing columns in idxlist
|
||||
# e.g. idxlist = [(0),(2,3),(0,1,2,3)] gives
|
||||
# column 0 for 1st dim,
|
||||
# columns 2 and 3 for 2nd dim,
|
||||
# columns 0,1,2 and 3 for 3rd dim.
|
||||
def subarray(a,idxlist,n=None) :
|
||||
if n is None: n = len(a.shape)-1
|
||||
sa = a[tuple(slice(x) for x in a.shape[:n]) + (idxlist[n],)]
|
||||
return subarray(sa,idxlist, n-1) if n > 0 else sa
|
||||
|
|
@ -1,335 +0,0 @@
|
|||
|
||||
################################################################################
|
||||
#
|
||||
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||
#
|
||||
# Copyright (C) 2011 by M. Ferrero, O. Parcollet
|
||||
#
|
||||
# TRIQS is free software: you can redistribute it and/or modify it under the
|
||||
# terms of the GNU General Public License as published by the Free Software
|
||||
# Foundation, either version 3 of the License, or (at your option) any later
|
||||
# version.
|
||||
#
|
||||
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
||||
# details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License along with
|
||||
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
||||
#
|
||||
################################################################################
|
||||
|
||||
from types import *
|
||||
#from triqs_dft_tools.U_matrix import *
|
||||
from U_matrix import *
|
||||
from pytriqs.gf import *
|
||||
#from hubbard_I import gf_hi_fullu, sigma_atomic_fullu
|
||||
import pytriqs.utility.mpi as mpi
|
||||
from itertools import izip
|
||||
import numpy
|
||||
import copy
|
||||
|
||||
class Solver:
|
||||
"""
|
||||
Hartree-Fock Solver
|
||||
"""
|
||||
|
||||
# initialisation:
|
||||
def __init__(self, beta, l, n_msb=1025, use_spin_orbit=False, Nmoments=5, dudarev=False):
|
||||
|
||||
self.name = "Hartree-Fock"
|
||||
self.beta = beta
|
||||
self.l = l
|
||||
self.Nmsb = n_msb
|
||||
self.UseSpinOrbit = use_spin_orbit
|
||||
self.Converged = False
|
||||
self.Nspin = 2
|
||||
self.Nmoments=Nmoments
|
||||
|
||||
self.dudarev = dudarev
|
||||
|
||||
assert use_spin_orbit == False, "Spin-orbit is not implemented"
|
||||
|
||||
self.Nlm = 2*l+1
|
||||
if (use_spin_orbit):
|
||||
# no blocks!
|
||||
self.gf_struct = [ ('ud', range(2*self.Nlm)) ]
|
||||
else:
|
||||
# up/down blocks:
|
||||
self.gf_struct = [ ('up', range(self.Nlm)), ('down', range(self.Nlm)) ]
|
||||
|
||||
# construct Greens functions:
|
||||
self.a_list = [a for a,al in self.gf_struct]
|
||||
glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
|
||||
self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
|
||||
self.G_Old = self.G.copy()
|
||||
self.G0 = self.G.copy()
|
||||
self.Sigma = self.G.copy()
|
||||
self.Sigma_Old = self.G.copy()
|
||||
|
||||
|
||||
|
||||
def solve(self, U_int, J_hund, T=None, verbosity=0, Iteration_Number=1, Test_Convergence=0.0001):
|
||||
"""Calculation of the impurity Greens function using Hubbard-I"""
|
||||
|
||||
if self.Converged :
|
||||
mpi.report("Solver %(name)s has already converged: SKIPPING"%self.__dict__)
|
||||
return
|
||||
|
||||
if mpi.is_master_node():
|
||||
self.verbosity = verbosity
|
||||
else:
|
||||
self.verbosity = 0
|
||||
|
||||
#self.Nmoments = 5
|
||||
|
||||
ur,ujmn,umn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
|
||||
|
||||
|
||||
M = [x for x in self.G.mesh]
|
||||
self.zmsb = numpy.array([x for x in M],numpy.complex_)
|
||||
|
||||
# # for the tails:
|
||||
# tailtempl={}
|
||||
# for sig,g in self.G:
|
||||
# tailtempl[sig] = copy.deepcopy(g.tail)
|
||||
# for i in range(9): tailtempl[sig][i] *= 0.0
|
||||
|
||||
# self.__save_eal('eal.dat',Iteration_Number)
|
||||
|
||||
# mpi.report( "Starting Fortran solver %(name)s"%self.__dict__)
|
||||
|
||||
self.Sigma_Old <<= self.Sigma
|
||||
self.G_Old <<= self.G
|
||||
|
||||
# # call the fortran solver:
|
||||
# temp = 1.0/self.beta
|
||||
# gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
|
||||
# zmsb=self.zmsb, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = self.verbosity)
|
||||
|
||||
#self.sig = sigma_atomic_fullu(gf=self.gf,e0f=self.eal,zmsb=self.zmsb,ns=self.Nspin,nlm=self.Nlm)
|
||||
def print_matrix(m):
|
||||
for row in m:
|
||||
print ''.join(map("{0:12.7f}".format, row))
|
||||
|
||||
|
||||
# Hartree-Fock solver
|
||||
self.Sigma.zero()
|
||||
dm = self.G.density()
|
||||
if mpi.is_master_node():
|
||||
# print
|
||||
# print " Reduced U-matrix:"
|
||||
# print " U:"
|
||||
# print_matrix(ujmn)
|
||||
# print " Up:"
|
||||
# print_matrix(umn)
|
||||
#
|
||||
## sig_test = {bl: numpy.zeros((self.Nlm, self.Nlm)) for bl in dm}
|
||||
# sig_test = {}
|
||||
# sig_test['up'] = numpy.dot(umn, dm['up'].real) + numpy.dot(ujmn, dm['down'].real)
|
||||
# sig_test['down'] = numpy.dot(umn, dm['down'].real) + numpy.dot(ujmn, dm['up'].real)
|
||||
# print " Sigma test:"
|
||||
# print_matrix(sig_test['up'])
|
||||
print
|
||||
print " Density matrix (up):"
|
||||
print_matrix(dm['up'])
|
||||
print
|
||||
print " Density matrix (down):"
|
||||
print_matrix(dm['down'])
|
||||
|
||||
if self.dudarev:
|
||||
Ueff = U_int - J_hund
|
||||
corr_energy = 0.0
|
||||
dft_dc = 0.0
|
||||
for bl1 in dm:
|
||||
# (U - J) * (1/2 - n)
|
||||
self.Sigma[bl1] << Ueff * (0.5 * numpy.identity(self.Nlm) - dm[bl1])
|
||||
# 1/2 (U - J) * \sum_{\sig} [\sum_{m} n_{m,m \sig} - \sum_{m1,m2} n_{m1,m2 \sig} n_{m2,m1 \sig}]
|
||||
corr_energy += 0.5 * Ueff * (dm[bl1].trace() - (dm[bl1]*dm[bl1].conj()).sum()).real
|
||||
# V[n] * n^{\dagger}
|
||||
dft_dc += (self.Sigma[bl1](0) * dm[bl1].conj()).sum().real
|
||||
|
||||
else:
|
||||
# !!!!!
|
||||
# !!!!! Mind the order of indices in the 4-index matrix!
|
||||
# !!!!!
|
||||
for il1 in xrange(self.Nlm):
|
||||
for il2 in xrange(self.Nlm):
|
||||
for il3 in xrange(self.Nlm):
|
||||
for il4 in xrange(self.Nlm):
|
||||
for bl1 in dm:
|
||||
for bl2 in dm:
|
||||
self.Sigma[bl1][il1, il2] += ur[il1, il3, il2, il4] * dm[bl2][il3, il4]
|
||||
if bl1 == bl2:
|
||||
self.Sigma[bl1][il1, il2] -= ur[il1, il3, il4, il2] * dm[bl1][il3, il4]
|
||||
|
||||
if mpi.is_master_node() and self.verbosity > 0:
|
||||
print
|
||||
print " Sigma (up):"
|
||||
print_matrix(self.Sigma['up'](0).real)
|
||||
print
|
||||
print " Sigma (down):"
|
||||
print_matrix(self.Sigma['down'](0).real)
|
||||
# if (self.verbosity==0):
|
||||
# # No fortran output, so give basic results here
|
||||
# mpi.report("Atomic occupancy in Hubbard I Solver : %s"%self.atocc)
|
||||
# mpi.report("Atomic magn. mom. in Hubbard I Solver : %s"%self.atmag)
|
||||
|
||||
# transfer the data to the GF class:
|
||||
if (self.UseSpinOrbit):
|
||||
nlmtot = self.Nlm*2 # only one block in this case!
|
||||
else:
|
||||
nlmtot = self.Nlm
|
||||
|
||||
# M={}
|
||||
# isp=-1
|
||||
# for a,al in self.gf_struct:
|
||||
# isp+=1
|
||||
# M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
|
||||
# for i in range(min(self.Nmoments,8)):
|
||||
# tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
|
||||
#
|
||||
# #glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb, data =M[a], tail =self.tailtempl[a])
|
||||
# # for a,al in self.gf_struct]
|
||||
# glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
|
||||
# self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
|
||||
#
|
||||
# self.__copy_Gf(self.G,M,tailtempl)
|
||||
#
|
||||
# # Self energy:
|
||||
# self.G0 <<= iOmega_n
|
||||
#
|
||||
# M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
|
||||
# self.G0 -= M
|
||||
# self.Sigma <<= self.G0 - inverse(self.G)
|
||||
#
|
||||
# # invert G0
|
||||
# self.G0.invert()
|
||||
|
||||
def test_distance(G1,G2, dist) :
|
||||
def f(G1,G2) :
|
||||
#print abs(G1.data - G2.data)
|
||||
dS = max(abs(G1.data - G2.data).flatten())
|
||||
aS = max(abs(G1.data).flatten())
|
||||
if mpi.is_master_node():
|
||||
print " Distances:", dS, " vs ", aS * dist
|
||||
return dS <= aS*dist
|
||||
return reduce(lambda x,y : x and y, [f(g1,g2) for (i1,g1),(i2,g2) in izip(G1,G2)])
|
||||
|
||||
mpi.report("\nChecking Sigma for convergence...\nUsing tolerance %s"%Test_Convergence)
|
||||
self.Converged = test_distance(self.Sigma,self.Sigma_Old,Test_Convergence)
|
||||
|
||||
if self.Converged :
|
||||
mpi.report("Solver HAS CONVERGED")
|
||||
else :
|
||||
mpi.report("Solver has not yet converged")
|
||||
|
||||
return corr_energy, dft_dc
|
||||
|
||||
|
||||
def GF_realomega(self, ommin, ommax, N_om, U_int, J_hund, T=None, verbosity=0, broadening=0.01):
|
||||
"""Calculates the GF and spectral function on the real axis."""
|
||||
|
||||
delta_om = (ommax-ommin)/(1.0*(N_om-1))
|
||||
|
||||
omega = numpy.zeros([N_om],numpy.complex_)
|
||||
|
||||
ur,umn,ujmn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
|
||||
|
||||
for i in range(N_om):
|
||||
omega[i] = ommin + delta_om * i + 1j * broadening
|
||||
|
||||
tailtempl={}
|
||||
for sig,g in self.G:
|
||||
tailtempl[sig] = copy.deepcopy(g.tail)
|
||||
for i in range(9): tailtempl[sig][i] *= 0.0
|
||||
|
||||
temp = 1.0/self.beta
|
||||
gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
|
||||
zmsb=omega, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = verbosity)
|
||||
|
||||
# transfer the data to the GF class:
|
||||
if (self.UseSpinOrbit):
|
||||
nlmtot = self.Nlm*2 # only one block in this case!
|
||||
else:
|
||||
nlmtot = self.Nlm
|
||||
|
||||
M={}
|
||||
isp=-1
|
||||
for a,al in self.gf_struct:
|
||||
isp+=1
|
||||
M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
|
||||
for i in range(min(self.Nmoments,8)):
|
||||
tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
|
||||
|
||||
#glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om, data = M[a], tail = self.tailtempl[a])
|
||||
# for a,al in self.gf_struct] # Indices for the upfolded G
|
||||
glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om) for a,al in self.gf_struct]
|
||||
self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
|
||||
|
||||
self.__copy_Gf(self.G,M,tailtempl)
|
||||
|
||||
# Self energy:
|
||||
self.G0 = self.G.copy()
|
||||
self.Sigma = self.G.copy()
|
||||
self.G0 <<= Omega + 1j*broadening
|
||||
|
||||
M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
|
||||
self.G0 -= M
|
||||
self.Sigma <<= self.G0 - inverse(self.G)
|
||||
self.Sigma.note='ReFreq' # This is important for the put_Sigma routine!!!
|
||||
|
||||
#sigmamat = sigma_atomic_fullu(gf=gf,e0f=self.ealmat,zmsb=omega,nlm=self.Nlm,ns=self.Nspin)
|
||||
|
||||
#return omega,gf,sigmamat
|
||||
|
||||
|
||||
def __save_eal(self,Filename,it):
|
||||
f=open(Filename,'a')
|
||||
f.write('\neff. atomic levels, Iteration %s\n'%it)
|
||||
for i in range(self.Nlm*self.Nspin):
|
||||
for j in range(self.Nlm*self.Nspin):
|
||||
f.write("%10.6f %10.6f "%(self.ealmat[i,j].real,self.ealmat[i,j].imag))
|
||||
f.write("\n")
|
||||
f.close()
|
||||
|
||||
def __copy_Gf(self,G,data,tail):
|
||||
""" Copies data and tail to Gf object GF """
|
||||
for s,g in G:
|
||||
g.data[:,:,:]=data[s][:,:,:]
|
||||
for imom in range(1,min(self.Nmoments,8)):
|
||||
g.tail.data[1+imom,:,:]=tail[s][imom]
|
||||
|
||||
|
||||
def set_atomic_levels(self,eal):
|
||||
""" Helps to set correctly the variables for the atomic levels from a dictionary."""
|
||||
|
||||
assert (type(eal)==DictType), "Give a dictionary to set_atomic_levels!"
|
||||
|
||||
cnt = 0
|
||||
self.ealmat[:,:] *= 0.0
|
||||
|
||||
for ind in eal:
|
||||
self.Eff_Atomic_Levels[ind] = copy.deepcopy(eal[ind])
|
||||
|
||||
if self.UseSpinOrbit:
|
||||
for ii in range(self.Nlm*2):
|
||||
for jj in range(self.Nlm*2):
|
||||
self.ealmat[ii,jj] = self.Eff_Atomic_Levels[ind][ii,jj]
|
||||
else:
|
||||
for ii in range(self.Nlm):
|
||||
for jj in range(self.Nlm):
|
||||
self.ealmat[cnt*self.Nlm + ii,cnt*self.Nlm + jj] = self.Eff_Atomic_Levels[ind][ii,jj]
|
||||
|
||||
cnt += 1
|
||||
|
||||
|
||||
def __set_umatrix(self,U,J,T=None):
|
||||
# U matrix:
|
||||
# l = (Nlm-1)/2
|
||||
# If T is specified, it is used to transform the Basis set
|
||||
Umat = U_matrix(l=self.l, U_int=U, J_hund=J, basis='cubic', T=T)
|
||||
U, Up = reduce_4index_to_2index(Umat)
|
||||
|
||||
return Umat, U, Up
|
||||
|
|
@ -1,15 +0,0 @@
|
|||
[General]
|
||||
#EFERMI =
|
||||
#BASENAME = plo_vasp
|
||||
#DOSMESH = -5.0 10.0 4001
|
||||
|
||||
[Group 1]
|
||||
SHELLS = 1
|
||||
NORMALIZE = False
|
||||
EWINDOW = -41.45 41.8
|
||||
|
||||
[Shell 1]
|
||||
LSHELL = 3
|
||||
IONS = 1
|
||||
|
||||
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
#!/bin/bash
|
||||
PYTHONPATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/python:$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/c:$PYTHONPATH pytriqs $@
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
PLOVASP_PATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src
|
||||
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@
|
||||
|
|
@ -1,164 +0,0 @@
|
|||
|
||||
import os
|
||||
import errno
|
||||
import signal
|
||||
import sys
|
||||
import time
|
||||
import pytriqs.utility.mpi as mpi
|
||||
import converters.plovasp.main as plovasp
|
||||
from test_ham_hf import dmft_cycle
|
||||
|
||||
debug = True
|
||||
#
|
||||
# Helper functions
|
||||
#
|
||||
def sigint_handler(signal, frame):
|
||||
raise SystemExit(1)
|
||||
|
||||
def is_vasp_lock_present():
|
||||
return os.path.isfile('./vasp.lock')
|
||||
|
||||
def is_vasp_running(vasp_pid):
|
||||
"""
|
||||
Tests if VASP initial process is still alive.
|
||||
"""
|
||||
pid_exists = False
|
||||
if mpi.is_master_node():
|
||||
try:
|
||||
os.kill(vasp_pid, 0)
|
||||
except OSError, e:
|
||||
pid_exists = e.errno == errno.EPERM
|
||||
else:
|
||||
pid_exists = True
|
||||
|
||||
pid_exists = mpi.bcast(pid_exists)
|
||||
return pid_exists
|
||||
|
||||
def get_dft_energy():
|
||||
"""
|
||||
Reads energy from the last line of OSZICAR.
|
||||
"""
|
||||
with open('OSZICAR', 'r') as f:
|
||||
nextline = f.readline()
|
||||
while nextline.strip():
|
||||
line = nextline
|
||||
nextline = f.readline()
|
||||
# print "OSZICAR: ", line[:-1]
|
||||
|
||||
try:
|
||||
dft_energy = float(line.split()[2])
|
||||
except ValueError:
|
||||
print "Cannot read energy from OSZICAR, setting it to zero"
|
||||
dft_energy = 0.0
|
||||
|
||||
return dft_energy
|
||||
|
||||
class bcolors:
|
||||
MAGENTA = '\033[95m'
|
||||
BLUE = '\033[94m'
|
||||
GREEN = '\033[92m'
|
||||
YELLOW = '\033[93m'
|
||||
RED = '\033[91m'
|
||||
ENDC = '\033[0m'
|
||||
#
|
||||
# Main self-consistent cycle
|
||||
#
|
||||
def run_all(vasp_pid):
|
||||
"""
|
||||
"""
|
||||
mpi.report(" Waiting for VASP lock to appear...")
|
||||
while not is_vasp_lock_present():
|
||||
time.sleep(1)
|
||||
|
||||
vasp_running = True
|
||||
|
||||
while vasp_running:
|
||||
if debug: print bcolors.RED + "rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
mpi.report(" Waiting for VASP lock to disappear...")
|
||||
mpi.barrier()
|
||||
while is_vasp_lock_present():
|
||||
time.sleep(1)
|
||||
# if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
if not is_vasp_running(vasp_pid):
|
||||
mpi.report(" VASP stopped")
|
||||
vasp_running = False
|
||||
break
|
||||
|
||||
if debug: print bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
err = 0
|
||||
exc = None
|
||||
try:
|
||||
if debug: print bcolors.BLUE + "plovasp: rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
if mpi.is_master_node():
|
||||
plovasp.generate_and_output_as_text('plo.cfg', vasp_dir='./')
|
||||
# Read energy from OSZICAR
|
||||
dft_energy = get_dft_energy()
|
||||
except Exception, exc:
|
||||
err = 1
|
||||
|
||||
err = mpi.bcast(err)
|
||||
if err:
|
||||
if mpi.is_master_node():
|
||||
raise exc
|
||||
else:
|
||||
raise SystemExit(1)
|
||||
|
||||
mpi.barrier()
|
||||
|
||||
try:
|
||||
if debug: print bcolors.GREEN + "rank %s"%(mpi.rank) + bcolors.ENDC
|
||||
corr_energy, dft_dc = dmft_cycle()
|
||||
except:
|
||||
if mpi.is_master_node():
|
||||
print " master forwarding the exception..."
|
||||
raise
|
||||
else:
|
||||
print " rank %i exiting..."%(mpi.rank)
|
||||
raise SystemExit(1)
|
||||
mpi.barrier()
|
||||
|
||||
if mpi.is_master_node():
|
||||
total_energy = dft_energy + corr_energy - dft_dc
|
||||
print
|
||||
print "="*80
|
||||
print " Total energy: ", total_energy
|
||||
print " DFT energy: ", dft_energy
|
||||
print " Corr. energy: ", corr_energy
|
||||
print " DFT DC: ", dft_dc
|
||||
print "="*80
|
||||
print
|
||||
|
||||
if mpi.is_master_node() and vasp_running:
|
||||
open('./vasp.lock', 'a').close()
|
||||
|
||||
if mpi.is_master_node():
|
||||
total_energy = dft_energy + corr_energy - dft_dc
|
||||
with open('TOTENERGY', 'w') as f:
|
||||
f.write(" Total energy: %s\n"%(total_energy))
|
||||
f.write(" DFT energy: %s\n"%(dft_energy))
|
||||
f.write(" Corr. energy: %s\n"%(corr_energy))
|
||||
f.write(" DFT DC: %s\n"%(dft_dc))
|
||||
f.write(" Energy correction: %s\n"%(corr_energy - dft_dc))
|
||||
|
||||
mpi.report("***Done")
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
try:
|
||||
vasp_pid = int(sys.argv[1])
|
||||
except (ValueError, KeyError):
|
||||
if mpi.is_master_node():
|
||||
print "VASP process pid must be provided as the first argument"
|
||||
raise
|
||||
# if len(sys.argv) > 1:
|
||||
# vasp_path = sys.argv[1]
|
||||
# else:
|
||||
# try:
|
||||
# vasp_path = os.environ['VASP_DIR']
|
||||
# except KeyError:
|
||||
# if mpi.is_master_node():
|
||||
# print "Path to VASP must be specified either as an argument of in VASP_DIR"
|
||||
# raise
|
||||
signal.signal(signal.SIGINT, sigint_handler)
|
||||
|
||||
run_all(vasp_pid)
|
||||
|
|
@ -1,9 +0,0 @@
|
|||
#/bin/bash
|
||||
|
||||
NPROC=8
|
||||
|
||||
rm -f vasp.lock
|
||||
stdbuf -o 0 mpirun -np $NPROC ~/Codes/vasp/build_20151202/vasp.5.4.2/build/std/vasp &
|
||||
|
||||
mpirun -np $NPROC ./run_build.sh sc_dmft.py $(jobs -p) || kill %1
|
||||
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
#!/bin/bash -l
|
||||
#
|
||||
#SBATCH --job-name="JOBNAME"
|
||||
#SBATCH --time=36:00:00
|
||||
##SBATCH --nodes=8
|
||||
#SBATCH --ntasks=8
|
||||
#SBATCH --partition=qmac
|
||||
##SBATCH --mem-per-cpu=1024
|
||||
#SBATCH --output=job.%j.o
|
||||
#SBATCH --error=job.%j.e
|
||||
|
||||
#======START=====
|
||||
# source $HOME/.profile
|
||||
# export KMP_AFFINITY=none
|
||||
|
||||
# module load slurm
|
||||
source /cm/shared/apps/intel/composer_xe/composer_xe_2015.0.090/mkl/bin/mklvars.sh intel64
|
||||
# module switch PrgEnv-gnu PrgEnv-intel
|
||||
# module load intel/14.0.2.144
|
||||
#module load intel/mkl/64/11.2/2015.0.090
|
||||
|
||||
export OMP_NUM_THREADS=1
|
||||
|
||||
VASP_DIR=$HOME/Codes/vasp/build_20151202/vasp.5.4.2/build/std
|
||||
#mpirun -np $SLURM_NTASKS $VASP_DIR/vasp
|
||||
./sc_dmft.sh
|
||||
|
||||
|
||||
|
|
@ -1,437 +0,0 @@
|
|||
#from triqs_dft_tools.sumk_dft import *
|
||||
from sumk_dft import *
|
||||
#from triqs_dft_tools.converters.wien2k_converter import *
|
||||
from converters.vasp_converter import *
|
||||
#from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
|
||||
from hf_solver import Solver
|
||||
import shutil
|
||||
|
||||
class TestSumkDFT(SumkDFT):
|
||||
# calculate and save occupancy matrix in the Bloch basis for VASP charge denity recalculation
|
||||
def calc_density_correction(self, filename='GAMMA', dm_type='wien2k'):
|
||||
r"""
|
||||
Calculates the charge density correction and stores it into a file.
|
||||
|
||||
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
|
||||
It represents a density matrix of the interacting system defined in Bloch basis
|
||||
and it is calculated from the sum over Matsubara frequecies of the full GF,
|
||||
|
||||
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
|
||||
|
||||
The density matrix for every `k`-point is stored into a file.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : string
|
||||
Name of the file to store the charge density correction.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(deltaN, dens) : tuple
|
||||
Returns a tuple containing the density matrix `deltaN` and
|
||||
the corresponing total charge `dens`.
|
||||
|
||||
"""
|
||||
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
|
||||
assert dm_type in ('vasp', 'wien2k'), "'type' must be either 'vasp' or 'wienk'"
|
||||
|
||||
ntoi = self.spin_names_to_ind[self.SO]
|
||||
spn = self.spin_block_names[self.SO]
|
||||
dens = {sp: 0.0 for sp in spn}
|
||||
|
||||
# Fetch Fermi weights and energy window band indices
|
||||
if dm_type == 'vasp':
|
||||
fermi_weights = 0
|
||||
band_window = 0
|
||||
if mpi.is_master_node():
|
||||
with HDFArchive(self.hdf_file,'r') as ar:
|
||||
fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
|
||||
band_window = ar['dft_misc_input']['band_window']
|
||||
fermi_weights = mpi.bcast(fermi_weights)
|
||||
band_window = mpi.bcast(band_window)
|
||||
|
||||
# Convert Fermi weights to a density matrix
|
||||
dens_mat_dft = {}
|
||||
for sp in spn:
|
||||
dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in xrange(self.n_k)]
|
||||
|
||||
# Set up deltaN:
|
||||
deltaN = {}
|
||||
for sp in spn:
|
||||
deltaN[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
|
||||
|
||||
ikarray = numpy.array(range(self.n_k))
|
||||
for ik in mpi.slice_array(ikarray):
|
||||
G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
|
||||
for bname,gf in G_latt_iw:
|
||||
deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
|
||||
dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
|
||||
if dm_type == 'vasp':
|
||||
# In 'vasp'-mode subtract the DFT density matrix
|
||||
diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
|
||||
deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
|
||||
dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
|
||||
|
||||
# mpi reduce:
|
||||
for bname in deltaN:
|
||||
for ik in range(self.n_k):
|
||||
deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y)
|
||||
dens[bname] = mpi.all_reduce(mpi.world, dens[bname], lambda x,y : x+y)
|
||||
mpi.barrier()
|
||||
|
||||
# now save to file:
|
||||
if dm_type == 'wien2k':
|
||||
if mpi.is_master_node():
|
||||
if self.SP == 0:
|
||||
f = open(filename,'w')
|
||||
else:
|
||||
f = open(filename+'up','w')
|
||||
f1 = open(filename+'dn','w')
|
||||
# write chemical potential (in Rydberg):
|
||||
f.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
|
||||
if self.SP != 0: f1.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
|
||||
# write beta in rydberg-1
|
||||
f.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
|
||||
if self.SP != 0: f1.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
|
||||
|
||||
if self.SP == 0: # no spin-polarization
|
||||
|
||||
for ik in range(self.n_k):
|
||||
f.write("%s\n"%self.n_orbitals[ik,0])
|
||||
for inu in range(self.n_orbitals[ik,0]):
|
||||
for imu in range(self.n_orbitals[ik,0]):
|
||||
valre = (deltaN['up'][ik][inu,imu].real + deltaN['down'][ik][inu,imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][inu,imu].imag + deltaN['down'][ik][inu,imu].imag) / 2.0
|
||||
f.write("%.14f %.14f "%(valre,valim))
|
||||
f.write("\n")
|
||||
f.write("\n")
|
||||
f.close()
|
||||
|
||||
elif self.SP == 1: # with spin-polarization
|
||||
|
||||
# dict of filename: (spin index, block_name)
|
||||
if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')}
|
||||
if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')}
|
||||
for fout in to_write.iterkeys():
|
||||
isp, sp = to_write[fout]
|
||||
for ik in range(self.n_k):
|
||||
fout.write("%s\n"%self.n_orbitals[ik,isp])
|
||||
for inu in range(self.n_orbitals[ik,isp]):
|
||||
for imu in range(self.n_orbitals[ik,isp]):
|
||||
fout.write("%.14f %.14f "%(deltaN[sp][ik][inu,imu].real,deltaN[sp][ik][inu,imu].imag))
|
||||
fout.write("\n")
|
||||
fout.write("\n")
|
||||
fout.close()
|
||||
elif dm_type == 'vasp':
|
||||
# assert self.SP == 0, "Spin-polarized density matrix is not implemented"
|
||||
|
||||
if mpi.is_master_node():
|
||||
with open(filename, 'w') as f:
|
||||
f.write(" %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
|
||||
for sp in spn:
|
||||
for ik in xrange(self.n_k):
|
||||
ib1 = band_window[0][ik, 0]
|
||||
ib2 = band_window[0][ik, 1]
|
||||
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
|
||||
for inu in xrange(self.n_orbitals[ik, 0]):
|
||||
for imu in xrange(self.n_orbitals[ik, 0]):
|
||||
if self.SP == 0:
|
||||
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
|
||||
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
|
||||
else:
|
||||
valre = deltaN[sp][ik][inu, imu].real
|
||||
valim = deltaN[sp][ik][inu, imu].imag
|
||||
|
||||
f.write(" %.14f %.14f"%(valre, valim))
|
||||
f.write("\n")
|
||||
else:
|
||||
raise NotImplementedError("Unknown density matrix type: '%s'"%(dm_type))
|
||||
|
||||
return deltaN, dens
|
||||
|
||||
def calc_hamiltonian_correction(self, filename='GAMMA'):
|
||||
r"""
|
||||
Calculates the charge density correction and stores it into a file.
|
||||
|
||||
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
|
||||
It represents a density matrix of the interacting system defined in Bloch basis
|
||||
and it is calculated from the sum over Matsubara frequecies of the full GF,
|
||||
|
||||
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
|
||||
|
||||
The density matrix for every `k`-point is stored into a file.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : string
|
||||
Name of the file to store the charge density correction.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(deltaN, dens) : tuple
|
||||
Returns a tuple containing the density matrix `deltaN` and
|
||||
the corresponing total charge `dens`.
|
||||
|
||||
"""
|
||||
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
|
||||
|
||||
ntoi = self.spin_names_to_ind[self.SO]
|
||||
spn = self.spin_block_names[self.SO]
|
||||
dens = {sp: 0.0 for sp in spn}
|
||||
|
||||
# Fetch Fermi weights and energy window band indices
|
||||
fermi_weights = 0
|
||||
band_window = 0
|
||||
if mpi.is_master_node():
|
||||
with HDFArchive(self.hdf_file,'r') as ar:
|
||||
fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
|
||||
band_window = ar['dft_misc_input']['band_window']
|
||||
fermi_weights = mpi.bcast(fermi_weights)
|
||||
band_window = mpi.bcast(band_window)
|
||||
|
||||
# Set up deltaH:
|
||||
deltaH = {}
|
||||
for sp in spn:
|
||||
deltaH[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
|
||||
|
||||
ikarray = numpy.array(range(self.n_k))
|
||||
for ik in mpi.slice_array(ikarray):
|
||||
sigma_minus_dc = [s.copy() for s in self.Sigma_imp_iw]
|
||||
sigma_minus_dc = self.add_dc('iw')
|
||||
|
||||
beta = self.Sigma_imp_iw[0].mesh.beta # override beta if Sigma_iw is present
|
||||
n_iw = len(self.Sigma_imp_iw[0].mesh)
|
||||
|
||||
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
|
||||
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
|
||||
block_ind_list = [block for block,inner in gf_struct]
|
||||
glist = lambda : [ GfImFreq(indices=inner,beta=beta,n_points=n_iw) for block,inner in gf_struct]
|
||||
G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
|
||||
G_latt.zero()
|
||||
|
||||
for icrsh in range(self.n_corr_shells):
|
||||
for bname, gf in G_latt:
|
||||
gf += self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf)
|
||||
|
||||
for sp in spn:
|
||||
deltaH[sp][ik][:, :] = G_latt[sp](0) # Any Matsubara frequency will do
|
||||
# G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
|
||||
# for bname,gf in G_latt_iw:
|
||||
# deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
|
||||
# dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
|
||||
# if dm_type == 'vasp':
|
||||
## In 'vasp'-mode subtract the DFT density matrix
|
||||
# diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
|
||||
# deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
|
||||
# dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
|
||||
|
||||
# mpi reduce:
|
||||
for bname in deltaH:
|
||||
for ik in range(self.n_k):
|
||||
deltaH[bname][ik] = mpi.all_reduce(mpi.world, deltaH[bname][ik], lambda x,y : x+y)
|
||||
mpi.barrier()
|
||||
|
||||
# now save to file:
|
||||
if mpi.is_master_node():
|
||||
with open(filename, 'w') as f:
|
||||
f.write("H %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
|
||||
for sp in spn:
|
||||
for ik in xrange(self.n_k):
|
||||
ib1 = band_window[0][ik, 0]
|
||||
ib2 = band_window[0][ik, 1]
|
||||
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
|
||||
for inu in xrange(self.n_orbitals[ik, 0]):
|
||||
for imu in xrange(self.n_orbitals[ik, 0]):
|
||||
if self.SP == 0:
|
||||
valre = (deltaH['up'][ik][inu, imu].real + deltaH['down'][ik][inu, imu].real) / 2.0
|
||||
valim = (deltaH['up'][ik][inu, imu].imag + deltaH['down'][ik][inu, imu].imag) / 2.0
|
||||
else:
|
||||
valre = deltaH[sp][ik][inu, imu].real
|
||||
valim = deltaH[sp][ik][inu, imu].imag
|
||||
|
||||
f.write(" %.14f %.14f"%(valre, valim))
|
||||
f.write("\n")
|
||||
return deltaH
|
||||
|
||||
def dmft_cycle():
|
||||
lda_filename = 'vasp'
|
||||
beta = 400
|
||||
U_int = 4.00
|
||||
J_hund = 0.70
|
||||
Loops = 1 # Number of DMFT sc-loops
|
||||
Mix = 1.0 # Mixing factor in QMC
|
||||
DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
|
||||
DC_Mix = 1.0 # 1.0 ... all from imp; 0.0 ... all from Gloc
|
||||
useBlocs = False # use bloc structure from LDA input
|
||||
useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
|
||||
chemical_potential_init=-0.0 # initial chemical potential
|
||||
|
||||
use_dudarev = True
|
||||
|
||||
HDFfilename = lda_filename+'.h5'
|
||||
|
||||
# Convert DMFT input:
|
||||
# Can be commented after the first run
|
||||
Converter = VaspConverter(filename=lda_filename)
|
||||
Converter.convert_dft_input()
|
||||
|
||||
#check if there are previous runs:
|
||||
previous_runs = 0
|
||||
previous_present = False
|
||||
|
||||
if mpi.is_master_node():
|
||||
with HDFArchive(HDFfilename,'a') as ar:
|
||||
if 'iterations' in ar:
|
||||
previous_present = True
|
||||
previous_runs = ar['iterations']
|
||||
else:
|
||||
previous_runs = 0
|
||||
previous_present = False
|
||||
|
||||
mpi.barrier()
|
||||
previous_runs = mpi.bcast(previous_runs)
|
||||
previous_present = mpi.bcast(previous_present)
|
||||
|
||||
# Init the SumK class
|
||||
#SK=SumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
|
||||
SK=TestSumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
|
||||
|
||||
Norb = SK.corr_shells[0]['dim']
|
||||
l = SK.corr_shells[0]['l']
|
||||
|
||||
# Init the Hubbard-I solver:
|
||||
S = Solver(beta = beta, l = l, dudarev=use_dudarev)
|
||||
|
||||
# DEBUG
|
||||
#SK.put_Sigma(Sigma_imp=[S.Sigma])
|
||||
#dH = SK.calc_hamiltonian_correction(filename='GAMMA')
|
||||
# END DEBUG
|
||||
|
||||
chemical_potential=chemical_potential_init
|
||||
# load previous data: old self-energy, chemical potential, DC correction
|
||||
if (previous_present):
|
||||
mpi.report("Using stored data for initialisation")
|
||||
if (mpi.is_master_node()):
|
||||
with HDFArchive(HDFfilename,'a') as ar:
|
||||
S.Sigma <<= ar['SigmaF']
|
||||
things_to_load=['chemical_potential','dc_imp']
|
||||
old_data=SK.load(things_to_load)
|
||||
chemical_potential=old_data[0]
|
||||
SK.dc_imp=old_data[1]
|
||||
S.Sigma = mpi.bcast(S.Sigma)
|
||||
chemical_potential=mpi.bcast(chemical_potential)
|
||||
SK.dc_imp=mpi.bcast(SK.dc_imp)
|
||||
|
||||
|
||||
# DMFT loop:
|
||||
for Iteration_Number in range(1,Loops+1):
|
||||
|
||||
itn = Iteration_Number + previous_runs
|
||||
|
||||
# put Sigma into the SumK class:
|
||||
S.Sigma.zero() # !!!!
|
||||
SK.put_Sigma(Sigma_imp = [ S.Sigma ])
|
||||
|
||||
# Compute the SumK, possibly fixing mu by dichotomy
|
||||
if SK.density_required and (Iteration_Number > 1):
|
||||
chemical_potential = SK.calc_mu( precision = 0.001 )
|
||||
else:
|
||||
mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=chemical_potential))
|
||||
|
||||
# Density:
|
||||
S.G <<= SK.extract_G_loc()[0]
|
||||
mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
|
||||
|
||||
# calculated DC at the first run to have reasonable initial non-interacting atomic level positions
|
||||
if ((Iteration_Number==1)and(previous_present==False)):
|
||||
if use_dudarev:
|
||||
dc_value_init = 0.0
|
||||
else:
|
||||
dc_value_init=U_int/2.0
|
||||
dm=S.G.density()
|
||||
SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_dc_value=dc_value_init)
|
||||
|
||||
# calculate non-interacting atomic level positions:
|
||||
# eal = SK.eff_atomic_levels()[0]
|
||||
# S.set_atomic_levels( eal = eal )
|
||||
|
||||
# solve it:
|
||||
corr_energy, dft_dc = S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1)
|
||||
|
||||
# Now mix Sigma and G:
|
||||
if ((itn>1)or(previous_present)):
|
||||
if (mpi.is_master_node()and (Mix<1.0)):
|
||||
with HDFArchive(HDFfilename,'r') as ar:
|
||||
mpi.report("Mixing Sigma and G with factor %s"%Mix)
|
||||
if ('SigmaF' in ar):
|
||||
S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF']
|
||||
if ('GF' in ar):
|
||||
S.G <<= Mix * S.G + (1.0-Mix) * ar['GF']
|
||||
S.G = mpi.bcast(S.G)
|
||||
S.Sigma = mpi.bcast(S.Sigma)
|
||||
|
||||
# after the Solver has finished, set new double counting:
|
||||
# dm = S.G.density()
|
||||
# if not use_dudarev:
|
||||
# SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type )
|
||||
|
||||
# correlation energy calculations:
|
||||
correnerg = 0.5 * (S.G * S.Sigma).total_density()
|
||||
mpi.report("Corr. energy = %s"%correnerg)
|
||||
|
||||
# store the impurity self-energy, GF as well as correlation energy in h5
|
||||
if (mpi.is_master_node()):
|
||||
with HDFArchive(HDFfilename,'a') as ar:
|
||||
ar['iterations'] = itn
|
||||
ar['chemical_cotential%s'%itn] = chemical_potential
|
||||
ar['SigmaF'] = S.Sigma
|
||||
ar['GF'] = S.G
|
||||
ar['correnerg%s'%itn] = correnerg
|
||||
ar['DCenerg%s'%itn] = SK.dc_energ
|
||||
|
||||
#Save essential SumkDFT data:
|
||||
things_to_save=['chemical_potential','dc_energ','dc_imp']
|
||||
SK.save(things_to_save)
|
||||
|
||||
# print out occupancy matrix of Ce 4f
|
||||
# mpi.report("Orbital densities of impurity Green function:")
|
||||
# for s in dm:
|
||||
# mpi.report("Block %s: "%s)
|
||||
# for ii in range(len(dm[s])):
|
||||
# str = ''
|
||||
# for jj in range(len(dm[s])):
|
||||
# if (dm[s][ii,jj].real>0):
|
||||
# str += " %.4f"%(dm[s][ii,jj].real)
|
||||
# else:
|
||||
# str += " %.4f"%(dm[s][ii,jj].real)
|
||||
# mpi.report(str)
|
||||
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
|
||||
|
||||
|
||||
# find exact chemical potential
|
||||
#if (SK.density_required):
|
||||
# SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
|
||||
# SK.chemical_potential = SK.calc_mu( precision = 0.01 )
|
||||
|
||||
#dN, d = SK.calc_density_correction(filename='GAMMA', dm_type='vasp')
|
||||
#mpi.report("Trace of Density Matrix: %s"%d)
|
||||
mpi.report("Storing Hamiltonian correction GAMMA...")
|
||||
SK.put_Sigma(Sigma_imp=[S.Sigma])
|
||||
dH = SK.calc_hamiltonian_correction(filename='GAMMA')
|
||||
# shutil.copyfile('GAMMA', 'it%i.GAMMA'%(itn))
|
||||
|
||||
# store correlation energy contribution to be read by Wien2ki and then included to DFT+DMFT total energy
|
||||
if (mpi.is_master_node()):
|
||||
with HDFArchive(HDFfilename) as ar:
|
||||
itn = ar['iterations']
|
||||
correnerg = ar['correnerg%s'%itn]
|
||||
DCenerg = ar['DCenerg%s'%itn]
|
||||
correnerg -= DCenerg[0]
|
||||
f=open(lda_filename+'.qdmft','a')
|
||||
f.write("%.16f\n"%correnerg)
|
||||
f.close()
|
||||
|
||||
return corr_energy, dft_dc
|
||||
|
||||
if __name__ == '__main__':
|
||||
dmft_cycle()
|
||||
|
|
@ -1,11 +0,0 @@
|
|||
*
|
||||
!.gitignore
|
||||
!INCAR
|
||||
!POSCAR
|
||||
!POTCAR
|
||||
!KPOINTS
|
||||
!plo.cfg
|
||||
!run_plovasp.sh
|
||||
!conv_example.py
|
||||
!rotations
|
||||
!rot_dz2_dx2
|
||||
|
|
@ -1,39 +0,0 @@
|
|||
ADDGRID = TRUE
|
||||
ALGO = NORMAL
|
||||
EDIFF = 0.00001
|
||||
EDIFFG = -0.01
|
||||
EMAX = 13.0
|
||||
EMIN = -3.0
|
||||
ENCUT = 500
|
||||
IBRION = -1
|
||||
ICHARG = 11
|
||||
ISIF = 2
|
||||
ISMEAR = -5
|
||||
ISPIN = 1
|
||||
ISYM = -1
|
||||
LASPH = FALSE
|
||||
LCHARG = TRUE
|
||||
LDAU = FALSE
|
||||
LDAUJ = 0 0 0 0 0
|
||||
LDAUL = 0 2 0 0 0
|
||||
LDAUTYPE = 2
|
||||
LDAUU = 0 6.2 0 0 0
|
||||
LMAXMIX = 6
|
||||
LORBIT = 0
|
||||
LREAL = Auto
|
||||
ROPT = 1e-3 1e-3 1e-3
|
||||
LWAVE = FALSE
|
||||
NEDOS = 1001
|
||||
NELM = 100
|
||||
NELMDL = -9
|
||||
NELMIN = 5
|
||||
NPAR = 4
|
||||
NSIM = 1
|
||||
NSW = 0
|
||||
POTIM = 0.1
|
||||
PREC = Accurate
|
||||
SIGMA = 0.05
|
||||
SMASS = 0.5
|
||||
SYMPREC = 1.0e-12
|
||||
|
||||
LOCPROJ = 5 6 7 8 : d : Hy
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
Automatic kpoint scheme
|
||||
0
|
||||
Gamma
|
||||
5 5 3
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
LuNiO3, exp. low-T
|
||||
1.0
|
||||
5.1234998703 0.0000000000 0.0000000000
|
||||
0.0000000000 5.5089001656 0.0000000000
|
||||
-0.0166880521 0.0000000000 7.3551808822
|
||||
Lu Ni O
|
||||
4 4 12
|
||||
Direct
|
||||
0.977199972 0.077000000 0.252999991
|
||||
0.022800028 0.922999978 0.746999979
|
||||
0.522800028 0.577000022 0.247000009
|
||||
0.477199972 0.423000008 0.753000021
|
||||
0.500000000 0.000000000 0.000000000
|
||||
0.000000000 0.500000000 0.500000000
|
||||
0.500000000 0.000000000 0.500000000
|
||||
0.000000000 0.500000000 0.000000000
|
||||
0.110100001 0.462700009 0.244100004
|
||||
0.889899969 0.537299991 0.755900025
|
||||
0.389899999 0.962700009 0.255899996
|
||||
0.610100031 0.037299991 0.744099975
|
||||
0.693300009 0.313699991 0.053900000
|
||||
0.306699991 0.686300039 0.946099997
|
||||
0.806699991 0.813699961 0.446099997
|
||||
0.193300009 0.186300009 0.553900003
|
||||
0.185100004 0.201600000 0.943799973
|
||||
0.814899981 0.798399985 0.056200027
|
||||
0.314899981 0.701600015 0.556200027
|
||||
0.685100019 0.298399985 0.443799973
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,12 +0,0 @@
|
|||
[General]
|
||||
DOSMESH = 401
|
||||
|
||||
[Shell 1]
|
||||
EWINDOW = -8.6 2.7
|
||||
#EWINDOW = -0.6 2.7
|
||||
LSHELL = 2
|
||||
IONS = 5..8
|
||||
NORMION = True
|
||||
TRANSFILE = rotations
|
||||
#TRANSFILE = rot_dz2_dx2
|
||||
NORMALIZE = True
|
||||
|
|
@ -1,18 +0,0 @@
|
|||
# Rotations to the local frames
|
||||
|
||||
# Atom 1, Euler angles: 44.57639968 16.85783381 -11.27701880
|
||||
0.02793493813 -0.09400347991 0.87385177595 -0.47142545676 0.06726141074
|
||||
-0.87830274777 -0.28356545426 0.00107688253 0.06072810673 0.38011294872
|
||||
|
||||
# Atom 2, Euler angles: 47.67919913 16.85783423 11.27702026
|
||||
-0.02793494291 0.09400349403 0.87385176987 -0.47142546468 0.06726141256
|
||||
-0.84565748799 -0.27808790765 -0.00680142052 -0.08187952499 -0.44808482755
|
||||
|
||||
# Atom 3, Euler angles:-133.92682764 17.32502138 165.73386489
|
||||
0.03668336697 0.12133769713 0.86698007031 0.47720500874 0.06747170631
|
||||
-0.85519672598 0.29102955535 -0.00287628125 -0.06301365234 0.42421853379
|
||||
|
||||
# Atom 4, Euler angles:-135.60182669 17.32502172 -165.73386378
|
||||
-0.03668337098 -0.12133770845 0.86698006522 0.47720501465 0.06747170748
|
||||
-0.84809293709 0.29001625337 0.00161325481 0.06758002057 -0.43824568573
|
||||
|
||||
|
|
@ -1,30 +0,0 @@
|
|||
# Rotations to the local frames
|
||||
|
||||
# Atom 1, Euler angles: 44.57639968 16.85783381 -11.27701880
|
||||
0.38045444863 0.05847209742 0.07282374934 0.27130935338 0.87916060116
|
||||
-0.26548421431 0.78270891327 0.33738843326 0.43922773606 -0.10066245585
|
||||
0.02793493813 -0.09400347991 0.87385177595 -0.47142545676 0.06726141074
|
||||
0.11214484811 -0.54286439359 0.34241445738 0.71241867354 -0.26064104853
|
||||
-0.87830274777 -0.28356545426 0.00107688253 0.06072810673 0.38011294872
|
||||
|
||||
# Atom 2, Euler angles: 47.67919913 16.85783423 11.27702026
|
||||
-0.44835779602 -0.08059449489 0.07251344244 0.26569233153 0.84653954413
|
||||
-0.10160916307 0.51163832768 0.35542791179 0.72915358208 -0.26440093208
|
||||
-0.02793494291 0.09400349403 0.87385176987 -0.47142546468 0.06726141256
|
||||
0.26969226294 -0.80346823694 0.32365046757 0.41084840677 -0.09032628101
|
||||
-0.84565748799 -0.27808790765 -0.00680142052 -0.08187952499 -0.44808482755
|
||||
|
||||
# Atom 3, Euler angles:-133.92682764 17.32502138 165.73386489
|
||||
0.42474364412 -0.05919623743 0.07674521204 -0.27806951500 0.85608186496
|
||||
0.27929719729 0.78784591591 -0.35463192502 0.41101628718 0.08120158747
|
||||
0.03668336697 0.12133769713 0.86698007031 0.47720500874 0.06747170631
|
||||
-0.09422739956 -0.52571244938 -0.34158989281 0.72252862625 0.27571062205
|
||||
-0.85519672598 0.29102955535 -0.00287628125 -0.06301365234 0.42421853379
|
||||
|
||||
# Atom 4, Euler angles:-135.60182669 17.32502172 -165.73386378
|
||||
-0.43875616246 0.06397634469 0.07678214913 -0.27699516110 0.84898659561
|
||||
0.09192658044 0.51921347838 -0.34449569877 0.72588538389 0.27636930122
|
||||
-0.03668337098 -0.12133770845 0.86698006522 0.47720501465 0.06747170748
|
||||
-0.28006291089 -0.79214400488 -0.35180986204 0.40505852046 0.07893071297
|
||||
-0.84809293709 0.29001625337 0.00161325481 0.06758002057 -0.43824568573
|
||||
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
PLOVASP_PATH=/path/to/triqs1.2/applications/dft_tools/src
|
||||
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
*
|
||||
!example.cfg
|
||||
!INCAR
|
||||
!KPOINTS
|
||||
!POSCAR
|
||||
!POTCAR
|
||||
!conv_example.py
|
||||
|
|
@ -1,42 +0,0 @@
|
|||
ADDGRID = TRUE
|
||||
ALGO = NORMAL
|
||||
EDIFF = 0.00001
|
||||
EDIFFG = -0.002
|
||||
EMAX = 13.0
|
||||
EMIN = -2.0
|
||||
ENCUT = 550
|
||||
ENAUG = 1100
|
||||
IBRION = -1
|
||||
ICHARG = 11
|
||||
ISIF = 4
|
||||
ISMEAR = -5
|
||||
ISPIN = 1
|
||||
ISYM = -1
|
||||
LASPH = TRUE
|
||||
LCHARG = TRUE
|
||||
LDAU = False
|
||||
LDAUJ = 0 1.0 0 0 0
|
||||
LDAUL = -1 2 -1
|
||||
LDAUTYPE = 1
|
||||
LDAUU = 0 5.0 0 0 0
|
||||
LDAUPRINT = 1
|
||||
LMAXMIX = 6
|
||||
LORBIT = 0
|
||||
LREAL = False
|
||||
ROPT = 1e-3 1e-3 1e-3
|
||||
LWAVE = FALSE
|
||||
NBANDS = 16
|
||||
NEDOS = 501
|
||||
NELM = 100
|
||||
NELMDL = -7
|
||||
NELMIN = 7
|
||||
NPAR = 4
|
||||
NSIM = 8
|
||||
NSW = 0
|
||||
POTIM = 0.4
|
||||
PREC = Accurate
|
||||
SIGMA = 0.1
|
||||
SMASS = 0.5
|
||||
SYMPREC = 1.0e-6
|
||||
|
||||
LOCPROJ = 1 : d : Hy
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
Automatic kpoint scheme
|
||||
0
|
||||
Gamma
|
||||
2 2 2
|
||||
|
|
@ -1,10 +0,0 @@
|
|||
V SF test
|
||||
2.9878
|
||||
-0.5 0.5 0.5
|
||||
0.5 -0.5 0.5
|
||||
0.5 0.5 -0.5
|
||||
V
|
||||
1
|
||||
Direct
|
||||
0.0 0.0 0.0
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
from vasp_converter import VaspConverter
|
||||
|
||||
if __name__ == '__main__':
|
||||
conv = VaspConverter('vasp')
|
||||
conv.convert_dft_input()
|
||||
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
[Shell 1]
|
||||
LSHELL = 2
|
||||
IONS = 1
|
||||
EWINDOW = -15.0 5.0
|
||||
|
||||
|
||||
|
|
@ -1,9 +0,0 @@
|
|||
*
|
||||
!.gitignore
|
||||
!INCAR
|
||||
!POSCAR
|
||||
!POTCAR
|
||||
!KPOINTS
|
||||
!plo.cfg
|
||||
!run_plovasp.sh
|
||||
!conv_example.py
|
||||
|
|
@ -1,43 +0,0 @@
|
|||
ADDGRID = TRUE
|
||||
ALGO = NORMAL
|
||||
EDIFF = 0.00001
|
||||
EDIFFG = -0.002
|
||||
EMAX = 13.0
|
||||
EMIN = -2.0
|
||||
ENCUT = 550
|
||||
ENAUG = 1100
|
||||
IBRION = -1
|
||||
ICHARG = 11
|
||||
ISIF = 4
|
||||
ISMEAR = -5
|
||||
ISPIN = 1
|
||||
ISYM = -1
|
||||
LASPH = TRUE
|
||||
LCHARG = TRUE
|
||||
LDAU = False
|
||||
LDAUJ = 0 1.0 0 0 0
|
||||
LDAUL = -1 2 -1
|
||||
LDAUTYPE = 1
|
||||
LDAUU = 0 5.0 0 0 0
|
||||
LDAUPRINT = 1
|
||||
LMAXMIX = 6
|
||||
LORBIT = 0
|
||||
LREAL = False
|
||||
ROPT = 1e-3 1e-3 1e-3
|
||||
LWAVE = FALSE
|
||||
NBANDS = 32
|
||||
NEDOS = 501
|
||||
NELM = 100
|
||||
NELMDL = -7
|
||||
NELMIN = 7
|
||||
NPAR = 1
|
||||
NSIM = 8
|
||||
NSW = 0
|
||||
POTIM = 0.4
|
||||
PREC = Accurate
|
||||
SIGMA = 0.1
|
||||
SMASS = 0.5
|
||||
SYMPREC = 1.0e-6
|
||||
|
||||
# Projectors
|
||||
LOCPROJ = 2 : d : Hy
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
Automatic kpoint scheme
|
||||
0
|
||||
Gamma
|
||||
3 3 3
|
||||
|
|
@ -1,13 +0,0 @@
|
|||
SrVO3
|
||||
3.841
|
||||
1.0000000000000000 0.0000000000000000 0.0000000000000000
|
||||
0.0000000000000000 1.0000000000000000 0.0000000000000000
|
||||
0.0000000000000000 0.0000000000000000 1.0000000000000000
|
||||
Sr V O
|
||||
1 1 3
|
||||
Direct
|
||||
0.5 0.5 0.5
|
||||
0.0 0.0 0.0
|
||||
0.5 0.0 0.0
|
||||
0.0 0.5 0.0
|
||||
0.0 0.0 0.5
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,51 +0,0 @@
|
|||
|
||||
import numpy as np
|
||||
from pytriqs.gf import *
|
||||
#from sumk_dft import SumkDFT
|
||||
from sumk_dft_tools import SumkDFTTools
|
||||
from converters.vasp_converter import VaspConverter
|
||||
|
||||
np.set_printoptions(suppress=True)
|
||||
|
||||
def density_matrix_and_overlap(sk):
|
||||
glist = [GfImFreq(indices=inds,beta=1000.0) for bl, inds in sk.gf_struct_sumk[0]]
|
||||
sigma = BlockGf(name_list=[bl for bl, inds in sk.gf_struct_sumk[0]], block_list=glist)
|
||||
sigma.zero()
|
||||
sk.put_Sigma(Sigma_imp=[sigma])
|
||||
|
||||
print "Overlap matrix:"
|
||||
dm_blocks = sk.check_projectors()
|
||||
print dm_blocks[0].real
|
||||
|
||||
sk.calc_mu(precision=0.001)
|
||||
|
||||
print
|
||||
print "Denisty matrix:"
|
||||
dm_blocks = sk.density_matrix(method='using_gf', beta=1000.0)
|
||||
ntot = 0.0
|
||||
for bl in dm_blocks[0]:
|
||||
print
|
||||
print bl
|
||||
print dm_blocks[0][bl].real
|
||||
ntot += dm_blocks[0][bl].real.trace()
|
||||
|
||||
print
|
||||
print " Impurity density:", ntot
|
||||
|
||||
def density_of_states(sk):
|
||||
glist = [GfReFreq(indices=inds, window=(-10.0, 10.0), n_points=2000) for bl, inds in sk.gf_struct_sumk[0]]
|
||||
sigma = BlockGf(name_list=[bl for bl, inds in sk.gf_struct_sumk[0]], block_list=glist)
|
||||
sigma.zero()
|
||||
sk.put_Sigma(Sigma_imp=[sigma])
|
||||
|
||||
print " Evaluating DOS..."
|
||||
sk.dos_wannier_basis(broadening=0.03, with_dc=False)
|
||||
|
||||
if __name__ == '__main__':
|
||||
conv = VaspConverter('vasp')
|
||||
conv.convert_dft_input()
|
||||
|
||||
sk = SumkDFTTools(hdf_file='vasp.h5')
|
||||
|
||||
density_matrix_and_overlap(sk)
|
||||
# density_of_states(sk)
|
||||
|
|
@ -1,16 +0,0 @@
|
|||
[General]
|
||||
DOSMESH = -3.0 5.0 4001
|
||||
|
||||
[Group 1]
|
||||
SHELLS = 1
|
||||
NORMALIZE = True
|
||||
EWINDOW = -1.45 1.8
|
||||
|
||||
[Shell 1]
|
||||
LSHELL = 2
|
||||
IONS = 2
|
||||
|
||||
TRANSFORM = 1.0 0.0 0.0 0.0 0.0
|
||||
0.0 1.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 1.0 0.0
|
||||
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
PLOVASP_PATH=/path/to/triqs1.2/applications/dft_tools/src
|
||||
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@
|
||||
|
|
@ -1,42 +0,0 @@
|
|||
ADDGRID = TRUE
|
||||
ALGO = NORMAL
|
||||
EDIFF = 0.00001
|
||||
EDIFFG = -0.002
|
||||
EMAX = 13.0
|
||||
EMIN = -2.0
|
||||
ENCUT = 550
|
||||
ENAUG = 1100
|
||||
IBRION = -1
|
||||
ICHARG = 11
|
||||
ISIF = 4
|
||||
ISMEAR = -5
|
||||
ISPIN = 1
|
||||
ISYM = -1
|
||||
LASPH = TRUE
|
||||
LCHARG = TRUE
|
||||
LDAU = False
|
||||
LDAUJ = 0 1.0 0 0 0
|
||||
LDAUL = -1 2 -1
|
||||
LDAUTYPE = 1
|
||||
LDAUU = 0 5.0 0 0 0
|
||||
LDAUPRINT = 1
|
||||
LMAXMIX = 6
|
||||
LORBIT = 0
|
||||
LREAL = False
|
||||
ROPT = 1e-3 1e-3 1e-3
|
||||
LWAVE = FALSE
|
||||
NBANDS = 32
|
||||
NEDOS = 501
|
||||
NELM = 100
|
||||
NELMDL = -7
|
||||
NELMIN = 7
|
||||
NPAR = 4
|
||||
NSIM = 8
|
||||
NSW = 0
|
||||
POTIM = 0.4
|
||||
PREC = Accurate
|
||||
SIGMA = 0.1
|
||||
SMASS = 0.5
|
||||
SYMPREC = 1.0e-6
|
||||
|
||||
LOCPROJ = 1 2 : d : Hy
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
Automatic kpoint scheme
|
||||
0
|
||||
Gamma
|
||||
2 2 2
|
||||
|
|
@ -1,11 +0,0 @@
|
|||
V SF test
|
||||
2.9878
|
||||
1.0 0.0 0.0
|
||||
0.0 1.0 0.0
|
||||
0.0 0.0 1.0
|
||||
V
|
||||
2
|
||||
Direct
|
||||
0.0 0.0 0.0
|
||||
0.5 0.5 0.5
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,7 +0,0 @@
|
|||
|
||||
[Shell 1]
|
||||
LSHELL = 2
|
||||
IONS = 1
|
||||
EWINDOW = -15.0 5.0
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user