Molecular dynamics is used to simulate the the physical motion of atoms and molecules. It can for example be used to simulate how lipids of a cell membrane behave, how a fictitious material would perform, if a protein would denaturalize in water or how a liposome fuses with a membrane. The simulation itself is based on Newton's laws of motion, that define classical mechanics, and force fields. These laws are as follows:
Template:Team:Groningen/CONTENT/MODELING/Supporting material
Introduction to molecular dynamics
If no force is applied, the object continues at it current speed.
The acceleration of the object is equal to its mass and the sum of the forces applied to the object.
When an object influences another object then the force on the second object is equal to the exerted force, but in the opposite direction.
These laws by themselves are not enough to simulate the motions of atoms and molecules. The forces that are mention in these laws have to be defined. The forces are defined by force fields in molecular dynamics simulations. Depending on the contents of the simulation some may be better than others, because one could model electrostatic interactions well while others are good for studying phase separations. Often these force fields are based on experimental data and have certain trade offs, therefore even for something as common as water different force fields exist.
There are however a few limitations to keep in mind. First of all, because these simulations are build from small elements, the simulations grows very fast. This in turn means that the simulations can take a long time even though the systems dimensions are usually on the scale of nanometers. Not only can the systems take long because they contain a lot of elements, but the effect of interest happens after a long time (the nanoseconds and microseconds range). Secondly because classical molecular dynamics is based on Newton’s laws of motion, it does not have chemical reactions. Lastly the behavior the simulations exhibit strongly depend on the input and the settings. For example the initial position might cause excessive forces and thus crashes the program.
Scripts
All the python scripts are made for python 2.7. There are scripts for the following things:
PGA topology creation (<a href="scpga" class="link">link</a>)
Cellulose topology creation (<a href="sccell" class="link">link</a>)
Cellulose phosphate topology creation (<a href="sccellp" class="link">link</a>)
PGA topology creation
This script can be used to generate a PGA topology file.
- /usr/bin/python
namef = 'pga{}.itp'
num = 100
def outputLine(f, line):
def p(*args):
f.write(line.format(*args))
return p
def printAtoms(f, num):
f.write(';nr\ttype\tresnr\tresidu\tatom\tcgnr\tcharge\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{0}\t{5: d}\n')
pa(1, 'Qd', 1, 'PGATN', 'NCC', 1)
pa(2, 'Qa', 1, 'PGATN', 'OCO', -1)
pa(3, 'Na', 1, 'PGATN', 'CCO', 0)
for x in range(0, num - 2):
i = 3 * x + 4
pa(i, 'Nd', x + 2, 'PGAM', 'NCHC', 0)
pa(i + 1, 'Qa', x + 2, 'PGAM', 'OCO', -1)
pa(i + 2, 'Na', x + 2, 'PGAM', 'CCO', 0)
i = 3 * num - 2
pa(i, 'Nd', num, 'PGATO', 'NCHC', 0)
pa(i + 1, 'Qa', num, 'PGATO', 'OCO', -1)
pa(i + 2, 'Qa', num, 'PGATO', 'OCO', -1)
def printBonds(f, num):
f.write(";i\tj\tfunct\tlength\tforce.c.\n")
pb = outputLine(f, '{0}\t{1}\t1\t{2}\t{3}\n')
pb(1, 2, 0.23, 18000)
pb(1, 3, 0.33, 4000)
for x in range(0, num - 2):
i = 3 * x + 3
pb(i, i + 1, 0.28, 12000)
pb(i + 1, i + 2, 0.24, 22000)
pb(i + 1, i + 3, 0.32, 2000)
i = 3 * num - 2
pb(i - 1, i, 0.28, 12000)
pb(i, i + 1, 0.24, 18000)
pb(i, i + 2, 0.36, 4000)
def printAngles(f, num):
f.write(';i\tj\tk\tfunct\tangle\tforce.c.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t2\t{3}\t{4}\n')
pa(2, 1, 3, 111.013, 100)
pa(1, 3, 4, 141.237, 55)
for x in range(0, num - 2):
i = 3 * x + 3
# next one is the same as the first one of the terminal end
pa(i, i + 1, i + 2, 114.646, 25)
pa(i, i + 1, i + 3, 124.643, 250)
pa(i + 1, i + 3, i + 4, 133.794, 250)
pa(i + 2, i + 1, i + 3, 113.393, 40)
i = 3 * num - 3
pa(i, i + 1, i + 2, 114.646, 25)
pa(i, i + 1, i + 3, 127.488, 100)
pa(i + 2, i + 1, i + 3, 116.572, 45)
def main():
with open(namef.format(num), 'w') as f:
f.write('[ moleculetype ]\n;molname\tnrexcl\nPGA{}\t\t1\n'.format(num))
f.write('\n[ atoms ]\n')
printAtoms(f, num)
f.write('\n[ bonds ]\n')
printBonds(f, num)
f.write('\n[ angles ]\n')
printAngles(f, num)
if __name__ == "__main__":
main()
Cellulose topology creation
This script can be used to generate a cellulose topology file and a initial coordinates file.
- !/usr/bin/python
import math
namef = 'cellulose{}.{}'
name = 'CELL{}'
num = 2
monomer = {"atoms": [("P1", "CELL", "B1", 0.0, 44.0534),
("P2", "CELL", "B2", 0.0, 75.0442),
("P4", "CELL", "B3", 0.0, 60.0528),
("P2", "CELL", "B4", 0.0, 58.0368),
("P1", "CELL", "B5", 0.0, 44.0534),
("P4", "CELL", "B6", 0.0, 60.0528)],
"bonds": [(1, 2, 1, 0.242, 30000),
(2, 3, 1, 0.284, 30000),
(2, 4, 1, 0.518, 30000),
(4, 5, 1, 0.234, 30000),
(4, 6, 1, 0.278, 30000)],
"constraints": [],
"angles": [(1, 2, 4, 2, 126, 50),
(3, 2, 4, 2, 120, 50),
(5, 4, 2, 2, 60, 100),
(6, 4, 2, 2, 65, 25)],
"dihedrals": [(1, 2, 4, 5, 1, 30, 8, 1),
(1, 2, 4, 6, 1, -150, 5, 1),
(3, 2, 4, 5, 1, -150, 5, 1)]}
def outputLine(f, line):
def p(*args):
f.write(line.format(*args))
return p
def printAtoms(f, num):
f.write('\n[ atoms ]\n')
f.write(';nr\ttype\tresnr\tresidu\tatom\tcgnr\tcharge\tmass\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{0}\t{5: .4f}\t{6: .4f}\n')
monolen = len(monomer["atoms"])
atomcnt = num * monolen
atoms = num * monomer["atoms"]
resnr = 1
for x in range(0, atomcnt):
atype, res, atom, charge, mass = atoms[x]
pa(x + 1, atype, resnr, res, atom, charge, mass)
if x % monolen == monolen - 1:
resnr += 1
def printBonds(f, num):
f.write('\n[ bonds ]\n')
f.write(';i\tj\tfunct\tlength\tforce.c.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\n')
numatoms = len(monomer["atoms"])
numbonds = len(monomer["bonds"])
if numbonds > 0:
for x in range(0, num):
for bond in monomer["bonds"]:
i, j, funct, leng, forc = bond
i = i + x * numatoms
j = j + x * numatoms
pa(i, j, funct, leng, forc)
f.write(';inter monomer bonds\n')
for x in range(0, num - 1):
i, j, funct, leng, forc = monomer["bonds"][2]
i = i + 2 + x * numatoms
j = j + 4 + x * numatoms
pa(i, j, funct, leng, forc)
def printConstraints(f, num):
f.write('\n[ constraints ]\n')
f.write(";i\tj\tfunct\tlength\tforce.c.\n")
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\n')
numatoms = len(monomer["atoms"])
numconstraints = len(monomer["constraints"])
if numconstraints > 0:
for x in range(0, num):
for bond in monomer["constraints"]:
i, j, funct, leng, forc = bond
i = i + x * numatoms
j = j + x * numatoms
pa(i, j, funct, leng, forc)
f.write(';inter monomer constraints\n')
for x in range(0, num - 1):
i, j, funct, leng, forc = monomer["constraints"][2]
i = i + 2 + x * numatoms
j = j + 4 + x * numatoms
pa(i, j, funct, leng, forc)
def printAngles(f, num):
f.write('\n[ angles ]\n')
f.write(';i\tj\tk\tfunct\tangle\tforce.c.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{5}\n')
numatoms = len(monomer["atoms"])
for x in range(0, num):
for angle in monomer["angles"]:
i, j, k, funct, leng, forc = angle
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
pa(i, j, k, funct, leng, forc)
f.write(';backbone angles between monomers\n')
angles = [(2, 4, 8, 2, 168, 50),
(4, 8, 10, 2, 168, 50)]
for x in range(0, num - 1):
for angle in angles:
i, j, k, funct, leng, forc = angle
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
pa(i, j, k, funct, leng, forc)
def printDihedrals(f, num):
f.write('\n[ dihedrals ]\n')
f.write(';i\tj\tk\tl\tfunct\tangle\tforce.c.\tmultiplic.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\n')
numatoms = len(monomer["atoms"])
for x in range(0, num):
for dihedral in monomer["dihedrals"]:
i, j, k, l, funct, leng, forc, multi = dihedral
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
l = l + x * numatoms
pa(i, j, k, l, funct, leng, forc, multi)
f.write(';inter monomer dihedrals\n')
dihedrals = [(7, 8, 4, 5, 1, 30, 5, 1),
(7, 8, 4, 6, 1, -150, 5, 1),
(9, 8, 4, 5, 1, -150, 5, 1)]
for x in range(0, num - 1):
for dihedral in dihedrals:
i, j, k, l, funct, leng, forc, multi = dihedral
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
l = l + x * numatoms
pa(i, j, k, l, funct, leng, forc, multi)
def main():
with open(namef.format(num, 'itp'), 'w') as f:
f.write('[ moleculetype ]\n;molname\tnrexcl\n' + name.format(num) + '\t\t1\n')
printAtoms(f, num)
printBonds(f, num)
printConstraints(f, num)
printAngles(f, num)
printDihedrals(f, num)
with open(namef.format(num, 'gro'), 'w') as f:
f.write('cellulose\n{}\n'.format(num * len(monomer["atoms"])))
y, z = (1, 1)
for i in range(0, num):
x = 0.23 + i * 2 * 0.518 * math.cos(math.radians(6))
andx = 1 + i * 6
resnr = 1 + i * 2
x1 = x - 0.242 * math.cos(math.radians(48))
y1 = y + 0.242 * math.sin(math.radians(48))
x3 = x - 0.284 * math.sin(math.radians(54))
y3 = y - 0.284 * math.cos(math.radians(54))
x4 = x + 0.518 * math.cos(math.radians(6))
y4 = y + 0.518 * math.sin(math.radians(6))
a5 = 0.234 * math.cos(math.radians(60))
s5 = 0.234 * math.sin(math.radians(60))
x5 = x4 - a5 * math.cos(math.radians(6))
y5 = y4 - a5 * math.sin(math.radians(6)) + s5 * math.cos(math.radians(30))
z5 = z + s5 * math.sin(math.radians(30))
a6 = 0.278 * math.cos(math.radians(65))
s6 = 0.278 * math.sin(math.radians(65))
x6 = x4 - a6 * math.cos(math.radians(6))
y6 = y4 + a6 * math.sin(math.radians(6)) + s6 * math.cos(math.radians(-150))
z6 = z + s6 * math.sin(math.radians(-150))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELL", "B1", andx, x1, y1, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELL", "B2", andx + 1, x, y, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELL", "B3", andx + 2, x3, y3, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELL", "B4", andx + 3, x4, y4, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELL", "B5", andx + 4, x5, y5, z5))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELL", "B6", andx + 5, x6, y6, z6))
f.write("%10.5f%10.5f%10.5f\n" % (x * 2, y * 2, z * 2))
if __name__ == "__main__":
main()
Cellulose phosphate topology creation
This script can be used to generate a cellulose phosphate topology file and a initial coordinates file.
- !/usr/bin/python
import math
namef = 'cellulosephos{}.{}'
name = 'CELLP{}'
num = 10
monomer = {"atoms": [("Qa", "CELLP", "B1", -2.0, 44.0534),
("P2", "CELLP", "B2", 0.0, 75.0442),
("P4", "CELLP", "B3", 0.0, 60.0528),
("P2", "CELLP", "B4", 0.0, 58.0368),
("Qa", "CELLP", "B5", -2.0, 44.0534),
("P4", "CELLP", "B6", 0.0, 60.0528)],
"bonds": [(1, 2, 1, 0.242, 30000),
(2, 3, 1, 0.284, 30000),
(2, 4, 1, 0.518, 30000),
(4, 5, 1, 0.234, 30000),
(4, 6, 1, 0.278, 30000)],
"constraints": [],
"angles": [(1, 2, 4, 2, 126, 50),
(3, 2, 4, 2, 120, 50),
(5, 4, 2, 2, 60, 100),
(6, 4, 2, 2, 65, 25)],
"dihedrals": [(1, 2, 4, 5, 1, 30, 8, 1),
(1, 2, 4, 6, 1, -150, 5, 1),
(3, 2, 4, 5, 1, -150, 5, 1)]}
def outputLine(f, line):
def p(*args):
f.write(line.format(*args))
return p
def printAtoms(f, num):
f.write('\n[ atoms ]\n')
f.write(';nr\ttype\tresnr\tresidu\tatom\tcgnr\tcharge\tmass\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{0}\t{5: .4f}\t{6: .4f}\n')
monolen = len(monomer["atoms"])
atomcnt = num * monolen
atoms = num * monomer["atoms"]
resnr = 1
for x in range(0, atomcnt):
atype, res, atom, charge, mass = atoms[x]
pa(x + 1, atype, resnr, res, atom, charge, mass)
if x % monolen == monolen - 1:
resnr += 1
def printBonds(f, num):
f.write('\n[ bonds ]\n')
f.write(';i\tj\tfunct\tlength\tforce.c.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\n')
numatoms = len(monomer["atoms"])
numbonds = len(monomer["bonds"])
if numbonds > 0:
for x in range(0, num):
for bond in monomer["bonds"]:
i, j, funct, leng, forc = bond
i = i + x * numatoms
j = j + x * numatoms
pa(i, j, funct, leng, forc)
f.write(';inter monomer bonds\n')
for x in range(0, num - 1):
i, j, funct, leng, forc = monomer["bonds"][2]
i = i + 2 + x * numatoms
j = j + 4 + x * numatoms
pa(i, j, funct, leng, forc)
def printConstraints(f, num):
f.write('\n[ constraints ]\n')
f.write(";i\tj\tfunct\tlength\tforce.c.\n")
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\n')
numatoms = len(monomer["atoms"])
numconstraints = len(monomer["constraints"])
if numconstraints > 0:
for x in range(0, num):
for bond in monomer["constraints"]:
i, j, funct, leng, forc = bond
i = i + x * numatoms
j = j + x * numatoms
pa(i, j, funct, leng, forc)
f.write(';inter monomer constraints\n')
for x in range(0, num - 1):
i, j, funct, leng, forc = monomer["constraints"][2]
i = i + 2 + x * numatoms
j = j + 4 + x * numatoms
pa(i, j, funct, leng, forc)
def printAngles(f, num):
f.write('\n[ angles ]\n')
f.write(';i\tj\tk\tfunct\tangle\tforce.c.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{5}\n')
numatoms = len(monomer["atoms"])
for x in range(0, num):
for angle in monomer["angles"]:
i, j, k, funct, leng, forc = angle
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
pa(i, j, k, funct, leng, forc)
f.write(';backbone angles between monomers\n')
angles = [(2, 4, 8, 2, 168, 50),
(4, 8, 10, 2, 168, 50)]
for x in range(0, num - 1):
for angle in angles:
i, j, k, funct, leng, forc = angle
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
pa(i, j, k, funct, leng, forc)
def printDihedrals(f, num):
f.write('\n[ dihedrals ]\n')
f.write(';i\tj\tk\tl\tfunct\tangle\tforce.c.\tmultiplic.\n')
pa = outputLine(f, '{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\n')
numatoms = len(monomer["atoms"])
for x in range(0, num):
for dihedral in monomer["dihedrals"]:
i, j, k, l, funct, leng, forc, multi = dihedral
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
l = l + x * numatoms
pa(i, j, k, l, funct, leng, forc, multi)
f.write(';inter monomer dihedrals\n')
dihedrals = [(7, 8, 4, 5, 1, 30, 5, 1),
(7, 8, 4, 6, 1, -150, 5, 1),
(9, 8, 4, 5, 1, -150, 5, 1)]
for x in range(0, num - 1):
for dihedral in dihedrals:
i, j, k, l, funct, leng, forc, multi = dihedral
i = i + x * numatoms
j = j + x * numatoms
k = k + x * numatoms
l = l + x * numatoms
pa(i, j, k, l, funct, leng, forc, multi)
def main():
with open(namef.format(num, 'itp'), 'w') as f:
f.write('[ moleculetype ]\n;molname\tnrexcl\n' + name.format(num) + '\t\t1\n')
printAtoms(f, num)
printBonds(f, num)
printConstraints(f, num)
printAngles(f, num)
printDihedrals(f, num)
with open(namef.format(num, 'gro'), 'w') as f:
f.write('cellulose\n{}\n'.format(num * len(monomer["atoms"])))
y, z = (1, 1)
for i in range(0, num):
x = 0.23 + i * 2 * 0.518 * math.cos(math.radians(6))
andx = 1 + i * 6
resnr = 1 + i * 2
x1 = x - 0.242 * math.cos(math.radians(48))
y1 = y + 0.242 * math.sin(math.radians(48))
x3 = x - 0.284 * math.sin(math.radians(54))
y3 = y - 0.284 * math.cos(math.radians(54))
x4 = x + 0.518 * math.cos(math.radians(6))
y4 = y + 0.518 * math.sin(math.radians(6))
x5 = x4 - 0.234 * math.cos(math.radians(54))
y5 = y4 + 0.234 * math.sin(math.radians(54))
z5 = z + 0.234 * math.sin(math.radians(60))
x6 = x4 - 0.278 * math.sin(math.radians(19))
y6 = y4 - 0.278 * math.cos(math.radians(19))
z6 = z - 0.278 * math.sin(math.radians(30))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELLP", "B1", andx, x1, y1, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELLP", "B2", andx + 1, x, y, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr, "CELLP", "B3", andx + 2, x3, y3, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELLP", "B4", andx + 3, x4, y4, z))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELLP", "B5", andx + 4, x5, y5, z5))
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnr + 1, "CELLP", "B6", andx + 5, x6, y6, z6))
f.write("%10.5f%10.5f%10.5f\n" % (x * 2, y * 2, z * 2))
if __name__ == "__main__":
main()