# YASARA MACRO # TOPIC: 3. Molecular Dynamics # TITLE: Analyzing a molecular dynamics trajectory # REQUIRES: Dynamics # AUTHOR: Elmar Krieger # LICENSE: GPL # DESCRIPTION: This macro analyzes a simulation and creates a table with energies and RMSDs from the starting structure, as well as the minimum energy structure and the time average structure with B-factors calculated from the root mean square fluctuations. If you want to do your own analysis, search for 'examples'. # The structure to analyze must be present with a .sce extension. # You can either set the target structure by clicking on Options > Macro > Set target, # or by uncommenting the line below and specifying it directly. #MacroTarget = 'c:\MyProject\1crn' # Number of the object whose RMSDs from the starting conformation will be calculated # If the protein is an oligomer, check the documentation of the 'Sup' command at 'analyzing a simulation' to avoid pitfalls. current = 1 # Forcefield to use (these are all YASARA commands, so no '=' used) # Use YASARA2 in YASARA Structure to include a quality Z-score ForceField AMBER03 # Cutoff Cutoff 7.86 # Cell boundary Boundary periodic # Use longrange coulomb forces (particle-mesh Ewald) Longrange Coulomb # The B-factors calculated from the root-mean-square fluctuations can be too large to fit them # into the PDB file's B-factor column. Replace e.g. 1.0 with 0.1 to scale them down to 10% bfactorscale=1.0 # Flag to save a PDB file of the solute snapshots for further analysis pdbsaved=0 # Selection of atoms to include for 'Calpha' RMSD calculation (also consider DNA/RNA) casel='CA Protein or C1* NucAcid' # No change required below this point # Do we have a target? if MacroTarget=='' RaiseError "This macro requires a target. Either edit the macro file or click Options > Macro > Set target to choose a target structure" Clear Console Off # Do we have a scene with water? scene = FileSize (MacroTarget)_water.sce if not scene RaiseError 'Could not find initial scene file (MacroTarget)_water.sce' # Load the scene LoadSce (MacroTarget)_water calphas = CountAtom (casel) if calphas>0 and calphas<3 # We cannot superpose 1 or 2 Calpha atoms casel='None' ShowMessage "Preparing analysis, please wait..." Wait 1 # See if structure validation checks should be done (require YASARA2 force field) checked=0 fof = ForceField if fof=='YASARA2' checked=1 # Duplicate the intial object for RMSD calculation initial = DuplicateObj (current) RemoveObj (initial) i=00000 emin=1e99 while 1 # See if next snapshot is present sim = FileSize (MacroTarget)(i).sim if not sim break # Yes, load it LoadSim (MacroTarget)(i) Sim Pause # Add time in picoseconds to table simtime = Time ShowMessage 'Analyzing snapshot (0+i) at (0+(simtime/1000)) ps' Wait 1 Tabulate (simtime/1000) # Add energy to table, first the total energy, then all components individually elist(i) = Energy ebndlist(i),eanglist(i),edihlist(i),eplnlist(i),ecoulist(i),evdwlist(i) = Energy All # The following examples provide a few hints for other things to analyze. # If you uncomment one of the examples, keep in mind that every value you # tabulate becomes an additional column in the table, so you must increase the # 'Columns' parameter of the SaveTab command further below (and the table # header will also no longer indicate the proper results unless you adapt it). # # Example: Measure the distance between the carboxyl group of Glu 123 (Cdelta) # and the guanidinium group of Arg 345 (Czeta): # Tabulate Distance CD Res Glu 123, CZ Res Arg 345 # # Example: Check if the carboxyl group of Glu 123 and the guanidinium group of # Arg 345 are currently hydrogen-bonded (1) or not (0): # hbolist() = ListHBoAtom Sidechain Res Glu 123, Sidechain Res Arg 345 # Tabulate (count hbolist>0) # # Example: Count the number of atoms in molecule A that are hydrogen bonded to B # hbolist() = ListHBoAtom Mol A,Mol B # Tabulate (count hbolist) # # Example: Calculate the current potential energy of residue Glu 123: # Tabulate EnergyRes Glu 123 # # Example: Calculate the backbone RMSD for residues 1-120: # Tabulate SupAtom Backbone Res 1-120 Obj (current),Backbone Res 1-120 Obj (initial) # # Save the minimum energy structure (ignoring solvent-solvent interactions) e = EnergyObj (current) if checked # Calculate structure validation Z-scores, using the formula from homology modeling for type in 'dihedrals','packing1d','packing3d' zscore(type) = CheckObj (current),(type) zscorelist(i)=zscoredihedrals*0.145+zscorepacking1d*0.390+zscorepacking3d*0.465 Sim Off if e