app.py

import os
import tempfile
import numpy as np
import gradio as gr
from ase.io import read, write
from ase.io.trajectory import Trajectory
import hashlib

# ==== Usar componente nativo Molecule3D ====
try:
    from gradio_molecule3d import Molecule3D
    HAVE_MOL3D = True
    print("✅ gradio_molecule3d loaded successfully")
except Exception as e:
    HAVE_MOL3D = False
    print(f"❌ gradio_molecule3d not available: {e}")

# ==== Función para convertir trayectoria a archivo temporal PDB ====
def traj_to_molecule3d_file(traj_path):
    """
    Convierte una trayectoria ASE (.traj) a un archivo PDB temporal para Molecule3D.
    Retorna el path del archivo temporal (.pdb).
    """
    if not traj_path or not os.path.exists(traj_path):
        return None
    try:
        traj = Trajectory(traj_path)
        if len(traj) == 0:
            return None

        # Usar último frame para visualización estática
        atoms = traj[-1]

        # Crear archivo PDB temporal (formato aceptado por Molecule3D)
        with tempfile.NamedTemporaryFile(suffix=".pdb", delete=False) as f:
            pdb_path = f.name
        write(pdb_path, atoms)  # extensión .pdb selecciona el formato automáticamente
        return pdb_path

    except Exception as e:
        print(f"Error converting trajectory: {e}")
        return None

# ==== Fallback HTML con 3Dmol.js (por si acaso) ====
def traj_to_html(traj_path, width=520, height=520, interval_ms=200):
    """
    Fallback HTML viewer (solo si Molecule3D no está disponible).
    """
    if not traj_path or not os.path.exists(traj_path):
        return "<div style='color:#b00; padding:20px;'>No trajectory file found</div>"

    viewer_id = f"viewer_{abs(hash(traj_path)) % 10000}"

    try:
        traj = Trajectory(traj_path)
        if len(traj) == 0:
            return "<div style='color:#555; padding:20px;'>Empty trajectory</div>"
    except Exception as e:
        return f"<div style='color:#b00; padding:20px;'>Error: {e}</div>"

    xyz_frames = []
    for atoms in traj:
        symbols = atoms.get_chemical_symbols()
        coords = atoms.get_positions()
        parts = [str(len(symbols)), "frame"]
        for s, (x, y, z) in zip(symbols, coords):
            parts.append(f"{s} {x:.6f} {y:.6f} {z:.6f}")
        xyz_frames.append("\n".join(parts))

    frames_json = str(xyz_frames).replace("'", '"')

    html = f"""
<div style="margin-bottom:10px; padding:10px; background:#f5f5f5; border-radius:5px;">
    <strong>🧬 3D Molecular Viewer</strong> - {len(xyz_frames)} frames
</div>
<div id="{viewer_id}" style="width:{width}px; height:{height}px; position:relative; border:2px solid #ddd; border-radius:8px; background:#fafafa;"></div>
<script>
if (typeof window.$3Dmol === 'undefined') {{
    var script = document.createElement('script');
    script.src = 'https://3dmol.org/build/3Dmol-min.js';
    script.onload = function() {{
        setTimeout(function() {{ initViewer_{viewer_id}(); }}, 100);
    }};
    document.head.appendChild(script);
}} else {{
    initViewer_{viewer_id}();
}}
function initViewer_{viewer_id}() {{
    var el = document.getElementById("{viewer_id}");
    if (!el || typeof $3Dmol === "undefined") return;

    var viewer = $3Dmol.createViewer(el, {{backgroundColor: 'white'}});
    var frames = {frames_json};
    var currentFrame = 0;

    function showFrame(frameIndex) {{
        viewer.clear();
        viewer.addModel(frames[frameIndex], "xyz");
        viewer.setStyle({{}}, {{stick: {{radius: 0.1}}, sphere: {{radius: 0.3}}}});
        viewer.zoomTo();
        viewer.render();
    }}

    showFrame(0);

    if (frames.length > 1) {{
        setInterval(function() {{
            currentFrame = (currentFrame + 1) % frames.length;
            showFrame(currentFrame);
        }}, {interval_ms});
    }}
}}
</script>
"""
    return html

# ==== OrbMol SPE directo (sin cambios) ====
from orb_models.forcefield import pretrained
from orb_models.forcefield.calculator import ORBCalculator

_MODEL_CALC = None
def _load_orbmol_calc():
    global _MODEL_CALC
    if _MODEL_CALC is None:
        orbff = pretrained.orb_v3_conservative_inf_omat(
            device="cpu", precision="float32-high"
        )
        _MODEL_CALC = ORBCalculator(orbff, device="cpu")
    return _MODEL_CALC

def predict_molecule(xyz_content, charge=0, spin_multiplicity=1):
    """
    Single Point Energy + fuerzas. No escribe nada salvo un .xyz temporal.
    """
    try:
        calc = _load_orbmol_calc()
        if not xyz_content or not xyz_content.strip():
            return "Error: Please enter XYZ coordinates", "Error"

        with tempfile.NamedTemporaryFile(mode="w", suffix=".xyz", delete=False) as f:
            f.write(xyz_content)
            xyz_file = f.name

        atoms = read(xyz_file)
        atoms.info = {"charge": int(charge), "spin": int(spin_multiplicity)}
        atoms.calc = calc

        energy = atoms.get_potential_energy()  # eV
        forces = atoms.get_forces()            # eV/Å

        lines = [f"Total Energy: {energy:.6f} eV", "", "Atomic Forces:"]
        for i, fc in enumerate(forces):
            lines.append(f"Atom {i+1}: [{fc[0]:.4f}, {fc[1]:.4f}, {fc[2]:.4f}] eV/Å")
        max_force = float(np.max(np.linalg.norm(forces, axis=1)))
        lines += ["", f"Max Force: {max_force:.4f} eV/Å"]

        try:
            os.unlink(xyz_file)
        except Exception:
            pass

        return "\n".join(lines), "Calculation completed with OrbMol"
    except Exception as e:
        return f"Error during calculation: {e}", "Error"

# ==== Simulaciones (sin cambios) ====
from simulation_scripts_orbmol import (
    run_md_simulation,
    run_relaxation_simulation,
)

def _string_looks_like_xyz(text: str) -> bool:
    try:
        first = (text or "").strip().splitlines()[0]
        int(first.split()[0])
        return True
    except Exception:
        return False

def _to_file_if_xyz(input_or_path: str):
    if isinstance(input_or_path, str) and _string_looks_like_xyz(input_or_path):
        tf = tempfile.NamedTemporaryFile(mode="w", suffix=".xyz", delete=False)
        tf.write(input_or_path)
        tf.flush(); tf.close()
        return tf.name, True
    return input_or_path, False

# Wrappers actualizados para devolver archivos para Molecule3D
def md_wrapper(xyz_content, charge, spin, steps, tempK, timestep_fs, ensemble):
    tmp_created = False
    path_or_str = xyz_content
    try:
        path_or_str, tmp_created = _to_file_if_xyz(xyz_content)

        traj_path, log_text, script_text, explanation = run_md_simulation(
            path_or_str,
            int(steps),
            20,                          # pre-relax fija
            float(timestep_fs),
            float(tempK),
            "NVT" if ensemble == "NVT" else "NVE",
            int(charge),
            int(spin),
        )
        status = f"MD completed: {int(steps)} steps at {int(tempK)} K ({ensemble})"

        # Usar Molecule3D si está disponible, sino HTML
        if HAVE_MOL3D:
            pdb_file = traj_to_molecule3d_file(traj_path)
            return (status, traj_path, log_text, script_text, explanation, pdb_file, "")
        else:
            html_value = traj_to_html(traj_path)
            return (status, traj_path, log_text, script_text, explanation, None, html_value)

    except Exception as e:
        return (f"Error: {e}", None, "", "", "", None, "")
    finally:
        if tmp_created and isinstance(path_or_str, str) and os.path.exists(path_or_str):
            try: os.remove(path_or_str)
            except Exception: pass

