Molecular Dynamics and Monte Carlo Simulation Opcodes

The opcodes in this section are linked below:

MDYN — Molecular DYNamics

MDSA — Molecular Dynamics structure SAmpling

MDMT — Molecular Dynamics MomenTum reset

MDAR — Molecular Dynamics ARea monitoring

MDAP — Molecular Dynamics Atom Position monitoring

MDDI — Molecular Dynamics DIstance monitoring

MDBA — Molecular Dynamics Bond Angle monitoring

MDDA — Molecular Dynamics Dihedral Angle monitoring

MHBD — Monitor the occurrence of a Hydrogen BonD

MDAV — Molecular Dynamics coordinate AVeraging

MDIT — Molecular Dynamics Initial Temperature

MDFT — Molecular Dynamics Final Temperature

MDMC — Molecular Dynamics / Monte Carlo mixed mode

MCSD — Monte Carlo / Stochastic Dynamics mixed mode

MDVE — Molecular Dynamics VElocity file

MDFR — Molecular Dynamics update FRequency

MDRE — Molecular Dynamics REset

IMPS — perform IMPortance Sampling

IMPO — IMPortance sampling Options

ITBS — Improper Torsion energy BiaSing

IMCC — IMportance sampling Chirality Check

TCMP — Torsional CoMParison

ZMAT — Z-MATrix

MDYN — Molecular DYNamics

Carry out molecular dynamics. Structures must be fully minimized before molecular dynamics is initiated and initial conditions (temperature, time step, and total time) must be selected. Temperature maintenance with a constant temperature bath is recommended. If too high a temperature or too large a time step is selected, the structure will eventually gain too much kinetic energy and blow up (program automatically aborts if the kinetic energy/atom exceeds 100 kJ/mol). The structure may also blow up if it is not sufficiently minimized and energy is added (e.g., MDIT). If this occurs, reduce the time step, turn on SHAKE (arg 2) and try again.

During the run, average temperatures (both over the most recent 0.5 ps and over the entire simulation so far), the instantaneous total energy (potential plus kinetic energy) and the average enthalpy over the entire simulation (taken as the average steric energy) are written to the log file. If you choose the full printing option (arg1=2), the information is written to the .mmo file as well. If you are using constant energy molecular dynamics (no temperature control), the energy may change over the course of the simulation: if the temperature listed after 0.5 ps or so is more than 5% too high, then it is likely that the starting geometry was not completely minimized or that a lower energy nearby minimum has been found.

If you are not using temperature control and wish to readjust the temperature, stop the batch job and set the desired temperature with arg 7. Start the run again and the temperature will adjust itself over the first several ps.

It is best to do the simulation in stages (e.g., to begin by equilibrating the system prior the actual run) having a MDYN command for each stage. If the MDRE (reset) command is not used, the run will continue from the previous end point. Equilibration runs are suggested for good thermodynamics especially if there is no reason to believe that the starting structure is the fully minimized global minimum. We suggest at least 5-50 ps of equilibration before data collection for all molecular dynamics operations other than conformation searching where equilibration is unnecessary. Experience shows that ca .2-.5 times the number of heavy atoms gives a reasonable equilibration time in picoseconds.

The time step (in femtoseconds (Default = 1.5 fs), arg5) is the time interval of each step of the molecular dynamics run. Energy stability is best at or below 1.0 fs but time steps up to 2.5 fs may be used to reduce the overall time required for the simulation. If too large a time step is used, the structure will gain energy each step and eventually blow up (simulation will abort). In general, higher temperatures (kinetic energy) require smaller timesteps. If mediocre energy stability is acceptable (e.g., for local conformation searching), then an initial temperature of 300-400 K with a timestep of 1.7-2.0 fs provides a good compromise of speed and stability if SHAKE (constrains H bond lengths) is used. If precise energies are important, a timestep of 1.0 fs with temperature control (see arg 7), or a timestep of 1.5 fs with SHAKE (arg 2) and temperature control (arg 7) is recommended. Molecules with explicit hydrogens (e.g., when using MM2) may require a smaller timestep and/or use of SHAKE (see arg 2) to constrain bonds to hydrogen.

