Skills › Scientific Agent Skills › rdkit › molecular_properties.py
molecular_properties.py
235 lines · skills/rdkit/scripts/molecular_properties.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Molecular Properties Calculator
5Calculate comprehensive molecular properties and descriptors for molecules.
6Supports single molecules or batch processing from files.
8Usage:
9 python molecular_properties.py "CCO"
10 python molecular_properties.py --file molecules.smi --output properties.csv
11"""
13import argparse
14import sys
16try:
17 from rdkit import Chem, rdBase
18 from rdkit.Chem import Descriptors, Lipinski
19except ImportError:
20 print("Error: RDKit not installed. Install with: uv pip install rdkit")
21 sys.exit(1)
24from _common import read_molecule_records, molecule_smiles
27def calculate_properties(mol):
28 """Calculate comprehensive molecular properties."""
29 if mol is None or mol.GetNumAtoms() == 0:
30 return None
32 properties = {
33 # Basic properties
34 'SMILES': molecule_smiles(mol),
35 'RDKit_Version': rdBase.rdkitVersion,
36 'Molecular_Formula': Chem.rdMolDescriptors.CalcMolFormula(mol),
38 # Molecular weight
39 'MW': Descriptors.MolWt(mol),
40 'ExactMW': Descriptors.ExactMolWt(mol),
42 # Lipophilicity
43 'LogP': Descriptors.MolLogP(mol),
44 'MR': Descriptors.MolMR(mol),
46 # Polar surface area
47 'TPSA': Descriptors.TPSA(mol),
48 'LabuteASA': Descriptors.LabuteASA(mol),
50 # Hydrogen bonding
51 'HBD': Descriptors.NumHDonors(mol),
52 'HBA': Descriptors.NumHAcceptors(mol),
54 # Atom counts
55 'Heavy_Atoms': Descriptors.HeavyAtomCount(mol),
56 'Heteroatoms': Descriptors.NumHeteroatoms(mol),
57 'Valence_Electrons': Descriptors.NumValenceElectrons(mol),
59 # Ring information
60 'Rings': Descriptors.RingCount(mol),
61 'Aromatic_Rings': Descriptors.NumAromaticRings(mol),
62 'Saturated_Rings': Descriptors.NumSaturatedRings(mol),
63 'Aliphatic_Rings': Descriptors.NumAliphaticRings(mol),
64 'Aromatic_Heterocycles': Descriptors.NumAromaticHeterocycles(mol),
66 # Flexibility
67 'Rotatable_Bonds': Descriptors.NumRotatableBonds(mol),
68 'Fraction_Csp3': Lipinski.FractionCSP3(mol),
70 # Complexity
71 'BertzCT': Descriptors.BertzCT(mol),
73 # Drug-likeness
74 'QED': Descriptors.qed(mol),
75 }
77 # Lipinski's Rule of Five
78 properties['Lipinski_Pass'] = (
79 properties['MW'] <= 500 and
80 properties['LogP'] <= 5 and
81 properties['HBD'] <= 5 and
82 properties['HBA'] <= 10
83 )
85 # Lead-likeness
86 properties['Lead-like'] = (
87 250 <= properties['MW'] <= 350 and
88 properties['LogP'] <= 3.5 and
89 properties['Rotatable_Bonds'] <= 7
90 )
92 return properties
95def process_single_molecule(smiles):
96 """Process a single SMILES string."""
97 mol = Chem.MolFromSmiles(smiles)
98 if mol is None or mol.GetNumAtoms() == 0:
99 print(f"Error: Failed to parse SMILES: {smiles}")
100 return None
102 props = calculate_properties(mol)
103 return props
106def process_file(input_file, output_file=None):
107 """Process molecules from a file."""
108 results = []
109 for index, mol in read_molecule_records(input_file):
110 props = calculate_properties(mol)
111 props['Index'] = index
112 props['Name'] = mol.GetProp('_Name') if mol.HasProp('_Name') else ''
113 results.append(props)
114 if not results:
115 raise ValueError("No valid, nonempty molecules in input")
117 # Output results
118 if output_file:
119 write_csv(results, output_file)
120 print(f"Results written to: {output_file}")
121 else:
122 # Print to console
123 for props in results:
124 print("\n" + "="*60)
125 for key, value in props.items():
126 print(f"{key:25s}: {value}")
128 return results
131def write_csv(results, output_file):
132 """Write results to CSV file."""
133 import csv
135 if not results:
136 print("No results to write")
137 return
139 with open(output_file, 'w', newline='') as f:
140 fieldnames = results[0].keys()
141 writer = csv.DictWriter(f, fieldnames=fieldnames)
142 writer.writeheader()
143 writer.writerows(results)
146def print_properties(props):
147 """Print properties in formatted output."""
148 print("\nMolecular Properties:")
149 print("="*60)
151 # Group related properties
152 print("\n[Basic Information]")
153 print(f" SMILES: {props['SMILES']}")
154 print(f" Formula: {props['Molecular_Formula']}")
156 print("\n[Size & Weight]")
157 print(f" Molecular Weight: {props['MW']:.2f}")
158 print(f" Exact MW: {props['ExactMW']:.4f}")
159 print(f" Heavy Atoms: {props['Heavy_Atoms']}")
160 print(f" Heteroatoms: {props['Heteroatoms']}")
162 print("\n[Lipophilicity]")
163 print(f" LogP: {props['LogP']:.2f}")
164 print(f" Molar Refractivity: {props['MR']:.2f}")
166 print("\n[Polarity]")
167 print(f" TPSA: {props['TPSA']:.2f}")
168 print(f" Labute ASA: {props['LabuteASA']:.2f}")
169 print(f" H-bond Donors: {props['HBD']}")
170 print(f" H-bond Acceptors: {props['HBA']}")
172 print("\n[Ring Systems]")
173 print(f" Total Rings: {props['Rings']}")
174 print(f" Aromatic Rings: {props['Aromatic_Rings']}")
175 print(f" Saturated Rings: {props['Saturated_Rings']}")
176 print(f" Aliphatic Rings: {props['Aliphatic_Rings']}")
177 print(f" Aromatic Heterocycles: {props['Aromatic_Heterocycles']}")
179 print("\n[Flexibility & Complexity]")
180 print(f" Rotatable Bonds: {props['Rotatable_Bonds']}")
181 print(f" Fraction Csp3: {props['Fraction_Csp3']:.3f}")
182 print(f" Bertz Complexity: {props['BertzCT']:.1f}")
184 print("\n[Drug-likeness]")
185 print(f" QED Score: {props['QED']:.3f}")
186 print(f" Lipinski Pass: {'Yes' if props['Lipinski_Pass'] else 'No'}")
187 print(f" Lead-like: {'Yes' if props['Lead-like'] else 'No'}")
188 print("="*60)
191def main():
192 parser = argparse.ArgumentParser(
193 description='Calculate molecular properties for molecules',
194 formatter_class=argparse.RawDescriptionHelpFormatter,
195 epilog="""
196Examples:
197 # Single molecule
198 python molecular_properties.py "CCO"
200 # From file
201 python molecular_properties.py --file molecules.smi
203 # Save to CSV
204 python molecular_properties.py --file molecules.sdf --output properties.csv
205 """
206 )
208 parser.add_argument('smiles', nargs='?', help='SMILES string to analyze')
209 parser.add_argument('--file', '-f', help='Input file (SDF or SMILES)')
210 parser.add_argument('--output', '-o', help='Output CSV file')
212 args = parser.parse_args()
214 if not args.smiles and not args.file:
215 parser.print_help()
216 sys.exit(1)
218 if args.smiles is not None:
219 # Process single molecule
220 props = process_single_molecule(args.smiles)
221 if props is None:
222 parser.error("Invalid or empty SMILES")
223 print_properties(props)
224 if args.output:
225 write_csv([props], args.output)
226 elif args.file:
227 # Process file
228 try:
229 process_file(args.file, args.output)
230 except (OSError, ValueError) as error:
231 parser.error(str(error))
234if __name__ == '__main__':
235 main()