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:
\(V^{\rm{cl}}(\mathbf{r}, t)\) in classical grid
\(V^{\rm{qm}}(\mathbf{r}, t)\) in quantum grid
\(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:
fromaseimport Atoms fromgpaw.fdtd.poisson_fdtdimport QSFDTD fromgpaw.fdtd.polarizable_materialimport (PermittivityPlus, PolarizableMaterial, PolarizableSphere) fromgpaw.tddftimport photoabsorption_spectrum importnumpyasnp # 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:
fromaseimport Atoms fromgpawimport GPAW fromgpaw.fdtd.poisson_fdtdimport FDTDPoissonSolver fromgpaw.fdtd.polarizable_materialimport (PermittivityPlus, PolarizableMaterial, PolarizableSphere) fromgpaw.tddftimport TDDFT, DipoleMomentWriter, photoabsorption_spectrum fromgpaw.inducedfield.inducedfield_tddftimport TDDFTInducedField fromgpaw.inducedfield.inducedfield_fdtdimport FDTDInducedField fromgpaw.mpiimport world importnumpyasnp # 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:
fromgpaw.tddftimport TDDFT fromgpaw.inducedfield.inducedfield_fdtdimport ( FDTDInducedField, calculate_hybrid_induced_field) fromgpaw.inducedfield.inducedfield_tddftimport 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 fromgpaw.mpiimport world importnumpyasnp importmatplotlib.pyplotasplt fromgpaw.inducedfield.inducedfield_baseimport BaseInducedField fromgpaw.tddft.unitsimport aufrequency_to_eV assert world.size == 1, 'This script should be run in serial mode.' # Helper function defdo_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:
