????

I’m currently working on an energy calculation for an atomistic structure (Graphene) and would appreciate some help in verifying my calculations. Specifically, I want to ensure that certain constants used in my code remain consistent throughout different computations.

To simplify my explanation:

I’ve defined a constant K = 63.9 which is integral to the energy calculation within the function calculate_energy(distances, bonds).

The distances between atoms are computed using several functions that plot the atoms and determine the distances between them. In the stable state, the energy should ideally be 0.

As part of testing my code, I’ve created a plot showing the relationship between energy growth and lattice constant growth (a). This plot helps me verify if I consistently get the same “K” constant (where the coefficient of “K” is the slope of the energy versus lattice constant plot) for each energy calculation, even when expanding the unit cell.

THE PROBLEM:
Despite my efforts, I’m encountering a multiplication factor of approximately 1.16667 on the K constant. Instead of the expected 63.9, I’m getting around 74.5520. I’ve thoroughly reviewed my script but cannot pinpoint the source of this factorization.

Additional Information:

The energy calculation involves calculating bond lengths and angles within the graphene lattice.
The unit cell is expanded by varying the lattice constant a and observing its effect on the total energy.
I’m particularly interested in understanding why there’s a deviation from the expected K constant and how to rectify it.
Any insights or suggestions on where I might look to resolve this issue would be greatly appreciated! Thank you in advance for your help.

SCRIPT SNIPPET:

import numpy as np
import matplotlib.pyplot as plt

# Constants
K = 63.9 # Force constant for bond elongation
xi = 1.56  # Unknown constant (not used in this code yet)
ideal_angle = 120  # Ideal bond angle in degrees
angle_constant = 1.0  # Constant for angle deviation energy contribution

# Define additional basis vectors
origin = np.array([0, 0])

# Function to plot lattice vectors
def plot_vectors(vectors, origin, color):
    for vector in vectors:
        plt.quiver(*origin, *vector, color=color, angles='xy', scale_units='xy', scale=1, width=0.003)

# Function to plot rectangle unit cell
def plot_rectangle(origin, V1, d2, color):
    rect_vertices = [
        origin,
        origin + V1,
        origin + 3 * d2,
        origin + V1 + 3 * d2
    ]
    rect_vertices = np.array(rect_vertices)
    rect_indices = [0, 2, 3, 1, 0]
    for i in range(len(rect_indices) - 1):
        plt.plot([rect_vertices[rect_indices[i]][0], rect_vertices[rect_indices[i + 1]][0]],
                 [rect_vertices[rect_indices[i]][1], rect_vertices[rect_indices[i + 1]][1]], color, linewidth=2)

# Function to plot atoms
def plot_atoms(offset, V1, V2, d1, d2, color1, color2):
    # Red atoms
    plt.plot(offset[0] + d1[0], offset[1] + d1[1], color1, markersize=10)
    plt.plot(offset[0] + d1[0] + V1[0], offset[1] + d1[1] + V1[1], color1, markersize=10)
    plt.plot(offset[0] + d1[0] + V2[0], offset[1] + d1[1] + V2[1], color1, markersize=10)
    plt.plot(offset[0] + d1[0] + V1[0] + V2[0], offset[1] + d1[1] + V1[1] + V2[1], color1, markersize=10)
    plt.plot(offset[0] + a * np.cos(np.radians(60)), offset[1] + a * np.sin(np.radians(60)), color1, markersize=10)

    # Green atoms
    plt.plot(offset[0] + d2[0], offset[1] + d2[1], color2, markersize=10)
    plt.plot(offset[0] + d2[0] + V1[0], offset[1] + d2[1] + V1[1], color2, markersize=10)
    plt.plot(offset[0] + a * np.cos(np.radians(60)), offset[1] + a * np.sin(np.radians(60)) + d2[1], color2, markersize=10)

