The simulate Stage

Simulation/Desmond
  • Overview
  • Examples
  • Default Configuration

The simulate stage is the most fundamental Desmond stage – its basic functionality is to run a Desmond simulation. The simulation is defined by the keywords of the stage, which can modify thermodynamic properties, apply restraints to atoms, activate plugins, and more.

Sample Simulate Stage

To understand the basics of a simulate stage, see the sample breakdown of a Desmond minimization stage below:

simulate 
      {
	title       = "Brownian Dynamics NVT, T = 10 K, small timesteps, and restraints on solute heavy atoms, 100ps"	
	time        = 100
	timestep    = [0.001 0.001 0.003]
	temperature = 10.0
	ensemble = 
		{
		class = "NVT"
		method = "Brownie"

		# The delta_max stops any particle ever moving >0.1 A regardless of
		# the force compelling it to do so.
		brownie = {delta_max = 0.1}
		}

	restraints.new = 
		[{
		name = posre_harm
		atoms = solute_heavy_atom
		force_constants = 50.0
		}]
	}

The keywords of this simulation stage define a “first-pass” minimization of the input system - a rough relaxation of a potentially unstable input structure. By modifying the given keywords, the stage could behave quite differently, highlighting how versatile the simulate stage is.

The first keyword is the title, which is a general-purpose keyword that can serve as a description for a customized stage. This is heavily used in template MSJs generated through Maestro to offer some insight into what the default generated stages are doing.

This simulation duration is 100ps, as set by the time keyword, and at a temperature of 10K. The integrator is configured with short timesteps (as mentioned in the title), and restraints are applied to the solute heavy atoms. Effectively, this stage is relaxing the solvent.

The ensemble keyword is used to specify the integrator that defines the time evolution of the system (method), as well as the ensemble class (class). In the example above, the Brownie integrator is used, which is able to handle very large potential energies/forces without having the system explode. This may be useful at the start of a simulation, where there may be extremely large forces due to misplaced solvents or other clashes.

Simulate Keywords

The keywords specified here work for all Desmond simulation stages (i.e. lambda_hopping, replica_exchange), not just the simulate stage. For a full list of the keywords in the simulation stage as well as their defaults, see the Default MSJ tab.

Table 1. Keywords for the simulate stage

Keyword

Description

cfg_file

This specifies the .cfg configuration file to be used for the given simulation, as an input to the Desmond engine. This keyword is not required; multisim is capable of creating and supplying the relevant .cfg file.

restraints

Specify internal coordinate restraints on bonds, angles, and torsions for specified atoms. Restraints can be a critical part of a simulation; see Simulation Restraints below for more information on this keyword and the Examples tab for examples.
Restraints set in all stages except the assign_forcefield stage apply only to the current stage. Restraints set in the assign_forcefield stage are “permanent”—they are inherited by all subsequent stages. Permanent restraints can be overridden in a particular stage, but will still be applied in subsequent stages. The restraints.existing keyword can be used to inherit restraints from a previous stage.

atom_group

Define atom groups within the .cms file (defined by the i_ffio_grp_name property). Atom groups can be restrained, or associated with particular thermostats. This keyword can take on the following values:

  • none—Remove all atom groups.

  • retain—Keep all atom groups from the previous stage.

  • { atom = atoms index = i name = name }—atom group block. Put the specified atoms in the atom group i that has the name name. The atom keyword accepts an atom list, an ASL expression prefixed by asl:, or keywords heavy atom, solute, solute_heavy_atom, solvent, solvent_heavy_atom. The index i is the value for the i_ffio_grp_name property. Desmond only supports index numbers from 0 to 7.

  • [ group1 group2 ... ]—Specify multiple atom groups. Each group can be specified as a block, in the above format.

  • {generator = "energy" residue_groups = group ligand_group = group cutoff = n }—Set up individual energy groups within a certain cut off. residue_groups and ligand_group can be set to "auto", which automatically selects the appropriate atoms for these groups. For example, atom_group = {generator="energy" residue_groups="asl:protein" ligand_group="ligand" cutoff=7.0} automatically sets up individual energy groups for all protein residues that are within 7Å of the "ligand_group".

Atom groups defined in the assign_forcefield stage are “permanent”—they are inherited by all subsequent stages. This can be overridden in a particular stage, but will still be applied in subsequent stages unless the retain keyword is used to keep atom groups from the previous stage.

Simulation Restraints

Table 2. Keywords accepted for simulation restraints

Keyword

Description

Syntax

restraints.new

