QMol_disc defines the Cartesian-grid discretization of the domain for all grid-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.

`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 results.

Differential operators

Differential-operator properties are dependent and hold meaningful values only after the QMol_disc 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 conditions, with ifft(obj.Tv.fft(f)).
Empty Tv indicates that the velocity-gauge laplacian has not been initialized using the setTv method.

Initialization status

Internally, QMol_disc 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 object, irrespective of prior initialization are guaranteed to define a meaningful initialization status.

`isInitialized (isInit)`

Initialization status of the QMol_disc 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 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 grid points in the domain discretization [ positive integer (default []) ]

At initialization, basisSize is defined as numel(obj.xspan).
basisSize is typically used for data preallocation (thus the name).

`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 methods -- such methods are identified below.
For developers: For practical reasons, QM is editable by QMol_DFT classes.

`isBasis`

Flag indicating that QMol_disc objects are not a basis-set discretization [ false ]

At run time, isBasis can be used to discriminate QMol_disc objects from QMol_disc_basis (which overloads the class)
For class development purposes, isBasis 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 object with all properties set to their default values

obj = QMol_disc;

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

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

Changing class properties

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

`set`

Update the name property of a QMol_disc 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, 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.

Initializing `QMol_disc` objects

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

`initialize`

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

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

This initializes the dependent properties D, T, 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 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 object has been initialized, set the velocity-gauge Laplacian operator Tv 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 with either of

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

Run-time documentation

`getMemoryProfile`

Get an estimate of the memory held by a QMol_disc 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 objects

`eq (==)`

Test if two QMol_disc objects describe the same domain discretization

obj1 == obj2

This is defined as all(abs(obj1.xspan(:)-obj2.xspan(:)) <= eqTol,'all')
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 returns a false result.
The equality test ignores any QMol-grid theory object specified in QM.

`ne (~=)`

Test if two QMol_disc 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 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'

`DFT_classPotential (static)`

Get the name of the class describing Kohn-Sham potentials

obj.DFT_classPotential;

Returns 'QMol_DFT_Vks'

`DFT_classPotentialGradient (static)`

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

obj.DFT_classPotentialGradient;

Returns 'QMol_DFT_Vks_grad'

Data allocation

Note: all DFT data/object allocation require the QMol_disc 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.

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 arease 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 object are allocated to zero.
The allocation also initializes the QMol_DFT_orbital 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 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 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 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 arease all data in the input QMol_DFT_Vks object V and then perform the (re)allocation.

`DFT_allocatePotentialGradient`

Allocate a new Kohn-Sham potential gradient object with either

DV = obj.DFT_allocatePotentialGradient;
DV = obj.DFT_allocatePotentialGradient([]);

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

Reallocate an existing Kohn-Sham potential gradient object

obj.DFT_allocatePotentialGradient(DV);

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

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
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
p = DFT_randomOrbital(randStr,0) is equivalent to p = DFT_randomOrbital(randStr) and produces a random orbital with no specific symmetry

Generate a random symmetric Kohn-Sham orbital

p = obj.DFT_randomOrbital(randStr,1);

Generate a random antisymmetric Kohn-Sham orbital

p = obj.DFT_randomOrbital(randStr,-1);

Random orbitals with a specific symmetry are typically used as a seed for ground-state DFT computations with enforced symmetry for the solution.

