# Setup a quartz cell.
lattice = Hexagonal(4.916*Angstrom, 5.4054*Angstrom)

elements = [Silicon, Silicon, Silicon, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen,
            Oxygen]

fractional_coordinates = [[ 0.4697        ,  0.            ,  0.            ],
                          [ 0.            ,  0.4697        ,  0.666666666667],
                          [ 0.5303        ,  0.5303        ,  0.333333333333],
                          [ 0.4135        ,  0.2669        ,  0.1191        ],
                          [ 0.2669        ,  0.4135        ,  0.547567      ],
                          [ 0.7331        ,  0.1466        ,  0.785767      ],
                          [ 0.5865        ,  0.8534        ,  0.214233      ],
                          [ 0.8534        ,  0.5865        ,  0.452433      ],
                          [ 0.1466        ,  0.7331        ,  0.8809        ]]

bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# Repeat the primitive cell to make a larger configuration.
bulk_configuration = bulk_configuration.repeat(3, 3, 3)

# Make the amorphous SiO2 cell with a density of 2.2 gram/cm**3.
bulk_configuration = amorphousPrebuilder(bulk_configuration, density=2.2*gram/cm**3)

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

potentialSet = Tersoff_SiO_2007()
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setVerletListsDelta(0.25*Angstrom)

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('amorphous_sio2.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = MaxwellBoltzmannDistribution(
    temperature=2000.0*Kelvin,
    remove_center_of_mass_momentum=True
)

method = NVTNoseHoover(
    time_step=1*femtoSecond,
    reservoir_temperature=2000*Kelvin,
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    chain_length=3,
    initial_velocity=initial_velocity,
)

constraints = [FixCenterOfMass()]

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=constraints,
    trajectory_filename='amorphous_sio2.hdf5',
    steps=10000,
    log_interval=100,
    method=method
)

bulk_configuration = md_trajectory.lastImage()

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = ConfigurationVelocities(
    remove_center_of_mass_momentum=True
)

method = NPTMartynaTobiasKlein(
    time_step=1*femtoSecond,
    reservoir_temperature=1500*Kelvin,
    reservoir_pressure=0*bar,
    thermostat_timescale=100*femtoSecond,
    barostat_timescale=500*femtoSecond,
    initial_velocity=initial_velocity,
    heating_rate=0*Kelvin/picoSecond,
)

constraints = [FixCenterOfMass()]

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=constraints,
    trajectory_filename='amorphous_sio2.hdf5',
    steps=50000,
    log_interval=100,
    method=method
)

bulk_configuration = md_trajectory.lastImage()