From 16d2ec58cbc1902216a374d9076a96e361a41030 Mon Sep 17 00:00:00 2001 From: Tomas Jira Date: Fri, 16 Feb 2024 15:42:26 +0100 Subject: [PATCH] small refactoring of educational scripts --- .gitignore | 1 + README.md | 2 +- TODO.md | 2 +- .../python/adiabatic_quantum_dynamics.py | 69 ----------- education/python/eqadet.py | 116 ++++++++++++++++++ education/python/resmet.py | 30 +++-- 6 files changed, 139 insertions(+), 81 deletions(-) delete mode 100644 education/python/adiabatic_quantum_dynamics.py create mode 100644 education/python/eqadet.py diff --git a/.gitignore b/.gitignore index 1254ce3..4fb53f5 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ trajectory.xyz molecule.xyz # exts +*.clean *.lock *.dat *.gif diff --git a/README.md b/README.md index 7d75a58..d8a9282 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,7 @@ Quantum Acorn, a dynamic collection of electronic structure methods, effortlessl ## Features -Below are all the important features of Acorn divided into categories. +Below are all the important features of Acorn divided into categories. If you are looking for the educational scripts, you can find them in the education folder. ### Quantum Mechanical Methods diff --git a/TODO.md b/TODO.md index 558b7db..65a89e5 100644 --- a/TODO.md +++ b/TODO.md @@ -2,7 +2,7 @@ - [ ] Add BAGEL interface with excited state dynamics - [ ] Add excitation specifications to FCI -- [ ] Create a compare script that compares results with ORCA +- [ ] Create a compare script that compares results with ORCA or Psi4 - [x] Add Mulliken analysis - [x] Add some tests - [x] Add the ability to export the atomic integrals diff --git a/education/python/adiabatic_quantum_dynamics.py b/education/python/adiabatic_quantum_dynamics.py deleted file mode 100644 index e3dfdd5..0000000 --- a/education/python/adiabatic_quantum_dynamics.py +++ /dev/null @@ -1,69 +0,0 @@ -#!/usr/bin/env python - -# import plotting and animating functions -import matplotlib.animation as anm -import matplotlib.pyplot as plt - -# import numpy library -import numpy as np - -# define constants and variables -xhalfrange, xpoints, dt, steps, thresh = 16, 1024, 0.1, 1000, 1e-12 -dx, nstates = 2 * xhalfrange / (xpoints - 1), 3 -legend, plotlim = False, 1e-8 - -# define x and k space -x = np.linspace(-xhalfrange, xhalfrange, xpoints) -k = 2 * np.pi * np.fft.fftfreq(len(x), dx) - -# define initial wavefunction and potential -psi0, V = np.exp(-(x - 0.5)**2), 0.5*x**2 - -def energy(wfn): - Ek = 0.5 * np.conj(wfn) * np.fft.ifft(k**2 * np.fft.fft(wfn)) - Er = np.conj(wfn) * V * wfn; return np.sum(Ek + Er).real * dx - -if __name__ == "__main__": - # define R and K operators - R, K, states = np.exp(-0.5 * V * dt), np.exp(-0.5 * k**2 * dt), [] - - for i in range(nstates): - # define initial wavefunction - psi = [0j + psi0 / np.sqrt(np.sum(np.abs(psi0)**2) * dx)] - - for _ in range(steps): - # apply R and K operators and append to wavefunction list - psi.append(R * np.fft.ifft(K * np.fft.fft(R * psi[-1]))); - - # apply Gram-Schmidt orthogonalization - for j in range(i): - psi[-1] -= np.sum(np.conj(states[j][-1]) * psi[-1]) * dx * states[j][-1] - - # normalize wavefunction - psi[-1] /= np.sqrt(np.sum(np.abs(psi[-1])**2) * dx) - - # break if wavefunction has converged - if np.sum(np.abs(psi[-1] - psi[-2])**2) < thresh or np.abs(energy(psi[-1]) - energy(psi[-2])) < thresh: break - - # append wavefunction to list of states and print energy - states.append(psi); print("E_{}:".format(i), energy(psi[-1])) - - # define probability density - D = [[energy(psi) + np.abs(psi)**2 for psi in state] for state in states] - - # define minimum and maximum x values for plotting - xmin = np.min([np.min([np.min(x[np.abs(psi)**2 > plotlim]) for psi in Si]) for Si in states]) - xmax = np.max([np.max([np.max(x[np.abs(psi)**2 > plotlim]) for psi in Si]) for Si in states]) - [fig, ax], xrange = plt.subplots(), np.max([np.abs(xmin), np.abs(xmax)]); plt.tight_layout() - - # set limits and plot initial wavefunction and potential - ax.set_xlim(-xrange, xrange); ax.set_ylim(0, 1.3 * np.max(D[-1][-1])) - ax.plot(x, V); psiplots = [ax.plot(x, Di[0])[0] for Di in D] - - def update(j): - labels = ["$\\psi_{{{},{}}}$".format(i, j if j < len(D[i]) else len(D[i]) - 1) for i in range(len(psiplots))] - for i in range(len(psiplots)): psiplots[i].set_ydata(D[i][j if j < len(D[i]) else -1]) - if legend: ax.legend(psiplots, labels, loc="upper right") - - # animate the wavefunction - ani = anm.FuncAnimation(fig, update, frames=np.max([len(state) for state in states]), interval=30); plt.show(); plt.close("all") diff --git a/education/python/eqadet.py b/education/python/eqadet.py new file mode 100644 index 0000000..c048914 --- /dev/null +++ b/education/python/eqadet.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python + +import argparse as ap, matplotlib.animation as anm, matplotlib.pyplot as plt, numpy as np + +def energy(wfn): + Ek = 0.5 * np.conj(wfn) * np.fft.ifft(k**2 * np.fft.fft(wfn)) + Er = np.conj(wfn) * V * wfn; return np.sum(Ek + Er).real * dx + +if __name__ == "__main__": + # SECTION FOR PARSING COMMAND LINE ARGUMENTS ======================================================================================================================================================= + + # create the parser + parser = ap.ArgumentParser( + prog="eqdet.py", description="Exact Quantum Adiabatic Dynamics Educational Toolkit", add_help=False, + formatter_class=lambda prog: ap.HelpFormatter(prog, max_help_position=128) + ) + + # add the arguments + parser.add_argument("-h", "--help", help="Print this help message.", action="store_true") + parser.add_argument("-i", "--iters", help="Maximum number of iterations. (default: %(default)s)", type=int, default=1000) + parser.add_argument("-n", "--nstate", help="Number of states to consider. (default: %(default)s)", type=int, default=3) + parser.add_argument("-p", "--points", help="Number of discretization points. (default: %(default)s)", type=int, default=1024) + parser.add_argument("-r", "--range", help="Range for the calculated wavefunction. (default: %(default)s)", nargs=2, type=float, default=[-16, 16]) + parser.add_argument("-s", "--tstep", help="Time step of the propagation. (default: %(default)s)", type=float, default=0.1) + parser.add_argument("-t", "--threshold", help="Convergence threshold for the wavefunction. (default: %(default)s)", type=float, default=1e-12) + + # add the flags + parser.add_argument("--real", help="Perform the real time propagation.", action="store_true") + + # parse the arguments + args = parser.parse_args() + + # print the help message if the flag is set + if args.help: parser.print_help(); exit() + + # VARIABLE INITIALIZATION ========================================================================================================================================================================== + + # define the discretization step and state array + dx, states = (args.range[1] - args.range[0]) / (args.points - 1), list() + + # define x and k space + x, k = np.linspace(args.range[0], args.range[1], args.points), 2 * np.pi * np.fft.fftfreq(args.points, dx) + + # define initial wavefunction and potential + psi0, V = np.exp(-(x - 0.5)**2), 0.5 * x**2 + + # define R and K operators for real and imaginary time propagation + Rr, Kr = np.exp(-0.5j * V * args.tstep), np.exp(-0.5j * k**2 * args.tstep) + Ri, Ki = np.exp(-0.5 * V * args.tstep), np.exp(-0.5 * k**2 * args.tstep) + + # WAVEFUNCTION PROPAGATION ========================================================================================================================================================================= + + # propagate each state + for i in range(args.nstate): + # define initial wavefunction + psi = [0j + psi0 / np.sqrt(np.sum(np.abs(psi0)**2) * dx)] + + # loop over imaginary and real time + for R, K in zip([Ri, Rr], [Ki, Kr]): + # break if real time propagation is not required + if not args.real and np.array_equal(R, Rr) and np.array_equal(K, Kr): break + + # reset the wavefunction guess for real time propagation + if np.array_equal(R, Rr) and np.array_equal(K, Kr): psi = [psi[-1]] + + # propagate the wavefunction + for _ in range(args.iters): + # apply R and K operators and append to wavefunction list + psi.append(R * np.fft.ifft(K * np.fft.fft(R * psi[-1]))); + + # apply Gram-Schmidt orthogonalization + for j in (j for j in range(i)): + psi[-1] -= np.sum(np.conj(states[j][-1]) * psi[-1]) * states[j][-1] * dx + + # normalize wavefunction + psi[-1] /= np.sqrt(np.sum(np.abs(psi[-1])**2) * dx) + + # break if wavefunction has converged + if np.sum(np.abs(psi[-1] - psi[-2])**2) < args.threshold or np.abs(energy(psi[-1]) - energy(psi[-2])) < args.threshold: break + + # append wavefunction to list of states and print energy + states.append(psi); print("E_{}:".format(i), energy(psi[-1])) + + # RESULTS AND PLOTTING ============================================================================================================================================================================= + + # define probability density + D = [[energy(psi) + np.abs(psi)**2 for psi in state] for state in states] + + # create the figure and definte tight layout + [fig, ax] = plt.subplots(); plt.tight_layout() + + # define minimum and maximum x values for plotting + xmin = np.min([np.min([np.min(x[np.abs(psi)**2 > 1e-8]) for psi in Si]) for Si in states]) + xmax = np.max([np.max([np.max(x[np.abs(psi)**2 > 1e-8]) for psi in Si]) for Si in states]) + + # define maximum y values for plotting + ymaxreal = max([max([energy(psi) + np.real(psi).max() for psi in state]) for state in states]) + ymaximag = max([max([energy(psi) + np.imag(psi).max() for psi in state]) for state in states]) + + # define the minimum y values for plotting + yminreal = min([min([energy(psi) + np.real(psi).min() for psi in state]) for state in states]) + yminimag = min([min([energy(psi) + np.imag(psi).min() for psi in state]) for state in states]) + + # set limits of the plot + ax.set_xlim(xmin, xmax); ax.set_ylim(np.block([V, yminreal, yminimag]).min(), max([ymaxreal, ymaximag])) + + # plot the potential and initial wavefunctions + ax.plot(x, V); plots = [[ax.plot(x, np.real(state[0]))[0], ax.plot(x, np.imag(state[0]))[0]] for state in states] + + # animation update function + def update(j): + for i in range(len(plots)): plots[i][0].set_ydata(energy(states[i][j if j < len(states[i]) else -1]) + np.real(states[i][j if j < len(states[i]) else -1])) + for i in range(len(plots)): plots[i][1].set_ydata(energy(states[i][j if j < len(states[i]) else -1]) + np.imag(states[i][j if j < len(states[i]) else -1])) + + # animate the wavefunction + ani = anm.FuncAnimation(fig, update, frames=np.max([len(state) for state in states]), interval=30); plt.show(); plt.close("all") # type: ignore diff --git a/education/python/resmet.py b/education/python/resmet.py index ee632c8..42f372c 100644 --- a/education/python/resmet.py +++ b/education/python/resmet.py @@ -5,16 +5,26 @@ A2BOHR = 1.8897259886 ATOM = { - "H": 1, - "O": 8, - "F": 9, + "H" : 1, "He": 2, + "Li": 3, "Be": 4, "B" : 5, "C" : 6, "N" : 7, "O" : 8, "F" : 9, "Ne": 10, + "Na": 11, "Mg": 12, "Al": 13, "Si": 14, "P" : 15, "S" : 16, "Cl": 17, "Ar": 18, + "K" : 19, "Ca": 20, "Sc": 21, "Ti": 22, "V" : 23, "Cr": 24, "Mn": 25, "Fe": 26, "Co": 27, "Ni": 28, "Cu": 29, "Zn": 30, "Ga": 31, "Ge": 32, "As": 33, "Se": 34, "Br": 35, "Kr": 36, + "Rb": 37, "Sr": 38, "Y" : 39, "Zr": 40, "Nb": 41, "Mo": 42, "Tc": 43, "Ru": 44, "Rh": 45, "Pd": 46, "Ag": 47, "Cd": 48, "In": 49, "Sn": 50, "Sb": 51, "Te": 52, "I" : 53, "Xe": 54, + "Cs": 55, "Ba": 56, "La": 57, "Hf": 72, "Ta": 73, "W" : 74, "Re": 75, "Os": 76, "Ir": 77, "Pt": 78, "Au": 79, "Hg": 80, "Tl": 81, "Pb": 82, "Bi": 83, "Po": 84, "At": 85, "Rn": 86, + "Fr": 87, "Ra": 88, "Ac": 89, "Rf":104, "Db":105, "Sg":106, "Bh":107, "Hs":108, "Mt":109, "Ds":110, "Rg":111, "Cn":112, "Nh":113, "Fl":114, "Mc":115, "Lv":116, "Ts":117, "Og":118, + + "Ce": 58, "Pr": 59, "Nd": 60, "Pm": 61, "Sm": 62, "Eu": 63, "Gd": 64, "Tb": 65, "Dy": 66, "Ho": 67, "Er": 68, "Tm": 69, "Yb": 70, "Lu": 71, + "Th": 90, "Pa": 91, "U" : 92, "Np": 93, "Pu": 94, "Am": 95, "Cm": 96, "Bk": 97, "Cf": 98, "Es": 99, "Fm":100, "Md":101, "No":102, "Lr":103, } if __name__ == "__main__": # SECTION FOR PARSING COMMAND LINE ARGUMENTS ======================================================================================================================================================= # create the parser - parser = ap.ArgumentParser(prog="resmet.py", description="Restricted Electronic Structure Methods Educational Toolkit", add_help=False, formatter_class=lambda prog: ap.HelpFormatter(prog, max_help_position=128)) + parser = ap.ArgumentParser( + prog="resmet.py", description="Restricted Electronic Structure Methods Educational Toolkit", + add_help=False, formatter_class=lambda prog: ap.HelpFormatter(prog, max_help_position=128) + ) # add the arguments parser.add_argument("-m", "--molecule", help="Molecule file in the .xyz format. (default: %(default)s)", type=str, default="molecule.xyz") @@ -35,7 +45,7 @@ # print the help message if the flag is set if args.help: parser.print_help(); exit() - # OBTAIN THE MOLECULE ANS ATOMIC INTEGRALS ========================================================================================================================================================= + # OBTAIN THE MOLECULE AND ATOMIC INTEGRALS ========================================================================================================================================================= # get the atomic numbers and coordinates of all atoms atoms = np.array([ATOM[line.split()[0]] for line in open(args.molecule).readlines()[2:]], dtype=int) @@ -77,7 +87,7 @@ VNN += 0.5 * atoms[i] * atoms[j] / np.linalg.norm(coords[i, :] - coords[j, :]) / A2BOHR if i != j else 0 # print the results - print("RHF ENERGY:", E_HF + VNN) + print("RHF ENERGY: {:.8f}".format(E_HF + VNN) + (" (PSI4: {:.8f})".format(psi4.energy("scf/{}".format(args.psi))) if args.psi else "")) # type: ignore # MOLLER-PLESSET PERTRUBATION THEORY OF 2ND ORDER ================================================================================================================================================== if args.mp2: @@ -90,7 +100,7 @@ E_MP2 += Jmo[i, a, j, b] * (2 * Jmo[i, a, j, b] - Jmo[i, b, j, a]) / (eps[i] + eps[j] - eps[a] - eps[b]); # print the results - print("MP2 ENERGY:", E_HF + E_MP2 + VNN) + print("MP2 ENERGY: {:.8f}".format(E_HF + E_MP2 + VNN) + (" (PSI4: {:.8f})".format(psi4.energy("mp2/{}".format(args.psi))) if args.psi else "")) # type: ignore # FULL CONFIGUIRATION INTERACTION ================================================================================================================================================================== if args.fci: @@ -113,7 +123,7 @@ # define the CI Hamiltonian Hci = np.zeros([len(dets), len(dets)]) - # define two-electron part of the Slater-Condon rules, so is array of unique and common spinorbitals [unique, common] + # define the Slater-Condon rules, "so" is an array of unique and common spinorbitals [unique, common] slater0 = lambda so: sum([Hms[m, m] for m in so]) + sum([0.5 * (Jms[m, m, n, n] - Jms[m, n, n, m]) for m, n in it.product(so, so)]) slater1 = lambda so: Hms[so[0], so[1]] + sum([Jms[so[0], so[1], m, m] - Jms[so[0], m, m, so[1]] for m in so[2:]]) slater2 = lambda so: Jms[so[0], so[2], so[1], so[3]] - Jms[so[0], so[3], so[1], so[2]] @@ -121,7 +131,7 @@ # fill the CI Hamiltonian for i in range(Hci.shape[0]): for j in range(Hci.shape[1]): - # copy the determinants so they are not modified + # copy the determinant and define sign aligned, sign = dets[i].copy(), 1 # align the first determinant to the second and calculate the sign @@ -145,4 +155,4 @@ eci, Cci = np.linalg.eigh(Hci); E_FCI = eci[0] - E_HF # print the results - print("FCI ENERGY:", E_HF + E_FCI + VNN) + print("FCI ENERGY: {:.8f}".format(E_HF + E_FCI + VNN) + (" (PSI4: {:.8f})".format(psi4.energy("fci/{}".format(args.psi))) if args.psi else "")) # type: ignore