Skip to content
Snippets Groups Projects
Commit 8d84d022 authored by cmaffeo2's avatar cmaffeo2
Browse files

Add tests to repository

parent b50a4efc
No related branches found
No related tags found
No related merge requests found
Showing
with 1592 additions and 0 deletions
BrownDyn.out*
BrownDyn.log
rho.dat
__pycache__
*.pdf
\ No newline at end of file
timestep 2e-6
steps 1000000
temperature 295
outputPeriod 1000
outputEnergyPeriod 1000
decompPeriod 100
outputFormat dcd
cutoff 10.0
# Potassium
particle POT
num 8000
gridFile pmf.dx
gridFileScale 1.0
diffusion 5
tabulatedPotential 1
tabulatedFile 0@0@vdw.dat
71.01% 47.987ms 11 4.3624ms 4.3024ms 4.4236ms createPairlists_debug(Vector3*, int, int, BaseGrid const *, CellDecomposition const *, int, int*, int2*, int, int const *, int*, Exclude const *, int2 const *, int, float)
60.08% 28.986ms 11 2.6351ms 2.6195ms 2.6664ms void createPairlists<int=64, int=64, int=8>(Vector3*, int, int, BaseGrid const *, CellDecomposition const *, int, int*, int2*, int, int const *, int*, Exclude const *, int2 const *, int, float)
# a(v) = 1.76375*[1 - (1-100/1.76375)*exp(-v/2.93) ]
object 1 class gridpositions counts 1 1 1
origin -20 -20 -20
delta 43.083 0 0
delta 0 43.083 0
delta 0 0 43.083
object 2 class gridconnections counts 1 1 1
object 3 class array type double rank 0 items 1 data follows
0
object "3D grid" class field
#! /bin/bash
mkdir -p output
../../src/arbd $@ BrownDyn.bd output/BrownDyn
0.000000 5328.199575
0.100000 5088.358871
0.200000 4848.518167
0.300000 4608.677462
0.400000 4368.836758
0.500000 4128.996053
0.600000 3889.155349
0.700000 3649.314644
0.800000 3409.473940
0.900000 3169.633235
1.000000 2929.792531
1.100000 2689.951826
1.200000 2450.111122
1.300000 2210.270417
1.400000 1970.429713
1.500000 1730.589008
1.600000 1490.748304
1.700000 1250.907600
1.800000 1011.066895
1.900000 771.226191
2.000000 531.385486
2.100000 291.544782
2.200000 163.668090
2.300000 93.681755
2.400000 54.481429
2.500000 32.073032
2.600000 19.034741
2.700000 11.333138
2.800000 6.727069
2.900000 3.945871
3.000000 2.255816
3.100000 1.226149
3.200000 0.600191
3.300000 0.222952
3.400000 -0.000256
3.500000 -0.127862
3.600000 -0.196223
3.700000 -0.228109
3.800000 -0.237841
3.900000 -0.234446
4.000000 -0.223608
4.100000 -0.208887
4.200000 -0.192491
4.300000 -0.175762
4.400000 -0.159494
4.500000 -0.144129
4.600000 -0.129888
4.700000 -0.116855
4.800000 -0.105032
4.900000 -0.094370
5.000000 -0.084796
5.100000 -0.076221
5.200000 -0.068557
5.300000 -0.061712
5.400000 -0.055603
5.500000 -0.050151
5.600000 -0.045285
5.700000 -0.040939
5.800000 -0.037054
5.900000 -0.033580
6.000000 -0.030470
6.100000 -0.027682
6.200000 -0.025181
6.300000 -0.022934
6.400000 -0.020914
6.500000 -0.019095
6.600000 -0.017456
6.700000 -0.015976
6.800000 -0.014639
6.900000 -0.013429
7.000000 -0.012333
7.100000 -0.011339
7.200000 -0.010436
7.300000 -0.009616
7.400000 -0.008869
7.500000 -0.008189
7.600000 -0.007569
7.700000 -0.007002
7.800000 -0.006484
7.900000 -0.006010
8.000000 -0.005575
8.100000 -0.005177
8.200000 -0.004812
8.300000 -0.004476
8.400000 -0.004167
8.500000 -0.003882
8.600000 -0.003620
8.700000 -0.003378
8.800000 -0.003155
8.900000 -0.002949
9.000000 -0.002758
9.100000 -0.002582
9.200000 -0.002418
9.300000 -0.002267
9.400000 -0.002126
9.500000 -0.001996
9.600000 -0.001875
9.700000 -0.001762
9.800000 -0.001657
9.900000 -0.001559
10.000000 -0.001468
10.100000 -0.001383
10.200000 -0.001304
10.300000 -0.001230
10.400000 -0.001160
10.500000 -0.001096
10.600000 -0.001035
10.700000 -0.000979
10.800000 -0.000926
10.900000 -0.000876
11.000000 -0.000829
11.100000 -0.000785
11.200000 -0.000744
11.300000 -0.000706
11.400000 -0.000669
11.500000 -0.000635
11.600000 -0.000603
11.700000 -0.000573
11.800000 -0.000544
11.900000 -0.000517
12.000000 -0.000492
12.100000 -0.000468
12.200000 -0.000446
12.300000 -0.000424
12.400000 -0.000404
12.500000 -0.000385
12.600000 -0.000367
12.700000 -0.000350
12.800000 -0.000334
12.900000 -0.000319
13.000000 -0.000304
13.100000 -0.000291
13.200000 -0.000278
13.300000 -0.000265
13.400000 -0.000254
13.500000 -0.000243
13.600000 -0.000232
13.700000 -0.000222
13.800000 -0.000213
13.900000 -0.000204
14.000000 -0.000195
14.100000 -0.000187
14.200000 -0.000179
14.300000 -0.000172
14.400000 -0.000165
14.500000 -0.000158
14.600000 -0.000152
14.700000 -0.000146
14.800000 -0.000140
14.900000 -0.000134
15.000000 -0.000129
15.100000 -0.000124
15.200000 -0.000119
15.300000 -0.000115
15.400000 -0.000110
15.500000 -0.000106
15.600000 -0.000102
15.700000 -0.000098
15.800000 -0.000094
15.900000 -0.000091
16.000000 -0.000088
16.100000 -0.000084
16.200000 -0.000081
16.300000 -0.000078
16.400000 -0.000075
16.500000 -0.000073
16.600000 -0.000070
16.700000 -0.000068
16.800000 -0.000065
16.900000 -0.000063
17.000000 -0.000061
17.100000 -0.000059
17.200000 -0.000057
17.300000 -0.000055
17.400000 -0.000053
17.500000 -0.000051
17.600000 -0.000049
17.700000 -0.000048
17.800000 -0.000046
17.900000 -0.000045
18.000000 -0.000043
18.100000 -0.000042
18.200000 -0.000040
18.300000 -0.000039
18.400000 -0.000038
18.500000 -0.000037
18.600000 -0.000035
18.700000 -0.000034
18.800000 -0.000033
18.900000 -0.000032
19.000000 -0.000031
19.100000 -0.000030
19.200000 -0.000029
19.300000 -0.000028
19.400000 -0.000027
19.500000 -0.000027
19.600000 -0.000026
19.700000 -0.000025
19.800000 -0.000024
19.900000 -0.000024
20.000000 -0.000023
20.100000 -0.000022
20.200000 -0.000022
20.300000 -0.000021
20.400000 -0.000020
20.500000 -0.000020
20.600000 -0.000019
20.700000 -0.000019
20.800000 -0.000018
20.900000 -0.000018
21.000000 -0.000017
21.100000 -0.000017
21.200000 -0.000016
21.300000 -0.000016
21.400000 -0.000015
21.500000 -0.000015
21.600000 -0.000014
21.700000 -0.000014
21.800000 -0.000014
21.900000 -0.000013
22.000000 -0.000013
22.100000 -0.000013
22.200000 -0.000012
22.300000 -0.000012
22.400000 -0.000012
22.500000 -0.000011
22.600000 -0.000011
22.700000 -0.000011
22.800000 -0.000010
22.900000 -0.000010
23.000000 -0.000010
23.100000 -0.000010
23.200000 -0.000009
23.300000 -0.000009
23.400000 -0.000009
23.500000 -0.000009
23.600000 -0.000008
23.700000 -0.000008
23.800000 -0.000008
23.900000 -0.000008
24.000000 -0.000008
24.100000 -0.000007
24.200000 -0.000007
24.300000 -0.000007
24.400000 -0.000007
24.500000 -0.000007
24.600000 -0.000007
24.700000 -0.000006
24.800000 -0.000006
24.900000 -0.000006
25.000000 -0.000006
25.100000 -0.000006
25.200000 -0.000006
25.300000 -0.000006
25.400000 -0.000005
25.500000 -0.000005
25.600000 -0.000005
25.700000 -0.000005
25.800000 -0.000005
25.900000 -0.000005
26.000000 -0.000005
26.100000 -0.000005
26.200000 -0.000004
26.300000 -0.000004
26.400000 -0.000004
26.500000 -0.000004
26.600000 -0.000004
26.700000 -0.000004
26.800000 -0.000004
26.900000 -0.000004
27.000000 -0.000004
27.100000 -0.000004
27.200000 -0.000004
27.300000 -0.000003
27.400000 -0.000003
27.500000 -0.000003
27.600000 -0.000003
27.700000 -0.000003
27.800000 -0.000003
27.900000 -0.000003
28.000000 -0.000003
28.100000 -0.000003
28.200000 -0.000003
28.300000 -0.000003
28.400000 -0.000003
28.500000 -0.000003
28.600000 -0.000003
28.700000 -0.000003
28.800000 -0.000002
28.900000 -0.000002
29.000000 -0.000002
29.100000 -0.000002
29.200000 -0.000002
29.300000 -0.000002
29.400000 -0.000002
29.500000 -0.000002
29.600000 -0.000002
29.700000 -0.000002
29.800000 -0.000002
29.900000 -0.000002
30.000000 -0.000002
30.100000 -0.000002
30.200000 -0.000002
30.300000 -0.000002
30.400000 -0.000002
30.500000 -0.000002
30.600000 -0.000002
30.700000 -0.000002
30.800000 -0.000002
30.900000 -0.000002
31.000000 -0.000002
31.100000 -0.000002
31.200000 -0.000001
31.300000 -0.000001
31.400000 -0.000001
31.500000 -0.000001
31.600000 -0.000001
31.700000 -0.000001
31.800000 -0.000001
31.900000 -0.000001
32.000000 -0.000001
32.100000 -0.000001
32.200000 -0.000001
32.300000 -0.000001
32.400000 -0.000001
32.500000 -0.000001
32.600000 -0.000001
32.700000 -0.000001
32.800000 -0.000001
32.900000 -0.000001
33.000000 -0.000001
33.100000 -0.000001
33.200000 -0.000001
33.300000 -0.000001
33.400000 -0.000001
33.500000 -0.000001
33.600000 -0.000001
33.700000 -0.000001
33.800000 -0.000001
33.900000 -0.000001
34.000000 -0.000001
34.100000 -0.000001
34.200000 -0.000001
34.300000 -0.000001
34.400000 -0.000001
34.500000 -0.000001
34.600000 -0.000001
34.700000 -0.000001
34.800000 -0.000001
34.900000 -0.000001
35.000000 -0.000001
35.100000 -0.000001
35.200000 -0.000001
35.300000 -0.000001
35.400000 -0.000001
35.500000 -0.000001
35.600000 -0.000001
35.700000 -0.000001
35.800000 -0.000001
35.900000 -0.000001
36.000000 -0.000001
36.100000 -0.000001
36.200000 -0.000001
36.300000 -0.000001
36.400000 -0.000001
36.500000 -0.000001
36.600000 -0.000001
36.700000 -0.000001
36.800000 -0.000000
36.900000 -0.000000
37.000000 -0.000000
37.100000 -0.000000
37.200000 -0.000000
37.300000 -0.000000
37.400000 -0.000000
37.500000 -0.000000
37.600000 -0.000000
37.700000 -0.000000
37.800000 -0.000000
37.900000 -0.000000
38.000000 -0.000000
38.100000 -0.000000
38.200000 -0.000000
38.300000 -0.000000
38.400000 -0.000000
38.500000 -0.000000
38.600000 -0.000000
38.700000 -0.000000
38.800000 -0.000000
38.900000 -0.000000
39.000000 -0.000000
39.100000 -0.000000
39.200000 -0.000000
39.300000 -0.000000
39.400000 -0.000000
39.500000 -0.000000
39.600000 -0.000000
39.700000 -0.000000
39.800000 -0.000000
39.900000 -0.000000
40.000000 -0.000000
40.100000 -0.000000
40.200000 -0.000000
40.300000 -0.000000
40.400000 -0.000000
40.500000 -0.000000
40.600000 -0.000000
40.700000 -0.000000
40.800000 -0.000000
40.900000 -0.000000
41.000000 -0.000000
41.100000 -0.000000
41.200000 -0.000000
41.300000 -0.000000
41.400000 -0.000000
41.500000 -0.000000
41.600000 -0.000000
41.700000 -0.000000
41.800000 -0.000000
41.900000 -0.000000
42.000000 -0.000000
42.100000 -0.000000
42.200000 -0.000000
42.300000 -0.000000
42.400000 -0.000000
42.500000 -0.000000
42.600000 -0.000000
42.700000 -0.000000
42.800000 -0.000000
42.900000 -0.000000
43.000000 -0.000000
43.100000 -0.000000
43.200000 -0.000000
43.300000 -0.000000
43.400000 -0.000000
43.500000 -0.000000
43.600000 -0.000000
43.700000 -0.000000
43.800000 -0.000000
43.900000 -0.000000
44.000000 -0.000000
44.100000 -0.000000
44.200000 -0.000000
44.300000 -0.000000
44.400000 -0.000000
44.500000 -0.000000
44.600000 -0.000000
44.700000 -0.000000
44.800000 -0.000000
44.900000 -0.000000
45.000000 -0.000000
45.100000 -0.000000
45.200000 -0.000000
45.300000 -0.000000
45.400000 -0.000000
45.500000 -0.000000
45.600000 -0.000000
45.700000 -0.000000
45.800000 -0.000000
45.900000 -0.000000
46.000000 -0.000000
46.100000 -0.000000
46.200000 -0.000000
46.300000 -0.000000
46.400000 -0.000000
46.500000 -0.000000
46.600000 -0.000000
46.700000 -0.000000
46.800000 -0.000000
46.900000 -0.000000
47.000000 -0.000000
47.100000 -0.000000
47.200000 -0.000000
47.300000 -0.000000
47.400000 -0.000000
47.500000 -0.000000
47.600000 -0.000000
47.700000 -0.000000
47.800000 -0.000000
47.900000 -0.000000
48.000000 -0.000000
48.100000 -0.000000
48.200000 -0.000000
48.300000 -0.000000
48.400000 -0.000000
48.500000 -0.000000
48.600000 -0.000000
48.700000 -0.000000
48.800000 -0.000000
48.900000 -0.000000
49.000000 -0.000000
49.100000 -0.000000
49.200000 -0.000000
49.300000 -0.000000
49.400000 -0.000000
49.500000 -0.000000
49.600000 -0.000000
49.700000 -0.000000
49.800000 -0.000000
49.900000 -0.000000
50.000000 0.000000
#! /usr/bin/env python
import numpy as np
eps=0.238000
rmin = 2*1.908100
sig = rmin/(2**(1/6))
r = np.linspace(0,50,501)
r[0] = 0.1
u = 4*eps*( (sig/r)**12 - (sig/r)**6 )
r[0] = 0
ids = np.where( u > 200 )
idMax = ids[0][-1]
du = u[idMax] - u[idMax-1]
u[ids] = u[idMax] -du * (np.arange(idMax,-1,-1) )
u = u-u[-1]
ch = open('dummy-pot.dat','w')
for data in zip(r,u):
ch.write("%f %f\n" % data)
timestep 2e-6
steps 1000000
temperature 295
outputPeriod 1000
outputEnergyPeriod 1000
decompPeriod 100
outputFormat dcd
cutoff 10.0
# Potassium
particle POT
num 64000
gridFile pmf.dx
diffusion 5
tabulatedPotential 1
tabulatedFile 0@0@vdw.dat
all:
$(MAKE) rho.dat
$(MAKE) $(patsubst %.plt.py,%.pdf, $(wildcard *.plt.py))
plot.pdf: plot.plt.py rho.dat plotTools.py
./$< $@
rho.dat: get-data.tcl ../BrownDyn.out.0.dcd
vmd -dispdev text < $<
.SECONDARY:
set ID [mol new ../BrownDyn.out.0.dcd first 50 last 100 waitfor all]
set sel [atomselect $ID all]
set out rho.dat
set rmax 10
set dr 0.1
## center on first atom
set size [expr 57*0.7017544]
set sel [atomselect $ID "index 1"]
set all [atomselect $ID "all"]
set dim [expr 1.5117018*57]
set dim "$dim $dim $dim"
set last [expr [molinfo $ID get numframes]-1]
for {set f 0} {$f <= $last} {incr f} {
animate goto $f
molinfo $ID set {a b c} $dim
molinfo $ID set {alpha beta gamma} {90 90 90}
}
# set result [measure gofr $all $all delta 0.2 usepbc 1 first 0 last $last]
set result [measure gofr $all $all delta 0.1 usepbc 1 first 0 last $last step 5]
lassign $result rs gofrs int hist frames
set ch [open $out w]
foreach r $rs g $gofrs {
puts $ch "$r $g"
}
#! /usr/bin/env python3
# /home/cmaffeo2/bin/python3.sh
# /software/python3-3.4.1/bin/python3
import matplotlib as mpl
mpl.use('agg')
import numpy as np
from scipy import exp
from scipy import sqrt
import matplotlib.pyplot as plt
# import scipy.optimize as opt
from glob import glob
# from natsort import natsorted
import re
from plotTools import *
import sys
def info(*obj):
print('INFO: ',*obj , file=sys.stderr)
def makeTitle(ax,text):
# ax.set_title(text)
ax.annotate(text, xy=(0.05,0.92),
fontsize=9, xycoords='axes fraction',
horizontalalignment='left', verticalalignment='top')
fig,axes = plt.subplots(1)
ax1 = axes
# makeTitle(ax1,'Pair interaction (rev. a158e71)')
ax1.set_ylabel('likelihood')
ax1.get_yaxis().set_label_coords(-0.085,0.5)
ax1.set_xlabel("Pair distribution g(r)")
xmin, xmax = [0,10]
setRange(ax1,[xmin,xmax, 0,5])
def loadData(fname):
d0 = np.loadtxt(fname)
x0 = d0[:,0]
y0 = d0[:,1]
# ids = np.where( (x0 > xmin) & (x0 < xmax) )
# x0 = x0[ids]
# y0 = y0[ids]
# y0 = y0 / np.sum(y0)
return x0,y0
def estimate():
eps = 2*10
sig = 3
x = np.linspace(xmin,xmax,100)
u = 4*eps*( (sig/x)**12 - (sig/x)**6 )
y = np.exp( - u ) * x**2
y = y / np.sum(y)
return x,y
def estimate():
d0 = np.loadtxt('../dummy-pot.dat')
x = d0[10:,0]
u = d0[10:,1]
ids = np.where( (x > xmin) & (x < xmax) )
x = x[ids]
u = u[ids]
y = np.exp( - u ) * x**2
y = y / np.sum(y)
return x,y
def plot(x,y,l,c):
return ax1.plot(x,y, ls='-', marker='', color=c, mec=c)
# return ax1.bar(x,y, label=l, lw=0, width=np.mean(x[1:]-x[:-1]), alpha=0.5, color=c)
x,y = loadData('rho.dat')
ex,ey = estimate()
plot(x,y,'simulated', c0)
# plot(ex,ey,'expected',c1)
plt.legend()
# plt.show()
plt.tight_layout()
# plt.savefig(sys.stdout, format='svg')
plt.savefig(sys.argv[1])
import numpy as np
import matplotlib as mpl
import matplotlib.ticker as plticker
import sys
def hsv2rgb(h,s,v):
# color=hsv_to_rgb( np.reshape( np.array([1./3, 0.5, 1]), [1,1,3] ) )[0,0,:],
hsv = np.reshape( [h,s,v], [1,1,3] )
return mpl.colors.hsv_to_rgb( hsv )[0,0,:]
# c0 = hsv2rgb(1/3, 0.5, 0.9)
# c1a = hsv2rgb(1/2, 0.25, 0.9)
# c1b = hsv2rgb(1/2, 0.5, 0.7)
# c1c = hsv2rgb(1/2, 0.75, 0.5)
# c2 = hsv2rgb(2/3, 0.5, 0.9)
# c3 = hsv2rgb(1/6, 0.5, 0.9)
# c3 = hsv2rgb(1/2, 0.5, 0.9)
def hsv2rgb256(h,s,v):
return hsv2rgb(h/360,s/255,v/255)
# http://mathematica.stackexchange.com/questions/20851/mathematica-color-schemes-for-the-colorblind
# c0 = np.array([51, 34, 136])/255
# c1a = np.array([136, 204, 238])/255
# c1b = np.array([68, 170, 153])/255
# c1c = np.array([17, 119, 51])/255
# c2 = np.array([153, 153, 51])/255
# # c3 = np.array([136, 34, 85])/255
# c3 = np.array([204, 121, 167])/255
# c0 = np.array([153, 34, 136])/255
# c1a = np.array([51, 102, 170])/255
# c1b = np.array([17, 170, 153])/255
# c1c = np.array([102, 170, 85])/255
# c2 = np.array([238, 51, 51])/255
# c3 = np.array([238, 119, 34])/255
# c0 = np.array([136, 46, 114])/255
# c1a = np.array([ 25, 101, 176])/255
# c1b = np.array([ 82, 137, 199])/255
# c1c = np.array([123, 175, 222])/255
# # c2 = np.array([144, 201, 135])/255
# # c2 = np.array([246, 193, 65])/255
# c2 = np.array([241, 147, 45])/255
# c3 = np.array([ 78, 178, 101])/255
# c0 = hsv2rgb256(314, 169, 136)
c0 = hsv2rgb256(280, 226, 94)
# c1a = hsv2rgb256(209, 219, 176)
# c1b = hsv2rgb256(211, 150, 199)
# c1c = hsv2rgb256(208, 114, 222)
c1a = hsv2rgb256(209, 219, 120)
c1b = hsv2rgb256(211, 150, 180)
c1c = hsv2rgb256(208, 114, 230)
c2 = hsv2rgb256( 31, 208, 241)
c3 = hsv2rgb256(133, 143, 178)
c1 = c1b
m0 = 's'
m1 = '<'
m2 = 'o'
m3 = '>'
# ----------------------------------------------------------------- #
# Make colormap based on Paul Tol's best visibility gradients. See #
# <http://www.sron.nl/~pault/> for more info on these colors. Also #
# see <http://matplotlib.sourceforge.net/api/colors_api.html> and #
# <http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps> on some #
# matplotlib examples #
# ----------------------------------------------------------------- #
# Deviation around zero colormap (blue--red)
cols = []
for x in np.linspace(0,1, 256):
rcol = 0.237 - 2.13*x + 26.92*x**2 - 65.5*x**3 + 63.5*x**4 - 22.36*x**5
gcol = ((0.572 + 1.524*x - 1.811*x**2)/(1 - 0.291*x + 0.1574*x**2))**2
bcol = 1/(1.579 - 4.03*x + 12.92*x**2 - 31.4*x**3 + 48.6*x**4 - 23.36*x**5)
cols.append((rcol, gcol, bcol))
cm_plusmin = mpl.colors.LinearSegmentedColormap.from_list("PaulT_plusmin", cols)
# Linear colormap (white--red)
from scipy.special import erf
cols = []
for x in np.linspace(0,1, 256):
rcol = (1 - 0.392*(1 + erf((x - 0.869)/ 0.255)))
gcol = (1.021 - 0.456*(1 + erf((x - 0.527)/ 0.376)))
bcol = (1 - 0.493*(1 + erf((x - 0.272)/ 0.309)))
cols.append((rcol, gcol, bcol))
cm_linear = mpl.colors.LinearSegmentedColormap.from_list("PaulT_linear", cols)
# Linear colormap (rainbow)
cols = [(0,0,0)]
for x in np.linspace(0,1, 254):
rcol = (0.472-0.567*x+4.05*x**2)/(1.+8.72*x-19.17*x**2+14.1*x**3)
gcol = 0.108932-1.22635*x+27.284*x**2-98.577*x**3+163.3*x**4-131.395*x**5+40.634*x**6
bcol = 1./(1.97+3.54*x-68.5*x**2+243*x**3-297*x**4+125*x**5)
cols.append((rcol, gcol, bcol))
cols.append((1,1,1))
cm_rainbow = mpl.colors.LinearSegmentedColormap.from_list("PaulT_rainbow", cols)
cm_default = cm_plusmin
def info(*obj):
print('INFO: ',*obj , file=sys.stderr)
def blkavg(a,blklength):
"""
a is (n*N+extra)x2 array
blklength (units of x-axis)
returns Nx2
"""
t = a[:,0]
dt = t[1]-t[0]
y = a[:,1]
# n is number of elements in a block, N is the number of blocks
n = int(np.floor(blklength/dt))
N = int(np.floor(len(t)/n))
tb = np.reshape(t[:n*N],(N,n)).mean(axis=1)
yb = np.reshape(y[:n*N],(N,n)).mean(axis=1)
r = np.vstack((tb,yb)).transpose()
# v = np.gradient(yb,dt*n)
# s = v.std()/np.sqrt(N)
# return r,v,s
return r
# def blockAvg(data, pointsPerBlock=0, nblocks=100):
# ## http://stackoverflow.com/questions/14229029/block-mean-of-numpy-2d-array
# if pointsPerBlock > 0:
# nblocks = len(data)/pointsPerBlock
# l = data.shape[0]
# pointsPerBlock = int(l/nblocks)
# data = data[
# data.reshape( (l*nblocks)
# info(nblocks)
# info(data.shape)
# blocks = data.reshape((data.shape[0], -1, nblocks))
# info(blocks)
# blockaverage = np.mean(blocks, axis=1)
# # blockstd = np.std(blocks])
# # return blockaverage, blockstd
# return blockaverage
def getTickDivisions(xmin,xmax):
dx = np.abs(xmax-xmin)
minticksize = dx / 2
magnitude = np.log(minticksize)/np.log(10)
tmp = magnitude
magnitude = 10**np.floor(magnitude)
residual = minticksize / magnitude
if np.abs(np.round(residual*2) - residual*2) < 1e-6 and magnitude > 0.01:
tick = residual * magnitude
elif residual > 5:
tick = 10*magnitude
elif residual > 2:
tick = 5*magnitude
elif residual > 1:
tick = 2*magnitude
else:
tick = magnitude
# info('range:',dx)
# info('minticksize:',minticksize)
# info('tmp:',tmp)
# info('magnitude:',magnitude)
# info('residual:',residual)
# info('tick:',tick)
# info(dx,magnitude,residual,tick)
return tick
# for f in np.linspace(0,1,11)[1:]:
# getTickDivisions(0,f)
def setRange(axes, minmax):
xmin,xmax,ymin,ymax = minmax
axes.axis(minmax)
dx = getTickDivisions(xmin,xmax)
dy = getTickDivisions(ymin,ymax)
setTicks(axes,dx,dy)
# def setTicks(axes, dx, dy, numXminor=1, numYminor=1, xoffset=0, yoffset=0):
def setTicks(axes, dx, dy, numXminor=1, numYminor=1):
axes.xaxis.set_major_locator( plticker.MultipleLocator( base=dx ))
axes.xaxis.set_minor_locator( plticker.MultipleLocator( base=dx/(numXminor+1) ))
axes.yaxis.set_major_locator( plticker.MultipleLocator( base=dy ))
axes.yaxis.set_minor_locator( plticker.MultipleLocator( base=dy/(numXminor+1) ))
62.57% 1.63078s 101 16.146ms 15.913ms 18.956ms void createPairlists<int=64, int=64, int=8>(Vector3*, int, int, BaseGrid const *, CellDecomposition const *, int, int*, int2*, int, int const *, int*, Exclude const *, int2 const *, int, float)
75.75% 3.10843s 101 30.777ms 30.282ms 36.200ms createPairlists_debug(Vector3*, int, int, BaseGrid const *, CellDecomposition const *, int, int*, int2*, int, int const *, int*, Exclude const *, int2 const *, int, float)
# a(v) = 1.76375*[1 - (1-100/1.76375)*exp(-v/2.93) ]
object 1 class gridpositions counts 1 1 1
origin -20 -20 -20
delta 86.167 0 0
delta 0 86.167 0
delta 0 0 86.167
object 2 class gridconnections counts 1 1 1
object 3 class array type double rank 0 items 1 data follows
0
object "3D grid" class field
#! /bin/bash
mkdir output
nvprof ../../src/arbd_lin_lin $@ BrownDyn.bd output/BrownDyn
0.000000 5328.199575
0.100000 5088.358871
0.200000 4848.518167
0.300000 4608.677462
0.400000 4368.836758
0.500000 4128.996053
0.600000 3889.155349
0.700000 3649.314644
0.800000 3409.473940
0.900000 3169.633235
1.000000 2929.792531
1.100000 2689.951826
1.200000 2450.111122
1.300000 2210.270417
1.400000 1970.429713
1.500000 1730.589008
1.600000 1490.748304
1.700000 1250.907600
1.800000 1011.066895
1.900000 771.226191
2.000000 531.385486
2.100000 291.544782
2.200000 163.668090
2.300000 93.681755
2.400000 54.481429
2.500000 32.073032
2.600000 19.034741
2.700000 11.333138
2.800000 6.727069
2.900000 3.945871
3.000000 2.255816
3.100000 1.226149
3.200000 0.600191
3.300000 0.222952
3.400000 -0.000256
3.500000 -0.127862
3.600000 -0.196223
3.700000 -0.228109
3.800000 -0.237841
3.900000 -0.234446
4.000000 -0.223608
4.100000 -0.208887
4.200000 -0.192491
4.300000 -0.175762
4.400000 -0.159494
4.500000 -0.144129
4.600000 -0.129888
4.700000 -0.116855
4.800000 -0.105032
4.900000 -0.094370
5.000000 -0.084796
5.100000 -0.076221
5.200000 -0.068557
5.300000 -0.061712
5.400000 -0.055603
5.500000 -0.050151
5.600000 -0.045285
5.700000 -0.040939
5.800000 -0.037054
5.900000 -0.033580
6.000000 -0.030470
6.100000 -0.027682
6.200000 -0.025181
6.300000 -0.022934
6.400000 -0.020914
6.500000 -0.019095
6.600000 -0.017456
6.700000 -0.015976
6.800000 -0.014639
6.900000 -0.013429
7.000000 -0.012333
7.100000 -0.011339
7.200000 -0.010436
7.300000 -0.009616
7.400000 -0.008869
7.500000 -0.008189
7.600000 -0.007569
7.700000 -0.007002
7.800000 -0.006484
7.900000 -0.006010
8.000000 -0.005575
8.100000 -0.005177
8.200000 -0.004812
8.300000 -0.004476
8.400000 -0.004167
8.500000 -0.003882
8.600000 -0.003620
8.700000 -0.003378
8.800000 -0.003155
8.900000 -0.002949
9.000000 -0.002758
9.100000 -0.002582
9.200000 -0.002418
9.300000 -0.002267
9.400000 -0.002126
9.500000 -0.001996
9.600000 -0.001875
9.700000 -0.001762
9.800000 -0.001657
9.900000 -0.001559
10.000000 -0.001468
10.100000 -0.001383
10.200000 -0.001304
10.300000 -0.001230
10.400000 -0.001160
10.500000 -0.001096
10.600000 -0.001035
10.700000 -0.000979
10.800000 -0.000926
10.900000 -0.000876
11.000000 -0.000829
11.100000 -0.000785
11.200000 -0.000744
11.300000 -0.000706
11.400000 -0.000669
11.500000 -0.000635
11.600000 -0.000603
11.700000 -0.000573
11.800000 -0.000544
11.900000 -0.000517
12.000000 -0.000492
12.100000 -0.000468
12.200000 -0.000446
12.300000 -0.000424
12.400000 -0.000404
12.500000 -0.000385
12.600000 -0.000367
12.700000 -0.000350
12.800000 -0.000334
12.900000 -0.000319
13.000000 -0.000304
13.100000 -0.000291
13.200000 -0.000278
13.300000 -0.000265
13.400000 -0.000254
13.500000 -0.000243
13.600000 -0.000232
13.700000 -0.000222
13.800000 -0.000213
13.900000 -0.000204
14.000000 -0.000195
14.100000 -0.000187
14.200000 -0.000179
14.300000 -0.000172
14.400000 -0.000165
14.500000 -0.000158
14.600000 -0.000152
14.700000 -0.000146
14.800000 -0.000140
14.900000 -0.000134
15.000000 -0.000129
15.100000 -0.000124
15.200000 -0.000119
15.300000 -0.000115
15.400000 -0.000110
15.500000 -0.000106
15.600000 -0.000102
15.700000 -0.000098
15.800000 -0.000094
15.900000 -0.000091
16.000000 -0.000088
16.100000 -0.000084
16.200000 -0.000081
16.300000 -0.000078
16.400000 -0.000075
16.500000 -0.000073
16.600000 -0.000070
16.700000 -0.000068
16.800000 -0.000065
16.900000 -0.000063
17.000000 -0.000061
17.100000 -0.000059
17.200000 -0.000057
17.300000 -0.000055
17.400000 -0.000053
17.500000 -0.000051
17.600000 -0.000049
17.700000 -0.000048
17.800000 -0.000046
17.900000 -0.000045
18.000000 -0.000043
18.100000 -0.000042
18.200000 -0.000040
18.300000 -0.000039
18.400000 -0.000038
18.500000 -0.000037
18.600000 -0.000035
18.700000 -0.000034
18.800000 -0.000033
18.900000 -0.000032
19.000000 -0.000031
19.100000 -0.000030
19.200000 -0.000029
19.300000 -0.000028
19.400000 -0.000027
19.500000 -0.000027
19.600000 -0.000026
19.700000 -0.000025
19.800000 -0.000024
19.900000 -0.000024
20.000000 -0.000023
20.100000 -0.000022
20.200000 -0.000022
20.300000 -0.000021
20.400000 -0.000020
20.500000 -0.000020
20.600000 -0.000019
20.700000 -0.000019
20.800000 -0.000018
20.900000 -0.000018
21.000000 -0.000017
21.100000 -0.000017
21.200000 -0.000016
21.300000 -0.000016
21.400000 -0.000015
21.500000 -0.000015
21.600000 -0.000014
21.700000 -0.000014
21.800000 -0.000014
21.900000 -0.000013
22.000000 -0.000013
22.100000 -0.000013
22.200000 -0.000012
22.300000 -0.000012
22.400000 -0.000012
22.500000 -0.000011
22.600000 -0.000011
22.700000 -0.000011
22.800000 -0.000010
22.900000 -0.000010
23.000000 -0.000010
23.100000 -0.000010
23.200000 -0.000009
23.300000 -0.000009
23.400000 -0.000009
23.500000 -0.000009
23.600000 -0.000008
23.700000 -0.000008
23.800000 -0.000008
23.900000 -0.000008
24.000000 -0.000008
24.100000 -0.000007
24.200000 -0.000007
24.300000 -0.000007
24.400000 -0.000007
24.500000 -0.000007
24.600000 -0.000007
24.700000 -0.000006
24.800000 -0.000006
24.900000 -0.000006
25.000000 -0.000006
25.100000 -0.000006
25.200000 -0.000006
25.300000 -0.000006
25.400000 -0.000005
25.500000 -0.000005
25.600000 -0.000005
25.700000 -0.000005
25.800000 -0.000005
25.900000 -0.000005
26.000000 -0.000005
26.100000 -0.000005
26.200000 -0.000004
26.300000 -0.000004
26.400000 -0.000004
26.500000 -0.000004
26.600000 -0.000004
26.700000 -0.000004
26.800000 -0.000004
26.900000 -0.000004
27.000000 -0.000004
27.100000 -0.000004
27.200000 -0.000004
27.300000 -0.000003
27.400000 -0.000003
27.500000 -0.000003
27.600000 -0.000003
27.700000 -0.000003
27.800000 -0.000003
27.900000 -0.000003
28.000000 -0.000003
28.100000 -0.000003
28.200000 -0.000003
28.300000 -0.000003
28.400000 -0.000003
28.500000 -0.000003
28.600000 -0.000003
28.700000 -0.000003
28.800000 -0.000002
28.900000 -0.000002
29.000000 -0.000002
29.100000 -0.000002
29.200000 -0.000002
29.300000 -0.000002
29.400000 -0.000002
29.500000 -0.000002
29.600000 -0.000002
29.700000 -0.000002
29.800000 -0.000002
29.900000 -0.000002
30.000000 -0.000002
30.100000 -0.000002
30.200000 -0.000002
30.300000 -0.000002
30.400000 -0.000002
30.500000 -0.000002
30.600000 -0.000002
30.700000 -0.000002
30.800000 -0.000002
30.900000 -0.000002
31.000000 -0.000002
31.100000 -0.000002
31.200000 -0.000001
31.300000 -0.000001
31.400000 -0.000001
31.500000 -0.000001
31.600000 -0.000001
31.700000 -0.000001
31.800000 -0.000001
31.900000 -0.000001
32.000000 -0.000001
32.100000 -0.000001
32.200000 -0.000001
32.300000 -0.000001
32.400000 -0.000001
32.500000 -0.000001
32.600000 -0.000001
32.700000 -0.000001
32.800000 -0.000001
32.900000 -0.000001
33.000000 -0.000001
33.100000 -0.000001
33.200000 -0.000001
33.300000 -0.000001
33.400000 -0.000001
33.500000 -0.000001
33.600000 -0.000001
33.700000 -0.000001
33.800000 -0.000001
33.900000 -0.000001
34.000000 -0.000001
34.100000 -0.000001
34.200000 -0.000001
34.300000 -0.000001
34.400000 -0.000001
34.500000 -0.000001
34.600000 -0.000001
34.700000 -0.000001
34.800000 -0.000001
34.900000 -0.000001
35.000000 -0.000001
35.100000 -0.000001
35.200000 -0.000001
35.300000 -0.000001
35.400000 -0.000001
35.500000 -0.000001
35.600000 -0.000001
35.700000 -0.000001
35.800000 -0.000001
35.900000 -0.000001
36.000000 -0.000001
36.100000 -0.000001
36.200000 -0.000001
36.300000 -0.000001
36.400000 -0.000001
36.500000 -0.000001
36.600000 -0.000001
36.700000 -0.000001
36.800000 -0.000000
36.900000 -0.000000
37.000000 -0.000000
37.100000 -0.000000
37.200000 -0.000000
37.300000 -0.000000
37.400000 -0.000000
37.500000 -0.000000
37.600000 -0.000000
37.700000 -0.000000
37.800000 -0.000000
37.900000 -0.000000
38.000000 -0.000000
38.100000 -0.000000
38.200000 -0.000000
38.300000 -0.000000
38.400000 -0.000000
38.500000 -0.000000
38.600000 -0.000000
38.700000 -0.000000
38.800000 -0.000000
38.900000 -0.000000
39.000000 -0.000000
39.100000 -0.000000
39.200000 -0.000000
39.300000 -0.000000
39.400000 -0.000000
39.500000 -0.000000
39.600000 -0.000000
39.700000 -0.000000
39.800000 -0.000000
39.900000 -0.000000
40.000000 -0.000000
40.100000 -0.000000
40.200000 -0.000000
40.300000 -0.000000
40.400000 -0.000000
40.500000 -0.000000
40.600000 -0.000000
40.700000 -0.000000
40.800000 -0.000000
40.900000 -0.000000
41.000000 -0.000000
41.100000 -0.000000
41.200000 -0.000000
41.300000 -0.000000
41.400000 -0.000000
41.500000 -0.000000
41.600000 -0.000000
41.700000 -0.000000
41.800000 -0.000000
41.900000 -0.000000
42.000000 -0.000000
42.100000 -0.000000
42.200000 -0.000000
42.300000 -0.000000
42.400000 -0.000000
42.500000 -0.000000
42.600000 -0.000000
42.700000 -0.000000
42.800000 -0.000000
42.900000 -0.000000
43.000000 -0.000000
43.100000 -0.000000
43.200000 -0.000000
43.300000 -0.000000
43.400000 -0.000000
43.500000 -0.000000
43.600000 -0.000000
43.700000 -0.000000
43.800000 -0.000000
43.900000 -0.000000
44.000000 -0.000000
44.100000 -0.000000
44.200000 -0.000000
44.300000 -0.000000
44.400000 -0.000000
44.500000 -0.000000
44.600000 -0.000000
44.700000 -0.000000
44.800000 -0.000000
44.900000 -0.000000
45.000000 -0.000000
45.100000 -0.000000
45.200000 -0.000000
45.300000 -0.000000
45.400000 -0.000000
45.500000 -0.000000
45.600000 -0.000000
45.700000 -0.000000
45.800000 -0.000000
45.900000 -0.000000
46.000000 -0.000000
46.100000 -0.000000
46.200000 -0.000000
46.300000 -0.000000
46.400000 -0.000000
46.500000 -0.000000
46.600000 -0.000000
46.700000 -0.000000
46.800000 -0.000000
46.900000 -0.000000
47.000000 -0.000000
47.100000 -0.000000
47.200000 -0.000000
47.300000 -0.000000
47.400000 -0.000000
47.500000 -0.000000
47.600000 -0.000000
47.700000 -0.000000
47.800000 -0.000000
47.900000 -0.000000
48.000000 -0.000000
48.100000 -0.000000
48.200000 -0.000000
48.300000 -0.000000
48.400000 -0.000000
48.500000 -0.000000
48.600000 -0.000000
48.700000 -0.000000
48.800000 -0.000000
48.900000 -0.000000
49.000000 -0.000000
49.100000 -0.000000
49.200000 -0.000000
49.300000 -0.000000
49.400000 -0.000000
49.500000 -0.000000
49.600000 -0.000000
49.700000 -0.000000
49.800000 -0.000000
49.900000 -0.000000
50.000000 0.000000
#! /usr/bin/env python
import numpy as np
eps=0.238000
rmin = 2*1.908100
sig = rmin/(2**(1/6))
r = np.linspace(0,50,501)
r[0] = 0.1
u = 4*eps*( (sig/r)**12 - (sig/r)**6 )
r[0] = 0
ids = np.where( u > 200 )
idMax = ids[0][-1]
du = u[idMax] - u[idMax-1]
u[ids] = u[idMax] -du * (np.arange(idMax,-1,-1) )
u = u-u[-1]
ch = open('dummy-pot.dat','w')
for data in zip(r,u):
ch.write("%f %f\n" % data)
seed 992
timestep 100e-6
steps 100000
temperature 295
outputPeriod 100
outputEnergyPeriod 10
outputFormat dcd
cutoff 100.0
# Potassium
particle POT
num 2
# num 9000
# gridFile null.dx
gridFile null2.dx
diffusionGridFile diffusion.POT.dx
inputCoordinates coord.txt
tabulatedPotential 1
tabulatedFile 0@0@dummy-pot.dat
tabulatedBondFile dummy-pot.dat
inputBonds bondfile.txt
all:
$(MAKE) $(patsubst %.plt.py,%.pdf, $(wildcard *.plt.py))
boltzmann-check.pdf: boltzmann-check.plt.py rho.dat plotTools.py
./$< $@
rho.dat: ../out.0.dcd get-data.tcl
vmd -dispdev text -args $< $@ < get-data.tcl
.SECONDARY:
#! /usr/bin/env python3
# /home/cmaffeo2/bin/python3.sh
# /software/python3-3.4.1/bin/python3
import numpy as np
from scipy import exp
from scipy import sqrt
import matplotlib.pyplot as plt
# import scipy.optimize as opt
from glob import glob
# from natsort import natsorted
import re
from plotTools import *
import sys
def info(*obj):
print('INFO: ',*obj , file=sys.stderr)
def makeTitle(ax,text):
# ax.set_title(text)
ax.annotate(text, xy=(0.05,0.92),
fontsize=9, xycoords='axes fraction',
horizontalalignment='left', verticalalignment='top')
fig,axes = plt.subplots(1)
ax1 = axes
makeTitle(ax1,'Pair interaction (rev. a158e71)')
ax1.set_ylabel('likelihood')
ax1.get_yaxis().set_label_coords(-0.085,0.5)
ax1.set_xlabel("distance between isolated pair of particles")
xmin, xmax = [2,8]
setRange(ax1,[xmin,xmax, 0,0.2])
def loadData(fname):
d0 = np.loadtxt(fname)
x0 = d0[:,0]
y0 = d0[:,1]
ids = np.where( (x0 > xmin) & (x0 < xmax) )
x0 = x0[ids]
y0 = y0[ids]
y0 = y0 / np.sum(y0)
return x0,y0
kT=0.58622592
def estimate():
eps = 2*10
sig = 3
x = np.linspace(xmin,xmax,100)
u = 4*eps*( (sig/x)**12 - (sig/x)**6 )
y = np.exp( - u/kT ) * x**2
y = y / np.sum(y)
return x,y
def estimate():
d0 = np.loadtxt('../dummy-pot.dat')
x = d0[10:,0]
u = d0[10:,1]
u = 2*u
ids = np.where( (x > xmin) & (x < xmax) )
x = x[ids]
u = u[ids]
y = np.exp( - u/kT ) * x**2
y = y / np.sum(y)
return x,y
def plot(x,y,l,c):
# return ax1.plot(x,y, ls='-', marker='o', color=c, mec=c)
return ax1.bar(x,y, label=l, lw=0, width=np.mean(x[1:]-x[:-1]), alpha=0.5, color=c)
x,y = loadData('rho.dat')
ex,ey = estimate()
plot(x,y,'simulated', c0)
plot(ex,ey,'expected',c1)
plt.legend()
# plt.show()
plt.tight_layout()
# plt.savefig(sys.stdout, format='svg')
plt.savefig(sys.argv[1])
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment