Der Pharma Chemica
Journal for Medicinal Chemistry, Pharmaceutical Chemistry, Pharmaceutical Sciences and Computational Chemistry

All submissions of the EM system will be redirected to Online Manuscript Submission System. Authors are requested to submit articles directly to Online Manuscript Submission System of respective journal.

Research Article - Der Pharma Chemica ( 2018) Volume 10, Issue 6

Molecular Structure, Nbo Charges, Vibrational Assignments, Homo-Lumo and Fukui Functions of Pyrazino [2, 3-D] Pyrimidine based on DFT Calculations

Ghamri M1, Harkati D1*, Saleh BA2*, Chikhaoui AR1, Belaidi S1

1Department of Chemistry, Faculty of Sciences, University of Biskra, Group of Computational and Pharmaceutical Chemistry, LMCE Laboratory, 07000 Biskra, Algeria

2Department of Chemistry, College of Science, University of Basrah, Basrah, Iraq

*Corresponding Author:
Harkati D
Department of Chemistry
Faculty of Sciences
University of Biskra
Group of Computational and Pharmaceutical Chemistry
LMCE Laboratory, 07000 Biskra, Algeria
Saleh BA
Department of Chemistry
College of Science
University of Basrah
Basrah, Iraq

Abstract

The optimized molecular structures, vibrational frequencies, corresponding vibrational assignments and atomic NBO charges of a biologically important molecule namely, pyrazino[2,3-d] pyrimidine have been carried out through density functional theory (B3LYP) level employing several basis sets. The results of vibrational frequencies and optimized geometric parameters were found in good agreement with the experimental data, which indicate that is planar and more stable. Moreover, the possible electrophilic and nucleophilic reactive sites were predicted by local reactivity descriptors and MEP map to describe the presence of intermolecular bonding and the hydrogen bonding interaction. In addition, the design format of our molecule present an opportunity to reveal the electron density profile by using the FMOs analysis and other electronic structure parameters such as ionization energy (I), Electron affinity (A), Global hardness (η), Chemical potential (μ), Electrophilicity index(ω), Electronegativity (X), Softness (S).

Keywords

Vibrational spectroscopy, Fukui function, NBO charge, DFT, FMO.

Introduction

One of the main objectives of organic and medicinal chemistry is to design, synthesize and produce molecules possessing value as human therapeutic agents. Compounds containing heterocyclic ring systems are of great importance receiving special attention as they belong to a class of compounds with proven utility in medicinal chemistry [1]. Pteridines, structurally the term “pteridine” describes the pyrazino [2,3-d] pyrimidine nucleus, which belong to a class of nitrogen heterocyclic compound present in a wide range of living systems. In the last century, Frederick Gowland Hopkins isolated a yellow pigment from the English Brimstone butterfly and a white pigment from the cabbage white butterfly [2].

These pigments in crystalline form and named them according to their colors, xanthopterin and leucopterin respectively [3]. There after isoxanthopterin and erythropterin were isolated from tropical butterflies and the isolation of pteridin europterin (later to be recognized as xanthopterin) from human urine [4]. In insects, they occur as metabolic end products and function as cofactors in hydroxylation reactions and as pigments. They are localized in the cuticle, wing scales, hypodermis, compound eyes, nervous system, light organ (of Lampyridae) and numerous other structures. The kind and quantity of pteridines found in insect tissues vary with developmental stages [5]. Although less symmetric, pteridine is of great intrinsic interest since it is the parent compound of a very large family of biologically important compounds [6]. So some pteridine derivatives are practically used in the chemotherapy or diagnosis of various diseases. Pteridines play an essential role in growth processes and the metabolism of one carbon unit as cofactors in enzyme catalysis and in biological coloration [7]. Many synthetic pteridines proved useful in medicine as anticancer, antiviral, antibacterial and diuretic drugs [8].

Conceptual Density Functional Theory [9-12] (DFT) is a very elegant way to describe chemical processes and provides a powerful theoretical framework for the study of both reactivity and selectivity [12]. In addition, Density functional theory methods offer an alternative [13] use of inexpensive computational methods which could handle relatively large molecules [14-16] and has been favorite one due to its great accuracy in reproducing the experimental values of molecular geometry, vibrational frequencies, atomic charges, etc. [17-29]. Thus, the present study is destined to study molecular geometry, electronic properties, NBO atomic charges, reactive sites by molecular electrostatic MEP/local reactivity descriptors, frontier molecular orbitals (FMOs) and the theoretical calculation of vibrational properties which is used to understand the spectra of large number of donor-acceptor systems [30-32]. All of these calculations have been carried out by Density Functional Theory (DFT).

