-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpyaimd.py
135 lines (110 loc) · 4.23 KB
/
pyaimd.py
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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
from pyscf import gto
from pyscf import scf
from pyscf import grad
import numpy as np
import md_func as aimd
import constants as const
# Read molecular structure from "H2O.xyz"
mol = gto.Mole()
mol.atom = "H2O.xyz"
mol.basis = "sto-3g"
mol.build()
# Calculate energy as potential
mf = scf.RHF(mol)
Epot = mf.kernel()
# Calculate gradient
g = mf.Gradients()
grad = g.kernel()
atom_mass = np.array(mol.atom_mass_list(), dtype=float)
mass_tot = np.sum(atom_mass) # Store mass in units of dalton
natm = np.size(atom_mass)
# Store elemental symbol
elem = [None]*natm
for i in range(natm):
elem[i] = mol.atom_symbol(i)
# Initialize velocity
# Generate random numbers with mean 0, variance 1, following a normal distribution
init_temp = 298.15
vel_old = np.random.normal(0, 1.0, size=np.prod((natm,3))).reshape(natm,3)
# Subracting the average velocity so the center of mass does not move
vel_avg = np.zeros(3, dtype=float)
for i in range(3):
vel_avg[i] = np.sum(vel_old[:,i])/natm
for i in range(natm):
vel_old[i,:] = vel_old[i,:] - vel_avg[i]
# Maxwell-Boltzmann
for i in range(natm):
vel_old[i,:] = np.sqrt((init_temp*const.kB*const.J2au) \
/(atom_mass[i]*const.dalton2au))*vel_old[i,:]
pos_old = np.array(mol.atom_coords(), dtype=float) # Store initial coordinates in units of Bohr
acc_old = np.zeros((natm,3), dtype=float)
pos_new = np.zeros((natm,3), dtype=float)
vel_new = np.zeros((natm,3), dtype=float)
acc_new = np.zeros((natm,3), dtype=float)
# Write coordinates in pos.xyz
with open('pos.xyz', 'w') as f:
f.write(f'{natm}\n')
f.write('nstep = 0\n')
for symbol, coords in zip(elem, pos_old):
f.write(f'{symbol} {coords[0]*const.ang2bo:> .6f} {coords[1]*const.ang2bo:> .6f} {coords[2]*const.ang2bo:> .6f}\n')
# Time step(fs) and number of steps
dt = 0.5 * const.fs2au
nstep = 1000
# Temperature of the thermostat (K), time constant (fs)
bath_temp = 298.15
con_time = 30.0 * const.fs2au
# Calculate temperature and kinetic
Ekin = aimd.cal_kin(natm,vel_old,atom_mass)
temp_cal = aimd.cal_temp(natm,Ekin)
# Write AIMD data in aimd.dat
with open('aimd.dat','w') as f:
f.write(f' step Ekin Epot Etot temp\n')
f.write(f'===== ============ ============== ============== ======\n')
f.write(f'{0:5d} {Ekin:.10f} {Epot:.10f} {Epot+Ekin:.10f} {temp_cal:.2f}\n')
for n in range(nstep):
# Update coordinates with Verlet velocity
for i in range(natm):
acc_old[i] = - grad[i]/(atom_mass[i]*const.dalton2au)
pos_new[i,:] = pos_old[i,:] + vel_old[i,:] * dt \
+ 0.5 * acc_old[i,:] * (dt**2.0)
# Update coordinates in PySCF
mol = gto.Mole()
mol.atom = ''
for sym, coords in zip(elem, pos_new):
x, y, z = coords[:3]
mol.atom += f'{sym} {x} {y} {z}; ' # Note the format: 'symbol x y z;'
mol.unit = 'B'
mol.build()
# Calculate energy and gradients
mf = scf.RHF(mol)
Epot = mf.kernel()
g = mf.Gradients()
grad = g.kernel()
# Update velocity with Verlet velocity
for i in range(natm):
acc_new[i] = - grad[i]/(atom_mass[i]*const.dalton2au)
vel_new[i,:] = vel_old[i,:] \
+ 0.5 * (acc_new[i,:]+acc_old[i,:]) * dt
# Remove translational motion
mv = np.zeros((3), dtype=float)
vel_cent = np.zeros((3), dtype=float)
for i in range(natm):
mv[:] = mv[:] + vel_new[i,:]*(atom_mass[i]*const.dalton2au)
vel_cent[:] = mv[:]/(mass_tot*const.dalton2au)
for i in range(natm):
vel_new[i,:] = vel_new[i,:] - vel_cent[:]
# Thermostat:
Ekin = aimd.cal_kin(natm,vel_old,atom_mass)
temp_cal = aimd.cal_temp(natm,Ekin)
f = aimd.berendsen(dt,bath_temp,con_time,temp_cal)
vel_new = vel_new * f
# Update pos.xyz and aimd.dat
with open('pos.xyz', 'a') as f:
f.write(f'{natm}\n')
f.write(f'nstep = {n+1}\n')
for symbol, coords in zip(elem, pos_new):
f.write(f'{symbol} {coords[0]*const.ang2bo:> .6f} {coords[1]*const.ang2bo:> .6f} {coords[2]*const.ang2bo:> .6f}\n')
with open('aimd.dat','a') as f:
f.write(f'{n+1:5d} {Ekin:.10f} {Epot:.10f} {Epot+Ekin:.10f} {temp_cal:.2f}\n')
pos_old = pos_new
vel_old = vel_new