acemd | manual | plugins | resources | gallery | contact


newACEMDlogo.png
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

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:

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

  • 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 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 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 2500

  • 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

3.3. Simulations parameters

  • celldimension { real real real}: Spatial dimension of the rectilinear periodic cell box.
    Ex.: celldimension 66 66 66

  • 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|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 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

  • run {integer}: Number of steps of the molecular dynamics simulation.
    Ex: run 10000000

  • minimize {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 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 (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 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).
    Ex: useflexiblecell on ;# for globular proteins

  • 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

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

  • getstep: 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.

Copyright 2008-2009. All rights reserved.