Team:Groningen/Modeling/Supporting material
Blue Bio Energy
Supporting material
Introduction to molecular dynamics
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:
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.
Visualization with VMD
To view the results of a simulation VMD is a great tool (version 1.9.2 was used in this text). While it has many advanced features, not many need to be known to effectively use it. Before one can use VMD one has to load the system topology and optionally the trajectory in it. This is usually done with the following command: “vmd -f topology trajectory”. With the “-f” option in the command one specifies that the following files belong together. Now that the topology and trajectory is loaded, one can visualize the data. Sometimes it is also useful to see the periodic boundary conditions. To show them one has to enter “pbc box” into VMD's console and to hide them “pbc box -off”.
Firstly one can change how groups of atoms/beads (selections) are represented. The current rules are shown in the Graphics -> Representations menu. The most useful drawing styles with the MARTINI systems are the VDW and the QuickSurf options. With VDW VMD shows the van der waals radii for each atom/bead and with QuickSurf shows the surface of the selection. The selected atoms/beads can be changed by typing a selection “string” in the Selected Atoms field or by clicking in the Selection tab. For more details one could consult the documentation of VMD or do some tutorials. Some starting keywords are “resname”, “index”, “name” and “within of”.
Secondly it is possible to orientate the system differently. This is done by dragging on the output screen of VMD after pressing one of the hotkeys mentioned in Table 1.
Key
Description
r
rotate the system
t
translate the system
s
scale the system
=
reset the system orientation
Hotkeys related to the orientation of the system. When VMD is started, it is in the rotation mode.
Thirdly it is possible to label particles and measure things. An overview of some useful hotkeys are given in Table 2. The information VMD shows in these cases can be changed in the Graphics -> Labels menu. For example in the properties tab of a label one can change the format to “%1i”, so that VMD shows the atom index starts from one, like the numbering in the Gromacs topology files. Also note that VMD shows distances in ångström.
Key
Description
1
label the atom/bead
2
measure the distance between two particles
3
measure the angle between three particles
4
measure the dihedral between four particles
Hotkeys related to the label and measurement modes.
Lastly sometimes it is useful to have the bond information in VMD, so that it can show the bonds between the coarse grained molecules. This can be done with the cg_bonds script provided on the MARTINI homepage (link, 2015).
Sometimes one needs better looking images than just screenshots from the program. This is possible with the rendering options in VMD. One has to edit the representation, change some options and then one can options to get a nice picture.
In general there are a three things that make the final image look a lot better. Firstly enabling shadows and ambient occlusion in the ray tracing options in the Display -> Display Settings menu. Secondly using a nice material for the selections like ''AOChalky''. And lastly using a ray tracer to render the image. VMD comes with the Tachyon ray tracer. This renderer can be selected in the the File -> Render menu, before rendering the picture by hitting the start rendering button in the same menu. Note that movies can also be made in a similar matter with the movie maker extension.
Scripts
All the python scripts are made for python 2.7. There are scripts for the following things:
PGA topology creation (link)
Cellulose topology creation (link)
Cellulose phosphate topology creation (link)
Changing the concentration and preparing the membranes (link)
Calculating the flow of ions through the membranes (link)
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()
Changing ion concentrations in membranes
This script can be used to prepare a system with two membranes for simualtions.
#!/usr/bin/python2
import sys
import argparse
import re
def calcAtomNumConc(conc, vol):
return int(conc*vol*0.6)
def minmaxz(sel, box):
zmin, zmax = (box[2], 0)
for elem in sel:
z = elem[6]
zmin = min(zmin, z)
zmax = max(zmax, z)
return (zmin, zmax)
def calcConc(sel):
vol = calcVol(sel)
na, cl = countNaCl(sel)
print("conc: {} NA+ and {} CL-".format(na/(vol*0.6), cl/(vol*0.6)))
def ispoly(elem):
resname = elem[1]
return resname in ["PGATN", "PGATO", "PGAM", "CELL", "CELLP"]
def isw(elem):
resname = elem[1]
return resname in ["W", "PW"]
def change_ion_to_w(sel, diff, data):
for elem in sel[0:diff]:
elem[1] = "W"
elem[2] = "W"
def change_w_to_ion(sel, diff, name, data):
for elem in sel[0:diff]:
elem[1] = "ION"
elem[2] = name
def change_ion_to_pw(sel, diff, data):
ions = sel[0:diff]
for ion in ions:
i = data.index(ion)
ion[1] = "PW"
ion[2] = "W"
if data[i + 1][1] == data[i + 2][1] == "RM":
data[i + 1][1] = "PW"
data[i + 2][1] = "PW"
else:
data.insert(i + 1, [ion[0], "PW", "WP", 0, ion[4] + 0.06, ion[5] + 0.09, ion[6] + 0.09])
data.insert(i + 2, [ion[0], "PW", "WM", 0, ion[4] + 0.07, ion[5] + 0.07, ion[6] + 0.10])
def change_pw_to_ion(sel, diff, name, data):
replaced = 0
for i in range(2, len(sel)):
if replaced == diff:
break
if sel[i][0] == sel[i - 1][0] == sel[i - 2][0]:
ndx = data.index(sel[i])
data[ndx - 2][1] = "RM"
data[ndx - 1][1] = "RM"
sel[i][1] = "ION"
sel[i][2] = name
replaced += 1
def countNaCl(sel):
na = len([x for x in sel if x[2] == "NA+"])
cl = len([x for x in sel if x[2] == "CL-"])
return (na, cl)
def changeTopology(data, box, chItoW, chWtoI, midConc, outerConc):
# select the ions between the two bounding boxes of the membranes
mid = box[2] / 2
membba = minmaxz([x for x in data if ispoly(x) and x[6] < mid], box)
membbb = minmaxz([x for x in data if ispoly(x) and x[6] > mid], box)
volmid = box[0] * box[1] * (membbb[0] - membba[1])
volouter = box[0] * box[1] * (membba[0] + box[2] - membbb[1])
# remove ions between mid membrane to bounding box of pga
memmida = (membba[1] - membba[0]) / 2 + membba[0]
memmidb = (membbb[1] - membbb[0]) / 2 + membbb[0]
ionsa = [x for x in data if x[1] == "ION" and memmida < x[6] and x[6] < membba[1]]
ionsb = [x for x in data if x[1] == "ION" and membbb[0] < x[6] and x[6] < memmidb]
rna, rcl = countNaCl(ionsa)
na, cl = countNaCl(ionsb)
rna += na
rcl += cl
print("removed {} na and {} cl from membrane\n".format(rna, rcl))
chItoW(ionsa, len(ionsa), data)
chItoW(ionsb, len(ionsb), data)
# function to correct the difference atom numbers
def correctdiff(diff, name, sel):
if diff > 0:
print("removing {} {}".format(diff, name))
chItoW([x for x in sel if x[2] == name], diff, data)
elif diff < 0:
print("adding {} {}".format(-diff, name))
chWtoI([x for x in sel if isw(x) or x[1] == "RM"], -diff, name, data)
# change the concentration in the middle
selmid = [x for x in data if membba[1] < x[6] and x[6] < membbb[0]]
na, cl = countNaCl(selmid)
target = calcAtomNumConc(midConc, volmid)
print("\nmiddle target ion number: {}".format(target))
print("middle has {} na and {} cl".format(na, cl))
cldiff = cl - target
nadiff = na - target
correctdiff(cldiff, "CL-", selmid)
correctdiff(nadiff, "NA+", selmid)
rna += nadiff
rcl += cldiff
# change the concentration in the outer part
selouter = [x for x in data if x[6] < membba[0] or membbb[1] < x[6]]
na, cl = countNaCl(selouter)
target = calcAtomNumConc(outerConc, volouter)
print("\nouter target ion number: {}".format(target))
print("outer has {} na and {} cl".format(na, cl))
cldiff = cl - target
nadiff = na - target
correctdiff(cldiff, "CL-", selouter)
correctdiff(nadiff, "NA+", selouter)
rna += nadiff
rcl += cldiff
# correct charge
diff = rcl - rna
print("\ncorrecting charge with {} na".format(diff))
hdiff = int(diff / 2)
correctdiff(hdiff, "NA+", [x for x in selouter if x[6] < mid])
correctdiff(diff - hdiff, "NA+", [x for x in selouter if x[6] > mid])
data[:] = [x for x in data if not x[1] == "RM"]
def readGro(file):
alist = []
linere = re.compile("(.{5})(.{5})(.{5})(.{5})(.{8})(.{8})(.{8})(?:(.{8})(.{8})(.{8}))?$")
floatre = '(\S+)'
sre = floatre + "\s*" + floatre + "\s*" + floatre
boxre = re.compile("\s*" + sre + "(?:\s*" + sre + "\s*" + sre + ")?$")
with open(file) as file:
name = file.next()
numatoms = int(file.next())
for i in range(numatoms):
line = file.next()
data = list(linere.match(line).groups())
data = data[:7]
alist.append([int(data[0]), data[1].strip(), data[2].strip(), int(data[3]), float(data[4]), float(data[5]), float(data[6])])
# only works with cubic boxes
box = map(float, boxre.match(file.next()).groups()[:3])
return (name, numatoms, alist, box)
def resnamecmp(x, y):
if ispoly(x) and ispoly(y):
return 0
if x[1] == "ION" and y[1] == "ION":
return cmp(x[2], y[2])
return cmp(x[1], y[1])
def writeGro(file, name, data, box):
data.sort(cmp=resnamecmp)
with open(file, "w") as f:
f.write(name)
f.write(str(len(data)) + "\n")
resnum, andx = (1, 1)
resname = data[0][1]
for elem in data:
if elem[1] != resname:
resname = elem[1]
resnum = 1
f.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (resnum, elem[1], elem[2], andx, elem[4], elem[5], elem[6]))
andx += 1
resnum += 1
f.write("%10.5f%10.5f%10.5f\n" % tuple(box))
def printTopology(data):
print("change your topology to:")
text = [("PGA", "PGA20\t{}".format(len([x for x in data if ispoly(x)])/60)),
("W", "W\t{}".format(len([x for x in data if isw(x)]))),
("WF", "WF\t{}".format(len([x for x in data if x[2] == "WF"]))),
("NA+", "NA+\t{}".format(len([x for x in data if x[2] == "NA+"]))),
("CL-", "CL-\t{}".format(len([x for x in data if x[2] == "CL-"])))]
text.sort(resnamecmp)
for line in text:
print(line[1])
def printPolTopology(data):
print("change your topology to:")
text = [("PGA", "PGA20\t{}".format(len([x for x in data if ispoly(x)])/60)),
("PW", "PW\t{}".format(len([x for x in data if isw(x)])/3)),
("NA+", "NA+\t{}".format(len([x for x in data if x[2] == "NA+"]))),
("CL-", "CL-\t{}".format(len([x for x in data if x[2] == "CL-"])))]
text.sort(resnamecmp)
for line in text:
print(line[1])
def main():
parser = argparse.ArgumentParser()
parser.add_argument("-p", "--polarizable", help="use polarizable water", action="store_true")
parser.add_argument("-d", "--double", help="double the universe (MDAnalysis 0.10 is required)", action="store_true")
parser.add_argument("-f", "--file", type=str, help="the input topology", default="zbox.gro")
parser.add_argument("-o", "--output", type=str, help="the output topology", default="sbox.gro")
parser.add_argument("-cm", "--concmid", type=float, help="the target concentration for the middle part of the box", default=0.017)
parser.add_argument("-co", "--concouter", type=float, help="the target concentration for the outer part of the box", default=0.51)
args = parser.parse_args()
name, a, data, box = readGro(args.file)
if args.polarizable:
changeTopology(data, box, change_ion_to_pw, change_pw_to_ion, args.concmid, args.concouter)
printPolTopology(data)
else:
changeTopology(data, box, change_ion_to_w, change_w_to_ion, args.concmid, args.concouter)
printTopology(data)
writeGro(args.output, name, data, box)
if __name__ == "__main__":
main()
Calculating the flow through membranes
This script can be used to calculate the flow of ions and water.
#!/usr/bin/python2
from __future__ import division
import MDAnalysis
import numpy
import sys
import math
def getZBounds(mem1, mem2):
bbox1 = mem1.bbox()
bbox2 = mem2.bbox()
return (bbox1[0][2], bbox1[1][2], bbox2[0][2], bbox2[1][2])
def partition(zcoord, bounds):
m1min, m1max, m2min, m2max = bounds
if zcoord < m1min:
return "outer"
if zcoord < m1max:
return "mem"
if zcoord < m2min:
return "mid"
if zcoord < m2max:
return "mem"
if m2max < zcoord:
return "outer"
def direction(frompart, topart):
dir = (0, "none")
if frompart == "mem" and topart == "mid":
dir = (1, "mem")
if frompart == "mem" and topart == "outer":
dir = (-1, "mem")
if frompart == "mid" and topart == "outer":
dir = (-1, "pass")
if frompart == "outer" and topart == "mid":
dir = (1, "pass")
return dir
def count(ind, orig, num, ts, bounds):
f, b, fm, bm = (0, 0, 0, 0)
if orig[0] == "init":
orig[:] = [partition(ts[ind[i]][2], bounds) for i in xrange(0, len(ind))]
return (f, b, fm, bm)
for offset in xrange(0, len(orig), num):
inds = [ind[i] for i in xrange(offset, offset + num)]
loc = [partition(ts[i][2], bounds) for i in inds]
if loc.count(loc[0]) != len(loc) or loc[0] == "mem":
continue
dire = [direction(orig[j], loc[i]) for i, j in enumerate(xrange(offset, offset + num))]
direp = sum([x[0] for x in dire if x[1] == "pass"])
if direp > 0:
f += 1
if direp < 0:
b += 1
direm = sum([x[0] for x in dire if x[1] == "mem"])
if direm > 0:
fm += 1
if direm < 0:
bm += 1
orig[offset:(offset + num)] = loc
return (f, b, fm, bm)
def norm(a, b, c):
return (a - b) / c
def main():
u = MDAnalysis.Universe(sys.argv[1], sys.argv[2])
prefix = sys.argv[3]
nai = u.selectAtoms("name NA+").indices()
nao = ["init"]
cli = u.selectAtoms("name CL-").indices()
clo = ["init"]
wi = u.selectAtoms("resname W or resname PW or resname WF").indices()
wo = ["init"]
middle = u.dimensions[2] / 2
mem = "not (resname W or resname WF or resname PW or name CL- or name NA+)"
mem1 = u.selectAtoms(mem + " and prop z < {}".format(middle))
mem2 = u.selectAtoms(mem + " and prop z > {}".format(middle))
nadc = [0, 0, 0, 0]
cldc = [0, 0, 0, 0]
wdc = [0, 0, 0, 0]
nal = len(nai)
cll = len(cli)
wl = len(wi)
with open(prefix + "_naclwflow.txt", "w") as f:
f.write("frame\tfna\tbna\tfnam\tbnam\tfcl\tbcl\tfclm\tbclm\tfw\tbw\tfwm\tbwm\t" +
"flna\tflnam\tflcl\tflclm\tflw\tflwm\t" +
"cflna\tcflnam\tcflcl\tcflclm\tcflw\tcflwm\n")
for ts in u.trajectory:
bounds = getZBounds(mem1, mem2)
nad = count(nai, nao, 1, ts, bounds)
cld = count(cli, clo, 1, ts, bounds)
if u.atoms[wi[0]].resname == "PW":
wd = count(wi, wo, 3, ts, bounds)
wlcor = wl/3
else:
wd = count(wi, wo, 1, ts, bounds)
wlcor = wl
flna = norm(nad[0], nad[1], nal)
flnam = norm(nad[2], nad[3], nal)
flcl = norm(cld[0], cld[1], cll)
flclm = norm(cld[2], cld[3], cll)
flw = norm(wd[0], wd[1], wlcor)
flwm = norm(wd[2], wd[3], wlcor)
nadc[:] = [x + y for x, y in zip(nadc, nad)]
cldc[:] = [x + y for x, y in zip(cldc, cld)]
wdc[:] = [x + y for x, y in zip(wdc, wd)]
cflna = norm(nadc[0], nadc[1], nal)
cflnam = norm(nadc[2], nadc[3], nal)
cflcl = norm(cldc[0], cldc[1], cll)
cflclm = norm(cldc[2], cldc[3], cll)
cflw = norm(wdc[0], wdc[1], wlcor)
cflwm = norm(wdc[2], wdc[3], wlcor)
f.write("{0}\t".format(ts.frame) +
"\t".join([str(x) for x in nad]) + "\t" +
"\t".join([str(x) for x in cld]) + "\t" +
"\t".join([str(x) for x in wd]) + "\t" +
"{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t".format(flna, flnam, flcl, flclm, flw, flwm) +
"{0}\t{1}\t{2}\t{3}\t{4}\t{5}\n".format(cflna, cflnam, cflcl, cflclm, cflw, cflwm))
if __name__ == "__main__":
main()
This script can be used to graph the flow of ions and water.
#!/usr/bin/gnuplot
set terminal pngcairo size 2666,776 enhanced font 'Verdana,20'
# define axis
# remove border on top and right and set color to gray
set style line 11 lc rgb '#808080' lt 1
set border 3 back ls 11
set tics nomirror
# define grid
set style line 12 lc rgb '#808080' lt 0 lw 1
set grid back ls 12
set ylabel "Proportion"
set xlabel "Time (ns)"
set key off
set output 'naclwflow.png'
set multiplot layout 1,3
set yrange [-0.05:0.5]
set autoscale x
set title "Net cumulative Na^{+} flow"
plot "cel_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fna") + column("fnam") - column("bna")) / 293) smooth cumulative with lines lt rgb "red" lw 2 title "Cellulose",\
"pga_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fna") + column("fnam") - column("bna")) / 1582) smooth cumulative with lines lt rgb "green" lw 2 title "Poly-γ-glutamic acid",\
"celp_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fna") + column("fnam") - column("bna")) / 3099) smooth cumulative with lines lt rgb "blue" lw 2 title "Cellulose phosphate"
set title "Net cumulative Cl^{-} flow"
plot "cel_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fcl") + column("fclm") - column("bcl")) / 293) smooth cumulative with lines lt rgb "red" lw 2 title "Cellulose",\
"pga_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fcl") + column("fclm") - column("bcl")) / 182) smooth cumulative with lines lt rgb "green" lw 2 title "Poly-γ-glutamic acid",\
"celp_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):((column("fcl") + column("fclm") - column("bcl")) / 299) smooth cumulative with lines lt rgb "blue" lw 2 title "Cellulose phosphate"
set autoscale y
set ylabel "Water count"
set title "Net cumulative water flow"
plot "cel_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):(column("fw") - column("bw")) smooth cumulative with lines lt rgb "red" lw 2 title "Cellulose",\
"pga_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):(column("fw") - column("bw")) smooth cumulative with lines lt rgb "green" lw 2 title "Poly-γ-glutamic acid",\
"celp_naclwflow.txt" u (20 * (column("frame") - 1) / 1000):(column("fw") - column("bw")) smooth cumulative with lines lt rgb "blue" lw 2 title "Cellulose phosphate"
unset multiplot