Molecular Dynamics Simulation of Poly(Ether Ether Ketone) (PEEK) Polymer to Analyze Intermolecular Ordering by Low Wavenumber Raman Spectroscopy and X-ray Diffraction

Poly(ether ether ketone) (PEEK) is an important engineering plastic and evaluation of its local crystallinity in composites is critical for producing strong and reliable mechanical parts. Low wavenumber Raman spectroscopy and X-ray diffraction are promising techniques for the analysis of crystal ordering but a detailed understanding of the spectra has not been established. Here, we use molecular dynamics combined with a newly developed approximation to simulate local vibrational features to understand the effect of intermolecular ordering in the Raman spectra. We found that intermolecular ordering does affect the low wavenumber Raman spectra and the X-ray diffraction as observed in the experiment. Raman spectroscopy of intermolecular vibration modes is a promising technique to evaluate the local crystallinity of PEEK and other engineering plastics, and the present technique offers an estimation without requiring heavy computational resources.


Introduction
Poly(ether ether ketone) (PEEK) is a high-performance thermoplastic and crystalline polymer with excellent mechanical properties, inertness in a harsh chemical environment, and thermal durability [1][2][3]. Reinforced composite materials of PEEK with inorganic or carbon fiber fillers have attracted attention as structural materials for aircraft and automobiles because of their light weight and excellent properties. It is established that the degree of crystallinity greatly affects the mechanical properties of crystalline polymers including PEEK. Since it has been reported that the crystallinity of fiber-reinforced polymers may be spatially non-uniform [4,5], it is important to evaluate local crystallinity less than the diameter of the filler fibers. Since the size of the fiber fillers is on the scale of several µm, evaluation of crystallinity in the sub-µm region is important.
There are various techniques for the evaluation of the crystallinity of polymers such as wide-angle X-ray diffraction (WAXRD), small-angle X-ray diffraction (SAXS), differential scanning calorimetry (DSC), Fourier transform infrared (FTIR) spectroscopy, and Raman microscopy [6,7]. Among them, Raman microscopy, which takes Raman spectra under the focus of an objective lens, is probably the only readily available technique [8]. X-ray diffraction microscopy using a synchrotron light source [9] or cryogenic transmission electron microscopy [10] is now possible, but the location and machine time for the measurement are very limited and preparation of the samples in very thin slabs is troublesome. Raman microscopy can be measured in the reflection geometry with the high resolution of optical microscopy that can easily reach sub-µm and possibly measure the local vibrational spectrum.

Extraction of Approximate Raman Signal from MD Simulation
The probability of Raman scattering is described by the Raman polarizability tensor, which is simplified by Placzek approximation [20,21] to Equation (1): in the case of Stokes Raman scattering. Here, the lefthand side of the equation is the Raman polarizability tensor corresponding to the vibrational transition from v n to v m , and n and j denote the ground state of electrons and excited states, respectively. ω jn is the energy difference between states n and j in angular frequency and ω the angular frequency of the incoming electromagnetic wave. µ α and µ β are dipole moment operators in the direction of α and β. Calculation of this quantity requires the calculation of precise electronic states of the molecule and its dependence on the position of the nuclei vibrating from the ground state. It becomes computationally very demanding in the case of polymer aggregates because of the many degrees of structural freedom as mentioned in the Introduction. Also, accurate treatment of van der Waals interactions becomes computationally very demanding but necessary for the estimation of intermolecular vibration. Our approach proposed here is to use an approximate local valence electron density. The motion of valence electrons by the external electric field is the cause of the dynamic dipole moment associated with the transition n|µ α |j and j|µ α |n . The fluctuation of the local valence electron density at the frequency of a vibration mode v m can be assumed proportional to them. Therefore, we can approximate the Raman intensity at vibrational frequency ν as proportional to the squared oscillation of valence electron density ρ caused by vibration: Since it is established that Raman scattering can occur in a single molecule [22], we can justify a crude approximation to neglect the spatial phase of the electric field of the incoming photons and consider the scattering locally at position x in the material: There will be various ways to calculate the local valence electron density. In this paper, we use a very simple approximation as follows: where ρ (x) is the valence electron density at position x and N i and R i are the valence electron number and characteristic atom radius of atom i, respectively. r i is the position of atom i. Here we used values for N i the nominal valence electron count of each atom (C 4, O 6, H 1), and for R i covalent radii were used (C 0.73 Å, O 0.66 Å, H 0.32 Å) [23]. In the following, we show the Fourier transform of ρ, not ρ 2 , to show the small features in the figures. The conclusion does not change if we use ρ 2 . This type of simplification of molecules has recently been used in machine learning for protein engineering and has been proven to work to extract the microchemical environment of biological polymers [24,25].
We must admit that the crude approximation in Equation (3) neglects the shape and symmetry of the wave functions of the molecule or its fragment, and the "spectrum" obtained by the Fourier transform of Equation (4) contains both infrared absorption and Raman scattering. The merit of doing this is the short computational time for large systems. We expect the usage of the present method is calculating large systems to find the local environment that causes a certain Raman signal (or infrared absorption). Then we can proceed with the first-principle calculation of extracted local structures to accurately evaluate it.