`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 = [numel(obj.x) 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 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 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 object KSO.
The input Vks is a QMol_DFT_Vks 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 with either

Hp = obj.DFT_operatorHamiltonian(Vks,p);
Hp = obj.DFT_operatorHamiltonian(Vks,p,0);

p and Hp are both numel(xspan)-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 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,true,0);
Hp = obj.DFT_operatorHamiltonian(Vks,p,false);  % spin down
Hp = obj.DFT_operatorHamiltonian(Vks,p,false,0);

p and Hp are both numel(xspan)-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 object describing the Kohn-Sham potential.
DFT_operatorHamiltonian automatically detect whether to use the length- or velocity-gauge kinetic operator.

Enforce symmetric or antisymmetric results respectively with

Hp = obj.DFT_operatorHamiltonian(___,1);
Hp = obj.DFT_operatorHamiltonian(___,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced DFT_operatorPotential on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state DFT computations with enforced symmetry.

`DFT_operatorKinetic`

Apply the kinetic operator to a wave function with either

Hp = obj.DFT_operatorKinetic(p);
Hp = obj.DFT_operatorKinetic(p,0);

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

Enforce symmetric or antisymmetric results respectively with

Hp = obj.DFT_operatorKinetic(p,1);
Hp = obj.DFT_operatorKinetic(p,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced DFT_operatorKinetic on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state DFT computations with enforced symmetry.

`DFT_operatorPotential`

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

Hp = obj.DFT_operatorPotential(Vks,p);
Hp = obj.DFT_operatorPotential(Vks,p,0);

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 numel(xspan)-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 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,true,0);
Hp = obj.DFT_operatorPotential(Vks,p,false);  % spin down
Hp = obj.DFT_operatorPotential(Vks,p,false,0);

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 numel(xspan)-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 object
This method interface is provided to completeness, instead Hp = Vks.applyPotential(p,true) and Hp = Vks.applyPotential(p,false) are recommended

Enforce symmetric or antisymmetric results respectively with

Hp = obj.DFT_operatorPotential(___,1);
Hp = obj.DFT_operatorPotential(___,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced DFT_operatorPotential on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state DFT computations with enforced symmetry.

Schrödinger-equation class methods

Component names (static)

`SE_classWaveFunction (static)`

Get the name of the class describing wave functions

obj.SE_classWaveFunction;

Returns 'QMol_SE_wfcn'

`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 object are allocated to zero.
The allocation also initializes the QMol_SE_wfcn 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 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 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

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
p = SE_randomWaveFunction(randStr,0) is equivalent to p = SE_randomWaveFunction(randStr) and produces a random orbital with no specific symmetry

Generate a random symmetric wave function

p = obj.SE_randomWaveFunction(randStr,1);

Generate a random antisymmetric wave function

p = obj.SE_randomWaveFunction(randStr,-1);

Random wave functions with a specific symmetry are typically used as a seed for ground-state Schrödinger-equation computations with enforced symmetry for the solution.

`SE_sizeWaveFunction`

Get the size of a wave function

s = obj.SE_sizeWaveFunction;

This is defined as s = [numel(obj.x) 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 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);
Hp = obj.SE_operatorHamiltonian(Vks,p,0);

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.

Enforce symmetric or antisymmetric results respectively with