Materials and Methods

Computational details

For the quantum chemical calculations, the Hyper Chem program package version 8.08, [33] was used to get the initial geometry optimization of pteridine which leading to energy minima and are achieved by using Molecular Mechanics Force Field (MM+) followed by semi-empirical PM3 method with a gradient norm limit of 0.1 kcal/Å for the geometry optimization. Then, a parallel study has been performed using Gaussian 09 program package, [34] gradient-corrected DFT with Becke 3 exchange [35] and Lee-Yang-Parr correlation functions (B3LYP) [36] with three different basis sets 6-311G, 6-311G (d,p) and 6-311G++(d,p), the overall view of pteridine with atom numbering is stimulated by MarvinSketch 6.2.1 software [37].

In the next step, the vibrational frequencies are calculated and scaled by the appropriate factor, the corresponding IR intensities, Raman activity of title molecule and the vibrational wavenumber assignments have been carried out by the GaussView 05, [38] and the Avogadro 1.0.1 software [39]. The electronic properties such as HOMO and LUMO energies were determined by the DFT approach. For obtaining chemical reactivity of the molecule, the molecular electrostatic potential (MEP) surface and the frontier molecular orbitals (FMOs) were visualized by the Gauss View 05 software, Fukui functions were calculated using Natural Bond Orbital (NBO) analysis by doing the single point energy calculation of N, N − 1 and N + 1 species of the molecule with same basis set 6-311G for analysis of molecular orbitals, the Gauss-Sum 3.0 program [40] was used to simulate the density of states (DOS) spectrum.

Results and Discussion

Molecular geometry analysis

The calculations of geometrical parameters of title compound were carried out calculated by the DFT/B3LYP method, in conjugation with three different basis sets 6-311G, 6-311G (d,p) and 6-311G++ (d,p). The results are listed in Table 1. Figure 1 depicts the overall view of pteridine with atom numbering.

Parameters EXP DFT/B3LYP
6-311G 6-311G(d,p) 6-311G++(d,p)
Bond length (A°)
N(1)–C(2) 1.304 1.3271 1.3133 1.3133
N(1)–C(9) 1.361 1.3712 1.3559 1.3558
C(2)–N(3) 1.364 1.3793 1.3633 1.3631
C(2)–H(11) 1.01 1.0776 1.0858 1.0858
N(3)–C(4) 1.35 1.3283 1.3126 1.3127
C(4)–C(10) 1.399 1.4171 1.4172 1.4176
C(4)–H(12) 1.1 1.0812 1.0875 1.0874
N(5)–C(6) 1.308 1.3281 1.3121 1.3118
N(5)–C(10) 1.384 1.3753 1.3594 1.3589
C(6)–C(7) 1.398 1.4231 1.4241 1.4245
C(6)–H(13) 1.09 1.0806 1.0866 1.0865
C(7)–N(8) 1.286 1.3289 1.3122 1.3119
C(7)–H(14) 0.87 1.0812 1.0876 1.0874
N(8)–C(9) 1.359 1.3744 1.3614 1.361
C(9)–C(10) 1.392 1.419 1.4189 1.4195
Bond angle(°)
C(2)–N(1)–C(9) 115.2 116.5942 116.1787 116.3381
N(1)–C(2)–N(3) 129 126.6008 127.9528 127.7701
N(1)–C(2)–H(11) 108 117.5376 116.7988 116.898
N(3)–C(2)–H(11) 115.8617 115.2484 115.3318
C(2)–N(3)–C(4) 115.7 117.0518 116.2098 116.3617
N(3)–C(4)–C(10) 119.8 121.4228 122.0316 121.9644
N(3)–C(4)–H(12) 116 118.7996 118.8683 118.7382
C(10)–C(4)–H(12) 119.7776 119.1002 119.2974
C(6)–N(5)–C(10) 114.5 116.1983 115.3681 115.5034
N(5)–C(6)–C(7) 123 121.9983 122.3063 122.2295
N(5)–C(6)–H(13) 124 117.2294 117.6232 117.7226
C(7)–C(6)–H(13) 120.7723 120.0705 120.048
C(6)–C(7)–N(8) 123.4 122.8506 123.3538 123.2931
C(6)–C(7)–H(14) 120.4038 119.5447 119.5018
N(8)–C(7)–H(14) 113 116.7456 117.1015 117.2051
C(7)–N(8)–C(9) 116.4 116.3418 115.6552 115.7876
N(1)–C(9)–N(8) 118.3447 118.225 118.3628
N(1)–C(9)–C(10) 121.2 120.916 121.0035 120.936
N(8)–C(9)–C(10) 120.9 120.7393 120.7715 120.7012
C(4)–C(10)–N(5) 120.7139 120.8313 120.8851
C(4)–C(10)–C(9) 119.1 117.4144 116.6236 116.6297
N(5)–C(10)–C(9) 121.8 121.8716 122.5451 122.4852

