Difference between revisions of "Template:Team:Groningen/CONTENT/MODELING/Supporting material"

(fix links)
(add chmemions and naclwflow)
Line 32: Line 32:
 
<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 class="item">Changing the concentration and preparing the membranes (<a href="#scchmem" class="link">link</a>)</div>
 +
<div class="item">Calculating the flow of ions through the membranes (<a href="#scnaclw" class="link">link</a>)</div>
 
</div>
 
</div>
  
Line 493: Line 495:
 
if __name__ == "__main__":
 
if __name__ == "__main__":
 
     main()
 
     main()
 +
</code></div></div>
 +
 +
 +
<div class="subtitle" id="scchmem">Changing ion concentrations in membranes</div>
 +
<div class="object code">
 +
<div class="caption">This script can be used to prepare a system with two membranes for simualtions.</div>
 +
 +
<div class="raw"><code class="python">
 +
#!/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()
 +
</code></div></div>
 +
 +
<div class="subtitle" id="scnaclw">Calculating the flow through membranes</div>
 +
<div class="object code">
 +
<div class="caption">This script can be used to calculate the flow of ions and water.</div>
 +
 +
<div class="raw"><code class="python">
 +
#!/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()
 +
</code></div></div>
 +
 +
<div class="object code">
 +
<div class="caption">This script can be used to graph the flow of ions and water.</div>
 +
 +
<div class="raw"><code>
 +
#!/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
 
</code></div></div>
 
</code></div></div>
  
 
</div>
 
</div>

Revision as of 00:20, 19 September 2015

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.

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>)
Changing the concentration and preparing the membranes (<a href="#scchmem" class="link">link</a>)
Calculating the flow of ions through the membranes (<a href="#scnaclw" class="link">link</a>)
PGA topology creation
This script can be used to generate a PGA topology file.
  1. /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.
  1. !/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.
  1. !/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.
  1. !/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.
  1. !/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.
  1. !/usr/bin/gnuplot

set terminal pngcairo size 2666,776 enhanced font 'Verdana,20'

  1. define axis
  2. 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

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