add magnetic option to bulk calculation

Author: sevyharris

adlib/bulk/calc.py

@@ -99,7 +99,13 @@ def make_scf_run_file_array(dest_dir, N_runs, job_name='bulk_energy', nproc=16):
         f.write(f'python calc.py\n')
 
 
-def make_scf_calc_file(calc_dir, lattice_constant, metal='Cu', crystal_structure='fcc', ecutwfc=1000, kpt=9, smear=0.1, nproc=16):
+def make_scf_calc_file(calc_dir, lattice_constant, metal='Cu', crystal_structure='fcc', ecutwfc=1000, kpt=9, smear=0.1, nproc=16, magnetism=None):
+
+    magnetism_line = ""
+    if str(magnetism).lower() == 'antiferromagnetic':
+        magnetism_line = "bulk.set_initial_magnetic_moments([1.0, -1.0])\nassert len(bulk) == 2"
+    elif str(magnetism).lower() == 'ferromagnetic':
+        magnetism_line = "bulk.set_initial_magnetic_moments([1.0, 1.0])\nassert len(bulk) == 2"
 
     python_file_lines = [
         "import os",
@@ -134,6 +140,7 @@ def make_scf_calc_file(calc_dir, lattice_constant, metal='Cu', crystal_structure
         "",
         f"bulk = ase.build.bulk('{metal}', crystalstructure='{crystal_structure}', a={lattice_constant}, cubic=True)",
         "",
+        magnetism_line,
         "pw_executable = os.environ['PW_EXECUTABLE']",
         "",
         "pseudopotentials = {",

adlib/bulk/convergence.py

@@ -46,7 +46,7 @@
 import adlib.bulk.calc
 
 
-def setup_kpts_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc'):
+def setup_kpts_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc', magnetism=None):
     """
     script to set up N jobs to check kpts convergence
     """
@@ -61,12 +61,12 @@ def setup_kpts_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_s
     for i, k in enumerate(kpts):
         for j, lattice_constant in enumerate(lattice_constants):
             calc_dir = os.path.join(kpts_dir, str(k), f'run_{j:04}')
-            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, ecutwfc=500, kpt=k, smear=0.1, nproc=16)
+            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, magnetism=magnetism, ecutwfc=500, kpt=k, smear=0.1, nproc=16)
 
         adlib.bulk.calc.make_scf_run_file_array(os.path.join(kpts_dir, str(k)), j, job_name='kpts_bulk_converge')
 
 
-def setup_ecutwfc_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc'):
+def setup_ecutwfc_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc', magnetism=None):
     """
     script to set up N jobs to check ecutwfc convergence
     """
@@ -80,12 +80,12 @@ def setup_ecutwfc_converge(bulk_dir, lattice_constant_center, metal='Cu', crysta
     for i, ecut in enumerate(ecuts):
         for j, lattice_constant in enumerate(lattice_constants):
             calc_dir = os.path.join(ecutwfc_dir, str(ecut), f'run_{j:04}')
-            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, ecutwfc=ecut, kpt=7, smear=0.1, nproc=16)
+            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, magnetism=magnetism, ecutwfc=ecut, kpt=7, smear=0.1, nproc=16)
 
         adlib.bulk.calc.make_scf_run_file_array(os.path.join(ecutwfc_dir, str(ecut)), j, job_name='ecut_bulk_converge')
 
 
-def setup_smear_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc'):
+def setup_smear_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_structure='fcc', magnetism=None):
     """
     script to set up N jobs to check MV smearing convergence
     """
@@ -98,7 +98,7 @@ def setup_smear_converge(bulk_dir, lattice_constant_center, metal='Cu', crystal_
     for i, smear in enumerate(smears):
         for j, lattice_constant in enumerate(lattice_constants):
             calc_dir = os.path.join(smear_dir, str(smear), f'run_{j:04}')
-            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, ecutwfc=1000, kpt=7, smear=smear, nproc=16)
+            adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, magnetism=magnetism, ecutwfc=1000, kpt=7, smear=smear, nproc=16)
 
         adlib.bulk.calc.make_scf_run_file_array(os.path.join(smear_dir, str(smear)), j, job_name='smear_bulk_converge')
 

adlib/bulk/eos.py

@@ -53,7 +53,7 @@
 import adlib.bulk.calc
 
 
