From 8c074a0d7d97d240820c5d389ab40c701b10c734 Mon Sep 17 00:00:00 2001 From: wolearyc Date: Thu, 8 Aug 2024 21:10:24 -0700 Subject: [PATCH] implemented ARTModel --- docs/coverage-badge.svg | 2 +- .../generated/ramannoodle.polarizability.rst | 8 + docs/tests-badge.svg | 2 +- ramannoodle/polarizability/art.py | 191 ++++++++++++++++ test/data/TiO2/known_art_spectrum.npz | Bin 0 -> 5714 bytes test/tests/test_art.py | 211 ++++++++++++++++++ test/tests/test_spectrum.py | 62 ++++- 7 files changed, 472 insertions(+), 4 deletions(-) create mode 100644 ramannoodle/polarizability/art.py create mode 100644 test/data/TiO2/known_art_spectrum.npz create mode 100644 test/tests/test_art.py diff --git a/docs/coverage-badge.svg b/docs/coverage-badge.svg index 8691f97..83b2e82 100644 --- a/docs/coverage-badge.svg +++ b/docs/coverage-badge.svg @@ -1 +1 @@ -coverage: 93.46%coverage93.46% +coverage: 93.22%coverage93.22% diff --git a/docs/source/generated/ramannoodle.polarizability.rst b/docs/source/generated/ramannoodle.polarizability.rst index 555192f..d27d7ab 100644 --- a/docs/source/generated/ramannoodle.polarizability.rst +++ b/docs/source/generated/ramannoodle.polarizability.rst @@ -4,6 +4,14 @@ ramannoodle.polarizability package Submodules ---------- +ramannoodle.polarizability.art module +------------------------------------- + +.. automodule:: ramannoodle.polarizability.art + :members: + :undoc-members: + :show-inheritance: + ramannoodle.polarizability.interpolation module ----------------------------------------------- diff --git a/docs/tests-badge.svg b/docs/tests-badge.svg index 9611b98..5ee71a0 100644 --- a/docs/tests-badge.svg +++ b/docs/tests-badge.svg @@ -1 +1 @@ -tests: 87tests87 +tests: 104tests104 diff --git a/ramannoodle/polarizability/art.py b/ramannoodle/polarizability/art.py new file mode 100644 index 0000000..e8c5a11 --- /dev/null +++ b/ramannoodle/polarizability/art.py @@ -0,0 +1,191 @@ +"""Polarizability model based on atomic Raman tensors (ARTs).""" + +from pathlib import Path + +import numpy as np +from numpy.typing import NDArray + +from .interpolation import InterpolationModel +from ..exceptions import ( + get_type_error, + get_shape_error, + verify_ndarray_shape, + InvalidDOFException, +) + + +class ARTModel(InterpolationModel): + """Polarizability model based on atomic Raman tensors (ARTs). + + This model is a flavor of InterpolationModel with added restrictions. The degrees + of freedom (DOF) in ARTModel are displacements of individual atoms. 1D + interpolations are used, with only two amplitudes allowed per DOF. Atomic + displacements can be in any direction, so long as all three directions for a given + atom are orthogonal. + + All assumptions from InterpolationModel apply. + + .. note:: + Degrees of freedom cannot (and should not) be added using `add_dof` and + `add_dof_from_files`. Instead, use `add_art` and `add_art_from_files`. + + .. warning:: + Due to the way InterpolationModel works, an ARTModel's predicted polarizability + tensor of the equilibrium structure may be offset slightly. This is a fixed + offset and therefore has no influence on Raman spectra calculations. + + Parameters + ---------- + structural_symmetry + equilibrium_polarizability + 2D array with shape (3,3) giving polarizability of system at equilibrium. This + would usually correspond to the minimum energy structure. + + """ + + def add_dof( # pylint: disable=too-many-arguments + self, + displacement: NDArray[np.float64], + amplitudes: NDArray[np.float64], + polarizabilities: NDArray[np.float64], + interpolation_order: int, + include_equilibrium_polarizability: bool = True, + ) -> None: + """Disable add_dof. + + :meta private: + """ + raise AttributeError("'ARTModel' object has no attribute 'add_dof'") + + def add_dof_from_files( + self, + filepaths: str | Path | list[str] | list[Path], + file_format: str, + interpolation_order: int, + ) -> None: + """Disable add_dof_from_files. + + :meta private: + """ + raise AttributeError("'ARTModel' object has no attribute 'add_dof_from_files'") + + def add_art( + self, + atom_index: int, + direction: NDArray[np.float64], + amplitudes: NDArray[np.float64], + polarizabilities: NDArray[np.float64], + ) -> None: + """Add an atomic Raman tensor (ART). + + Specification of an ART requires an atom index, displacement direction, + displacement amplitudes, and corresponding known polarizabilities for each + amplitude. Alongside ART's related to the specified ART by symmetry are added + automatically. + + Parameters + ---------- + atom_index + Index of atom, consistent with the StructuralSymmetry used to initialize + the model. + direction + 1D array with shape (3,). Must be orthogonal to all previously added ARTs + belonging to specified atom. + amplitudes + 1D array of length 1 or 2 containing amplitudes in angstroms. Duplicate + amplitudes are not allowed, including symmetrically equivalent + amplitudes. + polarizabilities + 3D array with shape (1,3,3) or (2,3,3) containing known polarizabilities for + each amplitude. + + Raises + ------ + InvalidDOFException + Provided ART was invalid. + + """ + if not isinstance(atom_index, int): + raise get_type_error("atom_index", atom_index, "int") + try: + if amplitudes.shape not in [(1,), (2,)]: + raise get_shape_error("amplitudes", amplitudes, "(1,) or (2,)") + except AttributeError as exc: + raise get_type_error("amplitudes", amplitudes, "ndarray") from exc + verify_ndarray_shape( + "polarizabilities", polarizabilities, (amplitudes.size, 3, 3) + ) + + displacement = self._structural_symmetry.get_fractional_positions() * 0 + try: + displacement[atom_index] = direction / np.linalg.norm(direction * 10.0) + except TypeError as exc: + raise get_type_error("direction", direction, "ndarray") from exc + except ValueError as exc: + raise get_shape_error("direction", direction, "(3,)") from exc + + super().add_dof( + displacement, + amplitudes, + polarizabilities, + 1, + include_equilibrium_polarizability=False, + ) + + def add_art_from_files( + self, + filepaths: str | Path | list[str] | list[Path], + file_format: str, + ) -> None: + """Add a atomic Raman tensor (ART) from file(s). + + Required directions, amplitudes, and polarizabilities are automatically + determined from provided files. Files should be chosen wisely such that the + resulting ARTs are valid under the restrictions set by `add_art`. + + Parameters + ---------- + filepaths + file_format + supports: "outcar" + + Raises + ------ + FileNotFoundError + File could not be found. + InvalidDOFException + ART assembled from supplied files was invalid (see get_art) + + """ + displacements, amplitudes, polarizabilities = super()._read_dof( + filepaths, file_format + ) + # Check whether only one atom is displaced. + _displacement = np.copy(displacements[0]) + atom_index = int(np.argmax(np.sum(_displacement**2, axis=1))) + _displacement[atom_index] = np.zeros(3) + if not np.isclose(_displacement, 0.0, atol=1e-6).all(): + raise InvalidDOFException("multiple atoms displaced simultaneously") + + basis_vectors_to_add, interpolation_xs, interpolation_ys = super()._get_dof( + displacements[0], amplitudes, polarizabilities, False + ) + + num_amplitudes = len(interpolation_xs[0]) + if num_amplitudes != 2: + raise InvalidDOFException( + f"wrong number of amplitudes: {num_amplitudes} != 2" + ) + + super()._construct_and_add_interpolations( + basis_vectors_to_add, interpolation_xs, interpolation_ys, 1 + ) + + def get_status_dict(self) -> dict[str, int]: + """Return information on model status. + + Relevant information includes which ARTs are required for each atom and whether + or not these ARTs have been specified. + + """ + raise NotImplementedError diff --git a/test/data/TiO2/known_art_spectrum.npz b/test/data/TiO2/known_art_spectrum.npz new file mode 100644 index 0000000000000000000000000000000000000000..0fc48c09ca9279f63340af659da24cb7e4934358 GIT binary patch literal 5714 zcmd7WX*8Ex+W_!CG7}}!k%J_ZA!F!V+fIWLg)(QJm3fvaQ#2AKq)9~+nT0YO86)#p zDf3k3x%^+P=dsrNob}va`t)3T-S=hk2j1xp`)Q`Ktb{%ofo#UcDHmBR^k`_%SJ(1 zj9=Ks+0Dbv+{x70&C2@M`VsSU?$+z;?sn!b*6VW-8EIKDQ8E4i{@)IUfu#q{k;YM& zNVfLMPI{IoOys55RbGasNc_9T1UC&j+C^d_r43beWK<(Dkx@BUl{0H2Fp=b}^D~X* z5tzs}$a+0iABs1mnJZ5o7YxN4-W&9FS>y=C8(L)*kgKOcFp<7>ZP&kzgCja|m5YWAmk{kFmp z16F;TJTZ|f2}feg6hRH?7lJ#}Ym`H0b?@&8CXG|o6RQzd2cPC8b)4h4DR$m#8Ltl9@b8!ysgdH;k1V6Zm6~1>UgQ zd~Z;Au{kF4#;w|(7*%sjOu{$CH%RIQzg&jnu6p$$?N8S#a2v0V*Kg?n4V!ldsVisUhq881P4PF7QkNfa4DE-8 z7y*heg&ok%KNIEb*$81J+KRV7*20b}hHVBsT&SdY`=9CME5O^UZxmR$2$ESRgSbT} z!R&k%W9YuO%At+3YZundX3A1l$=r+zCUy9Zxy1n0R{Dnq##KPCA9wxTf)(8Hz5QuKbpXT0}y0+xGSrw6_Trs{}5uYf-xaW1_qxCa83Tl zVd5evnu>7_3x>fc-4~e^-VA6@aZ(y*u!BpZkI06ihoF*_A}y$YAGrD6+x1Qy1kc-Z zZ#(Hy;d>5W`D%AA&~L4!(E4~4Y_8|$_?8#Kqg|avH*CXUP?YtZ)y0QEElPi9C0r4O zw)ExiY0XuxCEqPv=%7T+^xbrG{DR0!F?~R*uNb-vqOu8i{scx^ZQxEzuX4As}U8G=8+raJz?7v)1BRW~%;%$Eh3(5OsN9FH;jZOk5vtDfx_Y z1dlg&Uq6r(35neOY?VDPLD_P`v{jA^G44ySx&EySs`YDq@11-D!RjVy{I6%=Ih%;U z1D*~@+Yzr7mo@-L_i!@Ua+iTFLp#o8Z-zdWpeM)F;-OtGqt-^h6gZFL;?~X0Aj!zf zqaQvE6)zBdT*D-I><=(!v2F!*mEEpB`xhYH?1z9x=uQ39e__jKIOlbCc9Zll>@H(d)M@5L zjK@1eI}2tZNbFRzbk-EGMwljjFr-H(+%4)C#kZrDo{sH7CIZNZwel1zA01-5>FSbu zWiyhf;^MHCp+xKdV|Z+A6Dkk8<>*zE3pe?(TkrjH5TM9$|FWYKhK`hSoLL=))W}x$ zQ2h=VVG8>0<5v#ztM{jT%csF}>7;?{90ht9J+Q}?oe3pPYS?mInS&O`tW5oFEXWW^ z_NMLEM6bHSslAQak(4Og$9)KmYoe^6WoTjGSbpV;yOV%d{-zq$7uGWf+;O zjLq*rwnq0Y|Lx6CBR&mlhF%aD6LzsjF#~S(%t}{FHi2-2)$ow&JQ#S)TKBFk0OOUJ zm5IeMXt}i<9do06^Le>o*<4CqP$-8hm$UL`Vv9gX>u`+qHVaUu zshDi(Nd^~n(q#KcJy_H+*TrgYL2h-yQ;}yEfzxcHO|p}O`gy+=+xKyxdXdD1=qCq} znk%z@*i%z*&ZO*`J8%l@ZS8jtZ&L#|!LZa^iFjZ&4C&AbI|WaqBgvx{mO!1EC2+LW z72-dzKi+rt0*H(397-Q<09o&7F43BPSbIM{93(#q8Rp&Uhx58YQYxx^Y{z;YbMper z9raKb);F`8x)M6}nM|c~tig^PLk&@x3r;M9xOA4nbmhr}MP5;bn*BW7H=ca#b#%`SCNkIusPsNRD)pT0ui1dH!?cl!of2-4w8mEcRCI( z!OoDha2{4Uq;l0W_y9XU;yLYjN6%LYh4;-5r0dp#Aaj?#_Pz>GFquCZ>P(GH7EhWq z(S|{6fx?A2syJx7z1B=o6byxR&M9o;%Mh@I++N`D1*(pu|3zuH0=wSm{V-MHLk`KY z+x#OKk)}ghR^fhDX^iUk-8z?#$KSwCh zl3OA}trOO=ejkQq z;FQ-(>=ctgLBp!uy;@Xg{brUs-M0W@@%$x$7Y5)(7Y)_f{6RRuL|6LCcOFdB=CndB z`@z{qC6W76BmB*}R*_3n4n{rS9G5>HP<9gY`Fm;Irh*5I#>< z_*QIREHSJZ%Nla)srnv+bLTkJhYJJMfX&3hX@~JF`9-aLh+{+U>+19cMCvmLG zt{_@ezZ0}hzN15ne>|BJ+j|hn+3h&|r=$@waU4!Iy}TEtHt^dsA3hd^nWGDV`P5$v=LqIL*ZK#U&iG@nE_ z@aTGX4K~GrOe_}g&{2lnteQ^gKVHE1ErK@=xVk~ziTEjb?+NI6egAL%TPx5LW0*O{>WfU;L~#*QRqmqta4aqAljrI7OGxuh9tzRt*4935$AjZ zo-C6=E1G93RZ>NjqchVTJ)FAH=cUksmgPaDd7Oe$oa#S&xWAgSzaEx*@YcEo&a9rumhVa((!%Th9J#_J0@j5##^> literal 0 HcmV?d00001 diff --git a/test/tests/test_art.py b/test/tests/test_art.py new file mode 100644 index 0000000..d749ff4 --- /dev/null +++ b/test/tests/test_art.py @@ -0,0 +1,211 @@ +"""Testing for ARTModel.""" + +from typing import Type + +import numpy as np +from numpy.typing import NDArray + +import pytest + +from ramannoodle.polarizability.art import ARTModel +from ramannoodle.exceptions import InvalidDOFException +from ramannoodle.symmetry import StructuralSymmetry + +# pylint: disable=protected-access +# pylint: disable=too-many-arguments + + +@pytest.mark.parametrize( + "outcar_symmetry_fixture,atom_index, direction, amplitudes, known_dof_added", + [ + ( + "test/data/STO_RATTLED_OUTCAR", + 0, + np.array([1, 0, 0]), + np.array([-0.01, 0.01]), + 1, + ), + ("test/data/TiO2/phonons_OUTCAR", 0, np.array([1, 0, 0]), np.array([0.01]), 72), + ], + indirect=["outcar_symmetry_fixture"], +) +def test_add_art( + outcar_symmetry_fixture: StructuralSymmetry, + atom_index: int, + direction: NDArray[np.float64], + amplitudes: NDArray[np.float64], + known_dof_added: int, +) -> None: + """Test add_art (normal).""" + symmetry = outcar_symmetry_fixture + model = ARTModel(symmetry, np.zeros((3, 3))) + model.add_art(atom_index, direction, amplitudes, np.zeros((amplitudes.size, 3, 3))) + assert len(model._cartesian_basis_vectors) == known_dof_added + assert np.isclose(np.linalg.norm(model._cartesian_basis_vectors[0]), 1) + + +@pytest.mark.parametrize( + "outcar_symmetry_fixture,atom_indexes,directions,amplitudes,polarizabilities," + "exception_type,in_reason", + [ + ( + "test/data/STO_RATTLED_OUTCAR", + [0], + np.array([[1, 0, 0]]), + np.array([0.1]), + np.zeros((1, 3, 3)), + InvalidDOFException, + "insufficient points", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0], + np.array([[1, 0, 0]]), + np.array([-0.1, 0.1]), + np.zeros((2, 4, 3)), + ValueError, + "polarizabilities has wrong shape", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0], + np.array([[1, 0, 0]]), + np.array([-0.1, 0.1]), + np.zeros((4, 3, 3)), + ValueError, + "polarizabilities has wrong shape", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0, 0], + np.array([[1, 0, 0], [1, 0, 0]]), + np.array([-0.1, 0.1]), + np.zeros((2, 3, 3)), + InvalidDOFException, + "not orthogonal", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0, 0], + np.array([[1, 0, 0], [1, 0, 0]]), + [-0.1, 0.1], + np.zeros((2, 3, 3)), + TypeError, + "amplitudes should have type ndarray, not list", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0, 0], + np.array([[1, 0, 0], [1, 0, 0]]), + np.array([-0.1, 0.1]), + list(np.zeros((2, 3, 3))), + TypeError, + "polarizabilities should have type ndarray, not list", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [[0, 0]], + np.array([[1, 0, 0], [1, 0, 0]]), + np.array([-0.1, 0.1]), + np.zeros((2, 3, 3)), + TypeError, + "atom_index should have type int, not list", + ), + ( + "test/data/STO_RATTLED_OUTCAR", + [0], + np.array([[1, 0, 0], [1, 0, 0]]), + np.array([0.1, 0.1]), + np.zeros((2, 3, 3)), + InvalidDOFException, + "due to symmetry, amplitude 0.1 should not be specified", + ), + ], + indirect=["outcar_symmetry_fixture"], +) +def test_add_art_exception( + outcar_symmetry_fixture: StructuralSymmetry, + atom_indexes: list[int], + directions: NDArray[np.float64], + amplitudes: NDArray[np.float64], + polarizabilities: NDArray[np.float64], + exception_type: Type[Exception], + in_reason: str, +) -> None: + """Test add_art (exception).""" + symmetry = outcar_symmetry_fixture + model = ARTModel(symmetry, np.zeros((3, 3))) + with pytest.raises(exception_type) as error: + for atom_index, direction in zip(atom_indexes, directions): + model.add_art(atom_index, direction, amplitudes, polarizabilities) + + assert in_reason in str(error.value) + + +@pytest.mark.parametrize( + "outcar_symmetry_fixture,outcar_files,exception_type,in_reason", + [ + ( + "test/data/STO_RATTLED_OUTCAR", + ["test/data/TiO2/Ti5_0.1x_eps_OUTCAR"], + InvalidDOFException, + "incompatible outcar", + ), + ( + "test/data/TiO2/phonons_OUTCAR", + [ + "test/data/TiO2/Ti5_0.1x_eps_OUTCAR", + "test/data/TiO2/Ti5_0.2x_eps_OUTCAR", + ], + InvalidDOFException, + "wrong number of amplitudes: 4 != 2", + ), + ( + "test/data/TiO2/phonons_OUTCAR", + [ + "test/data/TiO2/Ti5_0.1x_eps_OUTCAR", + "test/data/TiO2/Ti5_m0.1x_eps_OUTCAR", + ], + InvalidDOFException, + "wrong number of amplitudes: 4 != 2", + ), + ( + "test/data/TiO2/phonons_OUTCAR", + [ + "this_outcar_does_not_exist", + ], + FileNotFoundError, + "No such file or directory", + ), + ( + "test/data/TiO2/phonons_OUTCAR", + [ + "test/data/TiO2/Ti5_0.1x_eps_OUTCAR", + "test/data/TiO2/Ti5_0.1y_eps_OUTCAR", + ], + InvalidDOFException, + "is not collinear", + ), + ( + "test/data/TiO2/phonons_OUTCAR", + [ + "test/data/TiO2/O43_0.1z_eps_OUTCAR", + ], + InvalidDOFException, + "wrong number of amplitudes: 1 != 2", + ), + ], + indirect=["outcar_symmetry_fixture"], +) +def test_add_art_from_files_exception( + outcar_symmetry_fixture: StructuralSymmetry, + outcar_files: list[str], + exception_type: Type[Exception], + in_reason: str, +) -> None: + """Test add_dof_from_files (exception).""" + symmetry = outcar_symmetry_fixture + model = ARTModel(symmetry, np.zeros((3, 3))) + with pytest.raises(exception_type) as error: + model.add_art_from_files(outcar_files, "outcar") + assert in_reason in str(error.value) diff --git a/test/tests/test_spectrum.py b/test/tests/test_spectrum.py index 125ba48..0e16126 100644 --- a/test/tests/test_spectrum.py +++ b/test/tests/test_spectrum.py @@ -9,6 +9,7 @@ import pytest from ramannoodle.polarizability.interpolation import InterpolationModel +from ramannoodle.polarizability.art import ARTModel from ramannoodle.symmetry import StructuralSymmetry from ramannoodle import io from ramannoodle.spectrum.spectrum_utils import ( @@ -64,12 +65,12 @@ def _validate_polarizabilities(model: InterpolationModel, data_directory: str) - ], indirect=["outcar_symmetry_fixture"], ) -def test_spectrum( +def test_interpolation_spectrum( outcar_symmetry_fixture: StructuralSymmetry, data_directory: str, dof_eps_outcars: list[str], ) -> None: - """Test a full spectrum calculation.""" + """Test a full spectrum calculation using InterpolationModel.""" # Setup model symmetry = outcar_symmetry_fixture _, polarizability = io.read_positions_and_polarizability( @@ -105,6 +106,63 @@ def test_spectrum( assert np.isclose(intensities, known_intensities, atol=1e-4).all() +@pytest.mark.parametrize( + "outcar_symmetry_fixture,data_directory,dof_eps_outcars", + [ + ( + "test/data/TiO2/phonons_OUTCAR", + "test/data/TiO2/", + [ + ["Ti5_0.1z_eps_OUTCAR"], + ["Ti5_0.1x_eps_OUTCAR"], + [ + "O43_0.1z_eps_OUTCAR", + "O43_m0.1z_eps_OUTCAR", + ], + ["O43_0.1x_eps_OUTCAR"], + ["O43_0.1y_eps_OUTCAR"], + ], + ), + ], + indirect=["outcar_symmetry_fixture"], +) +def test_art_spectrum( + outcar_symmetry_fixture: StructuralSymmetry, + data_directory: str, + dof_eps_outcars: list[str], +) -> None: + """Test a full spectrum calculation using ARTModel.""" + # Setup model + symmetry = outcar_symmetry_fixture + _, polarizability = io.read_positions_and_polarizability( + f"{data_directory}/ref_eps_OUTCAR", file_format="outcar" + ) + model = ARTModel(symmetry, polarizability) + for outcar_names in dof_eps_outcars: + model.add_art_from_files( + [f"{data_directory}/{name}" for name in outcar_names], file_format="outcar" + ) + + # Spectrum test + with np.load(f"{data_directory}/known_art_spectrum.npz") as known_spectrum: + phonons = io.read_phonons( + f"{data_directory}/phonons_OUTCAR", file_format="outcar" + ) + spectrum = phonons.get_raman_spectrum(model) + wavenumbers, intensities = spectrum.measure( + laser_correction=True, + laser_wavelength=532, + bose_einstein_correction=True, + temperature=300, + ) + + known_wavenumbers = known_spectrum["wavenumbers"] + known_intensities = known_spectrum["intensities"] + + assert np.isclose(wavenumbers, known_wavenumbers).all() + assert np.isclose(intensities, known_intensities, atol=1e-4).all() + + @pytest.mark.parametrize( "spectrum_path,known_gaussian_spectrum_path,known_lorentzian_spectrum_path", [