Table 1: The experimental and theoretical bond lengths (A°), bond angles (°) and dihedral angles (°) of pteridine

derpharmachemica-system-pteridine

Figure 1: The structural and numbering system of pteridine

A close look at Table 1 indicates that our theoretical results are agreed with the available experimental data [41]. Plotting of the theoretical results versus the experimental ones indicates robust linear dependencies (Figure 2). The values of the correlation coefficients (r) were 0.325(6- 311G), 0.236 (6-311G (d, p)) and 0.654(6-311G++ (d, p)) for bond lengths; and the bond angles were 0.325(6-311G), 0.236 (6-311G (d, p)) and 0.654(6-311G++ (d, p)). Hence, these figures suggest that the B3LYP/6–311G bond lengths and bond angles are well fitted to X-ray experimental data, with the correlation coefficients R=0.99635 and R=0.56786 respectively. Although there are no experimental bond dihedral angles, the calculated dihedrals are present within the same table. The dihedrals angles of this molecule vary between 0° and 180° degree which explain that the geometry of pteridine is planar, which make this conformation more stable.

derpharmachemica-bond-lengths

Figure 2: Correlation graphics of : a) calculated bond lengths vs. Experimental bond lengths (A˚), b) calculated bond angle vs. experimental bond angle (degree) at B3LYP/6-311G level

Atomic charges

Natural Bond Orbital (NBO) analysis employs a more sophisticated scheme that uses the concept of natural orbitals and occupation-based symmetric orthogonalization that leads to more robust results. The net result of NBO analysis is that the total density is broken into localized contributions associated with individual atoms giving rise to natural charges and further into 1-center (lone pairs) and two-center (bonds) pieces. This leads to a straight toward determination of Lewis structures. NBO also allows one to determine hybridization of the AOs contributing to a particular bond [42-43]. Atomic net charges were calculated applying the natural bond orbital (NBO) analysis [44-45] at the B3LYP/6–311G level of calculation.

Figure 3 shows the atomic NBO charges. The magnitudes of the carbon atomic charges are found positive and C (9) has maximum positive atomic charge values of about 0.325 electron. Whereas the magnitudes of charges on nitrogen atoms are found negative and the maximum negative atomic charge was obtained for N (3), -0.449 electron. In addition, the magnitudes of the hydrogen atomic charges are arranged in an order from 0.196 to 0.212 electrons.

derpharmachemica-Structure-pteridine

Figure 3: Structure of pteridine and the calculated NBO charges at B3LYP/6-311G level

Vibration analysis

The molecule pteridine consists of 14 atoms and it possesses thirty-six possible modes of vibration, hence the molecule being of low symmetry, most of these vibrations should be active in the infrared [46]. FTIR spectrum of the investigated compound was measured in the 4000-400 cm-1 region and the vibrational frequencies of the compound were calculated in region 4000-0 cm-1 by using B3LYP/6-311Gmethod (Figure 4).

derpharmachemica-spectrum

Figure 4: (a) FT-IR spectrum, (b) Simulated B3LYP level IR spectra of the title

Then Table 2 presents the detailed vibrational assignment scaled by the factors 1.06 with observed wavenumbers, Infrared Radiation (IR) intensities and Raman activity of title molecule, respectively. All compared data have been shown to be in a good agreement to each other for the most stable optimized structure of the pteridine molecule.

