Future Work

Hydrogen hydrate C2 has given research interest in the field of energy storage. In this material compound, hydrogen content can reach 10 % wt which is suitable for energy storage material. However, to synthesize this filled ice, the high-pressure condition is needed. Several published data revealed fast diffusion with less affected by pressurization is found under this highly compressed system. The material should be able to be generated under the more moderate thermodynamic condition to become available in a possible usage state. Also, it is the necessity that the product is stable for an extended time under ambient thermodynamic state concerning mobility and depository.

Kumar et al. [99] have synthesized hydrogen filled ice Ic from amorphous ice in

low-pressure condition less than 45 MPa. The process was performed by conditioning the amorphous ice in hydrogen gas pressure at low temperature. Subsequently, the annealing process to the higher temperature was implemented. The small amount of SII hydrate is then formed along with the crystalline structure of hydrate C2 as the primary product.

Additionally, other fascinating points are the availability and stability at ambient pres-sure and low 77 K temperature. Despite less H2 content to make this material under a possible condition to be the energy storage candidate, dynamic properties could be valu-able to study. For instance, diffusion process of guest hydrogen molecules is associating with the time needed to complete the synthesize. Also, conversion kinetic gives more properties to be considered in the study.

Quantum effect of hydrogen guest dynamics is expected to be evident under this low temperature. It has been reported in hexagonal ice that quantum effect gives noticeable contrast in lattice vibrations and potential energy in comparison to the classically cal-culated counterpart [100]. Also, a molecular interaction governed by covalent bond of water when interacting with others shifted the phase-space boundary [101]. This effect also explained the inter-site hopping dynamics in clathrate hydrate structure dependable on temperature and lattice site accommodation number [102].

To further investigate the process, the long calculation is needed while revealing the quantum effect. Concerning the necessary features to be studied, Path Integral Molecular Dynamics calculation is suitable to be used. The method is implemented by using the isomorphism between the real quantum particle and the set of classical ring polymers.

The formalism is as follows. The real quantum particle is represented by L replicas of particles in imaginary time sequences. These duplicates known as beads and connected to each other by a harmonic spring consecutively to form a ring polymer. The number of beads corresponds to the Trotter number while the ring to a cyclic quantum path.

The classical equation of motion is designed in a such a way that the resulted partition function can be formally equivalent to the quantum set. The method is proven influential to study many-body quantum system.

A

Mathematical Backgrounds of MD

In this section, molecular dynamics simulation is briefly explained. Classical MD method is chosen to be used in these two filled ice system’s calculation. Since the system under the thermodynamic condition where the quantum effect is negligibly small, the classical approach is sufficient to be used in the simulation of this condensed phase system. The molecular dynamics scheme is allowing one to proceed the extended period of simulation time with reasonable computational cost. Also, the method is proven to be statistically reliable in producing the required trajectory data for the study of a diffusion process.

LAMMPS molecular dynamics simulation package was used to simulate the system.

This MD package offers an algorithm that can be used for large bulky system enhanced by the ability to be done in multi-platform parallel processes. The program can also deliver fast and accurate computation for extensive MD work. Within this dissertation document, several algorithms that LAMMPS closely followed are briefly derived and ex-plained. The discussion is limited to only those that were used in this particular calcu-lation. Subsequently, post ”production run” analysis was done with additional self-build programs based on a commonly-used algorithm for the specific intention. The codes are only suitable for this particular purpose.

The derivation and discussion of the algorithm are started from periodic boundary condition scheme. The algorithm is suitable to be used in the case of orthogonal and non-orthogonal periodic cells calculation. In the next discussion section, the mathemati-cal background of long-range energy contributions is explained historimathemati-cally, starting from the first proposed algorithm until the one that closely-followed by LAMMPS recently.

Subsequently, the attached thermostat and barostat method are discussed. These algo-rithms are presented in historical sequences, in which similar to the previous algoalgo-rithms discussion. Finally, the last discussion section in this mathematical background is the constraints algorithm used to incorporate rigid molecule model in the calculation. This final section will include a summary of SHAKE and RATTLE algorithm complemented with general stages concerning the implementation.

A.1 Periodic Boundary Condition

In post-production stage, the periodic boundary condition is calculated by using the fol-lowing procedure. Suppose that the unit cell of the system is

a=axˆi+ayˆj+azkˆ (A.1) b=bxˆi+byˆj+bzkˆ (A.2) c=cxˆi+cyˆj+czk,ˆ (A.3) in which can be represented in the matrix form as

h=

⎡

⎢⎣

ax bx cx

ay by cy

az bz cz

⎤

⎥⎦. (A.4)

The particle’s displacement vector is written by

dr=dxˆi+dyˆj+dzk.ˆ (A.5) The periodic condition is applied by taking the matrix operation as follows. Let

ht =h^t, (A.6)

then the matrix of periodic boundary condition can be written as

Bpbc=h⁻¹_t

⎡

⎢⎣ dx dy dz

⎤

⎥⎦. (A.7)

By taking B_pbc = B_pbc one can obtain the displacement vector subtracted from the periodic boundary condition as

⎡

⎢⎣ dx dy dz

⎤

⎥⎦=h Bpbc. (A.8)