QMol_disc_basis defines a basis-set discretization built on a Cartesian-grid domain for all basis-based simulations in the QMol-grid package. It also provides a suite of functionalities for data (object) allocation, arithmetic, and various low-level computations that are used throughout the rest of the package.

Class properties

Domain definition

Properties associated with the domain definition can be accessed using the regular dot operator but, outside of the QMol-grid package, can only be edited using the set method.

`basis (v)`

One-electron wave function basis set [ matrix (default []) ]

basis is a numel(xspan)-by-N matrix where N is the basis size (number of vectors): each column in the matrix corresponds to the discretization of a basis vector on the domain grid.
The basis is assumed to be orthonormal. If not (e.g., when the basis is built from a collection of atomic orbitals centered around each atoms in the system) one may use the orthonormalizeBasis to perform the orthonormalization.
Calculations with a non-orthonormal basis set will be erroneous. The class does not check for or enforce the orthonormalization of the basis.

`xspan (x)`

Discretization grid of the simulation domain [ vector (default []) ]

xspan is assumed to be a Cartesian grid (all increasing, equally-spaced, vertices)
QMol_disc does not check whether xspan is indeed a Cartesian grid. Using a non-Cartesian grid will almost certainly produce wrong resutls.

Differential operators

Differential-operator properties are dependent and hold meaningful values only after the QMol_disc_basis object has been initialized. Access differential-operator properties using the regular dot operator. Differential-operator properties cannot be edited with the set method.

`D`

Discretization of the gradient operator $\nabla$ in Fourier space [ vector (default []) ]

At initialization, D is defined as 1i*fourierTool.fftGrid(obj.x).'
D is a column vector with the same number of element as xspan.
Compute the gradient of a function f, defined over xspan and with with vanishing or periodic boundary conditions, with ifft(obj.D.*fft(f)).

`T`

Discretization of the Laplacian operator $-\frac{\Delta }{2}$ in Fourier space [ vector (default []) ]

At initialization, T is defined as -0.5*D.^2
T is a column vector with the same number of element as xspan.
Compute the Laplacian of a function f, defined over xspan and with with vanishing or periodic boundary conditions, with ifft(obj.T.*fft(f)).

`Tv`

Discretization of the velocity-gauge Laplacian operator in Fourier space [ vector (default []) ]

Tv is independently initialized using the setTv method. For a potential vector A, Tv is defined as 0.5*(-1i*D+A).^2.
Tv is a column vector with the same number of element as xspan.
Compute the velocity-gauge Laplacian of a function f, defined over xspan and with with vanishing or periodic boundary consitions, with ifft(obj.Tv.*fft(f)).
Empty Tv indicates that the velocity-gauge laplacian has not been initialized using the setTv method.

`mD`

Matrix representation of the gradient operator D in the basis expansion [ square matrix (default []) ]

mD is a size(basis,2)-by-size(basis,2) matrix with components defined as ${\textrm{mD}}_{k,l} =\int v_k^* (x)\nabla v_l (x)~dx$ .
Compute the gradient of a function f, defined through its basis-expansion population coefficients (in a column vector) as mD*f.

`mT`

Matrix representation of the Laplacian operator T in the basis expansion [ square matrix (default []) ]

mT is a size(basis,2)-by-size(basis,2) matrix with components defined as ${\textrm{mT}}_{k,l} =-\int v_k^* (x)\frac{\Delta }{2}v_l (x)~dx$.
Compute the Laplacian of a function f, defined through its basis-expansion population coefficients (in a column vector) as mT*f.

`mTv`

Matrix representation of the velocity-gauge Laplacian operator Tv in the basis expansion [ square matrix (default []) ]

mTv is independently initialized using the setTv method.
mTv is a size(basis,2)-by-size(basis,2) matrix with components defined as ${\textrm{mTv}}_{k,l} =\int v_k^* (x)\frac{(-i\nabla +A)^2 }{2}v_l (x)~dx$ .
Compute the velocity-gauge Laplacian of a function f, defined through its basis-expansion population coefficients (in a column vector) as mTv*f.
Empty mTv indicates that the velocity-gauge Laplacian has not been initialized using the setTv method.

Initialization status

Internally, QMol_disc_basis objects keep track as to whether they have been initialized after creation. Access the initialization status using the regular dot operator. The initialization status cannot be edited with the set method. All QMol_disc_disc object, irrespective of prior initialization are guarantied to define a meaningful initialization status.

`isInitialized (isInit)`

Initialization status of the QMol_disc_basis object [ false (default) | true ]

isInitialized == false indicates that the object is not (fully) initialized. This is the status after creation, and using any of the set, reset, or clear methods of the class.
isInitialized == true indicates that the object is (fully) initialized. This is the status after running the initialize method.

Miscellaneous

Internally, QMol_disc_basis objects keeps track of some properties that are necessary for some of its methods. Access miscellaneous properties using the regular dot operator. Miscellaneous properties cannot be edited with the set method, they are defined at initialization of the object.

`dx`

Discretization stencil [ scalar (default []) ]

At initialization, dx is defined as xspan(2)-xspan(1).

`basisSize (nV)`

Number of basis vectors discretization [ positive integer (default []) ]

At initialization, basisSize is defined as size(obj.basis,2).
basisSize is typically used for data preallocation (thus the name).
For developers: Note that this definition is different from the parent QMol_disc class where basisSize = numel(xspan). This is deliberate such that in both cases one basisSize specifies the size of discretization objects like orbitals, etc.

`QM`

Parent QMol-grid theory object [ handle-class object (default []) ]

Knowledge of the parent QM object is required for running some of the QMol_disc_basis methods -- such methods are identified below.
For developers: For practical reasons, QM is editable by QMol_DFT classes.

