- Email: [email protected]
Intermolecular force models and the crystal structure of carbon disulphide
CHEMICAL PHYSICS LETTERS
Volume 103, number 2
INTERMOLECULAR
FORCE
MODELS
AND THE CRYSTAL
13 December 1983
STRUCTURE
OF CARBON
DISULPHIDE
Roger W. IMPEY and Michael L. KLEIN Division of Chemistry, National Research Council of Gmada. Ottawa*Ontario, Gmada KlA OR6 Received 16 August 1983; in final form 15 October 1983
Two recently proposed intermolecular potential mod& for CSe are tested by using the new molecular dynamics (N\lD) computer simulation technique that for a given temperature and pressure determines the stable crystal structure of each model. Comparison with the experimental structure of CSe, and its known variation with temperature, reveals deficiencies in both models.
Table 1 Potential parameters of CSZ. vap(R) = 4e,~[(euplR)” (u&)~ I
There has been considerable interest in the properties of solid and liquid CS2 [l-8] _Much of this interest has been focused on the construction of simple intermolecular force models that can explain both the observed crystal structure, that remains orthorhombic (Cmca, with four molecules per unit cell) up to its
crs
melting point, as well as the crystal dynamics_ In a recent article Grout and Leech [l] compared the calculated phonon frequencies of four specific models (called A, B, C and D). These calculations employed harmonic lattice dynamics and used as input the ex-
0 009-2614/83/S 03.00 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B-V.
A)
hiodel C (dcs = 1.57 X) e (k)
0 (A)
E (I;)
0 (.Q
cc
544.4
3.4225
51.10
3.35
cs
289.9 12.93
3.2025 42818
96.60 183.00
3.44 3.52
ss
perimental orthorhombic cell constants a, b, c appropriate to 150 K [7]. Model A [6] was found to be the worst and will not concern us further. Model B was designed to fit crystal data at 0 K [Z] , whereas model C was based in part on molecular dynamics studies of liquid CS, [3,4]. Grout and Leech [ I,?] found a strong correlation between the predictions of models B and C, which are both based upon Lermard-Jones 12-6 atom-atom potentials. The potential parameters for these two models differ widely, however (see table 1). The phonon frequencies, which fitted experiment with mean absolute errors of 8 and 9% respectively for models B and C, were unable to distinguish between these models. Model D, which included an electrostatic quadrupole moment, was able to fit the phonon data only slightly better than models B and C. This finding reinforced the view that harmonic lattice dynamics alone cannot be used to identify realistic potential models [ 11. In the calculations of Grout and Leech [l]
bfodel B (dcs = lJ46
-
the comparison of the four models was carried out at T = 150 K, a temperature very close to the normal melting point of 161.7 K. At such high temperatures anharmonic effects, that are ignored in traditional lattice dynamics calculations [ 1.21, will likely give important contributions to both structural and dynamic properties. Recent advances in computer simulation. in particular, the so called new molecular dynamics (NMD) technique, have made it possible to investigate directly the relationship between crystal structure and intermolecular potentials such at temperatures [9-l 11. The NMD method for molecular systems has been described in the literature so we give only the essential details of our own calculations. In the present application the MD cell consisted of 108 CS, molecules initially arranged as 3 X 3 X 3 unit cells of the observed orthorhombic crystal structure, periodic boundary conditions being used to simulate the infmite crystal. The equations of motion were integrated by standard tech143
CHEMICAL
Volume 163, number 2
PHYSICS LETFERS
tially P, was set at 3.0 kbar, a value which yields the observed crystal density at 150 K, but not the observed
niques using a time step of 2 fs. In the NMD method the MD cell, defined by the lengths Ln = 3a, Lb = gb, L, = 3c and the angles ob, MT,be, is ailowed to change both its size and shape in response to any net imbalance between an externally set pressure P, and the spori-taneously generated internal stresses. In models B and C the CS, molecules interact via Lennard-Jones potentials, the parameters of which are given in table 1. Ini-
0
4
23 December 1983
cell constants. The effect on the MD cell of reducing the temperature at constant P, for model C is shown in fig. I_ Several points should be noted. Firstly, the effect of including thermal pressure, due to translational and rotational motion of the molecules, completely alters the cell dimension compared with those
6
l2 time (&s>
20
Fig- l- Time evolution of quantities charaderising the basic MD cell of C&- in its orthorhombic p&se for model h~Dc&kngth&>L,>,?&.
144
C.(RW
that the
Volume 103. number 2
determined using the more traditional energy minimisation of the static lattice [ 1,2]. Secondly, the detailed temperature behaviour of the MD cell parameter L,, which experimentally shows a negative thermal expansion [S] is not correct. In the real crystal, the anomalous thermal expansion along c can be linked to the temperature dependence of the angle a,, between the CS2 molecular axis and the crystal axes (Y(ol = a,b,c). Experimentally, the angle 6, increases with temperature, while the corresponding angle Bb decreases. This results in the net thermal expansion along c being negative. In the NMD calculations, model C yields 8, = 42” at T = 130 K, which is somewhat smaller than the experimental value of 44” [7], but unfortunately ~9~exhibits the wrong temperature dependence and hence there is no negative expansion along c. Fig. 2 shows an instantaneous configuration generated for model C at T = 84 K, while details of the three molecular dynamics simulations are given in table 2. Accordingly, we conclude that model C, that was parameterised to data in the liquid, is not particularly successful for the solid. Similar calculations were attempted for model B. In
T
Fii. 2. An instantaneous confiwration generated for model C at T= 84 K, projected onto the bc plane. The shaded molecules are displaced by a/2. (La= 18.7 A. Lb = 17.3 -9, Lc = 26.3 A.)
1983
Table 2 Details of the NMD caJcuJations using model C [ 3]_ The last row gives the orientation angle between the molecular mis and the b direction T (W
Pex Ocbx) vol. (cm3 mof’) U (kJ mol-r) At/ (kJ mol-r) a) P &bar) Q (A) b (A) c (A) eb (deg)
130
84
57
3 47.5 -33.6 -3.6 1.47 6.28
3 46.6 -34.8 -3.7 1.41 6.23
3 46.1 -35.5 -3.7 1.38 6.21
5.67
5.65 8.81 41-2
5.64 8.73 40.8
8.87 41.9
a) Au is the correction due to the truncation of the potential at Rcut = 8.0 X. Both U and P have been corrected for this effect. this case starting
C -
23 December
CHEMICAL PHYSICS LETTERS
from
the experimental
lattice
constants
and at T = 150 K, the cell immediately transformed to a radically different structure, similar to that shown in fig. 3, indicating that this model is a poor one for solid CS2. Several additional attempts were made to assess whether or not the experimentally observed structure was metastable with respect to the structure shown in fig. 3. For example, we found that starting at P, = 0.5 kbar and T = 20 K and MD cell with dimensions close to that of the experimental structure remained stable for 1-2 ps, but invariably eventually transformed to a structure similar to that shown in fig. 3. This latter structure whichseems to be monoclinic or triclinic, probably arises because of the artifically enhanced CS interaction in model B with respect to that of model C. It should be noted in passing that traditional constant-volume (fured cell geometry) molecular dynamics cakulations carried out using both models B and C, produced very similar results for the orientation of the CS, molecules with respect to the b axis to those reported by Grout and Leech [2] _We conclude therefore that the structure reported by Grout and Leech [2] for their model B is metastable with respect to that shown in fig. 3. In summary, we fmd that the NMD method, unlike traditional lattice dynamics, is able to distinguish between two recently proposed models for CS2_ It is clear 145
Volume 103, number 2
CHEMICAL PHYSICS LETTERS
23 December 1983
finding is in agreement with the recent work of Burgos and Righini [S] who proposed that suitable models for CS, must not only include anisotropic electrostatic interactions, including hexadecapole-hexadecapole components, but also anisotropic repulsive forces. It remajns to be verified whether or not-the specific model they propose-is consistent with experiment once thermal effects are included in the cakulations.
We thank the referee for some helpful remarks.
References
Fig. 3. An instantamious configuration from the NMD simulation using model B at T= 80 K. projected onto the bc plane. Wructure appears to be monoclinic or triclinic, Ln = 19.1 A, Lb = 115 A&= 395iLL.nb = 9&L2°,nc= 90_3°,bc= 82x?_)
that both models suffer from deficiencies but those of model B are particularly severe. The inadequacy of an atom-atom model for CS, is thus demonstrated. This
146
[l] PJ. Grout and J.W. Leech, J. Phys. Cl5 (1982) L1083. [2] PJ. Grout and J-W. Leech, Mol. Phys. 45 (1982) 51. [3] DJ. Tildesley and P.A. Madden, Mol. Phys. 42 (1981) i137. [4] DJ. Ttidesley and P.A. Madden, Mol. Phys. 48 (1983) 129. [S] E. Burgos and R. Righini, Chem. Phys. Letters 96 (1983) 584. [6] GS. Pawley. G. Dolling. B_M_Powell and B.H. Torrie, Can. J. Phys. 59 (1981) 122. [7] B.M. Powell, G. Dolling and B.H. Torrie. Acta Cryst. B38 (1982) 28. [8] B.?+. Powell. G. Dolling, B_H. Torrie and G.S. Pawley. J. Phys. Cl5 (1932) 4265. [9] M. Paninello and A. Rahman. J. Appl. Phys. 52 (1981) 7182. [lo] S. No& and M.L. Klein, Phys. Rev. Letters 50 (1983) 12077. [ll] S. No& and M-L. Klein, J. Chem. Phys. 78 (1983) 6928.
Volume 103, number 2
INTERMOLECULAR
FORCE
MODELS
AND THE CRYSTAL
13 December 1983
STRUCTURE
OF CARBON
DISULPHIDE
Roger W. IMPEY and Michael L. KLEIN Division of Chemistry, National Research Council of Gmada. Ottawa*Ontario, Gmada KlA OR6 Received 16 August 1983; in final form 15 October 1983
Two recently proposed intermolecular potential mod& for CSe are tested by using the new molecular dynamics (N\lD) computer simulation technique that for a given temperature and pressure determines the stable crystal structure of each model. Comparison with the experimental structure of CSe, and its known variation with temperature, reveals deficiencies in both models.
Table 1 Potential parameters of CSZ. vap(R) = 4e,~[(euplR)” (u&)~ I
There has been considerable interest in the properties of solid and liquid CS2 [l-8] _Much of this interest has been focused on the construction of simple intermolecular force models that can explain both the observed crystal structure, that remains orthorhombic (Cmca, with four molecules per unit cell) up to its
crs
melting point, as well as the crystal dynamics_ In a recent article Grout and Leech [l] compared the calculated phonon frequencies of four specific models (called A, B, C and D). These calculations employed harmonic lattice dynamics and used as input the ex-
0 009-2614/83/S 03.00 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B-V.
A)
hiodel C (dcs = 1.57 X) e (k)
0 (A)
E (I;)
0 (.Q
cc
544.4
3.4225
51.10
3.35
cs
289.9 12.93
3.2025 42818
96.60 183.00
3.44 3.52
ss
perimental orthorhombic cell constants a, b, c appropriate to 150 K [7]. Model A [6] was found to be the worst and will not concern us further. Model B was designed to fit crystal data at 0 K [Z] , whereas model C was based in part on molecular dynamics studies of liquid CS, [3,4]. Grout and Leech [ I,?] found a strong correlation between the predictions of models B and C, which are both based upon Lermard-Jones 12-6 atom-atom potentials. The potential parameters for these two models differ widely, however (see table 1). The phonon frequencies, which fitted experiment with mean absolute errors of 8 and 9% respectively for models B and C, were unable to distinguish between these models. Model D, which included an electrostatic quadrupole moment, was able to fit the phonon data only slightly better than models B and C. This finding reinforced the view that harmonic lattice dynamics alone cannot be used to identify realistic potential models [ 11. In the calculations of Grout and Leech [l]
bfodel B (dcs = lJ46
-
the comparison of the four models was carried out at T = 150 K, a temperature very close to the normal melting point of 161.7 K. At such high temperatures anharmonic effects, that are ignored in traditional lattice dynamics calculations [ 1.21, will likely give important contributions to both structural and dynamic properties. Recent advances in computer simulation. in particular, the so called new molecular dynamics (NMD) technique, have made it possible to investigate directly the relationship between crystal structure and intermolecular potentials such at temperatures [9-l 11. The NMD method for molecular systems has been described in the literature so we give only the essential details of our own calculations. In the present application the MD cell consisted of 108 CS, molecules initially arranged as 3 X 3 X 3 unit cells of the observed orthorhombic crystal structure, periodic boundary conditions being used to simulate the infmite crystal. The equations of motion were integrated by standard tech143
CHEMICAL
Volume 163, number 2
PHYSICS LETFERS
tially P, was set at 3.0 kbar, a value which yields the observed crystal density at 150 K, but not the observed
niques using a time step of 2 fs. In the NMD method the MD cell, defined by the lengths Ln = 3a, Lb = gb, L, = 3c and the angles ob, MT,be, is ailowed to change both its size and shape in response to any net imbalance between an externally set pressure P, and the spori-taneously generated internal stresses. In models B and C the CS, molecules interact via Lennard-Jones potentials, the parameters of which are given in table 1. Ini-
0
4
23 December 1983
cell constants. The effect on the MD cell of reducing the temperature at constant P, for model C is shown in fig. I_ Several points should be noted. Firstly, the effect of including thermal pressure, due to translational and rotational motion of the molecules, completely alters the cell dimension compared with those
6
l2 time (&s>
20
Fig- l- Time evolution of quantities charaderising the basic MD cell of C&- in its orthorhombic p&se for model h~Dc&kngth&>L,>,?&.
144
C.(RW
that the
Volume 103. number 2
determined using the more traditional energy minimisation of the static lattice [ 1,2]. Secondly, the detailed temperature behaviour of the MD cell parameter L,, which experimentally shows a negative thermal expansion [S] is not correct. In the real crystal, the anomalous thermal expansion along c can be linked to the temperature dependence of the angle a,, between the CS2 molecular axis and the crystal axes (Y(ol = a,b,c). Experimentally, the angle 6, increases with temperature, while the corresponding angle Bb decreases. This results in the net thermal expansion along c being negative. In the NMD calculations, model C yields 8, = 42” at T = 130 K, which is somewhat smaller than the experimental value of 44” [7], but unfortunately ~9~exhibits the wrong temperature dependence and hence there is no negative expansion along c. Fig. 2 shows an instantaneous configuration generated for model C at T = 84 K, while details of the three molecular dynamics simulations are given in table 2. Accordingly, we conclude that model C, that was parameterised to data in the liquid, is not particularly successful for the solid. Similar calculations were attempted for model B. In
T
Fii. 2. An instantaneous confiwration generated for model C at T= 84 K, projected onto the bc plane. The shaded molecules are displaced by a/2. (La= 18.7 A. Lb = 17.3 -9, Lc = 26.3 A.)
1983
Table 2 Details of the NMD caJcuJations using model C [ 3]_ The last row gives the orientation angle between the molecular mis and the b direction T (W
Pex Ocbx) vol. (cm3 mof’) U (kJ mol-r) At/ (kJ mol-r) a) P &bar) Q (A) b (A) c (A) eb (deg)
130
84
57
3 47.5 -33.6 -3.6 1.47 6.28
3 46.6 -34.8 -3.7 1.41 6.23
3 46.1 -35.5 -3.7 1.38 6.21
5.67
5.65 8.81 41-2
5.64 8.73 40.8
8.87 41.9
a) Au is the correction due to the truncation of the potential at Rcut = 8.0 X. Both U and P have been corrected for this effect. this case starting
C -
23 December
CHEMICAL PHYSICS LETTERS
from
the experimental
lattice
constants
and at T = 150 K, the cell immediately transformed to a radically different structure, similar to that shown in fig. 3, indicating that this model is a poor one for solid CS2. Several additional attempts were made to assess whether or not the experimentally observed structure was metastable with respect to the structure shown in fig. 3. For example, we found that starting at P, = 0.5 kbar and T = 20 K and MD cell with dimensions close to that of the experimental structure remained stable for 1-2 ps, but invariably eventually transformed to a structure similar to that shown in fig. 3. This latter structure whichseems to be monoclinic or triclinic, probably arises because of the artifically enhanced CS interaction in model B with respect to that of model C. It should be noted in passing that traditional constant-volume (fured cell geometry) molecular dynamics cakulations carried out using both models B and C, produced very similar results for the orientation of the CS, molecules with respect to the b axis to those reported by Grout and Leech [2] _We conclude therefore that the structure reported by Grout and Leech [2] for their model B is metastable with respect to that shown in fig. 3. In summary, we fmd that the NMD method, unlike traditional lattice dynamics, is able to distinguish between two recently proposed models for CS2_ It is clear 145
Volume 103, number 2
CHEMICAL PHYSICS LETTERS
23 December 1983
finding is in agreement with the recent work of Burgos and Righini [S] who proposed that suitable models for CS, must not only include anisotropic electrostatic interactions, including hexadecapole-hexadecapole components, but also anisotropic repulsive forces. It remajns to be verified whether or not-the specific model they propose-is consistent with experiment once thermal effects are included in the cakulations.
We thank the referee for some helpful remarks.
References
Fig. 3. An instantamious configuration from the NMD simulation using model B at T= 80 K. projected onto the bc plane. Wructure appears to be monoclinic or triclinic, Ln = 19.1 A, Lb = 115 A&= 395iLL.nb = 9&L2°,nc= 90_3°,bc= 82x?_)
that both models suffer from deficiencies but those of model B are particularly severe. The inadequacy of an atom-atom model for CS, is thus demonstrated. This
146
[l] PJ. Grout and J.W. Leech, J. Phys. Cl5 (1982) L1083. [2] PJ. Grout and J-W. Leech, Mol. Phys. 45 (1982) 51. [3] DJ. Tildesley and P.A. Madden, Mol. Phys. 42 (1981) i137. [4] DJ. Ttidesley and P.A. Madden, Mol. Phys. 48 (1983) 129. [S] E. Burgos and R. Righini, Chem. Phys. Letters 96 (1983) 584. [6] GS. Pawley. G. Dolling. B_M_Powell and B.H. Torrie, Can. J. Phys. 59 (1981) 122. [7] B.M. Powell, G. Dolling and B.H. Torrie. Acta Cryst. B38 (1982) 28. [8] B.?+. Powell. G. Dolling, B_H. Torrie and G.S. Pawley. J. Phys. Cl5 (1932) 4265. [9] M. Paninello and A. Rahman. J. Appl. Phys. 52 (1981) 7182. [lo] S. No& and M.L. Klein, Phys. Rev. Letters 50 (1983) 12077. [ll] S. No& and M-L. Klein, J. Chem. Phys. 78 (1983) 6928.