Computers, Linux and compiled code

A machine that follows instructions

Every simulation, every neural network and every plot you will make runs on the same kind of machine, and it is simpler than it looks. A computer does exactly one thing, billions of times a second: it reads an instruction from memory, carries it out, and moves on to the next one. Everything else (operating systems, Python, GROMACS, PyTorch) is built from that.

The design almost every computer still follows was written down by John von Neumann in 1945, in a report on the EDVAC machine. Its central idea is the stored program: instructions are just numbers, kept in the same memory as the data they work on. A computer that can hold its own program can load a different program in seconds, which is what made general-purpose computing possible. The design has four parts:

• Memory holds numbers in numbered cells called addresses: both the program and its data.
• The processor (CPU) contains a control unit that steps through the program, an arithmetic logic unit (ALU) that adds, subtracts and compares, and a handful of very fast storage cells called registers.
• Buses carry addresses and numbers between them.
• Input and output connect it to disks, networks, screens and keyboards.

Watch one run

The machine below is a working computer with sixteen memory cells and six instructions. Its program, already in memory, adds up the numbers $1 + 2 + \dots + n$. Step through it one phase at a time and watch the three registers and the bus.

Every instruction goes through the same three-phase cycle.

1. Fetch. The program counter (PC) holds the address of the next instruction. The control unit reads that memory cell into the instruction register and adds one to the PC.
2. Decode. The control unit works out what the instruction asks for and which address it refers to.
3. Execute. The ALU does the arithmetic, or a number moves between memory and the accumulator register, or (for a jump) the PC is overwritten so that the next fetch comes from somewhere else.

That last trick, the conditional jump, is what makes loops and if statements possible. In the program above, JNZ 0 sends the machine back to address 0 as long as the counter is not zero. Every for loop in Python or do loop in Fortran becomes exactly this pattern after compilation.

The toy machine needs 3 cycles per instruction. A modern core finishes several instructions per cycle, by overlapping their phases like an assembly line (pipelining), running several at once (superscalar execution) and reordering them to avoid waiting (out-of-order execution). Its instructions are also wider: one SIMD (vector) instruction can multiply eight double-precision numbers at once.

From Python to machine code

Nobody writes these instructions by hand. A compiler (gfortran, gcc, nvcc) translates a whole program written in a high-level language into machine code before it runs; the statement sum = sum + n becomes something very like LOAD, ADD, STORE. An interpreter such as Python instead reads your program and acts on it statement by statement while it runs, which is flexible but far slower for arithmetic loops. That is why numerical Python hands its heavy work to compiled libraries such as NumPy, which you will see in the Python chapter.

How fast can a computer calculate?

Scientific computing counts floating-point operations per second (FLOP/s): additions and multiplications of real numbers. The peak rate of a processor is a simple product:

A 64-core server processor at 2.5 GHz doing 16 double-precision operations per cycle peaks at $64 \times 2.5\times10^9 \times 16 \approx 2.6\times10^{12}$ FLOP/s, 2.6 teraflops. A data-centre GPU, with thousands of simpler cores, reaches tens of teraflops in double precision and hundreds in the lower precisions used for machine learning. That is why both molecular dynamics and neural-network training moved to GPUs.

The memory wall

Here is the catch. Arithmetic is cheap; getting the numbers to the processor is not. Over the last forty years processors became far faster than memory, so a modern core can wait hundreds of cycles for a single number from main memory. Hardware hides this with a hierarchy of ever larger, ever slower memories:

Caches work because programs are predictable. They tend to reuse what they used recently (temporal locality) and to use neighbouring addresses next (spatial locality). Memory therefore moves in cache lines of 64 bytes, eight doubles at a time. A loop that walks through an array in storage order gets seven of every eight numbers for free; a loop that jumps around wastes most of each line. The Fortran chapter measures this: the same arithmetic runs 4.5 times slower with the wrong loop order.

Whether a calculation is limited by arithmetic or by memory depends on its arithmetic intensity, the number of operations it does per byte it moves. The roofline model puts it in one line:

Adding two long vectors ($c_i = a_i + b_i$) does 1 operation per 24 bytes moved, so it runs at a few percent of peak whatever the processor. Multiplying two large matrices reuses each number many times, has high intensity, and runs near peak. This one idea explains a great deal about scientific software: why libraries such as BLAS reorganise loops into blocks that fit in cache, why neighbour lists in molecular dynamics are sorted by space, and why neural networks are built from large matrix multiplications.

Many processors at once

Clock frequencies stopped rising around 2005 because of heat, so performance now comes from parallelism at every level:

Each level has a cost. Threads must not write to the same data at once; MPI processes must send messages over a network that is thousands of times slower than cache. Good parallel codes divide the work so that each processor mostly uses data it already has, a theme that returns in the chapters on clusters and molecular dynamics.

References and further reading

Everything is bits

Memory stores only two states per cell, written 0 and 1: a bit. Eight bits make a byte. Every number, letter, image and instruction is a pattern of bits, and the meaning of a pattern depends entirely on how a program agrees to read it. The same 32 bits can be an integer, a floating-point number, four letters of text or an instruction. Most of the numerical surprises in scientific computing come from forgetting how these patterns work.

Integers

Bits read as an ordinary binary number weigh powers of two, exactly as decimal digits weigh powers of ten:

To represent negative numbers, computers use two’s complement: the leftmost bit counts as $-2^{N-1}$ instead of $+2^{N-1}$. In 8 bits it is worth $-128$, so the range runs from $-128$ to $+127$. The beauty of the scheme is that the hardware adds signed and unsigned numbers with exactly the same circuit. The price is that arithmetic wraps around silently when it leaves the range: $127 + 1$ gives $-128$.

Real programs use 32-bit integers (up to about $2.1\times10^9$) or 64-bit ones (up to $9.2\times10^{18}$). The 32-bit limit is not academic: a simulation of more than two billion atoms, or a counter of steps × atoms, overflows it, and codes such as LAMMPS have build options for 64-bit indices for exactly this reason.

Floating-point numbers

