2.1. 안정된 삼차원 분자구조 계산

PC 내의 Chem-3D 프로그램을 사용하여 분자들을 스케치 한 뒤, Chem-3D에 내장된 반실험적 계산 엔진을 사용하여 기본적인 최적화 작업을 진행하고 본 계산에 필요한 입력자료를 작성하였다. 계산은 가우시안-03 프로그램을 사용하였으며,⁸ 안정된 구조들을 밀도함수이론(density functional theory, DFT)을 사용하여 B3LYP/6-31G^* 수준으로 계산하였다. 분자구조를 계산하는 과정에서는 어떠한 symmetry도 가하지 않았으며, 계산 최적화 과정도 cartesian coordinate 하에서 진행하였으며, Z-매트릭스는 사용하지 않았다.

2.2. 화학적 특성 예측 작업

가장 안정된 분자구조가 확인되면 이들 분자구조에서 화약성능 계산 작업에 필요한 밀도와 생성열을 예측하는 작업을 수행하였다. 생성열 예측은 B3LYP/6-31G*로 최적화된 분자구조에서 출발하여 DFT 계산의 BP86/6-31G** 수준으로 구조 최적화 및 에너지 계산을 재진행하였다. 본 DFT 계산결과를 이용하여 기존에 Politzer 교수가 개발한 방안에 따라 계산을 진행하였다.^9,10 Politzer 교수 계산법은 고에너지물질에 적용할 수 있도록 실험적인 보정치를 각 원소의 형태에 따라 기 진행하였으며, 이들 방안에 대해 국과연에서도 다양한 유형의 고에너지물질에 적용한 바 상당히 우수한 결과를 제공하고 있음을 확인할 수 있었다.¹¹

보편적으로 학계에서 화약 분자들의 밀도를 ‘그룹 영향력의 합’ 방법을 이용하여 계산하고 있지만 본 연구에서는 국과연과 인하대학교에서 공동으로 개발한 분자표면 정전기 전위를 이용하는 방법을 사용하였다.^12,13 일반적인 양자역학적인 분자 표면 계산 방법은 양자역학적으로 계산한 전자운(electron cloud)의 밀도가 특정치 이상 되는 지점을 연결하여 분자 표면으로 구성하고 있었다. 하지만 이들 김찬경 교수팀에서는 이들 분자표면을 사용하는 대신 분자역학적으로 널리 알려져 있는 구형의 반데르발스 직경을 사용하였으며, 이들을 봉합하는 분자 역학적 프로그램으로 삼차원적인 분자 표면을 구성한 뒤, 이들 위치에서 정전기 전위 값을 계산하였다.¹² 이들 분자 표면 정전기 전위가 계산되면 반데르발스 표면적과 분자 표면 정전기 전위 값의 평균 오차를 변수로 하는 새로운 QSPR 관계식을 사용하여 고체 화약 분자들의 밀도를 예측하는 작업을 수행하였는데 기존의 계산 결과에 비하여 우수한 결과를 제공하는 것을 확인하였다.¹³

2.3. 화약 성능 예측 작업

일단 각 화약 분자에 대한 고체 상태에서의 생성열과 밀도값이 정밀하게 예측되면, 화약 성능을 진행할 준비가 완료된 상태이다. 본 실험에서의 화약 성능 예측 작업은 미국 로렌스 리버모어 국립연구소의 프라이드 박사가 개발한 치타 프로그램 버전 2.0을 사용하여 진행하였으며,¹⁴ 그 결과를 기존에 널리 사용되는 화약들인 TNT, RDX, HMX와 비교하여 분석하였다.

3.1. 안정된 삼차원 분자구조 분석