Computation Details
Molecular vibrations of single-strand PEEK oligomers were calculated with DFT calculation using Gaussian16 [26] at B3LYP/6-31G level of theory. The MD simulation was performed using LAMMPS [27] with generalized AMBER FF (GAFF) [28]. The atomic charges for the calculation of Coulomb interactions were calculated by Gaussian16 of a PEEK oligomer at B3LYP/6-31G level. The initial structure of the polymer based on the reported crystal structure [29] was prepared using the moltemplate utility [30]. There are several crystal structure parameters for PEEK [31][32][33][34][35][36][37][38] but they agree with each other only with a small deviation of the lattice constant and atom positions [6]. The results were visualized with OVITO [39] and VESTA [40]. The atomic trajectory data, including Fourier transform and smoothing spline fitting, were analyzed using programs written with Python 3.7. The vibration mode at frequency ν (Hz) is expressed in the unit of wavenumber ν (cm −1 ) throughout this paper using where c is the light speed.

Experiments
To compare the calculation with the experiments, WAXRD patterns and Raman spectra were measured on the samples annealed at various temperatures up to 300 • C for 2 h. The sample material was a quenched PEEK thin film with a thickness of 50 mm (Victrex APTIV film, grade 1000). Details of the experiment are described in [7]. WAXRD was measured by using Rigaku Ultima IV (rotating anode, 40 kV 30 mA) operating with CuKa X-ray (wavelength 1.5406 Å). Raman spectroscopy was measured using a Horiba Jobin Yvon LablamHR spectrometer with a Nd:YAG laser (wavelength 1064 nm, 60 mW, spot size 1 µm) with a 50× magnification microscope. Signals from a 10 s measurement were accumulated four times to confirm no damage to the sample during the measurement.

Experimental Results
Before presenting computational analysis, the experimental findings are reviewed briefly. The experimental data have been reported previously [7] but here we newly processed the data in a form suitable for comparison with the computational results. Figure 1 shows Raman spectra of PEEK after annealing at various temperatures for two hours. It can be noticed that there are two peaks in this low wavenumber region at 97 cm −1 and 135 cm −1 , and the intensity of the 135 cm −1 peak increases whereas that at 97 cm −1 does not change significantly. Figure 2 shows WAXRD patterns (inset) and the peak ratios of the same samples. There are four distinct peaks (Peak A 2θ = 18.9 • , B 20.7 • , C 22.8 • , D 28.8 • ), and the peak ratio changes by annealing at different temperatures. These features in Raman and WAXRD can be used for the evaluation of the crystallinity of PEEK, but the mechanism of the change has not been understood so far. In the following sections, we use computational methods to understand these features.
with a small deviation of the lattice constant and atom positions [6]. The results w ualized with OVITO [39] and VESTA [40]. The atomic trajectory data, including transform and smoothing spline fitting, were analyzed using programs written thon 3.7. The vibration mode at frequency ν (Hz) is expressed in the unit of wave � (cm -1 ) throughout this paper using where c is the light speed.

Experiments
To compare the calculation with the experiments, WAXRD patterns and Ram tra were measured on the samples annealed at various temperatures up to 300 °C The sample material was a quenched PEEK thin film with a thickness of 50 mm APTIV film, grade 1000). Details of the experiment are described in [7]. WAX measured by using Rigaku Ultima IV (rotating anode, 40 kV 30 mA) operating wi X-ray (wavelength 1.5406 Å). Raman spectroscopy was measured using a Hori Yvon LablamHR spectrometer with a Nd:YAG laser (wavelength 1064 nm, 60 m size 1 μm) with a 50× magnification microscope. Signals from a 10 s measureme accumulated four times to confirm no damage to the sample during the measure

Experimental Results
Before presenting computational analysis, the experimental findings are r briefly. The experimental data have been reported previously [7] but here we ne cessed the data in a form suitable for comparison with the computational results Figure 1 shows Raman spectra of PEEK after annealing at various tempera two hours. It can be noticed that there are two peaks in this low wavenumber reg cm -1 and 135 cm -1 , and the intensity of the 135 cm -1 peak increases whereas that a does not change significantly. Figure 2 shows WAXRD patterns (inset) and the pe of the same samples. There are four distinct peaks (Peak A 2θ =18.9 °, B 20.7 °, C 28.8 °), and the peak ratio changes by annealing at different temperatures. These in Raman and WAXRD can be used for the evaluation of the crystallinity of PEEK mechanism of the change has not been understood so far. In the following sect use computational methods to understand these features.

Calculated Raman Spectroscopy of a Single Molecule PEEK Chain
In order to understand the vibration of a PEEK molecule, DFT (b3lyp/6-31G) optimization and calculation of Raman intensity of each vibrational mode of phenyl-group-terminated PEEK oligomers with different lengths were performed using the Gaussian 16 package. The results are shown in Figures 3 and 4. Figure 3 shows the Raman spectrum derived from the Raman intensity of each mode with 5 cm -1 broadening. The structures of the molecules are written in the right panel, in which "Ph", "-O-" and "-CO-" denote phenyl group (C6H5-or -C6H4-), ether, and ketone, respectively. The peaks around 85-98 cm -1 (labeled as "(i)") and 130-145 cm -1 ("(ii)") are those corresponding to the experimental observation (97 cm -1 and 135 cm -1 ). The scaling factor for vibrational analysis using DFT (b3lyp functional) is 0.961 (only ~4% difference) [41] and can be neglected because the shift of the modes due to the chain elongation is larger than 4%.

Calculated Raman Spectroscopy of a Single Molecule PEEK Chain
In order to understand the vibration of a PEEK molecule, DFT (b3lyp/6-31G) optimization and calculation of Raman intensity of each vibrational mode of phenyl-groupterminated PEEK oligomers with different lengths were performed using the Gaussian 16 package. The results are shown in Figures 3 and 4. Figure 3 shows the Raman spectrum derived from the Raman intensity of each mode with 5 cm −1 broadening. The structures of the molecules are written in the right panel, in which "Ph", "-O-" and "-CO-" denote phenyl group (C 6 H 5 -or -C 6 H 4 -), ether, and ketone, respectively. The peaks around 85-98 cm −1 (labeled as "(i)") and 130-145 cm −1 ("(ii)") are those corresponding to the experimental observation (97 cm −1 and 135 cm −1 ). The scaling factor for vibrational analysis using DFT (b3lyp functional) is 0.961 (only~4% difference) [41] and can be neglected because the shift of the modes due to the chain elongation is larger than 4%.