-def setup_eos(bulk_dir, lattice_constant_guess, metal='Cu', N=21, half_range=0.05):
+def setup_eos(bulk_dir, lattice_constant_guess, metal='Cu', N=21, half_range=0.05, magnetism=None):
     """
     script to set up N jobs
 
@@ -71,7 +71,7 @@ def setup_eos(bulk_dir, lattice_constant_guess, metal='Cu', N=21, half_range=0.0
 
     for i, lattice_constant in enumerate(lattice_constants):
         calc_dir = os.path.join(eos_dir, f'run_{i:04}')
-        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, ecutwfc=1000, kpt=9, smear=0.1, nproc=16)
+        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, magnetism=magnetism, ecutwfc=1000, kpt=9, smear=0.1, nproc=16)
 
     adlib.bulk.calc.make_scf_run_file_array(eos_dir, i, job_name='bulk_eos')
 
@@ -90,7 +90,7 @@ def get_lattice_constant_from_atoms(atoms):
     return lattice_constant
 
 
-def setup_eos_coarse(bulk_dir, lattice_constant_guess, metal='Cu'):
+def setup_eos_coarse(bulk_dir, lattice_constant_guess, metal='Cu', magnetism=None):
     """
     script to set up coarse calculation of E vs. lattice constant
     """
@@ -102,12 +102,12 @@ def setup_eos_coarse(bulk_dir, lattice_constant_guess, metal='Cu'):
 
     for i, lattice_constant in enumerate(lattice_constants):
         calc_dir = os.path.join(eos_dir, f'run_{i:04}')
-        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal)
+        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, magnetism=magnetism)
 
     adlib.bulk.calc.make_scf_run_file_array(eos_dir, i, job_name='bulk_eos_coarse')
 
 
-def setup_eos_fine(bulk_dir, lattice_constant_guess, metal='Cu', crystal_structure='fcc'):
+def setup_eos_fine(bulk_dir, lattice_constant_guess, metal='Cu', crystal_structure='fcc', magnetism=None):
     """
     script to set up fine calculation of E vs. lattice constant
     """
@@ -119,7 +119,7 @@ def setup_eos_fine(bulk_dir, lattice_constant_guess, metal='Cu', crystal_structu
 
     for i, lattice_constant in enumerate(lattice_constants):
         calc_dir = os.path.join(eos_dir, f'run_{i:04}')
-        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure)
+        adlib.bulk.calc.make_scf_calc_file(calc_dir, lattice_constant, metal=metal, crystal_structure=crystal_structure, magnetism=magnetism)
 
     adlib.bulk.calc.make_scf_run_file_array(eos_dir, i, job_name='bulk_eos_fine')
 
@@ -182,8 +182,11 @@ def plot_energy_vs_lattice(calc_dir, dest_dir=None):
     for pwo_file in pwo_files:
         with open(pwo_file, 'r') as f:
             traj = list(read_espresso_out(f, index=slice(None)))
-            atoms = traj[-1]
-            energies.append(atoms.get_potential_energy())
+            try:
+                atoms = traj[-1]
+                energies.append(atoms.get_potential_energy())
+            except IndexError:
+                energies.append(np.nan)
 
             lattice_constant = get_lattice_constant_from_atoms(atoms)
             lattice_constants.append(lattice_constant)

adlib/bulk/vcrelax.py

@@ -76,7 +76,7 @@ def make_run_vc_relax_script(calc_dir, nproc=16):
         f.write(f'python relax_bulk.py\n')
 
 
-def make_vc_relax_script(calc_dir, metal, lattice_constant, crystal_structure='fcc', nproc=16):
+def make_vc_relax_script(calc_dir, metal, lattice_constant, crystal_structure='fcc', magnetism=None, nproc=16):
     """
     Make a python script that uses ase to run quantum espresso
     """
@@ -85,12 +85,18 @@ def make_vc_relax_script(calc_dir, metal, lattice_constant, crystal_structure='f
     kpts = 5
     ecutwfc = 500
 
+    magnetism_line = ""
+    if str(magnetism).lower() == 'antiferromagnetic':
+        magnetism_line = "bulk.set_initial_magnetic_moments([1.0, -1.0])\nassert len(bulk) == 2"
+    elif str(magnetism).lower() == 'ferromagnetic':
+        magnetism_line = "bulk.set_initial_magnetic_moments([1.0, 1.0])\nassert len(bulk) == 2"
+
     python_file_lines = [
         "import os",
         "import sys",
         "import time",
         "from ase.calculators.espresso import Espresso, EspressoProfile",
-        "from ase.build import bulk",
+        "import ase.build",
         "",
         "",
         "T = time.localtime()",
@@ -123,7 +129,8 @@ def make_vc_relax_script(calc_dir, metal, lattice_constant, crystal_structure='f
         "}",
         "",
         "",
-        f"bulk_metal = bulk('{metal}', crystalstructure='{crystal_structure}', a={lattice_constant}, cubic=True)",
+        f"bulk_metal = ase.build.bulk('{metal}', crystalstructure='{crystal_structure}', a={lattice_constant}, cubic=True)",
+        magnetism_line,
         "",
         "pw_executable = os.environ['PW_EXECUTABLE']",
         "",
@@ -183,7 +190,7 @@ def run_vc_relax(bulk_dir):
     calc_dir = os.path.join(bulk_dir, 'vc_relax')
     cur_dir = os.getcwd()
     os.chdir(calc_dir)
-    if environment == 'DISCOVERY':
+    if environment in ['DISCOVERY', 'EXPLORER']:
         vc_relax_job = job_manager.SlurmJob()
         cmd = "sbatch run.sh"
     elif environment == 'SINGLE_NODE':