GROMACS-2018分子模拟教程系列：入门与进阶

需积分: 5 46 浏览量更新于2024-06-14 收藏 12.05MB PDF 举报

本文档是一份针对GROMACS-2018分子模拟软件的系列教程——"ALiveCoMSTutorial"。该教程由Justin A. Lemkul撰写，发表于2019年1月2日，旨在帮助用户深入了解并掌握分子动力学（MD）模拟技术在研究各种分子系统中的应用。MD模拟作为一种流行的原子级行为研究工具，对于想要从事这一领域的人来说，不仅需要对所选软件如GROMACS有深入理解，包括命令行操作、选项配置以及文件格式，还必须具备一定的物理学、数学知识，特别是计算机编程技能，并能够适应命令行环境。教程内容涵盖了七个核心主题，这些主题旨在降低入门门槛，解决新用户在开始MD模拟时可能遇到的复杂性。每个教程都围绕实际操作和理论知识相结合，通过实例演示如何设置模拟参数、处理系统初始构型、执行动力学模拟、分析结果等关键步骤。例如，"FromProteinstoPerturbed Hamiltonians"部分可能着重于如何模拟蛋白质的动态行为，以及如何引入扰动来研究系统响应。此外，文档强调了在线资源的重要性，提供了GitHub上的维护版本（<https://github.com/jalemkul/gmx_tutorials_livecoms>），用户可以通过这个平台提交反馈、建议或参与到教程的改进过程中。通过这种方式，作者鼓励互动与学习共享，确保教程内容的实用性和时效性。 GROMACS-2018的ALiveCoMSTutorial是一份全面且实用的指南，旨在帮助读者熟悉MD模拟的基础，并通过一系列系统的教学，提升他们在GROMACS软件上的操作能力和对分子动力学原理的理解。对于想要进入MD领域或提高现有技能的科研人员来说，这份教程是不可或缺的学习资料。

A LiveCoMS Tutorial

emstep = 0.01

nsteps = 50000

These instructions tell the GROMACS program

mdrun

(dis-

cussed below) to use the steepest descents method for en-

ergy minimization, and to terminate the process if the magni-

tude of the potential energy gradient is 1000.0 kJ mol

-1

or smaller. The maximum step size along the gradient is 0.01

nm, and a maximum of 50000 steps are allowed. The default

value of

emstep

in GROMACS of 0.01 nm is used, but it is often

necessary to use a smaller maximum step (e.g. 0.002 nm)

for systems that have diﬃculty converging; the smaller step

size allows for a more thorough walk over the potential en-

ergy surface, whereas a larger step may miss a path to the

true minimum. It is also possible to allow the process to go

on indeﬁnitely, stopping only when convergence is reached

(due to numerical precision or by actually achieving

emtol

) by

setting

nsteps = -1

. Note that

emstep

has units of distance,

not time. Energy minimization is not a dynamical process;

time does not elapse between each step and there are no

velocities. Energy minimization steps are expressed in terms

of maximum displacement along the vector indicated by the

force.

The remainder of the

.mdp

ﬁle contains instructions for

how to compute nonbonded interactions:

ns tlist = 1

cutoff - sche me = Verlet

ns _type = grid

cou l o m b type = PME

rc oulomb = 1.0

rvdw = 1.0

pbc = xyz

The value of

nstlist

controls how many steps (minimization

or MD integration) elapse between updating the neighbor list

for determining which atoms contribute to the short-range

forces. In MD simulations, this value is larger because atoms

exchange from the neighbor list in diﬀusion-limited time, but

for energy minimization, it needs to be set to 1 because

the conﬁguration may change considerably between each

step. The

cutoff-scheme = Verlet

setting uses a buﬀered

neighbor list, that is, atoms outside the longest cutoﬀ are still

tracked to improve energy conservation. Neighbor searching

is done by checking atoms in neighboring grids (

ns_type =

grid

) rather than checking every possible atom (

ns_type =

simple

). All short-range nonbonded interactions (electrostat-

ics and van der Waals) are truncated at 1.0 nm (

rcoulomb =

1.0

and

rvdw= 1.0

), and periodic boundary conditions are

applied in all three dimensions (

pbc = xyz

). The PME method

is used to calculate long-range electrostatic forces.

Note that in

