Phosphomimetic Bioisosteres.
from pathlib import Path
from datetime import datetime
import csv
import random
import traceback
from rdkit import Chem
from rdkit.Chem import AllChem
from dimorphite_dl import protonate_smiles
# ============================================================
# 1. CONFIGURATION
# ============================================================
INPUT_DIR = Path(
r"D:\encinar\tempjant\115-CD44-inhibitors\03-chemical-library"
)
OUTPUT_DIR = Path(
r"D:\encinar\tempjant\115-CD44-inhibitors\05-library-ASM-500-ion"
)
# Protonation conditions
PH_MIN = 3.9
PH_MAX = 4.1
PKA_PRECISION = 1.0
MAX_VARIANTS = 128
# Output
MOLECULES_PER_FILE = 500
RECURSIVE_SEARCH = False
# 3D generation and optimization
RANDOM_SEED = 12345
MAX_EMBED_ATTEMPTS = 3
MMFF_MAX_ITERATIONS = 500
# Reproducibility
random.seed(RANDOM_SEED)
# ============================================================
# 2. ZINC-BINDING GROUPS (ZBGs)
# ============================================================
ZBG_SMARTS = {
"Hydroxamic_Acid": "C(=O)NO",
"Primary_Sulfonamide": "S(=O)(=O)N",
"N_Hydroxyurea": "NC(=O)NO",
"1_2_Hydroxypyridinone": "O=C1C=CC=CN1O",
"Methylene_Phosphonate": "CP(=O)(O)O",
"Difluoromethylphosphonate": "FC(F)P(=O)(O)O",
"Bisphosphonate": "O=P(O)(O)CP(=O)(O)O",
"Phosphonamidate": "P(=O)(O)N",
"Thiophosphate": "SP(=O)(O)O",
"Sulfonate": "S(=O)(=O)O",
"Sulfate": "OS(=O)(=O)O",
"Carboxylate": "C(=O)O",
"Dicarboxylate_Malonate": "O=C(O)CC(=O)O",
"Boronic_Acid": "B(O)O",
"Tetrazole": "c1nnnn1",
}
COMPILED_ZBG_SMARTS = {}
for name, smarts in ZBG_SMARTS.items():
pattern = Chem.MolFromSmarts(smarts)
if pattern is not None:
COMPILED_ZBG_SMARTS[name] = pattern
else:
print(f"WARNING: Invalid SMARTS pattern: {name}")
# ============================================================
# 3. GLOBAL VARIABLES
# ============================================================
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
run_timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
output_prefix = f"library_ionized_pH4_{run_timestamp}"
seen_variants = set()
molecule_batch = []
output_file_number = 1
statistics = {
"input_files": 0,
"input_molecules": 0,
"generated_variants": 0,
"charged_variants": 0,
"uncharged_variants_discarded": 0,
"duplicate_variants": 0,
"embedding_failures": 0,
"optimization_failures": 0,
"output_molecules": 0,
}
# ============================================================
# 4. MOLECULE UTILITIES
# ============================================================
def get_formal_charge(mol):
"""Return total formal molecular charge."""
return sum(
atom.GetFormalCharge()
for atom in mol.GetAtoms()
)
def get_number_of_charged_atoms(mol):
"""Count atoms with non-zero formal charge."""
return sum(
1
for atom in mol.GetAtoms()
if atom.GetFormalCharge() != 0
)
def has_charged_atoms(mol):
"""Return True if at least one atom has formal charge."""
return any(
atom.GetFormalCharge() != 0
for atom in mol.GetAtoms()
)
def identify_zbgs(mol):
"""Identify zinc-binding group motifs using SMARTS."""
found_groups = []
for name, pattern in COMPILED_ZBG_SMARTS.items():
if mol.HasSubstructMatch(pattern):
found_groups.append(name)
return found_groups
# ============================================================
# 5. GENERATE 3D STRUCTURE AND OPTIMIZE
# ============================================================
def generate_3d_conformer(smiles, molecule_name):
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None, "Invalid_SMILES"
try:
# Add explicit hydrogens.
mol = Chem.AddHs(mol)
# Generate 3D coordinates.
embedded = False
for attempt in range(MAX_EMBED_ATTEMPTS):
params = AllChem.ETKDGv3()
params.randomSeed = RANDOM_SEED + attempt
params.useRandomCoords = False
params.maxIterations = 1000
result = AllChem.EmbedMolecule(
mol,
params
)
if result == 0:
embedded = True
break
if not embedded:
return None, "Embedding_Failed"
# Geometry optimization.
optimization_status = "Not_Optimized"
if AllChem.MMFFHasAllMoleculeParams(mol):
result = AllChem.MMFFOptimizeMolecule(
mol,
maxIters=MMFF_MAX_ITERATIONS,
mmffVariant="MMFF94"
)
optimization_status = "MMFF94"
if result < 0:
optimization_status = "MMFF_Failed"
elif AllChem.UFFHasAllMoleculeParams(mol):
result = AllChem.UFFOptimizeMolecule(
mol,
maxIters=MMFF_MAX_ITERATIONS
)
optimization_status = "UFF"
if result < 0:
optimization_status = "UFF_Failed"
# IMPORTANT:
# Do NOT remove hydrogens from the output molecule.
# Explicit hydrogen atoms and their 3D coordinates
# will be preserved in the final SDF.
# Remove internal RDKit properties.
for prop in mol.GetPropNames(
includePrivate=True,
includeComputed=True
):
mol.ClearProp(prop)
# Generate canonical SMILES without changing the
# original molecule containing explicit hydrogens.
canonical_smiles = Chem.MolToSmiles(
Chem.RemoveHs(mol),
canonical=True,
isomericSmiles=True
)
# Restore molecule name and SMILES property.
mol.SetProp("_Name", molecule_name)
mol.SetProp("SMILES", canonical_smiles)
return mol, optimization_status
except Exception as error:
return None, str(error)
# ============================================================
# 6. WRITE OUTPUT SDF FILES
# ============================================================
def write_batch(molecules, file_number):
if not molecules:
return
output_filename = (
f"{output_prefix}_{file_number:03d}.sdf"
)
output_path = OUTPUT_DIR / output_filename
writer = Chem.SDWriter(str(output_path))
for mol in molecules:
writer.write(mol)
writer.close()
print(
f"Written {len(molecules)} molecules: "
f"{output_filename}"
)
# ============================================================
# 7. SAVE A SUCCESSFULLY PROCESSED MOLECULE
# ============================================================
def save_molecule(
mol,
source_file,
source_index,
variant_index,
molecule_name,
net_charge,
charged_atoms,
zbg_list,
optimization_status,
canonical_smiles,
csv_writer
):
global molecule_batch
global output_file_number
molecule_batch.append(mol)
statistics["output_molecules"] += 1
csv_writer.writerow({
"Molecule_Name": molecule_name,
"Source_File": source_file,
"Source_Molecule_Index": source_index,
"Variant_Index": variant_index,
"Net_Formal_Charge": net_charge,
"Number_of_Charged_Atoms": charged_atoms,
"ZBGs": ";".join(zbg_list),
"Optimization_Method": optimization_status,
"SMILES": canonical_smiles,
})
if len(molecule_batch) >= MOLECULES_PER_FILE:
write_batch(
molecule_batch,
output_file_number
)
molecule_batch = []
output_file_number += 1
# ============================================================
# 8. PROCESS ONE INPUT MOLECULE
# ============================================================
def process_molecule(
original_mol,
source_file,
source_index,
csv_writer
):
statistics["input_molecules"] += 1
# Recover original molecule name.
if original_mol.HasProp("_Name"):
molecule_name = original_mol.GetProp("_Name").strip()
else:
molecule_name = ""
if not molecule_name:
molecule_name = (
f"{Path(source_file).stem}_"
f"Molecule_{source_index:06d}"
)
# Remove input hydrogens before protonation enumeration.
original_mol = Chem.RemoveHs(original_mol)
try:
original_smiles = Chem.MolToSmiles(
original_mol,
canonical=True,
isomericSmiles=True
)
except Exception as error:
print(
f"Could not generate SMILES for "
f"{molecule_name}: {error}"
)
return
if not original_smiles:
return
# Enumerate protonation states around pH 4.
try:
variants = protonate_smiles(
original_smiles,
ph_min=PH_MIN,
ph_max=PH_MAX,
precision=PKA_PRECISION,
max_variants=MAX_VARIANTS,
validate_output=True
)
except Exception as error:
print(
f"Protonation failed for {molecule_name}: {error}"
)
return
if not variants:
return
statistics["generated_variants"] += len(variants)
for variant_index, variant_smiles in enumerate(
variants,
start=1
):
variant_mol = Chem.MolFromSmiles(variant_smiles)
if variant_mol is None:
continue
# Retain molecules with at least one charged atom.
# Net charge may be zero (zwitterions are retained).
if not has_charged_atoms(variant_mol):
statistics["uncharged_variants_discarded"] += 1
continue
statistics["charged_variants"] += 1
net_charge = get_formal_charge(variant_mol)
charged_atoms = get_number_of_charged_atoms(
variant_mol
)
canonical_smiles = Chem.MolToSmiles(
variant_mol,
canonical=True,
isomericSmiles=True
)
# Remove duplicate variants only within the same
# source molecule, not between different input records.
duplicate_key = (
str(Path(source_file).resolve()),
source_index,
canonical_smiles
)
if duplicate_key in seen_variants:
statistics["duplicate_variants"] += 1
continue
# Identify ZBG motifs.
zbg_list = identify_zbgs(variant_mol)
# Generate a new 3D structure with explicit hydrogens.
final_mol, optimization_status = generate_3d_conformer(
canonical_smiles,
molecule_name
)
if final_mol is None:
if optimization_status == "Embedding_Failed":
statistics["embedding_failures"] += 1
print(
f"3D generation failed for "
f"{molecule_name}, variant {variant_index}: "
f"{optimization_status}"
)
continue
if "Failed" in optimization_status:
statistics["optimization_failures"] += 1
# Register only successfully generated variants.
seen_variants.add(duplicate_key)
save_molecule(
mol=final_mol,
source_file=source_file,
source_index=source_index,
variant_index=variant_index,
molecule_name=molecule_name,
net_charge=net_charge,
charged_atoms=charged_atoms,
zbg_list=zbg_list,
optimization_status=optimization_status,
canonical_smiles=canonical_smiles,
csv_writer=csv_writer
)
# ============================================================
# 9. MAIN PIPELINE
# ============================================================
def main():
if not INPUT_DIR.exists():
raise FileNotFoundError(
f"Input directory does not exist: {INPUT_DIR}"
)
OUTPUT_DIR.mkdir(
parents=True,
exist_ok=True
)
# Find input SDF files.
if RECURSIVE_SEARCH:
sdf_files = sorted(INPUT_DIR.rglob("*.sdf"))
else:
sdf_files = sorted(INPUT_DIR.glob("*.sdf"))
if not sdf_files:
print("No SDF files found.")
return
statistics["input_files"] = len(sdf_files)
print("=" * 65)
print("CD44 INHIBITOR LIBRARY PREPARATION")
print("=" * 65)
print(f"Input directory: {INPUT_DIR}")
print(f"Output directory: {OUTPUT_DIR}")
print(f"Input SDF files: {len(sdf_files)}")
print(f"Target pH: {PH_MIN}–{PH_MAX}")
print(f"Max molecules/file: {MOLECULES_PER_FILE}")
print("=" * 65)
csv_filename = (
OUTPUT_DIR /
f"{output_prefix}_report.csv"
)
csv_fields = [
"Molecule_Name",
"Source_File",
"Source_Molecule_Index",
"Variant_Index",
"Net_Formal_Charge",
"Number_of_Charged_Atoms",
"ZBGs",
"Optimization_Method",
"SMILES",
]
with open(
csv_filename,
"w",
newline="",
encoding="utf-8-sig"
) as csv_file:
csv_writer = csv.DictWriter(
csv_file,
fieldnames=csv_fields
)
csv_writer.writeheader()
for file_index, sdf_file in enumerate(
sdf_files,
start=1
):
print()
print(
f"[{file_index}/{len(sdf_files)}] "
f"Processing: {sdf_file.name}"
)
supplier = Chem.SDMolSupplier(
str(sdf_file),
removeHs=False,
sanitize=True
)
for molecule_index, mol in enumerate(
supplier,
start=1
):
if mol is None:
print(
f"Invalid molecule in "
f"{sdf_file.name}, record "
f"{molecule_index}"
)
continue
try:
process_molecule(
original_mol=mol,
source_file=sdf_file.name,
source_index=molecule_index,
csv_writer=csv_writer
)
except Exception:
print(
f"Unexpected error in "
f"{sdf_file.name}, record "
f"{molecule_index}"
)
traceback.print_exc()
# Write final batch if it contains fewer than 500 molecules.
if molecule_batch:
write_batch(
molecule_batch,
output_file_number
)
molecule_batch.clear()
# Final statistics.
print()
print("=" * 65)
print("PROCESSING COMPLETED")
print("=" * 65)
for key, value in statistics.items():
print(f"{key:35s}: {value}")
print()
print(f"CSV report: {csv_filename}")
print(f"Output directory: {OUTPUT_DIR}")
print("=" * 65)
# ============================================================
# 10. EXECUTION
# ============================================================
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
print("\nExecution interrupted by user.")
except Exception as error:
print(f"\nFATAL ERROR: {error}")
traceback.print_exc()
Functionality of the next Python script
- Input Parsing: Reads the molecule identifier from the first column and the SMILES string from the second column. It supports both tab- and space-delimited formats and handles input files encoded in UTF-8 or Windows character sets.
- 3D Conformer Generation: Adds explicit hydrogen atoms, generates 3D geometries using the ETKDGv3 algorithm, and optimizes structures via the MMFF94 force field (falling back to UFF when MMFF94 parameters are unavailable).
- Charge Preservation: Formal charges are strictly maintained without automatic neutralization. Validation confirmed charge retention across diverse species, including tetrazolate (N−), ammonium (N+), carboxylate (O−), and zwitterionic structures.
- File Segmentation: Outputs partitioned SDF files named
ASM-ion_0001.sdf,ASM-ion_0002.sdf, etc., each containing a maximum of 500 molecules. - Embedded Metadata: Each entry in the output SDF includes the compound identifier alongside three custom properties: the original SMILES string (
SMILES_original), the net formal charge (Carga_neta), and the force field optimization method applied (Optimizacion). - Error Logging: Processing failures (e.g., invalid SMILES syntax, malformed lines missing required columns, or failed 3D conformer generation) are logged in
ASM-ion_errores.txtwith their line numbers and error descriptions. Failed records do not reserve slots in the output files; SDF files are populated exclusively with valid structures, meaning the final file may contain fewer than 500 entries. Header lines in input files will be flagged in this error log as expected.
Considerations and Technical Notes
- Delimiter Sensitivity: If a molecule identifier contains spaces within a space-delimited file, the string will be truncated at the first space. Using tab-delimited files prevents this issue.
- Stereochemistry: If a SMILES string lacks stereochemical specification at chiral centers, the generator assigns a configuration automatically. Although fully reproducible due to a fixed random seed (
SEED = 42), this choice is arbitrary rather than an informed stereochemical selection. - File Overwriting: Any existing
ASM-ion_XXXX.sdffiles in the working directory will be overwritten during execution. - Execution Performance: The script processes molecules sequentially on a single thread. This approach is sufficient for datasets of a few thousand molecules; parallelization can be implemented if required for larger chemical libraries.
"""
Genera ficheros SDF 3D (con hidrógenos explícitos) a partir de un TXT con
nombre y SMILES, conservando las cargas formales del SMILES original.
Entrada : TXT con el nombre en la 1.ª columna y el SMILES en la 2.ª
(separados por tabulador o espacios; el nombre puede llevar espacios
solo si el separador es un tabulador).
Salida : ASM-ion_0001.sdf, ASM-ion_0002.sdf, ... con MOLS_POR_FICHERO moléculas
cada uno, más un fichero de errores con las moléculas que fallen.
Las moléculas que fallan NO ocupan hueco: cada fichero se rellena solo con
moléculas generadas correctamente. Las cargas no se neutralizan en ningún paso.
Requisito: pip install rdkit (o conda install -c conda-forge rdkit)
"""
import sys
import time
from pathlib import Path
from rdkit import Chem, RDLogger
from rdkit.Chem import AllChem, rdDistGeom
RDLogger.DisableLog('rdApp.*')
# --- CONFIGURACIÓN ---
DIR_SALIDA = Path(r'D:\encinar\tempjant\115-CD44-inhibitors\05-library-ASM-500-ion')
INPUT_FILE = DIR_SALIDA / 'asm-500-pH4.txt'
PREFIJO = 'ASM-ion'
MOLS_POR_FICHERO = 500
SEMILLA = 42 # semilla fija: mismos resultados en cada ejecución
MAX_ITER_OPT = 2000
PROGRESO_CADA = 100
def leer_lineas(ruta):
"""Lee el fichero probando UTF-8 y, si falla, la codificación de Windows."""
for codificacion in ('utf-8-sig', 'cp1252'):
try:
return ruta.read_text(encoding=codificacion).splitlines()
except UnicodeDecodeError:
continue
return ruta.read_text(encoding='latin-1').splitlines()
def separar_linea(linea):
"""Devuelve (nombre, smiles) o None si la línea no tiene dos columnas."""
if '\t' in linea:
partes = [p.strip() for p in linea.split('\t') if p.strip()]
else:
partes = linea.split()
if len(partes) < 2:
return None
return partes[0], partes[1]
def optimizar(mol):
"""Minimiza la geometría con MMFF94 (o UFF si faltan parámetros)."""
try:
if AllChem.MMFFHasAllMoleculeParams(mol):
AllChem.MMFFOptimizeMolecule(mol, maxIters=MAX_ITER_OPT)
return 'MMFF94'
if AllChem.UFFHasAllMoleculeParams(mol):
AllChem.UFFOptimizeMolecule(mol, maxIters=MAX_ITER_OPT)
return 'UFF'
except Exception:
pass
return 'ninguna (geometría ETKDG)'
def generar_3d(nombre, smiles):
"""Devuelve (mol_3D_con_H, None) o (None, motivo_del_error)."""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None, 'SMILES inválido (RDKit no pudo interpretarlo)'
carga_neta = Chem.GetFormalCharge(mol)
mol = Chem.AddHs(mol) # hidrógenos explícitos; no altera las cargas formales
params = rdDistGeom.ETKDGv3()
params.randomSeed = SEMILLA
if AllChem.EmbedMolecule(mol, params) == -1:
params.useRandomCoords = True # segundo intento, más robusto
if AllChem.EmbedMolecule(mol, params) == -1:
return None, 'No se pudo generar la geometría 3D'
metodo = optimizar(mol)
mol.SetProp('_Name', nombre)
mol.SetProp('SMILES_original', smiles)
mol.SetIntProp('Carga_neta', carga_neta)
mol.SetProp('Optimizacion', metodo)
return mol, None
def procesar(ruta_entrada, dir_salida, prefijo, por_fichero):
dir_salida.mkdir(parents=True, exist_ok=True)
lineas = leer_lineas(ruta_entrada)
print(f'{len(lineas)} líneas leídas de {ruta_entrada}', flush=True)
errores = []
escritas = 0
n_fichero = 0
en_fichero = 0
escritor = None
ficheros = []
inicio = time.time()
try:
for n_linea, linea in enumerate(lineas, start=1):
if not linea.strip():
continue
campos = separar_linea(linea)
if campos is None:
errores.append((n_linea, linea.strip(), '', 'La línea no tiene dos columnas'))
continue
nombre, smiles = campos
try:
mol, motivo = generar_3d(nombre, smiles)
except Exception as e: # una molécula problemática no debe parar el proceso
mol, motivo = None, f'Error inesperado: {e}'
if mol is None:
errores.append((n_linea, nombre, smiles, motivo))
continue
if escritor is None or en_fichero >= por_fichero:
if escritor is not None:
escritor.close()
n_fichero += 1
ruta = dir_salida / f'{prefijo}_{n_fichero:04d}.sdf'
ficheros.append(ruta)
escritor = Chem.SDWriter(str(ruta))
en_fichero = 0
escritor.write(mol)
en_fichero += 1
escritas += 1
if escritas % PROGRESO_CADA == 0:
print(
f' {escritas} moléculas generadas ({time.time() - inicio:.0f} s)',
flush=True,
)
finally:
if escritor is not None:
escritor.close()
ruta_err = dir_salida / f'{prefijo}_errores.txt'
if errores:
with open(ruta_err, 'w', encoding='utf-8') as f:
f.write('linea\tnombre\tsmiles\tmotivo\n')
for fila in errores:
f.write('\t'.join(str(x) for x in fila) + '\n')
print('=' * 60)
print(f'Moléculas generadas : {escritas}')
print(f'Moléculas con error : {len(errores)}' + (f' -> {ruta_err}' if errores else ''))
print(f'Ficheros SDF : {len(ficheros)} en {dir_salida}')
print(f'Tiempo total : {time.time() - inicio:.0f} s')
print('=' * 60)
def main():
if not INPUT_FILE.exists():
print(f'Error: no se encontró {INPUT_FILE}')
sys.exit(1)
procesar(INPUT_FILE, DIR_SALIDA, PREFIJO, MOLS_POR_FICHERO)
if __name__ == '__main__':
main()
Linux version for HPC cluster with 32 CPUs and 64 GB of RAM memory.
"""
Parallel (HPC / Linux) version of generar_sdf_3d.py.
Builds 3D SDF files (explicit hydrogens, formal charges kept exactly as in the
input SMILES) from a text file with the molecule name in column 1 and the SMILES
in column 2. Molecules are distributed over N worker processes and results are
written IN INPUT ORDER, so the output matches the serial version (same RDKit
version and same seed).
Usage:
python generar_sdf_3d_hpc.py INPUT.txt OUTPUT_DIR [--processes 32]
[--prefix ASM-ion] [--per-file 500] [--seed 42]
(The older Spanish option names --procesos, --prefijo, --por-fichero and
--semilla are still accepted as aliases.)
If --processes is not given, SLURM_CPUS_PER_TASK is used or, failing that, the
number of CPUs this process is allowed to use.
Output: ASM-ion_0001.sdf, ASM-ion_0002.sdf, ... and ASM-ion_errores.txt.
Failed molecules do not take a slot: every SDF file is filled only with
successfully generated molecules (the last file may hold fewer).
Note: this file is plain ASCII on purpose, so it works whatever encoding it is
saved or copied with.
"""
import argparse
import io
import multiprocessing as mp
import os
import sys
import time
from pathlib import Path
from rdkit import Chem, RDLogger
from rdkit.Chem import AllChem, rdDistGeom
MAX_OPT_ITERS = 2000
EMBED_TIMEOUT_S = 120 # only applied if your RDKit version supports it
PROGRESS_EVERY = 2000
_SEED = 42 # set in every worker process by _init_worker
# ----------------------------------------------------------------------------
# Input
# ----------------------------------------------------------------------------
def read_lines(path):
"""Read the file trying UTF-8 first, then the Windows code page."""
for encoding in ('utf-8-sig', 'cp1252'):
try:
return path.read_text(encoding=encoding).splitlines()
except UnicodeDecodeError:
continue
return path.read_text(encoding='latin-1').splitlines()
def split_line(line):
"""Return (name, smiles) or None if the line does not have two columns."""
if '\t' in line:
parts = [p.strip() for p in line.split('\t') if p.strip()]
else:
parts = line.split()
if len(parts) < 2:
return None
return parts[0], parts[1]
# ----------------------------------------------------------------------------
# Work done by each worker process
# ----------------------------------------------------------------------------
def _init_worker(seed):
global _SEED
_SEED = seed
RDLogger.DisableLog('rdApp.*')
def optimize(mol):
"""Minimise the geometry with MMFF94 (or UFF if parameters are missing)."""
try:
if AllChem.MMFFHasAllMoleculeParams(mol):
AllChem.MMFFOptimizeMolecule(mol, maxIters=MAX_OPT_ITERS)
return 'MMFF94'
if AllChem.UFFHasAllMoleculeParams(mol):
AllChem.UFFOptimizeMolecule(mol, maxIters=MAX_OPT_ITERS)
return 'UFF'
except Exception:
pass
return 'none (ETKDG geometry)'
def generate_3d(name, smiles):
"""Return (3D mol with H, None) or (None, error reason)."""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None, 'Invalid SMILES (RDKit could not parse it)'
net_charge = Chem.GetFormalCharge(mol)
mol = Chem.AddHs(mol) # explicit hydrogens; formal charges are untouched
params = rdDistGeom.ETKDGv3()
params.randomSeed = _SEED
params.numThreads = 1 # parallelism is done per molecule, by this script
if hasattr(params, 'timeout'):
params.timeout = EMBED_TIMEOUT_S
if AllChem.EmbedMolecule(mol, params) == -1:
params.useRandomCoords = True # second, more robust attempt
if AllChem.EmbedMolecule(mol, params) == -1:
return None, 'Could not generate 3D geometry'
method = optimize(mol)
# SDF property names are kept identical to the Windows version on purpose
mol.SetProp('_Name', name)
mol.SetProp('SMILES_original', smiles)
mol.SetIntProp('Carga_neta', net_charge)
mol.SetProp('Optimizacion', method)
return mol, None
def to_sdf_record(mol):
"""Full SDF record text for one molecule (ends with $$$$)."""
buffer = io.StringIO()
writer = Chem.SDWriter(buffer)
writer.write(mol)
writer.close()
return buffer.getvalue()
def work(task):
"""Process one molecule. Return (line_no, name, smiles, record, reason)."""
line_no, name, smiles = task
try:
mol, reason = generate_3d(name, smiles)
if mol is None:
return line_no, name, smiles, None, reason
return line_no, name, smiles, to_sdf_record(mol), None
except Exception as e: # one bad molecule must not stop the whole job
return line_no, name, smiles, None, f'Unexpected error: {e}'
# ----------------------------------------------------------------------------
# Main process
# ----------------------------------------------------------------------------
def available_cpus():
env = os.environ.get('SLURM_CPUS_PER_TASK')
if env and env.isdigit():
return int(env)
try:
return len(os.sched_getaffinity(0))
except AttributeError:
return os.cpu_count() or 1
def run(input_path, out_dir, prefix, per_file, processes, seed):
out_dir.mkdir(parents=True, exist_ok=True)
lines = read_lines(input_path)
tasks, errors = [], []
for line_no, line in enumerate(lines, start=1):
if not line.strip():
continue
fields = split_line(line)
if fields is None:
errors.append((line_no, line.strip(), '', 'Line does not have two columns'))
continue
tasks.append((line_no, fields[0], fields[1]))
total = len(tasks)
print(f'{len(lines)} lines read; {total} molecules to process with {processes} processes.', flush=True)
chunksize = max(1, min(25, total // (processes * 20) or 1))
written = 0
file_no = 0
in_file = 0
out = None
files = []
start = time.time()
try:
with mp.Pool(processes, initializer=_init_worker, initargs=(seed,)) as pool:
# imap (not imap_unordered): keeps the input order
for line_no, name, smiles, record, reason in pool.imap(work, tasks, chunksize):
if record is None:
errors.append((line_no, name, smiles, reason))
else:
if out is None or in_file >= per_file:
if out is not None:
out.close()
file_no += 1
path = out_dir / f'{prefix}_{file_no:04d}.sdf'
files.append(path)
out = open(path, 'w', encoding='utf-8', newline='\n')
in_file = 0
out.write(record)
in_file += 1
written += 1
done = written + len(errors)
if done % PROGRESS_EVERY == 0:
elapsed = time.time() - start
rate = done / elapsed if elapsed > 0 else 0
eta = (total - done) / rate / 60 if rate > 0 else 0
print(f' {done}/{total} ({rate:.1f} mol/s, ~{eta:.0f} min left)', flush=True)
finally:
if out is not None:
out.close()
err_path = out_dir / f'{prefix}_errores.txt'
if errors:
errors.sort(key=lambda e: e[0])
with open(err_path, 'w', encoding='utf-8') as f:
f.write('line\tname\tsmiles\treason\n')
for row in errors:
f.write('\t'.join(str(x) for x in row) + '\n')
print('=' * 60)
print(f'Molecules written : {written}')
print(f'Molecules failed : {len(errors)}' + (f' -> {err_path}' if errors else ''))
print(f'SDF files : {len(files)} in {out_dir}')
print(f'Total time : {(time.time() - start) / 60:.1f} min')
print('=' * 60)
def main():
ap = argparse.ArgumentParser(description='3D SDF files with H and charges, in parallel.')
ap.add_argument('input', type=Path, help='text file: name and SMILES')
ap.add_argument('output', type=Path, help='output directory')
ap.add_argument('--processes', '--procesos', type=int, default=None)
ap.add_argument('--prefix', '--prefijo', default='ASM-ion')
ap.add_argument('--per-file', '--por-fichero', type=int, default=500)
ap.add_argument('--seed', '--semilla', type=int, default=42)
args = ap.parse_args()
if not args.input.exists():
print(f'Error: file not found: {args.input}')
sys.exit(1)
processes = args.processes or available_cpus()
run(args.input, args.output, args.prefix, args.per_file, processes, args.seed)
if __name__ == '__main__':
main()
# el fichero sbatch usado para ejecutar el scipt anterior:
#!/bin/bash
#SBATCH --job-name=sdf3d
#SBATCH --qos=long
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --mem=16G
#SBATCH --time=02:00:00
#SBATCH --output=sdf3d_%j.log
# --- AJUSTA ESTAS LÍNEAS A TU CLÚSTER -------------------------------------
# Si algo falla (p. ej. la activación del entorno), el trabajo se detiene aquí
set -e
# Inicializa conda (conda.sh está en la instalación BASE de Miniconda, no en el entorno)
source /home/jant.encinar/miniconda3/etc/profile.d/conda.sh
conda activate vina
ENTRADA=/home/jant.encinar/name-SMILES_to_SDF-3D-H-500/asm-500-pH4.txt
SALIDA=/home/jant.encinar/name-SMILES_to_SDF-3D-H-500/05-library-ASM-500-ion
# Un solo hilo por proceso: el paralelismo lo gestiona el script
export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
# El trabajo se ejecuta desde la carpeta donde lanzaste sbatch (ahí debe estar el .py)
cd "$SLURM_SUBMIT_DIR"
# Comprobación que queda en el log: ¿es el Python del entorno y tiene RDKit?
echo "Python: $(which python)"
python -c "import rdkit; print('RDKit', rdkit.__version__)"
python3 generar_sdf_3d_hpc.py "$ENTRADA" "$SALIDA" --procesos "${SLURM_CPUS_PER_TASK}"
The ADMETsar 3.0 for the chemical library can be download from the next link as xls file.
Access to the chemical library in SDF format (strictly for scientific research purposes) is available at the following link.