#!/usr/bin/env python3"""Pareto-Aware Multi-Objective Molecular Optimization Generates candidates and identifies the Pareto front across multiple property objectives.Based on MultiMol (Yu 2025) and MOLLM (Ran 2025). Usage: python pareto_optimize.py --smiles "CCO" --objectives "LogP:minimize:3.0,QED:maximize:0.5" --output results.json python pareto_optimize.py --input candidates.csv --objectives "LogP:minimize:3.0,QED:maximize:0.5" --mode analyze""" import argparseimport jsonimport sys try: from rdkit import Chem from rdkit.Chem import Descriptors, QED, AllChem, DataStructs from rdkit.Chem.Scaffolds import MurckoScaffold from rdkit import RDLogger RDLogger.logger().setLevel(RDLogger.ERROR)except ImportError: print("Error: rdkit-pypi required.", file=sys.stderr) sys.exit(1) import numpy as np PROPERTY_FUNCS = { 'MW': lambda m: Descriptors.MolWt(m), 'LogP': lambda m: Descriptors.MolLogP(m), 'TPSA': lambda m: Descriptors.TPSA(m), 'HBA': lambda m: float(Descriptors.NumHAcceptors(m)), 'HBD': lambda m: float(Descriptors.NumHDonors(m)), 'RotBonds': lambda m: float(Descriptors.NumRotatableBonds(m)), 'QED': lambda m: QED.qed(m), 'AromaticRings': lambda m: float(Descriptors.NumAromaticRings(m)), 'FractionCSP3': lambda m: Descriptors.FractionCSP3(m), 'NumRings': lambda m: float(Descriptors.RingCount(m)),} BIOISOSTERES = [ ('[cH]', '[n]', 'C→N in ring'), ('Cl', 'F', 'Cl→F'), ('OC', 'O', 'OMe→OH'), ('CC', 'C', 'shorten chain'), ('c1ccccc1', 'c1ccncc1', 'phenyl→pyridine'), ('C(=O)N', 'S(=O)(=O)N', 'amide→sulfonamide'),] def parse_objectives(obj_str): """Parse objective string like 'LogP:minimize:3.0,QED:maximize:0.5'.""" objectives = [] for item in obj_str.split(','): parts = item.strip().split(':') if len(parts) == 3: prop, direction, target = parts objectives.append({ 'property': prop.strip(), 'direction': direction.strip(), 'target': float(target.strip()), }) elif len(parts) == 4: # range: TPSA:range:20:130 prop, _, lo, hi = parts objectives.append({ 'property': prop.strip(), 'direction': 'range', 'target_min': float(lo.strip()), 'target_max': float(hi.strip()), }) return objectives def evaluate_molecule(mol, objectives): """Evaluate a molecule against all objectives. Returns scores (lower = better).""" scores = {} props = {} for obj in objectives: prop = obj['property'] if prop not in PROPERTY_FUNCS: continue val = PROPERTY_FUNCS[prop](mol) props[prop] = round(val, 4) if obj['direction'] == 'minimize': score = max(0, val - obj['target']) elif obj['direction'] == 'maximize': score = max(0, obj['target'] - val) elif obj['direction'] == 'range': if val < obj['target_min']: score = obj['target_min'] - val elif val > obj['target_max']: score = val - obj['target_max'] else: score = 0 else: score = abs(val - obj.get('target', val)) scores[prop] = round(score, 4) return props, scores def dominates(scores_a, scores_b): """Check if solution A Pareto-dominates solution B (all scores lower = better).""" a_vals = list(scores_a.values()) b_vals = list(scores_b.values()) at_least_one_better = False for a, b in zip(a_vals, b_vals): if a > b: return False if a < b: at_least_one_better = True return at_least_one_better def find_pareto_front(candidates): """Identify the Pareto front from scored candidates.""" pareto = [] dominated = [] for i, cand_i in enumerate(candidates): is_dominated = False for j, cand_j in enumerate(candidates): if i == j: continue if dominates(cand_j['scores'], cand_i['scores']): is_dominated = True break if is_dominated: dominated.append(cand_i) else: pareto.append(cand_i) return pareto, dominated def generate_candidates(smiles, max_candidates=16): """Generate candidate molecules via bioisosteric replacements.""" mol = Chem.MolFromSmiles(smiles) if mol is None: return [] candidates = [] for old_pat, new_pat, rationale in BIOISOSTERES: if len(candidates) >= max_candidates: break try: old_mol = Chem.MolFromSmarts(old_pat) or Chem.MolFromSmiles(old_pat) new_mol = Chem.MolFromSmiles(new_pat) or Chem.MolFromSmarts(new_pat) if old_mol is None or new_mol is None: continue if mol.HasSubstructMatch(old_mol): products = AllChem.ReplaceSubstructs(mol, old_mol, new_mol) for prod in products: try: smi = Chem.MolToSmiles(prod) if Chem.MolFromSmiles(smi) is not None and smi != smiles: candidates.append({'smiles': smi, 'rationale': rationale}) except Exception: continue except Exception: continue # Deduplicate seen = set() unique = [] for c in candidates: if c['smiles'] not in seen: seen.add(c['smiles']) unique.append(c) return unique[:max_candidates] def main(): parser = argparse.ArgumentParser(description='Multi-objective molecular optimization') parser.add_argument('--smiles', type=str, help='Input SMILES to optimize') parser.add_argument('--input', type=str, help='CSV with candidate molecules') parser.add_argument('--smiles-col', type=str, default='SMILES') parser.add_argument('--objectives', type=str, required=True, help='Objectives (e.g., "LogP:minimize:3.0,QED:maximize:0.5")') parser.add_argument('--candidates', type=int, default=16, help='Max candidates to generate') parser.add_argument('--mode', choices=['optimize', 'analyze'], default='optimize') parser.add_argument('--output', type=str, default='pareto_results.json') args = parser.parse_args() objectives = parse_objectives(args.objectives) if not objectives: print("Error: No valid objectives parsed.", file=sys.stderr) return 1 molecules = [] if args.smiles and args.mode == 'optimize': ref_mol = Chem.MolFromSmiles(args.smiles) if ref_mol is None: print(f"Error: Invalid SMILES: {args.smiles}", file=sys.stderr) return 1 ref_props, ref_scores = evaluate_molecule(ref_mol, objectives) print(f"Reference: {Chem.MolToSmiles(ref_mol)}") print(f"Properties: {ref_props}") print(f"Objective gaps: {ref_scores}") candidates = generate_candidates(args.smiles, args.candidates) molecules = [{'smiles': args.smiles, 'is_reference': True}] + candidates elif args.input: import pandas as pd df = pd.read_csv(args.input) molecules = [{'smiles': smi} for smi in df[args.smiles_col]] else: print("Error: Provide --smiles or --input", file=sys.stderr) return 1 # Evaluate all molecules evaluated = [] for entry in molecules: mol = Chem.MolFromSmiles(entry['smiles']) if mol is None: continue props, scores = evaluate_molecule(mol, objectives) # Similarity to reference if available sim = None if args.smiles: ref_mol = Chem.MolFromSmiles(args.smiles) fp1 = AllChem.GetMorganFingerprintAsBitVect(ref_mol, 2, nBits=2048) fp2 = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=2048) sim = round(DataStructs.TanimotoSimilarity(fp1, fp2), 4) evaluated.append({ 'smiles': Chem.MolToSmiles(mol), 'properties': props, 'scores': scores, 'total_gap': round(sum(scores.values()), 4), 'similarity': sim, 'rationale': entry.get('rationale', ''), 'is_reference': entry.get('is_reference', False), }) # Find Pareto front pareto, dominated = find_pareto_front(evaluated) # Sort Pareto front by total gap pareto.sort(key=lambda x: x['total_gap']) result = { 'objectives': objectives, 'total_evaluated': len(evaluated), 'pareto_front_size': len(pareto), 'pareto_front': pareto, 'dominated': dominated[:10], # Top 10 dominated for reference } # Print summary print(f"\n{'='*70}") print(f"PARETO FRONT ({len(pareto)} solutions)") print(f"{'='*70}") for i, p in enumerate(pareto): ref_tag = ' [REFERENCE]' if p.get('is_reference') else '' sim_str = f" sim={p['similarity']}" if p['similarity'] is not None else '' print(f"\n#{i+1}{ref_tag}: {p['smiles']}{sim_str}") print(f" Properties: {p['properties']}") print(f" Gaps: {p['scores']} (total: {p['total_gap']})") if p.get('rationale'): print(f" Modification: {p['rationale']}") with open(args.output, 'w') as f: json.dump(result, f, indent=2) print(f"\nResults saved to {args.output}") return 0 if __name__ == '__main__': sys.exit(main() or 0)
#!/usr/bin/env python3
"""
Pareto-Aware Multi-Objective Molecular Optimization
Generates candidates and identifies the Pareto front across multiple property objectives.
Based on MultiMol (Yu 2025) and MOLLM (Ran 2025).
Usage:
python pareto_optimize.py --smiles "CCO" --objectives "LogP:minimize:3.0,QED:maximize:0.5" --output results.json
python pareto_optimize.py --input candidates.csv --objectives "LogP:minimize:3.0,QED:maximize:0.5" --mode analyze
"""
import argparse
import json
import sys
try:
from rdkit import Chem
from rdkit.Chem import Descriptors, QED, AllChem, DataStructs
from rdkit.Chem.Scaffolds import MurckoScaffold
from rdkit import RDLogger
RDLogger.logger().setLevel(RDLogger.ERROR)
except ImportError:
print("Error: rdkit-pypi required.", file=sys.stderr)
sys.exit(1)
import numpy as np
PROPERTY_FUNCS = {
'MW': lambda m: Descriptors.MolWt(m),
'LogP': lambda m: Descriptors.MolLogP(m),
'TPSA': lambda m: Descriptors.TPSA(m),
'HBA': lambda m: float(Descriptors.NumHAcceptors(m)),
'HBD': lambda m: float(Descriptors.NumHDonors(m)),
'RotBonds': lambda m: float(Descriptors.NumRotatableBonds(m)),
'QED': lambda m: QED.qed(m),
'AromaticRings': lambda m: float(Descriptors.NumAromaticRings(m)),
'FractionCSP3': lambda m: Descriptors.FractionCSP3(m),
'NumRings': lambda m: float(Descriptors.RingCount(m)),
}
BIOISOSTERES = [
('[cH]', '[n]', 'C→N in ring'),
('Cl', 'F', 'Cl→F'),
('OC', 'O', 'OMe→OH'),
('CC', 'C', 'shorten chain'),
('c1ccccc1', 'c1ccncc1', 'phenyl→pyridine'),
('C(=O)N', 'S(=O)(=O)N', 'amide→sulfonamide'),
]
def parse_objectives(obj_str):
"""Parse objective string like 'LogP:minimize:3.0,QED:maximize:0.5'."""
objectives = []
for item in obj_str.split(','):
parts = item.strip().split(':')
if len(parts) == 3:
prop, direction, target = parts
objectives.append({
'property': prop.strip(),
'direction': direction.strip(),
'target': float(target.strip()),
})
elif len(parts) == 4: # range: TPSA:range:20:130
prop, _, lo, hi = parts
objectives.append({
'property': prop.strip(),
'direction': 'range',
'target_min': float(lo.strip()),
'target_max': float(hi.strip()),
})
return objectives
def evaluate_molecule(mol, objectives):
"""Evaluate a molecule against all objectives. Returns scores (lower = better)."""
scores = {}
props = {}
for obj in objectives:
prop = obj['property']
if prop not in PROPERTY_FUNCS:
continue
val = PROPERTY_FUNCS[prop](mol)
props[prop] = round(val, 4)
if obj['direction'] == 'minimize':
score = max(0, val - obj['target'])
elif obj['direction'] == 'maximize':
score = max(0, obj['target'] - val)
elif obj['direction'] == 'range':
if val < obj['target_min']:
score = obj['target_min'] - val
elif val > obj['target_max']:
score = val - obj['target_max']
else:
score = 0
else:
score = abs(val - obj.get('target', val))
scores[prop] = round(score, 4)
return props, scores
def dominates(scores_a, scores_b):
"""Check if solution A Pareto-dominates solution B (all scores lower = better)."""
a_vals = list(scores_a.values())
b_vals = list(scores_b.values())
at_least_one_better = False
for a, b in zip(a_vals, b_vals):
if a > b:
return False
if a < b:
at_least_one_better = True
return at_least_one_better
def find_pareto_front(candidates):
"""Identify the Pareto front from scored candidates."""
pareto = []
dominated = []
for i, cand_i in enumerate(candidates):
is_dominated = False
for j, cand_j in enumerate(candidates):
if i == j:
continue
if dominates(cand_j['scores'], cand_i['scores']):
is_dominated = True
break
if is_dominated:
dominated.append(cand_i)
else:
pareto.append(cand_i)
return pareto, dominated
def generate_candidates(smiles, max_candidates=16):
"""Generate candidate molecules via bioisosteric replacements."""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return []
candidates = []
for old_pat, new_pat, rationale in BIOISOSTERES:
if len(candidates) >= max_candidates:
break
try:
old_mol = Chem.MolFromSmarts(old_pat) or Chem.MolFromSmiles(old_pat)
new_mol = Chem.MolFromSmiles(new_pat) or Chem.MolFromSmarts(new_pat)
if old_mol is None or new_mol is None:
continue
if mol.HasSubstructMatch(old_mol):
products = AllChem.ReplaceSubstructs(mol, old_mol, new_mol)
for prod in products:
try:
smi = Chem.MolToSmiles(prod)
if Chem.MolFromSmiles(smi) is not None and smi != smiles:
candidates.append({'smiles': smi, 'rationale': rationale})
except Exception:
continue
except Exception:
continue
# Deduplicate
seen = set()
unique = []
for c in candidates:
if c['smiles'] not in seen:
seen.add(c['smiles'])
unique.append(c)
return unique[:max_candidates]
def main():
parser = argparse.ArgumentParser(description='Multi-objective molecular optimization')
parser.add_argument('--smiles', type=str, help='Input SMILES to optimize')
parser.add_argument('--input', type=str, help='CSV with candidate molecules')
parser.add_argument('--smiles-col', type=str, default='SMILES')
parser.add_argument('--objectives', type=str, required=True, help='Objectives (e.g., "LogP:minimize:3.0,QED:maximize:0.5")')
parser.add_argument('--candidates', type=int, default=16, help='Max candidates to generate')
parser.add_argument('--mode', choices=['optimize', 'analyze'], default='optimize')
parser.add_argument('--output', type=str, default='pareto_results.json')
args = parser.parse_args()
objectives = parse_objectives(args.objectives)
if not objectives:
print("Error: No valid objectives parsed.", file=sys.stderr)
return 1
molecules = []
if args.smiles and args.mode == 'optimize':
ref_mol = Chem.MolFromSmiles(args.smiles)
if ref_mol is None:
print(f"Error: Invalid SMILES: {args.smiles}", file=sys.stderr)
return 1
ref_props, ref_scores = evaluate_molecule(ref_mol, objectives)
print(f"Reference: {Chem.MolToSmiles(ref_mol)}")
print(f"Properties: {ref_props}")
print(f"Objective gaps: {ref_scores}")
candidates = generate_candidates(args.smiles, args.candidates)
molecules = [{'smiles': args.smiles, 'is_reference': True}] + candidates
elif args.input:
import pandas as pd
df = pd.read_csv(args.input)
molecules = [{'smiles': smi} for smi in df[args.smiles_col]]
else:
print("Error: Provide --smiles or --input", file=sys.stderr)
return 1
# Evaluate all molecules
evaluated = []
for entry in molecules:
mol = Chem.MolFromSmiles(entry['smiles'])
if mol is None:
continue
props, scores = evaluate_molecule(mol, objectives)
# Similarity to reference if available
sim = None
if args.smiles:
ref_mol = Chem.MolFromSmiles(args.smiles)
fp1 = AllChem.GetMorganFingerprintAsBitVect(ref_mol, 2, nBits=2048)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=2048)
sim = round(DataStructs.TanimotoSimilarity(fp1, fp2), 4)
evaluated.append({
'smiles': Chem.MolToSmiles(mol),
'properties': props,
'scores': scores,
'total_gap': round(sum(scores.values()), 4),
'similarity': sim,
'rationale': entry.get('rationale', ''),
'is_reference': entry.get('is_reference', False),
})
# Find Pareto front
pareto, dominated = find_pareto_front(evaluated)
# Sort Pareto front by total gap
pareto.sort(key=lambda x: x['total_gap'])
result = {
'objectives': objectives,
'total_evaluated': len(evaluated),
'pareto_front_size': len(pareto),
'pareto_front': pareto,
'dominated': dominated[:10], # Top 10 dominated for reference
}
# Print summary
print(f"\n{'='*70}")
print(f"PARETO FRONT ({len(pareto)} solutions)")
print(f"{'='*70}")
for i, p in enumerate(pareto):
ref_tag = ' [REFERENCE]' if p.get('is_reference') else ''
sim_str = f" sim={p['similarity']}" if p['similarity'] is not None else ''
print(f"\n#{i+1}{ref_tag}: {p['smiles']}{sim_str}")
print(f" Properties: {p['properties']}")
print(f" Gaps: {p['scores']} (total: {p['total_gap']})")
if p.get('rationale'):
print(f" Modification: {p['rationale']}")
with open(args.output, 'w') as f:
json.dump(result, f, indent=2)
print(f"\nResults saved to {args.output}")
return 0
if __name__ == '__main__':
sys.exit(main() or 0)