At the end of each MDYN run, the final structure must be written to the output file using the WRIT command, otherwise no structure will be saved in the output file (unless MDSA is used). It is often useful to save final velocities using the MDVE command so that subsequent runs can start with equilibrated structure (read the final structure and the corresponding velocity file for the next run). Multiple MDYN commands within the same MacroModel run need not be accompanied with WRIT or MDVE commands.

DEBG 1 gives energy output to the log file ten times as often as it otherwise would.

nn fs

MDSA — Molecular Dynamics structure SAmpling

Also used to specify structure sampling for Monte Carlo procedures.

Sampling of structures during molecular dynamics or Monte Carlo run by time. MDSA commands must come before the MDYN, MCSD or MCLO command from which the samples are to be taken.

Initially, MacroModel is set to not clear the output structure file at the start of each MDYN run. If arg7 of MDSA is set to a value greater than zero, the default behavior is turned off and the output structure file is wiped clean at the start of each subsequent MDYN run. Setting arg7 of MDSA to a value less than zero “switches off” the effect of having set arg7 to a value greater than zero in a previous MDSA line in the command file.

If arg7 is set to zero, then the currently set MDSA behavior will be used. For example, if arg7 had been set to a value greater than zero at an earlier point in the .com file, setting arg7 to zero in subsequent line would cause the previously set behavior (wipe the output structure file) to continue. On the other hand, if arg7 had been set to a value less than zero at an earlier point in the .com file, setting arg7 to zero in subsequent line would cause the previously set behavior (do not wipe the output structure file) to continue.

If the first MDSA line in the .com file is set to 0, then the previously set behavior (MacroModel’s default behavior) of not wiping the output structure file would be used.

MDMT — Molecular Dynamics MomenTum reset

Used to control translational and rotational (angular) momentum which will be periodically (0.1 ps by default) reset to zero. MDMT is operational by default (i.e., MDMT is unnecessary unless it has previously been turned off), and its normal mode is constraint of the total molecular system to have zero translational momentum.

Note: Structures viewed graphically will not appear to translate or rotate even if momentum is not zeroed because they are translated/rotated to the original position and overall orientation before writing to any output file.

The frequency of momentum zeroing is set by arg 5.

MDAR — Molecular Dynamics ARea monitoring

Atomic area monitoring during the molecular dynamics run. Output consists of histogram-like compilations of individual atomic areas and the averages. The output will be written to the mmo output file. If no probe radius is specified, then the areas will reflect the van der Waals surface areas for the atoms specified. If all atoms are to have the same probe radius, use MDAR commands loaded with four atoms at a time using args1-4.

Data is accumulated every 10 time steps throughout the simulation and also when the structure is sampled. DEBG switch 1 lists each value to the .log file. A maximum of 250 atoms may be monitored simultaneously.

This command must precede the MDYN command to be monitored.

MDAP — Molecular Dynamics Atom Position monitoring

Atom Position monitoring during the molecular dynamics run. Output consists of average rms deviation (Å) from original position. The output is written to the log and .mmo output file.

Data is accumulated every 10 time steps throughout the simulation and also when the structure is sampled. DEBG switch 1 lists each value to the .log file. A maximum of 250 atoms may be monitored simultaneously.

Note: This command must precede the MDYN command to be monitored.

MDDI — Molecular Dynamics DIstance monitoring

Distance monitoring during a molecular dynamics run. Output consists of histogram-like compilations of internuclear distances/frequencies and the average distance. Histograms have 100 bins covering the range of arg6 to arg6+arg5*100. The default range is 0–10 Å in 0.1 Å increments. The output is written to the .mmo output file.

Data is accumulated every 10 time steps throughout the simulation and also when the structure is sampled. DEBG 1 lists each value to the log file. A maximum of 100 distances may be monitored simultaneously.

