Hybrid Quantum/Classical Scheme

The basic idea is to separate the calculation into two parts: the first one is the quantum subsystem, which is propagated using Time-propagation TDDFT scheme, and the second one is the classical subsystem that is treated using Quasistatic Finite-Difference Time-Domain method. The subsystems are propagated separately in their own real space grids, but they share a common electrostatic potential.

In the Time-propagation TDDFT part of the calculation the electrostatic potential is known as the Hartree potential \(V^{\rm{qm}}(\mathbf{r}, t)\) and it is solved from the Poisson equation \(\nabla^2 V^{\rm{qm}}(\mathbf{r}, t) = -4\pi\rho^{\rm{qm}}(\mathbf{r}, t)\)

In the Quasistatic Finite-Difference Time-Domain method the electrostatic potential is solved from the Poisson equation as well: \(V^{\rm{cl}}(\mathbf{r}, t) = -4\pi\rho^{\rm{cl}}(\mathbf{r}, t).\)

The hybrid scheme is created by replacing in both schemes the electrostatic (Hartree) potential by a common potential: \(\nabla^2 V^{\rm{tot}}(\mathbf{r}, t) = -4\pi\left[\rho^{\rm{cl}}(\mathbf{r}, t)+\rho^{\rm{qm}}(\mathbf{r}, t)\right].\)

Double grid

The observables of the quantum and classical subsystems are defined in their own grids, which are overlapping but can have different spacings. The following restrictions must hold:

  • The quantum grid must fit completely inside the classical grid

  • The spacing of the classical grid \(h_{\rm{cl}}\) must be equal to \(2^n h_{\rm{qm}}\), where \(h_{\rm{qm}}\) is the spacing of the quantum grid and n is an integer.

When these conditions hold, the potential from one subsystem can be transferred to the other one. The grids are automatically adjusted so that some grid points are common.

Transferring the potential between two grids

  • Transferring the potential from classical subsystem to the quantum grid is performed by interpolating the classical potential to the denser grid of the quantum subsystem. The interpolation only takes place in the small subgrid around the quantum mechanical region.

  • Transferring the potential from quantum subsystem to the classical one is done in another way: instead of the potential itself, it is the quantum mechanical electron density \(\rho^{\rm{qm}}(\mathbf{r}, t)\) that is copied to the coarser classical grid. Its contribution to the total electrostatic potential is then determined by solving the Poisson equation in that grid.

  • Altogether this means that although there is only one potential to be determined \((V^{\rm{tot}}(\mathbf{r}, t))\), three Poisson equations must be solved:

    1. \(V^{\rm{cl}}(\mathbf{r}, t)\) in classical grid

    2. \(V^{\rm{qm}}(\mathbf{r}, t)\) in quantum grid

    3. \(V^{\rm{qm}}(\mathbf{r}, t)\) in classical grid

    When these are ready and \(V^{\rm{cl}}(\mathbf{r}, t)\) is transferred to the quantum grid, \(V^{\rm{tot}}(\mathbf{r}, t)\) is determined in both grids.

Example: photoabsorption of Na2 near gold nanosphere

This example calculates the photoabsorption of \(\text{Na}_2\) molecule in (i) presence and (ii) absence of a gold nanosphere:

from ase import Atoms
from gpaw.fdtd.poisson_fdtd import QSFDTD
from gpaw.fdtd.polarizable_material import (PermittivityPlus,
                                            PolarizableMaterial,
                                            PolarizableSphere)
from gpaw.tddft import photoabsorption_spectrum
import numpy as np

# Nanosphere radius (Angstroms)
radius = 7.40

# Geometry
atom_center = np.array([30., 15., 15.])
sphere_center = np.array([15., 15., 15.])
simulation_cell = np.array([40., 30., 30.])

# Atoms object
atoms = Atoms('Na2', atom_center + np.array([[-1.5, 0.0, 0.0],
                                             [1.5, 0.0, 0.0]]))

# Permittivity of Gold
# J. Chem. Phys. 135, 084121 (2011); https://dx.doi.org/10.1063/1.3626549
eps_gold = PermittivityPlus(data=[[0.2350, 0.1551, 95.62],
                                  [0.4411, 0.1480, -12.55],
                                  [0.7603, 1.946, -40.89],
                                  [1.161, 1.396, 17.22],
                                  [2.946, 1.183, 15.76],
                                  [4.161, 1.964, 36.63],
                                  [5.747, 1.958, 22.55],
                                  [7.912, 1.361, 81.04]])

# 1) Nanosphere only
classical_material = PolarizableMaterial()
classical_material.add_component(PolarizableSphere(center=sphere_center,
                                                   radius=radius,
                                                   permittivity=eps_gold))

qsfdtd = QSFDTD(classical_material=classical_material,
                atoms=None,
                cells=simulation_cell,
                spacings=[2.0, 0.5],
                remove_moments=(1, 1))

energy = qsfdtd.ground_state('gs.gpw', mode='fd', nbands=1, symmetry='off')

qsfdtd.time_propagation('gs.gpw',
                        kick_strength=[0.001, 0.000, 0.000],
                        time_step=10,
                        iterations=1500,
                        dipole_moment_file='dm.dat')

photoabsorption_spectrum('dm.dat', 'spec.1.dat', width=0.15)


# 2) Na2 only (radius=0)
classical_material = PolarizableMaterial()
classical_material.add_component(PolarizableSphere(center=sphere_center,
                                                   radius=0.0,
                                                   permittivity=eps_gold))

qsfdtd = QSFDTD(classical_material=classical_material,
                atoms=atoms,
                cells=(simulation_cell, 4.0),  # vacuum = 4.0 Ang
                spacings=[2.0, 0.5],
                remove_moments=(1, 1))

energy = qsfdtd.ground_state('gs.gpw', mode='fd', nbands=-1, symmetry='off')

qsfdtd.time_propagation('gs.gpw',
                        kick_strength=[0.001, 0.000, 0.000],
                        time_step=10,
                        iterations=1500,
                        dipole_moment_file='dm.dat')

photoabsorption_spectrum('dm.dat', 'spec.2.dat', width=0.15)

# 3) Nanosphere + Na2
classical_material = PolarizableMaterial()
classical_material.add_component(PolarizableSphere(center=sphere_center,
                                                   radius=radius,
                                                   permittivity=eps_gold))

qsfdtd = QSFDTD(classical_material=classical_material,
                atoms=atoms,
                cells=(simulation_cell, 4.0),  # vacuum = 4.0 Ang
                spacings=[2.0, 0.5],
                remove_moments=(1, 1))

energy = qsfdtd.ground_state('gs.gpw', mode='fd', nbands=-1, symmetry='off')

qsfdtd.time_propagation('gs.gpw',
                        kick_strength=[0.001, 0.000, 0.000],
                        time_step=10,
                        iterations=1500,
                        dipole_moment_file='dm.dat')

photoabsorption_spectrum('dm.dat', 'spec.3.dat', width=0.15)

enhanced_absorption

The optical response of the molecule apparently enhances when it is located near the metallic nanoparticle, see Ref. [1] for more examples. The geometry and the distribution of the grid points are shown in the following figure (generated with this script):

geometry

Advanced example: Near field enhancement of hybrid system

In this example we calculate the same hybrid Na2 + gold nanoparticle system as above, but using the advanced syntax instead of the QSFDTD wrapper. This allows us to include InducedField observers in the calculation, see TDDFTInducedField module documentation:

from ase import Atoms
from gpaw import GPAW
from gpaw.fdtd.poisson_fdtd import FDTDPoissonSolver
from gpaw.fdtd.polarizable_material import (PermittivityPlus,
                                            PolarizableMaterial,
                                            PolarizableSphere)
from gpaw.tddft import TDDFT, DipoleMomentWriter, photoabsorption_spectrum
from gpaw.inducedfield.inducedfield_tddft import TDDFTInducedField
from gpaw.inducedfield.inducedfield_fdtd import FDTDInducedField
from gpaw.mpi import world
import numpy as np

# Nanosphere radius (Angstroms)
radius = 7.40

# Geometry
atom_center = np.array([30., 15., 15.])
sphere_center = np.array([15., 15., 15.])
simulation_cell = np.array([40., 30., 30.])

# Atoms object
atoms = Atoms('Na2', atom_center + np.array([[-1.5, 0.0, 0.0],
                                             [1.5, 0.0, 0.0]]))

# Permittivity of Gold
# J. Chem. Phys. 135, 084121 (2011); https://dx.doi.org/10.1063/1.3626549
eps_gold = PermittivityPlus(data=[[0.2350, 0.1551, 95.62],
                                  [0.4411, 0.1480, -12.55],
                                  [0.7603, 1.946, -40.89],
                                  [1.161, 1.396, 17.22],
                                  [2.946, 1.183, 15.76],
                                  [4.161, 1.964, 36.63],
                                  [5.747, 1.958, 22.55],
                                  [7.912, 1.361, 81.04]])

# 3) Nanosphere + Na2
classical_material = PolarizableMaterial()
classical_material.add_component(PolarizableSphere(center=sphere_center,
                                                   radius=radius,
                                                   permittivity=eps_gold))

# Combined Poisson solver
poissonsolver = FDTDPoissonSolver(classical_material=classical_material,
                                  qm_spacing=0.5,
                                  cl_spacing=2.0,
                                  cell=simulation_cell,
                                  communicator=world,
                                  remove_moments=(1, 1))
poissonsolver.set_calculation_mode('iterate')

# Combined system
atoms.set_cell(simulation_cell)
atoms, qm_spacing, gpts = poissonsolver.cut_cell(atoms, vacuum=4.0)

# Initialize GPAW
gs_calc = GPAW(mode='fd',
               gpts=gpts,
               nbands=-1,
               poissonsolver=poissonsolver,
               symmetry={'point_group': False})
atoms.calc = gs_calc

# Ground state
energy = atoms.get_potential_energy()

# Save state
gs_calc.write('gs.gpw', 'all')

# Initialize TDDFT and FDTD
kick = [0.001, 0.000, 0.000]
time_step = 10
iterations = 1500