`isBasis`

Flag indicating that QMol_disc_basis objects are indeed a basis-set discretization [ true ]

At run time, isBasis can be used to discriminate QMol_disc_basis objects from the parent QMol_disc, which model grid-based discretizations.
For class development purposes, isBasis is technically a static method. Practically, though, it can virtually almost always be treated as a constant property.

`isBasisPol`

Flag indicating that QMol_disc_basis objects do not describe a spin-polarized basis set [ false ]

This means that up- and down-spin components are described using the same basis set of one-electron wave functions
For class development purposes, isBasisPol is technically a static method. Practically, though, it can virtually almost always be treated as a constant property.

Generic class methods

Creation

default constructor

Create a QMol_disc_basis object with all properties set to their default values

obj = QMol_disc_basis;

Create a QMol_disc_basis object with name property set to the specified value. Several name-value pairs can be specified consecutively. Suitable name is any of the domain definition properties and is case insensitive.

obj = QMol_disc_basis(name,value);
obj = QMol_disc_basis(name1,value1,name2,value2,___);

Changing class properties

After creation, the domain definition properties of a QMol_disc_basis object can be edited.

`set`

Update the name property of a QMol_disc_basis object to the specified value. Several name-value pairs can be specified consecutively. Suitable name are specified in the domain definition properties and are case insensitive.

obj.set(name1,value1);
obj.set(name1,value1,name2,value2,___);

The set method also reset the class.
After running, the set property updates the isInitialized flag to a false value.

`reset`

Reset the object and update the isInitialized flag to a false.

obj.reset;

The reset method clears the dependent properties D, T, Tv, mD, mT, mTv, dx, basisSize, and QM.
After running, the reset property updates the isInitialized flag to a false value.

`clear`

Clear all the domain definition properties

obj.clear;

Clear a specific set of the domain definition properties. Suitable name are any of the domain definition properties and are case insensitive.

obj.clear(name1,name2,___);

The clear method also reset the class.
The clear method can be used to delete specific properties before saving a QMol_disc object.

`orthonormalizeBasis`

Orthonormalize the basis vectors

obj.orthonormalizeBasis;

Performs the basis-set orthonormalization using the default (Gram-Schmidt, see next) algorithm such that sum(conj(basis(:,k)).*basis(:,l))*dx = 1 for k == l, and 0 otherwise.
If the object was initialized (isInitialized == true ), after the orthonormalization it is reinitialized with the new basis set.

Use a Gram-Schmidt orthonormalization scheme with either

orthonormalizeBasis('Gram–Schmidt');
orthonormalizeBasis('Gram–Schmidt',[]);

Basis vectors are orthonormalized in the (row) order of the basis matrix.

Specify the order in which the basis would be orthonormalized with

orthonormalizeBasis('Gram–Schmidt',ind);

ind is a size(basis,2) vector containing the ordered indexes in which the orthonormalization should be performed (e.g., ind = [4 3 2 1 5 ___]).

Use a balanced orthonormalization scheme

orthonormalizeBasis('overlap');

The orthonormalization is performed by diagonalizing the overlap matrix between the basis vectors, therefore treating them on an equal footing.

Initializing `QMol_disc_basis` objects

To set dependent class properties that are required for most of the class properties, QMol_disc_basis objects must be initialized.

`initialize`

Initialize a QMol_disc_disc object without attaching a QMol-grid theory object

obj.initialize;
obj.initialize([]);

This initializes the dependent properties D, T, mD, mT, dx, and basisSize, and sets the isInitialized flag to true.
The QM property is left empty.
Methods requiring a QMol-grid theory object specified in QM will not work, most likely causing an error.

Initialize a QMol_disc_basis object with attaching a QMol-grid theory object

obj.initialize(QM);

This performs the same initialization as previously plus copying the input QM object into the QM property.
Methods requiring a QMol-grid theory object specified in QM will not work provided QM is of a compatible type.
Properties that require a QMol-grid theory object specified in QM and the compatibility criteria are specified in each property description below.

`setTv`

After the QMol_disc_basis object has been initialized, set the velocity-gauge Laplacian operator Tv and matrix mTv with

obj.setTv(A);

The input scalar A specifies the value of the vector potential for the velocity-gauge Laplacian.

Clear the velocity-gauge Laplacian operator Tv matrix mTv with either of

obj.setTv;
obj.setTv([]);

Run-time documentation

`getMemoryProfile`

Get an estimate of the memory held by a QMol_disc_basis object with either

mem = obj.getMemoryProfile;
mem = obj.getMemoryProfile(false);

The object must be properly initialized.
The estimate only includes the discretization of member components over the domain grid and ignores other (small) properties.
The output mem is the estimated size in bytes.

Additionally display the detail of the memory footprint with

mem = obj.getMemoryProfile(true);

Comparing domain discretizations

The class overloads the == and ~= operators to facilitate the comparison between QMol_disc_basis objects

`eq (==)`

Test if two QMol_disc_basis objects describe the same domain discretization

obj1 == obj2

a true result is defined as (i) both domain discretization grids xspan are identical within eqTol tolerance; (ii) both are QMol_disc_basis objects; (iii) both have the same basis-set vectors (in the same order) within eqTol tolerance.
eqTol = 1e-10 is a tolerance parameter to account for possible roundoff mismatch between two domain-discretization objects.
Any mismatch in the size of the two domain grids xspan or size of the basis vectors basis returns a false result.
The equality test ignores any QMol-grid theory object specified in QM.

`ne (~=)`

Test if two QMol_disc_basis objects describe different domain discretizations

obj1 ~= obj2

This is equivalent to ~(obj1 == obj2) .

DFT class methods

Component names (static)

