mirror of https://github.com/QMCPACK/qmcpack.git
107 lines
2.7 KiB
Python
Executable File
107 lines
2.7 KiB
Python
Executable File
#!/bin/env python
|
|
from numpy import *
|
|
import sys
|
|
|
|
|
|
def corr(i,x,mean,var):
|
|
N=len(x)
|
|
if var==0:#if the variance is 0 return an effectively infinity corr
|
|
return 1e100
|
|
corr=1.0/var*1.0/(N-i)*sum((x[0:N-i]-mean)*(x[i:N]-mean))
|
|
return corr
|
|
|
|
def Stats(x):
|
|
N=len(x)
|
|
mean=sum(x)/(N+0.0)
|
|
xSquared=x*x
|
|
var=sum(xSquared)/(N+0.0)-mean*mean
|
|
i=0
|
|
tempC=0.5
|
|
kappa=0.0
|
|
while (tempC>0 and i<(N-1)):
|
|
kappa=kappa+2.0*tempC
|
|
i=i+1
|
|
tempC=corr(i,x,mean,var)
|
|
if kappa == 0.0:
|
|
kappa = 1.0
|
|
Neff=(N+0.0)/(kappa+0.0)
|
|
error=sqrt(var/Neff)
|
|
return (mean,var,error,kappa)
|
|
|
|
def Block (x, blockfactor):
|
|
N = len(x)
|
|
Nblocks = N / blockfactor
|
|
xb = zeros(Nblocks)
|
|
for i in range(0,Nblocks):
|
|
start = i*blockfactor
|
|
end = (i+1)*blockfactor
|
|
xb[i] = sum(x[start:end])/(blockfactor+0.0)
|
|
return xb
|
|
|
|
def MeanErrorString (mean, error):
|
|
if (mean!=0.0):
|
|
meanDigits = math.floor(math.log(abs(mean))/math.log(10))
|
|
else:
|
|
meanDigits=2
|
|
if (error!=0.0):
|
|
rightDigits = -math.floor(math.log(error)/math.log(10))+1
|
|
else:
|
|
rightDigits=2
|
|
if (rightDigits < 0):
|
|
rightDigits = 0
|
|
formatstr = '%1.' + '%d' % rightDigits + 'f'
|
|
meanstr = formatstr % mean
|
|
errorstr = formatstr % error
|
|
return (meanstr, errorstr)
|
|
# return (meanstr, errorstr)
|
|
|
|
|
|
|
|
file = open (sys.argv[1], 'r')
|
|
names = file.readline().split()
|
|
file.close()
|
|
|
|
numat = 0
|
|
for n in names:
|
|
if (n[0] == 'F' and n[1] != 'W'):
|
|
d = n.split('_')
|
|
iat = int(d[-2])
|
|
dim = int(d[-1])
|
|
numat = max(iat+1,numat)
|
|
collist = [ [ [], [], [] ] ]
|
|
for iat in range(0,numat-1):
|
|
collist.append([ [],[],[] ])
|
|
|
|
col = 0
|
|
for n in names:
|
|
if (n[0] == 'F' and n[1] != 'W'):
|
|
d = n.split('_')
|
|
iat = int(d[-2])
|
|
dim = int(d[-1])
|
|
numat = max(iat+1,numat)
|
|
collist[iat][dim].append(col-1)
|
|
col = col + 1
|
|
|
|
s = loadtxt(sys.argv[1])
|
|
c = s.shape[0]
|
|
if len(sys.argv) > 2:
|
|
first = int(sys.argv[2])
|
|
else:
|
|
first = 20
|
|
|
|
print 'Atom Force +/- Error'
|
|
for iat in range(0,numat):
|
|
means = []
|
|
errors = []
|
|
for dim in range(0,3):
|
|
Ftot = 0.0*array(c-first)
|
|
for col in collist[iat][dim]:
|
|
# print 'iat = %d dim = %d col=%d' % (iat, dim, col)
|
|
Ftot = Ftot + s[first:c+1,col]
|
|
(avg, var, error, kapp) = Stats(Ftot)
|
|
means.append (avg)
|
|
errors.append(error)
|
|
# print '%-20s = %s' % (n, MeanErrorString(avg,error))
|
|
print '%3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f' \
|
|
% (iat+1,means[0], means[1], means[2], errors[0], errors[1], errors[2])
|