Run Jobs with Python Scripts
****************************

PLAMS Python Multiple jobs Multiple molecules

AMS includes the PLAMS Python Library for Automating Molecular
Simulations. It takes care of input preparation, job execution, file
management and output processing as well as helps with building more
advanced data workflows.

[image]

To run a job with Python, we will

* create a Python file "jobname.py" containing the calculation
  settings and system, and

* run it with the "amspython" Python interpreter that is included with
  AMS.

You can also run jobs with Jupyter notebooks if you prefer.

Important:

  This tutorial requires that you run commands from the command-line,
  so make sure to familiarize yourself with the command-line and the
  AMS input format first.


Set up the job in AMSinput
==========================

Let's convert the Getting started: Geometry optimization of ethanol to
a PLAMS Python script. It is easiest to do this from AMSinput:

   Start AMSInput
   Choose [image: ADFPanel][image] in the yellow drop-down if not already active
   Set **Task** to **Geometry optimization**
   Set **Basis** to **SZ**
   Draw an ethanol molecule, or copy-paste the below coordinates into the molecule view

   C  0.01247  0.02254  1.08262
   C -0.00894 -0.01624 -0.43421
   H -0.49334  0.93505  1.44716
   H  1.05522  0.04512  1.44808
   H -0.64695 -1.12346  2.54219
   H  0.50112 -0.91640 -0.80440
   H  0.49999  0.86726 -0.84481
   H -1.04310 -0.02739 -0.80544
   O -0.66442 -1.15471  1.56909


View the run script in AMSinput
===============================

   In the panel bar, select **Details → Run script**

This shows the command and input to AMS that will be executed when the
job is run.

Here is an annotated version with the Python equivalents (we will
create the Python script in the next step):

   Task GeometryOptimization  # settings.input.ams.Task = "GeometryOptimization"
   System                     # system = ChemicalSystem(....)
       Atoms
           C 0.01247 0.02254 1.08262
           C -0.00894 -0.01624 -0.43421
           H -0.49334 0.93505 1.44716
           H 1.05522 0.04512 1.44808
           H -0.64695 -1.12346 2.54219
           H 0.50112 -0.9164 -0.8044
           H 0.49999 0.86726 -0.8448099999999999
           H -1.0431 -0.02739 -0.80544
           O -0.66442 -1.15471 1.56909
       End
       BondOrders
            1 2 1.0
            1 3 1.0
            1 4 1.0
            1 9 1.0
            2 6 1.0
            2 7 1.0
            2 8 1.0
            5 9 1.0
       End
   End

   Engine ADF       # settings.input.ADF
       Basis        # settings.input.ADF.Basis
           Type SZ  # settings.input.ADF.Basis.Type = "SZ"
       End
   EndEngine
   eor


Export a PLAMS Python script from AMSinput
==========================================

   **File → Export PLAMS script...**
   Export with the name "ethanol.py"

Now open "ethanol.py" in a text editor. The full contents can be
expanded below, but we will focus on specific parts of the script.

For simple examples, it is often easy to write the script directly
yourself. However, **Export PLAMS script** becomes especially helpful
for more complicated jobs, where the AMS input contains many blocks
and subblocks and it is harder to see how to construct the
corresponding "Settings" object in Python. In that case, exporting
from AMSinput gives you a correct working starting point that you can
then simplify or modify.

-[ Reveal full contents of ethanol.py ]-

   #!/usr/bin/env amspython
   """ Run this script with this command: $AMSBIN/amspython scriptname.py """
   from scm.base import ChemicalSystem
   from scm.plams import AMSJob, Settings, init, log
   from typing import Dict, Union


   def main():
       settings = get_settings()
       system = get_system()
       job_name = "ethanol"
       jobs = [AMSJob(settings=settings, molecule=system, name=job_name)]

       init(folder="plams_workdir")

       #use_parallel_jobrunner() # Uncomment this line to run each job on 1 core, but run many simultaneously. (Only useful if you run more than 1 job)

       for job in jobs:
           job.run()

       # === Start accessing results here ===

       for job in jobs:
           log(f"Job {job.name} has finished.")

       return jobs


   def get_settings() -> Settings:
       """Returns Settings for the AMSJob."""

       s = Settings()
       s.input.ams.Task = 'GeometryOptimization'
       s.input.ADF = Settings()
       s.input.ADF.Basis.Type = 'SZ'
       s.runscript.preamble_lines = []
       s.runscript.postamble_lines = []

       return s


   def get_system() -> Union[ChemicalSystem, Dict[str, ChemicalSystem], None]:
       """Returns None, a ChemicalSystem, or a dictionary of string keys and ChemicalSystem values"""
       from scm.base import InputFile

       system_input = """
           System
             Atoms
               C 0.01247 0.02254 1.08262
               C -0.00894 -0.01624 -0.43421
               H -0.49334 0.93505 1.44716
               H 1.05522 0.04512 1.44808
               H -0.64695 -1.12346 2.54219
               H 0.50112 -0.9164 -0.8044
               H 0.49999 0.86726 -0.8448099999999999
               H -1.0431 -0.02739 -0.80544
               O -0.66442 -1.15471 1.56909
             End
             BondOrders
               1 2 1.0
               1 3 1.0
               1 4 1.0
               1 9 1.0
               2 6 1.0
               2 7 1.0
               2 8 1.0
               5 9 1.0
             End
           End
           """.strip()

       if not system_input:
           return None

       systems = ChemicalSystem.all_from_input(InputFile("chemical_system", system_input))

       if len(systems) == 1 and "" in systems:
           return systems[""]
       else:
           return systems


   def use_parallel_jobrunner(maxjobs=None):
       from scm.plams import config, JobRunner

       if maxjobs is None:
           import multiprocessing

           maxjobs = multiprocessing.cpu_count()
       log(f"Running up to {maxjobs} jobs in parallel simultaneously")
       config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)
       config.job.runscript.nproc = 1


   if __name__ == "__main__":
       main()