Scientific quantities span enormous ranges, from $10^{-23}$ J for $k_B$ to $10^{23}$ for Avogadro’s number, so computers store real numbers in a binary form of scientific notation: a sign, a few significant digits, and an exponent. The universal standard is IEEE 754 (first published 1985). A 32-bit single-precision number is three fields:

Two patterns of the exponent are reserved: all zeros encode zero and the subnormal numbers that fill the gap near zero, and all ones encode $\pm\infty$ and NaN (not a number, the result of $0/0$ or $\sqrt{-1}$). The 64-bit double precision used for most scientific work has the same structure with 11 exponent bits and 52 fraction bits:

Why 0.1 is not exactly 0.1

Only fractions whose denominator is a power of two can be written exactly in binary. One tenth is $0.000110011001100\ldots_2$, repeating forever, so it is rounded to the nearest representable number. In the lab, type 0.1 in floating-point mode: the stored value is 0.100000001490116…, a relative error of about $1.5\times10^{-8}$.

The spacing between neighbouring floats grows with their size. The key number is machine epsilon, the gap between 1 and the next larger number:

So above $2^{24} = 16\,777\,216$, single precision cannot even represent every integer: type 16777217 in the lab and it becomes 16777216.

Where rounding bites

Each single error is tiny; the trouble is how they combine. Run this in the Python chapter or on your machine:

import numpy as np

print(0.1 + 0.2 == 0.3, 0.1 + 0.2)                 # the most famous surprise
print(f"{0.1:.20f}")                               # what is really stored
print(np.finfo(np.float64).eps, np.finfo(np.float32).eps)   # machine epsilon

big, small = 1.0e16, 1.0
print((big + small) - big)                         # small is lost entirely: catastrophic cancellation

x = np.full(10_000_000, 0.1, dtype=np.float32)
print("float32 running sum:", np.cumsum(x)[-1], "  pairwise sum:", x.sum(), "  exact: 1000000")

print(np.isclose(0.1 + 0.2, 0.3))                  # compare floats with a tolerance
print(np.int32(2**31 - 1) + np.int32(1))           # integer overflow wraps around

False 0.30000000000000004
0.10000000000000000555
2.220446049250313e-16 1.1920929e-07
0.0
float32 running sum: 1.087937e+06   pairwise sum: 1.0000001e+06   exact: 1000000
True
-2147483648

Four lessons are hiding in that output.

1. Never test floats for equality. Compare with a tolerance: abs(a - b) < tol or np.isclose.
2. Subtracting nearly equal numbers destroys digits (catastrophic cancellation). $10^{16} + 1 - 10^{16}$ gives 0 because the 1 was rounded away before the subtraction. Rewrite formulas to avoid it: for example, compute $\ln(1+x)$ for small $x$ with log1p, not log(1 + x).
3. Summation order matters. Adding ten million single-precision 0.1s one after another gives 1,087,937, an error of 9%: once the running total is large, each small addition is rounded badly. NumPy’s sum adds in pairs, a tree, and is accurate to $10^{-7}$. Kahan’s compensated summation does better still. This is also why parallel sums give slightly different answers from serial ones (see the Fortran chapter): addition in floating point is not associative, $(a+b)+c \ne a+(b+c)$.
4. Integers overflow silently, and the sign can flip.

References and further reading

Why there is an operating system at all

The previous chapters showed a processor executing one program. A real computer runs hundreds at once: your shell, an editor, a simulation on 64 cores, a network daemon, someone else’s job. Each believes it has the machine to itself. Something must share out the processor time, keep the programs from overwriting each other’s memory, and give them a simple way to use disks, networks and GPUs without knowing how each device works. That something is the operating system, and its core, the part that is always running with full control of the hardware, is the kernel.

An operating system gives programs three illusions, its three great abstractions:

• a process: each program seems to have a processor of its own;
• virtual memory: each process seems to have a large, private memory of its own;
• files: storage, devices and even other programs look like named streams of bytes.

Unix, Linux and GNU

Unix was written at Bell Labs from 1969 by Ken Thompson and Dennis Ritchie, who also created the C language to rewrite it portably. Its ideas (a hierarchical file system, small programs connected by pipes, everything as a file) shaped every system since. In 1991 Linus Torvalds, then a student in Helsinki, wrote a new Unix-like kernel, Linux, and released it freely. Combined with the GNU project’s compiler, C library and command-line tools, it became a complete free operating system. Today Linux runs Android phones, most web servers, and every one of the world’s 500 fastest supercomputers. A distribution (Ubuntu, Rocky, Debian, SUSE) packages the kernel with tools and an installer; clusters usually run an enterprise distribution such as Rocky Linux or RHEL.

The layers

Linux is built in layers, each using only the one below it. Click any block to read what it does, then pick a command and step through what happens when you press Enter.

User mode and kernel mode. The processor itself enforces the division. Ordinary programs run in user mode, in which instructions that touch hardware or other programs’ memory are forbidden; if a program tries, the processor stops it and the kernel kills it (the familiar segmentation fault). Only the kernel runs in kernel mode. A program that needs something done, such as reading a file, asks the kernel through a system call: a special instruction that switches the processor into kernel mode at a fixed entry point, where the kernel checks the request, does the work and switches back. Linux has about 350 system calls; a handful (open, read, write, close, fork, execve, exit, mmap) do most of the work. You can watch them: on Linux, strace -c ls counts every system call ls makes.

Processes

A process is a running program: its code, its memory, its open files and the processor state (registers, program counter) needed to resume it. Each has a number, the PID. New processes are created in a two-step dance that is pure Unix:

1. fork() makes an almost exact copy of the calling process;
2. the copy calls execve() to replace its program with a different one.

That is how your shell starts every command, and why processes form a family tree (pstree shows it), all descended from process 1, systemd, which the kernel starts at boot. The scheduler decides which process runs on which core, switching between them many times a second, so that hundreds share a few cores fairly. A process can be running, waiting (for disk, network or a child), or stopped. Processes are controlled with signals: Ctrl-C sends SIGINT, kill PID sends SIGTERM (please stop), kill -9 PID sends SIGKILL (stop now, no clean-up).

