Chapter 5: Python in Practice: Polymer Data Analysis Workflow

5.1 Polymer Structure Generation with RDKit

RDKit is an open-source cheminformatics library that enables generation, visualization, and descriptor calculation of molecular structures. For polymer materials, monomer structures are expressed in SMILES notation and used to model repeat units.

5.1.1 Polymer SMILES Notation and Structure Generation

SMILES (Simplified Molecular Input Line Entry System) is a notation that represents molecular structures as strings. For polymers, repeat units are expressed with [*].

from rdkit import Chem
from rdkit.Chem import Draw, Descriptors, AllChem
import matplotlib.pyplot as plt
import numpy as np

# Polymer structure generation with RDKit
def generate_polymer_structure(monomer_smiles, repeat_units=5, polymer_name="Polymer"):
    """
    Generate polymer structure from SMILES notation

    Parameters:
    - monomer_smiles: Monomer SMILES notation ([*] specifies attachment points)
    - repeat_units: Number of repeat units
    - polymer_name: Polymer name

    Returns:
    - polymer_mol: RDKit molecule object
    - properties: Dictionary of molecular descriptors
    """
    print(f"=== {polymer_name} Structure Generation ===")
    print(f"Monomer SMILES: {monomer_smiles}")

    # Generate monomer structure
    monomer = Chem.MolFromSmiles(monomer_smiles)
    if monomer is None:
        print("Error: Invalid SMILES notation")
        return None, {}

    # Generate polymer SMILES (simple repetition)
    # Actual polymerization requires handling the attachment points [*]
    # For simplicity, we simulate linear concatenation here

    # Examples of common polymer SMILES
    polymer_smiles_examples = {
        "Polyethylene": "CC" * repeat_units,
        "Polystyrene": "CC(c1ccccc1)" * repeat_units,
        "PMMA": "CC(C)(C(=O)OC)" * repeat_units,
        "Nylon-6": "N(CCCCC)C(=O)" * repeat_units,
        "PET": "O=C(c1ccc(cc1)C(=O)O)OCCOC" * repeat_units
    }

    # From monomer to polymer (known SMILES are used here for illustration)
    if polymer_name in polymer_smiles_examples:
        polymer_smiles = polymer_smiles_examples[polymer_name]
    else:
        # Simple concatenation for custom monomers
        polymer_smiles = monomer_smiles.replace("[*]", "") * repeat_units

    # Generate polymer molecule object
    polymer_mol = Chem.MolFromSmiles(polymer_smiles)

    if polymer_mol is None:
        print("Error: Polymer SMILES generation failed")
        return None, {}

    # 3D structure optimization
    polymer_mol_3d = Chem.AddHs(polymer_mol)
    AllChem.EmbedMolecule(polymer_mol_3d, randomSeed=42)
    AllChem.UFFOptimizeMolecule(polymer_mol_3d)

    # Calculate molecular descriptors
    properties = {
        "Molecular Weight": Descriptors.MolWt(polymer_mol),
        "LogP": Descriptors.MolLogP(polymer_mol),
        "TPSA": Descriptors.TPSA(polymer_mol),
        "Rotatable Bonds": Descriptors.NumRotatableBonds(polymer_mol),
        "H-Bond Donors": Descriptors.NumHDonors(polymer_mol),
        "H-Bond Acceptors": Descriptors.NumHAcceptors(polymer_mol),
        "Heavy Atoms": Descriptors.HeavyAtomCount(polymer_mol)
    }

    print(f"\nNumber of repeat units: {repeat_units}")
    print("Molecular descriptors:")
    for key, value in properties.items():
        print(f"  {key}: {value:.2f}")

    # Structure visualization
    img = Draw.MolToImage(polymer_mol, size=(600, 400))

    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.axis('off')
    plt.title(f'{polymer_name} Structure (n={repeat_units})', fontsize=14, fontweight='bold')

    # Descriptor bar chart
    plt.subplot(1, 2, 2)
    descriptor_names = ['MW\n(Da)', 'LogP', 'TPSA\n(Ų)', 'Rot.\nBonds', 'HBD', 'HBA']
    descriptor_values = [
        properties["Molecular Weight"] / 100,  # Scale adjustment
        properties["LogP"],
        properties["TPSA"] / 10,  # Scale adjustment
        properties["Rotatable Bonds"],
        properties["H-Bond Donors"],
        properties["H-Bond Acceptors"]
    ]

    colors = ['#f093fb', '#f5576c', '#4A90E2', '#27ae60', '#f39c12', '#e74c3c']
    bars = plt.bar(descriptor_names, descriptor_values, color=colors, edgecolor='black', linewidth=2)
    plt.ylabel('Value (scaled)', fontsize=12)
    plt.title('Molecular Descriptors', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, axis='y')

    for bar, val in zip(bars, descriptor_values):
        height = bar.get_height()
        plt.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.1f}', ha='center', va='bottom', fontsize=9, fontweight='bold')

    plt.tight_layout()
    plt.savefig(f'{polymer_name}_structure.png', dpi=300, bbox_inches='tight')
    plt.show()

    return polymer_mol, properties

# Example 1: Polystyrene
generate_polymer_structure("CC(c1ccccc1)", repeat_units=5, polymer_name="Polystyrene")

# Example 2: PMMA
generate_polymer_structure("CC(C)(C(=O)OC)", repeat_units=5, polymer_name="PMMA")

5.1.2 Morgan Fingerprints and Structural Descriptors

