import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from IPython.display import HTML
import json
import random
from rdkit import RDLogger
RDLogger.DisableLog('rdApp.*') # This hides all RDKit warnings
# --- SETUP COLORS ---
COLORS = {
'neutral': {'text_primary': '#333333', 'text_secondary': '#666666', 'border': '#dee2e6', 'background': '#f5f5f5'},
'chart': {'curve': '#9C27B0', 'pka_line': '#FF9800', 'current': '#f44336'}
}
# --- 1. PREPARE MOLECULAR TEMPLATES ---
def prepare_templates():
# A. Aspirin System
aspirin = Chem.MolFromSmiles('CC(=O)Oc1ccccc1C(=O)O')
aspirin = Chem.AddHs(aspirin)
# Identify the acidic proton (H on COOH) and its connected Oxygen
patt = Chem.MolFromSmarts('C(=O)[OH]')
matches = aspirin.GetSubstructMatches(patt)
acidic_h_idx = -1
acidic_o_idx_orig = -1
if matches:
# matches[0] is (Carbon, Dbl_O, Sgl_O)
o_atom_idx = matches[0][2]
o_atom = aspirin.GetAtomWithIdx(o_atom_idx)
acidic_o_idx_orig = o_atom_idx
for neighbor in o_atom.GetNeighbors():
if neighbor.GetSymbol() == 'H':
neighbor.SetIsotope(2)
acidic_h_idx = neighbor.GetIdx()
AllChem.EmbedMolecule(aspirin, randomSeed=42)
AllChem.MMFFOptimizeMolecule(aspirin)
# Create Base (RCOO-) by removing the H
aspirin_base = Chem.RWMol(aspirin)
aspirin_base.RemoveAtom(acidic_h_idx)
# Calculate new O index in the base molecule (indices shift after removal)
asp_base_o_idx = acidic_o_idx_orig
if acidic_h_idx < acidic_o_idx_orig:
asp_base_o_idx -= 1
# Get H position relative to Aspirin
conf_asp = aspirin.GetConformer()
start_pos = np.array(conf_asp.GetAtomPosition(acidic_h_idx))
# B. Water System (Hydronium)
hydronium = Chem.MolFromSmiles('[OH3+]')
hydronium = Chem.AddHs(hydronium)
# Mark the 3rd H as the "accepted" one
target_h_idx = 3
hydronium.GetAtomWithIdx(target_h_idx).SetIsotope(2)
# Get the Oxygen index (usually 0 for [OH3+])
hyd_o_idx_orig = 0
for atom in hydronium.GetAtoms():
if atom.GetSymbol() == 'O':
hyd_o_idx_orig = atom.GetIdx()
break
AllChem.EmbedMolecule(hydronium, randomSeed=42)
AllChem.MMFFOptimizeMolecule(hydronium)
# Create Base (H2O) by removing the H
water_base = Chem.RWMol(hydronium)
water_base.RemoveAtom(target_h_idx)
# Calculate new O index in water base
wat_base_o_idx = hyd_o_idx_orig
if target_h_idx < hyd_o_idx_orig:
wat_base_o_idx -= 1
conf_hyd = hydronium.GetConformer()
end_pos = np.array(conf_hyd.GetAtomPosition(target_h_idx))
return aspirin, aspirin_base, start_pos, asp_base_o_idx, hydronium, water_base, end_pos, wat_base_o_idx
# Initialize global templates
ASPIRIN, ASPIRIN_BASE, ASP_H_POS, ASP_O_IDX, HYDRONIUM, WATER_BASE, HYD_H_POS, WAT_O_IDX = prepare_templates()
def create_proton():
# A single atom molecule for the flying proton
mol = Chem.MolFromSmiles('[2H]') # Initially Isotope 2
AllChem.EmbedMolecule(mol)
return mol
def add_thermal_motion(mol, seed, amplitude=0.15):
"""Add random thermal vibration to molecule"""
conf = mol.GetConformer()
np.random.seed(seed)
for i in range(mol.GetNumAtoms()):
pos = conf.GetAtomPosition(i)
noise = np.random.normal(0, amplitude, 3)
conf.SetAtomPosition(i, [pos.x + noise[0], pos.y + noise[1], pos.z + noise[2]])
def position_molecules_interpolated(fraction_deprot, frame_idx):
"""
Creates a scene where the proton position is interpolated.
Manages bonds and coloring dynamically.
"""
pairs_layout = [
{'asp': (-6, 5, 0), 'wat': (-1.5, 3.5, 0)},
{'asp': (6, 5, 0), 'wat': (1.5, 3.5, 0)},
{'asp': (-6, -5, 0), 'wat': (-1.5, -3.5, 0)},
{'asp': (6, -5, 0), 'wat': (1.5, -3.5, 0)}
]
combined_mol = None
# Store data for post-processing bonds and colors
# List of dicts: {'proton_idx': int, 'asp_o_idx': int, 'wat_o_idx': int, 't': float}
bond_info = []
current_atom_count = 0
val = fraction_deprot * 4.0
# 1. BUILD GEOMETRY
for i, layout in enumerate(pairs_layout):
asp_offset = np.array(layout['asp'])
wat_offset = np.array(layout['wat'])
# Calculate progress 't' for this pair
if val >= i + 1: t = 1.0
elif val <= i: t = 0.0
else: t = val - i
# A. Aspirin Base
mol_asp = Chem.Mol(ASPIRIN_BASE)
conf = mol_asp.GetConformer()
for idx in range(mol_asp.GetNumAtoms()):
p = conf.GetAtomPosition(idx)
conf.SetAtomPosition(idx, [p.x + asp_offset[0], p.y + asp_offset[1], p.z + asp_offset[2]])
add_thermal_motion(mol_asp, seed=frame_idx*1000 + i*100, amplitude=0.1)
if combined_mol is None: combined_mol = mol_asp
else: combined_mol = Chem.CombineMols(combined_mol, mol_asp)
# Track Aspirin Oxygen Global Index
global_asp_o = current_atom_count + ASP_O_IDX
current_atom_count += mol_asp.GetNumAtoms()
# B. Water Base
mol_wat = Chem.Mol(WATER_BASE)
conf = mol_wat.GetConformer()
for idx in range(mol_wat.GetNumAtoms()):
p = conf.GetAtomPosition(idx)
conf.SetAtomPosition(idx, [p.x + wat_offset[0], p.y + wat_offset[1], p.z + wat_offset[2]])
add_thermal_motion(mol_wat, seed=frame_idx*1000 + i*100 + 50, amplitude=0.15)
combined_mol = Chem.CombineMols(combined_mol, mol_wat)
# Track Water Oxygen Global Index
global_wat_o = current_atom_count + WAT_O_IDX
current_atom_count += mol_wat.GetNumAtoms()
# C. Proton
proton = create_proton()
start_g = asp_offset + ASP_H_POS
end_g = wat_offset + HYD_H_POS
curr_pos = start_g * (1 - t) + end_g * t
proton.GetConformer().SetAtomPosition(0, curr_pos)
add_thermal_motion(proton, seed=frame_idx*1000 + i*100 + 99, amplitude=0.1)
combined_mol = Chem.CombineMols(combined_mol, proton)
# Track Proton Global Index
global_prot = current_atom_count # It's the first atom of the single-atom molecule
current_atom_count += 1
bond_info.append({
'proton_idx': global_prot,
'asp_o_idx': global_asp_o,
'wat_o_idx': global_wat_o,
't': t
})
# 2. ADD BONDS & UPDATE ISOTOPES (Requires Editable Mol)
rw_mol = Chem.RWMol(combined_mol)
for info in bond_info:
t = info['t']
p_idx = info['proton_idx']
# --- COLORING LOGIC ---
# 0.05 to 0.95 = Active Transfer (Cyan/Isotope 2)
# Ends = Attached (White/Isotope 0)
atom = rw_mol.GetAtomWithIdx(p_idx)
if 0.05 < t < 0.95:
atom.SetIsotope(2) # Cyan
else:
atom.SetIsotope(0) # White
# --- BONDING LOGIC ---
# t < 0.4: Bonded to Aspirin
# t > 0.6: Bonded to Water
# 0.4 - 0.6: Bond Broken (Gap)
if t < 0.4:
rw_mol.AddBond(p_idx, info['asp_o_idx'], Chem.BondType.SINGLE)
elif t > 0.6:
rw_mol.AddBond(p_idx, info['wat_o_idx'], Chem.BondType.SINGLE)
return Chem.MolToMolBlock(rw_mol)
# --- CALCULATIONS ---
def calculate_fraction_deprotonated(pH, pKa):
return 1 / (1 + 10**(pKa - pH))
def generate_ph_trajectory(pKa=3.5):
trajectory_data = []
ph_values = []
fractions_deprotonated = []
status_messages = []
# 81 frames from pH 0 to 8
for i in range(81):
pH = i * 0.1
fraction = calculate_fraction_deprotonated(pH, pKa)
# Generate interpolated scene
scene_sdf = position_molecules_interpolated(fraction, i)
trajectory_data.append(scene_sdf)
ph_values.append(pH)
fractions_deprotonated.append(fraction * 100)
if abs(pH - pKa) < 0.2:
status_messages.append(f"pH {pH:.1f}: Equilibrium (~50% Transfer)")
elif pH < pKa:
status_messages.append(f"pH {pH:.1f}: Mostly Protonated (H on Aspirin)")
else:
status_messages.append(f"pH {pH:.1f}: Mostly Deprotonated (H on Water)")
return trajectory_data, ph_values, fractions_deprotonated, status_messages, pKa
# --- HTML VIEWER ---
def create_ph_viewer(width=1000):
trajectory_data, ph_values, fractions, status_messages, pKa = generate_ph_trajectory()
full_trajectory = "\n$$$$\n".join(trajectory_data) + "\n$$$$\n"
ids = {k: f"{k}_{np.random.randint(1e5, 9e5)}" for k in
['viewer', 'slider', 'ph', 'fraction', 'status', 'canvas']}
html = f"""
<div style="font-family: Arial, sans-serif; max-width: {width}px; margin: 20px auto;">
<div style="border: 3px solid {COLORS['chart']['curve']}; border-radius: 10px; padding: 20px; background: white;">
<div style="text-align: center; margin-bottom: 20px;">
<h2 style="color: {COLORS['chart']['curve']}; margin: 0;">Interactive Proton Transfer</h2>
<p style="color: #666; font-size: 14px;">Drag the slider to visualize the bond breaking and formation</p>
</div>
<div style="display: flex; gap: 20px;">
<div style="flex: 1.5; position: relative;">
<div id="{ids['viewer']}" style="width: 100%; height: 450px; border: 2px solid #dee2e6; border-radius: 8px; background: #f0f8ff;"></div>
<div style="position: absolute; top: 10px; right: 10px; background: rgba(255,255,255,0.9); padding: 8px; border-radius: 4px; font-size: 12px; border: 1px solid #ccc;">
<div style="margin-bottom: 4px;"><span style="color: cyan; font-weight: bold; font-size: 14px;">●</span> Transferring Proton</div>
<div style="margin-bottom: 4px;"><span style="color: red; font-weight: bold;">●</span> Oxygen</div>
<div><span style="color: grey; font-weight: bold;">●</span> Carbon</div>
</div>
</div>
<div style="flex: 1; display: flex; flex-direction: column;">
<canvas id="{ids['canvas']}" width="350" height="300" style="border: 1px solid #eee; border-radius: 8px; margin-bottom: 15px;"></canvas>
<div style="background: {COLORS['neutral']['background']}; padding: 15px; border-radius: 8px;">
<div style="display: flex; justify-content: space-between; margin-bottom: 10px;">
<div style="text-align: center;">
<div style="font-size: 12px; color: #666;">pH</div>
<div id="{ids['ph']}" style="font-size: 24px; font-weight: bold; color: {COLORS['chart']['curve']};">0.0</div>
</div>
<div style="text-align: center;">
<div style="font-size: 12px; color: #666;">% Deprotonated</div>
<div id="{ids['fraction']}" style="font-size: 24px; font-weight: bold; color: #FF6B6B;">0%</div>
</div>
</div>
<input type="range" id="{ids['slider']}" min="0" max="80" value="0"
style="width: 100%; height: 8px; cursor: pointer; background: #ddd; accent-color: {COLORS['chart']['curve']};">
<div id="{ids['status']}" style="text-align: center; margin-top: 10px; font-size: 13px; font-weight: bold; color: #444;"></div>
</div>
</div>
</div>
</div>
</div>
<script>
(function() {{
setTimeout(function() {{
const container = document.getElementById('{ids['viewer']}');
const config = {{backgroundColor: '#f0f8ff', antialias: true}};
const viewer = $3Dmol.createViewer(container, config);
const trajectory = `{full_trajectory}`;
const phValues = {json.dumps(ph_values)};
const fractions = {json.dumps(fractions)};
const statusMsgs = {json.dumps(status_messages)};
const pKa = {pKa};
viewer.addModelsAsFrames(trajectory, 'sdf');
// Standard Styles
viewer.setStyle({{}}, {{stick: {{radius: 0.15, color: 'lightgrey'}}, sphere: {{scale: 0.3, color: 'lightgrey'}} }});
viewer.setStyle({{elem: 'C'}}, {{stick: {{color: '#909090'}}, sphere: {{scale: 0.3, color: '#909090'}} }});
viewer.setStyle({{elem: 'O'}}, {{stick: {{color: '#FF0000'}}, sphere: {{scale: 0.35, color: '#FF0000'}} }});
// Default H Style (White) - applies to standard H and inactive transfer H
viewer.setStyle({{elem: 'H'}}, {{stick: {{color: '#FFFFFF'}}, sphere: {{scale: 0.25, color: '#FFFFFF'}} }});
// SPECIAL STYLE: ACTIVE PROTON (Isotope 2)
// Only applies when the proton is in flight (mid-transfer)
viewer.setStyle({{elem: 'H', isotope: 2}}, {{
sphere: {{color: 'cyan', scale: 0.5}},
stick: {{color: 'cyan', radius: 0.2}}
}});
viewer.zoomTo();
viewer.render();
const slider = document.getElementById('{ids['slider']}');
const phDisplay = document.getElementById('{ids['ph']}');
const fracDisplay = document.getElementById('{ids['fraction']}');
const statusDisplay = document.getElementById('{ids['status']}');
const canvas = document.getElementById('{ids['canvas']}');
const ctx = canvas.getContext('2d');
let currentFrame = 0;
function drawChart() {{
const w = canvas.width;
const h = canvas.height;
const pad = 40;
ctx.clearRect(0,0,w,h);
// Draw Axes
ctx.beginPath();
ctx.strokeStyle = '#333';
ctx.lineWidth = 2;
ctx.moveTo(pad, pad);
ctx.lineTo(pad, h-pad);
ctx.lineTo(w-pad, h-pad);
ctx.stroke();
// Draw Curve
ctx.beginPath();
ctx.strokeStyle = '{COLORS['chart']['curve']}';
ctx.lineWidth = 3;
for(let i=0; i<fractions.length; i++) {{
let x = pad + (i/(fractions.length-1)) * (w-2*pad);
let y = (h-pad) - (fractions[i]/100) * (h-2*pad);
if(i===0) ctx.moveTo(x,y);
else ctx.lineTo(x,y);
}}
ctx.stroke();
// Draw pKa line
let pKaX = pad + (pKa/8.0) * (w-2*pad);
ctx.beginPath();
ctx.strokeStyle = 'orange';
ctx.setLineDash([5,5]);
ctx.moveTo(pKaX, pad);
ctx.lineTo(pKaX, h-pad);
ctx.stroke();
ctx.setLineDash([]);
// Draw Current Position Dot
let currX = pad + (currentFrame/(fractions.length-1)) * (w-2*pad);
let currY = (h-pad) - (fractions[currentFrame]/100) * (h-2*pad);
ctx.beginPath();
ctx.fillStyle = '{COLORS['chart']['current']}';
ctx.arc(currX, currY, 6, 0, 2*Math.PI);
ctx.fill();
ctx.strokeStyle = 'white';
ctx.stroke();
}}
function update() {{
viewer.setFrame(currentFrame);
viewer.render();
phDisplay.innerText = phValues[currentFrame].toFixed(1);
fracDisplay.innerText = fractions[currentFrame].toFixed(0) + '%';
statusDisplay.innerText = statusMsgs[currentFrame];
drawChart();
}}
slider.addEventListener('input', function() {{
currentFrame = parseInt(this.value);
update();
}});
update();
}}, 100);
}})();
</script>
"""
return HTML(html)
create_ph_viewer()
