2. Equilibrium volume of FCC Si

Here we will follow the second part of the FCC Si example in the VASP tutorial, i.e. perform the calculations of the total energies at different volumes. Please make sure you have completed the previous tutorial before continuing.

  1. When the standard tutorial is completed we will now do the same in AiiDA-VASP using different strategies to you get a feel for how you can structurize a simple workflow like this. Let us now fetch the AiiDA-VASP run file for this example:

    %/$ wget https://github.com/aiida-vasp/aiida-vasp/raw/develop/tutorials/run_fcc_si_multiple_volumes.py
    
  2. Inspect the file, which has the following changes compared to the file used in the previous tutorial

    --- /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/latest/tutorials/run_fcc_si_one_volume.py
    +++ /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/latest/tutorials/run_fcc_si_multiple_volumes.py
    @@ -1,8 +1,8 @@
     """
    -Call script to calculate the total energies for one volume of standard silicon.
    +Call script to calculate the total energies for different volumes of the silicon structure.
     
    -This particular call script set up a standard calculation that execute a calculation for
    -the fcc silicon structure.
    +This particular call script set up a standard calculation for each structure and submits
    +each of them.
     """
     # pylint: disable=too-many-arguments
     import numpy as np
    @@ -14,12 +14,12 @@
     load_profile()
     
     
    -def get_structure():
    +def get_structure(alat):
         """
         Set up Si primitive cell
     
         fcc Si:
    -       3.9
    +       alat
            0.5000000000000000    0.5000000000000000    0.0000000000000000
            0.0000000000000000    0.5000000000000000    0.5000000000000000
            0.5000000000000000    0.0000000000000000    0.5000000000000000
    @@ -31,11 +31,9 @@
         """
     
         structure_data = DataFactory('structure')
    -    alat = 3.9
         lattice = np.array([[.5, .5, 0], [0, .5, .5], [.5, 0, .5]]) * alat
         structure = structure_data(cell=lattice)
    -    positions = [[0.0, 0.0, 0.0]]
    -    for pos_direct in positions:
    +    for pos_direct in ([[0.0, 0.0, 0.0]]):
             pos_cartesian = np.dot(pos_direct, lattice)
             structure.append_atom(position=pos_cartesian, symbols='Si')
         return structure
    @@ -117,8 +115,10 @@
         OPTIONS.max_wallclock_seconds = 3600
         OPTIONS.max_memory_kb = 1024000
     
    -    # POSCAR equivalent
    -    # Set the silicon structure
    -    STRUCTURE = get_structure()
    +    # Iterate over each lattice constant and pass it explicitly
    +    for lattice_constant in [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3]:
    +        # POSCAR equivalent
    +        # Set the silicon structure
    +        STRUCTURE = get_structure(lattice_constant)
     
    -    main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)
    +        main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)
    

    As you can see, we did only minimal changes, all which should be self explanatory.

  3. Change the options and code_string as you did in the previous tutorial.

    Warning

    Make sure you have activated your AiiDA virtual environment and that the AiiDA daemon is running before continuing.

  4. Save and execute the resulting run script by issuing:

    %/$ python run_fcc_si_multiple_volumes.py
    
  5. Let us check the progress:

    %/$ verdi process list -a
    PK  Created    Process label                 Process State      Process status
    ------  ---------  ----------------------------  -----------------  -----------------------------------------------------------
    101383  22s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101393
    101392  22s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101403
    101393  21s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101413
    101402  21s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101423
    101403  21s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101424
    101412  21s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101438
    101413  21s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101422  21s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101444
    101423  21s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101454
    101424  20s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101433  20s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101465
    101443  20s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101467
    101438  20s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101464
    101444  20s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101466
    101453  20s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101468
    101454  19s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101463  19s ago    VerifyWorkChain               ⏵ Waiting          Waiting for child processes: 101469
    101465  19s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101471
    101464  19s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101466  18s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101467  18s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101472
    101468  18s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101470
    101469  17s ago    VaspWorkChain                 ⏵ Waiting          Waiting for child processes: 101473
    101470  17s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101471  17s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: submit
    101472  16s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: upload
    101473  16s ago    VaspCalculation               ⏵ Waiting          Waiting for transport task: submit
    

    All calculations are submitted to the cluster in parallel, which is pretty convenient and certainly an improvement compared to the regular tutorial, where you do the calculations sequentially. Not a problem for a few calculations, but maybe not ideal if you want to execute a few thousand calculations.

  6. After a while, we execut verdi process list -a and get:

    %/$ verdi process list -a
    PK  Created    Process label                 Process State      Process status
    ------  ---------  ----------------------------  -----------------  -----------------------------------------------------------
    101383  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101392  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101393  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101402  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101403  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101412  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101413  2m ago     VaspCalculation               ⏹ Finished [0]
    101422  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101423  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101424  2m ago     VaspCalculation               ⏹ Finished [0]
    101433  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101443  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101438  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101444  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101453  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101454  2m ago     VaspCalculation               ⏹ Finished [0]
    101463  2m ago     VerifyWorkChain               ⏹ Finished [0]
    101465  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101464  2m ago     VaspCalculation               ⏹ Finished [0]
    101466  2m ago     VaspCalculation               ⏹ Finished [0]
    101467  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101468  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101469  2m ago     VaspWorkChain                 ⏹ Finished [0]
    101470  2m ago     VaspCalculation               ⏹ Finished [0]
    101471  2m ago     VaspCalculation               ⏹ Finished [0]
    101472  2m ago     VaspCalculation               ⏹ Finished [0]
    101473  2m ago     VaspCalculation               ⏹ Finished [0]
    

    All processes are in a finished state and we can extract the total energies for each step. However, it should be obvious that extracting the total energies from misc from each step manually seems to be a waste of time. One could query (AiiDA has shortcuts for this), but one would still need to look up the values in some way. Maybe it is easier to be able to access them directly when all the calculations are complete? Let us do another modification to the run script above, namely:

    --- /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/latest/tutorials/run_fcc_si_multiple_volumes.py
    +++ /home/docs/checkouts/readthedocs.org/user_builds/aiida-vasp/checkouts/latest/tutorials/run_fcc_si_multiple_volumes_eos.py
    @@ -2,14 +2,15 @@
     Call script to calculate the total energies for different volumes of the silicon structure.
     
     This particular call script set up a standard calculation for each structure and submits
    -each of them.
    +each of them. We, in addition collect the total energies and eject them to a file called
    +eos.
     """
     # pylint: disable=too-many-arguments
     import numpy as np
     from aiida.common.extendeddicts import AttributeDict
     from aiida.orm import Code, Bool, Str
     from aiida.plugins import DataFactory, WorkflowFactory
    -from aiida.engine import submit
    +from aiida.engine import run
     from aiida import load_profile
     load_profile()
     
    @@ -54,31 +55,32 @@
         settings = AttributeDict()
         inputs = AttributeDict()
     
    -    # organize settings
    +    # Organize settings
         settings.parser_settings = {'output_params': ['total_energies', 'maximum_force']}
     
    -    # set inputs for the following WorkChain execution
    -    # set code
    +    # Set inputs for the following WorkChain execution
    +    # Set code
         inputs.code = Code.get_from_string(code_string)
    -    # set structure
    +    # Set structure
         inputs.structure = structure
    -    # set k-points grid density
    +    # Set k-points grid density
         kpoints = kpoints_data()
         kpoints.set_kpoints_mesh(kmesh)
         inputs.kpoints = kpoints
    -    # set parameters
    +    # Set parameters
         inputs.parameters = dict_data(dict=incar)
    -    # set potentials and their mapping
    +    # Set potentials and their mapping
         inputs.potential_family = Str(potential_family)
         inputs.potential_mapping = dict_data(dict=potential_mapping)
    -    # set options
    +    # Set options
         inputs.options = dict_data(dict=options)
    -    # set settings
    +    # Set settings
         inputs.settings = dict_data(dict=settings)
    -    # set workchain related inputs, in this case, give more explicit output to report
    +    # Set workchain related inputs, in this case, give more explicit output to report
         inputs.verbose = Bool(True)
    -    # submit the requested workchain with the supplied inputs
    -    submit(workchain, **inputs)
    +    # Submit the requested workchain with the supplied inputs
    +    results = run(workchain, **inputs)
    +    return results
     
     
     if __name__ == '__main__':
    @@ -102,7 +104,7 @@
         # potential instead
         POTENTIAL_MAPPING = {'Si': 'Si'}
     
    -    # jobfile equivalent
    +    # Jobfile equivalent
         # In options, we typically set scheduler options.
         # See https://aiida.readthedocs.io/projects/aiida-core/en/latest/scheduler/index.html
         # AttributeDict is just a special dictionary with the extra benefit that
    @@ -115,10 +117,19 @@
         OPTIONS.max_wallclock_seconds = 3600
         OPTIONS.max_memory_kb = 1024000
     
    +    EOS = []
         # Iterate over each lattice constant and pass it explicitly
         for lattice_constant in [3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3]:
             # POSCAR equivalent
             # Set the silicon structure
             STRUCTURE = get_structure(lattice_constant)
     
    -        main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)
    +        output = main(CODE_STRING, INCAR, KMESH, STRUCTURE, POTENTIAL_FAMILY, POTENTIAL_MAPPING, OPTIONS)
    +        # The output are stored as AiiDA datatypes, which is Dict in this case. To obtain a regular
    +        # dictionary, we use get_dict
    +        misc = output['misc'].get_dict()
    +        EOS.append([lattice_constant, misc['total_energies']['energy_extrapolated']])
    +    # Write volume and total energies to file
    +    with open('eos', 'w') as file_object:
    +        for item in EOS:
    +            file_object.write('{} {}\n'.format(item[0], item[1]))
    

Save the new file as run_fcc_si_multiple_volumes_eos.py or fetch it with:

  wget https://github.com/aiida-vasp/aiida-vasp/raw/develop/tutorials/run_fcc_si_multiple_volumes_eos.py

This file should enable automatic extraction of total energies versus the different volumes.
  1. Execute run_fcc_si_multiple_volumes_eos.py:

    %/$ python run_fcc_si_multiple_volumes_eos.py
    

    Give the script time to complete. Contrary to before you will now see the output in the terminal (the similar output we previously got with verdi process report for a given process). More on this later. A file eos should be generated.

  2. Inspect the eos file that was generated:

    %/$ more eos
    3.5 -4.42341934
    3.6 -4.66006377
    3.7 -4.79595549
    3.8 -4.86303425
    3.9 -4.87588342
    4.0 -4.8481406
    4.1 -4.78451894
    4.2 -4.69228806
    4.3 -4.58122037
    
  3. Plot the data with your favorite plotting tool, for instance Gnuplot:

    %/$ gnuplot
    gnuplot> plot "eos" with lp
    

    which should give you something similar that is shown in the FCC Si tutorial.

We have now completed the FCC Si tutorial in AiiDA-VASP. However, maybe not in the best way. As you might already have realized, we had to change submit to run in order to make sure we got the output of each executed workchain. If we would have used submit we would only have access to a link to some output that we did not know was occupied or not at the point of inspection (unless we waited until we know the execution of the workchain was complete). In doing so, all the volume steps now run sequentially and thus it takes far longer for all the volume calculations to be finished (assuming your cluster started each volume calculations at the time of submission).

In the next tutorial we will instead develop a workchain to perform the different volume calculations (in fact, more general, namely different structure calculations) and call that from the run script. In doing so, you will realize that we can keep the benefits of automatic extraction of the data and parallel execution. In addition, we get the benefit of storing the result in the database.