First-principles DFT investigations of the vibrational spectra of chloro-quine and hydroxychloroquine

Since the COVID-19 pandemic began, two drugs, chloroquine (CQ) and hydroxychloroquine (HCQ), have received renewed attention. Using the density functional theory method in the CASTEP and DMol3 packages, we calculated both molecules’ infrared spectra and the partial phonon density of states of the hydroxyl group to identify the origin of the differences between the two spectra. Some characteristic vibrational modes of the hydroxyl group in HCQ were analysed individually. We also compared their Fukui functions and found that the oxygen atom in HCQ possesses electrophilic properties. This finding may be related to the large difference in toxicity between these two drugs. The method herein presents a new pathway to investigate organic molecules from the view of physics.


Introduction
The spread of the novel coronavirus disease 2019 (COVID- 19), caused by the new coronavirus SARS-CoV-2, presents a huge challenge to countries and public health institutions around the world [1]. By 6th Sep. 2021, more than 220 million cases and 4.6 million deaths had been reported [2]. The development of vaccination are imminent and the invetigations of COVID-19 are continuous updated [3]. It has been reported that some compounds have been discovered through CADD. They can be well bound with COVID-19 protease receptors [4]. 5-COI derivative of U may help to provide strongest complex of interacting ligand target system and prevent virus growth [5]. Some evidence suggests that chloroquine (CQ) is a highly effective inhibitor of the proliferation of the virus in vitro [6], leading to a renewed surge of interest in this anti-malarial drug.
CQ was first synthesised in 1934 [7]. As the malarial parasite became more resistant to CQ, the dosage of the drug was gradually reduced [8]. In recent decades, CQ and its structural analogue hydroxychloroquine (HCQ) have often been used as drugs to treat autoimmune diseases. In 1950s, the HCQ derivative of quinacrine was used in the treatment of systemic lupus erythematosus (SLE) [9]. After the success of that application, the effectiveness of anti-malarial drugs in the treatment of SLE began to be recognised. Due to its good safety and reliability, HCQ has been increasingly widely used in the clinic [10]. HCQ has also been recommended for combination therapy with other anti-rheumatic drugs [11]. In recent years, especially after the SARS epidemic in 2003, researchers have repeatedly demonstrated that this drug has the effect of inhibiting viral division [12]. However, the drug has uncertain effects on the virus that causes COVID-19. Many clinical trials have shown that the drug's anti-viral effectiveness is not promising [13]. Based on these investigations above, we used density functional theory (DFT) method to simulate these two drug molecules [14]. Therefore, we studied the differences in structure and toxicity of these two drugs from a view of physics.
HCQ is distinguished from CQ by a hydroxyl group at the end of the carbon chain. It is reported that the toxicity of CQ is twice that of HCQ [15,16]. In this work, we calculated the infrares (IR) spectra of these two molecules. We theoretically analysed the difference between their vibrational spectra and their Fukui functions from the view of physics. We hope that these insights will benefit the understanding of their antiviral effect.

Computational methods
The CQ and HCQ models seen in figure 1 were constructed in the Materials Studio platform. First-principles DFT method is capable to determine the most stable structure of these two drugs and is also able to present its IR spectrum. DFT calculations were first conducted using the DMol 3 package to calculate their IR spectra and Fukui functions after geometry optimization. The DMol 3 code can simulate the atomic vibrations for a nonperiodic structure based on a linear combination of atomic orbitals (LCAO) basis set. We chose the RPBE [17,18] exchange-correlation functional based on the generalised gradient approximation according to our experience [19]. The energy threshold and self-consistent field tolerance were set as 1×10 −6 and 1×10 −8 Hartree. Then, for phonon calculations, the HCQ molecule was used to construct a primitive cubic cell with P1 symmetry. To avoid inter-molecular interactions, the side length was set as 30 Å. The molecular properties of HCQ were recalculated with the CASTEP [20] code. Finally, the OH-related phonons were extracted from the total phonon density of states (PDOS), allowing the contributions of OH vibrational modes to be discussed.

