Skip to content

Charged Systems: Ions, Zwitterions & Electrochemistry

What You'll Learn

  • How to set molecular charge correctly for single molecules and batches
  • Working with anions, cations, and zwitterions
  • Dipole moment caveats for charged species
  • Long-range electrostatics for periodic charged systems
  • Understanding charge prediction behavior in molecular dynamics

Prerequisites

Relevant Models

Charge handling is available across the standard AIMNet2 model families (aimnet2, aimnet2-2025, aimnet2-nse, aimnet2-pd, aimnet2-b973c). The aimnet2-rxn family is the exception: it supports only net-neutral systems, including zwitterions, and rejects nonzero total charge. Choose the model based on both chemistry and charge state.

Setting the Charge

Direct Calculator

With AIMNet2Calculator, the charge is passed in the input dictionary as a scalar or tensor:

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")

# Neutral molecule
result = calc({
    "coord": coords,
    "numbers": numbers,
    "charge": torch.tensor(0.0),
})

# Anion (e.g., hydroxide OH-)
result = calc({
    "coord": coords_oh,
    "numbers": numbers_oh,
    "charge": torch.tensor(-1.0),
})

# Cation (e.g., ammonium NH4+)
result = calc({
    "coord": coords_nh4,
    "numbers": numbers_nh4,
    "charge": torch.tensor(1.0),
})

Batched Calculations

When computing multiple molecules in a single batch using mol_idx, provide the charge as a 1D tensor with one value per molecule:

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")

# Two molecules: first neutral, second with charge -1
result = calc({
    "coord": combined_coords,    # (N_total, 3)
    "numbers": combined_numbers,  # (N_total,)
    "charge": torch.tensor([0.0, -1.0]),
    "mol_idx": mol_idx,           # (N_total,) mapping atoms to molecules
})

ASE Interface

With AIMNet2ASE, set the charge at construction or update it later:

from aimnet.calculators import AIMNet2ASE, AIMNet2Calculator

# Set charge at construction
base_calc = AIMNet2Calculator("aimnet2")
ase_calc = AIMNet2ASE(base_calc, charge=-1)

# Update charge later
ase_calc.set_charge(-1)

# Or set charge for a cation
ase_calc = AIMNet2ASE("aimnet2", charge=1)

Always Set the Correct Charge

The charge parameter directly affects predicted partial atomic charges and total energy. An incorrect charge will produce incorrect results without raising an error. The model constrains partial charges to sum to the specified total charge.

Worked Examples

Hydroxide Anion (OH-)

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")

# OH- ion
oh_coords = torch.tensor([
    [0.0, 0.0, 0.0],     # O
    [0.0, 0.0, 0.9572],  # H
])
oh_numbers = torch.tensor([8, 1])

result = calc({
    "coord": oh_coords,
    "numbers": oh_numbers,
    "charge": torch.tensor(-1.0),
}, forces=True)

print(f"OH- energy: {result['energy'].item():.4f} eV")
print(f"Partial charges: {result['charges']}")
# Expect: negative charges summing to -1.0
print(f"Charge sum: {result['charges'].sum().item():.4f}")

Protonated Methylamine (CH3NH3+)

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")

# Methylammonium: CH3-NH3+
coords = torch.tensor([
    [ 0.000,  0.000,  0.000],   # C
    [ 0.000,  0.000,  1.520],   # N
    [ 1.030,  0.000, -0.380],   # H (on C)
    [-0.515,  0.893, -0.380],   # H (on C)
    [-0.515, -0.893, -0.380],   # H (on C)
    [ 0.943,  0.000,  1.920],   # H (on N)
    [-0.471,  0.817,  1.920],   # H (on N)
    [-0.471, -0.817,  1.920],   # H (on N)
])
numbers = torch.tensor([6, 7, 1, 1, 1, 1, 1, 1])

result = calc({
    "coord": coords,
    "numbers": numbers,
    "charge": torch.tensor(1.0),
}, forces=True)

print(f"CH3NH3+ energy: {result['energy'].item():.4f} eV")
print(f"Partial charges: {result['charges']}")
print(f"Charge sum: {result['charges'].sum().item():.4f}")
# Should sum to +1.0

Amino Acid Zwitterion (Glycine)