Virtual memory

Every process sees its own private address space, as if it had the machine’s whole memory to itself. The kernel and the processor’s memory management unit translate these virtual addresses to physical ones, in blocks called pages (usually 4 KB), using a page table per process.

The consequences are everywhere in practice. Processes cannot read each other’s memory, so a crashing program cannot corrupt another. Memory is assigned lazily: a program that allocates 100 GB but touches only 1 GB uses only 1 GB of RAM. Pages not used for a while can be moved to disk (swap), which is why an overloaded machine suddenly becomes thousands of times slower; clusters usually disable swap and kill jobs that exceed their memory request instead. And the kernel uses spare RAM as a page cache for files, which is why reading the same trajectory a second time is much faster than the first.

Files, and everything else that looks like one

To Linux, a file is a named sequence of bytes. Directories are files listing names; so are devices (/dev/nvidia0), pipes between processes, network sockets, and live information about the system (/proc/cpuinfo, /proc/meminfo, /proc/PID/status). One set of system calls (open, read, write, close) works on all of them. This is the famous Unix idea that everything is a file.

A running process refers to its open files by small integers, file descriptors. Three are open from the start: 0 standard input, 1 standard output, 2 standard error. Redirection (>) and pipes (|), which the next chapter dissects, are nothing more than the shell rearranging these three numbers before it starts a program.

Each file’s metadata (owner, permissions, size, times, where its data lives on disk) is kept in an inode; the name you see is just a directory entry pointing to it. Permissions are three sets of read, write and execute bits for the owner, the group and everyone else, and the special user root (the administrator) bypasses them all. On a shared cluster you are never root, which is why software comes from modules or from environments you install in your own directories.

The directory tree

Linux has a single tree rooted at /; disks and network file systems are mounted at points inside it rather than having drive letters. The standard layout:

Looking inside your own machine

These commands work on any Linux system (macOS has some under different names):

uname -a                 # kernel name and version
lscpu                    # processor: cores, threads, caches, vector extensions
free -h                  # memory in use, free and used as page cache
cat /proc/meminfo | head # the kernel's own memory accounting
df -h                    # mounted file systems and free space
ps -ef --forest | head   # the process family tree
top                      # live view of processes (q to quit)
strace -c python3 -c 'print(1)' 2>&1 | tail -15   # system calls of a tiny Python run
ldd $(which python3)     # the shared libraries a program depends on

References and further reading

Why the shell, in 2026

Nearly every scientific computer you will work on, from a cluster with ten thousand cores to a GPU workstation or a cloud machine, runs Linux and is reached through a text shell. Simulation codes are launched from it, job schedulers are driven from it, and the files they produce, often thousands of them, are far quicker to inspect, filter and move with a few commands than with any graphical tool. The shell also makes work reproducible: a command you typed can be saved in a script and run again next year, on another machine, with identical results.

A few terms. The terminal is the window; the shell is the program inside it that reads your commands and runs them. The common shells are bash (the default on most Linux systems) and zsh (the default on macOS); everything in this lesson works in both. macOS’s tools are BSD versions that occasionally take slightly different options from the GNU versions on Linux, which is noted where it matters.

The file system is a tree

Every file lives in one tree that starts at the root, written /. Your own files live in your home directory, written ~ (for example /home/student). A path names a location in the tree:

• an absolute path starts at the root: /home/student/project/input/md.mdp;
• a relative path starts at your current directory: input/md.mdp;
• . means “this directory” and .. means “the directory above”;
• names starting with a dot, such as .bashrc, are hidden from a plain ls.

Avoid spaces in file names (run 1/ must then be quoted everywhere); use run_1/ or run-1/.

The anatomy of a command

ls -l -h runs/T300

A command is a program name (ls), then options that change its behaviour (-l long format, -h human-readable sizes; single-letter options can be combined as -lh), then arguments, usually files or directories (runs/T300). Long options start with two dashes (--human-readable). Every command documents itself: man ls opens its manual (press q to leave, /word to search), and most accept --help.

Practise here

The terminal below runs entirely in this page, on a pretend file system containing a small simulation project: three GROMACS runs at different temperatures, one of which crashed. Nothing you type can affect your computer. Work through the missions; each turns green when you have done it, and each has a solution if you are stuck. Type help to see what the practice shell understands, use Tab to complete names and ↑ to recall previous commands.

Moving around and looking

Reading ls -l:

-rw-r--r-- 1 student lab   412 Oct 14 09:30 md.mdp
drwxr-xr-x 1 student lab  4096 Oct 14 09:30 runs/

The first character is - for a file and d for a directory. The next nine are permissions in three groups of three, for the owner, the group and everyone else: r read, w write, x execute (for a directory, x means you may enter it). Then the owner, group, size in bytes, modification time and name. chmod +x script.sh makes a script executable; chmod 600 ~/.ssh/id_ed25519 makes a private key readable only by you, which SSH insists on.

Creating, copying, moving and deleting

mkdir -p analysis/rdf        # -p creates parent directories as needed
cp input/md.mdp md_new.mdp   # copy a file
cp -r input/ input_backup/   # copy a directory: -r means recursive
mv md_new.mdp md_T320.mdp    # rename (or move to another directory)
rm md_T320.mdp               # delete a file
rm -r old_runs/              # delete a directory and everything in it

Wildcards

The shell expands wildcards into matching file names before the command runs:

