diff --git a/src/fmo/gradients/embedding.rs b/src/fmo/gradients/embedding.rs index 3dbd9b0..37b0018 100644 --- a/src/fmo/gradients/embedding.rs +++ b/src/fmo/gradients/embedding.rs @@ -7,42 +7,31 @@ use nalgebra::Vector3; use ndarray::prelude::*; use ndarray::RawData; use std::iter::FromIterator; -use std::ops::AddAssign; -use crate::utils::array_helper::ToOwnedF; +use std::ops::{AddAssign, SubAssign}; impl SuperSystem { + + /// Computes and returns the gradient of the embedding energy. pub fn embedding_gradient(&mut self) -> Array1 { - // initialize the gradient of the embedding energy array with zeros + + // The gradient of the embedding energy is initialized as an array with zeros. let mut gradient: Array1 = Array1::zeros([3 * self.atoms.len()]); + + // A reference to the charge differences and gamma matrix for all atoms is created. let dq: ArrayView1 = self.properties.dq().unwrap(); let gamma: ArrayView2 = self.properties.gamma().unwrap(); - // broadcast the charge differences into the shape of the gradients (n_atoms -> 3 * n_atoms) + // The charge differences are broadcast into the shape the gradients. let dq_f: Array1 = dq .broadcast([3, self.atoms.len()]) .unwrap() .reversed_axes() - .to_owned_f() + .to_owned() .into_shape([3 * self.atoms.len()]) .unwrap(); - // combine all derivatives of the charge differences into one Array that has the shape of - // the total gradient. Only the diagonal elements are used. - let mut grad_dq: Array1 = Array1::zeros([3 * self.atoms.len()]); - - // PARALLEL or this could also be done at once instead of getting the diag in the loop. However - // the disadvantage would be that we need to create the 3D Array and store it temporarily. - for mol in self.monomers.iter() { - let mol_grad_dq: ArrayView3 = mol - .properties - .grad_dq() - .unwrap() - .into_shape([3, mol.n_atoms, mol.n_atoms]) - .unwrap(); - grad_dq - .slice_mut(s![mol.slice.grad]) - .assign(&diag_of_last_dimensions(mol_grad_dq, mol.n_atoms)); - } + // The derivative of the charge differences is initialized as an array with zeros. + let mut grad_dq: ArrayView1 = self.properties.grad_dq_diag().unwrap(); // TODO: it is not neccessary to calculate the derivative of gamma two times. this should be // improved! it is already computed in the gradient of the monomer/pair @@ -50,26 +39,30 @@ impl SuperSystem { gamma_gradients_atomwise_2d(&self.gammafunction, &self.atoms, self.atoms.len()); let mut grad_gamma_dot_dq: Array1 = grad_gamma_sparse.dot(&dq); - // compute the gradient of the embedding energy for each pair, that is not treated within - // the ES-dimer approximation PARALLEL + // Begin of the loop to compute the gradient of the embedding energy for each pair. + // PARALLEL for pair in self.pairs.iter() { - // get references to the corresponding monomers + // References to the corresponding monomers. let m_i: &Monomer = &self.monomers[pair.i]; let m_j: &Monomer = &self.monomers[pair.j]; - // if the derivative is w.r.t to an atom that is within this pair: - // the first part of the equation reads: + // If the derivative is w.r.t to an atom that is within this pair: + // The first part of the equation reads: // dDeltaE_IJ^V/dR_a x = DDq_a^IJ sum_(K!=I,J)^(N) sum_(C in K) Dq_C^K dgamma_(a C)/dR_(a x) + + // Reference to the DDq_a^IJ (difference of charge difference between pair and monomer) let delta_dq: ArrayView1 = pair.properties.delta_dq().unwrap(); - // broadcast Delta dq into the shape of the gradients + + // DDq is broadcasted into the shape of the gradients. let delta_dq_f: Array1 = delta_dq .broadcast([3, pair.n_atoms]) .unwrap() .reversed_axes() - .to_owned_f() + .to_owned() .into_shape([3 * pair.n_atoms]) .unwrap(); + // The gradient for a in I is computed and assigned. gradient.slice_mut(s![m_i.slice.grad]).add_assign( &(&delta_dq_f.slice(s![..3 * m_i.n_atoms]) * &(&grad_gamma_dot_dq.slice(s![m_i.slice.grad]) @@ -78,6 +71,7 @@ impl SuperSystem { .dot(&dq.slice(s![m_i.slice.atom])))), ); + // The gradient for a in J is computed and assigned. gradient.slice_mut(s![m_j.slice.grad]).add_assign( &(&delta_dq_f.slice(s![3 * m_i.n_atoms..]) * &(&grad_gamma_dot_dq.slice(s![m_j.slice.grad]) @@ -86,51 +80,64 @@ impl SuperSystem { .dot(&dq.slice(s![m_j.slice.atom])))), ); - // the second term... - let grad_delta_dq: Array1 = get_grad_delta_dq(pair, m_i, m_j); - let mut esp_ij: Array1 = Array1::zeros([pair.n_atoms]); - esp_ij - .slice_mut(s![..m_i.n_atoms]) - .assign(&m_i.properties.esp_q().unwrap()); - esp_ij - .slice_mut(s![m_i.n_atoms..]) - .assign(&m_j.properties.esp_q().unwrap()); - let esp_ij: Array1 = esp_ij - .broadcast([3, pair.n_atoms]) - .unwrap() - .reversed_axes() - .to_owned_f() - .into_shape([3 * pair.n_atoms]) - .unwrap(); - let gddq_esp: Array1 = &grad_delta_dq * &esp_ij; - gradient - .slice_mut(s![m_i.slice.grad]) - .add_assign(&gddq_esp.slice(s![..3 * m_i.n_atoms])); - - gradient - .slice_mut(s![m_j.slice.grad]) - .add_assign(&gddq_esp.slice(s![3 * m_i.n_atoms..])); - - // if the derivative is w.r.t to an atom that is not in this pair, a -> K where K != I,J - // left hand side - - // A in monomer I + // Right hand side of Eq. 24, but a is still in I,J + // The difference between the derivative of the charge differences between the monomer + // and the dimer is computed. + // let grad_delta_dq: Array1 = get_grad_delta_dq(pair, m_i, m_j); + // + // // The electrostatic potential (ESP) is collected from the corresponding monomers. + // let mut esp_ij: Array1 = Array1::zeros([pair.n_atoms]); + // esp_ij + // .slice_mut(s![..m_i.n_atoms]) + // .assign(&m_i.properties.esp_q().unwrap()); + // esp_ij + // .slice_mut(s![m_i.n_atoms..]) + // .assign(&m_j.properties.esp_q().unwrap()); + // + // // The ESP is transformed into the shape of the gradient. + // let esp_ij: Array1 = esp_ij + // .broadcast([3, pair.n_atoms]) + // .unwrap() + // .reversed_axes() + // .to_owned_f() + // .into_shape([3 * pair.n_atoms]) + // .unwrap(); + // + // // The (elementwise) product of the ESP with the derivative of the pair charge + // // differences is computed. + // let gddq_esp: Array1 = &grad_delta_dq * &esp_ij; + // + // // The gradient of the rhs for a in I is assigned. + // gradient + // .slice_mut(s![m_i.slice.grad]) + // .add_assign(&gddq_esp.slice(s![..3 * m_i.n_atoms])); + // + // // The gradient of the rhs for a in J is assigned. + // gradient + // .slice_mut(s![m_j.slice.grad]) + // .add_assign(&gddq_esp.slice(s![3 * m_i.n_atoms..])); + + // Start of the computation if the derivative is w.r.t to an atom that is not in + // this pair. So that a in K where K != I,J. + + // The matrix vector product of the gamma matrix derivative with DDq is computed for a in I let mut dg_ddq: Array1 = grad_gamma_sparse .slice(s![0.., m_i.slice.atom]) .dot(&delta_dq.slice(s![..m_i.n_atoms])); - // A in monomer J + // and for a in J. dg_ddq += &grad_gamma_sparse .slice(s![0.., m_j.slice.atom]) .dot(&delta_dq.slice(s![m_i.n_atoms..])); - // since K != I,J => the elements were K = I,J are set to zero + // Since K != I,J => the elements where K = I,J are set to zero. dg_ddq.slice_mut(s![m_i.slice.grad]).assign(&Array1::zeros([3 * m_i.n_atoms])); dg_ddq.slice_mut(s![m_j.slice.grad]).assign(&Array1::zeros([3 * m_j.n_atoms])); + // The (elementwise) product with the charge differences is computed and assigned. gradient += &(&dg_ddq * &dq_f); - // right hand side + // Start of the computation of the right hand side of Eq. 25. // A in monomer I let mut ddq_gamma: Array1 = delta_dq .slice(s![..m_i.n_atoms]) @@ -141,7 +148,7 @@ impl SuperSystem { .slice(s![m_i.n_atoms..]) .dot(&gamma.slice(s![m_j.slice.atom, 0..])); - // since K != I,J => the elements were K = I,J are set to zero + // Since K != I,J => the elements were K = I,J are set to zero. ddq_gamma.slice_mut(s![m_i.slice.atom]).assign(&Array1::zeros([m_i.n_atoms])); ddq_gamma.slice_mut(s![m_j.slice.atom]).assign(&Array1::zeros([m_j.n_atoms])); @@ -152,7 +159,7 @@ impl SuperSystem { .broadcast([3, self.atoms.len()]) .unwrap() .reversed_axes() - .to_owned_f() + .to_owned() .into_shape([3 * self.atoms.len()]) .unwrap()); } @@ -167,38 +174,43 @@ fn get_grad_delta_dq(pair: &Pair, m_i: &Monomer, m_j: &Monomer) -> Array1 { let grad_dq_i: ArrayView2 = m_i.properties.grad_dq().unwrap(); let grad_dq_j: ArrayView2 = m_j.properties.grad_dq().unwrap(); + // compute the difference between dimers and monomers and take the diagonal values - let mut grad_delta_dq_2d: Array2 = Array2::zeros([3 * pair.n_atoms, pair.n_atoms]); + let mut grad_delta_dq_2d: Array2 = grad_dq.to_owned(); // difference for monomer i grad_delta_dq_2d .slice_mut(s![..(3 * m_i.n_atoms), ..m_i.n_atoms]) - .assign(&(&grad_dq.slice(s![..(3 * m_i.n_atoms), ..m_i.n_atoms]) - &grad_dq_i)); + .sub_assign(&grad_dq_i); // difference for monomer j grad_delta_dq_2d .slice_mut(s![(3 * m_i.n_atoms).., m_i.n_atoms..]) - .assign(&(&grad_dq.slice(s![(3 * m_i.n_atoms).., m_i.n_atoms..]) - &grad_dq_j)); + .sub_assign(&grad_dq_j); + let grad_delta_dq_3d: Array3 = grad_delta_dq_2d .into_shape([3, pair.n_atoms, pair.n_atoms]) .unwrap(); - diag_of_last_dimensions(grad_delta_dq_3d, pair.n_atoms) + diag_of_last_dimensions(grad_delta_dq_3d) } -fn diag_of_last_dimensions(data: ArrayBase, n_atoms: usize) -> Array1 +pub fn diag_of_last_dimensions(data: ArrayBase) -> Array1 where S: ndarray::Data, { - // create a temporary array to store the values it will be flattened afterwards - let mut grad_charge: Array2 = Array2::zeros([3, n_atoms]); + let (a, b, c): (usize, usize, usize) = data.dim(); + assert_eq!(b, c, "The last two dimension should have the same length"); + + // A temporary array to store the values is created. + let mut grad_charge: Array2 = Array2::zeros([a, b]); - // take the diagonal of each of the three dimensions - for i in 0..3 { + // The diagonal of each of the three dimensions is saved. + for i in 0..a { grad_charge .slice_mut(s![i, ..]) .assign(&data.slice(s![i, .., ..]).diag()); } - // reshape it into a one dimensional array - grad_charge.into_shape([3 * n_atoms]).unwrap() + // The gradient of the charges are reshaped into a one dimensional array. + grad_charge.into_shape([a * b]).unwrap() } diff --git a/src/fmo/gradients/es_dimer.rs b/src/fmo/gradients/es_dimer.rs new file mode 100644 index 0000000..03899d2 --- /dev/null +++ b/src/fmo/gradients/es_dimer.rs @@ -0,0 +1,88 @@ +use crate::fmo::*; +use crate::initialization::*; +use crate::scc::gamma_approximation::gamma_gradients_atomwise_2d; +use crate::scc::*; +use crate::utils::array_helper::ToOwnedF; +use crate::utils::Timer; +use nalgebra::Vector3; +use ndarray::prelude::*; +use ndarray::RawData; +use std::iter::FromIterator; +use std::ops::{AddAssign, SubAssign}; + +impl SuperSystem { + /// Computes and returns the gradient of the embedding energy. + pub fn es_dimer_gradient(&mut self) -> Array1 { + let mut gradient: Array1 = Array1::zeros([3 * self.atoms.len()]); + + // A reference to the charge differences and gamma matrix for all atoms is created. + let dq: ArrayView1 = self.properties.dq().unwrap(); + let grad_dq: ArrayView1 = self.properties.grad_dq_diag().unwrap(); + + let gamma: ArrayView2 = self.properties.gamma().unwrap(); + + // The charge differences are broadcast into the shape the gradients. + let dq_f: Array1 = dq + .broadcast([3, self.atoms.len()]) + .unwrap() + .reversed_axes() + .to_owned() + .into_shape([3 * self.atoms.len()]) + .unwrap(); + + // TODO: it is not neccessary to calculate the derivative of gamma two times. this should be + // improved! it is already computed in the gradient of the monomer/pair + let grad_gamma_sparse: Array2 = + gamma_gradients_atomwise_2d(&self.gammafunction, &self.atoms, self.atoms.len()); + + for pair in self.esd_pairs.iter() { + // References to the corresponding monomers + let m_i: &Monomer = &self.monomers[pair.i]; + let m_j: &Monomer = &self.monomers[pair.j]; + + // a in I + let mut gradient_slice: ArrayViewMut1 = gradient.slice_mut(s![m_i.slice.grad]); + + // lhs. of Eq. 29 in Ref. [1] + gradient_slice += &(&dq_f.slice(s![m_i.slice.grad]) + * &grad_gamma_sparse + .slice(s![m_i.slice.grad, m_j.slice.atom]) + .dot(&dq.slice(s![m_j.slice.atom]))); + + // rhs of Eq. 29 in Ref. [1] + gradient_slice += &(&grad_dq.slice(s![m_i.slice.grad]) + * &gamma + .slice(s![m_i.slice.atom, m_j.slice.atom]) + .dot(&dq.slice(s![m_j.slice.atom])) + .broadcast([3, m_i.n_atoms]) + .unwrap() + .reversed_axes() + .to_owned_f() + .into_shape([3 * m_i.n_atoms]) + .unwrap()); + + // a in J + let mut gradient_slice: ArrayViewMut1 = gradient.slice_mut(s![m_j.slice.grad]); + + // lhs of Eq. 30 in Ref. [1] + gradient_slice += &(&dq_f.slice(s![m_j.slice.grad]) + * &grad_gamma_sparse + .slice(s![m_j.slice.grad, m_i.slice.atom]) + .dot(&dq.slice(s![m_i.slice.atom]))); + + // rhs of Eq. 30 in Ref [1] + gradient_slice += &(&grad_dq.slice(s![m_j.slice.grad]) + * &gamma + .slice(s![m_j.slice.atom, m_i.slice.atom]) + .dot(&dq.slice(s![m_i.slice.atom])) + .broadcast([3, m_j.n_atoms]) + .unwrap() + .reversed_axes() + .to_owned_f() + .into_shape([3 * m_j.n_atoms]) + .unwrap()); + } + + gradient + } +} diff --git a/src/fmo/gradients/mod.rs b/src/fmo/gradients/mod.rs index 0ff64a5..a866f00 100644 --- a/src/fmo/gradients/mod.rs +++ b/src/fmo/gradients/mod.rs @@ -11,14 +11,16 @@ mod embedding; mod monomer; mod pair; mod numerical; +mod es_dimer; // mod numerical; pub use monomer::*; pub use pair::*; use crate::fmo::helpers::get_pair_slice; +use crate::fmo::gradients::embedding::diag_of_last_dimensions; pub trait GroundStateGradient { - fn get_grad_dq(&self, atoms: &[Atom], grad_s: ArrayView3, p: ArrayView2) -> Array2; + fn get_grad_dq(&self, atoms: &[Atom], s: ArrayView2, grad_s: ArrayView3, p: ArrayView2) -> Array2; fn scc_gradient(&mut self, atoms: &[Atom]) -> Array1; } @@ -36,11 +38,30 @@ impl SuperSystem { fn monomer_gradients(&mut self) -> Array1 { let mut gradient: Array1 = Array1::zeros([3 * self.atoms.len()]); let atoms: &[Atom] = &self.atoms[..]; + // The derivative of the charge differences is initialized as an array with zeros. + let mut grad_dq: Array1 = Array1::zeros([3 * self.atoms.len()]); + + // PARALLEL or this could also be done at once instead of getting the diag in the loop. However + // the disadvantage would be that we need to create the 3D Array and store it temporarily. + // The derivative is collected from the monomers and only the diagonal elements are used. for mol in self.monomers.iter_mut() { gradient .slice_mut(s![mol.slice.grad]) .assign(&mol.scc_gradient(&atoms[mol.slice.atom_as_range()])); + + let mol_grad_dq: ArrayView3 = mol + .properties + .grad_dq() + .unwrap() + .into_shape([3, mol.n_atoms, mol.n_atoms]) + .unwrap(); + + grad_dq + .slice_mut(s![mol.slice.grad]) + .assign(&diag_of_last_dimensions(mol_grad_dq)); } + self.properties.set_grad_dq_diag(grad_dq); + return gradient; } diff --git a/src/fmo/gradients/monomer.rs b/src/fmo/gradients/monomer.rs index d1fbad5..1ef4a21 100644 --- a/src/fmo/gradients/monomer.rs +++ b/src/fmo/gradients/monomer.rs @@ -46,7 +46,9 @@ impl GroundStateGradient for Monomer { // the derivatives of the charge (difference)s are computed at this point, since they depend // on the derivative of S and this is available here at no additional cost. let p: ArrayView2 = self.properties.p().unwrap(); - let grad_dq: Array2 = self.get_grad_dq(&atoms, grad_s.view(), p.view()); + let s: ArrayView2 = self.properties.s().unwrap(); + let grad_dq: Array2 = self.get_grad_dq(&atoms, s.view(), grad_s.view(), p.view()); + drop(s); self.properties.set_grad_dq(grad_dq); // and reshape them into a 2D array. the last two dimension (number of orbitals) are compressed @@ -121,7 +123,7 @@ impl GroundStateGradient for Monomer { /// monomer. This means that the first dimension of the Array is the degree of freedom and the /// second dimension is the atom on which the charge resides. /// [1]: [J. Chem. Theory Comput. 2014, 10, 4801−4812](https://pubs.acs.org/doi/pdf/10.1021/ct500489d) - fn get_grad_dq(&self, atoms: &[Atom], grad_s: ArrayView3, p: ArrayView2) -> Array2 { + fn get_grad_dq(&self, atoms: &[Atom], s: ArrayView2, grad_s: ArrayView3, p: ArrayView2) -> Array2 { // get the shape of the derivative of S, it should be [f, n_orb, n_orb], where f = 3 * n_atoms let (f, n_orb, _): (usize, usize, usize) = grad_s.dim(); @@ -149,11 +151,15 @@ impl GroundStateGradient for Monomer { .into_shape([f, n_orb, n_orb]) .unwrap(); + // Compute W . S, and contract their last dimension + let w_s: Array2 = w.into_shape([f * n_orb, n_orb]).unwrap().dot(&s); + // compute P . S'; it is necessary to broadcast P into the shape of S' - let d_grad_s: Array3 = &grad_s * &p.broadcast([f, n_orb, n_orb]).unwrap(); + let d_grad_s: Array2 = grad_s_2d.dot(&p); + //println!("pair d_grad_s {}", w); // do the sum of both terms and sum over nu - let w_plus_ps: Array2 = (w + d_grad_s).sum_axis(Axis(2)); + let w_plus_ps: Array3 = (&w_s + &d_grad_s).into_shape([f, n_orb, n_orb]).unwrap(); // sum over mu where mu is on atom a let mut grad_dq: Array2 = Array2::zeros([f, self.n_atoms]); @@ -162,11 +168,12 @@ impl GroundStateGradient for Monomer { for _ in atom.valorbs.iter() { grad_dq .slice_mut(s![.., idx]) - .add_assign(&w_plus_ps.slice(s![.., mu])); + .add_assign(&w_plus_ps.slice(s![.., mu, mu])); mu += 1; } } + println!("grad dq {}", grad_dq); // Shape of returned Array: [f, n_atoms], f = 3 * n_atoms return grad_dq; } diff --git a/src/fmo/gradients/numerical.rs b/src/fmo/gradients/numerical.rs index 6601b15..c56679b 100644 --- a/src/fmo/gradients/numerical.rs +++ b/src/fmo/gradients/numerical.rs @@ -30,6 +30,7 @@ impl SuperSystem { } pub fn embedding_energy_wrapper(&mut self, geometry: Array1) -> f64 { + self.properties.reset(); self.update_xyz(geometry); self.prepare_scc(); let maxiter: usize = self.config.scf.scf_max_cycles; @@ -40,10 +41,10 @@ impl SuperSystem { } pub fn esd_energy_wrapper(&mut self, geometry: Array1) -> f64 { + self.properties.reset(); self.update_xyz(geometry); self.prepare_scc(); - let maxiter: usize = self.config.scf.scf_max_cycles; - let (_energy, _dq): (f64, Array1) = self.monomer_scc(maxiter); + //println!("{}", self.properties.gamma().unwrap()); self.esd_pair_energy() } @@ -71,5 +72,17 @@ impl SuperSystem { assert_deriv(self, SuperSystem::embedding_energy_wrapper, SuperSystem::embedding_gradient, self.get_xyz(), 0.01, 1e-6); } + pub fn test_esd_gradient(&mut self) { + self.prepare_scc(); + let maxiter: usize = self.config.scf.scf_max_cycles; + let (_energy, dq): (f64, Array1) = self.monomer_scc(maxiter); + + println!("ESD ENERGY {}", self.esd_pair_energy()); + self.properties.set_dq(dq); + let m_gradients: Array1 = self.monomer_gradients(); + + assert_deriv(self, SuperSystem::esd_energy_wrapper, SuperSystem::es_dimer_gradient, self.get_xyz(), 0.01, 1e-6); + } + } diff --git a/src/fmo/gradients/pair.rs b/src/fmo/gradients/pair.rs index a151774..68ad219 100644 --- a/src/fmo/gradients/pair.rs +++ b/src/fmo/gradients/pair.rs @@ -28,6 +28,8 @@ use rayon::iter::ParallelIterator; use rayon::prelude::IntoParallelRefMutIterator; use std::iter::FromIterator; use std::ops::{AddAssign, SubAssign}; +use crate::utils::array_helper::ToOwnedF; + impl GroundStateGradient for Pair { fn scc_gradient(&mut self, atoms: &[Atom]) -> Array1 { @@ -45,7 +47,9 @@ impl GroundStateGradient for Pair { // the derivatives of the charge (difference)s are computed at this point, since they depend // on the derivative of S and this is available here at no additional cost. let p: ArrayView2 = self.properties.p().unwrap(); - let grad_dq: Array2 = self.get_grad_dq(&atoms, grad_s.view(), p.view()); + let s: ArrayView2 = self.properties.s().unwrap(); + let grad_dq: Array2 = self.get_grad_dq(&atoms, s.view(), grad_s.view(), p.view()); + drop(s); self.properties.set_grad_dq(grad_dq); // and reshape them into a 2D array. the last two dimension (number of orbitals) are compressed @@ -120,7 +124,7 @@ impl GroundStateGradient for Pair { /// monomer. This means that the first dimension of the Array is the degree of freedom and the /// second dimension is the atom on which the charge resides. /// [1]: [J. Chem. Theory Comput. 2014, 10, 4801−4812](https://pubs.acs.org/doi/pdf/10.1021/ct500489d) - fn get_grad_dq(&self, atoms: &[Atom], grad_s: ArrayView3, p: ArrayView2) -> Array2 { + fn get_grad_dq(&self, atoms: &[Atom], s: ArrayView2, grad_s: ArrayView3, p: ArrayView2) -> Array2 { // get the shape of the derivative of S, it should be [f, n_orb, n_orb], where f = 3 * n_atoms let (f, n_orb, _): (usize, usize, usize) = grad_s.dim(); @@ -148,11 +152,15 @@ impl GroundStateGradient for Pair { .into_shape([f, n_orb, n_orb]) .unwrap(); + // Compute W . S, and contract their last dimension + let w_s: Array2 = w.into_shape([f * n_orb, n_orb]).unwrap().dot(&s); + // compute P . S'; it is necessary to broadcast P into the shape of S' - let d_grad_s: Array3 = &grad_s * &p.broadcast([f, n_orb, n_orb]).unwrap(); + let d_grad_s: Array2 = grad_s_2d.dot(&p); + //println!("pair d_grad_s {}", w); // do the sum of both terms and sum over nu - let w_plus_ps: Array2 = (w + d_grad_s).sum_axis(Axis(2)); + let w_plus_ps: Array3 = (&w_s + &d_grad_s).into_shape([f, n_orb, n_orb]).unwrap(); // sum over mu where mu is on atom a let mut grad_dq: Array2 = Array2::zeros([f, self.n_atoms]); @@ -161,7 +169,7 @@ impl GroundStateGradient for Pair { for _ in atom.valorbs.iter() { grad_dq .slice_mut(s![.., idx]) - .add_assign(&w_plus_ps.slice(s![.., mu])); + .add_assign(&w_plus_ps.slice(s![.., mu, mu])); mu += 1; } } diff --git a/src/fmo/scc/mod.rs b/src/fmo/scc/mod.rs index bfb7631..76d6475 100644 --- a/src/fmo/scc/mod.rs +++ b/src/fmo/scc/mod.rs @@ -44,6 +44,7 @@ impl RestrictedSCC for SuperSystem { fn prepare_scc(&mut self) { // prepare all individual monomers let atoms: &[Atom] = &self.atoms; + self.properties.set_gamma(gamma_atomwise(&self.gammafunction, &atoms, self.atoms.len())); self.monomers .par_iter_mut() .for_each(|mol: &mut Monomer| { @@ -59,6 +60,8 @@ impl RestrictedSCC for SuperSystem { let max_iter: usize = self.config.scf.scf_max_cycles; logging::fmo_scc_init(max_iter); + + // Assembling of the energy following Eq. 11 in // https://pubs.acs.org/doi/pdf/10.1021/ct500489d // E = sum_I^N E_I^ + sum_(I>J)^N ( E_(IJ) - E_I - E_J ) + sum_(I>J)^(N) DeltaE_(IJ)^V diff --git a/src/fmo/scc/monomer.rs b/src/fmo/scc/monomer.rs index 3cc1c79..5ee3cbd 100644 --- a/src/fmo/scc/monomer.rs +++ b/src/fmo/scc/monomer.rs @@ -28,6 +28,7 @@ impl Monomer { let atomic_numbers: Vec = atoms.iter().map(|atom| atom.number).collect(); self.properties.set_atomic_numbers(atomic_numbers); // get the gamma matrix + let gamma: Array2 = gamma_atomwise(&self.gammafunction, &atoms, self.n_atoms); // and save it as a `Property` self.properties.set_gamma(gamma); diff --git a/src/gradients/numerical.rs b/src/gradients/numerical.rs index dcc100c..131b37f 100644 --- a/src/gradients/numerical.rs +++ b/src/gradients/numerical.rs @@ -127,7 +127,7 @@ pub fn assert_deriv( // Ridder's method let (numerical_deriv, deriv_error): (f64, f64) = ridders_method(system, &function, origin.clone(), i, stepsize, con, safe, maxiter); let diff: f64 = (numerical_deriv - analytic_deriv).abs(); - let correct: bool = if diff >= deriv_error {false} else {true}; + let correct: bool = if diff >= deriv_error && diff > 1e-10 {false} else {true}; errors.push(correct); println!( diff --git a/src/initialization/properties.rs b/src/initialization/properties.rs index 355ae8f..b5284ac 100644 --- a/src/initialization/properties.rs +++ b/src/initialization/properties.rs @@ -288,6 +288,16 @@ impl Properties { } } + /// Returns a reference to the diagonal of the derivative of charge (differences) w.r.t. to the + /// degrees of freedom per atom. The first dimension is dof and the second one is the atom + /// where the charge resides. + pub fn grad_dq_diag(&self) -> Option> { + match self.get("grad_dq_diag") { + Some(value) => Some(value.as_array1().unwrap().view()), + _=> None, + } + } + /// Returns a reference to the differences of charges differences between the pair and the /// corresponding monomers per atom. pub fn delta_dq(&self) -> Option> { @@ -443,6 +453,9 @@ impl Properties { /// Set the derivative of the charge (differences) w.r.t. to degrees of freedom per atom pub fn set_grad_dq(&mut self, grad_dq: Array2) { self.set("grad_dq", Property::from(grad_dq))} + /// Set the diagonal of the derivative of the charge (differences) w.r.t. to degrees of freedom per atom + pub fn set_grad_dq_diag(&mut self, grad_dq_diag: Array1) { self.set("grad_dq_diag", Property::from(grad_dq_diag))} + /// Set the difference of charge differences between the pair dq and the dq's from the /// corresponding monomers per atom. pub fn set_delta_dq(&mut self, delta_dq: Array1) { diff --git a/src/main.rs b/src/main.rs index cbb8b88..7ec2e32 100644 --- a/src/main.rs +++ b/src/main.rs @@ -99,8 +99,9 @@ fn main() { let mut system = SuperSystem::from((frame, config)); //gamma_atomwise(&system.gammafunction, &system.atoms, system.atoms.len()); system.prepare_scc(); - let result = system.test_monomer_gradient(); + //let result = system.test_monomer_gradient(); let result = system.test_embedding_gradient(); + //let result = system.test_esd_gradient(); println!("Grad {:?}", system.ground_state_gradient()); info!("{}", timer);