1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
# fun collision of: 2 H2 + O2 -> 2 H2O
import os
from ase.build import molecule
from ase.calculators.dftb import Dftb
from ase.io import read, write
from ase.io.dftb import read_dftb_velocities, write_dftb_velocities
o2 = molecule('O2')
h2_1 = molecule('H2')
h2_2 = molecule('H2')
o2.translate([0, 0.01, 0])
h2_1.translate([0, 0, 3])
h2_1.euler_rotate(center='COP', theta=90)
h2_2.translate([0, 0, -3])
h2_2.euler_rotate(center='COP', theta=90)
o2.set_velocities(([0, 0, 0], [0, 0, 0]))
h2_1.set_velocities(([0, 0, -3.00], [0, 0, -3.000]))
h2_2.set_velocities(([0, 0, 3.000], [0, 0, 3.000]))
atoms = o2 + h2_1 + h2_2
# 1fs = 41.3 au
# 1000K = 0.0031668 au
calculator_NVE = Dftb(
label='h2o',
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_O='p',
Hamiltonian_MaxAngularMomentum_H='s',
Driver_='VelocityVerlet',
Driver_MDRestartFrequency=10,
Driver_Velocities_='',
Driver_Velocities_empty='<<+ "velocities.txt"',
Driver_Steps=1000,
Driver_KeepStationary='Yes',
Driver_TimeStep=4.13,
Driver_Thermostat_='None',
Driver_Thermostat_empty='',
)
# 1fs = 41.3 au
# 1000K = 0.0031668 au
calculator_NVT = Dftb(
label='h2o',
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_O='p',
Hamiltonian_MaxAngularMomentum_H='s',
Driver_='VelocityVerlet',
Driver_MDRestartFrequency=5,
Driver_Velocities_='',
Driver_Velocities_empty='<<+ "velocities.txt"',
Driver_Steps=500,
Driver_KeepStationary='Yes',
Driver_TimeStep=8.26,
Driver_Thermostat_='Berendsen',
Driver_Thermostat_Temperature=0.00339845142, # 800 degC
Driver_Thermostat_CouplingStrength=0.01,
)
write_dftb_velocities(atoms, 'velocities.txt')
atoms.calc = calculator_NVE
atoms.get_potential_energy() # run NVE ensemble using DFTB+'s own driver
atoms = read('geo_end.gen')
write('after_NVE.xyz', atoms)
read_dftb_velocities(atoms, filename='geo_end.xyz')
write_dftb_velocities(atoms, 'velocities.txt')
os.system('mv geo_end.xyz geo_end_NVE.xyz')
atoms.calc = calculator_NVT
atoms.get_potential_energy() # run NVT ensemble using DFTB+'s own driver
atoms = read('geo_end.gen')
write('after_NVT.xyz', atoms)
read_dftb_velocities(atoms, filename='geo_end.xyz')
write_dftb_velocities(atoms, 'velocities.txt')
os.system('mv geo_end.xyz geo_end_NVT.xyz')
# to watch:
# ase gui geo_end_NVE.xyz geo_end_NVT.xyz
|
