/// Compute the dynamical matrix at a q-point.
///
/// The force constants do not need to contain data for rows outside the
/// designated cell. (but if they do, it won't hurt)
Compute the dynamical matrix at a q-point.
The force constants do not need to contain data for rows outside the
designated cell. (but if they do, it won't hurt)
Compute the dynamical matrix at a q - point . The force constants do not need to contain data for rows outside the designated cell . ( but if they do it won ' t hurt )
pub fn dynmat_at_cart_q(
        &self,
        super_coords: &Coords,
        qpoint_cart: V3,
        sc: &SupercellToken,
        masses: &[f64],
    ) -> DynamicalMatrix {
        assert_eq!(masses.len(), sc.num_primitive_atoms());
        let masses: &Indexed<PrimI, _> = Indexed::from_raw_ref(&masses[..]);
        let sc = SupercellWrapper::new(sc);
        let primitive_atoms = sc.atom_primitive_atoms();
        // Since rsp2 pays so much attention to image vectors (see FracBond) in places,
        // one might wonder why it is now doing brute-force searches for nearest images.
        //
        // Well... our current representation of force_sets doesn't contain image information.
        // All we know is that our supercell should be large enough that each site in the
        // designated cell only interacts with one image of any other site in the supercell.
        let image_finder = rsp2_structure::NearestImageFinder::new(super_coords.lattice()).unwrap();
        // Reducing the carts is a requirement for NearestImageFinder::shortest_images_cart_fast
        let reduced_carts = Indexed::<SuperI, _>::from_raw({
            let mut reduced = super_coords.clone();
            reduced.reduce_positions();
            reduced.to_carts()
        });
        let mut shortest_images_buf = vec![];
        let coo = self.0.to_coo();
        let iter = zip_eq!(&coo.row, &coo.col, &coo.val)
            // each column of the dynamical matrix sums over columns for images in
            // the force constants matrix, with phase factors.
            .map(|(&prim_r, &super_c, &mat)| {
                let super_r = sc.designated_super(prim_r);
                let prim_c = primitive_atoms[super_c];
                // mass-normalizing scale factor
                let scale = 1.0 / f64::sqrt(masses[prim_r] * masses[prim_c]);
                // vector from atom in primitive cell to atom in supercell.
                //
                // NOTE: Strictly speaking, rsp2 does not currently need to consider ties for
                // length. If necessary, we can optimize this by only getting one vector.
                //
                // I believe that the reason that phonopy must do this is because phonopy's
                // forces are dense; phonopy assumes that the equilibrium structure has zero
                // force, so the small forces that are present in the equilibrium structure
                // will also appear between distant, noninteracting sites in the force sets.
                // Considering ties thus helps cancel out some terms that would otherwise
                // violate symmetry.
                //
                // In rsp2, noninteracting sites in the supercell truly have zero force
                // between them. There should not be any ties for any of the `(r, c)` pairs
                // that appear in our sparse force constants.
                //
                // For now, however, we *do* find all shortest images, simply because I don't want
                // this to come back to bite me if things change in the future.
                // `_fast` is okay because the carts are reduced.
                image_finder.shortest_images_cart_fast(
                    &mut shortest_images_buf,
                    reduced_carts[super_c] - reduced_carts[super_r],
                    1e-4,
                );
                if shortest_images_buf.len() > 1 {
                    // If this warning is generated spuriously, then it probably means a potential
                    // was added to rsp2 that is incapable of producing sparse force sets.
                    //
                    // If that is indeed the case, the existing code should still handle it properly
                    // (though it hasn't been tested under these conditions), and either the
                    // potential should be fixed if possible or the warning should be removed.
                    warn_once!("\
                        Multiple shortest images found for a vector in the force constants! \
                        This could mean that your supercell is not large enough.\n\
                        (or, the force sets have somehow become dense!)\n\
                        Lattice: {:?}\n\
                        Vectors: {:?}\n\
                          Force: {:?}\n\
                    ", super_coords.lattice(), shortest_images_buf, mat);
                }
                let (phase_real, phase_imag) = {
                    shortest_images_buf.iter()
                        .map(|cart_diff| V3::dot(&qpoint_cart, cart_diff))
                        .map(|arg| 2.0 * std::f64::consts::PI * arg)
                        .map(|arg| (arg.cos(), arg.sin()))
                        .fold((0.0, 0.0), |(ar, ai), (br, bi)| (ar + br, ai + bi))
                };
                // NOTE: dividing by the multiplicity here is something that phonopy does.
                //       I'm not sure precisely *why* it does so. (simply summing over the images
                //       seems more natural than averaging over them).  But for our sanity,
                //       we'll do what they do.
                //
                //       (I don't think it matters much anyways in the end, because multiplicity
                //        should only ever exceed 1 for distant atoms (about half of the supercell),
                //        so the forces are small to begin with)
                //
                //       Notably, this also makes the graphene_denseforce_111 test succeed...
                //       though I am still suspicious as to how that test is even correct!
                let real = scale * phase_real * mat / shortest_images_buf.len() as f64;
                let imag = scale * phase_imag * mat / shortest_images_buf.len() as f64;
                ((prim_r, prim_c), Complex33(real, imag))
            });
        let matrix = {
            let (pos, val): (Vec<_>, Vec<_>) = iter.unzip();
            let (row, col) = pos.into_iter().unzip();
            let dim = (sc.raw.num_primitive_atoms(), sc.raw.num_primitive_atoms());
            if val != val {
                panic!("Dynamical matrix contains NaN!");
            }
            RawCoo { dim, val, row, col }.into_csr()
        };
        DynamicalMatrix(matrix)
    }
pub fn dynmat_at_cart_q ( & self , super_coords : & Coords , qpoint_cart : V3 , sc : & SupercellToken , masses : & [ f64 ] , ) -> DynamicalMatrix { assert_eq ! ( masses len ( ) sc num_primitive_atoms ( ) ) ; let masses : & Indexed < PrimI , _ > = Indexed :: from_raw_ref ( & masses [ .. ] ) ; let sc = SupercellWrapper :: new ( sc ) ; let primitive_atoms = sc . atom_primitive_atoms ( ) ; // Since rsp2 pays so much attention to image vectors (see FracBond) in places, // one might wonder why it is now doing brute-force searches for nearest images. // // Well... our current representation of force_sets doesn't contain image information. // All we know is that our supercell should be large enough that each site in the // designated cell only interacts with one image of any other site in the supercell. let image_finder = rsp2_structure :: NearestImageFinder :: new ( super_coords . lattice ( ) ) . unwrap ( ) ; // Reducing the carts is a requirement for NearestImageFinder::shortest_images_cart_fast let reduced_carts = Indexed :: < SuperI , _ > :: from_raw ( { let mut reduced = super_coords . clone ( ) ; reduced . reduce_positions ( ) ; reduced . to_carts ( ) } ) ; let mut shortest_images_buf = vec ! [ ] ; let coo = self . 0 . to_coo ( ) ; let iter = zip_eq ! ( coo row coo col coo val ) // each column of the dynamical matrix sums over columns for images in // the force constants matrix, with phase factors. . map ( | ( & prim_r , & super_c , & mat ) | { let super_r = sc . designated_super ( prim_r ) ; let prim_c = primitive_atoms [ super_c ] ; // mass-normalizing scale factor let scale = 1.0 / f64 :: sqrt ( masses [ prim_r ] * masses [ prim_c ] ) ; // vector from atom in primitive cell to atom in supercell. // // NOTE: Strictly speaking, rsp2 does not currently need to consider ties for // length. If necessary, we can optimize this by only getting one vector. // // I believe that the reason that phonopy must do this is because phonopy's // forces are dense; phonopy assumes that the equilibrium structure has zero // force, so the small forces that are present in the equilibrium structure // will also appear between distant, noninteracting sites in the force sets. // Considering ties thus helps cancel out some terms that would otherwise // violate symmetry. // // In rsp2, noninteracting sites in the supercell truly have zero force // between them. There should not be any ties for any of the `(r, c)` pairs // that appear in our sparse force constants. // // For now, however, we *do* find all shortest images, simply because I don't want // this to come back to bite me if things change in the future. // `_fast` is okay because the carts are reduced. image_finder . shortest_images_cart_fast ( & mut shortest_images_buf , reduced_carts [ super_c ] - reduced_carts [ super_r ] , 1e-4 , ) ; if shortest_images_buf . len ( ) > 1 { // If this warning is generated spuriously, then it probably means a potential // was added to rsp2 that is incapable of producing sparse force sets. // // If that is indeed the case, the existing code should still handle it properly // (though it hasn't been tested under these conditions), and either the // potential should be fixed if possible or the warning should be removed. warn_once ! ( " \
 \
 \n \
 \n \
 \n \
 \n \
 \n \
 " super_coords lattice ( ) shortest_images_buf mat ) ; } let ( phase_real , phase_imag ) = { shortest_images_buf . iter ( ) . map ( | cart_diff | V3 :: dot ( & qpoint_cart , cart_diff ) ) . map ( | arg | 2.0 * std :: f64 :: consts :: PI * arg ) . map ( | arg | ( arg . cos ( ) , arg . sin ( ) ) ) . fold ( ( 0.0 , 0.0 ) , | ( ar , ai ) , ( br , bi ) | ( ar + br , ai + bi ) ) } ; // NOTE: dividing by the multiplicity here is something that phonopy does. //       I'm not sure precisely *why* it does so. (simply summing over the images //       seems more natural than averaging over them).  But for our sanity, //       we'll do what they do. // //       (I don't think it matters much anyways in the end, because multiplicity //        should only ever exceed 1 for distant atoms (about half of the supercell), //        so the forces are small to begin with) // //       Notably, this also makes the graphene_denseforce_111 test succeed... //       though I am still suspicious as to how that test is even correct! let real = scale * phase_real * mat / shortest_images_buf . len ( ) as f64 ; let imag = scale * phase_imag * mat / shortest_images_buf . len ( ) as f64 ; ( ( prim_r , prim_c ) , Complex33 ( real , imag ) ) } ) ; let matrix = { let ( pos , val ) : ( Vec < _ > , Vec < _ > ) = iter . unzip ( ) ; let ( row , col ) = pos . into_iter ( ) . unzip ( ) ; let dim = ( sc . raw . num_primitive_atoms ( ) , sc . raw . num_primitive_atoms ( ) ) ; if val != val { panic ! ( " " ) ; } RawCoo { dim , val , row , col } . into_csr ( ) } ; DynamicalMatrix ( matrix ) }
Compute the dynamical matrix at a q-point.
Compute the dynamical matrix at a q - point .
&self
super_coords
&Coords
qpoint_cart
V3
sc
&SupercellToken
masses
&[f64]