밀도함수 이론인 B3LYP/6-31G* 방법으로 Fig. 1에 제시된 분자들을 최적화한 결과 Fig. 2와 같이 안정된 분자구조들을 획득할 수 있었다. 양자역학적으로 계산된 구조들은 큰 이상적 징후없이 안정된 형태로 계산이 되고 있다. 하지만 이미다졸은 방향성 구조로 6개의 추가적인 파이전자들이 고리 전체에 잘 비편재화되어 있어 평면 구조를 형성하고 있다. 본 계산에서도 모든 구조들이 평면으로 유지되고 있음을 확인할 수 있다. 또한 치환된 나이트로기들은 이미다졸 링과 파이전자들을 나누기 위해 서로 평면으로 유지되는 것이 가장 안정된 구조이며, 이런 형태는 주변의 다른 복잡한 치환체가 없는 경우에 잘 유지된다. 그림에서 보는 바와 같이 모노나이트로 치환체들인 1A, 1B, 1C에서의 나이트로기는 이미다졸 링과의 이면각이 비틀림없이 거의 0에 가까운 것을 알 수 있다. 또한 다이나이트로에서도 나이트로기가 떨어져 있는 2A에서도 나이트로기의 이면각이 거의 0에 가까운 수준이다. 하지만 2B의 경우는 메틸기를 사이에 두고 나이트로기가 양 옆에 위치하여 구조적인 복잡성을 유발하고 2C의 경우에는 4,5위치에 이웃하는 나이트로기들 간의 상호 척력이 작용하여 심하게는 42°까지의 높은 비틀림이 있음을 볼 수 있다. 트라이나이트로 유도체인 3에서는 모든 나이트로기들이 이미다졸 고리로부터 적게는 20°, 많게는 36°까지의 비틀림이 있음을 확인할 수 있었다. 이들 비틀림으로 인해 분자의 내부 응력이 높아져 상대적으로 높은 에너지를 가지게 되어 생성열이 높이지며, 분자 간의 배열에도 영향을 미쳐 밀도에 영향이 있을 것으로 분석된다.

Fig. 2.

Optimized molecular structures of methylnitroimidzole derivatives calculated at B3LYP/6-31G*. The values in parentheses are dihedral angles of nitro groups twisted from the imidazole ring.

3.2. 화학적 특성 분석

안정된 삼차원 분자구조들을 획득하면 이들 분자구조에서 화학적 특성치들을 계산한다. 화약적 특성을 정확하게 예측하려면 고체상태에서의 생성열과 밀도를 정확하게 예측하여야 한다. 상기 기술한 방법에 따라 분석 작업을 진행하였으며 각 결과는 Table 1에 정리한 바와 같다. 그리고 생성열과 밀도의 나이트로기가 치환된 위치와 개수에 따른 변화를 확인하기 위해 Fig. 3과 Fig. 4에 도시하였다.

Table1.

Chemical formula, oxygen balance, heat of formation (in solid phase), and density of methylnitroimidazole derivatives^a

^aValues in parentheser are from the literatures. ^bRef. 14. ^cFrom Cambridge crystal structure database, refcodes RUGJUZ for 2C, VAHFIU for 3.

Fig. 3.

Heats of formation of methylnitroimidazole derivatives.

Fig. 4.

Densities of methylnitroimidazole derivatives.

본 계산으로 진행한 생성열은 모노나이트로 유도체들에서는 12.1-14.9 kcal/mol로 계산이 되며, 다이나이트로 유도체에서는 13.2-18.5 kcal/mol로 계산이 된다. 트라이나아트로 유도체인 3에서는 24.8 kcal/mol로 분석된다. 각 나이트로기가 모노나이트로에서 다이나이트로, 다이나이트로에서 트라이나이트로로 추가적으로 치환됨에 따라 약간씩 증가하는 현상을 Fig. 3에서 볼 수 있으나 증가폭은 미미한 편으로 나타나고 있다. 같은 나이트로기가 치환된 이성질체들의 생성열 변화도 매우 미미한 것으로 나타나고 있다. 하지만 다이나이트로 유도체들의 이성질체에서는 나이트로기의 비틀림 현상에 따른 응력 증가가 그대로 생성열 값에 반영되고 있는 것으로 보여진다. 즉 Fig. 2에서 보는 바와 같이 나이트로기의 비틀림 현상이 있는 2B, 2C의 생성열 값이 비틀림 현상이 전혀 없는 2A보다 3-5 kcal/mol 정도 높은 것으로 분석된다. 본 계산 결과는 Su 등이 분석한 수치에 비해 개략적으로 15-25 kcal/mol 낮은 것으로 분석되었다.¹⁵ 이는 Su 등이 유사한 계산 작업을 진행해서 보고한 내용으로 기체 상태에서의 생성열을 산출한 것으로 판단된다. 본 계산에서 산출한 승화열의 예측 결과가 18-27 kcal/mol로 이들 값을 기체 상태에서의 생성열 수치에서 제외하였을 경우 본 계산 결과와 거의 일치하는 것으로 추정되므로 본 연구에서의 생성열 수치는 매우 적절한 것으로 판단된다.