This command must precede the MDYN command for monitoring to take place.

MDBA — Molecular Dynamics Bond Angle monitoring

Bond Angle monitoring during the molecular dynamics run. Output consists of histogram-like compilations of angles/frequencies and the average angle. Default resolution for the histograms is 5 degrees (see arg5). The output will be written to the .mmo output file.

Data is accumulated every 10 time steps throughout the simulation and also when the structure is sampled. DEBG 1 lists each value to the .log file. A maximum of 100 angles may be monitored simultaneously.

This command must precede the MDYN command for monitoring to take place.

MDDA — Molecular Dynamics Dihedral Angle monitoring

Dihedral Angle monitoring during the molecular dynamics run. Output consists of histogram-like compilations of angles/frequencies and the average angle. Default resolution for the histograms is 10 degrees (see arg5). The output will be written to the .mmo output file (averages and histogram, printing must be on —see MDYN arg 1) and to the job log file (averages only). The average torsion angle cosine and the average cosine squared are also provided for use in Karplus-like coupling constant calculations.

Data is accumulated every 10 time steps throughout the simulation and also when the structure is sampled. DEBG switch 1 lists each value to the .log file. A maximum of 100 angles may be monitored simultaneously.

This command must precede the MDYN command to be monitored.

MHBD — Monitor the occurrence of a Hydrogen BonD

This command allows monitoring of the number of times that distance and angle criteria are met during a molecular or stochastic dynamics simulation. The normal use of this command is to monitor the population of hydrogen bonds during a simulation. At each sampling point a distance and one or two angles are determined, and if these values meet criteria defined in the MHBD command line then a hydrogen bond is considered to exist. At the end of the simulation the occurrence of the specified hydrogen bond is reported to the log file as the percentage of structures in which the bond occurs.

Example:

MDAV — Molecular Dynamics coordinate AVeraging

Averaging of atomic Cartesian coordinates during the molecular dynamics run. The average atomic coordinates will be used as the structure coordinates at the end of the dynamics run. The coordinates will be written to the output structure file. To save these coordinates in the output file, it is necessary to use a WRIT command after the MDYN command.

This command must precede the MDYN run to be averaged.

MDIT — Molecular Dynamics Initial Temperature

Initial Temperature (K, Absolute) for the start of the run. This temperature is converted to the corresponding kinetic energy, and that energy is used to initialize the velocities of all atoms. A random Gaussian distribution of velocities is generated such that the total energy is distributed equally along the three Cartesian axes. 3RT is also supplied for translation/rotation for each molecule beyond the first. The total energy added is (3N-6M)*RT (M is the number of translation/rotation constraints, normally 1). Note that the MDIT command replaces any current velocities with random new ones and assumes that the structure is completely minimized (if it is not, the resulting temperature will be higher than that given in this command unless temperature control is used).

If no initial temperature is specified, the simulation keeps the kinetic energy from a previous run. If a structure is fully minimized and no initial temperature (or initial energy) is specified, the simulation will not run. If a structure is not fully minimized, the difference between the current strain energy and the minimum energy builds up in the kinetic energy term.

The standard dynamics protocol is to fully minimize a structure, then start the simulation either with the desired temperature, or slowly increase the kinetic energy (set an initial (MDIT) and warm the molecule over ca 5 ps (MDYN args 7 and 8)) prior to an equilibration run. If an equilibration period is used, the actual run will occur in a separate step without specification of new initial or final temperatures.

If a precise (harmonic) temperature is desired, use the MDYN arg7 to set constant temperature dynamics. This will insure that the final temperature attained will be that desired regardless of the relative energy of the starting structure and any nearby (local) minima after equilibration.

MDFT — Molecular Dynamics Final Temperature