.mdp

ﬁles, there is no diﬀerence between

a hyphen (-) and underscore (_); hence in the above exam-

ple,

ns-type

and

ns_type

would be equivalent, as would

cutoff-scheme and cutoff_scheme.

As above, invoke

grompp

to assemble the instructions for

energy minimization (-f minim.mdp), coordinates

(

-c 1AKI_solv_ions.gro

), and topology (

-p topol.top

) to

create the run input ﬁle for energy minimization (

-o em.tpr

$ gmx gromp p -f minim . mdp -c

1 A K I _ s o l v _ i o n s . gro -p topol . top -o

em . tpr

The GROMACS

mdrun

program is responsible for perform-

ing all minimization and dynamics processes. To run energy

minimization, use the following syntax:

$ gmx mdrun -v - deffnm em

The

-v

option invokes "verbose" mode, in which an estimate

of time remaining is printed to the screen, along with the

current energy minimization step, the step size, potential

energy, and maximum force. It is not a required option, and if

printed to a ﬁle, may result in a large ﬁle that is saved to disk.

It is, however, instructive to watch the progress of energy

minimization to understand what is going on. The

-deffnm

option deﬁnes the base ﬁle name for all input and output

ﬁles and eliminates the need for using individual input and

output ﬂags to specify names.

mdrun

requires the

.tpr

ﬁle as its only input, normally

passed to the

-s

ﬂag.

mdrun

produces several ﬁle types, in-

cluding an ASCII text ﬁle with a

.log

extension, a binary ﬁle

with all energy values with a

.edr

extension, and trajectory

ﬁles with either

.trr

(full-precision coordinates, velocities,

and/or forces) or

.xtc

(reduced precision coordinates only)

extensions. When using

-deffnm

, these ﬁles will be called

with a corresponding ﬁle extension. It is useful to name ﬁles

in this manner to avoid relying on default ﬁle names (

md.log

ener.edr

traj.trr

traj_comp.xtc

), which will be the same

for any process carried out by mdrun.

When

mdrun

is done, the user will see something similar

to the following:

Stee p est D e scen t s co n verg e d to Fmax < 1000 i n 726 step s

Pote n tia l Ene r g y = -5. 8 4 4 8 7 69 e +05

Maxim u m force = 9.69 5 7593 e +02 on atom 736

Norm of f o r c e = 2 . 3290 9 28 e +01

For a system such as this, it is expected that the ﬁnal po-

tential energy will be a negative value, on the order of 10

- 10

kJ mol

-1

. Potential energy is an extrinsic quantity, so

larger systems will have larger magnitudes. For a solvated

protein, the value will be negative, as favorable electrostatic

interactions between water molecules dominate the poten-

tial energy terms. Smaller systems may even have positive

potential energy, arising from the fact that in such cases, the

intramolecular bonded terms have a larger magnitude than

9 of 53

https://doi.org/10.33011/livecoms.1.1.5068

Living J. Comp. Mol. Sci. 2019, 1(1), 5068

A LiveCoMS Tutorial

nonbonded terms. Bonded components of the potential en-

ergy function have positive values, by deﬁnition. So for small

systems, these may prevail and the ﬁnal, potential energy is

positive. Such an outcome is not necessarily an indication

that anything is wrong. It would, however, be unusual for

a fully solvated protein system such as this one to have a

positive potential energy.

The various energy terms saved during energy minimiza-

tion are stored in the

em.edr

ﬁle, and can be extracted with

the

energy

program. For example, to extract the potential

energy of a system as a function of energy minimization step:

$ gmx energ y -f em . edr -o potential . xvg

By default, GROMACS analysis programs write output for-

matted for use with the XmGrace plotting program, though

this formatting can be changed with the

-xvgr

option. The

output

potential.xvg

ﬁle can be plotted to produce Figure 5.

Figure 5.

Potential energy of the solvated lysozyme system as a

function of energy minimization step.

The solvated system is now at a reasonable energy min-

imum. Note that the steepest descents method is very eﬃ-

cient but also may have diﬃculty ﬁnding a true energy mini-

mum. The settings used in this tutorial are fairly permissive

but work reasonably well for most similar systems. There are

applications, however, that require more exhaustive energy

minimization (for instance, free energy calculations such as

those in section 3.7), either through multiple rounds of min-

imization with diﬀerent algorithms that may more carefully

sample energy minima, or with stricter convergence criteria.

It is also important to note that the large changes in poten-

tial energy are often confusing to new users. During dynamics,

it is expected that energy will be conserved, but the purpose

of energy minimization is to change the conﬁguration of the

system in a manner that leads to a decrease in energy. As

such, obtaining an average energy during energy minimiza-

tion is not a useful metric, nor is any examination of energy

conservation.

3.2.5 Equilibration

Following energy minimization, the system is now ready for

the start of dynamics. Equilibration is carried out to obtain

a stable thermodynamic ensemble for whatever conditions

are desired. Most biomolecular simulations are carried out

under a canonical (NVT) ensemble, in which the number of

particles (N), volume (V), and temperature (T) of the system

are conserved, or an isothermal-isobaric (NPT) ensemble, in

which the number of particles (N), pressure (P), and tempera-

ture (T) are conserved. There are other statistical mechanical

ensembles that can be employed, but this tutorial will focus

on those that are most relevant to biomolecular systems. Be-

fore carrying out a simulation under a given ensemble, the

system must be equilibrated under the same conditions. For

this tutorial, the production MD simulation will be carried out

under an NPT ensemble. It is possible to simply initiate an

equilibration simulation under the same conditions, but it is

often more robust to ﬁrst equilibrate the temperature of the

system under an NVT ensemble before applying a barostat

to control the temperature. The simultaneous combination

of velocity generation and coordinate scaling under the inﬂu-

ence of the barostat can introduce instabilities in a system

that may be far from equilibrium.

It is also important to note that rather than initiating the

simulation at the desired temperature, it is also possible to

slowly warm the system up. This technique will not be applied

in this tutorial, as there is some debate as to whether or not it

is necessary in most cases. Interested readers are directed to

the GROMACS manual section on "simulated annealing" for

information on how this process is carried out in GROMACS.

To initialize dynamics, velocities are generated by

grompp

and written to the

.tpr

ﬁle. The following settings in the

nvt.mdp ﬁle control the generation of velocities:

ge n_vel = yes

ge n_temp = 300

ge n_seed = -1

A new concept is the application of position restraints.

A restraint is a biasing potential that is applied to selected

atoms to disfavor motion. It should be noted that posi-

tion restraints do not prevent motion, rather they penalize

it. Position restraints are typically applied to solute non-

hydrogen atoms during equilibration to prevent deformations

that could be induced by the random nature of the starting

velocities and subsequent collisions between the solute and

10 of 53

https://doi.org/10.33011/livecoms.1.1.5068

Living J. Comp. Mol. Sci. 2019, 1(1), 5068

A LiveCoMS Tutorial

solvent. Position restraints are speciﬁed with the following

keyword in the nvt.mdp ﬁle:

define = - DPOSR ES

Recall the generation of the

posre.itp

ﬁle by

pdb2gmx

in sec-

tion 3.2.1, which is only invoked as part of an

#ifdef

block in

topol.top

. This construct means restraints are only applied

when the user wants them, with invocation being depen-

dent upon the use of the

define

keyword in the

.mdp

ﬁle.

The presence of

define = -DPOSRES

means that the

#ifdef

POSRES

condition is evaluated as true by

grompp

and the posi-

tion restraint topology (

posre.itp

) is applied to the system.

Note that the general syntax of

#ifdef KEYWORD

and

define

= -DKEYWORD

; failing to preﬁx the keyword with

-D

will not lead

to any error, but the #ifdef condition will not be matched.

Some of the settings in the

.mdp

ﬁle are the same as in

energy minimization and will not be repeated here. New

settings that are introduced during NVT are:

in t e g rator = md

nsteps = 50000

dt = 0.002

The above lines specify that this process is an MD simulation,

and the

setting speciﬁes the leap-frog integrator in GRO-

MACS. A total of 50,000 integration steps will be performed,

and each step corresponds to an integration time step,

t, of

0.002 ps, or 2 fs. The simulation will therefore last for 100 ps

(50,000 steps × 0.002 ps per step).

Output is speciﬁed in the next section of nvt.mdp:

ns txout = 500

ns tvout = 500

ns t energy = 500

nstlog = 500

The

nstxout

and

nstvout

settings control the interval of time

steps between writing coordinates (x) and velocities (v) to the

trajectory (

.trr

) ﬁle. With these settings, coordinates and

velocities are saved every 1 ps. Similarly, energy terms are

written to the

.edr

ﬁle (

nstenergy

) and

.log

ﬁle (

nstlog

) at

the same interval.

Bond constraints are applied to allow for the 2-fs integra-

tion time step speciﬁed above. The value of

t is limited by the

fastest vibrational modes in the system, which are typically

bonds to hydrogen atoms. As these bonds will rarely exist in

excited states, it is reasonable to model them as rigid, replac-

ing the harmonic oscillator term with a holonomic constraint.

Doing so allows the value of

t to be increased from roughly

0.5 fs up to 2 fs for typical dynamics runs. In the

.mdp

ﬁle, the

relevant settings are:

cont i n u a t i o n = no

constrain t _ a l g orithm = lincs

con s t r a ints = h- bonds

li n c s _iter = 1

lin c s _ o rder = 4

The

continuation

keyword tells

mdrun

that the coordinates

are not the output of a previously constrained simulation,

so the constraints need to be solved at the ﬁrst time step.

The method for applying constraints is LINCS, the Linear Con-

straint Solver [

]. Constraints are applied to bonds to

hydrogen atoms

constraints = h-bonds

, and the remaining

settings are the default LINCS parameters for accuracy of

the constraint equations. It is important to note that "con-

straints" and "restraints" are diﬀerent. As described above, a

restraint is a biasing potential that disfavors motion, whereas

a constraint is used to ﬁx a distance or angle in a simulation.

Older literature often uses these terms interchangeably, but

in modern simulation practice, these two terms have distinct

meaning.

The next relevant section is related to nonbonded interac-

tions, which were described above in section 3.2.4. After the

nonbonded settings are the keywords related to temperature

control:

tcoupl = V- re scale

tc - grps = Prot ein Non - Pr otein

tau_t = 0.1 0.1

ref_t = 300 300

The

tcoupl

keyword speciﬁes the algorithm that is used for

temperature coupling, also known as the thermostat. In this

example, the velocity rescaling method of Bussi et al. is ap-

plied [

]. Another thermostat commonly applied during

equilibration is the weak coupling method of Berendsen et

al. [

] (speciﬁed with

tcoupl = berendsen

in the

.mdp

ﬁle).

The Berendsen weak coupling method has fallen out of fa-

vor, as it has been demonstrated that it does not sample a

correct canonical velocity distribution [

]. The Bussi velocity

rescaling method, however, does sample a correct velocity

distribution, and has the same advantages as the Berendsen

method in that it quickly relaxes the system to the target

temperature, making it practical and eﬃcient for use in both

equilibration and production MD simulations. The

tc-grps

keyword speciﬁes which groups of atoms are coupled to the

speciﬁed thermostat. All atoms should be coupled, but they

can be separated into diﬀerent groups. Strictly speaking, any

separation of atoms into diﬀerent groups violates the Equipar-

tition Theorem, in that thermal energy can ﬂow into and out

of the thermal reservoir (thermostat) independently between

the diﬀerent groups, rather than being exchanged between

the atoms in the system. The practice of coupling solute and

solvent atoms separately arises from the "hot solvent/cold

solute" problem in MD simulations [

], in which the solvent

heats while the solute cools, due to a host of factors related

11 of 53

https://doi.org/10.33011/livecoms.1.1.5068

Living J. Comp. Mol. Sci. 2019, 1(1), 5068

剩余52页未读，继续阅读

FrontScience

粉丝: 0
资源: 74

GROMACS-2018分子模拟教程系列：入门与进阶

Tutorials.pdf

GROMACS使用教程.doc

Tutorial.pdf

2003-Tutorials-I.-Galton-Delta-Sigma-Fractional-N-Phase-Locked-Loops.pdf

tutorials-flow1-8.0b9.zip

gromacs-tutorials.pptx

Informix-Tutorials-2 Configuring-Informix.pdf

Tutorials - Application Timer .pdf

Informix-Tutorials-4 Managing-Informix-Logs.pdf

Informix-Tutorials-6 Connecting-to-Informix.pdf

最新资源