# Function to calculate angles between bonds
def calculate_angles(bonds):
    angles = []
    for i in range(len(bonds)):
        for j in range(i + 1, len(bonds)):
            common_atom = None
            if np.array_equal(bonds[i][0], bonds[j][0]):
                common_atom = bonds[i][0]
                vec1 = bonds[i][1] - common_atom
                vec2 = bonds[j][1] - common_atom
            elif np.array_equal(bonds[i][0], bonds[j][1]):
                common_atom = bonds[i][0]
                vec1 = bonds[i][1] - common_atom
                vec2 = bonds[j][0] - common_atom
            elif np.array_equal(bonds[i][1], bonds[j][0]):
                common_atom = bonds[i][1]
                vec1 = bonds[i][0] - common_atom
                vec2 = bonds[j][1] - common_atom
            elif np.array_equal(bonds[i][1], bonds[j][1]):
                common_atom = bonds[i][1]
                vec1 = bonds[i][0] - common_atom
                vec2 = bonds[j][0] - common_atom

            if common_atom is not None:
                cos_angle = np.dot(vec1, vec2) / (np.linalg.norm(vec1) * np.linalg.norm(vec2))
                cos_angle = np.clip(cos_angle, -1, 1)  # Ensure cos_angle is within the valid range
                angle = np.arccos(cos_angle) * (180 / np.pi)
                # Validate angle deviation
                if np.abs(angle) > 0.1:  # Include only realistic angles
                    angles.append(angle)
    return angles

# Function to plot carbon-carbon bonds and calculate distances
def plot_carbon_bonds(basis_atoms_red, basis_atoms_green):
    distances = []
    bonds = []  # Store the bonds to calculate angles later
    for red_atom in basis_atoms_red:
        for green_atom in basis_atoms_green:
            red_atom_np = np.array(red_atom)
            green_atom_np = np.array(green_atom)
            distance = np.linalg.norm(red_atom_np - green_atom_np)
            if distance < a + 0.1:
                distances.append(distance)
                bonds.append((red_atom_np, green_atom_np))
                plt.plot([red_atom[0], green_atom[0]], [red_atom[1], green_atom[1]], 'k-')

    return distances, bonds

def generate_basis_atoms(offset, V1, V2, d1, d2):
    basis_atoms_red = [
        offset + origin,
        offset + V1,
        offset + 3 * d2,
        offset + V1 + 3 * d2,
        offset + np.array([a * np.cos(np.radians(60)), a * np.sin(np.radians(60))])
    ]

    basis_atoms_green = [
        offset + d2,
        offset + d2 + V1,
        offset + np.array([a * np.cos(np.radians(60)), a * np.sin(np.radians(60)) + d2[1]])
    ]

    return basis_atoms_red, basis_atoms_green

# Function to plot atoms and lattice points
def plot_atoms_and_lattice(supercell_counter, vertical_extension_counter, V1, V2, d1, d2):
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid(True)
    plt.axis('equal')

    all_red_atoms = set()
    all_green_atoms = set()

    # Generate and plot lattice points for the supercell
    for i in range(-supercell_counter, supercell_counter + 1):
        for j in range(-vertical_extension_counter, vertical_extension_counter + 1):
            offset = np.array([i * V1[0] + j * V2[0], i * V1[1] + j * V2[1]])
            plt.plot(offset[0], offset[1], 'ko')  # Lattice points
            plot_atoms(offset, V1, V2, d1, d2, 'ro', 'go')
            red_atoms, green_atoms = generate_basis_atoms(offset, V1, V2, d1, d2)
            all_red_atoms.update(map(tuple, red_atoms))
            all_green_atoms.update(map(tuple, green_atoms))

    # Plot the central unit cell
    plot_rectangle(origin, V1, d2, 'blue')
    plot_atoms(origin, V1, V2, d1, d2, 'ro', 'go')

    # Plot the carbon-carbon bonds considering the entire supercell
    distances, bonds = plot_carbon_bonds(list(all_red_atoms), list(all_green_atoms))

    return distances, bonds, list(all_red_atoms), list(all_green_atoms)