Final Temperature (K, Absolute) to be achieved at end of run. The difference between the original bath temperature (MDYN arg7) and the final desired temperature (MDFT arg5) is slowly subtracted from the bath temperature over the course of the simulation. This temperature changing mechanism can be used to either heat or cool systems linearly during the simulation. The actual final temperature will be close but not equal to arg5. If a precise final temperature is desired, a subsequent MDYN run with the desired temperature in MDYN arg7 should be used. MDYN arg7 and arg8 can also be used to heat or cool molecules over the course of a simulation. This alternative would be effected by having an initial temperature as desired and the final temp as MDYN arg7 and weak coupling to the temperature bath via MDYN arg8.

MDMC — Molecular Dynamics / Monte Carlo mixed mode

MCSD — Monte Carlo / Stochastic Dynamics mixed mode

These commands are synonyms; MCSD is the preferred form, since its name reminds the user that the mixed-mode procedure should always be carried out using stochastic dynamics. This is because stochastic dynamics and the Monte Carlo procedure both generate the canonical ensemble, whereas (constant-energy) molecular dynamics generates the microcanonical ensemble.

This procedure mixes stochastic dynamics and Monte Carlo (MC) internal coordinate movements. This makes it possible to explore conformational (phase, configurational) space more effectively than dynamics alone, while using the dynamics paradigm to obtain free energies. Internal coordinate Monte Carlo movements must be specified using MOLS and TORS commands as is done in MCMM conformational searching. TORS commands involving ring bonds and RCA4 commands are allowed with MCSD, but the acceptance ratio is very low. If this combination of commands is used, DEBG 104 should be given.

This command should only be used with stochastic dynamics (MDYN arg3=1 and arg7 !=0).

We do not recommend using SHAKE (controlled via arg2 of MDYN) during MC/SD simulations, as its use can result in less than optimal temperature control. Using MCSD without SHAKE sometimes requires a smaller simulation timestep than might otherwise be used. But since MCSD is very efficient at exploring conformational space relative to dynamics alone, this should not be an issue.

If the MC/SD simulation is unstable and fails, consider turning off long-range constant derivatives, by setting EXNB arg1 to 1, and increasing the non-bonded cutoffs to a large number, such as 100 Å.

MDVE — Molecular Dynamics VElocity file

This command directs MacroModel to write or read a velocity file to or from disk. The file has double precision x,y,z velocities for each atom. The filename is out-filename.vel when writing a velocity file and in-filename.vel when reading one. This option is useful for listing the velocities, for saving the velocities to continue a molecular dynamics run at some later time and as a mechanism by which externally generated velocities may be used as molecular dynamics starting points. Velocity files are stored as formatted disk files. Each line contains x,y,z velocities (meters/sec) for each substructure atom. The Fortran format for each line in the velocity data in the file is (3D20.12).

A velocity file read (arg1=1) must come after READ or SUBS commands and before any MDYN command. A velocity file write (arg1=0) must come after the MDYN command.

MDYN does not write final structures to the output file unless the WRIT command follows it. We generally save both final velocities and structures at the end of a molecular dynamics run.

MDFR — Molecular Dynamics update FRequency

Sets the frequency in picoseconds for updates of the nonbonded array.

For molecules where large movements of groups are expected, the nonbonded array should be periodically updated (e.g., every 0.2-0.5 ps). A better alternative is to use very long cutoffs (using EXNB command) and a force field like OPLS, which does not have explicit hydrogen bonding potentials. The more frequently the nonbonded array is updated, the slower the dynamics will run.

MDRE — Molecular Dynamics REset

Reset or reinitialize all molecular dynamics variables, set initial velocities to zero.

IMPS — perform IMPortance Sampling

Note: IMPS affects the course of MCLO and MCSD simulations. For the full description of possible combinations and what they do, see the MacroModel Technical Manual.

This command activates the importance sampling option. Use of this command implies that the input file contains the output from a conformational search. When this command is encountered, all the structures in the input file are read and the values of the specified torsions and bond angles calculated for each conformation. These are used later in the importance sampling. The energy and number of times the structure was found in the conformational search are also read from the input file. These quantities may be used to weight the importance sampling (this last option is not recommended, see below).