Assignment Experimental frequencies [46] Calculated frequencies (cm−1)
B3LYP 6–311G
Theoretical frequency IR intensity (km/mol) R activity (Å4/amu)
γ NCH - 159.76 0.011 0.154
γ (Pym-Pyz) - 195.07 1.803 0.452
γ CH (Pym) - 412.52 3.462 4.02
β CH (Pym) 434 426.42 4.096 0.265
γ CH (Pym-pyz) + γ CCH ( Pym) 455 478.63 15.621 0.336
γ NCN (Pym-Pyz) + γ CH ( pym) 497 529.34 2.914 0.238
β CNC (Pym) + β CCN (Pyz) 533 540.7 0.769 12.18
β NCC (Pyz) + β CNC (Pyz) 616 567.08 0.324 9.089
β CNC (Pym-Pyz) + β CCN (Pym) 653 631.48 7.327 2.082
γ CH (Pym-Pyz) 777 690.19 13.785 0.407
β CH (Pym-Pyz) + υ C=C (Pym-Pyz) 795 772.9 5.476 31.956
γ CH (Pym-Pyz) 822 853.96 1.565 1.003
β CNC (Pym) + β CCN (Pyz) 860 868.95 6.843 0.282
γ CH (Pyz) 882 891.45 15.743 0.256
β NCC (Pyz) + β NCH (pym) 933 951.55 24.477 0.397
γ CH (Pym) 1015 956.62 24.095 0.23
γ CH (Pym) 1056 1003.39 0.377 0.53
β CH (Pyz) 1085 1008.29 0.003 2.461
β (Pym-Pyz) + υ C=C (Pyz) 1092 1048.04 17.735 3.926
β (Pym-Pyz) + υ C-N (Pym) 1109 1090.91 7.812 12.828
β (Pym-Pyz) + υ C-N (Pym-Pyz) 1158 1181.19 0.022 0.606
β NCN (Pym-Pyz) + β NCC (Pym-Pyz) 1197 1198.23 13.273 0.833
υ C=C (Pym-Pyz) + β (Pym-Pyz) 1286 1300.82 2.906 11.858
β CH (Pym-Pyz) 1348 1310.92 4.28 0.73
β CH (Pym) + υ C=C (Pym-Pyz) 1368 1362.42 7.83 68.004
β CH (Pym-Pyz) 1391 1392.24 16.222 47.435
β CH (Pym-Pyz) 1400 1423.32 11.406 25.691
β (Pym-Pyz) + υ C=C (Pyz) 1435 1457.74 61.776 0.247
β CH (Pym) 1461 1476.84 4.262 4.519
β (Pym-Pyz) + υ C=C (Pym-Pyz) 1558 1542.08 27.248 39.991
υ C=N (Pyz) + β (Pyz) 1573 1575.12 9.463 5.151
υ C=N (Pym) + β (Pym) 1580 1592.98 57.405 2.825
υas CH (Pyz) 3035 3178.71 5.539 96.235
υas CH (Pym) 3055 3190.47 7.84 113.199
υs CH (Pyz) 3090 3198.32 34.478 265.26
υs CH (Pym) 3117 3227.93 16.817 173.514

Table 2: Comparison of the experimental and calculated vibrational frequencies (cm-1)

C–H vibrations

The C–H stretching vibrations of heteroaromatic structure are expected in the region 31000-3000 cm-1 [47-48]. In this study, the C–H stretching vibrations were observed at 3178.71, 3190.47, 3198.32, 3227.93 cm−1 theoretically, the corresponding modes were showed at 3035, 3055, 3090, 3117 cm−1 and 96.235, 113.1994, 265.2597, 173.514 cm−1 in FTIR and Raman, respectively.

The CH in-plane bending vibrations interact somewhat with C–C stretching vibrations occurs in the region 1000-1300 cm-1 [49]. The C-H in-plane bending vibrations were predicted at 1310.92, 1362.42, 1392.24, 1423.32, 1476.84 cm-1 and the corresponding modes were observed as strong peaks at 1348, 1368, 1391, 1400, 1461 cm-1 in Fourier Transform Infrared (FTIR) spectroscopy and in Raman at 0.7304, 68.0035, 47.435, 25.6914, 4.5186 cm-1.

The CH out-of-plane bending vibrations appear in the region 1050-700 cm-1 [46]. Accordingly, The FTIR bands observed as weak to strong peaks at 690.19, 853.96, 891.45 and 956.62 cm-1 are assigned to C-H out-of-plane bending vibrations, which are in good agreement with the experimental values 777, 822, 882 and 1015 cm-1 at 0.4071, 1.0029, 0.2556 and 0.2301 cm-1 in the Raman spectrum.