• * matches any sequence of characters: runs/*/md.log is every md.log one level below runs/;
• ? matches exactly one character: T3?0 matches T300, T320 and T340;
• [ab] matches one of the listed characters; {a,b} expands to both words (bash and zsh).

Because the shell does the expansion, grep never sees the *; it receives a list of file names. Put a pattern in quotes (find . -name '*.py') when you want the command to see it instead.

Searching

grep 'Fatal error' runs/T340/md.log       # lines containing the text
grep -n -i 'lincs' runs/*/md.log          # with line numbers, ignoring case, across files
grep -l 'Fatal error' runs/*/md.log       # only the names of matching files
grep -r 'ref_t' .                         # recursively through all files below .
grep -v '^[#@]' energy.xvg                # invert: lines NOT starting with # or @
find . -name '*.xtc' -size +1G            # files by name and size
find . -type d -name 'T3*'                # directories by name

grep patterns are regular expressions: ^ anchors the start of a line, $ the end, . matches any character, [0-9] any digit, * repeats the previous item. They are worth learning; they appear in every programming language.

Pipes and redirection: the real power

Every command reads from standard input and writes to standard output (and error messages to standard error). You can redirect them:

grep -v '^[#@]' energy.xvg > potential.dat     # > writes output to a file (overwriting it)
echo "T300 done" >> progress.log                # >> appends
gmx mdrun -deffnm md > run.out 2>&1             # 2>&1 sends errors to the same file

and connect commands with a pipe, |, so that the output of one becomes the input of the next:

grep -v '^[#@]' energy.xvg | wc -l                          # how many data lines?
grep -h 'Performance' runs/*/md.log | sort -k2 -n | head -1  # slowest run
cut -c18-20 protein.pdb | sort | uniq -c | sort -rn         # residue counts in a PDB file
awk '!/^[#@]/ {print $1, $2}' energy.xvg > energy.dat       # strip headers, keep two columns
awk '!/^[#@]/ {s += $2; n++} END {print s / n}' energy.xvg  # mean of column 2

This is the Unix philosophy, set out by Doug McIlroy: small programs that each do one thing well, combined through a universal interface, text. A pipeline of four small commands often replaces a page of code.

Your first script

A shell script is a file of commands. This one prepares three simulations at different temperatures from one template:

#!/bin/bash
# make_runs.sh: one directory and one parameter file per temperature
set -euo pipefail                      # stop on errors, unset variables and failed pipes

for T in 300 320 340; do
    dir="T${T}"
    mkdir -p "$dir"
    # replace the reference temperature, writing a new file in each directory
    sed "s/^ref_t .*/ref_t      = ${T}/" md.mdp > "$dir/md.mdp"
    echo "prepared $dir: $(grep ref_t "$dir/md.mdp")"
done

chmod +x make_runs.sh
./make_runs.sh

prepared T300: ref_t      = 300
prepared T320: ref_t      = 320
prepared T340: ref_t      = 340

Points to notice: the first line (the shebang) says which program runs the script; variables are set with T=300 (no spaces) and used as $T or ${T}; double quotes around "$dir" keep the value as one word even if it contains spaces; $(…) inserts the output of a command; and set -euo pipefail makes the script stop at the first error instead of carrying on with broken state. sed "s/old/new/" substitutes text; on Linux sed -i edits a file in place, while macOS’s sed needs sed -i ''.

References and further reading

The shell is a program like any other

bash is an ordinary user program (it lives at /bin/bash) whose whole job is a loop: print a prompt, read a line, work out what it means, start the programs it names, wait for them, and repeat. Understanding the middle step, working out what a line means, removes almost all of the mystery from the command line: why spaces in file names break scripts, why * sometimes stays a star, why cd cannot be a separate program, and why 2>&1 must come after >.

Four stages for every line

1. Splitting into tokens

bash first cuts the line into words and operators. Operators are the characters with special meaning: |, &, ;, &&, ||, <, >, >>, parentheses. Quotes stop splitting: "my run" is one word, my run is two. Single quotes switch off every special meaning inside them; double quotes switch off splitting and wildcards but still allow $ expansions.

2. Expansions, in a fixed order

Each word is then rewritten, always in the same order:

Two consequences are worth remembering. Word splitting and globbing happen after variable expansion, and only to unquoted results, which is why variables should almost always be quoted: "$file". And a glob that matches nothing is passed on unchanged, so ls *.xyz complains about a file literally called *.xyz.

3. Redirections, pipes and starting programs

Before running a command, bash sets up its file descriptors (the chapter on the Linux system): 0 for input, 1 for output, 2 for errors.

Redirections are processed left to right. cmd > log 2>&1 first points 1 at log, then makes 2 point wherever 1 now points: both go to the file. cmd 2>&1 > log makes 2 point at the terminal (where 1 pointed at that moment) and only then moves 1 to the file, so errors still appear on screen.

A pipe a | b asks the kernel for a small buffer with two ends, connects a’s descriptor 1 to the writing end and b’s descriptor 0 to the reading end, and starts both programs at once. Data streams through while both run: grep can be filtering the first lines of a huge file before cat has finished reading it, and the file never has to fit in memory. When b stops reading (as head does after ten lines), a is told to stop.

Then bash runs the command. If the first word is a built-in such as cd, export, echo or source, bash does the work itself. cd must be a built-in: it changes the current directory of the shell process itself, and a separate program could only change its own. Otherwise bash searches the directories listed in $PATH, in order, for an executable with that name (type gmx or which gmx shows which one wins) and starts it with fork and exec.

4. Exit status

Every program ends by returning a small number, its exit status: 0 for success, anything from 1 to 255 for failure. bash keeps the last one in $?. The operators && and || act on it: make && ./run only runs if the build succeeded; grep -q Fatal md.log || echo "no errors" prints only if grep found nothing.

Variables and the environment

T=320 sets a shell variable, visible only to this shell. export T makes it an environment variable, copied into every program the shell starts; that is how settings such as OMP_NUM_THREADS, CUDA_VISIBLE_DEVICES and PATH reach your programs. Prefixing a single command, OMP_NUM_THREADS=8 ./md, sets the variable for that command alone.

At start-up an interactive bash reads ~/.bashrc (login shells on some systems read ~/.bash_profile instead). That is where aliases, module load lines and PATH additions go, which also means a broken .bashrc can break every job you submit.

Traps that catch everyone

• Spaces around =. T = 320 runs a command called T with arguments = and 320. Assignments have no spaces: T=320.
• Single against double quotes. echo '$HOME' prints $HOME; echo "$HOME" prints /home/student.
• A pipeline’s status is its last command’s. false | true succeeds. set -o pipefail makes a pipeline fail if any part fails.
• ls output is not for parsing. Loop over a glob instead: for f in *.log; do …; done.
• Windows line endings. A script edited on Windows may end lines with \r\n; bash then reports $'\r': command not found. Convert with dos2unix.

References and further reading

How a cluster is organised

A high-performance computing (HPC) cluster is many computers (nodes) connected by a fast network and sharing file systems. You log in to a login node, which is shared by everyone and meant for editing, compiling and submitting work. The calculations run on compute nodes, which you reach only through the scheduler (usually Slurm, sometimes PBS or LSF): you describe the resources you need in a job script, submit it, and the scheduler starts it when the resources are free.

Logging in with SSH, properly

ssh user@cluster.university.edu opens a remote shell. Typing a password every time is slow and less secure than a key pair: a private key that never leaves your computer and a public key that you install on the server.

ssh-keygen -t ed25519 -C "laptop key"        # creates ~/.ssh/id_ed25519 and id_ed25519.pub
ssh-copy-id user@cluster.university.edu       # appends the public key to the server's authorized_keys

Then give the cluster a short name in ~/.ssh/config:

Host cluster
    HostName login.cluster.university.edu
    User yourname
    IdentityFile ~/.ssh/id_ed25519
    ProxyJump gateway.university.edu      # only if you must go through a gateway first

Now ssh cluster, scp file cluster: and rsync … cluster: all work, and so does the remote mode of editors such as VS Code. Protect the private key with a passphrase and keep it readable only by you (chmod 600).

Port forwarding lets you use a remote service in your local browser, for example a Jupyter notebook running on the cluster:

ssh -L 8888:localhost:8888 cluster    # then open http://localhost:8888 on your laptop

Moving data

scp copies files over SSH; rsync synchronises directories and transfers only what has changed, which makes it the right tool for simulation data:

rsync -avzP cluster:~/project/runs/ runs/                     # fetch, preserving times; resumable
rsync -avzP --exclude '*.trr' --exclude '*.cpt' cluster:~/project/runs/ runs/   # skip huge files
rsync -avn --delete runs/ cluster:~/backup/runs/              # -n: dry run first, ALWAYS with --delete

A trailing slash on the source (runs/) copies the contents of the directory; without it (runs) it copies the directory itself. Pack many small files before moving or archiving them, and check that the copy is intact:

tar -czf runs_2026-10.tar.gz runs/          # c create, z gzip, f file name
sha256sum runs_2026-10.tar.gz > runs_2026-10.sha256   # on macOS: shasum -a 256
sha256sum -c runs_2026-10.sha256            # verify after transfer
tar -xzf runs_2026-10.tar.gz                # x extract

Sessions that survive a dropped connection

If your SSH connection drops, programs started in it are killed. tmux (or the older screen) keeps a session running on the server:

tmux new -s analysis       # start a named session
# ... work; press Ctrl-b then d to detach ...
tmux ls                    # list sessions after logging back in
tmux attach -t analysis    # reattach

For a single long command, nohup command > out.log 2>&1 & also survives logout, but on a cluster long runs should be batch jobs.

Software: modules and environments

Clusters provide many versions of compilers, MPI libraries and simulation codes through environment modules:

module avail gromacs        # which versions exist?
module load gromacs/2023.3  # put one on your PATH
module list                 # what is loaded
module purge                # unload everything (start job scripts with this)

For Python, keep each project in its own environment so that package versions do not collide:

conda create -n mlip python=3.11 numpy scipy matplotlib
conda activate mlip
pip install mdanalysis ase
conda env export > environment.yml       # record exact versions for reproducibility

Job scripts

A Slurm job script is a shell script whose #SBATCH comment lines request resources. The builder below writes one for common cases and explains every line; adjust the module names and partition to your cluster (sinfo -s lists partitions, module avail the software).

The essentials of resource requests:

• Nodes, tasks and cores. --ntasks is the number of MPI processes; --cpus-per-task is the number of threads each one runs (OpenMP); their product per node should not exceed the cores in a node. Set OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK so programs use exactly what was allocated.
• GPUs. Request them explicitly (--gres=gpu:1 or --gpus-per-node=1, depending on the site). GPU codes usually want one to a few MPI ranks per GPU.
• Time. The job is killed at the limit, so ask for a little more than you need, and write checkpoints: GROMACS -cpi, LAMMPS write_restart. Shorter jobs start sooner because the scheduler can fit them into gaps (backfilling).
• Memory. Too little and the job is killed; far too much and it waits longer. After a job, sacct -j JOBID --format=JobID,Elapsed,MaxRSS,State shows what it actually used.

Monitoring and controlling jobs:

sbatch job.sh                                  # submit; prints the job ID
squeue -u $USER                                # your jobs: PD pending, R running
squeue -u $USER --start                        # estimated start times
scontrol show job 123456                       # full details, including why it is pending
scancel 123456                                 # cancel
sbatch --dependency=afterok:123456 analyse.sh  # start only after job 123456 succeeds
sacct -u $USER -S 2026-10-01 --format=JobID,JobName,Elapsed,State,ExitCode   # history

Job arrays run the same script many times with different inputs, with the index in $SLURM_ARRAY_TASK_ID: the natural way to run a temperature series or a set of independent replicas (tick the option in the builder).

Where to keep files

Clusters usually have several file systems with different purposes:

Check your usage with quota, du -sh * or the site’s own tool, and copy anything you need to keep off scratch promptly.

Explain a command

When you meet an unfamiliar command line in a tutorial or a colleague’s script, paste it here to see what every part does.

Reference

References and further reading

Why learn Fortran now?

Fortran (FORmula TRANslation, 1957) is the oldest language still in heavy scientific use, and for a good reason: it was designed for one job, fast numerical work on arrays, and compilers have had seventy years to get good at it. Much of the code a computational materials scientist runs every day is written in it: VASP, Quantum ESPRESSO, CP2K, CASTEP, SIESTA and ABINIT for electronic structure, and many weather, climate and engineering codes. Reading, modifying and compiling Fortran is part of the job, even if your own analysis is in Python.