Results and discussion
To distinguish the effect of hydroxyl group of the two drugs, we compared the vibrational spectra of the two structures. Under the harmonic approximations, the calculated phonons show good agreement with experimental spectra. The peaks in an IR spectrum can be assigned by analyzing the vibrational dynamics of corresponding normal modes. Figure 2 depicts the simulated partial PDOS of -OH in HCQ, and the IR spectra of CQ and HCQ. To facilitate the comparison, we separated the spectra into four parts labelled (a)-(d). The phonons of -OH below 200 cm −1 are mainly from acoustic branches that cannot be detected by IR absorption and were thus excluded from the three displayed spectra. As there are no frequencies from 1650 to 2900 cm −1 , this band was also omitted from the spectra. With the help of vibrational analysis, we discuss the IR spectra in terms of the normal modes. Under the harmonic approximation, the total number of vibrational normal modes is 3N − 3, where N is the number of atoms in a primitive cell. Many of the modes may be so close to each other that the peaks are observed as combinations. For the simulated IR spectra, we selected one representative mode with the highest intensity around each peak. It is obvious that the phonon spectra of CQ and HCQ are almost identical except for the absence of the -OH phonons from CQ. Because the PDOS were calculated with CASTEP whilst the IR spectra were calculated by DMol 3 , there are some discrepancies, especially in part (a) and part (d).
The H-O bending modes in HCQ are all below 531 cm −1 as shown in part (a). The peaks at 310 cm −1 in CQ and 298 cm −1 in HCQ are the N-H vibrations, which are shown in figure 3. Frosch et al observed a peak of CQ was at 275 cm −1 [21], and a peak of HCQ was measured at 340 cm −1 reported by Ejuh et al [22]. For HCQ, we attribute the peak at 276 cm −1 in the partial PDOS of -OH to be the peak at 339 cm −1 in the IR spectrum due to systematic error.  The O-C stretching peak is at around 991 cm −1 and the H-O-C-C dihedral angle bending modes are spread from 1081 to 1359 cm −1 . As can be seen from figure 2(b), three relatively strong peaks appear at 991, 1018 and 1081 cm −1 in the -OH PDOS of HCQ. As there are no phonons in this area in the CQ spectrum, the difference here between CQ and HCQ is remarkable. In figure 2(c), the peak at 1359 cm −1 in the -OH PDOS can be seen clearly at 1358 cm −1 in the IR spectrum of HCQ. It can be found that both CQ and HCQ have three strong N-H bending vibrational modes from 1500 to 1600 cm −1 . For CQ, the normal modes are at 1529, 1560 and 1574 cm −1 , which correspond to 1556 cm −1 , 1576 cm −1 and 1597 cm −1 in the experiment [21], whilst the equivalent modes of HCQ are at 1523, 1550 and 1571 cm −1 respectively and the data observed by Ejuh et al are at 1560 cm −1 , 1600 cm −1 and 1640 cm −1 , respectively [22]. The mode at 990 cm −1 is shown in figure 3, and the modes at 1081, 1357 and 1550 cm −1 are shown in figure 4. The corresponding movies of their dynamic processes can be found in the supplementary files.
For HCQ, the O-H stretching vibration is at 3723 cm −1 corresponds to the experimental data at 3800 cm −1 [22]. The corresponding mode in the PDOS of -OH is at 3636 cm −1 . Thus, the peak in the PDOS from CASTEP is red-shifted by 87 cm −1 relative to the IR spectrum. This is a characteristic peak of HCQ, which is clearly absent from the IR spectrum of CQ as shown in part (d). The N-H stretching vibrations are at 3597 and 3572 cm −1 in Figure 2. Comparison of simulated IR spectra between CQ and HCQ. Due to the wide variation of the intensity in different regions, we adjusted the intensity scales appropriately.
HCQ and CQ, respectively. The C-H stretching vibrations in the carbon chain are widely distributed from 2980 to 3140 cm −1 , which are from 2900 cm −1 to 3200 cm −1 in the experiment [22][23][24]. For HCQ, mixing with some C-H stretching vibrations of the benzene ring, they are at 3084, 3115, 3145, 3176 and 3228 cm −1 , which were observed at 3020, 3080, 3120 and 3280 cm −1 [22]. For CQ, the corresponding modes are from 2895 to 3259 cm −1 , which were observed from 2876 to 2993 cm −1 [23]. From figure 2(d), it is difficult to distinguish between CQ and HCQ in this region of the spectrum except for the -OH stretching mode at 3723 cm −1 .
To predict the reactivity and selectivity of atomic sites in a molecule, we calculated the Fukui functions of CQ and HCQ to predict the reactive sites of the molecules [25]. The Fukui function was proposed by Parr and Yang in 1984 [26] and aims to predict the sensitivity of molecular atoms to radical, nucleophilic and electrophilic attacks [27]. The condensed Fukui function is designed to characterise the chemical reactivity and selectivity of specific atomic sites in a molecule [28]. The results for CQ and HCQ calculated herein are presented in table 1. F(+), F(−) and F(0) represent nucleophilic attack, electrophilic attack and radical attack, respectively. In addition, the dual descriptor is represented as δF [29]. A positive number means that the atom possesses nucleophilic properties, whereas a negative number represents electrophilic properties. We can identify various sites for radical attack on the molecule [30]. Cl7, N18, N11 and C6 are the preferred radical attack sites of CQ, and Cl7, N11, C6 and C3 are the corresponding sites of HCQ. Meanwhile, the atoms Cl7, C3, C6, C10 and C2 in  CQ and Cl7, C3, C6, C10 and C8 in HCQ are predicted to be preferred sites for nucleophilic attack. We can clearly see that both CQ and HCQ are susceptible to nucleophilic attacks at similar sites. The preferred sites for electrophilic attack are predicted to be N18 for both CQ and HCQ. Meanwhile, some sites, like C9, N11 and N12, possess both nucleophilic and electrophilic properties. O23, which is present in HCQ but not CQ, possesses electrophilic properties. This is an obvious difference between the two molecules in nucleophilic and electrophilic properties. We speculate that this may be the reason for the toxicity difference between CQ and HCQ.

Conclusions
In this study, we calculated the normal modes of CQ and HCQ to examine the differences between their IR spectra. Notably, the contributions of the -OH group in HCQ seem to be almost independent of the rest of the molecule. If the phonons of the -OH group are subtracted from the total PDOS of HCQ, the two spectra seem identical. These theoretical results represent a new method to distinguish similar molecules.
To investigate the toxicity difference between these two drugs, we compared the nucleophilic and electrophilic properties of the two molecules and found that the oxygen atom in HCQ possesses electrophilic properties. Although we cannot demonstrate a direct relationship of the oxygen atom with the toxicity of these two drugs, the investigation method represents a new pathway to study organic molecules.
Based on the work above, we shed light on a method to explore drug molecules from the view of physics. We confirmed that the two vibrational spectra of CQ and HCQ are almost identical except for the -OH stretching mode at 3723 cm −1 . Besides, the oxygen atom possesses electrophilic propertie and maybe responseible for its higher toxicity.

Acknowledgments
The numerical calculations were performed on the supercomputing system in the Supercomputing Center, Shandong University, Weihai.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.