To facilitate housekeeping of a specific domain discretization compared with overloading classes for a different flavors of discretization, QMol_disc_basis defines static methods returning the names of the classes that describe one-body densities, and Kohn-Sham orbitals and potentials

`DFT_classDensity (static)`

Get the name of the class describing one-body densities

obj.DFT_classDensity;

Returns 'QMol_DFT_density'

`DFT_classOrbital (static)`

Get the name of the class describing Kohn-Sham orbitals

obj.DFT_classOrbital;

Returns 'QMol_DFT_orbital_basis'

`DFT_classPotential (static)`

Get the name of the class describing Kohn-Sham potentials

obj.DFT_classPotential;

Returns 'QMol_DFT_Vks_basis'

`DFT_classPotentialGradient (static)`

Get the name of the class describing Kohn-Sham potential gradients

obj.DFT_classPotentialGradient;

The most appropriate implementation of the Kohn-Sham potential derivative for basis-set discretization is not yet clear. Until then, the method returns 'N/A'

Data allocation

Note: all DFT data/object allocation require the QMol_disc_basis object to have been initialized with a DFT-theory object (QMol_DFT_spinRes or QMol_DFT_spinPol object).

`DFT_allocateDensity`

Allocate a new one-body density object with either

rho = obj.DFT_allocateDensity;
rho = obj.DFT_allocateDensity([]);
rho = obj.DFT_allocateDensity([],false);

Only the density object is created and no memory is preallocated for the discretization of the one-body density itself.
The allocation also initializes the QMol_DFT_density object, i.e., the isInitialized flag of the output rho is set to true.
Note that the one-body density is discretized on the domain grid and not using the basis vectors (the same class is used for one-body density objects as with grid-based models).

Allocate a new one-body density object and force the discretization of the one-body density to be preallocated with zero values

rho = obj.DFT_allocateDensity([],true);

Reallocate an existing one-body density object

obj.DFT_allocateDensity(rho);
obj.DFT_allocateDensity(rho,false);
obj.DFT_allocateDensity(rho,true);

This will first erase all data in the input QMol_DFT_density object rho and then perform the (re)allocation.

`DFT_allocateOrbital`

Allocate a new Kohn-Sham-orbital object with either

KSO = obj.DFT_allocateOrbital(N);
KSO = obj.DFT_allocateOrbital(N,[]);

For spin restricted DFT models, N is a positive integer specifying the number of orbitals to be allocated.
For spin polarized DFT models, N is a two vector of positive integers specifying the number of orbitals to be allocated in the up- and down-spin channels, respectively.
In all cases, the relevant properties in the output QMol_DFT_orbital_basis object are allocated to zero.
The allocation also initializes the QMol_DFT_orbital_basis object, i.e., the isInitialized flag of the output KSO is set to true.

Reallocate and existing Kohn-Sham-orbital object

obj.DFT_allocateOrbital(N,KSO);

If the input QMol_DFT_orbital_basis object KSO has a different spin polarization than QM, or the size of the member orbital properties does not match the domain discretization, all content in KSO is wiped out and a new, zero-valued, orbital object is allocated.
Otherwise, the content of the orbital properties is resized to match the requested number of orbitals: The final orbitals are removed for smaller N or zero-value orbitals are added at the end for larger N.
Orbital resizing is typically used to restart (TD)DFT computations from a previous result instead of scratch.
To force the reallocation of zero-valued orbitals irrespective of matching spin polarization of domain size, first clear the orbital object with KSO.clear before obj.DFT_allocateOrbital(N,KSO).
The allocation also initializes the QMol_DFT_orbital_basis object, i.e., the isInitialized flag of the input KSO is set to true.

Allocate with random values (instead of the default zeros) with

KSO = obj.DFT_allocateOrbital(N,[],randStr);        % Allocate new orbitals
KSO = obj.DFT_allocateOrbital(N,[],randStr,false);
obj.DFT_allocateOrbital(N,KSO,randStr);             % Reallocate existing ones
obj.DFT_allocateOrbital(N,KSO,randStr,false);

randStr is a random number stream from which the random values are taken.
The (random) orbitals are complex valued and normalized (but not orthogonal).

Allocate with random real values with

KSO = obj.DFT_allocateOrbital(N,[],randStr,true);   % Allocate new orbitals
obj.DFT_allocateOrbital(N,KSO,randStr,true);        % Reallocate existing ones

`DFT_allocatePotential`

Allocate a new Kohn-Sham potential object with either

V = obj.DFT_allocatePotential;
V = obj.DFT_allocatePotential([]);

Only the potential object is created and no memory is preallocated for the discretization of the Kohn-Sham potential itself.
The allocation does not initialize the QMol_DFT_Vks_basis object, i.e., the isInitialized flag of the output V is left to false (unlike for the density and orbitals).

Reallocate an existing Kohn-Sham potential object

obj.DFT_allocatePotential(V);

This will first erase all data in the input QMol_DFT_Vks_basis object V and then perform the (re)allocation.

`DFT_allocatePotentialGradient`

The most appropriate implementation of the Kohn-sham potential gradient for basis-set discretization is not yet clear. Until then, accessing the DFT_allocatePotentialGradient method throws an error.

Kohn-Sham orbitals

QMol_disc_basis provides facilities to navigate between grid-based and basis-set based discretizations of the Kohn-Sham orbitals

`DFT_projectOrbital`

Convert a grid-discretization Kohn-Sham orbital object QMol_DFT_orbital into its corresponding basis discretization

proj = obj.DFT_projectOrbital(KSO);

The input KSO is a QMol_DFT_orbital object and the output proj is a QMol_DFT_orbital_basis
If KSO is initialized KSO.isInitialized == true , so does the output proj
The projected orbitals in proj may not be orthonormalized (if discretized orbitals have components outside of the hull of the projection basis).