Calculated Raman Spectroscopy of a Single Molecule PEEK Chain
In order to understand the vibration of a PEEK molecule, DFT (b3lyp/6-31G) optimization and calculation of Raman intensity of each vibrational mode of phenyl-group-terminated PEEK oligomers with different lengths were performed using the Gaussian 16 package. The results are shown in Figures 3 and 4. Figure 3 shows the Raman spectrum derived from the Raman intensity of each mode with 5 cm -1 broadening. The structures of the molecules are written in the right panel, in which "Ph", "-O-" and "-CO-" denote phenyl group (C6H5-or -C6H4-), ether, and ketone, respectively. The peaks around 85-98 cm -1 (labeled as "(i)") and 130-145 cm -1 ("(ii)") are those corresponding to the experimental observation (97 cm -1 and 135 cm -1 ). The scaling factor for vibrational analysis using DFT (b3lyp functional) is 0.961 (only ~4% difference) [41] and can be neglected because the shift of the modes due to the chain elongation is larger than 4%.   modes, i.e., {83, 86, 89 cm -1 } and {128,131,135 cm -1 }, due to the lowered symmetry of the chain after optimization as a single molecule. The Raman active vibration at 190 cm -1 appearing as a weak peak (iii) in Figure 3b is also displayed in Figure 4c, because it is observed in the simulation as explained later. This mode is related to the rotation around ketone and becomes lowered (164-189 cm -1 ) when the molecule becomes longer ( Figure  3c).

MD Simulations
The MD simulation was performed using LAMMPS with GAFF and a 0.1 fs time step. The atomic trajectory was recorded at 5 fs intervals. It was performed under NPT constraint, in which the number of particles (N) was constant and the pressure (P) and temperature (T) were controlled to be certain values with feedback. The following numerical analysis including temporal Fourier transform can be performed in several minutes. It is much shorter than the vibrational analysis of a polymer system using other techniques. Figure 5 shows snapshots from the MD trajectory at 300 K and 1 atm. The video can be found in supplement, Video S1(time interval corresponding to 90 cm-1) and Video S2(135 cm -1 ). Gaussian peaks with 5 cm −1 broadening were summed up after multiplying the calculated Raman intensities.
By scrutinizing the atomic motion of each mode (Figure 4), these modes were roughly identified as follows: The 85~98 cm −1 (mode (i)) is scissoring of ether (Ph-O-Ph) moiety and the 130~145 cm −1 (mode (ii)) is twisting of ketone and connected benzene rings (Ph-CO-Ph). In Figure 4, the molecule corresponding to Figure 3b is shown. The shorter molecule in Figure 3a has the vibration mode corresponding to Figure 3a,b at 90 cm −1 and 138 cm −1 , respectively. For the longer molecule in Figure 3c, these modes are split into several modes, i.e., {83, 86, 89 cm −1 } and {128,131,135 cm −1 }, due to the lowered symmetry of the chain after optimization as a single molecule. The Raman active vibration at 190 cm −1 appearing as a weak peak (iii) in Figure 3b is also displayed in Figure 4c, because it is observed in the simulation as explained later. This mode is related to the rotation around ketone and becomes lowered (164-189 cm −1 ) when the molecule becomes longer (Figure 3c).

MD Simulations
The MD simulation was performed using LAMMPS with GAFF and a 0.1 fs time step. The atomic trajectory was recorded at 5 fs intervals. It was performed under NPT constraint, in which the number of particles (N) was constant and the pressure (P) and temperature (T) were controlled to be certain values with feedback. The following numerical analysis including temporal Fourier transform can be performed in several minutes. It is much shorter than the vibrational analysis of a polymer system using other techniques. Figure 5 shows snapshots from the MD trajectory at 300 K and 1 atm. The video can be found in supplement, Video S1 (time interval corresponding to 90 cm −1 ) and Video S2 (135 cm −1 ).

Local Vibrational Spectra Simulated from MD
The trajectory of atoms from the MD simulation was used to calculate the local valence electron densities (ρ) at a 5 fs interval. The 2 × 3 × 3 × 12 (=216) points were selected in the unit cell and ρ of each point was calculated as a function of time using Equation (4). The points are identified using two characters as shown in Figure 6: circled numbers (1)~(18) in the ab plane and letters A~L and M~X in the c direction.

Local Vibrational Spectra Simulated from MD
The trajectory of atoms from the MD simulation was used to calculate the local valence electron densities (ρ) at a 5 fs interval. The 2 × 3 × 3 × 12 (=216) points were selected in the unit cell and ρ of each point was calculated as a function of time using Equation (4). The points are identified using two characters as shown in Figure 6: circled numbers (1)~(18) in the ab plane and letters A~L and M~X in the c direction.