The Morgan Fingerprint is a technique that represents molecular structures as fixed-length bit vectors, used for structural similarity evaluation and as features for machine learning.

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Morgan fingerprint calculation
def calculate_morgan_fingerprint(smiles_list, radius=2, n_bits=2048):
    """
    Calculate Morgan fingerprints and generate a similarity matrix

    Parameters:
    - smiles_list: List of SMILES notations (dictionary form {name: SMILES})
    - radius: Morgan fingerprint radius (default: 2)
    - n_bits: Bit length (default: 2048)

    Returns:
    - fingerprints: List of fingerprints
    - similarity_matrix: Tanimoto similarity matrix
    """
    print("=== Morgan Fingerprint Calculation ===")

    # Generate molecule objects
    mols = {}
    fingerprints = {}

    for name, smiles in smiles_list.items():
        mol = Chem.MolFromSmiles(smiles)
        if mol is not None:
            mols[name] = mol
            # Generate Morgan fingerprint
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=n_bits)
            fingerprints[name] = fp
            print(f"{name}: Fingerprint generated ({n_bits} bits)")
        else:
            print(f"Error: SMILES conversion failed for {name}")

    # Calculate Tanimoto similarity matrix
    polymer_names = list(fingerprints.keys())
    n = len(polymer_names)
    similarity_matrix = np.zeros((n, n))

    for i, name1 in enumerate(polymer_names):
        for j, name2 in enumerate(polymer_names):
            fp1 = fingerprints[name1]
            fp2 = fingerprints[name2]
            similarity = DataStructs.TanimotoSimilarity(fp1, fp2)
            similarity_matrix[i, j] = similarity

    # Visualization
    plt.figure(figsize=(14, 6))

    # Subplot 1: Similarity heatmap
    plt.subplot(1, 2, 1)
    sns.heatmap(similarity_matrix, annot=True, fmt='.3f', cmap='RdYlGn',
                xticklabels=polymer_names, yticklabels=polymer_names,
                vmin=0, vmax=1, cbar_kws={'label': 'Tanimoto Similarity'})
    plt.title('Morgan Fingerprint Similarity Matrix', fontsize=14, fontweight='bold')

    # Subplot 2: Fingerprint bit distribution
    plt.subplot(1, 2, 2)
    for name, fp in fingerprints.items():
        # Convert fingerprint to numpy array
        arr = np.zeros((n_bits,))
        DataStructs.ConvertToNumpyArray(fp, arr)
        # Plot the number of set bits
        bit_count = np.sum(arr)
        plt.bar(name, bit_count, edgecolor='black', linewidth=2)

    plt.ylabel('Number of Set Bits', fontsize=12)
    plt.title('Fingerprint Bit Density', fontsize=14, fontweight='bold')
    plt.xticks(rotation=45, ha='right')
    plt.grid(alpha=0.3, axis='y')

    plt.tight_layout()
    plt.savefig('morgan_fingerprint_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Output results
    print("\n=== Tanimoto Similarity Matrix ===")
    for i, name1 in enumerate(polymer_names):
        for j, name2 in enumerate(polymer_names):
            if i < j:  # Show upper triangle only
                print(f"{name1} - {name2}: {similarity_matrix[i, j]:.3f}")

    return fingerprints, similarity_matrix

# Example: structural similarity of representative polymers
polymer_smiles = {
    "Polyethylene": "CC",
    "Polypropylene": "CC(C)",
    "Polystyrene": "CC(c1ccccc1)",
    "PMMA": "CC(C)(C(=O)OC)",
    "PVC": "CC(Cl)",
    "PVDF": "CC(F)(F)"
}

calculate_morgan_fingerprint(polymer_smiles, radius=2, n_bits=2048)

5.2 Property Prediction with Machine Learning

We predict properties of polymer materials such as the glass transition temperature (T_g) and density using machine learning models. Here, we build a regression model using Random Forest.

5.2.1 T_g Prediction Model (scikit-learn Random Forest)

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import matplotlib.pyplot as plt

# Build Tg prediction model
def build_tg_prediction_model():
    """
    Polymer Tg prediction model using Random Forest

    Returns:
    - model: Trained model
    - X_test, y_test: Test data
    - metrics: Dictionary of evaluation metrics
    """
    print("=== Building Tg Prediction Model (Random Forest) ===")

    # Sample dataset (in practice, obtained from PolyInfo etc.)
    # Features: [MW, LogP, TPSA, RotBonds, HBD, HBA, Aromatic_Ratio]
    data = {
        'Polymer': ['PS', 'PMMA', 'PC', 'PET', 'Nylon-6', 'PE', 'PP', 'PVC',
                    'PVDF', 'PTFE', 'PAN', 'PVA', 'Cellulose', 'PLA', 'PEEK'],
        'MW': [104, 100, 254, 192, 113, 28, 42, 62.5, 64, 100, 53, 44, 162, 72, 288],
        'LogP': [3.2, 1.5, 2.8, 1.2, -0.5, 1.9, 2.1, 1.4, 1.0, 3.5, -0.3, -0.8, -1.2, 0.5, 3.8],
        'TPSA': [0, 26.3, 40.5, 52.6, 43.1, 0, 0, 0, 0, 0, 23.8, 20.2, 90.2, 26.3, 25.8],
        'RotBonds': [2, 3, 8, 6, 10, 1, 2, 1, 1, 0, 1, 1, 4, 2, 6],
        'HBD': [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0],
        'HBA': [0, 2, 2, 4, 1, 0, 0, 0, 2, 0, 1, 1, 5, 2, 2],
        'Aromatic_Ratio': [100, 0, 80, 60, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 90],
        'Tg_K': [373, 378, 423, 343, 323, 195, 253, 354, 233, 115, 358, 358, 503, 333, 416]
    }

    df = pd.DataFrame(data)

    # Features and target
    feature_cols = ['MW', 'LogP', 'TPSA', 'RotBonds', 'HBD', 'HBA', 'Aromatic_Ratio']
    X = df[feature_cols].values
    y = df['Tg_K'].values

    # Data split (train:test = 80:20)
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42
    )

    # Build Random Forest model
    model = RandomForestRegressor(
        n_estimators=100,
        max_depth=10,
        min_samples_split=2,
        random_state=42
    )

    # Training
    model.fit(X_train, y_train)

    # Prediction
    y_train_pred = model.predict(X_train)
    y_test_pred = model.predict(X_test)

    # Evaluation metrics
    train_r2 = r2_score(y_train, y_train_pred)
    test_r2 = r2_score(y_test, y_test_pred)
    train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
    test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
    test_mae = mean_absolute_error(y_test, y_test_pred)

    # Cross-validation
    cv_scores = cross_val_score(model, X, y, cv=5, scoring='r2')

    metrics = {
        'Train R²': train_r2,
        'Test R²': test_r2,
        'Train RMSE': train_rmse,
        'Test RMSE': test_rmse,
        'Test MAE': test_mae,
        'CV R² Mean': cv_scores.mean(),
        'CV R² Std': cv_scores.std()
    }

    # Visualization
    plt.figure(figsize=(14, 5))

    # Subplot 1: Predicted vs actual
    plt.subplot(1, 3, 1)
    plt.scatter(y_train, y_train_pred, alpha=0.6, s=100, label='Train', edgecolors='black')
    plt.scatter(y_test, y_test_pred, alpha=0.8, s=150, label='Test',
                edgecolors='black', linewidths=2)

    # Diagonal (perfect prediction)
    min_val = min(y.min(), y_train_pred.min(), y_test_pred.min())
    max_val = max(y.max(), y_train_pred.max(), y_test_pred.max())
    plt.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Prediction')

    plt.xlabel('Actual Tg (K)', fontsize=12)
    plt.ylabel('Predicted Tg (K)', fontsize=12)
    plt.title(f'Tg Prediction (R² = {test_r2:.3f})', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 2: Residual plot
    plt.subplot(1, 3, 2)
    residuals = y_test - y_test_pred
    plt.scatter(y_test_pred, residuals, s=150, alpha=0.6, edgecolors='black', linewidths=2)
    plt.axhline(0, color='red', linestyle='--', linewidth=2)
    plt.xlabel('Predicted Tg (K)', fontsize=12)
    plt.ylabel('Residuals (K)', fontsize=12)
    plt.title('Residual Plot', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3)

    # Subplot 3: Evaluation metrics
    plt.subplot(1, 3, 3)
    metric_names = ['R² Train', 'R² Test', 'RMSE\nTrain', 'RMSE\nTest', 'MAE\nTest']
    metric_values = [train_r2, test_r2, train_rmse/100, test_rmse/100, test_mae/100]  # Scale adjustment

    colors = ['#27ae60', '#3498db', '#f39c12', '#e74c3c', '#9b59b6']
    bars = plt.bar(metric_names, metric_values, color=colors, edgecolor='black', linewidth=2)
    plt.ylabel('Value (R² or RMSE/100)', fontsize=12)
    plt.title('Model Performance Metrics', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, axis='y')

    for bar, val, orig in zip(bars, metric_values, [train_r2, test_r2, train_rmse, test_rmse, test_mae]):
        height = bar.get_height()
        if 'RMSE' in metric_names[bars.index(bar)] or 'MAE' in metric_names[bars.index(bar)]:
            label_text = f'{orig:.1f}K'
        else:
            label_text = f'{orig:.3f}'
        plt.text(bar.get_x() + bar.get_width()/2., height,
                label_text, ha='center', va='bottom', fontsize=9, fontweight='bold')

    plt.tight_layout()
    plt.savefig('tg_prediction_model.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Output results
    print("\n=== Model Evaluation Metrics ===")
    for key, value in metrics.items():
        print(f"{key}: {value:.4f}")

    print("\n=== Test Set Prediction Results ===")
    test_indices = df.index[len(y_train):]
    for idx, (actual, pred) in zip(test_indices, zip(y_test, y_test_pred)):
        polymer_name = df.loc[idx, 'Polymer']
        print(f"{polymer_name}: Actual {actual:.1f}K, Predicted {pred:.1f}K, Error {abs(actual-pred):.1f}K")

    return model, X_test, y_test, metrics, feature_cols

# Execute
model, X_test, y_test, metrics, feature_cols = build_tg_prediction_model()

5.2.2 Feature Importance and Model Interpretation

import numpy as np
import matplotlib.pyplot as plt
from sklearn.inspection import permutation_importance

# Feature importance analysis
def analyze_feature_importance(model, X_test, y_test, feature_names):
    """
    Analyze Random Forest feature importance

    Parameters:
    - model: Trained model
    - X_test: Test data features
    - y_test: Test data target
    - feature_names: List of feature names

    Returns:
    - importances: Dictionary of feature importances
    """
    print("=== Feature Importance Analysis ===")

    # Gini importance (built-in)
    gini_importances = model.feature_importances_

    # Permutation importance (more reliable)
    perm_importance = permutation_importance(
        model, X_test, y_test, n_repeats=10, random_state=42
    )
    perm_importances = perm_importance.importances_mean

    # Visualization
    plt.figure(figsize=(14, 6))

    # Subplot 1: Gini importance
    plt.subplot(1, 2, 1)
    sorted_idx = np.argsort(gini_importances)[::-1]

    colors = plt.cm.RdYlGn(np.linspace(0.3, 0.9, len(feature_names)))
    bars = plt.barh(range(len(feature_names)),
                    gini_importances[sorted_idx],
                    color=colors, edgecolor='black', linewidth=2)
    plt.yticks(range(len(feature_names)), np.array(feature_names)[sorted_idx])
    plt.xlabel('Gini Importance', fontsize=12)
    plt.title('Feature Importance (Gini)', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, axis='x')

    for bar, val in zip(bars, gini_importances[sorted_idx]):
        width = bar.get_width()
        plt.text(width, bar.get_y() + bar.get_height()/2.,
                f'{val:.3f}', ha='left', va='center', fontsize=10, fontweight='bold')

    # Subplot 2: Permutation importance
    plt.subplot(1, 2, 2)
    sorted_idx_perm = np.argsort(perm_importances)[::-1]

    bars2 = plt.barh(range(len(feature_names)),
                     perm_importances[sorted_idx_perm],
                     color=colors, edgecolor='black', linewidth=2)
    plt.yticks(range(len(feature_names)), np.array(feature_names)[sorted_idx_perm])
    plt.xlabel('Permutation Importance', fontsize=12)
    plt.title('Feature Importance (Permutation)', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, axis='x')

    for bar, val in zip(bars2, perm_importances[sorted_idx_perm]):
        width = bar.get_width()
        plt.text(width, bar.get_y() + bar.get_height()/2.,
                f'{val:.3f}', ha='left', va='center', fontsize=10, fontweight='bold')

    plt.tight_layout()
    plt.savefig('feature_importance.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Output results
    print("\n=== Gini Importance Ranking ===")
    for idx in sorted_idx:
        print(f"{feature_names[idx]}: {gini_importances[idx]:.4f}")

    print("\n=== Permutation Importance Ranking ===")
    for idx in sorted_idx_perm:
        print(f"{feature_names[idx]}: {perm_importances[idx]:.4f}")

    importances = {
        'gini': dict(zip(feature_names, gini_importances)),
        'permutation': dict(zip(feature_names, perm_importances))
    }

    return importances

# Execute
analyze_feature_importance(model, X_test, y_test, feature_cols)

5.3 MD Simulation Data Analysis

We analyze trajectory data obtained from molecular dynamics (MD) simulations with MDAnalysis to calculate diffusion coefficients and structural parameters.

5.3.1 Trajectory Analysis with MDAnalysis

import numpy as np
import matplotlib.pyplot as plt

# MD simulation data analysis (simplified)
# With actual MDAnalysis: import MDAnalysis as mda
# Here we simulate the behavior with synthetic data

def analyze_md_trajectory_simplified():
    """
    Simplified analysis of MD simulation trajectories
    (use MDAnalysis in real environments)

    Returns:
    - time: Time (ps)
    - rg: Radius of gyration (Å)
    - end_to_end: End-to-end distance (Å)
    """
    print("=== MD Trajectory Analysis (Simplified) ===")

    # Generate synthetic MD data
    # In practice: u = mda.Universe('topology.pdb', 'trajectory.dcd')
    time = np.linspace(0, 10000, 1000)  # ps

    # Radius of gyration
    # Measure of polymer compactness
    rg_mean = 15.0  # Å
    rg = rg_mean + np.random.normal(0, 1.5, len(time))

    # End-to-end distance
    # Measure of polymer chain extension
    ete_mean = 40.0  # Å
    ete = ete_mean + np.random.normal(0, 5.0, len(time))

    # Visualization
    plt.figure(figsize=(14, 10))

    # Subplot 1: Rg vs time
    plt.subplot(3, 2, 1)
    plt.plot(time, rg, linewidth=1, alpha=0.7)
    plt.axhline(rg.mean(), color='red', linestyle='--', linewidth=2,
                label=f'Mean = {rg.mean():.2f} Å')
    plt.xlabel('Time (ps)', fontsize=12)
    plt.ylabel('Radius of Gyration Rg (Å)', fontsize=12)
    plt.title('Rg vs Time', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 2: Rg distribution
    plt.subplot(3, 2, 2)
    plt.hist(rg, bins=30, edgecolor='black', alpha=0.7, color='#3498db')
    plt.axvline(rg.mean(), color='red', linestyle='--', linewidth=2,
                label=f'Mean = {rg.mean():.2f} Å')
    plt.xlabel('Radius of Gyration Rg (Å)', fontsize=12)
    plt.ylabel('Frequency', fontsize=12)
    plt.title('Rg Distribution', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 3: End-to-end distance vs time
    plt.subplot(3, 2, 3)
    plt.plot(time, ete, linewidth=1, alpha=0.7, color='green')
    plt.axhline(ete.mean(), color='red', linestyle='--', linewidth=2,
                label=f'Mean = {ete.mean():.2f} Å')
    plt.xlabel('Time (ps)', fontsize=12)
    plt.ylabel('End-to-End Distance (Å)', fontsize=12)
    plt.title('End-to-End Distance vs Time', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 4: End-to-end distance distribution
    plt.subplot(3, 2, 4)
    plt.hist(ete, bins=30, edgecolor='black', alpha=0.7, color='#27ae60')
    plt.axvline(ete.mean(), color='red', linestyle='--', linewidth=2,
                label=f'Mean = {ete.mean():.2f} Å')
    plt.xlabel('End-to-End Distance (Å)', fontsize=12)
    plt.ylabel('Frequency', fontsize=12)
    plt.title('End-to-End Distance Distribution', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 5: Rg vs end-to-end correlation
    plt.subplot(3, 2, 5)
    plt.scatter(rg, ete, alpha=0.3, s=20, edgecolors='none')

    # Linear regression
    z = np.polyfit(rg, ete, 1)
    p = np.poly1d(z)
    plt.plot(rg, p(rg), "r--", linewidth=2, label=f'y = {z[0]:.2f}x + {z[1]:.2f}')

    plt.xlabel('Radius of Gyration Rg (Å)', fontsize=12)
    plt.ylabel('End-to-End Distance (Å)', fontsize=12)
    plt.title('Rg vs End-to-End Correlation', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 6: Statistical summary
    plt.subplot(3, 2, 6)
    stats = {
        'Rg Mean (Å)': rg.mean(),
        'Rg Std (Å)': rg.std(),
        'EtE Mean (Å)': ete.mean(),
        'EtE Std (Å)': ete.std(),
        'Correlation': np.corrcoef(rg, ete)[0, 1]
    }

    stat_names = list(stats.keys())
    stat_values = list(stats.values())

    colors = ['#3498db', '#3498db', '#27ae60', '#27ae60', '#9b59b6']
    bars = plt.barh(stat_names, stat_values, color=colors, edgecolor='black', linewidth=2)
    plt.xlabel('Value', fontsize=12)
    plt.title('Statistical Summary', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, axis='x')

    for bar, val in zip(bars, stat_values):
        width = bar.get_width()
        plt.text(width, bar.get_y() + bar.get_height()/2.,
                f'{val:.2f}', ha='left', va='center', fontsize=10, fontweight='bold')

    plt.tight_layout()
    plt.savefig('md_trajectory_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Output results
    print("\n=== Trajectory Statistics ===")
    for key, value in stats.items():
        print(f"{key}: {value:.3f}")

    return time, rg, ete

# Execute
time, rg, ete = analyze_md_trajectory_simplified()

5.3.2 MSD (Mean Square Displacement) Calculation and Diffusion Coefficient

import numpy as np
import matplotlib.pyplot as plt

# MSD calculation and diffusion coefficient derivation
def calculate_msd_diffusion_coefficient(temperature=300):
    """
    Calculate diffusion coefficient from mean square displacement (MSD)

    Parameters:
    - temperature: Temperature (K)

    Returns:
    - time: Time (ps)
    - msd: MSD (Ų)
    - diffusion_coeff: Diffusion coefficient (cm²/s)
    """
    print(f"=== MSD Calculation and Diffusion Coefficient (T = {temperature}K) ===")

    # Generate synthetic data
    # In real MD analysis, this is computed from atomic coordinates
    time = np.linspace(0, 1000, 500)  # ps

    # MSD = 6Dt (3D diffusion)
    # D: diffusion coefficient (Ų/ps)
    D_true = 5e-3  # Ų/ps (assumed value)
    msd = 6 * D_true * time + np.random.normal(0, 0.5, len(time))
    msd = np.maximum(msd, 0)  # Prevent negative values

    # Linear fitting (use only the later linear region)
    # Exclude the initial nonlinear region
    linear_region = time > 200  # ps
    z = np.polyfit(time[linear_region], msd[linear_region], 1)
    slope = z[0]  # Ų/ps

    # Diffusion coefficient D = slope / 6
    D_fitted = slope / 6  # Ų/ps

    # Unit conversion: Ų/ps to cm²/s
    # 1 Ų = 10⁻¹⁶ cm², 1 ps = 10⁻¹² s
    D_cm2_s = D_fitted * 1e-16 / 1e-12  # cm²/s
    D_cm2_s = D_fitted * 1e-4  # cm²/s

    # Visualization
    plt.figure(figsize=(14, 5))

    # Subplot 1: MSD vs Time
    plt.subplot(1, 3, 1)
    plt.plot(time, msd, 'b-', linewidth=2, label='MSD Data')

    # Fitted line
    fit_line = z[0] * time + z[1]
    plt.plot(time, fit_line, 'r--', linewidth=2,
             label=f'Fit: MSD = {z[0]:.4f}t + {z[1]:.2f}')

    # Mark the linear region
    plt.axvline(200, color='green', linestyle=':', linewidth=1.5,
                label='Linear Region Start')

    plt.xlabel('Time (ps)', fontsize=12)
    plt.ylabel('MSD (Ų)', fontsize=12)
    plt.title('Mean Square Displacement', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3)

    # Subplot 2: Log-log plot (power-law check)
    plt.subplot(1, 3, 2)
    plt.loglog(time[time > 0], msd[time > 0], 'b-', linewidth=2)

    # Power law MSD ∝ t^α (α=1: normal diffusion)
    alpha = np.polyfit(np.log(time[linear_region]), np.log(msd[linear_region]), 1)[0]

    plt.xlabel('Time (ps)', fontsize=12)
    plt.ylabel('MSD (Ų)', fontsize=12)
    plt.title(f'Log-Log Plot (α = {alpha:.2f})', fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3, which='both')

    # Subplot 3: Temperature dependence of the diffusion coefficient (Arrhenius)
    plt.subplot(1, 3, 3)
    temperatures = np.linspace(250, 400, 20)  # K

    # Arrhenius equation: D = D0 * exp(-Ea/RT)
    D0 = 1e-2  # cm²/s
    Ea = 30000  # J/mol
    R = 8.314  # J/mol·K

    D_temps = D0 * np.exp(-Ea / (R * temperatures))

    plt.semilogy(temperatures, D_temps, 'purple', linewidth=2, label='Arrhenius Model')
    plt.scatter([temperature], [D_cm2_s], s=200, c='red', edgecolors='black',
                linewidths=2, zorder=5, label=f'Current ({temperature}K)')

    plt.xlabel('Temperature (K)', fontsize=12)
    plt.ylabel('Diffusion Coefficient (cm²/s)', fontsize=12)
    plt.title('Temperature Dependence of D', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(alpha=0.3, which='both')

    plt.tight_layout()
    plt.savefig('msd_diffusion.png', dpi=300, bbox_inches='tight')
    plt.show()

    # Output results
    print(f"\n=== MSD Analysis Results ===")
    print(f"Linear fit: MSD = {z[0]:.4f} * t + {z[1]:.2f}")
    print(f"Diffusion coefficient D: {D_fitted:.2e} Ų/ps")
    print(f"Diffusion coefficient D: {D_cm2_s:.2e} cm²/s")
    print(f"Power exponent α: {alpha:.3f} (α=1: normal diffusion)")

    # Validation via the Einstein relation
    # D = kT / (6πηr) (Stokes-Einstein equation)
    # Here we simply check plausibility
    print(f"\nExperimental diffusion coefficient range: 10⁻⁷ - 10⁻⁵ cm²/s (typical polymer melts)")

    return time, msd, D_cm2_s

# Execute
calculate_msd_diffusion_coefficient(temperature=300)

5.4 Building an Integrated Workflow

We now integrate all the methods learned so far and implement them as the PolymerAnalysis class. This class can handle structure generation, property prediction, and MD analysis in one place.

5.4.1 Implementing the PolymerAnalysis Class

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem import Descriptors, AllChem
from sklearn.ensemble import RandomForestRegressor
import matplotlib.pyplot as plt

# Integrated workflow class
class PolymerAnalysis:
    """
    Integrated analysis workflow for polymer materials

    Features:
    - Structure generation and descriptor calculation with RDKit
    - Tg prediction with machine learning
    - MD data analysis
    - Automated report generation
    """

    def __init__(self, polymer_name, smiles):
        """
        Initialization

        Parameters:
        - polymer_name: Polymer name
        - smiles: SMILES notation
        """
        self.polymer_name = polymer_name
        self.smiles = smiles
        self.mol = Chem.MolFromSmiles(smiles)
        self.descriptors = {}
        self.predicted_tg = None

        print(f"=== PolymerAnalysis Initialization: {polymer_name} ===")
        print(f"SMILES: {smiles}")

    def calculate_descriptors(self):
        """Calculate RDKit descriptors"""
        if self.mol is None:
            print("Error: Invalid molecule object")
            return None

        self.descriptors = {
            'MW': Descriptors.MolWt(self.mol),
            'LogP': Descriptors.MolLogP(self.mol),
            'TPSA': Descriptors.TPSA(self.mol),
            'RotBonds': Descriptors.NumRotatableBonds(self.mol),
            'HBD': Descriptors.NumHDonors(self.mol),
            'HBA': Descriptors.NumHAcceptors(self.mol),
            'AromaticRings': Descriptors.NumAromaticRings(self.mol)
        }

        print("\n=== Molecular Descriptors ===")
        for key, value in self.descriptors.items():
            print(f"{key}: {value:.2f}")

        return self.descriptors

    def predict_tg(self, model=None):
        """
        Tg prediction (using a simplified model)

        Parameters:
        - model: Trained model (simple estimation if None)

        Returns:
        - predicted_tg: Predicted Tg (K)
        """
        if not self.descriptors:
            self.calculate_descriptors()

        # Simple estimation formula (use a trained model in practice)
        # Tg ≈ 250 + 0.5*MW + 10*HBD + 5*HBA - 20*RotBonds
        tg_estimate = (250 +
                      0.5 * self.descriptors['MW'] +
                      10 * self.descriptors['HBD'] +
                      5 * self.descriptors['HBA'] -
                      20 * self.descriptors['RotBonds'])

        self.predicted_tg = tg_estimate

        print(f"\n=== Tg Prediction ===")
        print(f"Predicted Tg: {tg_estimate:.1f} K ({tg_estimate - 273.15:.1f}°C)")

        return tg_estimate

    def analyze_md_data(self, time_range=1000):
        """
        MD analysis (synthetic data)

        Parameters:
        - time_range: Simulation time (ps)

        Returns:
        - md_results: Dictionary of analysis results
        """
        print(f"\n=== MD Analysis ({time_range} ps) ===")

        # Synthetic MD data
        time = np.linspace(0, time_range, 500)
        rg = 15 + np.random.normal(0, 1.5, len(time))
        msd = 6 * 5e-3 * time + np.random.normal(0, 0.5, len(time))

        # Calculate diffusion coefficient
        linear_region = time > time_range * 0.2
        slope = np.polyfit(time[linear_region], msd[linear_region], 1)[0]
        D = slope / 6 * 1e-4  # cm²/s

        md_results = {
            'Rg_mean': rg.mean(),
            'Rg_std': rg.std(),
            'Diffusion_coeff': D
        }

        print(f"Radius of gyration Rg: {md_results['Rg_mean']:.2f} ± {md_results['Rg_std']:.2f} Å")
        print(f"Diffusion coefficient D: {md_results['Diffusion_coeff']:.2e} cm²/s")

        return md_results

    def generate_report(self):
        """Generate integrated report"""
        print(f"\n{'='*60}")
        print(f"Polymer Material Analysis Report: {self.polymer_name}")
        print(f"{'='*60}")

        print(f"\n[Structure Information]")
        print(f"SMILES: {self.smiles}")

        if self.descriptors:
            print(f"\n[Molecular Descriptors]")
            for key, value in self.descriptors.items():
                print(f"  {key}: {value:.2f}")

        if self.predicted_tg:
            print(f"\n[Property Prediction]")
            print(f"  Glass transition temperature Tg: {self.predicted_tg:.1f} K ({self.predicted_tg - 273.15:.1f}°C)")

        print(f"\n[Recommended Applications]")
        if self.predicted_tg and self.predicted_tg > 350:
            print("  - High-temperature engineering plastics")
            print("  - Aerospace materials")
        elif self.predicted_tg and self.predicted_tg > 300:
            print("  - General-purpose engineering plastics")
            print("  - Automotive components")
        else:
            print("  - Commodity plastics")
            print("  - Packaging materials")

        print(f"\n{'='*60}\n")

    def visualize_summary(self):
        """Visualize summary"""
        if not self.descriptors:
            self.calculate_descriptors()

        fig, axes = plt.subplots(2, 2, figsize=(14, 10))

        # Subplot 1: Descriptor radar chart
        ax1 = axes[0, 0]
        categories = list(self.descriptors.keys())[:6]
        values = [self.descriptors[cat] for cat in categories]

        # Normalization
        max_vals = [500, 5, 150, 20, 5, 10]  # Assumed maximum for each descriptor
        values_norm = [v / m for v, m in zip(values, max_vals)]

        angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist()
        values_norm += values_norm[:1]
        angles += angles[:1]

        ax1 = plt.subplot(2, 2, 1, projection='polar')
        ax1.plot(angles, values_norm, 'o-', linewidth=2, color='#f093fb')
        ax1.fill(angles, values_norm, alpha=0.25, color='#f5576c')
        ax1.set_xticks(angles[:-1])
        ax1.set_xticklabels(categories, fontsize=10)
        ax1.set_ylim(0, 1)
        ax1.set_title(f'{self.polymer_name}: Descriptor Profile', fontsize=12, fontweight='bold', pad=20)
        ax1.grid(True)

        # Subplot 2: Chemical structure
        ax2 = axes[0, 1]
        from rdkit.Chem import Draw
        img = Draw.MolToImage(self.mol, size=(400, 300))
        ax2.imshow(img)
        ax2.axis('off')
        ax2.set_title('Chemical Structure', fontsize=12, fontweight='bold')

        # Subplot 3: Tg prediction
        ax3 = axes[1, 0]
        if self.predicted_tg:
            bars = ax3.bar(['Predicted Tg'], [self.predicted_tg],
                          color='#27ae60', edgecolor='black', linewidth=2)
            ax3.axhline(273.15, color='blue', linestyle='--', linewidth=1.5, label='0°C')
            ax3.axhline(373.15, color='red', linestyle='--', linewidth=1.5, label='100°C')
            ax3.set_ylabel('Temperature (K)', fontsize=12)
            ax3.set_title('Glass Transition Temperature Prediction', fontsize=12, fontweight='bold')
            ax3.legend()
            ax3.grid(alpha=0.3, axis='y')

            for bar in bars:
                height = bar.get_height()
                ax3.text(bar.get_x() + bar.get_width()/2., height,
                        f'{height:.1f}K\n({height-273.15:.1f}°C)',
                        ha='center', va='bottom', fontsize=11, fontweight='bold')

        # Subplot 4: Application classification
        ax4 = axes[1, 1]
        applications = ['Engineering\nPlastic', 'General\nPurpose', 'Commodity', 'High-Temp\nSpecialty']
        scores = [0, 0, 0, 0]

        if self.predicted_tg:
            if self.predicted_tg > 400:
                scores[3] = 1
            elif self.predicted_tg > 350:
                scores[0] = 1
            elif self.predicted_tg > 300:
                scores[1] = 1
            else:
                scores[2] = 1

        colors_app = ['#3498db', '#27ae60', '#f39c12', '#e74c3c']
        bars = ax4.barh(applications, scores, color=colors_app, edgecolor='black', linewidth=2)
        ax4.set_xlim(0, 1.2)
        ax4.set_xlabel('Suitability Score', fontsize=12)
        ax4.set_title('Application Classification', fontsize=12, fontweight='bold')
        ax4.grid(alpha=0.3, axis='x')

        plt.tight_layout()
        plt.savefig(f'{self.polymer_name}_analysis_summary.png', dpi=300, bbox_inches='tight')
        plt.show()

# Usage example
polymer = PolymerAnalysis("PMMA", "CC(C)(C(=O)OC)")
polymer.calculate_descriptors()
polymer.predict_tg()
polymer.analyze_md_data(time_range=1000)
polymer.generate_report()
polymer.visualize_summary()

5.4.2 PolyInfo API Database Integration

import requests
import json
import pandas as pd

# PolyInfo API integration (simplified)
def fetch_polyinfo_data(polymer_name, property_type='Tg'):
    """
    Fetch polymer data from the PolyInfo API (mock implementation)

    See the official PolyInfo documentation for the actual API specification:
    https://polymer.nims.go.jp/

    Parameters:
    - polymer_name: Polymer name
    - property_type: Property type ('Tg', 'density', 'modulus')

    Returns:
    - data: Retrieved data (DataFrame)
    """
    print(f"=== PolyInfo API Data Retrieval: {polymer_name} ===")

    # Example of an actual API call (pseudocode)
    # api_url = "https://polymer.nims.go.jp/api/v1/polymers"
    # params = {'name': polymer_name, 'property': property_type}
    # response = requests.get(api_url, params=params)
    # data = response.json()

    # Return mock data here
    sample_data = {
        'Polymer': [polymer_name] * 5,
        'Measurement_Condition': ['DSC-10K/min', 'DSC-20K/min', 'DMA-1Hz', 'DSC-10K/min', 'TMA'],
        'Tg_K': [378, 375, 380, 377, 379],
        'Reference': ['Smith2020', 'Jones2019', 'Lee2021', 'Kim2018', 'Wang2022']
    }

    df = pd.DataFrame(sample_data)

    print(f"\nNumber of records retrieved: {len(df)}")
    print("\nData sample:")
    print(df.head())

    # Statistics
    print(f"\n=== {property_type} Statistics ===")
    print(f"Mean: {df['Tg_K'].mean():.1f} K")
    print(f"Standard deviation: {df['Tg_K'].std():.1f} K")
    print(f"Min-Max: {df['Tg_K'].min():.1f} - {df['Tg_K'].max():.1f} K")

    return df

# Usage example
polyinfo_data = fetch_polyinfo_data("PMMA", property_type='Tg')

5.5 Practical Project Example

Finally, we put into practice a complete workflow from new polymer material design to property optimization.

5.5.1 New Material Design Workflow

import numpy as np
from rdkit import Chem
from rdkit.Chem import Descriptors

# New material design workflow
def design_new_polymer_material(target_tg=423, max_iterations=10):
    """
    Design new polymers based on target specifications

    Parameters:
    - target_tg: Target Tg (K)
    - max_iterations: Maximum number of iterations

    Returns:
    - best_candidate: Best candidate
    """
    print(f"=== New Polymer Material Design ===")
    print(f"Target Tg: {target_tg} K ({target_tg - 273.15}°C)")

    # Candidate monomer list
    monomers = {
        'Styrene': 'CC(c1ccccc1)',
        'MMA': 'CC(C)(C(=O)OC)',
        'Acrylonitrile': 'CC(C#N)',
        'Vinyl_Acetate': 'CC(OC(=O)C)',
        'Butadiene': 'C=CC=C',
        'Isoprene': 'CC(=C)C=C'
    }

    # Functional group library
    functional_groups = {
        'Phenyl': 'c1ccccc1',
        'Methyl': 'C',
        'Cyano': 'C#N',
        'Carbonyl': 'C(=O)',
        'Hydroxyl': 'O'
    }

    candidates = []

    for iteration in range(max_iterations):
        # Randomly combine monomers and functional groups
        monomer_name = np.random.choice(list(monomers.keys()))
        monomer_smiles = monomers[monomer_name]

        # Simple Tg estimation
        mol = Chem.MolFromSmiles(monomer_smiles)
        if mol is None:
            continue

        descriptors = {
            'MW': Descriptors.MolWt(mol),
            'HBD': Descriptors.NumHDonors(mol),
            'HBA': Descriptors.NumHAcceptors(mol),
            'RotBonds': Descriptors.NumRotatableBonds(mol),
            'AromaticRings': Descriptors.NumAromaticRings(mol)
        }

        # Simple Tg estimation formula
        tg_estimate = (250 +
                      0.5 * descriptors['MW'] +
                      10 * descriptors['HBD'] +
                      5 * descriptors['HBA'] -
                      20 * descriptors['RotBonds'] +
                      30 * descriptors['AromaticRings'])

        error = abs(tg_estimate - target_tg)

        candidates.append({
            'Monomer': monomer_name,
            'SMILES': monomer_smiles,
            'Predicted_Tg': tg_estimate,
            'Error': error,
            'Descriptors': descriptors
        })

        print(f"\nIteration {iteration + 1}:")
        print(f"  Monomer: {monomer_name}")
        print(f"  Predicted Tg: {tg_estimate:.1f} K ({tg_estimate - 273.15:.1f}°C)")
        print(f"  Error: {error:.1f} K")

    # Select best candidate
    candidates_sorted = sorted(candidates, key=lambda x: x['Error'])
    best_candidate = candidates_sorted[0]

    print(f"\n{'='*60}")
    print(f"Best candidate: {best_candidate['Monomer']}")
    print(f"SMILES: {best_candidate['SMILES']}")
    print(f"Predicted Tg: {best_candidate['Predicted_Tg']:.1f} K ({best_candidate['Predicted_Tg'] - 273.15:.1f}°C)")
    print(f"Error from target: {best_candidate['Error']:.1f} K")
    print(f"{'='*60}")

    # Visualize the top candidates
    import matplotlib.pyplot as plt

    top_candidates = candidates_sorted[:5]
    names = [c['Monomer'] for c in top_candidates]
    tgs = [c['Predicted_Tg'] for c in top_candidates]
    errors = [c['Error'] for c in top_candidates]

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Subplot 1: Tg comparison
    ax1 = axes[0]
    colors = ['#27ae60' if i == 0 else '#3498db' for i in range(len(names))]
    bars = ax1.bar(names, tgs, color=colors, edgecolor='black', linewidth=2)
    ax1.axhline(target_tg, color='red', linestyle='--', linewidth=2,
                label=f'Target Tg = {target_tg}K')
    ax1.set_ylabel('Predicted Tg (K)', fontsize=12)
    ax1.set_title('Top Candidate Comparison', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(alpha=0.3, axis='y')
    plt.setp(ax1.xaxis.get_majorticklabels(), rotation=45, ha='right')

    for bar, val in zip(bars, tgs):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.0f}K', ha='center', va='bottom', fontsize=9, fontweight='bold')

    # Subplot 2: Error
    ax2 = axes[1]
    bars2 = ax2.bar(names, errors, color=colors, edgecolor='black', linewidth=2)
    ax2.set_ylabel('Error from Target (K)', fontsize=12)
    ax2.set_title('Prediction Error', fontsize=14, fontweight='bold')
    ax2.grid(alpha=0.3, axis='y')
    plt.setp(ax2.xaxis.get_majorticklabels(), rotation=45, ha='right')

    for bar, val in zip(bars2, errors):
        height = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.0f}K', ha='center', va='bottom', fontsize=9, fontweight='bold')

    plt.tight_layout()
    plt.savefig('polymer_design_optimization.png', dpi=300, bbox_inches='tight')
    plt.show()

    return best_candidate

# Example: material design targeting Tg = 150°C (423 K)
design_new_polymer_material(target_tg=423, max_iterations=10)

Exercises

# Polypropylene monomer (propylene)
smiles = "CC(C)"  # or "C=CC" (double bond notation)
print(f"Polypropylene SMILES: {smiles}")

from rdkit import Chem
from rdkit.Chem import Descriptors

smiles = "CC(C)(C(=O)OC)"
mol = Chem.MolFromSmiles(smiles)
mw = Descriptors.MolWt(mol)
print(f"PMMA molecular weight: {mw:.2f} g/mol")
# Output: 100.12 g/mol

slope = 0.03  # Ų/ps
D = slope / 6  # 3D: MSD = 6Dt
print(f"Diffusion coefficient: {D:.4f} Ų/ps")
# Output: 0.0050 Ų/ps

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs

smiles1 = "CC"
smiles2 = "CCC"

mol1 = Chem.MolFromSmiles(smiles1)
mol2 = Chem.MolFromSmiles(smiles2)

fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, 2, nBits=1024)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, 2, nBits=1024)

similarity = DataStructs.TanimotoSimilarity(fp1, fp2)
print(f"Tanimoto similarity: {similarity:.3f}")
# Output: around 0.5-0.7 (similar structures)

import numpy as np
from sklearn.metrics import mean_squared_error

y_true = np.array([345, 390, 395])
y_pred = np.array([350, 380, 400])

rmse = np.sqrt(mean_squared_error(y_true, y_pred))
print(f"RMSE: {rmse:.2f} K")
# Output: 7.45 K

import numpy as np

rg_data = np.array([14.5, 15.2, 14.8, 15.5, 14.9])
rg_mean = rg_data.mean()
rg_std = rg_data.std()

print(f"Rg mean: {rg_mean:.2f} Å")
print(f"Rg standard deviation: {rg_std:.2f} Å")
# Output: mean 14.98 Å, standard deviation 0.35 Å

from rdkit import Chem

def calculate_aromatic_ratio(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return 0

    aromatic_atoms = sum([1 for atom in mol.GetAtoms() if atom.GetIsAromatic()])
    total_atoms = mol.GetNumAtoms()

    aromatic_ratio = aromatic_atoms / total_atoms if total_atoms > 0 else 0
    return aromatic_ratio

# Test: polystyrene
smiles_ps = "CC(c1ccccc1)"
ratio = calculate_aromatic_ratio(smiles_ps)
print(f"Aromatic ratio: {ratio:.3f}")
# Output: 0.667 (6 aromatic atoms / 9 total atoms)

class PolymerAnalysisExtended(PolymerAnalysis):
    @staticmethod
    def compare_polymers(polymer_list):
        """
        Compare multiple polymers

        Parameters:
        - polymer_list: List of [(name, smiles), ...]

        Returns:
        - comparison_df: Comparison result DataFrame
        """
        import pandas as pd

        results = []
        for name, smiles in polymer_list:
            poly = PolymerAnalysis(name, smiles)
            poly.calculate_descriptors()
            poly.predict_tg()

            results.append({
                'Name': name,
                'MW': poly.descriptors['MW'],
                'LogP': poly.descriptors['LogP'],
                'Predicted_Tg': poly.predicted_tg
            })

        df = pd.DataFrame(results)
        print(df)
        return df

# Usage example
polymers = [
    ("PS", "CC(c1ccccc1)"),
    ("PMMA", "CC(C)(C(=O)OC)"),
    ("PVC", "CC(Cl)")
]

PolymerAnalysisExtended.compare_polymers(polymers)

import numpy as np
from rdkit import Chem
from rdkit.Chem import Descriptors

def genetic_algorithm_tg_optimization(target_tg=400, generations=5, population_size=10):
    """
    Tg optimization using a genetic algorithm
    """
    # Initial population (SMILES variations)
    monomers = ["CC", "CC(C)", "CC(c1ccccc1)", "CC(C)(C(=O)OC)", "CC(Cl)", "CC(C#N)"]

    population = np.random.choice(monomers, size=population_size)

    for gen in range(generations):
        # Fitness evaluation
        fitness_scores = []
        for smiles in population:
            mol = Chem.MolFromSmiles(smiles)
            if mol is None:
                fitness_scores.append(1e6)  # Penalty
                continue

            mw = Descriptors.MolWt(mol)
            hbd = Descriptors.NumHDonors(mol)
            hba = Descriptors.NumHAcceptors(mol)
            rot = Descriptors.NumRotatableBonds(mol)

            tg_est = 250 + 0.5*mw + 10*hbd + 5*hba - 20*rot
            fitness = abs(tg_est - target_tg)  # Smaller error is better
            fitness_scores.append(fitness)

        # Selection (elitism)
        elite_idx = np.argmin(fitness_scores)
        elite = population[elite_idx]

        print(f"Generation {gen+1}: Best = {elite}, Fitness = {fitness_scores[elite_idx]:.1f}K")

        # Generate next generation (simplified: random mutation)
        population = np.random.choice(monomers, size=population_size)
        population[0] = elite  # Preserve elite

    return elite

best_smiles = genetic_algorithm_tg_optimization(target_tg=400, generations=5, population_size=10)
print(f"\nBest SMILES: {best_smiles}")

from reportlab.lib.pagesizes import A4
from reportlab.pdfgen import canvas
from reportlab.lib.utils import ImageReader

def generate_pdf_report(polymer_analysis, filename='polymer_report.pdf'):
    """
    Generate PDF report

    Parameters:
    - polymer_analysis: PolymerAnalysis instance
    - filename: Output file name
    """
    c = canvas.Canvas(filename, pagesize=A4)
    width, height = A4

    # Title
    c.setFont("Helvetica-Bold", 20)
    c.drawString(50, height - 50, f"Polymer Analysis Report: {polymer_analysis.polymer_name}")

    # Basic information
    c.setFont("Helvetica", 12)
    y = height - 100
    c.drawString(50, y, f"SMILES: {polymer_analysis.smiles}")

    # Descriptors
    y -= 40
    c.drawString(50, y, "Molecular Descriptors:")
    y -= 20
    for key, value in polymer_analysis.descriptors.items():
        c.drawString(70, y, f"{key}: {value:.2f}")
        y -= 15

    # Tg prediction
    y -= 20
    c.drawString(50, y, f"Predicted Tg: {polymer_analysis.predicted_tg:.1f} K")

    # Insert graph (previously generated image)
    # img_path = f'{polymer_analysis.polymer_name}_structure.png'
    # c.drawImage(img_path, 50, y-300, width=400, height=250)

    c.save()
    print(f"PDF report generated: {filename}")

# Usage example
polymer = PolymerAnalysis("PMMA", "CC(C)(C(=O)OC)")
polymer.calculate_descriptors()
polymer.predict_tg()
generate_pdf_report(polymer, filename='PMMA_report.pdf')