Importance sampling also activates calculation of conformational populations which can be used to obtain relative free energies (see MDSA command).

Importance sampling can be used either as a replacement or in conjunction with Metropolis Monte Carlo (MMC) in both Monte Carlo and mixed mode simulations. The method is described in detail in the technical manual and examples of command files are provided in MacroModel Dynamics Calculations.

IMPO — IMPortance sampling Options

This command controls options available for the IMPS procedure.

At the beginning of an IMPS run, the program checks the internal coordinates by making all possible n² conformational interconversions between the n input conformers. The program produces an n by n matrix whose (i,j) element gives the difference between the energy of conformer j obtained from conformer i and its actual energy as read from the input file. Ideally, these energetic differences should all be zero, but in practice they are nonzero, since the interconversion is accomplished by means of torsional and bond-angle transformations alone. Energetic differences larger than a predetermined threshold lead to the termination of the run, unless an override is specified. The n² structures obtained by these interconversions may be written to the output file. Superimposing high energy structures (those for which the energy difference exceeds the threshold) on their original counterparts may help in identifying the cause of large energy differences. In particular, the deletion of a crucial torsion or bond angle leads to these symptoms, and the procedure just described allows such missing degrees of freedom to be identified.

Arguments 1, 2, 3, and 5 control the process described above. Argument 6 controls the criterion for declaring a conformation encountered during a simulation to be an instance of one of the IMPS input conformations, and Argument 7 controls the action to be taken when an importance-sampling step is scheduled, but the starting conformation does not correspond to any IMPS input conformation.

ITBS — Improper Torsion energy BiaSing

This command is used to bias the energy of a conformational family with a proper or improper torsion (also with a proper torsion) in a specific range by a constant amount. If, for example, the steric binding energies of conformations of L-family ligands are much lower than those of the D-family ligands, then attempts to calculate the enantioselectivity by jumping between the L and D families is likely to fail, because most L → D jumps will be rejected. In such cases, it is possible to artificially stabilize the D family by subtracting from the steric energies of its members a constant amount. This number must then be added back to the free-energy difference at the end of the simulation. If, for example, an energy value of 1.5 kcal/mol is subtracted from the members of the D family during the run, leading to a D/L free-energy difference of 0.5 kcal/mol in favor of the latter, then the real free-energy difference is 0.5 + 1.5 = 2 kcal/mol.

IMCC — IMportance sampling Chirality Check

This command is used to define chiral centers to be checked when assigning the simulated structure to one of the input structures. The current structure will not be assigned to an input structure if the chirality of one of its centers is different from that of the input structure. This command must come before the IMPS command.

TCMP — Torsional CoMParison

Used during importance sampling to define the torsions for torsional RMS calculation and for the torsional check during the assignment of the simulated structure to one of the input structures. Can be used either with torsional RMS calculations (both functions) or with Cartesian RMS calculations (second function only).

ZMAT — Z-MATrix

Defines the internal degrees of freedom (Z-matrix variables) to be sampled during all degrees of freedom MMC or MC(JBW) runs (also used in MCSD and MC(JBW)/SD). In a Z-matrix, each atom (except the first three specified) is defined with a bond length, bond angle, and torsional angle with respect to atoms previously defined. The first atom is defined with respect to itself, the second by a bond length to the first, and the third by a bond length to the second and a bond angle to the first.

Also defined in this command are the ranges of change/randomization for each degree of freedom. Typical ranges are 0-0.1 Å for bond lengths and 5° for bond angles. In course of an MMC simulation, backbone torsions should be changed by large ranges to allow for a complete sampling of the entire potential energy surface. With the above definition of Z‑matrix, torsions to branched atoms should be changed by a small amount (0-5°) to avoid distorted bond angles. For MC(JBW) and MC(JBW)/SD simulations, typical randomization ranges are 0-0.1A for bond angles and 0-5 for bond and torsion angles. Large randomizations reduce the acceptance rate and the number of conformational interconversions, and have the effect of reducing the JBW procedure to an MMC one.