Understand the exported script
------------------------------

If you are not used to Python, it helps to read the script from top to
bottom:

* import the PLAMS tools we need

* define some Python functions

* call "main()" at the very end

The exported script is written in a fairly general style, so it can
also handle more complicated jobs. For a first script, some parts look
more advanced than they really are.


Imports
-------

At the top of the script you will see lines like:

   from scm.base import ChemicalSystem
   from scm.plams import AMSJob, Settings, init, log

These lines make PLAMS classes and functions available in the script:

* "ChemicalSystem" stores the molecular structure

* "Settings" stores the AMS and engine input options

* "AMSJob" defines a calculation to run

* "init()" prepares the PLAMS working directory

* "log()" prints a message in the PLAMS output

The line

   from typing import Dict, Union

is only there for type hints. Type hints can help readability, but
they are not essential for understanding the script.


The "main()" function
---------------------

The central part of the script is:

   def main():
       settings = get_settings()
       system = get_system()
       job_name = "ethanol"
       jobs = [AMSJob(settings=settings, molecule=system, name=job_name)]

       init(folder="plams_workdir")

       for job in jobs:
           job.run()

       for job in jobs:
           log(f"Job {job.name} has finished.")

       return jobs

This is already a complete workflow:

* make the settings

* make the system

* create the job

* run the job

* print a message afterwards

Why is there a "main()" function at all? Only because it is good
practice for Python scripts. In a short script, you do not strictly
need it.


The settings part
-----------------

The exported "get_settings()" function contains:

   s = Settings()
   s.input.ams.Task = "GeometryOptimization"
   s.input.ADF = Settings()
   s.input.ADF.Basis.Type = "SZ"

This is the Python version of the AMS input blocks. You can often read
it almost literally:

* "s.input.ams.Task = "GeometryOptimization"" corresponds to "Task
  GeometryOptimization"

* "s.input.ADF.Basis.Type = "SZ"" corresponds to the ADF block with
  "Basis" and "Type SZ"

The nice thing about "Settings" is that you usually do not need to
create every block manually. For example, this also works:

   s = Settings()
   s.input.ams.Task = "GeometryOptimization"
   s.input.ADF.Basis.Type = "SZ"

PLAMS will create missing subblocks automatically when possible.

This is another reason why **Export PLAMS script** is so useful for
advanced jobs. For small examples, writing the "Settings" object by
hand is often easy. For larger jobs with many engine options,
constraints, properties, or multiple systems, the exported script
shows exactly how AMSinput translated the GUI setup into Python
assignments.

The two lines

   s.runscript.preamble_lines = []
   s.runscript.postamble_lines = []

are generated by AMSinput for completeness. For this example they are
empty, so you can simply remove them.


The system part
---------------

The exported script constructs the structure in a fairly general way:

   def get_system():
       from scm.base import InputFile
       system_input = """
       System
         ...
       End
       """.strip()
       systems = ChemicalSystem.all_from_input(InputFile("chemical_system", system_input))
       ...

This works for one system or for multiple named systems, but it is
more complicated than needed here.

For a single molecule, a simpler version is:

   system = ChemicalSystem("""
   System
     Atoms
       C 0.01247 0.02254 1.08262
       C -0.00894 -0.01624 -0.43421
       H -0.49334 0.93505 1.44716
       H 1.05522 0.04512 1.44808
       H -0.64695 -1.12346 2.54219
       H 0.50112 -0.91640 -0.80440
       H 0.49999 0.86726 -0.84481
       H -1.04310 -0.02739 -0.80544
       O -0.66442 -1.15471 1.56909
     End
     BondOrders
       1 2 1.0
       1 3 1.0
       1 4 1.0
       1 9 1.0
       2 6 1.0
       2 7 1.0
       2 8 1.0
       5 9 1.0
     End
   End
   """)

That is often the easiest form to read if you already know the AMS
"System" block.

If you prefer, you can also read the structure from an external file
instead of putting the coordinates directly in the script. That can
make the script shorter when the system is large.

