Commit cf4e21f1 authored by Thorbjørn Skovhus's avatar Thorbjørn Skovhus

Merge branch 'get-general-component' into 'versatile-chi0-integration'

Get general component

See merge request !598
parents cd3fb490 a5deba64
Pipeline #97565757 passed with stage
in 3 minutes and 16 seconds
import sys
from time import ctime
from pathlib import Path
import pickle
import numpy as np
......@@ -108,7 +109,7 @@ class FourComponentSusceptibilityTensor:
frequencies,
txt=txt)
write_macroscopic_component(omega_w / Hartree, chiks_w, chi_w,
write_macroscopic_component(omega_w, chiks_w, chi_w,
filename, self.world)
return omega_w, chiks_w, chi_w
......@@ -128,7 +129,7 @@ class FourComponentSusceptibilityTensor:
omega_w, chiks_w, chi_w : nd.array, nd.array, nd.array
omega_w: frequencies in eV
chiks_w: macroscopic dynamic susceptibility (Kohn-Sham system)
chi_w: macroscopic(to be generalized?) dynamic susceptibility
chi_w: macroscopic dynamic susceptibility
"""
(pd, wd,
chiks_wGG, chi_wGG) = self.calculate_component(spincomponent, q_c,
......@@ -145,6 +146,66 @@ class FourComponentSusceptibilityTensor:
return omega_w, chiks_w, chi_w
def get_component(self, spincomponent, q_c, frequencies,
store_ecut=50, filename=None, txt=None):
"""Calculates a specific spin component of the
susceptibility tensor and writes it to a file.
Parameters
----------
spincomponent, q_c,
frequencies : see gpaw.response.chiks, gpaw.response.kslrf
store_ecut : float
Energy cutoff for the reduced plane wave representation.
The susceptibility is returned/saved in the reduced representation.
filename : str
Save chiks_w and chi_w to pickle file of given name.
Defaults to:
'chi%sGG_q«%+d-%+d-%+d».pckl' % (spincomponent,
*tuple((q_c * kd.N_c).round()))
txt : str
Save output of the calculation of this specific component into
a file with the filename of the given input.
Returns
-------
omega_w, G_Gc, chiks_wGG, chi_wGG : nd.array(s)
omega_w: frequencies in eV
G_Gc : plane wave repr. as coordinates on the reciprocal lattice
chiks_wGG: dynamic susceptibility (Kohn-Sham system)
chi_wGG: dynamic susceptibility
"""
if filename is None:
tup = (spincomponent,
*tuple((q_c * self.calc.wfs.kd.N_c).round()))
filename = 'chi%sGG_q«%+d-%+d-%+d».pckl' % tup
(pd, wd,
chiks_wGG, chi_wGG) = self.calculate_component(spincomponent, q_c,
frequencies, txt=txt)
# Get frequencies in eV
omega_w = wd.get_data() * Hartree
# Get susceptibility in a reduced plane wave representation
mask_G = get_pw_reduction_map(pd, store_ecut)
chiks_wGG = np.ascontiguousarray(chiks_wGG[:, mask_G, :][:, :, mask_G])
chi_wGG = np.ascontiguousarray(chi_wGG[:, mask_G, :][:, :, mask_G])
# Get reduced plane wave repr. as coordinates on the reciprocal lattice
G_Gc = get_pw_coordinates(pd)[mask_G]
# Gather susceptibilities for all frequencies
chiks_wGG = self.gather(chiks_wGG, wd)
chi_wGG = self.gather(chi_wGG, wd)
# Write susceptibilities to a pickle file
write_component(omega_w, G_Gc, chiks_wGG, chi_wGG,
filename, self.world)
return omega_w, G_Gc, chiks_wGG, chi_wGG
def calculate_component(self, spincomponent, q_c, frequencies, txt=None):
"""Calculate a single component of the susceptibility tensor.
......@@ -287,6 +348,29 @@ class FourComponentSusceptibilityTensor:
world.all_gather(b_w, A_w)
return A_w[:nw]
def gather(self, A_wGG, wd):
"""Gather a full susceptibility array to root."""
# Allocate arrays to gather (all need to be the same shape)
shape = (self.mynw,) + A_wGG.shape[1:]
tmp_wGG = np.empty(shape, dtype=A_wGG.dtype)
tmp_wGG[:self.w2 - self.w1] = A_wGG
# Allocate array for the gathered data
if self.world.rank == 0:
# Make room for all frequencies
shape = (self.mynw * self.world.size,) + A_wGG.shape[1:]
allA_wGG = np.empty(shape, dtype=A_wGG.dtype)
else:
allA_wGG = None
self.world.gather(tmp_wGG, 0, allA_wGG)
# Return array for w indeces on frequency grid
if allA_wGG is not None:
allA_wGG = allA_wGG[:len(wd), :, :]
return allA_wGG
@timer('Distribute frequencies')
def distribute_frequencies(self, chiks_wGG):
"""Distribute frequencies to all cores."""
......@@ -335,12 +419,87 @@ class FourComponentSusceptibilityTensor:
self.fd.close()
def get_pw_reduction_map(pd, ecut):
"""Get a mask to reduce the plane wave representation.
Please remark, that the response code currently works with one q-vector
at a time, at thus only a single plane wave representation at a time.
Returns
-------
mask_G : nd.array (dtype=bool)
Mask which reduces the representation
"""
assert ecut is not None
ecut /= Hartree
assert ecut <= pd.ecut
# List of all plane waves
G_Gv = np.array([pd.G_Qv[Q] for Q in pd.Q_qG[0]])
if pd.gammacentered:
mask_G = ((G_Gv ** 2).sum(axis=1) <= 2 * ecut)
else:
mask_G = (((G_Gv + pd.K_qv[0]) ** 2).sum(axis=1) <= 2 * ecut)
return mask_G
def get_pw_coordinates(pd):
"""Get the reciprocal lattice vector coordinates corresponding to a
givne plane wave basis.
Please remark, that the response code currently works with one q-vector
at a time, at thus only a single plane wave representation at a time.
Returns
-------
G_Gc : nd.array (dtype=int)
Coordinates on the reciprocal lattice
"""
# List of all plane waves
G_Gv = np.array([pd.G_Qv[Q] for Q in pd.Q_qG[0]])
# Use cell to get coordinates
B_cv = 2.0 * np.pi * pd.gd.icell_cv
return np.round(np.dot(G_Gv, np.linalg.inv(B_cv))).astype(int)
def write_macroscopic_component(omega_w, chiks_w, chi_w, filename, world):
"""Write the spatially averaged dynamic susceptibility."""
assert isinstance(filename, str)
if world.rank == 0:
with Path(filename).open('w') as fd:
for omega, chiks, chi in zip(omega_w * Hartree, chiks_w, chi_w):
for omega, chiks, chi in zip(omega_w, chiks_w, chi_w):
print('%.6f, %.6f, %.6f, %.6f, %.6f' %
(omega, chiks.real, chiks.imag, chi.real, chi.imag),
file=fd)
def read_macroscopic_component(filename):
"""Read a stored macroscopic susceptibility file"""
d = np.loadtxt(filename, delimiter=', ')
omega_w = d[:, 0]
chiks_w = np.array(d[:, 1], complex)
chiks_w.imag = d[:, 2]
chi_w = np.array(d[:, 3], complex)
chi_w.imag = d[:, 4]
return omega_w, chiks_w, chi_w
def write_component(omega_w, G_Gc, chiks_wGG, chi_wGG, filename, world):
"""Write the dynamic susceptibility as a pickle file."""
assert isinstance(filename, str)
if world.rank == 0:
with open(filename, 'wb') as fd:
pickle.dump((omega_w, G_Gc, chiks_wGG, chi_wGG), fd)
def read_component(filename):
"""Read a stored susceptibility component file"""
assert isinstance(filename, str)
with open(filename, 'rb') as fd:
omega_w, G_Gc, chiks_wGG, chi_wGG = pickle.load(fd)
return omega_w, G_Gc, chiks_wGG, chi_wGG
......@@ -42,6 +42,28 @@ class TransverseMagneticSusceptibility(FCST):
frequencies, filename=filename,
txt=txt)
def get_component(self, spincomponent, q_c, frequencies,
store_ecut=50, filename=None, txt=None):
"""Calculates a specific spin component of the
transverse magnetic susceptibility and writes it to a file.
Parameters
----------
spincomponent : str
'+-': calculate chi+-, '-+: calculate chi-+
q_c, frequencies,
store_ecut, filename, txt : see gpaw.response.susceptibility
Returns
-------
see gpaw.response.susceptibility
"""
assert spincomponent in ['+-', '-+']
return FCST.get_component(self, spincomponent, q_c,
frequencies, store_ecut=store_ecut,
filename=filename, txt=txt)
def _calculate_component(self, spincomponent, pd, wd):
"""Calculate a transverse magnetic susceptibility element.
......
......@@ -349,6 +349,7 @@ tests = [
'generic/Cu.py', # ~21s
'vdw/ts09.py', # ~21s
'response/na_plasmon.py', # ~22s
'response/two-aluminum_chi_RPA.py', # ~23s
'lcao/kpts_many_combinations.py', # ~23s
'fermilevel.py', # ~23s
'ralda/ralda_energy_H2.py', # ~23s
......
......@@ -19,6 +19,7 @@ from ase.parallel import parprint
from gpaw import GPAW, PW
from gpaw.response.tms import TransverseMagneticSusceptibility
from gpaw.response.susceptibility import read_macroscopic_component
from gpaw.test import findpeak, equal
from gpaw.mpi import world
......@@ -107,32 +108,32 @@ parprint('Excited state calculations took', (t3 - t2) / 60, 'minutes')
world.barrier()
# Part 3: identify magnon peak in scattering functions
d1 = np.loadtxt('iron_dsus_G1.csv', delimiter=', ')
d2 = np.loadtxt('iron_dsus_G2.csv', delimiter=', ')
d3 = np.loadtxt('iron_dsus_G3.csv', delimiter=', ')
d4 = np.loadtxt('iron_dsus_G4.csv', delimiter=', ')
d5 = np.loadtxt('iron_dsus_G5.csv', delimiter=', ')
d6 = np.loadtxt('iron_dsus_G6.csv', delimiter=', ')
d7 = np.loadtxt('iron_dsus_G7.csv', delimiter=', ')
d8 = np.loadtxt('iron_dsus_G8.csv', delimiter=', ')
wpeak1, Ipeak1 = findpeak(d1[:, 0], d1[:, 4])
wpeak2, Ipeak2 = findpeak(d2[:, 0], d2[:, 4])
wpeak3, Ipeak3 = findpeak(d3[:, 0], d3[:, 4])
wpeak4, Ipeak4 = findpeak(d4[:, 0], d4[:, 4])
wpeak5, Ipeak5 = findpeak(d5[:, 0], d5[:, 4])
wpeak6, Ipeak6 = findpeak(d6[:, 0], d6[:, 4])
wpeak7, Ipeak7 = findpeak(d7[:, 0], d7[:, 4])
wpeak8, Ipeak8 = findpeak(d8[:, 0], d8[:, 4])
mw1 = (wpeak1 + d1[0, 0]) * 1000
mw2 = (wpeak2 + d2[0, 0]) * 1000
mw3 = (wpeak3 + d3[0, 0]) * 1000
mw4 = (wpeak4 + d4[0, 0]) * 1000
mw5 = (wpeak5 + d5[0, 0]) * 1000
mw6 = (wpeak6 + d6[0, 0]) * 1000
mw7 = (wpeak7 + d7[0, 0]) * 1000
mw8 = (wpeak8 + d8[0, 0]) * 1000
w1_w, chiks1_w, chi1_w = read_macroscopic_component('iron_dsus_G1.csv')
w2_w, chiks2_w, chi2_w = read_macroscopic_component('iron_dsus_G2.csv')
w3_w, chiks3_w, chi3_w = read_macroscopic_component('iron_dsus_G3.csv')
w4_w, chiks4_w, chi4_w = read_macroscopic_component('iron_dsus_G4.csv')
w5_w, chiks5_w, chi5_w = read_macroscopic_component('iron_dsus_G5.csv')
w6_w, chiks6_w, chi6_w = read_macroscopic_component('iron_dsus_G6.csv')
w7_w, chiks7_w, chi7_w = read_macroscopic_component('iron_dsus_G7.csv')
w8_w, chiks8_w, chi8_w = read_macroscopic_component('iron_dsus_G8.csv')
wpeak1, Ipeak1 = findpeak(w1_w, chi1_w.imag)
wpeak2, Ipeak2 = findpeak(w2_w, chi2_w.imag)
wpeak3, Ipeak3 = findpeak(w3_w, chi3_w.imag)
wpeak4, Ipeak4 = findpeak(w4_w, chi4_w.imag)
wpeak5, Ipeak5 = findpeak(w5_w, chi5_w.imag)
wpeak6, Ipeak6 = findpeak(w6_w, chi6_w.imag)
wpeak7, Ipeak7 = findpeak(w7_w, chi7_w.imag)
wpeak8, Ipeak8 = findpeak(w8_w, chi8_w.imag)
mw1 = (wpeak1 + w1_w[0]) * 1000
mw2 = (wpeak2 + w2_w[0]) * 1000
mw3 = (wpeak3 + w3_w[0]) * 1000
mw4 = (wpeak4 + w4_w[0]) * 1000
mw5 = (wpeak5 + w5_w[0]) * 1000
mw6 = (wpeak6 + w6_w[0]) * 1000
mw7 = (wpeak7 + w7_w[0]) * 1000
mw8 = (wpeak8 + w8_w[0]) * 1000
# Part 4: compare new results to test values
test_mw1 = 234.63 # meV
......
......@@ -16,6 +16,7 @@ from ase.parallel import parprint
from gpaw import GPAW, PW
from gpaw.response.tms import TransverseMagneticSusceptibility
from gpaw.response.susceptibility import read_macroscopic_component
from gpaw.test import findpeak, equal
from gpaw.mpi import world
......@@ -81,14 +82,14 @@ parprint('Excited state calculation took', (t3 - t2) / 60, 'minutes')
world.barrier()
# Part 3: identify magnon peaks in scattering function
d1 = np.loadtxt('iron_dsus_1.csv', delimiter=', ')
d2 = np.loadtxt('iron_dsus_2.csv', delimiter=', ')
w1_w, chiks1_w, chi1_w = read_macroscopic_component('iron_dsus_1.csv')
w2_w, chiks2_w, chi2_w = read_macroscopic_component('iron_dsus_2.csv')
wpeak1, Ipeak1 = findpeak(d1[:, 0], d1[:, 4])
wpeak2, Ipeak2 = findpeak(d2[:, 0], d2[:, 4])
wpeak1, Ipeak1 = findpeak(w1_w, chi1_w.imag)
wpeak2, Ipeak2 = findpeak(w2_w, chi2_w.imag)
mw1 = (wpeak1 + d1[0, 0]) * 1000
mw2 = (wpeak2 + d2[0, 0]) * 1000
mw1 = (wpeak1 + w1_w[0]) * 1000
mw2 = (wpeak2 + w2_w[0]) * 1000
# Part 4: compare new results to test values
test_fxcs = 1.033
......
......@@ -8,6 +8,7 @@ from gpaw import GPAW, PW, FermiDirac
from gpaw.test import findpeak, equal
from gpaw.response.df import DielectricFunction
from gpaw.response.susceptibility import FourComponentSusceptibilityTensor
from gpaw.response.susceptibility import read_macroscopic_component
from gpaw.mpi import size, world
assert size <= 4**3
......@@ -63,8 +64,8 @@ parprint('For excited state calc 2, it took', (t4 - t3) / 60, 'minutes')
# The two response codes should hold identical results
d1 = np.loadtxt('Si_chi1.csv', delimiter=',')
wpeak1, Ipeak1 = findpeak(d1[:, 0], -d1[:, 4])
d2 = np.loadtxt('Si_chi2.csv', delimiter=',')
wpeak2, Ipeak2 = findpeak(d2[:, 0], d2[:, 4])
w_w, chiks_w, chi_w = read_macroscopic_component('Si_chi2.csv')
wpeak2, Ipeak2 = findpeak(w_w, chi_w.imag)
equal(wpeak1, wpeak2, 0.02)
equal(Ipeak1, Ipeak2, 1.0)
import numpy as np
import time
from ase.build import bulk
from ase.parallel import parprint
from gpaw import GPAW, PW
from gpaw.test import findpeak, equal
from gpaw.response.susceptibility import FourComponentSusceptibilityTensor
from gpaw.response.susceptibility import read_component
from gpaw.mpi import size, world
assert size <= 4**3
# Ground state calculation
t1 = time.time()
a = 4.043
atoms1 = bulk('Al', 'fcc', a=a)
atoms2 = atoms1.repeat((2, 1, 1))
calc1 = GPAW(mode=PW(200),
nbands=4,
kpts=(8, 8, 8),
parallel={'domain': 1},
idiotproof=False, # allow uneven distribution of k-points
xc='LDA')
atoms1.set_calculator(calc1)
atoms1.get_potential_energy()
t2 = time.time()
calc2 = GPAW(mode=PW(200),
nbands=8,
kpts=(4, 8, 8),
parallel={'domain': 1},
idiotproof=False, # allow uneven distribution of k-points
xc='LDA')
atoms2.set_calculator(calc2)
atoms2.get_potential_energy()
t3 = time.time()
# Excited state calculation
q1_qc = [np.array([1 / 8., 0., 0.]), np.array([3 / 8., 0., 0.])]
q2_qc = [np.array([1 / 4., 0., 0.]), np.array([- 1 / 4., 0., 0.])]
w = np.linspace(0, 24, 241)
# Calculate susceptibility using Al
fcst = FourComponentSusceptibilityTensor(calc1, fxc='RPA',
eta=0.2, ecut=50,
disable_point_group=False,
disable_time_reversal=False)
for q, q_c in enumerate(q1_qc):
fcst.get_component('00', q_c, w, store_ecut=25,
filename='Al1_chiGG_q%d.pckl' % (q + 1))
world.barrier()
t4 = time.time()
# Calculate susceptibility using Al2
fcst = FourComponentSusceptibilityTensor(calc2, fxc='RPA',
eta=0.2, ecut=50,
disable_point_group=False,
disable_time_reversal=False)
for q, q_c in enumerate(q2_qc):
fcst.get_component('00', q_c, w, store_ecut=25,
filename='Al2_chiGG_q%d.pckl' % (q + 1))
world.barrier()
t5 = time.time()
parprint('')
parprint('Ground state calc 1 took', (t2 - t1), 'seconds')
parprint('Ground state calc 2 took', (t3 - t2), 'seconds')
parprint('Susceptibility calc 1 took', (t4 - t3), 'seconds')
parprint('Susceptibility calc 2 took', (t5 - t4), 'seconds')
# Check that results are consistent, when structure is simply repeated
# Read results
w11_w, G11_Gc, chiks11_wGG, chi11_wGG = read_component('Al1_chiGG_q1.pckl')
w21_w, G21_Gc, chiks21_wGG, chi21_wGG = read_component('Al2_chiGG_q1.pckl')
w12_w, G12_Gc, chiks12_wGG, chi12_wGG = read_component('Al1_chiGG_q2.pckl')
w22_w, G22_Gc, chiks22_wGG, chi22_wGG = read_component('Al2_chiGG_q2.pckl')
# Check that reciprocal lattice vectors remain as assumed in check below
equal(np.linalg.norm(G11_Gc[0]), 0., 1e-6)
equal(np.linalg.norm(G21_Gc[0]), 0., 1e-6)
equal(np.linalg.norm(G12_Gc[0]), 0., 1e-6)
equal(np.linalg.norm(G22_Gc[1] - np.array([1, 0, 0])), 0., 1e-6)
# Check plasmon peaks remain the same
wpeak11, Ipeak11 = findpeak(w11_w, chi11_wGG[:, 0, 0].imag)
wpeak21, Ipeak21 = findpeak(w21_w, chi21_wGG[:, 0, 0].imag)
equal(wpeak11, wpeak21, 0.02)
equal(Ipeak11, Ipeak21, 1.0)
wpeak12, Ipeak12 = findpeak(w12_w, chi12_wGG[:, 0, 0].imag)
wpeak22, Ipeak22 = findpeak(w22_w, chi22_wGG[:, 1, 1].imag)
equal(wpeak12, wpeak22, 0.02)
equal(Ipeak12, Ipeak22, 1.0)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment