acemd | manual | plugins | resources | gallery | contact
v.1.3.0
G. De Fabritiis, I. Buch
Research Group of biomedical Informatics (GRIB-IMIM),
Universitat Pompeu Fabra, Barcelona Biomedical Research Park (PRBB),
C/ Dr. Aiguader 88, 08003, Barcelona, Spain.
g.defabritiis@gmail.com
Contents
1. General
1.1. Concept
The concept behind ACEMD is to create a new generation [ACEMD], simple and fast all-atom molecular dynamics (MD) simulation software using the same commands and input files of NAMD [NAMD]. The three design target features of ACEMD are:
- Speed
- Usability
- Support (via www.acellera.com)
ACEMD is fast, as it was designed from scratch to be an optimized MD engine that exploits the computational power of graphical processing units (GPUs). It is simple to use as it uses a subset of the syntax of one of the most user friendly MD packages around like VMD/NAMD. ACEMD provides the equivalent performance of parallel CPU simulations using many tens to hundreds of processors, depending on the number and type of GPU cards installed.
1.2. GPU hardware supported
We support all the latest Geforce graphics cards and Tesla units from Nvidia. Specifically, we suggest you have at least a GTX200 series cards or Tesla C1060 or S1060. We are now also supporting Nvidia Fermi-based cards, GTX470, GTX480 and Tesla C2050/C2070. Full list of compatible cards is available in the FAQ of the GPUGRID.net project.
1.3. How to install Nvidia drivers
Follow the simple instruction at the Nvidia website. If you want to make sure that you can run ACEMD on your hardware just follow the instruction for the GPUGRID.net project.
1.4. How to install ACEMD
ACEMD is distributed in a tar-compressed directory. You need a license to get the software. Currently it is available for Windows 32, 64 bit, Vista, XP and Linux 64 operating systems. In order to install ACEMD, unpackage the software. In Windows-like operating systems, just double click on the ACEMD-latest.tar.gz file. For Linux
$tar -xzvf ACEMD-latest.tar.gz
where ACEMD-latest.tar.gz is the ACEMD distribution. The ACEMD binary is located into the ACEMD/bin directory. The parallel version depends from the openmpi package, install it on your machine using yum for instance for Fedora distributions.
1.5. How to run
ACEMD is run with a single command line parameter, which is the path of the input file. For the case of running ACEMD from the same directory where the input file resides, type
acemd input.conf
If you have multiple NVIDIA devices in your machine, you can select device n by adding the --device n at the command line
acemd input.conf --device 1
If you posses the parallel version of the code, you can increase ACEMD speed by running it concurrently on 3 GPUs by a factor 2.3
mpirun -np 3 acemd --device 1 input
This command will start ACEMD on devices 1,2,3 (skipping device 0). Note you need mpirun installed on your system. Please contact support if you do not know how to enable mpirun on your system (see getting support).
The software prints a log file to standard output which includes the parameters used for the running simulation and the energies at the specified output frequency:
....
INFO ON PARAMETERS
....
# Step Bond Angle Dihed Impr Elec
VDW PE KE Total Efield Temp
0 9196.72 12058.89 3336.60 93.25 -77160.03
193.45 -52281.14 25807.91 -26473.23 0 298.13
1.6. Fundamental units
Ångström (Å=10-10meters) for lengths,
- kcal/mol for energies,
- Kelvin (K) for temperature,
- g/mol for mass.
- Derived units are:
time unit: 1[t] = 48.88821 femto seconds = 48.88821 10-15 seconds
- velocity unit: 1[v] = Å/[t]
- kB (Boltzmann constant): 0.001987191
2. Atomistic force fields
The default force-field formats used by ACEMD are CHARMM (with cross-term support) and Amber. The code can also simulate OPLS force-fields in CHARMM format via translated force-fields (see below).
2.1. CHARMM force field
The topology and parameters force fields for CHARMM are available at http://mackerell.umaryland.edu/CHARMM_ff_params.html. Please refer to this link for the complete and latest information. We also provide it here for convenience:
CHARMM27r (C27r, revised C27 for aliphatic chain of lipids): Par/Top files
CHARMM27rn (C27rn, force field for nitro compounds): Par/Top files
2.2. AMBER force field
For AMBER there are two options. If you are used to AMBER input files, then you can just keep doing that. Use in your ACEMD configuration files the following commands:
coordinates your_pdb_file amber on parmfile your_topology_file 1-4scaling 0.83333
We suggest that you keep switching on when using AMBER force fields in ACEMD to avoid force discontinuities at the cutoff distance as you do for CHARMM. Special parameters for the Amber forcefields can be found at http://www.pharmacy.manchester.ac.uk/bryce/amber
To generate a system with tleap using ff99SB and new ions, load the following settings
source leaprc.ff99SB # FF99 forcefield with Simmerling modifications source leaprc.gaff # If custom ligands/lipids: GAFF forcefield # Modified ions and ion names Joung/Cheatham ion parameters for TIP3P water loadamberparams frcmod.ionsjc_tip3p loadoff ions08.lib
To use the new parameter files for the D.E. Shaw "ILDN" force field, follow this link.
Alternatively, it is possible to use the common VMD psfgen utility. The AMBER force field in CHARMM format is available at http://www.glue.umd.edu/~jbklauda/research/download.html. When using this force field you should also set 1-4scaling 0.83333 in the input file. We provide the files here for convenience:
AMBER99SB for use in CHARMM (Hornak, V., et al. Proteins: S.F.B. 65:712 (2006)): Par/Top files
AMBER99-bs0 for use in CHARMM (Perez, A. et al. BJ. 92: 3817 (2007)): Par/Top files
2.3. OPLS force field
The use of OPLS force field parameters is possible using the version in CHARMM format, but only for the protein portion. See http://brooks.scripps.edu/charmm_docs/Data/oplsaa-toppar.tgz. This requires to use the configuration parameter "vdwgeometricsigma on".
3. Input commands
3.1. Input file parameters
coordinates {string}: name of the PDB file of the molecular structure.
Ex.: coordinate 1JN0.pdbstructure {string}: name of the PSF file generated by VMD/NAMD containing the structural information of the molecular system.
Ex.: structure 1JN0.psfbincoordinates {string}: name of the binary coordinate file generated by previous ACEMD runs or VMD/NAMD. Units for coordinates are in Ångström.
Ex.: bincoordinates restart.coorbinvelocities {string}: name of the binary coordinate file generated by previous ACEMD runs or VMD/NAMD. Units for velocities are in derived units Å/(48.88821 fs).
Ex.: binvelocities restart.velparameters {string}: CHARMM input parameter file name.
Ex.: parameters par_all27_prot_lipid.inpparmfile {string}: AMBER input parameter/topology file name. Only used if the command amber on is also used.
Ex.: parmfile protein.prmtop
3.2. Output file parameters
restart {on|off}: set to on to causes ACEMD to attempt to restart from restartfile at start-up. Default is off
restartsave {on|off}: set to on to append to restart filenames the iteration number. Default is off
restartname {string}: name of the restart files. This will dump 3 files: coordinates (.coor), velocities (.vel) and time step (.idx).
Ex.: restartname restartrestartfreq {integer}: Number of iterations between dumps of restart files.
Ex.: restartfreq 10000dcdfile {string}: Name of the dcd trajectory file for the atom coordinates.
Ex: dcdfile out.dcddcdfreq {integer}: Number of iterations between dumps of trajectory coordinates.
Ex: dcdfreq 2500veldcdfile {string}: Name of the dcd trajectory file for the atom velocities.
Ex: veldcdfile out.vel.dcdveldcdfreq {integer}: Number of iterations between dumps of trajectory velocities.
Ex: veldcdfreq 2500energyfreq {integer}: Number of iterations between prints of the energies of the system. (Default 1). Warning: too few iterations (low energyfreq values) may decrease overall performance of the code.
Ex: energyfreq 1000outputname {integer}: Name of the final output files. This will dump 2 files the coordinates (.coor) and velocities (.vel).
Ex: outputname output
3.3. Simulations parameters
celldimension { real real real}: Spatial dimension of the rectilinear periodic cell box.
Ex.: celldimension 66 66 66cutoff {real}: Length of the cut-off radius.
Ex: cutoff 12switching {on|off}: Use switching to zero of the van der Walls potential.
Ex: switching onswitchdist {real}: Length of the switching distance.
Ex: switchdist 101-4scaling {real}: Scaling factory to apply to electrostatic interactions between 1-4 bonded atoms
Ex: 1-4scaling 1.0exclude {none|1-2|1-3|1-4|scaled1-4}: Exclude from non-bonded interactions 1-2 or 1-3 or 1-4 bonded atoms. "none" computes all the non-bonded interactions. "scaled1-4" applies scales van der Walls parameters to 1-4 bonded atoms.
Ex: exclude scaled1-4amber {on|off}: Use amber topology file as input. Use parmfile to read the topology file. (Default off).
dielectric {real}: Set the relative dielectric constant. (Default 1).
vdwgeometricsigma {on|off}: Use a geometric mean to compute vdw sigma for two atom types instead of the arithmetic mean. (Default off).
hydrogenscale {integer}: Multiply the hydrogen mass by a factor. The total mass of the systems remains the same because mass is removed from heavy atoms bonded to the hydrogens. Suggested is factor 4 for which diffusion and viscosity of water are only marginally affected. Note that individual atom masses do not appear in the form of the equilibrium distribution, therefore changing them only affects the dynamic properties of the system (marginally) but not the equilibrium distribution. On the other hands, it allows longer time-steps. (Default 1).
Ex: hydrogenscale 4
rigidbonds {none|water|all}: Constrain bond lengths formed by hydrogen atoms with heavy atoms at the equilibrium length. None uses harmonic bonds. "water" applies rigid bonds only to water molecules. "all" to all hydrogen bonds. This is required for timesteps larger than 1 fs. (Default none).
Ex: rigidbonds alltimestep {integer}: Set the time step for the molecular dynamics simulation in fs. ACEMD supports a timestep of 4 fs used with hydrogenscale 4 and rigidbonds all. This is the suggested configuration.
seed {integer}: Initialize the random number generator. If not initialized, it will be set using system time.
Ex: seed 123456run {integer}: Number of steps of the molecular dynamics simulation.
Ex: run 10000000minimize {integer}: Number of minimization steps to perform. (Default 0).
Ex: minimize 1000
3.3.1. Thermostats
The Langevin thermostat is needed to keep the system in the NVT ensemble. This is the suggested ensemble for production runs. The langevindamping should be as small as possible in order to thermalize the system without affecting the transport parameters (diffusion). We suggest to use langevindamping 0.1 for all production runs in NVT. A langevindamping 1 is necessary when equilibrating also the pressure (see Barostat).
langevin {on|off}: Use the Langevin thermostat to control temperature. (Default off).
Ex: langevin onlangevintemp {real}: Equilibrium temperature of the thermostat.
Ex: langevintemp 298langevindamping {real}: Damping coefficient of the thermostat. Units are 1/ps. Suggested damping coefficient for production runs is 0.1.
Ex: langevindamping 0.1
3.3.2. Barostat (Testing!)
ACEMD implements a Berendsen barostat designed to be sufficient to equilibrate molecular systems (globular and in a membrane) to then start NVT production runs. NVT have to preferred for productions runs in any case as they introduce less noise into the dynamics. With the size of systems which are achievable nowadays it is not necessary to run the production runs with a pressure control unless you really know what you are doing (for large number of atoms all ensembles are equivalent statistically). Use langevindamping 1 to maintain the temperature during pressure equilibration. For molecular systems up to 100,000 atoms in a membrane allow for an equilibration of 20 ns. for globular proteins 5 ns are sufficient.
The Berendsen barostat is currently under tests and will be released soon. Current limitations are: rigidbonds have to be set to none, automatic restart does not work as the cell size is not available and it works only on one GPU.
berendsenpressure {on|off}: Use the Berendsen barostat to control the pressure (Default off).
Ex: berendsenpressure onberendsenpressuretarget {real}: Set the target pressure in units of bar. The default is 1 atm. (Default 1.01325).
Ex: berendsenpressuretarget 1.01325berendsenpressurerelaxationtime {real}: Set the relaxation time of equilibration of the pressure in units of fs. (Default 400).
Ex: berendsenpressurerelaxationtime 400berendsenpressurecompressibility {real}: Set the pressure compressibility for the barostat in units of 1/bar. Default value is set to water compressibility (Default 4.57E-5).
berendsenpressurefreq {integer}: Set the frequency of the rescaling of the volume in terms of time steps. It must be much more frequent than the relaxation time. (Default 1).
useflexiblecell {on|off}: If on it allows every dimension of the periodic cell to change independently. (Default off).
Ex: useflexiblecell on ;# for globular proteinsuseconstantratio {on|off}: If on it fixes the ratio of the XY dimensions of the periodic cell, but still allows them to change. The z dimension is changed independently. (Default off).
Ex: useconstantratio on; # for membranes
3.3.3. Particle-Mesh Ewald parameters
The particle-mesh-Ewald method is the most widely adopted methodology to resolve long range electrostatics. In ACEMD, it is suggested to be used with cutoff 9, and switchdist 7.5, and at least 1 grid point per A of the box size, e.g. for a 62 A box in the x direction use pmegridsizex 64. With these settings the approximation of the forces is good enough and the performance is maximum.
pme {on|off}: Select whether to use PME. Note, a cutoff less than 9.5 must be used with PME.(Default off).
pmegridsizex {integer}: The PME grid dimensions. (Default 0).
pmegridsizey {integer}: The PME grid dimensions. (Default 0).
pmegridsizez {integer}: The PME grid dimensions. (Default 0).
3.3.4. Generalized reaction field
The generalized reaction field method resolves long range electrostatics assuming an uniform medium faraway from the atom [Tironi]. It appears as correction shift of the electrostatic, so it is fast to compute but it then requires a cutoff of 12 at least. Overall, it is not faster than PME for ACEMD. This is an approximation, but it showed to be fine for proteins in solutions. Less so for membranes,for which PME should be used. I
grf {on|off}: Select whether to use generalized reaction field. Maximum cutoff 12. (Default off).
grfdielectric1 {real}: Set the dielectric within the cutoff radius. (Default 1).
grfdielectric2 {real}: Set the dielectric within the cutoff radius. (Default 0), where 0 means infinite dielectric which is the suggested choice to have no discontinuities at the cutoff.
4. Tcl Scripting
The entire input file is seen by ACEMD as a Tcl script. You can interleave Tcl command with the commands shown in this manual. Tcl is also useful to manipulate the molecular systems by reading coordinates, velocities and forces on-the-fly while the simulation is running.
ACEMD calls two tcl functions, calc_forces_init at startup and calc_forces at every step. Note that calc_forces is processed on the CPU.
4.1. Intrinsic tcl functions
tclforces {on|off}: Specifies usage of the Tcl interface. If it is set to on, the procedure calcforces is executed.
tclforcesfreq {integer}: Number of iterations between every execution of calcforces. (Default 1).
Ex: tclforcesfreq 5calcforces_init {}: Procedure where addatom and addgroup to be used in order to call atom coordinates within calcforces.
Ex: proc calcforces_init {} { ... }calcforces {}: Procedure where custom manipulation of the system must be specified.
Ex: proc calcforces {} { ... }addatom {integer}: Requests coordinates of the specified atom, and returns its atom ID. The ID may then be used to find coordinates and apply forces. Used within calcfoces_init.
Ex: set atom [addatom 42]addgroup {integer list}: Requests center of mass coordinates of the atoms list, and returns a group ID of the following form: gN (being N a small integer). This group ID may then be used to find coordinates and apply forces. Used within calcfoces_init.
Ex: set group [addgroup { 4 8 15 }]loadcoords {}: Loads the group coordinates from the selected atoms into an array. It must be called from within calcforces.
Ex: loadcoords listreadpdb {}: Loads the coordinates and beta value from a pdb file into an array. No need to be called within calcforces.
Ex: readpdb list "structure.pdb"addforce: Applies a given force in kcal·mol-1·Å-1 to the selected atom or group (i.e. group).
Ex: addforce $group { 5.2 6.0 7.1 }loadmasses {}: Loads the masses from the system into an array. It must be called from within calcforces.
Ex: loadmasses listgetstep: Returns the current simulation timestep value.
Ex: set step [ getstep ]
4.2. Extended tcl functions
veclength2 {v}: returns length to the power of two of vector v
veclength {v}: returns length of vector v
vecdot {x y}: returns cross product of vectors x y
vec_normalise {v}: returns the normal vector of v
get_groups { _idlist _grouplist }: builds a list of atom groups from a list of atom ids
loadsystem { _pdb value _id _coor }: loads a selection of atoms and its coordinates by the beta-column value in a PDB file loaded with the readpdb command
center_of_mass { _atoms _coords }: computes the center of mass of a list of atomic coordinates
restrain1_1D { coor O group Dir K d log }: restrains 1 group of atoms relative to a point in 1 dimension
restrain1_1D_flatbottom { coor O group Dir K d log }: restrains 1 group of atoms relative to a point in 1 dimension with a flatbottom potential of size 'd'
restrain1_3D { _coor O group K log }: restrains 1 group of atoms relative to a point in 3 dimensions
smd1 { _coor group refcoor k dir smd_vel }: Steered Molecular Dynamics. Displacement of an atomic selection on a certain direction at a certain velocity.
proc veclength2 {v} {
set retval 0
foreach term $v {
set retval [expr $retval + $term * $term]
}
return $retval
}
proc veclength {v} {
return [expr {sqrt([veclength2 $v])}]
}
proc vecdot {x y} {
if {[llength $x] != [llength $y]} {
error "vecdot needs vectors of the same size: $x : $y"
}
set ret 0
foreach t1 $x t2 $y {
set ret [expr $ret + $t1 * $t2]
}
return $ret
}
proc vec_normalise {v} {
set ret {}
set len [ veclength $v ]
foreach i $v {
lappend ret [ expr $i / $len ]
}
return $ret
}
proc get_groups { _idlist _grouplist } {
# USAGE: builds a list of atom groups from a list of atom ids
upvar $_idlist idlist $_grouplist grouplist
foreach i $idlist {
lappend grouplist [ addgroup $i ]
}
}
proc loadsystem { _pdb value _id _coor } {
# USAGE: loads a selection of atoms and its coordinates by the
# beta-column value in a PDB file loaded with the readpdb command
upvar $_pdb pdb $_id id $_coor coor
foreach atom [array names pdb] {
if {[lindex $pdb($atom) 3] == $value} {
lappend id [expr $atom - 1]
set c {}
lappend c [lindex $pdb($atom) 0]
lappend c [lindex $pdb($atom) 1]
lappend c [lindex $pdb($atom) 2]
lappend coor $c
}
}
}
proc center_of_mass { _atoms _coords } {
# USAGE: computes the center of mass of a list of atomic coordinates
upvar $_atoms atoms $_coords coords
loadmasses masses
set comsum { 0 0 0 }
set totalmass 0
for {set i 0} {$i<[llength $atoms]} {incr i} {
set atom [lindex $atoms $i]
set coo [lindex $coords $i]
set comsum [vecadd $comsum [vecscale $masses($atom) $coo]]
set totalmass [expr $totalmass + $masses($atom)]
}
return [ vecscale [expr 1.0/$totalmass] $comsum ]
}
proc restrain1_1D { _coor O group Dir K d log } {
# USAGE: restrains 1 group of atoms at at distance $d from $O in the direction $Dir with constant $K
upvar $_coor coor
global logfreq
set dr [expr [vecdot $coor($group) $Dir] - [vecdot $O $Dir] ]
set DR [expr $dr - $d]
set f [vecscale [expr -$K*$DR] $Dir]
addforce $group $f
set step [ getstep ]
if {$step % $logfreq == 0} {
set E [ expr 0.5*$K*$DR*$DR ]
print "$log $step $dr"
print "$log\_DR\_E $step $DR $E"
}
}
proc restrain1_1D_flatbottom { _coor O group Dir K d log } {
# USAGE: restrains 1 group of atoms with a flat bottom potential
# at a distance $d from $O in the direction $Dir with constant $K
upvar $_coor coor
global logfreq
set dr [expr [vecdot $coor($group) $Dir] - [vecdot $O $Dir] ]
set DR [expr abs($dr) - $d]
if {$DR > 0} {
if {$dr > 0} {
set DR [expr $dr - $d]
} else {
set DR [expr $dr + $d]
}
set f [vecscale [expr -$K*$DR] $Dir]
addforce $group $f
set step [ getstep ]
if {$step % $logfreq == 0} {
set E [ expr 0.5*$K*$DR*$DR ]
print "$log $step $dr $DR $E"
}
}
}
proc restrain1_3D { _coor O group K log } {
# USAGE: restrains 1 group of atoms relative to a point $O with constant $K
upvar $_coor coor
global logfreq
set dr [vecsub $coor($group) $O]
set f [vecscale [expr -$K] $dr]
addforce $group $f
set step [ getstep ]
if {$step == 0 || ($step % $logfreq == 0 && $log != "") } {
set DR [veclength $dr]
set E [ expr 0.5*$K*$DR*$DR ]
print "$log $step $DR $E"
}
}
proc smd1 { _coor group refcoor k dir smd_vel } {
# USAGE: Steered Molecular Dynamics. Displacement of 1 group of atoms
# on a certain direction $dir at a certain velocity $smd_vel with harmonic constant $k.
upvar $_coor coor
global logfreq
set step [ getstep ]
set cur_pos $coor($group)
set delta_r [ vecsub $cur_pos $refcoor ]
set diff [ vecdot $delta_r $dir ]
set c_1 [ expr $k*($smd_vel*$step - $diff) ]
set v_1 [ vecscale $c_1 $dir ]
addforce $group $v_1
if {$step % $logfreq==0} {
print "SMD $step $delta_r $v_1"
}
}
4.3. Examples of usage
4.3.1. Applying a force restrain to a single atom group in 1D
Note: Beware of adding to your input file the function definitions specified above
# Add function descriptions here
set structure mystructure.pdb
set Krestrain 10
set axis {1 0 0}
set logfreq 1000
tclforces on
proc calcforces_init {} {
global structure coor1 glist1 gcom1
# Extraction of atoms and coordinates from pdb file.
readpdb pdb $structure
# Loading of atoms with beta column equal to '1' and storage of
# selected atoms index values and coordinates.
loadsystem pdb 1 id1 coor1
# Creation of group corresponding to the center of mass of
# selected atoms to which we will apply the restrain.
set gcom1 [ addgroup $id1 ]
# Creation of a list of groups for every selected atom.
# We need this to obtain later the reference position for the center of mass.
get_groups id1 glist1
}
proc calcforces {} {
global coor1 glist1 gcom1 axis Krestrain logfreq
## Loading of system's current coordinates
loadcoords coords
## Getting atom selection coordinates
# Current
set gcom1_pos $coords($gcom1)
# Reference
set gcom1_start [ center_of_mass glist1 coor1 ]
## Applying a 1D restrain to the center of mass of the atom selection
## using as reference coordinates those extracted from the pdb.
restrain1_1D coords $gcom1_start $gcom1 $axis $Krestrain 0 "log_text"
return;
}
4.3.2. Applying a force restrain to a list of atom groups in 3D
Note: Beware of adding to your input file the function definitions specified above
# Add function descriptions here
set structure mystructure.pdb
set Krestrain 10
set logfreq 1000
tclforces on
proc calcforces_init {} {
global structure coor1 glist1
# Extraction of atoms and coordinates from pdb file.
readpdb pdb $structure
# Loading of atoms with beta column equal to '1' and storage of
# selected atoms index values and coordinates.
loadsystem pdb 1 id1 coor1
# Creation of a list of groups from the selected atoms
# to which we will apply the restrain.
get_groups id1 glist1
}
proc calcforces {} {
global coor1 glist1 Krestrain logfreq
## Loading of system's current coordinates
loadcoords coords
## Applying a 3D restrain to every atom in the selection
## using as reference coordinates those extracted from the pdb.
for {set i 0} {$i<[llength $glist1]} {incr i} {
set group [lindex $glist1 $i]
set initcoor [lindex $coor1 $i]
restrain1_3D coords $initcoor $group $Krestrain "log_text"
}
return;
}
5. Plugin interface
ACEMD is easily extended by adding plugin modules written in C and dynamically loaded by the application. Current plugins include metadynamics, a power biased free energy calculation method. It is easy to think of ways on which users might want to customize ACEMD for their needs. In practice, the plugin interface give access to the position,velocities and forces at each iteration from a C interface. Please visit the plugin web page for more information.
6. Tabular potentials
User defined potentials may be supplied for the electrostatic and VdW potentials. The command for enabling these is:
tabular {elec|vdw} {filename}
The file must have the following formats. For electrostatic:
r=<cutoff> r_0 U(r_0) dU(r_0)/dr ... r_N U(r_N) dU(r_N)/dr
For VdW:
r=<cutoff> r_0 C1(r_0) dC1(r_0)/dr C2(r_0) dC2(r_0)/dr ... r_N C1(r_N) dC1(r_N)/dr C2(r_N) dC2(r_N)/dr
The table must always have 4096 entries. The bin width is cutoff/4096 and the values of r are bin-centred. Tabular #potentials cannot be used in conjunction with GRF or PME. A runtime check is made to compare the cutoff valued defined in the simulation input with that specified in the tabular file.
7. Examples
7.1. A globular protein solvated in water
This demo is contained in the directory demo/dhfr (dihydrofolate reductase). This protein represents a common benchmark for molecular dynamics codes [Gromacs08]. The input files are provided here to allow users to benchmark the nanoseconds per day deliverd by ACEMD, to compare to other architectures, GPU versions and codes. Example input configuration file:
# CHARMM/AMBER joint benchmark # initial config coordinates 5dhfr_cube.pdb temperature 300 # integrator params timestep 4 hydrogenscale 4 rigidbonds all # force field params structure 5dhfr_cube.psf parameters par_all22_prot.inp exclude scaled1-4 1-4scaling 1.0 switching on switchdist 7.5 cutoff 9 # output params outputname dhfr_out dcdfreq 1000 # periodic cell celldimension 62.23 62.23 62.23 # full electrostatics pme on pmegridsizex 64 pmegridsizey 64 pmegridsizez 64 langevin on langevintemp 300 langevindamping 0.1 energyfreq 500 run 1500
7.2. Steered molecular dynamics on k+ ions for a protein pore embedded in a membrane
This demo is contained in the directory demo/ga_smd. It provides the full input files and input script to perform steered molecular dynamics simulations pulling a potassium ion within the Gramicidin A channel embedded into a lipid bilayer and solvated in water and ions at 150mM (taken from [DeFabritiis08GA]). This demo provides an example on how to use ACEMD Tcl interface for steered molecular dynamics simulations.
8. Performance tips
The most important production parameters are the one related to the cutoff. Please use the following settings in productions runs, unless you have specific reasons not to:
#... # WITH PME timestep 4 rigidbonds all hydrogenscale 4 switching on switchdist 7.5 cutoff 9 #...
#... # WITHOUT PME timestep 4 rigidbonds all hydrogenscale 4 switching on switchdist 10 cutoff 12 #...
These parameters have been tested on thousands of run and give you the maximum performance,accuracy and stability. Also ACEMD is optimized for these two cases. If you use PME, you cannot use cutoff larger that 9.5 due to the fixed box size needed to optimize occupancy of the GPU.
For the langevin thermostat, it is suggested to use a very low damping frequency
langevin on langevindamping 0.1
This will be sufficient to thermalize the system without reducing the diffusion coefficient of the solvent.
Outputing energies and coordinates is expensive, as the data must be moved in and out the GPU. For production runs, use frequency of the order of thousands.
9. Getting support
Users can receive individual and confidential support at Acellera Ltd.
10. References
[ACEMD] M. Harvey, G. Giupponi and G. De Fabritiis, ACEMD: Accelerated molecular dynamics simulations in the microseconds timescale, J. Chem. Theory and Comput. 5, 1632 (2009).
[ACEMDPME] M. J. Harvey and G. De Fabritiis, An implementation of the smooth particle-mesh Ewald (PME) method on GPU hardware, J. Chem. Theory Comput., 5, 2371–2377 (2009).
[DeFabritiis08GA] G. De Fabritiis, P. V. Coveney and J. Villa-Freixa, Energetics of K+ permeability through Gramicidin A by forward-reverse steered molecular dynamics, 73, 184 Proteins (2008).
[Hardy07] David J. Hardy, NAMD-Lite, http://www.ks.uiuc.edu/Development/MDTools/namdlite/, University of Illinois at Urbana-Champaign, 2007.
[hydrogenscale] Feenstra, K. A., Hess, B., Berendsen, H. J. C., Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems, J. Comp. Chem. 20, 786(1999).
[M-SHAKE]V. Krautler, W. F. Van Gunsteren, P. H. Hunenberger, A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations, J. Comp. Chem. 22, 501 (2001).
[NAMD] James C. Phillips, Rosemary Braun, Wei Wang, James Gumbart, Emad Tajkhorshid, Elizabeth Villa, Christophe Chipot, Robert D. Skeel, Laxmikant Kale, and Klaus Schulten. Scalable molecular dynamics with NAMD. J. Comp. Chem., 26, 1781-1802 (2005).
[PME] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth Particle Mesh Ewald Method, J. Chem. Phys. 103, 8577 (1995).
[RATTLE] H. C. Andersen, Rattle: A velocity version of the shake algorithm for molecular dynamics calculations, J. Comp. Phys. 52, 24 (1983).
[Tironi95] Ilario G. Tironi, Rene Sperb, Paul E. Smith, and Wilfred F. van Gunsteren, A generalized reaction field method for molecular dynamics simulations, J. Chem. Phys. 102, 5451 (1995).
[VMD] Humphrey, W., Dalke, A. and Schulten, K., "VMD - Visual Molecular Dynamics", J. Molec. Graphics, 1996, vol. 14, pp. 33-38.