svats73 / md-sfa-msm

7 stars 2 forks source link

Weights are not writing in the output PDB #4

Closed sbhakat closed 2 months ago

sbhakat commented 2 months ago

@svats73 md-sfa plumed-bfactor executes properly but the output PDB is not saving the weights. Adding the standalone plumed-bfactor.py which does it correctly for comparison

#!/usr/bin/env python3

import argparse
from Bio import PDB
import re
import math

def parse_plumed_input(plumed_file):
    with open(plumed_file, 'r') as f:
        plumed_content = f.read()

    # Extract variable names and coefficients
    var_pattern = r'ARG=(.*?)\s+COEFFICIENTS=(.*?)\s+PERIODIC'
    match = re.search(var_pattern, plumed_content, re.DOTALL)
    if not match:
        raise ValueError("Could not parse PLUMED input")

    variables = match.group(1).split(',')
    coefficients = [float(c) for c in match.group(2).split(',')]

    return dict(zip(variables, coefficients))

def combine_weights(weights):
    combined_weights = {}
    angle_types = set()
    specific_angles = set()

    # Determine which angle types are present
    for var in weights.keys():
        parts = var.split('_')
        if len(parts) >= 3 and parts[0] == 'meanfree':
            angle_type = parts[1]
            specific_angle = parts[2]
            residue = int(parts[-1])
            angle_types.add(angle_type)
            specific_angles.add(specific_angle)

            if residue not in combined_weights:
                combined_weights[residue] = {}

            if specific_angle not in combined_weights[residue]:
                combined_weights[residue][specific_angle] = 0

            combined_weights[residue][specific_angle] += weights[var] ** 2

    # Calculate the combined weight for each residue
    for residue in combined_weights:
        total_weight = sum(combined_weights[residue].values())
        combined_weights[residue] = math.sqrt(total_weight)

    return combined_weights, angle_types, specific_angles

def apply_weights_to_pdb(pdb_file, weights, output_file):
    parser = PDB.PDBParser()
    structure = parser.get_structure("protein", pdb_file)

    for model in structure:
        for chain in model:
            for residue in chain:
                res_id = residue.id[1]
                if res_id in weights:
                    for atom in residue:
                        atom.set_bfactor(weights[res_id])

    io = PDB.PDBIO()
    io.set_structure(structure)
    io.save(output_file)

def main():
    parser = argparse.ArgumentParser(description="Apply PLUMED dihedral angle weights to PDB B-factors")
    parser.add_argument("plumed_file", help="Input PLUMED file")
    parser.add_argument("pdb_file", help="Input PDB file")
    parser.add_argument("output_file", help="Output PDB file")
    args = parser.parse_args()

    try:
        weights = parse_plumed_input(args.plumed_file)
        combined_weights, angle_types, specific_angles = combine_weights(weights)
        apply_weights_to_pdb(args.pdb_file, combined_weights, args.output_file)
        print(f"PDB file updated with new B-factors based on PLUMED dihedral angle weights.")
        print(f"Angle components processed: {', '.join(sorted(angle_types))}")
        print(f"Specific angles found: {', '.join(sorted(specific_angles))}")
        print(f"Output saved to {args.output_file}")
    except Exception as e:
        print(f"An error occurred: {str(e)}")

if __name__ == "__main__":
    main()
sbhakat commented 2 months ago

@svats73 this issue is still present.

sbhakat commented 2 months ago

Functionally added.

sbhakat commented 2 months ago

Functionally added.