aboutsummaryrefslogtreecommitdiff
path: root/schroedinger/schrodinger_solve.py
blob: 2eb83e862ada94b3f9c9da83bf01264ddf6ac790 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import sys
import argparse
from pathlib import Path

import numpy as np
from numpy.typing import NDArray

from schroedinger import (
    Config, potential_interp, build_potential, solve_schroedinger
)

DESCRIPTION='Solve time independent Schrödinger\'s equation for a given system.'

def save_wavefuncs(
        filename: Path,
        x: NDArray[np.float64],
        v: NDArray[np.float64]
) -> None:
    wavefuncs = np.zeros((x.shape[0], v.shape[1] + 1))
    wavefuncs[:, 0] = x
    for i in range(v.shape[1]):
        wavefuncs[:, 1 + i] = v[:, i]
    np.savetxt(filename, wavefuncs)


def save_expvalues(
        filename: Path,
        x: NDArray[np.float64],
        v: NDArray[np.float64]
) -> None:
    n = v.shape[1]
    delta = np.abs(x[1] - x[0])
    expvalues = np.zeros((n, 2))
    for i in range(n):
        exp_x = delta * np.sum(v[:, i] * x * v[:, i])
        exp_xsq = delta * np.sum(v[:, i] * x**2 * v[:, i])
        expvalues[i, 0] = exp_x
        expvalues[i, 1] = np.sqrt(exp_xsq - exp_x ** 2)
    np.savetxt(filename, expvalues)


def main() -> None:
    parser = argparse.ArgumentParser(
        prog='schrodinger_solve',
        description=DESCRIPTION,
        epilog='')

    parser.add_argument('filename',
                        help='File describing the system to solve')
    parser.add_argument('-o', '--output-dir',
                        help='Output directory for the results')

    args = parser.parse_args()

    conf = Config(args.filename)
    output_path = Path(args.output_dir) if args.output_dir else Path.cwd()

    potential, delta = build_potential(conf)

    e, v = solve_schroedinger(conf.mass, potential[:, 1], delta,
                              conf.eig_interval)

    try:
        np.savetxt(output_path / 'potential.dat', potential)
        np.savetxt(output_path / 'energies.dat', e)
        save_wavefuncs(output_path / 'wavefuncs.dat', potential[:, 0], v)
        save_expvalues(output_path / 'expvalues.dat', potential[:, 0], v)
    except FileNotFoundError:
        print('Output files could not be saved.'
              ' Are you sure the output directory exists?')
        sys.exit(1)


if __name__ == '__main__':
    main()