Amino acids exist as zwitterions at physiological pH, with a protonated amine (NH3+) and deprotonated carboxylate (COO-). The overall charge is zero, but the internal charge separation is large.

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")

# Glycine zwitterion: +H3N-CH2-COO-
# Overall neutral, but with internal charge separation
coords = torch.tensor([
    [ 0.000,  0.000,  0.000],   # N (NH3+)
    [ 1.475,  0.000,  0.000],   # C (alpha)
    [ 2.011,  1.420,  0.000],   # C (carboxyl)
    [ 1.349,  2.399,  0.337],   # O
    [ 3.191,  1.557, -0.336],   # O
    [-0.380,  1.000,  0.200],   # H (on N)
    [-0.380, -0.500,  0.850],   # H (on N)
    [-0.380, -0.500, -0.850],   # H (on N)
    [ 1.830, -0.500,  0.920],   # H (on C-alpha)
    [ 1.830, -0.500, -0.920],   # H (on C-alpha)
])
numbers = torch.tensor([7, 6, 6, 8, 8, 1, 1, 1, 1, 1])

result = calc({
    "coord": coords,
    "numbers": numbers,
    "charge": torch.tensor(0.0),  # Overall neutral
}, forces=True)

print(f"Glycine zwitterion energy: {result['energy'].item():.4f} eV")
print(f"Partial charges:")
for i, (z, q) in enumerate(zip(numbers, result['charges'].flatten())):
    elem = {1: "H", 6: "C", 7: "N", 8: "O"}[z.item()]
    print(f"  {elem}({i}): {q.item():+.4f}")
# Expect: N has positive charge, O atoms have negative charges

Dipole Moment Caveats

AIMNet2 computes the dipole moment classically from the predicted partial atomic charges:

mu = sum_i(q_i * r_i)

where q_i is the partial charge on atom i and r_i is its position vector.

Origin Dependence for Charged Species

For a neutral molecule (total charge = 0), the dipole moment is independent of the coordinate origin. This is the standard physical definition.

For a charged species (total charge != 0), the dipole moment depends on the choice of coordinate origin. Shifting all coordinates by a constant vector d changes the dipole by Q_total * d, where Q_total is the total charge. This is a fundamental property of multipole expansions, not a limitation of AIMNet2.

Practical consequence: Dipole moments are only physically meaningful for neutral molecules. For charged species, compare dipoles only if all geometries use the same coordinate convention (e.g., center of mass at the origin).

from ase import Atoms
from aimnet.calculators import AIMNet2ASE, AIMNet2Calculator
import numpy as np

# Example: dipole of neutral water (origin-independent)
base_calc = AIMNet2Calculator("aimnet2")
water = Atoms("OH2", positions=[
    [0.0, 0.0, 0.1173],
    [0.0, 0.7572, -0.4692],
    [0.0, -0.7572, -0.4692],
])
water.calc = AIMNet2ASE(base_calc, charge=0)

dipole = water.calc.get_dipole_moment(water)
print(f"Water dipole: {np.linalg.norm(dipole):.3f} e*A")
# Approximately 0.38 e*A = 1.85 D

Charge Behavior in Molecular Dynamics

Atomic Charges Are Not Conserved Across MD Frames

In AIMNet2, atomic partial charges are predicted independently for each frame during molecular dynamics. The model predicts charges that sum to the specified total charge, but the distribution among atoms can fluctuate from frame to frame.

This means:

  • Charge on any individual atom may change between MD steps
  • Only the total charge is constrained to be constant
  • There is no continuity constraint on per-atom charges between frames
  • This is physically reasonable (charge redistribution happens in real molecules) but differs from fixed-charge force fields

For applications requiring smooth charge evolution (e.g., computing current from charge fluxes), consider post-processing to smooth the charge trajectory.

Ewald and PME for Periodic Charged Systems

Single Molecule in a Periodic Box

For charged molecules in a periodic box, Ewald summation (or PME for larger cells) handles the long-range electrostatics with controllable accuracy. Both methods share the ewald_accuracy parameter (default 1e-6, matching the nvalchemiops default) and use a per-call neighbor list sized to the estimated real-space cutoff:

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")
calc.set_lrcoulomb_method("ewald", ewald_accuracy=1e-6)