Store the projected/converted Kohn-Sham orbitals in an existing QMol_DFT_orbital_basis object

obj.DFT_projectOrbital(KSO,proj);

Any content in the QMol_DFT_orbital_basis object proj is the erased before the conversion/projection

`DFT_reconstructOrbital`

Convert a basis-set discretization object into its corresponding grid-based object QMol_DFT_orbital

KSO = obj.DFT_reconstructOrbital(proj);

The input proj is a QMol_DFT_orbital_basis and the output KSO is a QMol_DFT_orbital object
If proj is initialized proj.isInitialized == true , so does the output KSO

Store the reconstructed Kohn-Sham orbitals in an existing QMol_DFT_orbital object

obj.DFT_reconstructOrbital(proj,KSO);

Any content in the QMol_DFT_orbital object KSO is the erased before the reconstruction

Miscellaneous

`DFT_dist`

Compute the distance between two one-body densities

d = obj.DFT_dist(rho1,rho2);

For spin restricted DFT models, the output scalar d contains the numerical evaluation of $\sqrt{\int (\rho_1 (x)-\rho_2 (x))^2 ~dx}$
For spin polarized DFT models, the output two vector d contains the numerical evaluation of $\sqrt{\int (\rho_1^{\uparrow } (x)-\rho_2^{\uparrow } (x))^2 ~dx}$ and $\sqrt{\int (\rho_1^{\downarrow } (x)-\rho_2^{\downarrow } (x))^2 ~dx}$ , respectively
Both one-body density objects must be allocated and initialized beforehand, and have matching spin polarized or restricted nature.

`DFT_mix`

Linearly mix pairs of spin-restricted one-body densities

rho = obj.DFT_mix(rho,alpha1,rho1,alpha2,rho2,___);

Performs the linear mixing of densities $\rho =\alpha_1 \rho_1 +\alpha_2 \rho_2 +\cdots$
Each of the rho is a spin-restricted one-body density object and each of the alpha is a scalar
Note that this is an in-place function and the content of the first input one-body density object is changed by DFT_mix
obj.DFT_mix(rho1,alpha1,rho1,alpha2,rho2,___), with the one-body density object repeated in the first mixing pair, formally performs in place mixing $\rho_1 \gets \alpha_1 \rho_1 +\alpha_2 \rho_2 +\cdots$

Linearly mix pairs of spin-polarized one-body densities

rho = obj.DFT_mix(rho,[a11 a12],rho1,[a21 a22],rho2,___);

Performs the linear mixings of densities $\rho^{\uparrow } =a_{11} \rho_1^{\uparrow } +a_{21} \rho_2^{\uparrow } +\cdots$ and $\rho^{\downarrow } =a_{21} \rho_1^{\downarrow } +a_{22} \rho_2^{\downarrow } +\cdots$
Each of the rho is a spin-polarized one-body density object and each of the a is a scalar
Note that this is an in-place function and the content of the first input one-body density object is changed by DFT_mix
obj.DFT_mix(rho,[a11 a12],rho1,[a21 a22],rho2,___), with the one-body density object repeated in the first mixing pair, formally performs in place mixing

Linearly mix pairs of Kohn-Sham potentials

Vks = obj.DFT_mix(Vks,alpha1,Vks1,alpha2,Vks2,___);

Performs the linear mixing of potentials ${\mathcal{V}}_{{\mathrm{K}\mathrm{S}}} =\alpha_1 {\mathcal{V}}_{{\mathrm{K}\mathrm{S},1}} +\alpha_2 {\mathcal{V}}_{{\mathrm{K}\mathrm{S},2}} +\cdots$
Each of the Vks is a Kohn-Sham potential object
For spin-restricted DFT models, each of the alpha is a scalar
For spin-polarized DFT models, each of the alpha is a two vector indicating the mixing coefficients for the up- and down-spin channels
In-place mixing can be performed by repeating the output Kohn-Sham potential object in the first mixing pair obj.DFT_mix(Vks1,alpha1,Vks1,alpha2,Vks2,___)
Note that only the explicit part of the potential is mixed: any implicit components of input ${\mathcal{V}}_{{\mathrm{K}\mathrm{S}}}$ are left unchanged in the output while all implicit components of ${\mathcal{V}}_{{\mathrm{K}\mathrm{S},1}} ,{\mathcal{V}}_{{\mathrm{K}\mathrm{S},2}} , \cdots$ are ignored.
Note that this is an in-place function and the content of the first input Kohn-Sham potential object is changed by DFT_mix

DFT_mix does not perform any density/potential allocation or initialization; all input objects must be created and initialized beforehand. DFT_mix is typically used by self-consistent-field solvers in ground-state DFT computations

`DFT_normalizeOrbital`

Normalize Kohn-Sham orbitals to

p = obj.DFT_normalizeOrbital(p);

The input p is a matrix whose columns contain the Kohn-Sham orbitals to be normalized
The normalization is such that $\int |\phi_k (x)|^2 ~dx=1$ for all the orbitals $\phi_k$ in p (via their basis-expansion coefficients)
Note that this is an in-place function which avoids creating a copy of the matrix p when performing the normalization

`DFT_randomOrbital`

Generate a random Kohn-Sham orbital from random number stream

p = obj.DFT_randomOrbital(randStr);

The input random number stream randStr enables reproducibility across runs
The output vector p is a single normalized Kohn-Sham orbital

`DFT_SCF_AndersonMixCoeff`

Compute the mixing coefficient for self-consistent-field computations with an Anderson mixing scheme from one-body densities

b = obj.DFT_SCF_AndersonMixCoeff(rho1i,rho2i,rho1o,rho2o);