The language you will meet is modern Fortran (Fortran 90 and later; the current standard is Fortran 2023). It has modules, dynamic arrays, whole-array arithmetic, derived types and built-in parallelism, and it bears little resemblance to the fixed-column FORTRAN 77 of punched cards, which you will still find in older libraries (files ending .f or .F instead of .f90).

A first program

Save this as hello.f90, then compile and run it:

program hello
   implicit none
   print *, "Hello from Fortran"
end program hello

gfortran hello.f90 -o hello     # compile and link into an executable called hello
./hello

 Hello from Fortran

implicit none belongs in every program and module. Without it, Fortran silently declares any undeclared variable whose name starts with i–n as an integer and anything else as a real, so a typo creates a new variable instead of an error.

Types, kinds and the precision trap

Numerical code needs to say how many bits a number has. The portable way is to take kind parameters from the standard module iso_fortran_env:

program basics
   use, intrinsic :: iso_fortran_env, only: dp => real64
   implicit none
   real(dp), parameter :: kB = 1.380649e-23_dp     ! J/K, exact since 2019
   real(dp) :: T, E, a, b
   integer  :: ndof
   character(len=*), parameter :: name = "argon"

   T = 300.0_dp
   ndof = 3
   E = 0.5_dp * ndof * kB * T                      ! equipartition: (3/2) kB T
   print '(a, es12.5, a)', "mean kinetic energy of one " // name // " atom: ", E, " J"

   a = 0.1                                         ! a single-precision literal...
   b = 0.1_dp                                      ! ...and a double-precision one
   print '(a, f20.17)', "0.1     stored in a real(dp): ", a
   print '(a, f20.17)', "0.1_dp  stored in a real(dp): ", b
   print '(a, i0)', "integer division 7/2 = ", 7 / 2
end program basics

basics.f90:18:48:
   18 |    print '(a, i0)', "integer division 7/2 = ", 7 / 2
      |                                                1
Warning: Integer division truncated to constant '3' at (1) [-Winteger-division]

mean kinetic energy of one argon atom:  6.21292E-21 J
0.1     stored in a real(dp):  0.10000000149011612
0.1_dp  stored in a real(dp):  0.10000000000000001
integer division 7/2 = 3

Two classic bugs are visible here.

• Literal precision. 0.1 is a single-precision constant, accurate to about 7 digits, even when assigned to a double-precision variable. Write 0.1_dp. A force field parameter typed without its kind suffix silently loses half its digits.
• Integer division. 7/2 is 3. Write 7.0_dp/2 when you mean 3.5. The compiler warned us because -Wall was on, one of many reasons to always compile with warnings.

Arrays: the reason Fortran exists

Arrays are first-class. Whole-array expressions, slices and a large library of intrinsic functions replace most explicit loops:

program arrays
   use, intrinsic :: iso_fortran_env, only: dp => real64
   implicit none
   real(dp), allocatable :: x(:), y(:), A(:, :)
   integer :: i, n

   n = 5
   allocate(x(n), y(n), A(3, 3))
   x = [(real(i, dp), i = 1, n)]                   ! array constructor: 1, 2, ..., 5
   y = x**2                                        ! element-wise: no loop needed
   print '(a, 5f8.2)', "x            =", x
   print '(a, 5f8.2)', "x**2         =", y
   print '(a, f8.2)',  "sum(x**2)    =", sum(y)
   print '(a, f8.2)',  "dot(x, x)    =", dot_product(x, x)
   print '(a, 3f8.2)', "x(2:4)       =", x(2:4)    ! a slice

   A = reshape([(real(i, dp), i = 1, 9)], [3, 3])  ! filled column by column
   print '(a)', "A, printed row by row:"
   do i = 1, 3
      print '(3f6.1)', A(i, :)
   end do
   print '(a, 3f8.2)', "A @ [1,1,1]  =", matmul(A, [1.0_dp, 1.0_dp, 1.0_dp])

   where (y > 10.0_dp) y = 10.0_dp                 ! masked assignment
   print '(a, 5f8.2)', "capped at 10 =", y
   deallocate(x, y, A)
end program arrays

x            =    1.00    2.00    3.00    4.00    5.00
x**2         =    1.00    4.00    9.00   16.00   25.00
sum(x**2)    =   55.00
dot(x, x)    =   55.00
x(2:4)       =    2.00    3.00    4.00
A, printed row by row:
   1.0   4.0   7.0
   2.0   5.0   8.0
   3.0   6.0   9.0
A @ [1,1,1]  =   12.00   15.00   18.00
capped at 10 =    1.00    4.00    9.00   10.00   10.00

Notice two things about A. Indices start at 1 by default (you can declare other bounds, such as real :: psi(-10:10)). And reshape filled it column by column: the numbers 1, 2, 3 went down the first column. That is not a quirk of reshape; it is how Fortran stores every multidimensional array.

Memory layout and loop order

Memory is one-dimensional, so a 2D array must be flattened. Fortran stores arrays in column-major order: A(1,1), A(2,1), A(3,1), …, A(1,2), …, so the first index varies fastest. C, C++ and NumPy (by default) use row-major order, where the last index varies fastest.

This matters for speed. The processor never fetches a single number; it loads a whole cache line (64 bytes, eight doubles) at once, and it is fastest when consecutive iterations touch consecutive addresses. The lab shows which loop order does that.

The rule: in Fortran, make the innermost loop run over the first index. On a real machine the effect is large. This program performs exactly the same arithmetic twice, on a 4000 × 4000 array:

program loop_order
   use, intrinsic :: iso_fortran_env, only: dp => real64, int64
   implicit none
   integer, parameter :: n = 4000
   real(dp), allocatable :: a(:, :)
   integer :: i, j, rep
   integer(int64) :: t0, t1, rate

   allocate(a(n, n))
   a = 1.0_dp

   call system_clock(t0, rate)
   do rep = 1, 5
      do j = 1, n                  ! outer loop over columns
         do i = 1, n               ! inner loop over rows: consecutive addresses
            a(i, j) = a(i, j) * 1.000001_dp
         end do
      end do
   end do
   call system_clock(t1)
   print '(a, f7.3, a)', "inner loop over i (contiguous): ", real(t1 - t0, dp) / rate, " s"

   call system_clock(t0)
   do rep = 1, 5
      do i = 1, n                  ! outer loop over rows
         do j = 1, n               ! inner loop over columns: jumps of n*8 bytes
            a(i, j) = a(i, j) * 1.000001_dp
         end do
      end do
   end do
   call system_clock(t1)
   print '(a, f7.3, a)', "inner loop over j (strided):    ", real(t1 - t0, dp) / rate, " s"
   print '(a, f12.6)', "checksum: ", sum(a) / n**2
end program loop_order

gfortran -O2 loop_order.f90 -o loop_order && ./loop_order

inner loop over i (contiguous):   0.083 s
inner loop over j (strided):      0.390 s
checksum:     1.000010

The strided version is about 4.5 times slower on this laptop (an Apple M1), only because of the order in which memory is touched. Whole-array syntax (a = a * 1.000001_dp) lets the compiler pick the right order for you.

Procedures and modules

Real programs are split into subroutines (which act through their arguments) and functions (which return a value), placed inside modules. A module gives the compiler an explicit interface: it can check every call for the right number, type and shape of arguments. The intent attribute documents and enforces which arguments are inputs and outputs.

Here is a Lennard-Jones energy and force routine with periodic boundaries, in a module, the same physics as the molecular dynamics lesson:

module lj_mod
   use, intrinsic :: iso_fortran_env, only: dp => real64
   implicit none
   private
   public :: dp, lj_energy_forces

contains

   subroutine lj_energy_forces(x, box, rc, energy, f)
      real(dp), intent(in)  :: x(:, :)       ! positions, shape (3, n), reduced units
      real(dp), intent(in)  :: box, rc       ! cubic box length and cut-off
      real(dp), intent(out) :: energy
      real(dp), intent(out) :: f(:, :)       ! forces, shape (3, n)
      integer  :: i, j, n
      real(dp) :: d(3), r2, ir6, ff, ec

      n = size(x, 2)
      energy = 0.0_dp
      f = 0.0_dp
      ec = 4.0_dp * (rc**(-12) - rc**(-6))   ! shift so that the energy is zero at rc
      do i = 1, n - 1
         do j = i + 1, n
            d = x(:, i) - x(:, j)
            d = d - box * anint(d / box)     ! minimum-image convention
            r2 = sum(d * d)
            if (r2 < rc * rc) then
               ir6 = 1.0_dp / r2**3
               energy = energy + 4.0_dp * ir6 * (ir6 - 1.0_dp) - ec
               ff = 24.0_dp * ir6 * (2.0_dp * ir6 - 1.0_dp) / r2
               f(:, i) = f(:, i) + ff * d
               f(:, j) = f(:, j) - ff * d
            end if
         end do
      end do
   end subroutine lj_energy_forces

end module lj_mod

and a main program that builds a 256-atom fcc crystal and calls it:

program lj
   use lj_mod, only: dp, lj_energy_forces
   implicit none
   integer, parameter :: m = 4, n = 4 * m**3          ! 4 x 4 x 4 fcc cells: 256 atoms
   real(dp), parameter :: rho = 0.8442_dp
   real(dp), parameter :: basis(3, 4) = reshape([0.0_dp, 0.0_dp, 0.0_dp,  0.5_dp, 0.5_dp, 0.0_dp, &
                                                 0.5_dp, 0.0_dp, 0.5_dp,  0.0_dp, 0.5_dp, 0.5_dp], [3, 4])
   real(dp) :: x(3, n), f(3, n), a, box, e
   integer  :: i, j, k, b, p

   a = (4.0_dp / rho)**(1.0_dp / 3.0_dp)               ! lattice constant for this density
   box = m * a
   p = 0
   do k = 0, m - 1
      do j = 0, m - 1
         do i = 0, m - 1
            do b = 1, 4
               p = p + 1
               x(:, p) = a * ([real(i, dp), real(j, dp), real(k, dp)] + basis(:, b))
            end do
         end do
      end do
   end do

   call lj_energy_forces(x, box, 2.5_dp, e, f)
   print '(a, i0, a, f6.4)', "atoms: ", n, ", fcc lattice at density ", rho
   print '(a, f10.5)', "potential energy per atom: ", e / n
   print '(a, es10.2)', "largest force component (zero by symmetry): ", maxval(abs(f))
end program lj

gfortran -O2 -Wall -c lj_mod.f90      # produces lj_mod.o and lj_mod.mod
gfortran -O2 -Wall -c lj_main.f90     # needs lj_mod.mod, so compile the module first
gfortran -o lj lj_mod.o lj_main.o     # link
./lj

atoms: 256, fcc lattice at density 0.8442
potential energy per atom:   -6.33281
largest force component (zero by symmetry):   2.42E-14

On a perfect lattice every atom’s force cancels by symmetry, so the largest force is rounding noise: a free correctness test of the force code.

Compiling: the flags that matter

Develop with the debug set and run production with the optimised set. Most “the code gives garbage” problems in scientific Fortran are out-of-bounds array accesses that -fcheck=all turns into an immediate, located error message.

A Makefile records how to build a program, rebuilds only what changed, and handles the module order. This one builds the Lennard-Jones example (the indented lines must start with a tab):

# Makefile for the Lennard-Jones example. Usage: make, make debug, make clean
FC      = gfortran
FFLAGS  = -O2 -march=native -Wall
DEBUG   = -O0 -g -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow -Wall -Wextra
OBJ     = lj_mod.o lj_main.o

lj: $(OBJ)
	$(FC) $(FFLAGS) -o $@ $(OBJ)

# a module must be compiled before any file that uses it
lj_main.o: lj_mod.o

