# YASARA MACRO # TOPIC: 3. Molecular Dynamics # TITLE: Running a steered molecular dynamics simulation in water # REQUIRES: Dynamics # AUTHOR: Elmar Krieger # LICENSE: GPL # DESCRIPTION: This macro sets up and runs a steered simulation to pull a ligand from its binding site. As soon as they have been pulled apart completely and atoms start crossing periodic boundaries, the macro works no longer correctly. You can set the pulling 'acceleration' below, which is not changed during the simulation. So this macro can currently not tell you how strong you need to pull, it can only tell you whether the chosen acceleration is strong enough or not to pull the ligand off. # Parameter section - adjust as needed, but NOTE that some changes only take # effect if you start an entirely new simulation, not if you continue an existing one. # ==================================================================================== # The structure to simulate must be present with a .pdb or .sce extension. # If a .sce (=YASARA scene) file is present, the cell must have been added. # You can either set the target structure by clicking on Options > Macro > Set target, # by providing it as command line argument (see docs at Essentials > The command line), # or by uncommenting the line below and specifying it directly. #MacroTarget '/home/myself/test' # Selection of ligand residue(s), adapt as needed. For PDB file 1BFB, you would need # ligandsel='Res UAB SGN IDS'. If the ligand is a protein, you can simply use its molecule name, # e.g. 'Mol B'. If you leave ligandsel empty, YASARA will choose selected atoms as the ligand or, # if none are selected, try to identify the ligand automatically. ligandsel='' # Selection of receptor residues. By default, the ligand will be pulled away from the center # of mass of the selected residues. If the receptor has an unusual shape, select residues # immediately 'below' the ligand, e.g. the active site 'Res Ser 360 Glu 361'. If the ligand is # also a protein, add the receptor's molecule name, e.g. 'Mol A'. # If you provide an arrow to indicate the pulling direction, then select all receptor atoms, # YASARA will pick the right subset to update the pulling direction automatically. To verify that # your selection is right, just look at the red pulling arrow shown by YASARA during the steered MD. receptorsel='Protein' # Simulation temperature # If you run at a temperature that differs from 298K, you also need # to adapt the pressure control below, look in the PressureCtrl documentation. temperature='298K' # Solvent density in g/ml (0.997 is water at 298K) # If you run at a significantly different temperature, this density needs to # be adapted. Look at the PressureCtrl documentation. # For solvents other than water, you have to create your own solvent box # as described at the FillCellObj documentation and save it as .._solvent.sce density=0.997 # Name of solvent molecules used to measure solvent density solvent='HOH' # The starting and minimum pulling acceleration in picometers/picosecond^2 acceleration=500. # The steps in picometers/picosecond^2 in which acceleration will be increased, set 0 to get a constant acceleration accdelta=100 # Time in picoseconds when the pulling starts # (everything before can be considered equilibration) pullstart=3 # pH at which the simulation should be run, by default physiological pH ph=7.4 # The ion concentration as a mass fraction, here we use 0.9% NaCl (physiological solution) ions='Na,Cl,0.9' # Multiple timestep: 1.25 femtoseconds for intramolecular and 2*1.25 fs for intermolecular forces tslist=2,1.25 # Save simulation snapshots every 1000 femtoseconds = 1 picoseconds saveinterval=1000 # The format used to save the trajectories: YASARA 'sim', GROMACS 'xtc' or AMBER 'mdcrd'. format='sim' # The simulation ends when the distance between the centers of mass of receptor and ligand exceeds this value [A] enddis=15 # Use longrange coulomb forces (particle-mesh Ewald) longrange='Coulomb' # Forcefield to use (these are all YASARA commands, so no '=' used) ForceField AMBER14 # Cutoff Cutoff 8 # Number of simulation steps per screen update and pairlist update. # Since the macro can apply pulling forces only once per screen update, screen and pairlist are updated each # simulation step to ensure that forces are applied continuously. This slows down the simulation considerably # to about 40% of the normal speed. If speed is an issue, you can change the command below to 'SimSteps 2,2' # (50% of normal speed), 'SimSteps 3,3' (60% of normal speed) or 'SimSteps 5,5' (80% of normal speed). Going # higher is not recommended, since acceleration*N will be applied every Nth step, and no acceleration in between. SimSteps Screen=1,Pairlist=1 # Normally no change required below this point # ============================================ # Requires 'with distance from local point' RequireVersion 19.1.27 # Keep the solute from diffusing around and crossing periodic boundaries CorrectDrift On # Treat all simulation warnings as errors that stop the macro WarnIsError On # Normally no change required below this point # Temperature Temp (temperature) # Calculate total timestep, we want a float, so tslist2 is on the left side ts=tslist2*tslist1 # Snapshots are saved every 'savesteps' savesteps=saveinterval/ts Console Off # 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 # Do we already have a scene with water or other solvent? waterscene = FileSize (MacroTarget)_water.sce solventscene = FileSize (MacroTarget)_solvent.sce if waterscene LoadSce (MacroTarget)_water elif solventscene LoadSce (MacroTarget)_solvent else # No scene with water present yet if solvent!='HOH' RaiseError "When using a solvent other than water, the cell cannot be filled automatically. Look at the documentation of the FillCellObj command" # Do we have a scene at all? scene = FileSize (MacroTarget).sce if scene LoadSce (MacroTarget) else # No scene present, assume it's a PDB or YOB file for type in 'yob','pdb' size = FileSize (MacroTarget).(type) if size break if !size RaiseError 'Initial structure not found, expected (MacroTarget).pdb or .yob. Make sure to create a project directory and place the structure there' # Load structure Load(type) (MacroTarget) # Align model with major axes to minimize cell size NiceOriAll CleanAll # Create the simulation cell, make sure that there is enough space to reach 'enddis' Cell Auto,Extension=(enddis*2) SaveSce (MacroTarget) # Fill with water, predict pKas, place counter ions Boundary periodic Experiment Neutralization WaterDensity (density) pH (ph) Ions (ions) pKaFile (MacroTarget).pka Speed Fast Experiment On Wait ExpEnd # Save scene with water SaveSce (MacroTarget)_water # Identify the ligand: choose the largest hetgroup with >6 atoms if any # Changes here also in md_analyze if ligandsel=='' # Check if atoms are selected selatomlist()=ListAtom selected if count selatomlist ligandsel='(join selatomlist)' if ligandsel=='' mols=CountMol Obj !Water if mols>1 reslist()=ListRes Hetgroup !Water with >0 bonds to all reslenmax=6 ligandname='' for res in reslist reslen=CountAtom (res) if reslen>reslenmax reslenmax=reslen ligandsel=res for type in 'ligand','receptor' # Make sure that ligand and receptor are really present if (type)sel=='' (type)atoms=0 else (type)atoms = CountAtom ((type)sel) if !(type)atoms RaiseError 'This macro requires that the (type) is selected explicitly. The current (type) ((type)sel) was not found' # Get central ligand/receptor atom for arrow visualization cen()=GroupCenter ((type)sel) for atomsel in 'CA','!H' (type)atom=ListAtom (atomsel) and ((type)sel) with minimum distance from local point (join cen) if (type)atom break # Get ligand mass ligmasslist()=MassAtom (ligandsel) ligmass=sum ligmasslist # If an arrow object is present it will be used to determine the steering direction arrowobj=ListObj 'Arrow' if arrowobj cellsize1,cellsize2,cellsize3,angle1,angle2,angle3=Cell if angle1!=90 or angle2!=90 or angle3!=90 RaiseError "Steered MD with a user-provided arrow requires a cuboid cell (with all three cell angles at 90 degrees)" diag=norm cellsize # Create a chain of dummy atoms chainobj=BuildGrid 0,(2*diag) # Copy global position/orientation of arrow to chain of dummy atoms TransferObj (chainobj),(arrowobj),Local=Keep # Determine number of atoms that will be used to place the arrow neighboratoms = CountAtom (receptorsel) if neighboratoms>20 neighboratoms=20 # Select all receptor atoms near the chain cutoff=0.5 do arrowatomlist()=ListAtom (receptorsel) with distance<(cutoff) from Obj (chainobj) cutoff=cutoff+0.5 while count arrowatomlistpullstart # Steer the simulation ligpos()=GroupCenter (ligandsel) if arrowobj # Calculate pulling 'dir'ection from arrowatomlist, and 'dis'tance of the ligand from its starting position recpos()=GroupCenter (join arrowatomlist) else # Calculate pulling direction from centers of mass, and 'dis'tance between receptor and ligand recpos()=GroupCenter (receptorsel) disvec()=ligpos-recpos dis=norm disvec if disstart==-1 disstart=dis if !arrowobj dir()=disvec/dis dis=dis-disstart HideArrowAtom (ligandatom) ShowArrow AtAtom,(ligandatom),StepFromAtom,(ligandatom),(dir*10) if dis>enddis ShowMessage 'The distance has reached (enddis) A in (0.00+t) ps, simulation has been stopped' break # force[pN,kg*pm/s^2] = acceleration[pm/ps^2]*1e24[ps^2/s^2]*ligmass[g/mol]*1e-3[kg/g]*(1/AvoConst)[mol] force=1e21/AvoConst*acceleration*ligmass if dis>dismax progtime=t if (dismax>0) # work[EnergyUnit] = dis[A]*1e-10[m/A]*force[pN]*1e-12[N/pN]*JToUnit[EnergyUnit/J] work=JToUnit*1e-22*((dis-dismax)*force)+work Tabulate (dis),(t),(acceleration),(force) dismax=dis # Get current steps per screen update to adjust the acceleration steps=SimSteps if !(step%400) # Check if there is any progress if dismax<=progdismax # No, pull more strongly. # We use additive addaption instead multiplicative because the acceleration step # would get too large at the point where the ligand breaks free. acceleration=acceleration+accdelta progdismax=dismax if !(step%20) # Tabulate (dismax),(t),(acceleration),(dis) # Check if we just snapped out and are too fast. Calculate speed in m/s speed=(dismax-speeddismax)*100000/(steps*tslist1*tslist2) if speed>speedallowed # Too fast, reduce acceleration factor=1.-speedallowed/speed acceleration=acceleration*(1.-factor*factor) speeddismax=dismax if acceleration