td_calc = TDDFT('gs.gpw')
DipoleMomentWriter(td_calc, 'dm.dat')

# Attach InducedFields to the calculation
frequencies = [2.05, 2.60]
width = 0.15
cl_ind = FDTDInducedField(paw=td_calc,
                          frequencies=frequencies,
                          width=width)
qm_ind = TDDFTInducedField(paw=td_calc,
                           frequencies=frequencies,
                           width=width)

# Propagate TDDFT and FDTD
td_calc.absorption_kick(kick_strength=kick)
td_calc.propagate(time_step, iterations)

# Save results
td_calc.write('td.gpw', 'all')
cl_ind.write('cl.ind')
qm_ind.write('qm.ind')

photoabsorption_spectrum('dm.dat', 'spec.3.dat', width=width)

The TDDFTInducedField records the quantum part of the calculation and the FDTDInducedField records the classical part. We can calculate the individual and the total induced field by the following script:

from gpaw.tddft import TDDFT
from gpaw.inducedfield.inducedfield_fdtd import (
    FDTDInducedField, calculate_hybrid_induced_field)
from gpaw.inducedfield.inducedfield_tddft import TDDFTInducedField

td_calc = TDDFT('td.gpw')

# Classical subsystem
cl_ind = FDTDInducedField(filename='cl.ind', paw=td_calc)
cl_ind.calculate_induced_field(gridrefinement=2)
cl_ind.write('cl_field.ind', mode='all')

# Quantum subsystem
qm_ind = TDDFTInducedField(filename='qm.ind', paw=td_calc)
qm_ind.calculate_induced_field(gridrefinement=2)
qm_ind.write('qm_field.ind', mode='all')

# Total system, interpolate/extrapolate to a grid with spacing h
tot_ind = calculate_hybrid_induced_field(cl_ind, qm_ind, h=1.0)
tot_ind.write('tot_field.ind', mode='all')

All the InducedField objects can be analyzed in the same way as described in TDDFTInducedField module documentation. Here we show an example script for plotting (run in serial mode, i.e., with one process):

# web-page: cl_field.ind_Ffe.png, qm_field.ind_Ffe.png, tot_field.ind_Ffe.png
from gpaw.mpi import world

import numpy as np
import matplotlib.pyplot as plt

from gpaw.inducedfield.inducedfield_base import BaseInducedField
from gpaw.tddft.units import aufrequency_to_eV

assert world.size == 1, 'This script should be run in serial mode.'


# Helper function
def do_plot(d_g, ng, box, atoms):
    # Take slice of data array
    d_yx = d_g[:, :, ng[2] // 2]
    y = np.linspace(0, box[0], ng[0] + 1)[:-1]
    dy = box[0] / (ng[0] + 1)
    y += dy * 0.5
    ylabel = u'x / Å'
    x = np.linspace(0, box[1], ng[1] + 1)[:-1]
    dx = box[1] / (ng[1] + 1)
    x += dx * 0.5
    xlabel = u'y / Å'

    # Plot
    plt.figure()
    ax = plt.subplot(1, 1, 1)
    X, Y = np.meshgrid(x, y)
    plt.contourf(X, Y, d_yx, 40)
    plt.colorbar()
    for atom in atoms:
        pos = atom.position
        plt.scatter(pos[1], pos[0], s=50, c='k', marker='o')
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.xlim([x[0], x[-1]])
    plt.ylim([y[0], y[-1]])
    ax.set_aspect('equal')


for fname, name in zip(['cl_field.ind', 'qm_field.ind', 'tot_field.ind'],
                       ['Classical subsystem', 'Quantum subsystem',
                        'Total hybrid system']):
    # Read InducedField object
    ind = BaseInducedField(fname, readmode='all')

    # Choose array
    w = 0                                      # Frequency index
    freq = ind.omega_w[w] * aufrequency_to_eV  # Frequency
    box = np.diag(ind.atoms.get_cell())        # Calculation box
    d_g = ind.Ffe_wg[w]                        # Data array
    ng = d_g.shape                             # Size of grid
    atoms = ind.atoms                          # Atoms

    do_plot(d_g, ng, box, atoms)
    plt.title(f'{name}\nField enhancement @ {freq:.2f} eV')
    plt.savefig(fname + '_Ffe.png', bbox_inches='tight')

    # Imaginary part of density
    d_g = ind.Frho_wg[w].imag
    ng = d_g.shape
    do_plot(d_g, ng, box, atoms)
    plt.title('%s\nImaginary part of induced charge density @ %.2f eV' %
              (name, freq))
    plt.savefig(fname + '_Frho.png', bbox_inches='tight')

    # Imaginary part of potential
    d_g = ind.Fphi_wg[w].imag
    ng = d_g.shape
    do_plot(d_g, ng, box, atoms)
    plt.title(f'{name}\nImaginary part of induced potential @ {freq:.2f} eV')
    plt.savefig(fname + '_Fphi.png', bbox_inches='tight')

This produces the following figures for the electric near field:

cl_fe qm_fe tot_fe

References