For spin restricted DFT models, the output scalar b contains the numerical evaluation of the mixing coefficient for the Anderson mixing scheme
For spin polarized DFT models, the output two vector contains the numerical evaluation of the mixing coefficient for the Anderson mixing scheme
All one-body density objects must be allocated and initialized beforehand, and have matching spin polarized or restricted nature.

Compute the mixing coefficient for self-consistent-field computations with an Anderson mixing scheme from Kohn-Sham potentials

b = obj.DFT_SCF_AndersonMixCoeff(Vks1i,Vks2i,Vks1o,Vks2o);

For spin restricted DFT models, the output scalar b contains the numerical evaluation of the mixing coefficient for the Anderson mixing scheme
For spin polarized DFT models, the output two vector contains the numerical evaluation of the mixing coefficient for the Anderson mixing scheme
Note that only the explicit part of the potential is used to define the mixing coefficient and any implicit component is ignored.
All Kohn-Sham potentials must be allocated and initialized beforehand, and have matching spin polarized or restricted nature.

DFT_SCF_AndersonMixCoeff is used in the QMol_DFT_SCF_Anderson class and is not expected to have any use for end users.

`DFT_sizeOrbital`

Get the size of a Kohn-Sham orbital

s = obj.DFT_sizeOrbital;

This is defined as s = [size(obj.v,2) 1]
This is typically used for data preallocation

Operators and functionals

`DFT_energyKinetic`

Compute the DFT kinetic energy for specifies occupation coefficients and Kohn-Sham orbitals

E = obj.DFT_energyKinetic(occ,KSO);

For spin-restricted models, the output positive real E corresponds to the numerical evaluations of the DFT kinetic energy defined as $E=\sum_k n_k \langle \phi_k |\hat{\mathcal{T}} |\phi_k \rangle$ . The population coefficients $n_k$ are specified in the input vector occ and the Kohn-Sham orbitals $|\phi_k \rangle$ in the QMol_DFT_orbital_basis object KSO.
For spin-polarized models, the ouput positive vector [Eup Edown] corresponds to the numerical evaluation of the DFT kinetic energy in the up- and down-spin channels, respectively. The population coefficients for each spin channel are grouped in a cell occ = {occUp,occDown}.
DFT_energyKinetic automatically detect whether to use the length- or velocity-gauge kinetic operator $\hat{\mathcal{T}}$ to compute the energy.

`DFT_energyOrbital`

Compute the energy of Kohn-Sham orbitals for a specific potential

E = obj.DFT_energyOrbital(Vks,KSO);

For an orbital $\left|\phi \right\rangle$ , the orbital energy is defined as the real $E=\langle \phi |{\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}} |\phi \rangle$ where ${\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}}$ is the total DFT Hamiltonian operator.
For spin-restricted models, E is a N-by-1 vector containing the energies of each of the N Kohn-Sham orbitals in the QMol_DFT_orbital_basis object KSO.
For spin-polarized models, E is a cell containing a Nup-by-1 and Ndown-by-1 vectors respectively with energies of each of the Nup and Ndown Kohn-Sham orbitals in the up- and down-spin channels of the QMol_DFT_orbital_basis object KSO.
The input Vks is a QMol_DFT_Vks_basis object describing the Kohn-Sham potential to be used in the Hamiltonian operator

Additionally compute the errors for the Kohn-Sham orbitals

[E,DE] = obj.DFT_energyOrbital(Vks,KSO);

For an orbital $\left|\phi \right\rangle$ , the error is defined as the positive real $\left\|({\hat{\mathcal{H}} }_{{\mathrm{D}\mathrm{F}\mathrm{T}}} -E)|\phi \rangle \right\|$ with the orbital energy $E$ defined as above.
For spin-restricted or polarized models, the additional output DE is a vector or cell of vectors with the errors matching each orbital energy.

`DFT_operatorHamiltonian`

For spin-restricted models, apply the total, kinetic plus potential, DFT Hamiltonian operator to a wave function

Hp = obj.DFT_operatorHamiltonian(Vks,p);

p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization of the orbital to which the potential should be applied and its result. Vks is a QMol_DFT_Vks_basis object describing the Kohn-Sham potential.
DFT_operatorHamiltonian automatically detect whether to use the length- or velocity-gauge kinetic operator.

For spin-polarized models, apply the total, kinetic plus potential, DFT Hamiltonian operator to a wave function respectively with

Hp = obj.DFT_operatorHamiltonian(Vks,p,true);   % spin up
Hp = obj.DFT_operatorHamiltonian(Vks,p,false);  % spin down

p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization of the orbital to which the potential should be applied and its result. Vks is a QMol_DFT_Vks_basis object describing the Kohn-Sham potential.
DFT_operatorHamiltonian automatically detect whether to use the length- or velocity-gauge kinetic operator.

`DFT_operatorKinetic`

Apply the kinetic operator to a wave function

Hp = obj.DFT_operatorKinetic(p);

In the velocity gauge, after setting the Laplacian with setTv, numerically computes $\mathcal{P}\frac{(-i\nabla +A)^2 }{2}\left|\phi \right\rangle$ on the input wave function p, where $\mathcal{P}$ is the orthonormal projection onto the hull of hte basis set.
In the length gauge, numerically computes $-\frac{\Delta }{2}\left|\phi \right\rangle$ on the input wave function p
p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization of the orbital to which the potential should be applied and its result.

`DFT_operatorPotential`

For spin-restricted models, apply the Kohn-Sham potential operator to a wave function

Hp = obj.DFT_operatorPotential(Vks,p);

This computes the action of the entire, including both explicit and implicit components, Kohn-Sham potential operator described by the Vks object on the input wave function p.
p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization of the orbital to which the potential should be applied and its result. Vks is a QMol_DFT_Vks_basis object
This method interface is provided to completeness, instead Hp = Vks.applyPotential(p) is recommended