Hp = obj.SE_operatorHamiltonian(___,1);
Hp = obj.SE_operatorHamiltonian(___,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced SE_operatorHamiltonian on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state Schrödinger-equation computations with enforced symmetry.

`SE_operatorKinetic`

Apply the kinetic operator to a wave function with either

Hp = obj.SE_operatorKinetic(p);
Hp = obj.SE_operatorKinetic(p,0);

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 numel(xspan)-by-1 vectors respectively containing the discretization of the wave functions to which the potential should be applied and its result.

Enforce symmetric or antisymmetric results respectively with

Hp = obj.SE_operatorKinetic(p,1);
Hp = obj.SE_operatorKinetic(p,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced SE_operatorKinetic on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state Schrödinger-equation computations with enforced symmetry.

`SE_operatorPotential`

Apply the potential operator to a wave function with either

Hp = obj.SE_operatorPotential(V,p);
Hp = obj.SE_operatorPotential(V,p,0);

This computes the action of the potential operator described by the V object on the input wave function 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

Enforce symmetric or antisymmetric results respectively with

Hp = obj.SE_operatorPotential(___,1);
Hp = obj.SE_operatorPotential(___,-1);

The symmetric constraint enforces that Hp matches flip(Hp,1) within roundoff error
The antisymmetric constraint enforces that Hp matches -flip(Hp,1) within roundoff error
In both cases, the input wave function should match the same symmetry condition, but no check is performed for this (i.e., using symmetry enforced SE_operatorPotential on a wave function that does not match the same symmetry condition will produce erroneous results)
This is typically used in ground-state Schrödinger-equation computations with enforced symmetry.

Examples

Create a discretization

disc = QMol_disc('x',-10:.1:15);
disc.initialize;

Display the estimated memory footprint for the object

disc.getMemoryProfile(true);

  * Domain discretization                                        
    > domain                                                        2.0 KB
    > gradient                                                      3.9 KB
    > Laplacian                                                     2.0 KB
    > Laplacian (velocity gauge)                                    3.9 KB

Compute the gradient and Laplacian of a function

f   =   3/sqrt(2*pi*3)*exp(-(disc.x(:)-3.5).^2/2/3);
df  =   ifft(disc.D.*fft(f));
ddf =   ifft(disc.T.*fft(f));

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 compare the numerical differentiations to closed-form formula.

% Analytical derivatives
Df  =  -(disc.x(:)-3.5)/3 .* f;
Ddf =   (f/3 - (disc.x(:)-3.5).^2/9 .*f)/2;

% Plot differences
figure; hold on
plot(disc.x,Df-df,'-','LineWidth',2','DisplayName','\nabla{f}')
plot(disc.x,Ddf-ddf,'-','LineWidth',2','DisplayName','-\Delta{f}/2')
xlabel('position (a.u.)'); xlim(disc.x([1 end]));
ylabel('error')
legend show

and see that we have excellent agreement between the two (the error can be further driven down by increasing the box size).

Test suite

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

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

For developers

Class definition

QMol_disc defines a concrete handle class overloading QMol_suite.

QMol_disc is a cornerstone component to the QMol-grid package and, as such, contains many properties. For enhanced readability it implements most of it methods in separate files (inside the @QMol_disc 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.m file class definition and move the file outside of @QMol_disc folder (and delete the folder).

Run-time documentation

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

obj.showDocumentation;

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

Private methods

QMol_disc defines a handful of private methods (in the QMol_disc.m file itself).

`DFT_applySymmetry`

Apply the symmetrization operator to a Kohn-Sham orbital

p = obj.DFT_applySymmetry(1,p);

The symmetrized orbital verifies p = flip(p,1) within roundoff error

Apply the antisymmetrization operator to a Kohn-Sham orbital

p = obj.DFT_applySymmetry(-1,p);

The antisymmetrized orbital verifies p = -flip(p,1) within roundoff error

Note that DFT_applySymmetry is an in-place function. To provide a universal interface, p = obj.DFT_applySymmetry(0,p) is the identity operator, imposing no (anti)symmetrization to the orbital.

`SE_applySymmetry`

Apply the symmetrization operator to a wave function

p = obj.SE_applySymmetry(1,p);

The symmetrized wave function verifies p = flip(p,1) within roundoff error

Apply the antisymmetrization operator to a wave function

p = obj.SE_applySymmetry(-1,p);

The antisymmetrized wave function verifies p = -flip(p,1) within roundoff error

Note that SE_applySymmetry is an in-place function. To provide a universal interface, p = obj.SE_applySymmetry(0,p) is the identity operator, imposing no (anti)symmetrization to the orbital.

Notes

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

QMol_disc 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
getMemoryProfile was introduced in version 01.10
Schrödinger-equation methods were added in version 01.20

QMol_disc - fmauger1/QMol-grid GitHub Wiki

QMol_disc

Table of Contents

Description

Class properties

Domain definition

xspan (x)

Differential operators

D

T

Tv

Initialization status

isInitialized (isInit)

Miscellaneous

dx

basisSize (nV)

QM

isBasis

Generic class methods

Creation

default constructor

Changing class properties

set

reset

clear

Initializing QMol_disc 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

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

Miscellaneous

SE_normalizeWaveFunction

SE_randomWaveFunction

SE_sizeWaveFunction

Operators and functionals

SE_energyKinetic

SE_energyWaveFunction

SE_operatorHamiltonian

SE_operatorKinetic

SE_operatorPotential

Examples

Test suite

For developers

Class definition

Run-time documentation

Private methods

DFT_applySymmetry

`QMol_disc`

`xspan (x)`

`D`

`T`

`Tv`

`isInitialized (isInit)`

`dx`

`basisSize (nV)`

`QM`

`isBasis`

`set`

`reset`

`clear`

Initializing `QMol_disc` 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_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_normalizeWaveFunction`

`SE_randomWaveFunction`

`SE_sizeWaveFunction`

`SE_energyKinetic`

`SE_energyWaveFunction`

`SE_operatorHamiltonian`

`SE_operatorKinetic`

`SE_operatorPotential`

`DFT_applySymmetry`

`SE_applySymmetry`