# Function to calculate total energy
def calculate_energy(distances, bonds):
    total_energy = 0
    eq_distance = 0.824  # Equilibrium distance for bond length
    for distance in distances:
        total_energy += 0.5 * K * round((distance - eq_distance), 3) ** 2  # Energy contribution from bond elongation

    angles = calculate_angles(bonds)
    for angle in angles:
         angle_deviation = angle - ideal_angle
         total_energy += 0.5 * angle_constant * angle_deviation ** 2  # Energy contribution from angle deviation

    return total_energy

# Function to dynamically update lattice vectors and dependent vectors
def update_lattice_vectors(a):
    d2 = np.array([0, (1 / 3 ** 0.5) * a])
    V1 = np.array([a, 0])
    V2 = np.array([0, 3 * d2[1]])
    d1 = np.array([0, 0])
    return V1, V2, d1, d2

# Example usage: plot with horizontal supercell counter and vertical extension counter
horizontal_supercell_counter = 0
vertical_extension_counter = 0
a = 1.428  # Initial lattice constant
V1, V2, d1, d2 = update_lattice_vectors(a)
distances, bonds, all_red_atoms, all_green_atoms = plot_atoms_and_lattice(horizontal_supercell_counter, vertical_extension_counter, V1, V2, d1, d2)
plot_vectors([3 * d2, V1], [0, 0], 'blue')

# Generate lattice points for the central unit cell
lattice_points = []
for i in range(-horizontal_supercell_counter - 1, horizontal_supercell_counter + 3):
    for j in range(-vertical_extension_counter - 1, vertical_extension_counter + 3):
        lattice_points.append((i * V1[0] + j * V2[0], i * V1[1] + j * V2[1]))

for point in lattice_points:
    plt.plot(point[0], point[1], 'ko')

angles = calculate_angles(bonds)
energy = calculate_energy(distances, bonds)
plt.title(f'Graphene Lattice - Supercell VisualizationnTotal Energy: {energy:.4f} (only bond elongation for now)')
print("Distances:", len(distances))
print("Angles:", len(angles))

# Print coordinates of all atoms
print("Atom Coordinates (X, Y):")

all_atom_coords = list(all_red_atoms) + list(all_green_atoms)
all_atom_coords = [list(coord) for coord in all_atom_coords]

print("Number of atoms:", len(all_atom_coords))
print("Atom Coordinates:", all_atom_coords)
plt.show()


# Calculation and plotting of energy versus lattice constant
energy_values = []
eq_distance_values = []
a_values = np.arange(2, 10, 0.1)
for a in a_values:
    a = np.round(a,3)
    print(f"Lattice constant: {a}")
    V1, V2, d1, d2 = update_lattice_vectors(a)
    distances, bonds, all_red_atoms, all_green_atoms = plot_atoms_and_lattice(horizontal_supercell_counter, vertical_extension_counter, V1, V2, d1, d2)
    energy = calculate_energy(distances, bonds)
    print(energy)
    energy_values.append(energy)
    eq_distance_values.append(np.mean(distances))  # Calculate the average distance for equilibrium distance

# Plotting the energy versus lattice constant
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(a_values, energy_values, 'ro', label='Calculated Energies')
p_fit = np.poly1d(np.polyfit(a_values, energy_values, 2))
ax.plot(a_values, p_fit(a_values), label='Quadratic Fit')
K_squared_constant = np.round((np.polyfit(a_values, energy_values, 2)[0]), 3)
ax.set_title(f'Total Energy vs Lattice Constant nQuadratic Fit: Coefficient of K^2 = {K_squared_constant:.4f} [eV*A^-2]')
ax.set_xlabel('Lattice Constant (a)')
ax.set_ylabel('Total Energy (eV)')
ax.legend()
ax.grid(True)

plt.show()

I attempted to simulate energy calculations for a graphene structure using a predefined constant (K = 63.9) in my code. I expected the energy calculations to yield consistent results with K remaining unchanged. However, upon plotting Energy versus Lattice Constant (a), I found that the calculated constant K deviated unexpectedly, showing a multiplication factor of approximately 1.1667 (resulting in K = 74.5520). This discrepancy is puzzling as I can’t identify any source in my script that accounts for this factorization.