Manual ver.1.5.2384 ACEMD ver. 2384 or greater

Contents

General
Atomistic force fields
Input commands
Tcl Scripting
Plugin interface
Getting support
Appendices
References

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. Currently ACEMD runs on Windows 32, 64 bit, Vista, XP and Linux 64 operating systems. The basic version is for Linux only. In order to install ACEMD, unpack the software. 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.

The latest version of ACEMD is compiled with the Nvidia toolkit cuda3.1, therefore to run you need a driver which is cuda3.1 compatible (version 256.35 or greater). Download the new driver here: http://developer.nvidia.com/object/cuda_3_1_downloads.html

First set the ACEMD_HOME directory to the directory where the acemd binary is placed. For a bash shell

#replace /home/acellera/acemd with the complete path to your ACEMD directory
DIR=/home/acellera/acemd
export ACEMD_HOME=$DIR/bin
export PATH=${PATH}:$ACEMD_HOME
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${ACEMD_HOME}

You might like to add these commands to your .bashrc file.

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

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. Low CPU load mode

ACEMD runs most of the calculations on the GPU. The CPU is used to poll the termination of the kernel in the GPUs, for running TCL scripts and plugins. It is possible to change the syncronization with the GPU and reduce the load of the CPU to just few per cent by setting the shell variable SWAN_SYNC:

>export SWAN_SYNC=4
>acemd input

will run in low CPU load mode. Unset the variable to run normally. In low CPU load mode, ACEMD is just few per cent slower than with the CPU fully occupied.

1.7. 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:

CHARMM36a2 - this is a combined parameter file which uses the new lipid naming scheme used by CHARMM-GUI (cfr e.g. P1 instead of previously-used P for the phosporus atom in lipids). Note that matching, updated topology files are needed to generate PSFs. Requires CMAP and NBFIX support.
CHARMMC31b1
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. It is necessary to provide a PDB instead of the CRD file (load it and save it from VMD). Use the following commands in your ACEMD configuration files :

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
                     # use .ff99SBildn for D.E. Shaw force field 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

For more information on the D.E. Shaw "ILDN" force field, see this link.

Alternatively, but rather more complicated, 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:

AMBER99 for use in CHARMM: Par/Top files
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.pdb
structure {string}: name of the PSF file generated by VMD/NAMD containing the structural information of the molecular system.
Ex.: structure 1JN0.psf
bincoordinates {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.coor
binvelocities {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.vel
parameters {string}: CHARMM input parameter file name.
Ex.: parameters par_all27_prot_lipid.inp
parmfile {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 cause ACEMD to restart from restartfile once executed. DCD files will be automatically appended. With this command, it is very simple to increase the simulation length or restart a stopped one. Default is off
restartsave {on|off}: Set to on to append the iteration number to restart filenames. 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 restart
restartfreq {integer}: Number of iterations between dumps of restart files.
Ex.: restartfreq 10000
dcdfile {string}: Name of the dcd trajectory file for the atom coordinates.
Ex: dcdfile out.dcd
dcdfreq {integer}: Number of iterations between dumps of trajectory coordinates.
Ex: dcdfreq 25000

veldcdfile {string}: Name of the dcd trajectory file for the atom velocities.
Ex: veldcdfile out.vel.dcd
veldcdfreq {integer}: Number of iterations between dumps of trajectory velocities.
Ex: veldcdfreq 2500
energyfreq {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 1000
outputname {integer}: Name of the final output files. This will dump 2 files the coordinates (.coor) and velocities (.vel).
Ex: outputname output
wrap {none|all}: If set with wrap all, all the atom coordinates are wrapped around the periodic box in the output files. Water molecules are wrapped such that they are not split across the periodic boundary, so hydrogens are always wrapped close to the oxygen. Default is none.

3.3. Simulation parameters

celldimension { real real real}: Spatial dimension of the rectilinear periodic cell box. Note ACEMD supports only orthogonal boxes, the reason being that you can often get much better saving by restraining the rotational diffusion of the protein instead
Ex.: celldimension 66 66 66
extendedsystem {string}: Read from an xsc file the cell dimension of the periodic box overriding celldimension. From acemd ver. 2024.
Ex.: extendedsystem input.xsc
cutoff {real}: Length of the cut-off radius.
Ex: cutoff 12
switching {on|off}: Use switching to zero of the van der Walls potential.
Ex: switching on
switchdist {real}: Length of the switching distance.
Ex: switchdist 10
1-4scaling {real}: Scaling factory to apply to electrostatic interactions between 1-4 bonded atoms
Ex: 1-4scaling 1.0
exclude {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-4
amber {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|all}: Constrain bond lengths formed by hydrogen atoms with heavy atoms at the equilibrium length. "None" uses harmonic bonds, "all" to all hydrogen bonds. This is required for timesteps larger than 1 fs. (Default none).
Ex: rigidbonds all
timestep {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 123456
minimize {integer}: Number of minimization steps to perform. During minimization energies are plot at every time step and no dcd file is produced. (Default 0).
Ex: minimize 1000
run {integer}: Number of steps of the molecular dynamics simulation. This is the last command of the Tcl configuration file, commands after "run" are ignored. To execute complex protocols during the run, please look at the Tcl scripting section.
Ex: run 10000000

3.3.1. Langevin thermostat

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 better during equilibration.

langevin {on|off}: Use the Langevin thermostat to control temperature. (Default off).
Ex: langevin on
langevintemp {real}: Equilibrium temperature of the thermostat.
Ex: langevintemp 298
langevindamping {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

ACEMD implements a Berendsen barostat designed for the equilibration of molecular systems (globular and in a membrane) to then start NVT production runs. With the system sizes which are achievable nowadays it is not necessary to have a pressure control in the production run, unless you really know what you are doing (for large number of atoms all ensembles are equivalent statistically). For molecular systems up to 100,000 atoms in a membrane allow for an equilibration of 20 ns, for globular proteins 1 to 5 ns are sufficient.

berendsenpressure {on|off}: Use the Berendsen barostat to control the pressure (Default off).
Ex: berendsenpressure on
berendsenpressuretarget {real}: Set the target pressure in units of bar. The default is 1 atm. (Default 1.01325).
Ex: berendsenpressuretarget 1.01325
berendsenpressurerelaxationtime {real}: Set the relaxation time of equilibration of the pressure in units of fs. (Default 400).
Ex: berendsenpressurerelaxationtime 400
berendsenpressurecompressibility {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).
useconstantratio {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
computepressure {on|off}: Compute the pressure even if the barostat is off for output in tclforces. (Default off).

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 and the performance is maximum.

pme {on|off}: Select whether to use PME. (Default off).
pmegridspacing {real}: Set the PME grid automatically such that the bin size between grid points is at least pmegridspacing A. This is overridden by setting the explicit grid size. (Default 1.
Ex: pmegridspacing 1
fullelectfrequency {integer}: Set every how many time steps the PME calculations are performed. Suggested value is 2 always. It is likely that we will change the default in the next releases. (Default 1).
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.

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.

3.4. Harmonic positional constraints

Harmonic positional constraints can be applied on selected atoms. This is useful during equilibration but also in production runs, so ACEMD implements it in the fastest possible way: directly computed on the GPU. There are no limitations on the number of atoms to which the constraints are applied. Similar and more flexible constraints can also be applied using Tcl scripting.

constraints {on|off}: Apply harmonic positional constraints. (Default off).
consref {pdb filename}: Name of the PDB file from which to take the reference position for each atom and the harmonic constant in the beta column. A beta column with zero for one atom means no constraint. (For ACEMD ver 1862 and above) The occupancy column is used to specify which direction is restrained. With the following values: 0 for XYZ, 1 for X, 2 for Y, 3 for Z, 4 for XY, 5 for YZ and 6 for XZ.
constraintscaling {real}: Multiply the harmonic constant as read from the PDB file. Useful to increase or reduce the constraints without changing the PDB file. Differently from NAMD this cannot be interleaved with the run command. It must be used before it.

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. Tcl scripts are executed in the CPU, so they can be expensive if the number of atoms involved is large (depending on system size, but target for less than 100 atoms if possible). For simple harmonic positional constraints use the constraints command instead.

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 5
calcforces_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 list
readpdb {}: 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 list
loadcharges {}: Loads the charges from the system into an array. It must be called from within calcforces. Since version 2385.
Ex: loadmasses list
getstep: Returns the current simulation timestep value.
Ex: set step [ getstep ]
pressure: Returns the pressure tensor in the form of a 9 element list XX XY XZ YX YY YZ ZX ZY ZZ. It should be symmetric, but won't be because of numerical inaccuracy.
Ex: set P [ pressure ]
constraintscaling {real}: Multiply the harmonic constant as read from the PDB file. Useful in equilibration protocols to change the contraint forces within calcforces after some steps.

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. Getting support

Users can receive individual and confidential support at Acellera Ltd.
There is also a new ACEMD user mailing list. Email is acemd@googlegroups.com. Web interface is: http://groups.google.com/group/acemd

7. Appendices

7.1. 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.2. Configuration file example

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.3. Production run settings

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
fullelectfrequency 2
#...

#...
# 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.

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.

7.4. Equilibrium run settings

During equilibrium runs, like during the equilibration of the pressure, use a smaller timesteps to account for higher forces and a stronger Langevin damping:

#...
timestep 2
rigidbonds      all
hydrogenscale   1
#...

langevin on
langevindamping 1

7.5. Apply periodic boundary conditions

ACEMD by default does not wrap in the output the coordinates around the periodic box, so when the molecular system is visualized water molecules appear to diffuse away. This is done such that you can measure diffusion coefficients. It is just a visualization choice. In the newest versions of ACEMD there is the command "wrap all" which instead wraps the coordinates of every atoms.

Best and most flexible way to do the wrapping is after the simulation with VMD. From VMD 1.8.7 (http://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/) you can wrap using the information in the DCD unitcell (load the dcd directly on a structure file rather than the pdb):

pbc set [pbc get -all] -all
pbc wrap -center bb -centersel protein -compound res -all

NB: In the case that you want to align the protein as well to some structure, be careful to align only after wrapping, otherwise the wrapping is wrong as the periodic box is rotated by the alignment. With the older version of VMD, 1.8.6 use instead: pbc wrap -sel "not protein" -center "protein" -all

8. 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.