|
| 1 | +""" |
| 2 | +Module for setting up and running diffusion barrier calculations in ase/Quantum Espresso |
| 3 | +
|
| 4 | +diffusion |
| 5 | +____Pt111_top_top |
| 6 | +________O |
| 7 | +____________site1.pwo |
| 8 | +____________espresso.pwo |
| 9 | +____________calc.py |
| 10 | +________H |
| 11 | +____________site1.pwo |
| 12 | +____________espresso.pwo |
| 13 | +____________calc.py |
| 14 | +____Pt111_top_fcc |
| 15 | +________O |
| 16 | +____________site1.pwo |
| 17 | +____________espresso.pwo |
| 18 | +____________calc.py |
| 19 | +________H |
| 20 | +____________site1.pwo |
| 21 | +____________espresso.pwo |
| 22 | +____________calc.py |
| 23 | +____Au111_top_brg |
| 24 | +________O |
| 25 | +____________site1.pwo |
| 26 | +____________espresso.pwo |
| 27 | +____________calc.py |
| 28 | +________H |
| 29 | +____________site1.pwo |
| 30 | +____________espresso.pwo |
| 31 | +____________calc.py |
| 32 | +""" |
| 33 | + |
| 34 | +import os |
| 35 | +import adlib.env |
| 36 | + |
| 37 | + |
| 38 | +environment = adlib.env.load_environment() |
| 39 | + |
| 40 | + |
| 41 | +def make_run_relax_script(calc_dir, nproc=48, job_name='relax_system'): |
| 42 | + bash_filename = os.path.join(calc_dir, 'run.sh') |
| 43 | + # write the array job file |
| 44 | + with open(bash_filename, 'w') as f: |
| 45 | + f.write('#!/bin/bash\n\n') |
| 46 | + if environment == 'DISCOVERY': |
| 47 | + f.write('#SBATCH --time=24:00:00\n') |
| 48 | + f.write(f'#SBATCH --job-name={job_name}\n') |
| 49 | + f.write('#SBATCH --mem=40Gb\n') |
| 50 | + f.write('#SBATCH --cpus-per-task=1\n') |
| 51 | + f.write(f'#SBATCH --ntasks={nproc}\n') |
| 52 | + f.write('#SBATCH --partition=short\n') |
| 53 | + f.write('#SBATCH --constraint=cascadelake\n\n') |
| 54 | + f.write('module load gcc/10.1.0\n') |
| 55 | + f.write('module load openmpi/4.0.5-skylake-gcc10.1\n') |
| 56 | + f.write('module load scalapack/2.1.0-skylake\n\n') |
| 57 | + # f.write(f'cd {calc_dir}\n') |
| 58 | + f.write(f'python relax_system.py\n') |
| 59 | + |
| 60 | + |
| 61 | +def make_relax_script(calc_dir, ecutwfc=60, kpt=5, smear=0.1, nproc=48): |
| 62 | + """Function to make a python script to relax the slab-adsorption system |
| 63 | + """ |
| 64 | + fmax = 0.01 |
| 65 | + vacuum = 7.5 |
| 66 | + |
| 67 | + python_file_lines = [ |
| 68 | + "import os", |
| 69 | + "import sys", |
| 70 | + "import numpy as np", |
| 71 | + "from ase.build import bulk, fcc111, add_adsorbate", |
| 72 | + "from ase.io import read, write", |
| 73 | + "from ase import Atoms", |
| 74 | + "from ase.optimize import BFGS", |
| 75 | + "from ase.calculators.espresso import Espresso, EspressoProfile", |
| 76 | + "from ase.constraints import FixAtoms", |
| 77 | + "from ase.io.espresso import read_espresso_out", |
| 78 | + "from ase.io.trajectory import Trajectory", |
| 79 | + "from ase.io.ulm import InvalidULMFileError", |
| 80 | + "import time", |
| 81 | + "", |
| 82 | + "", |
| 83 | + "def place_adsorbate_top(metal_slab, adsorbate, height=1.0, ads_index=0):", |
| 84 | + " # remove the adsorbate cell", |
| 85 | + " adsorbate.cell = [0, 0, 0]", |
| 86 | + " top_layer_z = np.max([pos[2] for pos in metal_slab.get_positions()])", |
| 87 | + " for i, pos in enumerate(metal_slab.get_positions()):", |
| 88 | + " if pos[2] == top_layer_z:", |
| 89 | + " # place the atom height above here", |
| 90 | + " ads_origin = pos + [0, 0, height]", |
| 91 | + " pos_difference = ads_origin - adsorbate.get_positions()[ads_index]", |
| 92 | + " adsorbate.translate(pos_difference)", |
| 93 | + " metal_slab += adsorbate", |
| 94 | + " print(f'Using {adsorbate[ads_index]} for binding atom')", |
| 95 | + " return", |
| 96 | + " print('Failed to place adsorbate')", |
| 97 | + " exit(-1)", |
| 98 | + "", |
| 99 | + "", |
| 100 | + "# TODO accept height guess", |
| 101 | + "# if len(sys.argv) < 2:", |
| 102 | + "# raise IndexError('Must specify starting height for fine system run')", |
| 103 | + "# height = float(sys.argv[1])", |
| 104 | + "height = 4.0", |
| 105 | + f"fmax = {fmax}", |
| 106 | + "T = time.localtime()", |
| 107 | + "start = time.time()", |
| 108 | + "", |
| 109 | + "logfile = 'ase.log'", |
| 110 | + "with open(logfile, 'a') as f:", |
| 111 | + " f.write(f'Start: {T[3]:02}:{T[4]:02}:{T[5]:02}\\n')", |
| 112 | + "", |
| 113 | + "slab_file = 'slab.pwo'", |
| 114 | + "with open(slab_file, 'r') as f:", |
| 115 | + " traj = list(read_espresso_out(f, index=slice(None)))", |
| 116 | + "metal_slab = traj[-1]", |
| 117 | + "", |
| 118 | + "adsorbate_file = 'adsorbate.pwo'", |
| 119 | + "with open(adsorbate_file, 'r') as f:", |
| 120 | + " traj = list(read_espresso_out(f, index=slice(None)))", |
| 121 | + "adsorbate = traj[-1]", |
| 122 | + "", |
| 123 | + "# Fix the bottom two layers", |
| 124 | + "bottom_layer = []", |
| 125 | + "second_layer = []", |
| 126 | + "fixed_indices = []", |
| 127 | + "# Round to the nearest 0.1 Angstrom in determining layers", |
| 128 | + "z_values = list(set([np.round(pos[2], 1) for pos in metal_slab.get_positions()]))", |
| 129 | + "z_values.sort()", |
| 130 | + "", |
| 131 | + "", |
| 132 | + "for i, pos in enumerate(metal_slab.get_positions()):", |
| 133 | + " if np.round(pos[2], 1) == z_values[0]:", |
| 134 | + " bottom_layer.append(metal_slab[i].index)", |
| 135 | + " if np.round(pos[2], 1) == z_values[1]:", |
| 136 | + " second_layer.append(metal_slab[i].index)", |
| 137 | + "fixed_indicies = bottom_layer + second_layer", |
| 138 | + "fix_bottom_layers = FixAtoms(indices=fixed_indicies)", |
| 139 | + "metal_slab.set_constraint(fix_bottom_layers)", |
| 140 | + "", |
| 141 | + "", |
| 142 | + "# place the adsorbate", |
| 143 | + "element_priority = ['C', 'O', 'H'] # where does N fit?", |
| 144 | + "bond_atom_index = -1", |
| 145 | + "for element in element_priority:", |
| 146 | + " for atom in adsorbate:", |
| 147 | + " if atom.symbol == element:", |
| 148 | + " bond_atom_index = atom.index", |
| 149 | + " break", |
| 150 | + " if bond_atom_index > -1:", |
| 151 | + " break", |
| 152 | + "print(f'bond atom is {adsorbate[bond_atom_index]}')", |
| 153 | + "", |
| 154 | + "# 'ontop', 'bridge', 'fcc', 'hcp'", |
| 155 | + "# https://wiki.fysik.dtu.dk/ase/ase/build/surface.html#ase.build.add_adsorbate", |
| 156 | + "# add_adsorbate(metal_slab, adsorbate, height=height, position='ontop', mol_index=bond_atom_index)", |
| 157 | + "# can't use the add_adsorbate function if you load the atoms from a pwo file...", |
| 158 | + "place_adsorbate_top(metal_slab, adsorbate, height=height, ads_index=bond_atom_index)", |
| 159 | + "", |
| 160 | + "", |
| 161 | + "# Restart if available", |
| 162 | + "traj_file = 'system.traj'", |
| 163 | + "restart = False", |
| 164 | + "if os.path.exists(traj_file):", |
| 165 | + " try:", |
| 166 | + " traj = Trajectory(traj_file)", |
| 167 | + " if len(traj) > 0:", |
| 168 | + " metal_slab = traj[-1]", |
| 169 | + " restart = True", |
| 170 | + " with open(logfile, 'a') as f:", |
| 171 | + " f.write(f'Restarting from trajectory file {traj_file}\\n')", |
| 172 | + " except InvalidULMFileError:", |
| 173 | + " # traj file is empty. delete it.", |
| 174 | + " os.remove(traj_file)", |
| 175 | + "", |
| 176 | + "if not restart:", |
| 177 | + " write(f'initial_system.xyz', metal_slab)", |
| 178 | + "", |
| 179 | + "espresso_settings = {", |
| 180 | + " 'control': {", |
| 181 | + " 'verbosity': 'high',", |
| 182 | + " 'calculation': 'scf',", |
| 183 | + " 'disk_io': 'none',", |
| 184 | + " },", |
| 185 | + " 'system': {", |
| 186 | + " 'input_dft': 'BEEF-VDW',", |
| 187 | + " 'occupations': 'smearing',", |
| 188 | + " 'smearing': 'mv',", |
| 189 | + f" 'degauss': {smear},", |
| 190 | + f" 'ecutwfc': {ecutwfc},", |
| 191 | + " },", |
| 192 | + "}", |
| 193 | + "", |
| 194 | + "pseudopotentials = {", |
| 195 | + " 'C': 'C_ONCV_PBE-1.2.upf',", |
| 196 | + " 'Cu': 'Cu_ONCV_PBE-1.2.upf',", |
| 197 | + " 'O': 'O_ONCV_PBE-1.2.upf',", |
| 198 | + " 'N': 'N_ONCV_PBE-1.2.upf',", |
| 199 | + " 'H': 'H_ONCV_PBE-1.2.upf',", |
| 200 | + " 'Pt': 'Pt_ONCV_PBE-1.2.upf',", |
| 201 | + " 'Pd': 'Pd_ONCV_PBE-1.2.upf',", |
| 202 | + " 'Au': 'Au_ONCV_PBE-1.2.upf',", |
| 203 | + " 'Ag': 'Ag_ONCV_PBE-1.2.upf',", |
| 204 | + " 'Al': 'Al_ONCV_PBE-1.2.upf',", |
| 205 | + " 'Ni': 'Ni_ONCV_PBE-1.2.upf',", |
| 206 | + " 'Fe': 'Fe_ONCV_PBE-1.2.upf',", |
| 207 | + " 'Cr': 'Cr_ONCV_PBE-1.2.upf',", |
| 208 | + "}", |
| 209 | + "", |
| 210 | + "pw_executable = os.environ['PW_EXECUTABLE']", |
| 211 | + f"command = f'mpirun -np {nproc} " + "{pw_executable} -in PREFIX.pwi > PREFIX.pwo'", |
| 212 | + "profile = EspressoProfile(argv=command.split())", |
| 213 | + "calc = Espresso(", |
| 214 | + " # command=command,", |
| 215 | + " profile=profile,", |
| 216 | + " pseudopotentials=pseudopotentials,", |
| 217 | + " tstress=True,", |
| 218 | + " tprnfor=True,", |
| 219 | + f" kpts=({kpt}, {kpt}, 1),", |
| 220 | + " pseudo_dir=os.environ['PSEUDO_DIR'],", |
| 221 | + " input_data=espresso_settings,", |
| 222 | + ")", |
| 223 | + "", |
| 224 | + "metal_slab.calc = calc", |
| 225 | + "opt = BFGS(metal_slab, logfile=logfile, trajectory='system.traj')", |
| 226 | + "opt.run(fmax=fmax)", |
| 227 | + "", |
| 228 | + "", |
| 229 | + "# Read the energy back in", |
| 230 | + "energy = 0", |
| 231 | + "with open('espresso.pwo', 'r') as f:", |
| 232 | + " traj = list(read_espresso_out(f, index=slice(None)))", |
| 233 | + " energy = traj[-1].get_potential_energy()", |
| 234 | + "", |
| 235 | + "", |
| 236 | + "T = time.localtime()", |
| 237 | + "end = time.time()", |
| 238 | + "duration = end - start", |
| 239 | + "", |
| 240 | + "with open(logfile, 'a') as f:", |
| 241 | + " f.write(f'Energy: {energy}\\n')", |
| 242 | + " f.write(f'End: {T[3]:02}:{T[4]:02}:{T[5]:02}\\n')", |
| 243 | + " f.write(f'Completed in {duration} seconds\\n')", |
| 244 | + ] |
| 245 | + |
| 246 | + os.makedirs(calc_dir, exist_ok=True) |
| 247 | + calc_filename = os.path.join(calc_dir, f'relax_system.py') |
| 248 | + with open(calc_filename, 'w') as f: |
| 249 | + f.writelines([line + '\n' for line in python_file_lines]) |
| 250 | + |
| 251 | + |
| 252 | +def run_relax_system(system_dir): |
| 253 | + import job_manager |
| 254 | + cur_dir = os.getcwd() |
| 255 | + os.chdir(system_dir) |
| 256 | + if environment == 'DISCOVERY': |
| 257 | + relax_system_job = job_manager.SlurmJob() |
| 258 | + cmd = "sbatch run.sh" |
| 259 | + elif environment == 'SINGLE_NODE': |
| 260 | + relax_system_job = job_manager.DefaultJob() |
| 261 | + cmd = "/bin/bash run.sh" |
| 262 | + elif environment == 'THETA': |
| 263 | + raise NotImplementedError |
| 264 | + else: |
| 265 | + raise NotImplementedError |
| 266 | + |
| 267 | + relax_system_job.submit(cmd) |
| 268 | + os.chdir(cur_dir) |
0 commit comments