분자 표면 정전기 전위 계산으로부터 산출한 밀도 값들은 모노나이트로 유도체의 경우는 1.48-1.58 g/cc로 계산이 되며, 다이나이트로 유도체에서는 1.66-1.73 g/cc로 증가한다. 나이트로기가 이미다졸 링이 1개 치환된 경우와 2개 치환된 차이에서 밀도가 상당히 증가하는 것으로 확인할 수 있었다. 또한 나이트로기가 3개 치환된 트라이나이트로이미다졸에서는 1.83 g/cc로 상당히 높은 밀도를 가지게 된다. 본 계산을 수행한 후 캠브리지 결정구조를 조사하여 기존에 2C와 3에 대해서는 결정 밀도가 알려져 있음을 확인하였다. 이들 값들과 본 연구에서의 결과치를 비교한 결과 차이가 0.1 g/cc 이내로 매우 잘 일치하고 있음을 확인할 수 있었다.

3.3. 화약적 성능 분석

밀도함수 이론으로 분석된 밀도와 생성열을 입력자료로 사용하여 분석한 화약 특성치인 C-J 압력과 폭발속도를 Table 2에 정리하였다. 계산에는 결정 밀도의 97%를 실질적인 밀도로 재산정하여 치타 프로그램에서 분석 작업을 수행하였다. 97%의 실질적인 밀도를 사용함으로써 밀도에 매우 민감하게 변화하는 화약 분자의 성능을 과다하게 추정하는 것을 막을 수 있으며, 다양한 화약 물질들을 추가하여 고체 입자들을 압축 제조하는 화약의 성능 추정에도 유사한 결과를 유지할 수 있다. C-J 압력과 폭발속도 계산 결과들은 Fig. 5에 도시하였으며 군용 화약으로 널리 사용되는 TNT, RDX, HMX 등의 값들과 비교하였다. 치타 프로그램으로 분석된 모노나이트로 유도체들의 C-J 압력은 12.6-14.9 GPa 정도로 추정된다. 이 수치는 TNT의 C-J 압력인 21.0 GPa에 상당히 떨어지는 수치이다. 반면 다이나이트로 유도체의 C-J 압력은 12.4-23.5 GPa로 TNT 보다는 약간 높은 수치이지만 예측오차를 생각한다면 거의 유사한 수준으로 판단된다. 메틸트라이나이트로이미다졸의 경우 31.6 GPa로 분석되며 이 수치는 RDX 수치인 33.8 GPa에 거의 도달하는 수치로 나이트로기가 추가할수록 성능이 높이짐을 알 수 있다. 폭발속도를 사용하여 화약 성능들을 분석하면 모노나이트로 유도체들의 폭발속도는 16.12-6.48 km/sec로 추정되며 TNT의 폭발속도인 6.93 km/sec에 훨씬 못 미치는 수치로 분석된다. 반면 다이나이트로 유도체들의 폭발속도는 7.33-7.53 km/sec로 TNT 성능 보다는 우수한 것으로 예측된다. 3은 8.47 km/sec로 RDX의 폭발속도인 8.70 km/sec에 거의 육박하는 수치이다. TNT, RDX 등과의 비교적 분석은 폭발속도로 추정한 비교치가 약간 높은 것으로 나타나지만 전반적인 경향성은 매우 유사하게 나타난다. Ravi 등이 추정한 3의 C-J 압력과 폭발속도는 33.9 GPa과 8.66 km/sec로 본 연구에서 수행한 결과와 잘 일치하는 수치이다.¹⁶ 이들의 수치가 약간 과다하게 추정되는 것은 밀도 추정의 문제이거나 화약 성능 추정을 단순한 Kamlet-Jacobs 식을 이용하여 추정한 것이므로 오차가 발생했을 것으로 판단된다.

Table2.

C-J pressure^a and detonation velocity^a of methylnitroimidazole derivatives^b

^aEstimated at Cheetah 2.0 program. Estimation performed with the density values assumed to be packed with 97% of crystal density. ^bValues in parentheses are from the literatures. ^cRef. 16

Fig. 5.

C-J pressures (top) and detonation velocities (bottom) of methylnitroimidazole derivatives.