For spin-polarized models, apply the up- and down-spin Kohn-Sham potential operators to a wave function respectively with

Hp = obj.DFT_operatorPotential(Vks,p,true);   % spin up
Hp = obj.DFT_operatorPotential(Vks,p,false);  % spin down

This computes the action of the entire -- including both explicit and implicit components -- Kohn-Sham potential operator described by the Vks object on the input wave function p.
p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization of the orbital to which the potential should be applied and its result. Vks is a QMol_DFT_Vks_basis object
This method interface is provided to completeness, instead Hp = Vks.applyPotential(p,true) and Hp = Vks.applyPotential(p,false) are recommended

Schrödinger-equation class methods

Component names (static)

To facilitate housekeeping of a specific domain discretization compared with overloading classes for a different flavors of discretization, QMol_disc defines static methods returning the names of the classes that describe wave functions and potentials.

`SE_classWaveFunction (static)`

Get the name of the class describing wave functions

obj.SE_classWaveFunction;

Returns 'QMol_SE_wfcn_basis'

`SE_classPotential (static)`

Get the name of the class describing potentials

obj.DFT_classPotential;

Returns 'QMol_SE_V'

Data allocation

Note: all Schrödinger-equation data/object allocation require the QMol_disc object to have been initialized with a Schrödinger-equation object (QMol_SE object).

`SE_allocateWaveFunction`

Allocate a new wave function object with either

wfcn = obj.SE_allocateWaveFunction(N);
wfcn = obj.SE_allocateWaveFunction(N,[]);

N is a positive integer specifying the number of wave functions to be allocated.
The relevant properties in the output QMol_SE_wfcn_basis object are allocated to zero.
The allocation also initializes the QMol_SE_wfcn_basis object, i.e., the isInitialized flag of the output wfcn is set to true.

Reallocate and existing wave function object

wfcn = obj.SE_allocateWaveFunction(N,wfcn);

If the size of the member wave functions property does not match the domain discretization, all content of the wave function in the input QMol_SE_wfcn_basis object wfcn is wiped out and a new, zero-valued, wave function object is allocated.
Otherwise, the content of the wave function property is resized to match the requested number of orbitals: The final wave functions are removed for smaller N or zero-value wave functions are added at the end for larger N.
Wave function resizing is typically used to restart (TD)SE computations from a previous result instead of from scratch.
To force the reallocation of zero-valued wave functions irrespective of matching domain size, first clear the wave function object with wfcn.clear before obj.SE_allocateWaveFunction(N,wfcn).
The allocation also initializes the QMol_SE_wfcn_basis object, i.e., the isInitialized flag of the input wfcn is set to true.

Allocate with random values (instead of the default zeros) with

wfcn = obj.SE_allocateWaveFunction(N,[],randStr);        % Allocate new wave functions
wfcn = obj.SE_allocateWaveFunction(N,[],randStr,false);
obj.SE_allocateWaveFunction(N,wfcn,randStr);             % Reallocate existing ones
obj.SE_allocateWaveFunction(N,wfcn,randStr,false);

randStr is a random number stream from which the random values are taken.
The (random) wave functions are complex valued and normalized (but not orthogonal).

Allocate with random real values with

wfcn = obj.SE_allocateWaveFunction(N,[],randStr,true);   % Allocate new wave functions
obj.SE_allocateWaveFunction(N,wfcn,randStr,true);        % Reallocate existing ones

Wave functions

QMol_disc_basis provides facilities to navigate between grid-based and basis-set based discretizations of the wave functions

`SE_projectWaveFunction`

Convert a grid-discretization wave function object QMol_SE_wfcn into its corresponding basis discretization

proj = obj.SE_projectWaveFunction(wfcn);

The input wfcn is a QMol_SE_wfcn object and the output proj is a QMol_SE_wfcn_basis
If wfcn is initialized wfcn.isInitialized == true , so does the output proj
The projected wave functions in proj may not be orthonormal (if discretized orbitals have components outside of the hull of the projection basis).

Store the projected/converted wave functions in an existing QMol_SE_wfcn_basis object

obj.SE_projectWaveFunction(wfcn,proj);

Any content in the QMol_SE_wfcn_basis object proj is the erased before the conversion/projection

`SE_reconstructWaveFunction`

Convert a basis-set discretization object into its corresponding grid-based object QMol_SE_wfcn

wfcn = obj.SE_reconstructWaveFunction(proj);

The input proj is a QMol_SE_wfcn_basis and the output wfcn is a QMol_SE_wfcn object
If proj is initialized proj.isInitialized == true , so does the output wfcn

Store the reconstructed wave functions in an existing QMol_SE_wfcn object

obj.SE_reconstructWaveFunction(proj,wfcn);

Any content in the QMol_SE_wfcn object wfcn is the erased before the reconstruction

Miscellaneous

`SE_normalizeWaveFunction`

Normalize wave functions

p = obj.SE_normalizeWaveFunction(p);

The input p is a matrix whose columns contain the wave functions to be normalized
The normalization is such that $\int |\psi_k (x)|^2 ~dx=1$ for all the wave functions $\psi_k$ in p
Note that this is an in-place function which avoids creating a copy of the matrix p when performing the normalization

`SE_randomWaveFunction`

Generate a random wave function from random number stream

p = obj.SE_randomWaveFunction(randStr);

The input random number stream randStr enables reproducibility across runs
The output vector p is a single normalized wave function

`SE_sizeWaveFunction`

Get the size of a wave function

s = obj.SE_sizeWaveFunction;

This is defined as s = [size(obj.v,2) 1]
This is typically used for data preallocation

