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:
Difference between revisions of "Template:Team:Groningen/CONTENT/MODELING/Supporting material"
(added links) |
(fix links) |
||
Line 29: | Line 29: | ||
<div class="listing"> | <div class="listing"> | ||
− | <div class="item">PGA topology creation (<a href="scpga" class="link">link</a>)</div> | + | <div class="item">PGA topology creation (<a href="#scpga" class="link">link</a>)</div> |
− | <div class="item">Cellulose topology creation (<a href="sccell" class="link">link</a>)</div> | + | <div class="item">Cellulose topology creation (<a href="#sccell" class="link">link</a>)</div> |
− | <div class="item">Cellulose phosphate topology creation (<a href="sccellp" class="link">link</a>)</div> | + | <div class="item">Cellulose phosphate topology creation (<a href="#sccellp" class="link">link</a>)</div> |
</div> | </div> | ||
Revision as of 00:08, 19 September 2015
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.
All the python scripts are made for python 2.7. There are scripts for the following things:
- /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()
- !/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()
- !/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()