C=N, C-N vibrations

The identification of the C–N stretching frequency in the side chains is a rather difficult task since there are problems in identifying these frequencies from other vibrations [50]. In the present study, the ring C-N stretching vibrational modes occur in the region 1600-1080 cm-1 where the weak to strong intense bands observed at 1109, 1158, 1573, 1580 cm-1 in FTIR and the bands at 12.8276, 0.6055, 5.1509, 2.8246 cm-1 in Raman are allotted to C-N, C=N stretching vibration. The observed results are well matching with calculated values.

C=C Vibrations

The aromatic ring C=C stretching vibrations usually present in the region 1635-1080 cm-1 [48]. In the present study, the bands were of variable intensity and were found at 1092, 1286, 1368, 1435, 1558 cm-1 in FT-IR which are assigned to C-C stretching vibrations in pteridine. The same vibrations were observed in Raman at 3.9258, 11.8581, 68.0035, 0.2466 and 39.9912 cm-1. The calculated values at 1048.04, 1300.82, 1362.42, 1457.74 and 1542.08 cm-1, respectively, these assignments are in quite agreement with experimental data.

Ring Vibrations

Generally, ring vibrational modes are quite sensitive to substitution [51-54] where the bands were found at 1048.04, 1090.91, 1181.19, 1300.82, 1457.74, 1542.08, 1575.12 and 1592.98 cm-1 in FTIR. These bands are attributed to ring in plane bending modes, while experimentally found modes were found at 1092, 1109, 1158, 1286, 1435, 1558, 1573 and 1580 cm-1, whereas in Raman were found at 3.9258, 12.8276, 0.6055, 11.8581, 0.2466, 39.9912, 5.1509 and 2.8246 cm-1.

The ring out-of-plane bending mode frequencies were at 195.07 cm-1 and in Raman at 0.4515 cm-1, while no vibration was identified in experimental data.

The bands at 933, 455 cm-1 in the infrared and 0.3971, 0.3362, 0.1543 cm-1 in Raman were assigned to the N-C-H in plane bending vibrations and C-C-H/N-C-H out plane bending vibrations, respectively, whereas the same bands were calculated at 951.55, 478.63, 159.76 cm-1. In addition to these bending vibrations, C-N-C, C-C-N and N-C-N bending vibrations were assigned in the characteristic range in pteridine molecule.

Frontier Molecular Orbitals (FMOs) analysis

The highest occupied molecular orbital (HOMO) represents the outermost orbital filled by electrons and behaves as an electron donor, while the lowest unoccupied molecular orbital (LUMO) considers as the first empty innermost orbital unfilled by electron and behaves as an electron acceptor [55]. These orbitals are also called as frontier molecule orbital (FMOs) which analysis predicts the molecular interactions with other species and has great importance in modern biochemistry and molecular biology [56]. According to B3LYP/6–311G calculation, The energy gap between HOMO and LUMO is -3.85839 eV which indicates the molecular chemical stability and it is a critical parameter to determine molecular electrical transport properties [55]. The calculate ionization energy and electron affinity can express through HOMO and LUMO orbitals energies as I= -EHOMO and A = -ELUMO. The term ‘η’ has come to be used to refer to the global hardness, η =½(ELUMO - EHOMO), the hardness has been associated with the stability of chemical system. The electron affinity can be used in combination with ionization energy to give electronic chemical potential, μ= ½(ELUMO + EHOMO). The term ‘ω’ refer to the global electrophilicity index, image [57]. Electronegativity,image, is the tendency of an atom to attract electrons. The term ‘S’ is generally understood to mean the chemical softness, image [58]. These values of Electronic structure parameters are given in Table 3.

Molecularproperties Energy
EHOMO(eV) -7.02695
ELUMO(eV) -3.16857
Energie gap (EHoMO -ELuMO)(eV) -3.85839
Ionizationenergy (eV) 7.02695
Electron affinity (eV) 3.16857
Global hardness (eV) 1.92933
Chemical potential (eV) -5.0979
Electrophilicity index (eV) 6.73523
Electronegativity (eV) 5.0979
Softness (eV-1) 0.51831