If you want to learn more about what a "ChemicalSystem" can do, see
the ChemicalSystem overview.

If you want to see more examples of how structures and settings can be
written in Python, have a look at the PythonExamples section.


Functions are optional
----------------------

The exported script uses functions such as "get_settings()" and
"get_system()". They are useful, but not required.

For a short script, it is completely reasonable to write everything
directly in one place:

   #!/usr/bin/env amspython
   from scm.base import ChemicalSystem
   from scm.plams import AMSJob, Settings, init

   init(folder="plams_workdir")

   settings = Settings()
   settings.input.ams.Task = "GeometryOptimization"
   settings.input.ADF.Basis.Type = "SZ"

   system = ChemicalSystem("""
   System
     Atoms
       C 0.01247 0.02254 1.08262
       C -0.00894 -0.01624 -0.43421
       H -0.49334 0.93505 1.44716
       H 1.05522 0.04512 1.44808
       H -0.64695 -1.12346 2.54219
       H 0.50112 -0.91640 -0.80440
       H 0.49999 0.86726 -0.84481
       H -1.04310 -0.02739 -0.80544
       O -0.66442 -1.15471 1.56909
     End
     BondOrders
       1 2 1.0
       1 3 1.0
       1 4 1.0
       1 9 1.0
       2 6 1.0
       2 7 1.0
       2 8 1.0
       5 9 1.0
     End
   End
   """)

   job = AMSJob(settings=settings, molecule=system, name="ethanol")
   job.run()

This version does exactly the same calculation, but with less Python
syntax around it.

   Save the above script with the name "simplified_ethanol.py"

So when should you use functions?

* If the script is short, skipping functions is fine.

* If the script starts to grow, functions help split it into logical
  pieces.

* If you want to reuse the same setup many times, functions become
  very convenient.

In other words: functions are a tool for readability and reuse, not a
requirement of PLAMS.


The final two lines
-------------------

At the bottom of the exported script you will see:

   if __name__ == "__main__":
       main()

For now, you can read this simply as: "when this file is run as a
script, execute "main()"".

This is standard Python style for executable scripts.


Running the script
------------------

   Open a bash command-line
   Run the script with the AMS Python interpreter:

   $AMSBIN/amspython ethanol.py   # or simplified_ethanol.py

PLAMS will create a working directory called "plams_workdir" unless
you choose another name in "init(...)". Inside it, you will find a
subdirectory for the job, containing the input, output, and result
files. You can open those files in the GUI or text editors.


Running multiple systems
========================

One advantage of writing jobs in Python is that it becomes very easy
to repeat the same calculation for several molecules.

For example, you can create a list of systems from SMILES strings:

   from scm.base import ChemicalSystem

   systems = [
       ChemicalSystem.from_smiles("O"),
       ChemicalSystem.from_smiles("CO"),
       ChemicalSystem.from_smiles("CCO"),
   ]

Then run the same settings for all of them:

   from scm.plams import AMSJob, Settings, init

   init(folder="plams_workdir")

   settings = Settings()
   settings.input.ams.Task = "SinglePoint"
   settings.input.ADF.Basis.Type = "SZ"

   jobs = []
   for i, system in enumerate(systems, start=1):
       job = AMSJob(settings=settings, molecule=system, name=f"molecule_{i}")
       job.run()
       jobs.append(job)

After the jobs have finished, you can extract the energy in hartree:

   for job in jobs:
       energy = job.results.get_energy()
       print(job.name, energy)

Here is the full script:

-[ Reveal full multiple-molecule script ]-

   #!/usr/bin/env amspython
   from scm.base import ChemicalSystem
   from scm.plams import AMSJob, Settings, init

   init(folder="plams_workdir")

   systems = [
       ChemicalSystem.from_smiles("O"),
       ChemicalSystem.from_smiles("CO"),
       ChemicalSystem.from_smiles("CCO"),
   ]

   settings = Settings()
   settings.input.ams.Task = "SinglePoint"
   settings.input.ADF.Basis.Type = "SZ"

   jobs = []
   for i, system in enumerate(systems, start=1):
       job = AMSJob(settings=settings, molecule=system, name=f"molecule_{i}")
       job.run()
       jobs.append(job)

   for job in jobs:
       energy = job.results.get_energy()
       print(job.name, energy)

This kind of loop is one of the main reasons to use Python: once you
can run one job, running many similar jobs is only a small extra step.


What to remember
================

For a first PLAMS script, the most important ideas are:

* "Settings" is the Python version of the AMS input blocks

* "ChemicalSystem" stores the structure

* "AMSJob" combines settings and structure into a calculation

* helper functions are optional and mainly help keep larger scripts
  organized

Once you are comfortable with the short version above, you can
gradually add more Python features when they become useful.


Next steps
==========

* Run Jobs with Jupyter Notebooks

* To learn more about working with structures, see the ChemicalSystem
  overview.

* For many more ideas and ready-to-run examples, see Python Examples.