Operators and functionals

`SE_energyKinetic`

Compute the Schrödinger-equation kinetic energy for the input wave functions

E = obj.SE_energyKinetic(occ,wfcn);

The output positive real E corresponds to the numerical evaluations of the Schrödinger-equation kinetic energy, summed over all the wave functions, defined as $E=\sum_k \langle \psi_k |\hat{\mathcal{T}} |\psi_k \rangle$ .
SE_energyKinetic automatically detect whether to use the length- or velocity-gauge kinetic operator $\hat{\mathcal{T}}$ to compute the energy.

`SE_energyWaveFunction`

Compute the energy of wave functions for a specific potential

E = obj.SE_energyWaveFunction(V,wfcn);

For a wave function $\left|\psi \right\rangle$ , the wave function energy is defined as the real $E=\langle \psi |\hat{\mathcal{H}} |\psi \rangle$ where $\hat{\mathcal{H}}$ is the total Schrödinger-equation Hamiltonian operator.
E is a N-by-1 vector containing the energies of each of the N wave functions orbitals in the QMol_SE_wfcn_basis object wfcn.
The input V is a QMol_SE_V object describing the potential to be used in the Hamiltonian operator

Additionally compute the errors for the wave functions

[E,DE] = obj.SE_energyWaveFunction(Vks,wfcn);

For a wave function $\left|\psi \right\rangle$ , the error is defined as the positive real $\left\|(\hat{\mathcal{H}} -E)|\psi \rangle \right\|$ with the wave function energy $E$ defined as above.

`SE_operatorHamiltonian`

Apply the total, kinetic plus potential, Schrödinger Hamiltonian operator to a wave function with either

Hp = obj.SE_operatorHamiltonian(Vks,p);

p and Hp are both numel(xspan)-by-1 vectors respectively containing the discretization of the wave function to which the potential should be applied and its result. V is a QMol_SE_V object describing the potential.
SE_operatorHamiltonian automatically detect whether to use the length- or velocity-gauge kinetic operator.

`SE_operatorKinetic`

Apply the kinetic operator to a wave function with either

Hp = obj.SE_operatorKinetic(p);

In the velocity gauge, after setting the Laplacian with setTv, numerically computes $\frac{(-i\nabla +A)^2 }{2}\left|\psi \right\rangle$ on the input wave function p
In the length gauge, numerically computes $-\frac{\Delta }{2}\left|\psi \right\rangle$ on the input wave function p
p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization, on the basis, of the wave functions to which the potential should be applied and its result.

`SE_operatorPotential`

Apply the potential operator to a wave function with either

Hp = obj.SE_operatorPotential(V,p);

This computes the action of the potential operator described by the V object on the input wave function p.
p and Hp are both size(basis,2)-by-1 vectors respectively containing the discretization, on the basis, of the wave function to which the potential should be applied and its result. V is a QMol_SE_V object

Examples

Define the domain grid and set of vectors from which to build the basis set

% Domain grid
x   =   (-15:.1:20).';

% vector set
v   =  [exp(-(x-2).^2  ), (x-2).*exp(-(x-2).^2/2), ...
        exp(-(x+3).^2/2), (x+3).*exp(-(x+3).^2/4)];

Create a discretization object and orthonormalize the basis set

