# YASARA MACRO # TOPIC: 3. Molecular Dynamics # TITLE: Running an accurate molecular dynamics simulation in water # REQUIRES: Dynamics # AUTHOR: Elmar Krieger # LICENSE: GPL # DESCRIPTION: This macro sets up and runs a simulation. It can also continue a simulation that got interrupted. # 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, # or by uncommenting the line below and specifying it directly. #MacroTarget = 'c:\MyProject\1crn' # Extension of the cell on each side of the protein # '10' means that the cell will be 20 A larger than the protein extension=10 # pH at which the simulation should be run, by default physiological pH ph=7.4 # NaCl concentration in mass percent (0.9% is a physiological solution) nacl=0.9 # 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' # Pressure control mode # Default: Rescale the cell such that residues named HOH reach a density of 0.997 g/l. # For solvents other than water, you have to create your own solvent box # as described in the FillCellObj documentation and save it as .._solvent.sce. density=0.997 pressurectrl='SolventProbe,Name=HOH,Density=(density)' # Alternative: Uncomment below to calculate the pressure from the virial and # rescale the cell to reach a pressure of 1 bar. Use this method if you do not # know the correct density. #pressurectrl='Manometer,Pressure=1' # Constrain bond lengths to hydrogens and water bond angles to allow a larger timestep # If set to 'yes', the MD will run faster, but will be a bit less accurate constrain='no' # The format used to save the trajectories, sim or xtc format='sim' # Duration of the simulation, alternatively use e.g. 'duration=5000' to simulate for 5000 picoseconds duration='forever' # Flag to use a cubic simulation cell. This makes sure that also elongated # molecules can rotate freely during very long simulations. If set to 0, # the simulation cell will fit the solute more tightly, speeding up the simulation. cubic=1 # Forcefield to use (these are all YASARA commands, so no '=' used) ForceField AMBER03 # Cutoff Cutoff 7.86 # Cell boundary Boundary periodic # Use longrange coulomb forces (particle-mesh Ewald) Longrange Coulomb # Normally no change required below this point # ============================================ RequireVersion 9.9.25 # Keep the solute from diffusing around and crossing periodic boundaries CorrectDrift On # Treat all simulation warnings as errors that stop the macro WarnIsError On # 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 solvent present yet # 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 'pdb','yob','Error' size = FileSize (MacroTarget).(type) if size obj = Load(type) (MacroTarget) # Align object with major axes to minimize cell size NiceOriObj (obj) break if type=='Error' RaiseError "Initial structure not found. Make sure to create a project directory and place the structure there" # Prepare the structure for simulation CleanAll if Structure # Optimize the hydrogen-bonding network (more stable trajectories) OptHydAll # Create the simulation cell Cell Auto,Extension=(extension) if cubic # And make it cubic, taking the length of the X-axis, which is always the longest l = Cell Cell (l),(l),(l) SaveSce (MacroTarget) # Fill with water (always needs periodic boundaries), predict pKas, place counter ions Boundary periodic Experiment Neutralization WaterDensity (density) pH (ph) NaCl (nacl) pKaFile (MacroTarget).pka Speed Fast Experiment On Wait ExpEnd # Save scene with water SaveSce (MacroTarget)_water # Choose timestep and activate constraints if constrain=='yes' # Constrain bond lengths to hydrogens FixBond all,Element H # Constrain bond angles in water FixAngle Water,Water,Water # Multiple timestep: 1.3333 femtoseconds for intramolecular and 3*1.3333 = 4 fs for intermolecular forces TimeStep 3,1.3333 ts=4 # Save simulation snapshots every 6250 simulation steps # (with a timestep of 4 femtoseconds, that's 6540*4 fs = 25 picoseconds). savesteps=6250 else # Remove any constraints FreeBond all,all FreeAngle all,all,all # Smaller timestep, since we don't use constraints: 2*1.25 = 2.5 fs TimeStep 2,1.25 ts=2.5 # Save simulation snapshots every 10000 simulation steps # (with a timestep of 2.5 femtoseconds, that's 10000*2.5 fs = 25 picoseconds). savesteps=10000 # Temperature Temp (temperature) # Make sure all atoms are free to move FreeAll # Alread a snapshot/trajectory present? i=00000 filename='(MacroTarget)(i).(format)' running = FileSize (filename) if not running # Perform energy minimization Experiment Minimization Experiment On Wait ExpEnd # And now start the real simulation Sim On else # Simulation has been running before ShowMessage "Simulation has been running before, loading last snapshot..." Wait 1 # Switch console off to load the snapshots quickly Console Off if format=='sim' # Find and load the last 'sim' snapshot do i = i+1 found = FileSize (MacroTarget)(i).sim while found i=i-1 LoadSim (MacroTarget)(i) else # XTC format requires that the entire trajectory is read in to find the last one do i = i+1 eof,time = LoadXTC (filename),(i) ShowMessage 'Searching XTC trajectory for last snapshot, showing snapshot (i) at (0+time) fs' Sim Pause Wait 1 while !eof Sim Continue # Adjust savesteps to save snapshots in the same interval as previously if i>0 t = Time savesteps=0+t/(ts*i) #print 'Time=(t), TimeStep=(ts), savesteps=(savesteps)' HideMessage # Set temperature and pressure control TempCtrl Rescale PressureCtrl (pressurectrl) # Uncomment to add distance constraints # AddSpring O Res Lys 80,H Res Glu 84,Len=1.9 # And finally, make sure that future snapshots are saved Save(format) (filename),(savesteps) if duration=='forever' Console On Wait forever else Console Off # Wait for given number of picoseconds do Wait 10 t = Time while t