cell = torch.tensor([
    [15.0, 0.0, 0.0],
    [0.0, 15.0, 0.0],
    [0.0, 0.0, 15.0],
])

result = calc({
    "coord": coords_ion,
    "numbers": numbers_ion,
    "charge": torch.tensor(-1.0),
    "cell": cell,
}, forces=True, stress=True)

calc.set_lrcoulomb_method("pme", ewald_accuracy=1e-6)
result_pme = calc({
    "coord": coords_ion,
    "numbers": numbers_ion,
    "charge": torch.tensor(-1.0),
    "cell": cell,
}, forces=True, stress=True)

Ewald and PME support force/stress losses (train=True), but not calculator Hessian requests. DSF only supports inference forces/stress. See Long-Range Methods → Derivative Support.

For non-neutral cells, Ewald and PME use the standard uniform-background convention. DSF is cutoff-based and remains finite for non-neutral inputs, but it does not apply the same Ewald/PME background correction.

Automatic Method Switching for PBC

When periodic boundary conditions are detected (a cell is provided), the calculator automatically switches from the "simple" Coulomb method to DSF:

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")
# Default method is "simple"

# Providing a cell triggers automatic switch to DSF
cell = torch.tensor([
    [10.0, 0.0, 0.0],
    [0.0, 10.0, 0.0],
    [0.0, 0.0, 10.0],
])

result = calc({
    "coord": coords,
    "numbers": numbers,
    "charge": torch.tensor(-1.0),
    "cell": cell,
}, forces=True)
# A warning is issued: "Switching to DSF Coulomb for PBC"

To avoid the warning, set the method explicitly before the first periodic calculation:

calc.set_lrcoulomb_method("dsf", cutoff=15.0)

Protonation and Deprotonation Energetics

A common application is computing protonation/deprotonation energies. Since AIMNet2 is gas-phase only, these are gas-phase proton affinities:

import torch
from aimnet.calculators import AIMNet2Calculator

calc = AIMNet2Calculator("aimnet2")
ev_to_kcal = 23.0609

# Example: gas-phase proton affinity of acetate
# CH3COO- + H+ -> CH3COOH

# Acetate anion (CH3COO-)
acetate_coords = torch.tensor([
    [ 0.000,  0.000,  0.000],   # C (methyl)
    [ 1.520,  0.000,  0.000],   # C (carboxyl)
    [ 2.100,  1.090,  0.000],   # O
    [ 2.100, -1.090,  0.000],   # O
    [-0.400,  1.000,  0.200],   # H
    [-0.400, -0.500,  0.900],   # H
    [-0.400, -0.500, -0.900],   # H
])
acetate_numbers = torch.tensor([6, 6, 8, 8, 1, 1, 1])

# Acetic acid (CH3COOH) - add proton to one oxygen
acetic_coords = torch.tensor([
    [ 0.000,  0.000,  0.000],   # C (methyl)
    [ 1.520,  0.000,  0.000],   # C (carboxyl)
    [ 2.100,  1.090,  0.000],   # O (C=O)
    [ 2.100, -1.090,  0.000],   # O (C-OH)
    [-0.400,  1.000,  0.200],   # H
    [-0.400, -0.500,  0.900],   # H
    [-0.400, -0.500, -0.900],   # H
    [ 3.060, -1.090,  0.000],   # H (on O-H)
])
acetic_numbers = torch.tensor([6, 6, 8, 8, 1, 1, 1, 1])

e_acetate = calc({
    "coord": acetate_coords,
    "numbers": acetate_numbers,
    "charge": torch.tensor(-1.0),
})["energy"].item()

e_acetic = calc({
    "coord": acetic_coords,
    "numbers": acetic_numbers,
    "charge": torch.tensor(0.0),
})["energy"].item()

# Gas-phase proton affinity (no proton energy needed for relative comparison)
# For absolute PA, the proton energy is 0 by convention in gas phase
delta_e = (e_acetic - e_acetate) * ev_to_kcal
print(f"Deprotonation energy (gas phase): {-delta_e:.1f} kcal/mol")

Gas-Phase vs Solution

Gas-phase proton affinities differ substantially from solution-phase pKa values. Solvation stabilization of charged species is typically 50-80 kcal/mol for ions in water. Do not directly compare AIMNet2 gas-phase energetics with experimental solution-phase data without solvation corrections.

What's Next