Add new restraints to any simulate stage in the MSJ file.

 
restraints.new = [
{
 name = name 
 atoms = [atoms]
 force_constants = [value]
 r0 = value
 theta0 = value
 phi0 = value
 sigma = sigmavalue
}
]

restraints.existing

Specifies how to deal with existing restraints from a previous stage. Default value is ignore, which overwrites old restraints. This keyword cannot be used to control permanent restraints set in the assign_forcefield stage. To override permanent restraints for a stage, you must use restraints.new to define new restraints for that stage.

 
restraints.existing = retain

print_restraint

Prints restraints to the Multisim log file.

 
print_restraint = true

There are some important considerations when applying restraints to a simulation.

The number of selections defined in atoms has to match the type of restraint (i.e., two for distance, three for angle, and four for dihedral/improper restraints). Any number of atoms may be used for positional restraints. Each selection has to specify a single atom. Individual atoms can be selected using the ASL atom.num or a.num syntax. For example, to select atom number 1, the ASL would be "atom.num 1" or "a.num 1".

The r0, theta0 and phi0 settings respectively specify the equilibrium distance, angle, and dihedral/improper torsion. If they are not specified, the atomic coordinates from the input structures are used to compute them.

Force constants should be specified in units of kcal mol-1 Å-2. You can include anisotropic restraints on a positional restraint by assigning the force constant as follows:

force_constant = [x-valuey-valuez-value]

See the Examples tab for examples of various restraints that can be applied.

Note: the print_restraint = true setting prints the restraints in the multisim log file. restraints.existing = retain adds the new terms defined in restraints.new to existing restraints.

simulate {
 print_restraint = true
 restraints.existing = retain
 restraints.new = [
 { 
   name = stretch_harm
   atoms = ["atom.num 23" "atom.num 5"]
   force_constants = [25.0]
   r0 = 4.4
 }
 ]
}

Note: the print_restraint = true setting prints the restraints in the multisim log file. restraints.existing = retain adds the new terms defined in restraints.new to existing restraints.

simulate {
 print_restraint = true
 restraints.existing = retain
 restraints.new = [					
 {
   name = angle_harm
   atoms = ["atom.num 23" "atom.num 5" "atom.num 6"]
   force_constants = [35.0]
   theta0 = 4.4
 }
 ]
}

Note: the print_restraint = true setting prints the restraints in the multisim log file. restraints.existing = retain adds the new terms defined in restraints.new to existing restraints.

simulate {
 print_restraint = true
 restraints.existing = retain
 restraints.new = [
 {
  name = improper_harm
  atoms = ["atom.num 1" "atom.num 2" "atom.num 3" "atom.num 5"]
  force_constants = [50.0]
  phi0 = 120.0
 }
 ]
}

Note: the three numbers specified in force_constants = [10 10 10] are the force constants in the x,y, and z direction.

simulate {
 time = 2.0
 timestep = [ 0.002 0.002 0.006 ]
 ensemble = NVT
 temperature = 300.0
 print_restraint = true
 restraints.new = [
 {
   name = posre_harm
   atoms = ["atom.num 1, 2"]
   force_constants = [10 10 10]
 }
 ]
}

Note: the sigma value defines the range of the flat-bottom well symmetrically centered around r0, theta0, phi0 where no restraining force acts on the system.

simulate {
 print_restraint = true
 restraints.existing = retain
 restraints.new = [
 { 
   name = stretch_fbhw
   atoms = ["atom.num 1713" "atom.num 599"]
   force_constants = [25.0]
   r0 = 2.5
   sigma = 0.25
 }
 {
   name = angle_fbhw
   atoms = ["atom.num 1568" "atom.num 609" "atom.num 1569"]
   force_constants = [25.0]
   theta0 = 120.0
   sigma = 2.0
 }
 {
   name = improper_fbhw
   atoms = ["atom.num 101" "atom.num 102" "atom.num 103" "atom.num 104"]
   force_constants = [25.0]
   phi0 = 65.0
   sigma = 5.0
 }
 ]
}

Note: COM-based restraints are set with the group_restraints keyword. The following msj stage defines two groups using VMD selection language and then adds a flat-bottom harmonic well restraint between them. You will need to add the “backend.force.term = {}” section to your normal msj file.