Table 3: Electronic structure parameters of the title compound

The Gauss Sum program was used to obtain the density of states DOS diagram which was convoluted from the molecular orbital data as shown in Figure 5.

derpharmachemica-molecular-orbitals

Figure 5: The (a) frontier molecular orbitals and (b) density of states spectrum of pteridine molecule

Local reactivity descriptors and molecular electrostatic potential

The local reactivity descriptors such as Fukui functions (fk+, fk, fk0) were calculated at the B3LYP/6–311G for the pteridine molecule using Natural Bond Orbital (NBO) analysis to predict the possible reactive sites of the molecule where a chemical species will change its density when the number of electrons is modified Thus, the propensity of the electronic density to deform at a given position upon accepting or donating electrons. Itis an important in designing a pharmaceutical compound [59-62]. The corresponding condensed or atomic Fukui functions (FF) can be evaluated through the following equations [63]:

f k+ = qk(N +1) − qk(N) for nucleophilic attack

f k = qk(N) − qk(N −1) for electrophilic attack

f ko= [qk(N + 1) − qk(N −1)] /2 for radical attack

and the importance of MEP, it displays molecular size, shape as well as positive, negative and neutral electrostatic potential regions in terms of color grading red < orange < yellow < green < blue. Where red, blue and green represent the regions of most negative, most positive and zero electrostatic potential respectively [64-66]. As depicted in (Figure 6), regions of negative potential are over the nitrogen atoms and the regions having the positive potential are over the hydrogen atoms.

derpharmachemica-electrostatic-potential

Figure 6: 3D plot of molecular electrostatic potential (MEP) for pteridine

On the other hand, It was noted in Table 4 that the highest value of fukui functions for the LUMO (fk+=0.1900) present around N (1), which indicates electrophilic site and hence possible site for nucleophilic attacks. The highest of fukui functions for the HOMO (fk-=0.1410) present around N (5) and N (8), which indicate nucleophilic site and these are the possible sites for electrophilic attacks. While, the atom N (1) characterized by the highest value of fk0 (fk0=0.1500) indicates that this atom is the most important site for radical attacks.

B3LYP/6-311G
Atom P (N) P (N+1) P (N-1) fk+ fk- fk0
1N -0.424 -0.234 -0.534 0.19 0.11 0.15
2C 0.26 0.3 0.208 0.04 0.052 0.046
3N -0.449 -0.304 -0.538 0.145 0.089 0.117
4C 0.113 0.143 -0.017 0.03 0.13 0.08
5N -0.385 -0.256 -0.526 0.129 0.141 0.135
6C 0.026 0.078 -0.055 0.052 0.081 0.0665
7C 0.038 0.089 -0.069 0.051 0.107 0.079
8N -0.374 -0.273 -0.515 0.101 0.141 0.121
9C 0.325 0.329 0.338 0.004 -0.013 -0.0045
10C 0.054 0.099 0.056 0.045 -0.002 0.0215
11H 0.196 0.247 0.153 0.051 0.043 0.047
12H 0.212 0.271 0.174 0.059 0.038 0.0485
13H 0.204 0.254 0.163 0.05 0.041 0.0455
14H 0.204 0.257 0.163 0.053 0.041 0.047

Table 4: Values of the condensed Fukui function in pteridine computed from NBO charges

Conclusion

The present study deals with the comparison between the calculated and experimental structural parameters indicate that B3LYP/6–311G is in good agreement with experimental observations. In addition, theoretically calculated and experimentally observed vibrational wavenumbers were compared and assigned on the same level of calculations. Our work was mainly based on the study of both reactivity and selectivity also oriented towards the search for the most stable molecular conformation and which are largely responsible for binding of small molecule to its active site. Thus, at the B3LYP/6-311G level, which provide the nature of reactivity, the energies of HOMO/LUMO and their orbital energy gaps were calculated as well as electronic structure parameters of the molecules. Fukui function with NBO charge was positive value around N (1) at places which are better at accepting electron than they are at donating electrons. Conversely, zones with negative values around N (5) and N (8) which means that these sites donate electrons rather than accept electrons.

References

FIND ARTICLES
SCImago Journal & Country Rank
Recommended Conferences
bornova escort karşıyaka escort osmangazi escort buca escort anne porno fetiş porno anal porno bartın escort burdur escort eskişehir escort escort izmir bursa escort porno escort bayan