/// Use symmetry to expand the number of displacement-force equations for
/// a row belonging to a symmetry star representative.
//
// Inspired by phonopy, we build the complete forces in such a way that
// ensures that the forces at each affected atom are described by the same
// displacements (in the same order). This way we only need to compute one
// pseudoinverse per symmetry star.
//
// I think this function vaguely corresponds to the first half of phonopy's
// `_solve_force_constants_svd` (though to be entirely honest, I just picked
// a sizable looking block of code and did `Extract Method`)
Use symmetry to expand the number of displacement-force equations for
a row belonging to a symmetry star representative.

Inspired by phonopy, we build the complete forces in such a way that
ensures that the forces at each affected atom are described by the same
displacements (in the same order). This way we only need to compute one
pseudoinverse per symmetry star.

I think this function vaguely corresponds to the first half of phonopy's
`_solve_force_constants_svd` (though to be entirely honest, I just picked
a sizable looking block of code and did `Extract Method`)
Use symmetry to expand the number of displacement - force equations for a row belonging to a symmetry star representative . Inspired by phonopy we build the complete forces in such a way that ensures that the forces at each affected atom are described by the same displacements ( in the same order ) . This way we only need to compute one pseudoinverse per symmetry star . I think this function vaguely corresponds to the first half of phonopy ' s ` _solve_force_constants_svd ` ( though to be entirely honest I just picked a sizable looking block of code and did ` Extract Method ` )
fn build_all_equations_for_representative_row(
        &self,
        // The row we are generating, which should be the representative atom
        // of a site-symmetry star. (needed to determine lattice point corrections
        // for each operator)
        displaced_prim: PrimI,
        // All displacements that displace this atom.
        disp_indices: &[DispI],
        // Operators that map `displaced_prim` to itself in the primitive structure.
        invariant_opers: &[OperI],
    ) -> (
        Indexed<EqnI, Vec<V3>>,
        BTreeMap<SuperI, Indexed<EqnI, Vec<V3>>>,
    ) {
        assert!(
            disp_indices.len() <= 6,
            "(BUG) got {} displacements. That's a lot! Are you sure these are all for the same atom?",
            disp_indices.len(),
        );
        let mut all_displacements = Indexed::<EqnI, Vec<V3>>::new();
        // Initially a BTreeMap<EqnI, _> is used for forces in case the set of affected atoms in
        // a force set is somehow not symmetric (in which case some rotations will have a different
        // set of affected atoms than others)
        let mut all_sparse_forces = BTreeMap::<SuperI, BTreeMap<EqnI, V3>>::new();
        for &oper in invariant_opers {
            // get a composite rotation + lattice-point translation that maps the
            // displaced supercell atom directly into itself (not just an image)
            let (rotated_prim, rotate_and_translate_atom) = self.get_corrected_rotate(oper, displaced_prim);
            assert_eq!(rotated_prim, displaced_prim, "(BUG) prim not invariant under supplied oper!");
            let rotate_vector = |v: V3| v * self.cart_rots[oper].t();
            for &disp in disp_indices {
                assert_eq!(self.displacements[disp].0, displaced_prim, "(BUG) disp for wrong atom");
                let eqn_i = all_displacements.push(rotate_vector(self.displacements[disp].1));
                for (&affected_atom, &cart_force) in &self.force_sets[disp] {
                    if cart_force != cart_force {
                        panic!("Force sets contain NaN!");
                    }
                    let new_affected_atom = rotate_and_translate_atom(affected_atom);
                    let new_cart_force = rotate_vector(cart_force);
                    if let Some(_) = {
                        all_sparse_forces
                            .entry(new_affected_atom)
                            .or_insert_with(BTreeMap::new)
                            .insert(eqn_i, new_cart_force)
                    } {
                        panic!("BUG! multiple affected atoms rotated to same one?!");
                    }
                }
            }
        }
        // done mutating
        let all_displacements = all_displacements;
        let all_sparse_forces = all_sparse_forces;
        // now that all equations have been built, we can "densify" the force matrices.
        // (with respect to equations; not the entire supercell, of course!)
        let all_forces = {
            all_sparse_forces.into_iter()
                .map(|(atom, force_map)| {
                    let mut force_vec = Indexed::from_elem(V3::zero(), &all_displacements);
                    for (eqn, v3) in force_map {
                        force_vec[eqn] = v3;
                    }
                    (atom, force_vec)
                }).collect()
        };
        (all_displacements, all_forces)
    }
fn build_all_equations_for_representative_row ( & self , // The row we are generating, which should be the representative atom // of a site-symmetry star. (needed to determine lattice point corrections // for each operator) displaced_prim : PrimI , // All displacements that displace this atom. disp_indices : & [ DispI ] , // Operators that map `displaced_prim` to itself in the primitive structure. invariant_opers : & [ OperI ] , ) -> ( Indexed < EqnI , Vec < V3 > > , BTreeMap < SuperI , Indexed < EqnI , Vec < V3 > > > , ) { assert ! ( disp_indices len ( ) 6 " " disp_indices len ( ) ) ; let mut all_displacements = Indexed :: < EqnI , Vec < V3 > > :: new ( ) ; // Initially a BTreeMap<EqnI, _> is used for forces in case the set of affected atoms in // a force set is somehow not symmetric (in which case some rotations will have a different // set of affected atoms than others) let mut all_sparse_forces = BTreeMap :: < SuperI , BTreeMap < EqnI , V3 > > :: new ( ) ; for & oper in invariant_opers { // get a composite rotation + lattice-point translation that maps the // displaced supercell atom directly into itself (not just an image) let ( rotated_prim , rotate_and_translate_atom ) = self . get_corrected_rotate ( oper , displaced_prim ) ; assert_eq ! ( rotated_prim displaced_prim " " ) ; let rotate_vector = | v : V3 | v * self . cart_rots [ oper ] . t ( ) ; for & disp in disp_indices { assert_eq ! ( self displacements [ disp ] 0 displaced_prim " " ) ; let eqn_i = all_displacements . push ( rotate_vector ( self . displacements [ disp ] . 1 ) ) ; for ( & affected_atom , & cart_force ) in & self . force_sets [ disp ] { if cart_force != cart_force { panic ! ( " " ) ; } let new_affected_atom = rotate_and_translate_atom ( affected_atom ) ; let new_cart_force = rotate_vector ( cart_force ) ; if let Some ( _ ) = { all_sparse_forces . entry ( new_affected_atom ) . or_insert_with ( BTreeMap :: new ) . insert ( eqn_i , new_cart_force ) } { panic ! ( " " ) ; } } } } // done mutating let all_displacements = all_displacements ; let all_sparse_forces = all_sparse_forces ; // now that all equations have been built, we can "densify" the force matrices. // (with respect to equations; not the entire supercell, of course!) let all_forces = { all_sparse_forces . into_iter ( ) . map ( | ( atom , force_map ) | { let mut force_vec = Indexed :: from_elem ( V3 :: zero ( ) , & all_displacements ) ; for ( eqn , v3 ) in force_map { force_vec [ eqn ] = v3 ; } ( atom , force_vec ) } ) . collect ( ) } ; ( all_displacements , all_forces ) }
Use symmetry to expand the number of displacement-force equations for
a row belonging to a symmetry star representative.
Use symmetry to expand the number of displacement - force equations for a row belonging to a symmetry star representative .
&self
displaced_prim
PrimI
disp_indices
&[DispI]
invariant_opers
&[OperI]