Commit 9940178f authored by Tuomas Rossi's avatar Tuomas Rossi

Calculate density directly from LCAO matrix

parent cb710ef2
......@@ -312,9 +312,10 @@ Induced density
---------------
The density matrix gives access to any other quantities.
For instance, the induced density can be conveniently obtained:
For instance, the induced density can be conveniently obtained
from the density matrix:
.. literalinclude:: lcaotddft_Na8/ksd_ind.py
.. literalinclude:: lcaotddft_Na8/fdm_ind.py
The resulting cube files can be visualized, for example, with
:download:`this script <lcaotddft_Na8/ind_plot.py>` producing
......
import numpy as np
from ase.io import write
from gpaw import GPAW
from gpaw.tddft.units import au_to_eV
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
# Load the objects
calc = GPAW('unocc.gpw', txt=None)
calc.initialize_positions() # Initialize in order to calculate density
dmat = DensityMatrix(calc)
fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
def do(w):
# Select the frequency and the density matrix
rho_MM = fdm.FReDrho_wuMM[w][0]
freq = fdm.freq_w[w]
frequency = freq.freq * au_to_eV
print('Frequency: %.2f eV' % frequency)
print('Folding: %s' % freq.folding)
# Induced density
rho_g = dmat.get_density([rho_MM.imag])
# Save as a cube file
write('ind_%.2f.cube' % frequency, calc.atoms, data=rho_g)
# Calculate dipole moment for reference
dm_v = dmat.density.finegd.calculate_dipole_moment(rho_g, center=True)
absorption = 2 * freq.freq / np.pi * dm_v[0] / au_to_eV * 1e5
print('Total absorption: %.2f eV^-1' % absorption)
do(0)
do(1)
......@@ -6,7 +6,7 @@ def agts(queue):
spec = queue.add('spectrum.py', deps=[td], ncpus=1, walltime=2)
fdm = queue.add('td_fdm_replay.py', deps=[td], ncpus=1, walltime=5)
ksd = queue.add('ksd_init.py', deps=[gs], ncpus=1, walltime=5)
ind = queue.add('ksd_ind.py', deps=[ksd, fdm], ncpus=1, walltime=2)
ind = queue.add('fdm_ind.py', deps=[fdm], ncpus=1, walltime=2)
queue.add('spec_plot.py', deps=[spec], ncpus=1, walltime=2,
creates=['spec.png'])
queue.add('tcm_plot.py', deps=[ksd, fdm, spec], ncpus=1, walltime=2,
......
import numpy as np
from gpaw.utilities import pack
def get_density(rho_MM, wfs, density, density_type='comp', u=0):
rho_G = density.gd.zeros()
kpt = wfs.kpt_u[u]
assert kpt.q == 0
rho_MM = rho_MM.astype(wfs.dtype)
wfs.basis_functions.construct_density(rho_MM, rho_G, kpt.q)
# Uncomment this if you want to add the static part
# rho_G += density.nct_G
if density_type == 'pseudocoarse':
return rho_G
rho_g = density.finegd.zeros()
density.distribute_and_interpolate(rho_G, rho_g)
rho_G = None
if density_type == 'pseudo':
return rho_g
if density_type == 'comp':
D_asp = density.atom_partition.arraydict(density.D_asp.shapes_a)
Q_aL = {}
for a, D_sp in D_asp.items():
P_Mi = wfs.P_aqMi[a][kpt.q]
assert np.max(np.absolute(P_Mi.imag)) == 0
P_Mi = P_Mi.real
assert P_Mi.dtype == float
D_ii = np.dot(np.dot(P_Mi.T.conj(), rho_MM), P_Mi)
D_sp[:] = pack(D_ii)[np.newaxis, :]
Q_aL[a] = np.dot(D_sp.sum(axis=0), wfs.setups[a].Delta_pL)
tmp_g = density.finegd.zeros()
density.ghat.add(tmp_g, Q_aL)
rho_g += tmp_g
return rho_g
raise RuntimeError('Unknown density type: %s' % density_type)
class DensityMatrix(object):
def __init__(self, paw):
self.wfs = paw.wfs
self.density = paw.density
self.using_blacs = self.wfs.ksl.using_blacs
self.tag = None
......@@ -34,3 +76,10 @@ class DensityMatrix(object):
self.rho_uMM.append(rho_MM)
self.tag = tag
return self.rho_uMM
def get_density(self, rho_uMM=None, density_type='comp'):
assert len(self.wfs.kpt_u) == 1, 'K-points not implemented'
u = 0
if rho_uMM is None:
rho_uMM = self.rho_uMM
return get_density(rho_uMM[u], self.wfs, self.density, density_type, u)
......@@ -9,7 +9,6 @@ from gpaw.external import ConstantElectricField
from gpaw.lcaotddft.hamiltonian import KickHamiltonian
from gpaw.lcaotddft.utilities import read_uMM
from gpaw.lcaotddft.utilities import write_uMM
from gpaw.utilities import pack
from gpaw.utilities.tools import tri2full
......@@ -260,10 +259,11 @@ class KohnShamDecomposition(object):
return dm_v
def get_density(self, rho_up, density='comp'):
from gpaw.lcaotddft.densitymatrix import get_density
density_type = density
assert len(rho_up) == 1, 'K-points not implemented'
u = 0
kpt = self.wfs.kpt_u[u]
rho_p = rho_up[u]
C0_nM = self.C0_unM[u]
......@@ -275,43 +275,7 @@ class KohnShamDecomposition(object):
rho_MM = np.dot(C0_iM.T, np.dot(rho_ia, C0_aM.conj()))
rho_MM = 0.5 * (rho_MM + rho_MM.T)
rho_G = self.density.gd.zeros()
assert kpt.q == 0
rho_MM = rho_MM.astype(self.wfs.dtype)
self.wfs.basis_functions.construct_density(rho_MM, rho_G, kpt.q)
# Uncomment this if you want to add the static part
# rho_G += self.density.nct_G
if density_type == 'pseudocoarse':
return rho_G
rho_g = self.density.finegd.zeros()
self.density.distribute_and_interpolate(rho_G, rho_g)
rho_G = None
if density_type == 'pseudo':
return rho_g
if density_type == 'comp':
D_asp = self.density.atom_partition.arraydict(
self.density.D_asp.shapes_a)
Q_aL = {}
for a, D_sp in D_asp.items():
P_Mi = self.wfs.P_aqMi[a][kpt.q]
assert np.max(np.absolute(P_Mi.imag)) == 0
P_Mi = P_Mi.real
assert P_Mi.dtype == float
D_ii = np.dot(np.dot(P_Mi.T.conj(), rho_MM), P_Mi)
D_sp[:] = pack(D_ii)[np.newaxis, :]
Q_aL[a] = np.dot(D_sp.sum(axis=0),
self.wfs.setups[a].Delta_pL)
tmp_g = self.density.finegd.zeros()
self.density.ghat.add(tmp_g, Q_aL)
rho_g += tmp_g
return rho_g
raise RuntimeError('Unknown density type: %s' % density_type)
return get_density(rho_MM, self.wfs, self.density, density_type, u)
def get_contributions_table(self, weight_p, minweight=0.01,
zero_fermilevel=True):
......
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