% Domain discretization
disc = QMol_disc_basis('x',x.','basis',v);
disc.orthonormalizeBasis('overlap');
disc.initialize([]);

% Plot basis set
figure; hold on
plot(disc.x,disc.basis(:,1),'-','LineWidth',2,'DisplayName','v_1')
plot(disc.x,disc.basis(:,2),'-','LineWidth',2,'DisplayName','v_2')
plot(disc.x,disc.basis(:,3),'-','LineWidth',2,'DisplayName','v_3')
plot(disc.x,disc.basis(:,4),'-','LineWidth',2,'DisplayName','v_4')
xlabel('position (a.u.)'); xlim(disc.x([1 end]));
ylabel('basis set')
legend show

Display the estimated memory footprint for the object

disc.getMemoryProfile(true);

  * Domain discretization                                        
    > domain                                                        2.7 KB
    > gradient                                                      5.5 KB
    > Laplacian                                                     2.7 KB
    > Laplacian (velocity gauge)                                    5.5 KB
    > basis                                                        11.0 KB
    > gradient matrix                                               128  B
    > Laplacian matrix                                              128  B
    > Laplacian matrix (velocity gauge)                             256  B

Compute the projected ( $\mathcal{P}$ -- on the basis set) gradient and Laplacian of a function

% Basis-set expansion
f   =   [1; -.5; -.5; 1];
df  =   disc.mD*f;
ddf =   disc.mT*f;

% Reconstruct the functions
F   =     f(1)*disc.v(:,1)+  f(2)*disc.v(:,2)+ ...
          f(3)*disc.v(:,3)+  f(4)*disc.v(:,4);
DF  =    df(1)*disc.v(:,1)+ df(2)*disc.v(:,2)+ ...
         df(3)*disc.v(:,3)+ df(4)*disc.v(:,4);
DDF =   ddf(1)*disc.v(:,1)+ddf(2)*disc.v(:,2)+ ...
        ddf(3)*disc.v(:,3)+ddf(4)*disc.v(:,4);

Plot the three functions

figure; hold on
plot(disc.x,F,'-','LineWidth',2,'DisplayName','f')
plot(disc.x,DF,'-','LineWidth',2','DisplayName','\nabla{f}')
plot(disc.x,DDF,'-','LineWidth',2','DisplayName','-\Delta{f}/2')
xlabel('position (a.u.)'); xlim(disc.x([1 end]));
ylabel('function/gradient/Laplacian')
legend show

Finally, we can check the accuracy of the projected derivatives

% Compute derivatives
Df  =  ifft(disc.D.*fft(F));
Ddf =  ifft(disc.T.*fft(F));

% Project derivatives on basis
for k = 1:disc.basisSize
    df(k)   =    df(k) - sum(Df .*disc.v(:,k))*disc.dx;
    ddf(k)  =   ddf(k) - sum(Ddf.*disc.v(:,k))*disc.dx;
end

% Display errors
fprintf('Projected gradient  max. error = %g\n',max(abs(df )));
fprintf('Projected Laplacian max. error = %g\n',max(abs(ddf)));

yielding

Projected gradient  max. error = 4.44089e-16
Projected Laplacian max. error = 1.11022e-16

Test suite

Run the test suite for the QMol_disc_basis class in normal or summary mode respectively with

QMol_test.test('disc_basis');
QMol_test.test('-summary','disc_basis');

For developers

Class definition

QMol_disc_basis defines a concrete handle class overloading QMol_disc (and thus QMol_suite).

QMol_disc_basis is a cornerstone component to the QMol-grid package for basis-set models and, as such, contains many properties. For enhanced readability it implements most of it methods in separate files (inside the @QMol_disc_basis folder). Its efficiency as compared to a regular definition of the class and all its methods in a single MATLAB file are untested. To make a single-file class, copy all the QMol_disc_basis.m file class definition and move the file outside of @QMol_disc_basis folder (and delete the folder).

Run-time documentation

After it has been initialized, display the run-time documentation associated with a specific QMol_disc_basis object

obj.showDocumentation;

The run-time documentation is only accessible to the QMol_suite (and overloading) class.

Notes

The results displayed in this documentation page were generated using version 01.21 of the QMol-grid package.

QMol_disc_basis was introduced in version 01.00
Version 01.10 moved DFT_getCharge and DFT_getDensity to QMol_DFT_density
Version 01.10 changed the names of DFT_orbitalEnergy to DFT_energyOrbital and DFT_kineticEnergy to DFT_energyKinetic
Version 01.10 changed the names of DFT_projectKSO to DFT_projectOrbital and DFT_getKSO to DFT_reconstructOrbital
getMemoryProfile was introduced in version 01.10
Schrödinger-equation methods were added in version 01.20

QMol_disc_basis - fmauger1/QMol-grid GitHub Wiki

QMol_disc_basis

Table of Contents

Description

Class properties

Domain definition

basis (v)

xspan (x)

Differential operators

D

T

Tv

mD

mT

mTv

Initialization status

isInitialized (isInit)

Miscellaneous

dx

basisSize (nV)

QM

isBasis

isBasisPol

Generic class methods

Creation

default constructor

Changing class properties

set

reset

clear

orthonormalizeBasis

Initializing QMol_disc_basis objects

initialize

setTv

Run-time documentation

getMemoryProfile

Comparing domain discretizations

eq (==)

ne (~=)

DFT class methods

Component names (static)

DFT_classDensity (static)

DFT_classOrbital (static)

DFT_classPotential (static)

DFT_classPotentialGradient (static)

Data allocation

DFT_allocateDensity

DFT_allocateOrbital

DFT_allocatePotential

DFT_allocatePotentialGradient

Kohn-Sham orbitals

DFT_projectOrbital

DFT_reconstructOrbital

Miscellaneous

DFT_dist

DFT_mix

DFT_normalizeOrbital

DFT_randomOrbital

DFT_SCF_AndersonMixCoeff

DFT_sizeOrbital

Operators and functionals

DFT_energyKinetic

DFT_energyOrbital

DFT_operatorHamiltonian

DFT_operatorKinetic

DFT_operatorPotential

Schrödinger-equation class methods

Component names (static)

SE_classWaveFunction (static)

SE_classPotential (static)

Data allocation

SE_allocateWaveFunction

Wave functions

SE_projectWaveFunction

SE_reconstructWaveFunction

Miscellaneous

SE_normalizeWaveFunction

SE_randomWaveFunction

SE_sizeWaveFunction

Operators and functionals

`QMol_disc_basis`

`basis (v)`

`xspan (x)`

`D`

`T`

`Tv`

`mD`

`mT`

`mTv`

`isInitialized (isInit)`

`dx`

`basisSize (nV)`

`QM`

`isBasis`

`isBasisPol`

`set`

`reset`

`clear`

`orthonormalizeBasis`

Initializing `QMol_disc_basis` objects

`initialize`

`setTv`

`getMemoryProfile`

`eq (==)`

`ne (~=)`

`DFT_classDensity (static)`

`DFT_classOrbital (static)`

`DFT_classPotential (static)`

`DFT_classPotentialGradient (static)`

`DFT_allocateDensity`

`DFT_allocateOrbital`

`DFT_allocatePotential`

`DFT_allocatePotentialGradient`

`DFT_projectOrbital`

`DFT_reconstructOrbital`

`DFT_dist`

`DFT_mix`

`DFT_normalizeOrbital`

`DFT_randomOrbital`

`DFT_SCF_AndersonMixCoeff`

`DFT_sizeOrbital`

`DFT_energyKinetic`

`DFT_energyOrbital`

`DFT_operatorHamiltonian`

`DFT_operatorKinetic`

`DFT_operatorPotential`

`SE_classWaveFunction (static)`

`SE_classPotential (static)`

`SE_allocateWaveFunction`

`SE_projectWaveFunction`

`SE_reconstructWaveFunction`

`SE_normalizeWaveFunction`

`SE_randomWaveFunction`

`SE_sizeWaveFunction`

`SE_energyKinetic`

`SE_energyWaveFunction`

`SE_operatorHamiltonian`

`SE_operatorKinetic`

`SE_operatorPotential`