-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbasin_hopping.py
218 lines (186 loc) · 9.25 KB
/
basin_hopping.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
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
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
############################################################
# Implementation of the Basin Hopping (BH) algorithm for
# structure optimization.
#
# Version 1.2 November 2016
#
# BH algorithm has been implemented using Python3.4;
# coupled to Quantum Espresso 5.2 (DFT code as calculator)
# It also works with newest QE version, 6.0
#
#
# Note: Output folders will be generated in current directory
##############################################################
import glob
import os
import math
import random
import re
import sys
import datetime
from subprocess import*
# Variable declarations
init = ""
input_dir = ""
programs_dir = ""
substrate_nat = ""
cluster_nat = ""
cluster_ntyp = ""
temperature_K = 0
step_width = ""
iterations = ""
x_range = ""
y_range = ""
z_range = ""
z_vacuum = ""
time_initial = 0
time_final = 0
# Initial time
time_initial = datetime.datetime.now()
file = str(sys.argv[1])
# Variable initialization with information from command line
f = open(file, 'r')
lines = f.readlines()
f.close()
for line in lines:
l = line.strip()
if l.startswith("initialization_file"):
init = l.split("initialization_file = ", 1)[1]
if l.startswith("input_dir"):
input_dir = l.split("input_dir = ", 1)[1]
if l.startswith("programs_dir"):
programs_dir = l.split("programs_dir = ", 1)[1]
if l.startswith("substrate_nat"):
substrate_nat = l.split("substrate_nat = ", 1)[1]
if l.startswith("cluster_ntyp"):
cluster_ntyp = l.split("cluster_ntyp = ", 1)[1]
if l.startswith("temperature_K"):
temperature_K = float(l.split("temperature_K = ", 1)[1])
if l.startswith("step_width"):
step_width = l.split("step_width = ", 1)[1]
if l.startswith("iterations"):
iterations = int(l.split("iterations = ", 1)[1])
if l.startswith("x_range"):
x_range = l.split("x_range = ", 1)[1]
if l.startswith("y_range"):
y_range = l.split("y_range = ", 1)[1]
if l.startswith("z_range"):
z_range = l.split("z_range = ", 1)[1]
if l.startswith("z_vacuum"):
z_vacuum = l.split("z_vacuum = ", 1)[1]
sys.path.append(programs_dir)
from input_parser import*
# Cluster total number of atoms
cluster_nat = str(get_nAtoms_str(cluster_ntyp))
# Temperature
kBT = float(0.00008617 * temperature_K)
# Output folder name
output_dir_name = get_output_name(cluster_ntyp)
# Output folder is created
call(['cp', '-r', input_dir, './' + output_dir_name])
# The input file gets edited to remain consistent with the number of atoms
call(['python3.4', programs_dir + '/replace_line.py', output_dir_name + '/input.in', 'nat=n,', 'nat=' + str(int(cluster_nat) + int(substrate_nat)) + ','])
# The cluster gets randomly generated only if there is no initialization file
if init == "False" or init == "false":
call(['python3.4', programs_dir + '/RandomGenerator.py', output_dir_name + '/input.in', cluster_ntyp, x_range, y_range, z_range, z_vacuum])
else:
call(['python3.4', programs_dir + '/initializer.py', init, output_dir_name + '/input.in'])
print("=========================================================================================================")
print("Basin Hopping Monte Carlo global optimization program coupled to (DFT-Quantum Espresso) ")
print("for the analysis of supported clusters, and (eventually) surface/bulk alloys ")
print("Instituto de Fisica - UNAM (2016) ")
print("=========================================================================================================")
print(" ")
print(" --> Starting DFT relaxation of initial/random configuration ")
print("Iteration 1 of " + output_dir_name)
print(" ")
# The "run" script gets executed, with particular emphasis on its completion before continuing with
# the subsequent commands in the script
subproc = Popen([output_dir_name + '/run.sh'], stdout = PIPE, stderr = PIPE)
(out, err) = subproc.communicate()
while call('grep -q \"Begin final coordinates\" ' + output_dir_name + '/output.out', shell = True) == 1:
print(" --> SCF failed. Starting again from randomly generated structure! " + output_dir_name)
# The cluster gets randomly generated again; ie, the original input gets tossed out
call(['python3.4', programs_dir + '/RandomGenerator.py', output_dir_name + '/input.in', cluster_ntyp, x_range, y_range, z_range, z_vacuum])
# The failing output file gets deleted, for the next one to take its place
call(['rm', output_dir_name + '/output.out'])
subproc = Popen([output_dir_name + '/run.sh'], stdout = PIPE, stderr = PIPE)
(out, err) = subproc.communicate()
# Naming of corresponding input and output
call(['cp', output_dir_name + '/input.in', output_dir_name + '/input1.in'])
call(['cp', output_dir_name + '/output.out', output_dir_name + '/output1.out'])
#Naming of final coordinates and energy file
call('grep \"! \" ' + output_dir_name + '/output1.out | tail -1 > ' + output_dir_name + '/coord1.xyz', shell = True)
# -A to take the necessary subsequent lines, -m 1 to take only the first ocurrence of the argument
call('grep -A ' + str(int(substrate_nat) + 2 + int(cluster_nat)) + ' \"Begin final coordinates\" ' + output_dir_name + '/output1.out | head -' + str(int(substrate_nat) + 3 + int(cluster_nat)) + ' >> ' + output_dir_name + '/coord1.xyz', shell = True)
# Deletion of unnecessary files
call(['rm', '-r', output_dir_name + '/pwscf.save', output_dir_name + '/output.out'])
for fl in glob.glob(output_dir_name + '/pwscf.*'):
os.remove(fl)
print(" ")
print(" --> Relaxation of initial configuration: DONE! ")
print(" ")
print("=========================================================================================================")
# Steps for the Basin hopping
print("BH-DFT routine starts here! ")
print("Note: ")
print("For monometallic clusters: only random xyz moves will be applied ")
print("For bimetallic clusters : 1 atomic swap will be performed after 10 moves ")
print("=========================================================================================================")
print(" ")
i = 1
# In case swap energy fails
swap_fail_flag = False
while i < iterations:
# Every 10th iteration a swap routine gets executed
if i%10 == 0 and swap_fail_flag == False:
call(['python3.4', programs_dir + '/swap.py', output_dir_name + '/coord' + str(i) + '.xyz', output_dir_name + '/input.in', cluster_ntyp])
else:
swap_fail_flag = False
call(['python3.4', programs_dir + '/move.py', output_dir_name + '/coord' + str(i) + '.xyz', output_dir_name + '/input.in', step_width, cluster_ntyp])
print("Iteration " + str(i + 1)+ " of structure: " + output_dir_name)
# The "run" script gets executed, with particular emphasis on its completion before continuing with
# the subsequent commands in the script
subproc = Popen([output_dir_name + '/run.sh'], stdout = PIPE, stderr = PIPE)
(out, err) = subproc.communicate()
while call('grep -q \"Begin final coordinates\" ' + output_dir_name + '/output.out', shell = True) == 1:
print(" --> SCF failed. Starting again from previous configuration " + str(i + 1) + " of " + output_dir_name)
# The Pt coordinates get moved again; ie, the original input gets tossed out
call(['python3.4', programs_dir + '/move.py', output_dir_name + '/coord' + str(i) + '.xyz', output_dir_name + '/input.in', step_width, cluster_ntyp])
# The failing output file gets deleted, for the next one to take its place
call(['rm', output_dir_name + '/output.out'])
subproc = Popen([output_dir_name + '/run.sh'], stdout = PIPE, stderr = PIPE)
(out, err) = subproc.communicate()
# Naming of corresponding input and output
call(['cp', output_dir_name + '/input.in', output_dir_name + '/input' + str(i + 1) + '.in'])
call(['cp', output_dir_name + '/output.out', output_dir_name + '/output' + str(i + 1) + '.out'])
# Creation of final energy and coordinates file
call('grep \"! \" ' + output_dir_name + '/output' + str(i + 1) + '.out | tail -1 > ' + output_dir_name + '/coord' + str(i + 1) + '.xyz', shell = True)
call('grep -A ' + str(int(substrate_nat) + 2 + int(cluster_nat)) + ' \"Begin final coordinates\" ' + output_dir_name + '/output' + str(i + 1) + '.out | head -' + str(int(substrate_nat) + 3 + int(cluster_nat)) + ' >> ' + output_dir_name + '/coord' + str(i + 1) + '.xyz', shell = True)
# Deletion of unnecessary files
call(['rm', '-r', output_dir_name + '/pwscf.save', output_dir_name + '/output.out'])
for fl in glob.glob(output_dir_name + '/pwscf.*'):
os.remove(fl)
# Metropolis condition for accepting
f = open(output_dir_name + '/coord' + str(i) + '.xyz', 'r')
E0 = float(re.findall("\d+\.\d+", f.readline())[0]) * -13.6
f.close()
f = open(output_dir_name + '/coord' + str(i + 1) + '.xyz', 'r')
En = float(re.findall("\d+\.\d+", f.readline())[0]) * -13.6
f.close()
if math.exp((E0 - En)/kBT) > random.uniform(0,1):
print(" --> Basin Hopping MC criteria: Energy accepted! ")
print(" --> Finished iteration # " + str(i + 1))
i = i + 1
continue
else:
print(" --> Basin Hopping MC criteria: Energy rejected!")
call(['rm', '-r', output_dir_name + '/input' + str(i + 1) + '.in', output_dir_name + '/output' + str(i + 1) + '.out', output_dir_name + '/coord' + str(i + 1) + '.xyz'])
print("Iteration " + str(i + 1) + " failed!. Starting again from previous configuration + random moves ! ")
# Every 10th iteration swap gets involved
if i%10 == 0:
print(" --> Swap failed to converge.")
swap_fail_flag = True
# Print elapsed time
time_final = datetime.datetime.now()
print("\nTotal execution time: " + str(time_final - time_initial))