Class FlucDens
Defined in File FlucDens.h
Class Documentation
-
class FlucDens
To construct a new FlucDens force object, define the object with.
This class impliments the fluctuating monopole density interaction. Similar to a fluctuating charge model, this force uses atom centered monopole densities described by slater functions for the electrons and point-particles for the nuclei. As such, FlucDens includes effects of electron overlap between atomic sites.
To use this force, create a FlucDens force object by specifying the number of charge sites and the parameters for all sites. These parameters can be altered after the creation of the force object. Next, each site must be assigned a molecular “fragment” Each fragment, defined by it’s site indicies, are groups of charge sites whose dynamic desity coefficients will sum to zero when the electrostatic energy is minized. As a result, charge transfer is not allowed between fragments, but each fragment will however reduce to it’s frozen electron populations when complety separated from all other framgents. Each fragment is also used to define which frozen density - dynamic density interactions are included. Frozen densities on each site do not interact with the dynamic densities on each fragment, allowing for this limmiting behavior to occur.
Next, set the model’s pther settings by calling their respective setter function. This includes the polarization dampening strengths, charge-transfer coefficients, and additional hardness values to be applied to the Coulomb matrix. Polarization dampening is applied between a polarized density site and it’s interacing frozen density site and decays exponentially as a function of their distance. The strength and rate of this polarization can be adjusted, as well as whether this is applied to the linear
(default) or the quadratic matrix terms (old method). Charge transfer is not explicitly allowed in this force, and instead the energy associated with this interation is estimated by a scalar multiple of the total polarization energy.
Additional hardness values can be applied to the Coulomb matrix. These values are not overrides, but instead add a constant value to the self-interaction between two densities and hinder (or enhance) the ability of that site to gain or loose electrons. These values can be either positive or negative, with positive values increasing the hardness. Care must be taken that the system remains positive-definite if negative values are used.
Aconstant electric field can also be applied to each system. Simply call set_external_field() with the desired field components to apply this field. Periodic boundary conditions are also possible, and can be set by calling set_use_PBC(). If PBC is used, then a long-range cutoff distance must be used (the default is to use half the minimum periodic box lengths). Turning this off and using PBC will have undefined behavior.
This module also computes all interaction energies, including minimizing the total energy w.r.t. the dynamic density coefficients. Do run this computation, first call calc_energy() to compute all frozen-frozen energy components and to create the polarization matrix components. Then, call solve_minimization() to perform the actual matrix algebra associated with the minimization procedures. Energy components can then be returned calling get_energies() or any of the particular energy getter functions. Dipoles (including the polarized density) can also then be called with the updated electron populations.
All values in this module use atomic units unless otherwise specified.
- Param n_sites
number of particles (including cirtual sites)
- Param frozen_charges_in
the frozen charges of all sites
- Param nuclei_in
the nuclei of each site. Core electrons will be remove automatically
- Param frozen_exp
the Slater exponents of all frozen densities
- Param dynamic_exp
the Slater exponents of all dynamic densities
Public Types
Public Functions
-
FlucDens(const int n_sites, const double *frozen_charges_in, const double *nuclei_in, const double *frozen_exp, const double *dynamic_exp)
Construct a new FlucDens force object.
- Parameters
n_sites – number of particles (including cirtual sites)
frozen_charges_in – array of frozen charges of all sites
nuclei_in – array of nuclei of each site. Core electrons will be remove automatically
frozen_exp – array of Slater exponents of all frozen densities
dynamic_exp – array of Slater exponents of all dynamic densities
-
~FlucDens()
-
vec_d calc_density(const vec_d &points, const vec_d &pos, DensityType density_type)
Compute the electron density at the specified coordinates.
- Parameters
points – A 3K array of K number of points to evaluate the density at
pos – The 3N array of positions for the charge sites of the system
density_type – The type of density to compute
- Returns
vec_d A K-length array of density evaluations
-
void set_site_params(const int index, const double frz_chg, const double frz_exp, const double dyn_exp)
Change the parameters for one charge site.
- Parameters
index – The index of the site to adjust the parameters
frz_chg – The frozen charge to set
frz_exp – The Slater exponent of the frozen density
dyn_exp – The Slater exponent of the dynamic density
-
void get_site_params(const int index, double &frz_chg, double &frz_exp, double &dyn_exp)
Get the parameters of an atomic site.
- Parameters
index – The index of the site to adjust the parameters
frz_chg – The frozen charge of the site
frz_exp – The Slater exponent of the frozen density
dyn_exp – The Slater exponent of the dynamic density
-
void set_dyn_exp(const int index, const double value)
Set the Slater exponent of a particular dynamic density site.
- Parameters
index – The index of the site
value – The value of the Slater exponent
-
void set_dyn_exp(vec_d exponents)
Set the Slater exponents of all dynamic density sites at once.
- Parameters
exponents – An array of all exponents in order of their indicies
-
void set_frz_exp(const int index, const double value)
Set the Slater exponent of a particular frozen density site.
- Parameters
index – The index of the site
value – The value of the Slater exponent
-
void set_ct_coeff(const double coeff)
Set the charge-transfer coefficient. Charge transfer is extimated as a multiple of the polarization energy between two inducable sites. This does not apply to the polarization due to external fields.
- Parameters
coeff – The scalar multiple of the polarization energy used for the charge transfer energy
-
double get_ct_coeff()
Get the charge-transfer coefficient.
- Returns
double The coefficient
-
void set_dampening(double coeff, double exponent, DampType damp = DampType::Linear)
Set the polarization dampening parameters.
- Parameters
coeff – The strength of the dampening
exponent – The exponential decay of the dampening
damp – Tye type of dampening, enumerated by DampType
-
void get_dampening(double &coeff, double &exponent)
Return the polarization dampening parameters.
- Parameters
coeff – The strength of the dampening
exponent – The exponential decay of the dampening
-
void set_additional_hardness(const int index, const double value)
Add additional atomic hardness to the diagonal J-matrix elements. This is added to the value already computed by the self interaction in the matrix, not an override.
- Parameters
index – The particle to add the hardness to.
value – The value of the additional hardness
-
void set_additional_hardness(vec_d values)
Return the additional hardness values of all sites.
- Parameters
values – An array of the hardness values
-
void add_fragment(const std::vector<int> site_idx_list)
Create a new fragment whoes dynamic density whould sum to zero. Every site must be assigned to a fragment, else an exception will be thrown.
- Parameters
site_idx_list – A list of site indicies that make up this fragment
-
std::vector<vec_i> get_fragments()
Return the partitioned fragments and their site indicies for each.
- Returns
An array of fragment indicies. Each element correpsonds to a fragment and contains a list of indicies that make up that fragment
-
int get_num_fragments()
Get the num fragments.
- Returns
The number of fragments
-
void add_del_frz_exclusion(int delta_i, int frz_j)
Exclude the frozen electrostatic interactions between a dynamic and frozen site.
- Parameters
delta_i – The index of the dynamic density site
frz_j – The index of the frozen density site
-
std::set<int> get_del_frz_exclusions(const int particle1) const
Get the indicies of all frozen sites that do not interact with a single dynamic density.
- Parameters
particle1 – The index of the dynamic density
- Returns
std::set<int> The indicies of the frozen sites that do not interact with this dynamic density
-
void set_frag_constraints(const bool constr_frags)
Choose whether constraints are used on all fragments or not. NOT FULLY IMPLIMENTED YET!
- Parameters
constr_frags – If true, fragment constraints will be used. If false, then only the total dynamic density of the entire system will be summed to zero, not per fragment.
-
int get_num_constraints()
Get the number of total constrains applied to the energy minimization.
- Returns
The total number of constraints
-
std::vector<std::vector<int>> get_constraints()
Get constraint matrix applied to the dyanmic densities.
- Returns
std::vector<std::vector<int>>
-
void add_frz_frz_exclusion(int frz_i, int frz_j)
Exclude the interaction of two frozen densities.
- Parameters
frz_i – The index of the first frozen site
frz_j – The index of the second frozen site
-
std::set<int> get_frz_frz_exclusions(const int particle1) const
Get the all sites that are excluded from interacting from a particular site.
- Parameters
particle1 – The frozen site whoes exclusions are being inquired
- Returns
std::set<int> All indicies that are excluded
-
int get_num_frz_frz_exclusions() const
Get the number of frz-frz exclusions.
- Returns
int The total number of exclusions in the system
-
void create_frz_exclusions_from_bonds(const std::vector<std::pair<int, int>> bonds, int bond_cutoff)
Create frozen-frozen exclusions based on each site’s bonds.
- Parameters
bonds – An vector of pair of indicies for each bond
bond_cutoff – All indicies that are this number of bonds away will be excluded
-
void set_external_field(double field_x, double field_y, double field_z)
Apply a constant electric field to the entire system.
- Parameters
field_x – The x-direction Cartesian component of the field
field_y – The y-direction Cartesian component of the field
field_z – The z-direction Cartesian component of the field
-
void apply_field_to_system(const vec_d &coords)
Apply the field to the system. This must be called before any minimization of the total energy is performed.
- Parameters
coords – 1D array of positions of the all sites
-
std::vector<Vec3> get_dipoles(const vec_d &coords)
Return the electric dipoles of the system. This includes nuclear, frozen electrons, dynamic electrons, and total dipoles.
- Parameters
coords – 1D array of positions of the all sites
- Returns
std::vector<Vec3> The dipoles that are computed. Indicies are determind by DensityType enumeration
-
vec_d get_dipole(const vec_d &coords, DensityType density_dype)
Compute the dipole for a particule type of charge.
- Parameters
coords – 1D array of positions of the all sites
density_dype – DensityType: The type of charges to be used in the computation
- Returns
The dipole components
-
double calc_frz_ext_field_energy(const vec_d &coords, std::vector<Vec3> &forces)
Compute the energy and forces of the frozen charges (nuclei and frozen electrons) with the external electric field.
- Parameters
coords – 1D array of positions of the all sites
forces – Array of forces to be updated
- Returns
The value of the interaction energy
-
void set_use_PBC(const bool is_periodic)
Turn on or off periodic boundary conditions.
- Parameters
is_periodic – If true, use PBC, if false, do not use them
-
void set_use_PBC(const bool is_periodic, const double x, const double y, const double z)
Set the periodic boundaries.
- Parameters
is_periodic – If true, use PBC, if false, do not use them
x – The x-direction of the periodic box
y – The y-direction of the periodic box
z – The z-direction of the periodic box
-
bool get_use_PBC()
Inquire if PBC is used or not.
- Returns
true PBC is being used
- Returns
false PBC is not being used
-
void set_use_cutoff(bool useCutoff)
Set whether or not to use a nonbonded cutoff. Care must be taked to ensure that the polarization J-matrix remains is positive definite.
- Parameters
useCutoff – If true, enable the use of a cutoff distance between all interactions. If false, do not use any cutoff distance.
-
bool get_use_cutoff()
Inquire of nonbonded cutoffs are being used.
- Returns
true cutoffs are being used
- Returns
false cutoffs are not being used
-
void set_cutoff_distance(double distance_in_nm)
Set the cutoff distance (in bohr)
- Parameters
distance_in_nm – The cutoff distance (in bohr)
-
double get_cutoff_distance()
Get the cutoff distance (in bohr)
- Returns
distance_in_nm The cutoff distance (in bohr)
-
void set_use_SR_cutoff(bool use_cutoff)
Set the use SR cutoff. If enabled, electron overlap will be ignored beyond a particular distance.
- Parameters
use_cutoff – If true, enable the use this apprximation
-
bool get_use_SR_cutoff()
Inquire of SR cutoffs are being used.
- Returns
true SR cutoffs are being used
- Returns
false SR cutoffs are not being used
-
double calc_overlap(const vec_d &coords)
Calculate the total electron overlap.
- Parameters
coords – the positions of all sites
- Returns
the total overlap in e^2 /bohr^3
-
double calc_energy(const vec_d &coords, bool calc_frz = true, bool calc_pol = true)
Calculate the energy and forces of the system.
- Parameters
coords – A 1D array of all the positions (in bohr)
calc_frz – Whether or not to compute frozen electrostatic energies
calc_pol – Whether or not to compute polarization
- Returns
The total frozen energy computed. Must cal solve_minimization() in order to get polarization energy.
-
void calc_one_electro(DeltaR &deltaR, int i, int j, bool calc_pol, bool calc_frz, Energies &energies, std::vector<Vec3> &forces, int thread_num = 0)
Compute one pair of interactions between two atomic sites.
- Parameters
deltaR – The distance object between the two sites
i – The index of the first site
j – The index of the second site
calc_pol – Whether or not to include polarization
calc_frz – Whether or not to include frozen electrostatic energies
energies – Returns the energies computed
forces – Returns the forces between both sites
thread_num – the index of the OpenMP thread to compute the energies and forces on
-
Energies calc_one_frozen(const vec_d &coords, int i, int j)
Compute one pair of frozen interactions between two atomic sites (does not include polarization)
- Parameters
coords – 1D array of all the positions (in bohr)
i – The index of the first site
j – The index of the second site
- Returns
Energies energies of the computed pair
-
void initialize_calculation()
initizlize arrays before an energy computation is performed. This is called automatically by calc_energy(), but must be called manually to clear all energy and force arrays if calling calc_one_electro() or calc_one_frozen().
-
void solve_minimization(std::vector<Vec3> &forces)
Minimize the total energy of the system w.r.t. the fluctuating densities. This function solve the matrix equations that minimizes the quadratic polarization energy, under the constraint that all fragment dynamic densities sum to zero.
- Parameters
forces – array of forces to be updated based on the polarization interactions
-
double get_frozen_energy()
Get the energy associated with frozen-frozen dsitity interactions.
- Returns
double frozen-frozen energy
-
double get_polarization_energy()
Get the polarization energy.
- Returns
double polarization energy
-
double get_ct_energy()
Get the estimated CT energy as a multiple of the polarization energy.
- Returns
double
-
FlucDensEnergies get_energies()
Return all energy components computed.
- Returns
FlucDensEnergies Object that constains the energy components
-
std::vector<vec_d> get_forces()
Get the forces exerted on all particles after an energy computation is performed.
- Returns
std::vector<vec_d> Nx3 Array of forces on all sites
-
double frz_frz_overlap(const double inv_r, const double a, const double b, const double exp_ar, const double exp_br)
Overlap between frozen-frozen dsnsity.
- Parameters
inv_r – inverse distance between the two sites
a – frozen exponent of the first site
b – frozen exponent of the second site
exp_ar – pre-computed exp(-a*r) using the first exponent
exp_br – pre-computed exp(-a*r) using the secpnd exponent
- Returns
double the overlap in e^2/bohr^3
-
double elec_elec_energy(const double inv_r, const double a, const double b, const double exp_ar, const double exp_br, double &dEdR)
Compute the density-density interaction, either frozen or dynamic.
- Parameters
inv_r – inverse distance between the two sites
a – frozen exponent of the first site
b – frozen exponent of the second site
exp_ar – pre-computed exp(-a*r) using the first exponent
exp_br – pre-computed exp(-a*r) using the secpnd exponent
dEdR – The derivative of the energy w.r.t. the distance
- Returns
The elec-elec energy in a.u.
-
double elec_nuclei_energy(const double inv_r, const double a, const double exp_ar, double &dEdR)
Compute the energy between a nuclei and an electron density.
- Parameters
inv_r – inverse distance between the two sites
a – frozen exponent of the density site
exp_ar – pre-computed exp(-a*r) using the density exponent
dEdR – The derivative of the energy w.r.t. the distance
- Returns
The nuclei-density energy in q.u.
-
vec_d get_rho_coulomb_mat()
Get the delta-rho Coumomb J-matrix.
- Returns
1D array of the symmetric matrix elements. Contains NxN elements
-
vec_d get_delta_rho()
Get the delta rhos of each site that are computed after minimizing the energy.
- Returns
vec_d Nx1 array
-
double get_total_time()
Not implimented yet! Return the total wall time used for computations.
- Parameters
calculate_forces –
-
void set_calc_forces(bool calculate_forces)
Set wether or not forces should also be computed. If set to false, some computation time can be saved.
- Parameters
calculate_forces – If true, then forces will be computed
-
void print_params(const std::string message, const std::string param_name)
Print site parameters to the console.
- Parameters
message – string: A comment message to be printed
param_name – string: the name of the parameter to be printed
-
std::vector<std::string> get_param_names()
Get the name of the parameters that can be return.
- Returns
std::vector<std::string> A list of possible parameters to enquire about