simulate {
   cfg_file = "desmond_md_job_1.cfg"
   backend.force.term = {
      group_restraints = {
         groups = [
              { name=protein1   center=mass  query="protein and resid 57" } #query definitions for protein1
              { name=protein2   center=mass  query="protein and resid 136" } #query definitions for protein2
          ]
          restraints = [ { form            = harmonic_flat_bottom
                           groups          = [ protein2 protein1 ]
                           spring_constant = 0.35 # spring constant in kcal/mol/A^2
                           equilibria      = [ 8 14 ] # equilibrium distances in angstroms 
                         } ]
         type = group_restraints
      }
      list = [group_restraints]
   }
   jobname  = "$MAINJOBNAME"
   dir      = "."
   compress = ""
}

Change following lines:
- query definitions for your protein (remember this uses VMD selection language not ASL)
- spring_constant (kcal/mol/A^2)
- equilibria (currently set to exert no force if the distance is between 8 and 14 Angstroms, but outside of this the distance is harmonically restrained)

Then, run the job using the normal multisim command.

Default configurations for the simulate stage. A {type} value is shown for keywords which do not have defaults set. See General multisim Stage Keywords for descriptions of general keywords

{
   annealing = false
   atom_group = none
   backend = {
   }
   bigger_rclone = false
   box = {}
   cfg_file = ""
   checkpt = {
      first = 0.0
      interval = 240.06
      name = "$JOBNAME.cpt"
      write_last_step = true
   }
   compress = "$MAINJOBNAME_$STAGENO-out.tgz"
   coulomb_method = useries
   cpu = 1
   cutoff_radius = 9.0
   dipole_moment = false
   dir = "$[$JOBPREFIX/$]$[$PREFIX/$]$MAINJOBNAME_$STAGENO$[_lambda$LAMBDA$]"
   dryrun = false
   ebias_force = false
   effect_if = {list}
   elapsed_time = 0.0
   energy_group = false
   eneseq = {
      first = 0.0
      interval = 1.2
      name = "$JOBNAME$[_replica$REPLICA$].ene"
   }
   ensemble = NPT
   fep = {
      i_window = {int}
      lambda = "default:12"
      output = {
         first = 0.0
         interval = 1.2
         name = "$JOBNAME$[_replica$REPLICA$].dE"
      }
      trajectory = {
         record_windows = [0 -1 ]
      }
      type = small_molecule
   }
   gaussian_force = false
   gcmc = {
      ene_name = "$JOBNAME$[_replica$REPLICA$]_gcmc.ene"
      first = 0.0
      gcmc_region = {
         cell_size = 0.22
         exclusion_radius = 2.2
         global_switching = {
            frequency = 0.2
            move_factor = 3.0
            spacing_factor = 2.0
         }
         region_buffer = 4.0
         track_voids = true
      }
      interval = 4.8
      ligand_file = {str}
      moves = {
         moves_per_cycle = 5000
      }
      mu_excess = -6.18
      seed = random
      solvent = {
         s_file = ""
      }
      solvent_density = 0.03262
      verbose = 0
   }
   glue = solute
   host = "$SUBHOST"
   jin_file = []
   jin_must_transfer_file = []
   jlaunch_opt = ["" ]
   jobname = "$MAINJOBNAME_$STAGENO$[_lambda$LAMBDA$]"
   jout = ""
   lambda_dynamics = false
   maeff_output = {
      center_atoms = solute
      first = 0.0
      interval = 120.0
      name = "$JOBNAME$[_replica$REPLICA$]-out.cms"
      periodicfix = true
      trjdir = "$JOBNAME$[_replica$REPLICA$]_trj"
   }
   meta = false
   meta_file = {str}
   msd = false
   prefix = ""
   pressure = 1.01325
   pressure_tensor = false
   print_expected_memory = false
   print_restraint = false
   randomize_velocity = {
      first = 0.0
      interval = inf
      seed = 2007
      temperature = "@*.temperature"
   }
   restrain = none
   restraints = {
      existing = ignore
      new = []
   }
   rnemd = false
   should_skip = false
   should_sync = true
   simbox = {
      first = 0.0
      interval = 1.2
      name = "$JOBNAME$[_replica$REPLICA$]_simbox.dat"
   }
   spatial_temperature = false
   struct_output = ""
   surface_tension = 0.0
   taper = false
   temperature = 300.0
   time = 1200.0
   timestep = [0.002 0.002 0.006 ]
   title = {str}
   trajectory = {
      center = []
      first = 0.0
      format = dtr
      frames_per_file = 250
      interval = 4.8
      name = "$JOBNAME$[_replica$REPLICA$]_trj"
      periodicfix = true
      write_last_vel = false
      write_velocity = false
   }
   transfer_asap = false
   wall_force = false
   window = {}
}