%.o: %.f90
	$(FC) $(FFLAGS) -c $<

debug: FFLAGS = $(DEBUG)
debug: clean lj

clean:
	rm -f *.o *.mod lj

.PHONY: debug clean

make builds an optimised lj; make debug rebuilds it with all runtime checks; make clean removes the build files. Larger projects use CMake or the Fortran Package Manager (fpm), which find dependencies and module order automatically.

Parallel Fortran

OpenMP parallelises loops over the cores of one node with compiler directives, comment lines that start with !$omp:

program omp_sum
   use, intrinsic :: iso_fortran_env, only: dp => real64
   use omp_lib
   implicit none
   integer, parameter :: n = 200000000
   integer  :: i
   real(dp) :: s, t0

   s = 0.0_dp
   t0 = omp_get_wtime()
   !$omp parallel do reduction(+:s)
   do i = 1, n
      s = s + 1.0_dp / real(i, dp)**2       ! each thread sums its share; reduction combines them
   end do
   !$omp end parallel do
   print '(a, i0, a, f6.3, a)', "threads: ", omp_get_max_threads(), ",  time: ", omp_get_wtime() - t0, " s"
   print '(a, f12.9, a, f12.9)', "sum 1/i^2 = ", s, "   pi^2/6 = ", acos(-1.0_dp)**2 / 6
end program omp_sum

gfortran -O2 -fopenmp omp_sum.f90 -o omp_sum
OMP_NUM_THREADS=1 ./omp_sum
OMP_NUM_THREADS=8 ./omp_sum

threads: 1,  time:  0.569 s
sum 1/i^2 =  1.644934058   pi^2/6 =  1.644934067
threads: 8,  time:  0.091 s
sum 1/i^2 =  1.644934062   pi^2/6 =  1.644934067

Eight threads give a 6.2-fold speed-up. Look closely at the sums: they differ in the ninth digit. Each thread adds its own share and the partial sums are combined at the end, so the order of additions changes, and floating-point addition is not associative. Parallel results differing in the last digits from serial ones is normal; differing in the third digit is a bug, usually a variable that should have been private or a missing reduction.

MPI (message passing) parallelises across many nodes. Each process owns part of the data and exchanges messages explicitly; this is how VASP, Quantum ESPRESSO and CP2K run on thousands of cores. Programs are compiled with a wrapper and launched with mpirun or, inside a Slurm job, srun:

mpif90 -O2 my_mpi_code.f90 -o my_mpi_code
mpirun -np 64 ./my_mpi_code

Do not write your own linear algebra

BLAS (vector and matrix operations) and LAPACK (linear systems, eigenvalues, factorisations) are decades-tested and heavily optimised (OpenBLAS, Intel MKL, Apple Accelerate). Calling them is easy:

program solve
   use, intrinsic :: iso_fortran_env, only: dp => real64
   implicit none
   integer, parameter :: n = 3
   real(dp) :: A(n, n), b(n)
   integer  :: ipiv(n), info
   external :: dgesv                                  ! from LAPACK

   A = reshape([4.0_dp, 1.0_dp, 2.0_dp,  1.0_dp, 5.0_dp, 1.0_dp,  2.0_dp, 1.0_dp, 6.0_dp], [n, n])
   b = [1.0_dp, 2.0_dp, 3.0_dp]
   call dgesv(n, 1, A, n, ipiv, b, n, info)           ! solves A x = b; b is overwritten by x
   if (info /= 0) error stop "dgesv failed"
   print '(a, 3f10.6)', "x = ", b
end program solve

gfortran -O2 solve.f90 -o solve -llapack     # on HPC systems often -lopenblas, or MKL flags
./solve

x =  -0.063830  0.319149  0.468085

The same answer as numpy.linalg.solve, which calls the same LAPACK routine underneath.

Calling Fortran from Python

f2py, part of NumPy, wraps Fortran routines as a Python module, the best of both worlds: Python for workflow and analysis, Fortran for the hot loop.

! ljf.f90: a Fortran routine to call from Python
subroutine lj_energy(x, n, box, rc, energy)
   implicit none
   integer, intent(in) :: n
   double precision, intent(in) :: x(3, n), box, rc
   double precision, intent(out) :: energy
   integer :: i, j
   double precision :: d(3), r2, ir6
   energy = 0.0d0
   do i = 1, n - 1
      do j = i + 1, n
         d = x(:, i) - x(:, j)
         d = d - box * anint(d / box)
         r2 = sum(d * d)
         if (r2 < rc * rc) then
            ir6 = 1.0d0 / r2**3
            energy = energy + 4.0d0 * ir6 * (ir6 - 1.0d0)
         end if
      end do
   end do
end subroutine lj_energy

pip install meson ninja                 # f2py's build backend on Python >= 3.12
python -m numpy.f2py -c ljf.f90 -m ljf  # builds ljf.cpython-*.so

import numpy as np
import ljf                                    # the extension module built by f2py

# 216 atoms on a simple cubic lattice (spacing 1.12) in a periodic box
g = np.arange(6) * 1.12
x = np.array(np.meshgrid(g, g, g, indexing="ij")).reshape(3, -1)   # shape (3, n), as Fortran expects
box, rc = 6 * 1.12, 2.5

e = ljf.lj_energy(np.asfortranarray(x), box, rc)   # n is inferred from the array shape

d = x[:, :, None] - x[:, None, :]                  # the same sum in NumPy, for comparison
d -= box * np.round(d / box)
r2 = (d ** 2).sum(0)[np.triu_indices(x.shape[1], 1)]
r2 = r2[r2 < rc ** 2]
print(f"Fortran: {e:.8f}   NumPy: {np.sum(4 * (r2 ** -6 - r2 ** -3)):.8f}")
print(ljf.lj_energy.__doc__.splitlines()[0])

Fortran: -1039.38039581   NumPy: -1039.38039581
energy = lj_energy(x,box,rc,[n])

f2py inferred that n can be taken from the shape of x, so the Python call has one fewer argument; the docstring line shows the generated signature.

References and further reading