Local Vibrational Spectra Simulated from MD
The trajectory of atoms from the MD simulation was used to calculate the local valence electron densities (ρ) at a 5 fs interval. The 2 × 3 × 3 × 12 (=216) points were selected in the unit cell and ρ of each point was calculated as a function of time using Equation (4). The points are identified using two characters as shown in Figure 6: circled numbers (1)~(18) in the ab plane and letters A~L and M~X in the c direction.  Circled (1)~(9) in the ab plane were coupled with A~L in the c direction to make 9 × 12 points, and circled (10)~(18) were coupled with M~X to make another 9 × 12 points. The local vibrational features were calculated as follows. Approximated "valence electron density" was obtained by Equation (4) using an atomic coordinate snapshot of 300 K 1 atm NPT simulation recorded at 5 fs intervals. The sum of Equation (4) at each position (3 × 3 × 12 (=108) positions ( 1 A~9 L)) was taken over the atoms in the same unit cell and nearest-neighboring unit cells (total 3 × 3 × 3 = 27) using a cut-off distance of 10 Å. A time duration of 45 ps after reaching thermal equilibrium was used for the calculation. The approximated "valence electron density" had 3 × 3 × 12 array data with 9000 points. They were Fourier-transformed and the absolute values were smoothed by smoothingcubic spline. The results are shown in Figure 7 after converting the frequency axis to the wavenumber unit using Equation (5). The scaling factor for the vibration frequency derived from the MD using GAFF is reported to be 1.0 [41] and there is no need for scaling the frequency. Circled (1)~(9) in the ab plane were coupled with A~L in the c direction to make 9 × 12 points, and circled (10)~(18) were coupled with M~X to make another 9 × 12 points. The local vibrational features were calculated as follows. Approximated "valence electron density" was obtained by Equation (4) using an atomic coordinate snapshot of 300 K 1 atm NPT simulation recorded at 5 fs intervals. The sum of Equation (4) at each position (3 × 3 × 12 (=108) positions (①A~⑨L)) was taken over the atoms in the same unit cell and nearest-neighboring unit cells (total 3 × 3 × 3 = 27) using a cut-off distance of 10 Å. A time duration of 45 ps after reaching thermal equilibrium was used for the calculation. The approximated "valence electron density" had 3 × 3 × 12 array data with 9000 points. They were Fourier-transformed and the absolute values were smoothed by smoothing-cubic spline. The results are shown in Figure 7 after converting the frequency axis to the wavenumber unit using Equation (5). The scaling factor for the vibration frequency derived from the MD using GAFF is reported to be 1.0 [41] and there is no need for scaling the frequency. We notice that several prominent peaks are observed in the low-frequency region in Figure 7. Typical ones are picked up and magnified in Figure 8. It is noted that 175 cm -1 is very local in some positions in column ① (Figure 8 top left), which coincides with the We notice that several prominent peaks are observed in the low-frequency region in Figure 7. Typical ones are picked up and magnified in Figure 8. It is noted that 175 cm −1 is very local in some positions in column 1 (Figure 8   The results presented here show that Fourier transform of the "valence electron density" (Equation (4)) reproduces the precise vibrational analysis by DFT and the experimental observation. Considering that the 135 cm -1 Raman signal is much weaker than that of 90 cm -1 in DFT of single molecules, and the MD combined with Equation (4) gives prominent peaks at 135 cm -1 at the intermolecular positions, we can presume that the 135 cm -1 mode is strongly affected by intermolecular interaction. This identification is consistent with the different behavior of 90 cm -1 and 135 cm -1 peaks during the improvement of crystallinity by the annealing experiment.
We must admit that there are some insufficient points in the present analysis of Raman spectra-the 175 cm -1 peak found in the analysis was not experimentally observed. It may be explained by the selection rule considering the phase factor of the electric field of light. However, the result can be obtained with much fewer computational resources compared with first-principle techniques that are difficult to apply to complicated crystalline polymers and composites.

WAXRD Pattern from MD Simulation
X-ray diffraction calculation was implemented in LAMMPS. We used this feature to explain the experimental results of WAXRD. First, we used large supercells of infinite polymers, but the structure did not change even when the simulation temperature was raised to 1000 K. We considered the periodic boundary condition (PBC) with connected The results presented here show that Fourier transform of the "valence electron density" (Equation (4)) reproduces the precise vibrational analysis by DFT and the experimental observation. Considering that the 135 cm −1 Raman signal is much weaker than that of 90 cm −1 in DFT of single molecules, and the MD combined with Equation (4) gives prominent peaks at 135 cm −1 at the intermolecular positions, we can presume that the 135 cm −1 mode is strongly affected by intermolecular interaction. This identification is consistent with the different behavior of 90 cm −1 and 135 cm −1 peaks during the improvement of crystallinity by the annealing experiment.
We must admit that there are some insufficient points in the present analysis of Raman spectra-the 175 cm −1 peak found in the analysis was not experimentally observed. It may be explained by the selection rule considering the phase factor of the electric field of light. However, the result can be obtained with much fewer computational resources compared with first-principle techniques that are difficult to apply to complicated crystalline polymers and composites.

WAXRD Pattern from MD Simulation
X-ray diffraction calculation was implemented in LAMMPS. We used this feature to explain the experimental results of WAXRD. First, we used large supercells of infinite polymers, but the structure did not change even when the simulation temperature was raised to 1000 K. We considered the periodic boundary condition (PBC) with connected covalent bonding to be too strong a restriction. Next, we tried to remove the PBC and run an NPT simulation. The structure soon became amorphous when the temperature was raised beyond 400 K, and crystallization was not observed even though the system was cooled on the scale of tens of ns. This is because the in-silico crystallization requires time beyond ordinary computation in the case of polymer simulation [42]. We therefore prepared crystalline phenyl-terminated oligomers corresponding to 9 × 9 × 5 unit cells and set PBC to the whole supercell to run the NPT calculation. Then the WAXRD patterns from the simulation changed according to the temperature. The results are shown in Figure 9. Still the PBC has an influence and the crystalline order is kept beyond the melting point of PEEK. However, the peak intensities diminish when the temperature is raised. We believe that this result can be used to study the effect of the degree of disorder on the WAXRD. covalent bonding to be too strong a restriction. Next, we tried to remove the PBC an NPT simulation. The structure soon became amorphous when the temperat raised beyond 400 K, and crystallization was not observed even though the syst cooled on the scale of tens of ns. This is because the in-silico crystallization requi beyond ordinary computation in the case of polymer simulation [42]. We theref pared crystalline phenyl-terminated oligomers corresponding to 9 × 9 × 5 unit c set PBC to the whole supercell to run the NPT calculation. Then the WAXRD from the simulation changed according to the temperature. The results are show ure 9. Still the PBC has an influence and the crystalline order is kept beyond the point of PEEK. However, the peak intensities diminish when the temperature i We believe that this result can be used to study the effect of the degree of disorde WAXRD. Figure 9. WAXRD patterns simulated using LAMMPS for a PEEK oligomer supercell corres to 9 × 9 × 5 unit cells.
In Figure 2, the initial state of the sample was mainly amorphous obtained by ing and the crystallinity became improved as the annealing temperature increase certain time step. Since the simulation temperature corresponds to thermal disord supercell, it is reasonable to compare the peak area as a function of temperatur reverse direction as shown in Figure 10. It is noted that peak ratio C/A increases w crystallinity is increased. This agrees well with the experimental result shown in F Peaks A and C correspond to 110 and 200 diffractions, respectively. The result ca plained by considering the order of the nearest neighbor (110) to be first establis that of the next nearest neighbor (200) to follow to improve when the crystallini hanced during annealing. Peaks B and D correspond to 111 and 211 diffractions lieve the effect of PBC works strongly in the c-direction (the third coordinate) bec supercell did not exist in this direction and thus discrepancy with the experime served for these peaks. . WAXRD patterns simulated using LAMMPS for a PEEK oligomer supercell corresponding to 9 × 9 × 5 unit cells.
In Figure 2, the initial state of the sample was mainly amorphous obtained by quenching and the crystallinity became improved as the annealing temperature increased with a certain time step. Since the simulation temperature corresponds to thermal disorder in the supercell, it is reasonable to compare the peak area as a function of temperature in the reverse direction as shown in Figure 10. It is noted that peak ratio C/A increases when the crystallinity is increased. This agrees well with the experimental result shown in Figure 2. Peaks A and C correspond to 110 and 200 diffractions, respectively. The result can be explained by considering the order of the nearest neighbor (110) to be first established and that of the next nearest neighbor (200) to follow to improve when the crystallinity is enhanced during annealing. Peaks B and D correspond to 111 and 211 diffractions. We believe the effect of PBC works strongly in the c-direction (the third coordinate) because the supercell did not exist in this direction and thus discrepancy with the experiment is observed for these peaks.

Future Improvement of the Present Method
We have tried to extract information on intermolecular vibration using a con used in machine learning of microchemical environments. It is applied to the classical trajectory of PEEK crystal and it shows agreement with the experiment. However, we mit that the estimation of the Raman intensity of each mode requires a more accurat timation. It is desired to combine the MD with quantum chemical techniques, and s research in this direction has been reported on quantum MD [43,44]. The reported m ods still require heavy computational resources and cannot be applied to large syste For example, the use of machine learning techniques to identify the position and direc of Raman-sensitive chemical moieties will be the direct extension of the present work

Conclusions
We have proposed and demonstrated a new analysis procedure for MD time tra tory to derive approximated Raman spectra in the low-frequency region. It involves local fluctuation of "valence electron density" as a function of time, which is converte a vibrational spectrum by the Fourier transform. From the calculation of PEEK polym 90 cm -1 , 135 cm -1 , and 175 cm -1 peaks were observed, among which 90 cm -1 and 135 peaks were experimentally observed. The 135 cm -1 signal was strongly observed at intermolecular point, and it is suggested that it is strongly influenced by intermolec interaction. This explanation agreed well with the experimental behavior of the peak tensity during crystallization. The WAXRD experiment was compared with the MD w different crystallinity at different temperatures, and nearest-neighbor ordering is m quickly established than that involving the next-nearest neighbors. Although ther room for improvement, the present technique provides a convenient way to estimate vibrational spectra of polymers and their composites with much fewer computationa sources than ordinary techniques.

Future Improvement of the Present Method
We have tried to extract information on intermolecular vibration using a concept used in machine learning of microchemical environments. It is applied to the classical MD trajectory of PEEK crystal and it shows agreement with the experiment. However, we admit that the estimation of the Raman intensity of each mode requires a more accurate estimation. It is desired to combine the MD with quantum chemical techniques, and some research in this direction has been reported on quantum MD [43,44]. The reported methods still require heavy computational resources and cannot be applied to large systems. For example, the use of machine learning techniques to identify the position and direction of Raman-sensitive chemical moieties will be the direct extension of the present work.

Conclusions
We have proposed and demonstrated a new analysis procedure for MD time trajectory to derive approximated Raman spectra in the low-frequency region. It involves the local fluctuation of "valence electron density" as a function of time, which is converted to a vibrational spectrum by the Fourier transform. From the calculation of PEEK polymer, 90 cm −1 , 135 cm −1 , and 175 cm −1 peaks were observed, among which 90 cm −1 and 135 cm −1 peaks were experimentally observed. The 135 cm −1 signal was strongly observed at the intermolecular point, and it is suggested that it is strongly influenced by intermolecular interaction. This explanation agreed well with the experimental behavior of the peak intensity during crystallization. The WAXRD experiment was compared with the MD with different crystallinity at different temperatures, and nearest-neighbor ordering is more quickly established than that involving the next-nearest neighbors. Although there is room for improvement, the present technique provides a convenient way to estimate the vibrational spectra of polymers and their composites with much fewer computational resources than ordinary techniques.
Author Contributions: Conceptualization, M.Y. and T.S.; software, X.Y. and T.S.; validation, S.Y. and T.N. All authors have read and agreed to the published version of the manuscript.