-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpimcase_table_npt
executable file
·186 lines (175 loc) · 6.96 KB
/
pimcase_table_npt
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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
#!/usr/bin/env python
import numpy as np
import mcase.particle_moves as mcp
import mcase.cell_moves as mcc
import mcase.quantum_moves as mcq
from mcase.table import read_table, evaluate
from ase.io import read, write
from ase.units import kB, GPa
from mcase.reader import read_input
from sys import argv
settings = read_input(argv[1])
rng = np.random.default_rng()
# setup system
walker = read(settings['xyz'], index=':')
if not isinstance(walker, list):
print('walker must contain more than one slice')
exit()
# determine what the calculator is
if settings['calculator'] == 'nep':
from pynep.calculate import NEP
calc = NEP(settings['model'])
elif settings['calculator'] == 'dp':
from deepmd.calculator import DP
calc = DP(model=settings['model'])
elif settings['calculator'] == 'nequip':
from nequip.ase import NequIPCalculator
calc = NequIPCalculator.from_deployed_model(settings['model'],\
device='cuda')
elif settings['calculator'] == 'mace':
from mace.calculators import MACECalculator
calc = MACECalculator(model_path=settings['model'],\
device='cuda', default_dtype='float32')
else:
print('calculator not recognized')
exit()
# setup tabulated pair potential
pair_energy, pair_force = read_table(settings['table'])
# setup simulation
temperature = float(settings['temperature'])
pressure = float(settings['pressure']) * GPa
timestep = float(settings['timestep'])
cellstep = float(settings['cellstep'])
number_samples = int(settings['samples'])
beta = 1.0/(kB*temperature)
scale = np.sqrt(2*timestep)
save_interval = int(settings['save_every'])
# path-integral scale
num_beads = len(walker)
if 'mass' in settings:
mass = float(settings['mass'])
else:
mass = walker[0].get_masses()[0]
path_time = beta / num_beads / mass
path_time *= 0.00418015929
# successful read, print out stuff
print('read inputs from ' + argv[1])
print('particles read from ' + settings['xyz'])
print('they have mass %s (amu)' % mass)
print(settings['calculator'] + ' model read from ' + settings['model'])
print('tabulated pair potential read from %s' % settings['table'])
print('temperature: %s K (beta = %s eV^-1)' % (temperature, beta))
print('pressure: %s eV/A^3 (%s GPa)' % (pressure, settings['pressure']))
print('timestep: %s (length scale = %s)' % (timestep, scale))
print('cell step: %s A' % cellstep)
print('quantum scale: %s (%d beads)' % (path_time, num_beads))
print('number of samples to calculate: %d' % number_samples)
print('writing samples from bead 0 to ' + settings['save_slice'])
print('will write every %d steps' % save_interval)
print('will write final path to ' + settings['save_xyz'])
# setup
N = len(walker[0])
r = np.zeros((num_beads, N, 3))
energy_c = np.zeros(num_beads)
forces_c = np.zeros_like(r)
for j in range(num_beads):
r[j] = walker[j].get_positions()
walker[j].calc = calc
energy_c[j] = walker[j].get_potential_energy()
forces_c[j] = walker[j].get_forces()
# add table
pairs = walker[j].get_all_distances(mic=True, vector=True)
x, y = evaluate(pairs, pair_energy, pair_force)
energy_c[j] += x
forces_c[j] += y
cell = walker[0].cell.cellpar()[:3] # CELL MUST BE IDENTICAL OVER PATH
volume = np.prod(cell)
rc = mcq.get_rc(r)
energy_q = mcq.spring_potential(r, beta, path_time)
forces_q = mcq.spring_force(r, rc, beta, path_time)
# begin sampling
total_accepted = 0
cell_accepted = 0
energy_c_new = np.zeros_like(energy_c)
forces_c_new = np.zeros_like(forces_c) # keeping these arrays better for memory?
for i in range(number_samples):
path_accepted = 0
# particle move
for j in range(num_beads):
r_old = np.copy(r[j])
rc_old = np.copy(rc)
r[j] += rng.normal(size=(N,3), scale=scale)
r[j] += beta * timestep * forces_q[j]
r[j] += beta * timestep * forces_c[j] / num_beads
# calculate new quantities
rc = mcq.get_rc(r)
energy_q_new = mcq.spring_potential(r, beta, path_time)
forces_q_new = mcq.spring_force(r, rc, beta, path_time)
walker[j].set_positions(r[j])
energy_c_new[j] = walker[j].get_potential_energy()
forces_c_new[j] = walker[j].get_forces()
pairs = walker[j].get_all_distances(mic=True, vector=True)
x, y = evaluate(pairs, pair_energy, pair_force)
energy_c_new[j] += x
forces_c_new[j] += y
# accept / reject
if rng.random() < mcp.acceptance_ratio(r_old, r[j],\
energy_q + energy_c[j] / num_beads,\
energy_q_new + energy_c_new[j] / num_beads,\
forces_q[j] + forces_c[j] / num_beads,\
forces_q_new[j] + forces_c_new[j] / num_beads,\
beta, timestep):
energy_q = energy_q_new
forces_q = forces_q_new
energy_c[j] = energy_c_new[j]
forces_c[j] = forces_c_new[j]
path_accepted += 1
else:
r[j] = np.copy(r_old)
rc = np.copy(rc_old)
walker[j].set_positions(r[j])
# cell move
cell_new = cell + cellstep * (rng.random(size=3) - 0.5)
volume_new = np.prod(cell_new)
r_new = mcq.rescale_path(r, cell, cell_new)
energy_q_new = mcq.spring_potential(r_new, beta, path_time)
rc_new = mcq.get_rc(r_new)
forces_q_new = mcq.spring_force(r_new, rc_new, beta, path_time)
# have to set all slices
for j in range(num_beads):
walker[j].set_cell(cell_new, scale_atoms=True)
energy_c_new[j] = walker[j].get_potential_energy()
forces_c_new[j] = walker[j].get_forces()
pairs = walker[j].get_all_distances(mic=True, vector=True)
x, y = evaluate(pairs, pair_energy, pair_force)
energy_c_new[j] += x
forces_c_new[j] += y
if rng.random() < mcc.old_acceptance_ratio(\
energy_q + np.mean(energy_c),\
energy_q_new + np.mean(energy_c_new),\
volume, volume_new, N*num_beads, pressure, beta):
cell = np.copy(cell_new)
volume = volume_new
r = np.copy(r_new)
energy_q = energy_q_new
rc = np.copy(rc_new)
forces_q = np.copy(forces_q_new)
energy_c = np.copy(energy_c_new)
forces_c = np.copy(forces_c_new)
cell_accepted += 1
else:
for j in range(num_beads):
# need to reset walkers
walker[j].set_cell(cell, scale_atoms=True)
print(np.mean(energy_c), energy_q, cell[0], cell[1], cell[2])
total_accepted += path_accepted / num_beads
if i % save_interval == (save_interval - 1):
write(settings['save_slice'], walker[0], format='extxyz', append=True,
write_info=False, write_results=False)
write(settings['save_xyz'], walker, format='extxyz', append=False,
write_info=False, write_results=False)
# finished
print('particle move acceptance rate: %.2f' % (total_accepted / number_samples))
print('cell move acceptance rate: %.2f' % (cell_accepted / number_samples))
write(settings['save_xyz'], walker, format='extxyz', append=False,
write_info=False, write_results=False)