def relax_wrapper(xyz_content, steps, fmax, charge, spin, relax_cell):
    tmp_created = False
    path_or_str = xyz_content
    try:
        path_or_str, tmp_created = _to_file_if_xyz(xyz_content)

        traj_path, log_text, script_text, explanation = run_relaxation_simulation(
            path_or_str,
            int(steps),
            float(fmax),
            int(charge),
            int(spin),
            bool(relax_cell),
        )
        status = f"Relaxation finished (≤ {int(steps)} steps, fmax={float(fmax)} eV/Å)"

        # Usar Molecule3D si está disponible, sino HTML
        if HAVE_MOL3D:
            pdb_file = traj_to_molecule3d_file(traj_path)
            return (status, traj_path, log_text, script_text, explanation, pdb_file, "")
        else:
            html_value = traj_to_html(traj_path)
            return (status, traj_path, log_text, script_text, explanation, None, html_value)

    except Exception as e:
        return (f"Error: {e}", None, "", "", "", None, "")
    finally:
        if tmp_created and isinstance(path_or_str, str) and os.path.exists(path_or_str):
            try:
                os.remove(path_or_str)
            except Exception:
                pass

# ==== Ejemplos (sin cambios) ====
examples = [
    ["""2
Hydrogen molecule
H 0.0 0.0 0.0
H 0.0 0.0 0.74""", 0, 1],
    ["""3
Water molecule
O 0.0000 0.0000 0.0000
H 0.7571 0.0000 0.5864
H -0.7571 0.0000 0.5864""", 0, 1],
    ["""5
Methane
C 0.0000 0.0000 0.0000
H 1.0890 0.0000 0.0000
H -0.3630 1.0267 0.0000
H -0.3630 -0.5133 0.8887
H -0.3630 -0.5133 -0.8887""", 0, 1],
]

# ==== UI actualizada ====
with gr.Blocks(theme=gr.themes.Ocean(), title="OrbMol Demo") as demo:
    with gr.Tabs():
        # -------- SPE (sin cambios) --------
        with gr.Tab("Single Point Energy"):
            with gr.Row():
                with gr.Column(scale=2):
                    gr.Markdown("# OrbMol — Quantum-Accurate Molecular Predictions")
                    gr.Markdown("Energías y fuerzas con **charge** y **spin multiplicity**.")
                    xyz_input = gr.Textbox(label="XYZ Coordinates", lines=12, placeholder="Paste XYZ here…")
                    with gr.Row():
                        charge_input = gr.Slider(minimum=-10, maximum=10, value=0, step=1, label="Charge")
                        spin_input = gr.Slider(minimum=1, maximum=11, value=1, step=1, label="Spin Multiplicity")
                    run_spe = gr.Button("Run OrbMol Prediction", variant="primary")
                with gr.Column(variant="panel", min_width=500):
                    spe_out = gr.Textbox(label="Energy & Forces", lines=15, interactive=False)
                    spe_status = gr.Textbox(label="Status", interactive=False, max_lines=1)

            gr.Examples(examples=examples, inputs=[xyz_input, charge_input, spin_input])
            run_spe.click(predict_molecule, [xyz_input, charge_input, spin_input], [spe_out, spe_status])

        # -------- MD (actualizada con Molecule3D) --------
        with gr.Tab("Molecular Dynamics"):
            with gr.Row():
                with gr.Column(scale=2):
                    xyz_md = gr.Textbox(label="XYZ Coordinates or path (.xyz/.traj/.pdb/.cif)", lines=12, placeholder="Paste XYZ or path here…")
                    with gr.Row():
                        charge_md = gr.Slider(minimum=-10, maximum=10, value=0, step=1, label="Charge")
                        spin_md = gr.Slider(minimum=1, maximum=11, value=1, step=1, label="Spin Multiplicity")
                    with gr.Row():
                        steps_md = gr.Slider(minimum=10, maximum=2000, value=100, step=10, label="Steps")
                        temp_md = gr.Slider(minimum=10, maximum=1500, value=300, step=10, label="Temperature (K)")
                    with gr.Row():
                        timestep_md = gr.Slider(minimum=0.1, maximum=5.0, value=1.0, step=0.1, label="Timestep (fs)")
                        ensemble_md = gr.Radio(["NVE", "NVT"], value="NVE", label="Ensemble")
                    run_md_btn = gr.Button("Run MD Simulation", variant="primary")

                with gr.Column(variant="panel", min_width=520):
                    md_status = gr.Textbox(label="MD Status", interactive=False)
                    md_traj = gr.File(label="Trajectory (.traj)", interactive=False)

                    # Usar Molecule3D si está disponible
                    if HAVE_MOL3D:
                        md_viewer = Molecule3D(label="3D Molecular Viewer")
                        md_html = gr.HTML(visible=False)  # Oculto cuando usamos Molecule3D
                    else:
                        md_viewer = gr.HTML(visible=False)  # Placeholder
                        md_html = gr.HTML(label="Trajectory Viewer")

                    md_log = gr.Textbox(label="Log", interactive=False, lines=15, max_lines=25)
                    md_script = gr.Code(label="Reproduction Script", language="python", interactive=False, lines=20, max_lines=30)
                    md_explain = gr.Markdown()

            run_md_btn.click(
                md_wrapper,
                inputs=[xyz_md, charge_md, spin_md, steps_md, temp_md, timestep_md, ensemble_md],
                outputs=[md_status, md_traj, md_log, md_script, md_explain, md_viewer, md_html],
            )

        # -------- Relax (actualizada con Molecule3D) --------
        with gr.Tab("Relaxation / Optimization"):
            with gr.Row():
                with gr.Column(scale=2):
                    xyz_rlx = gr.Textbox(label="XYZ Coordinates or path (.xyz/.traj/.pdb/.cif)", lines=12, placeholder="Paste XYZ or path here…")
                    steps_rlx = gr.Slider(minimum=1, maximum=2000, value=300, step=1, label="Max Steps")
                    fmax_rlx = gr.Slider(minimum=0.001, maximum=0.5, value=0.05, step=0.001, label="Fmax (eV/Å)")
                    with gr.Row():
                        charge_rlx = gr.Slider(minimum=-10, maximum=10, value=0, step=1, label="Charge")
                        spin_rlx = gr.Slider(minimum=1, maximum=11, value=1, step=1, label="Spin")
                    relax_cell = gr.Checkbox(False, label="Relax Unit Cell")
                    run_rlx_btn = gr.Button("Run Optimization", variant="primary")

                with gr.Column(variant="panel", min_width=520):
                    rlx_status = gr.Textbox(label="Status", interactive=False)
                    rlx_traj = gr.File(label="Trajectory (.traj)", interactive=False)

                    # Usar Molecule3D si está disponible
                    if HAVE_MOL3D:
                        rlx_viewer = Molecule3D(label="Final Structure")
                        rlx_html = gr.HTML(visible=False)  # Oculto cuando usamos Molecule3D
                    else:
                        rlx_viewer = gr.HTML(visible=False)  # Placeholder
                        rlx_html = gr.HTML(label="Final Structure")

                    rlx_log = gr.Textbox(label="Log", interactive=False, lines=15, max_lines=25)
                    rlx_script = gr.Code(label="Reproduction Script", language="python", interactive=False, lines=20, max_lines=30)
                    rlx_explain = gr.Markdown()

            run_rlx_btn.click(
                relax_wrapper,
                inputs=[xyz_rlx, steps_rlx, fmax_rlx, charge_rlx, spin_rlx, relax_cell],
                outputs=[rlx_status, rlx_traj, rlx_log, rlx_script, rlx_explain, rlx_viewer, rlx_html],
            )

print("Starting OrbMol model loading…")
_ = _load_orbmol_calc()

if __name__ == "__main__":
    demo.launch(server_name="0.0.0.0", server_port=7860, show_error=True)