10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
quantum_package/plugins/QMC/qp_convert_qmcpack_to_ezfio.py

370 lines
8.6 KiB
Python
Raw Normal View History

2016-01-19 17:06:57 +01:00
#!/usr/bin/python
print "#QP -> QMCPACK"
2016-01-27 10:48:10 +01:00
# ___
# | ._ o _|_
# _|_ | | | |_
#
2016-01-27 10:48:10 +01:00
2016-01-19 17:06:57 +01:00
from ezfio import ezfio
2016-03-07 10:51:11 +01:00
import os
2016-01-19 17:06:57 +01:00
import sys
ezfio_path = sys.argv[1]
ezfio.set_file(ezfio_path)
do_pseudo = ezfio.get_pseudo_do_pseudo()
if do_pseudo:
2016-01-27 10:48:10 +01:00
print "do_pseudo True"
2016-03-07 10:51:11 +01:00
from qp_path import QP_ROOT
l_ele_path = os.path.join(QP_ROOT,"data","list_element.txt")
with open(l_ele_path, "r") as f:
data_raw = f.read()
l_element_raw = data_raw.split("\n")
l_element = [element_raw.split() for element_raw in l_element_raw]
d_z = dict((abr, z) for (z, abr, ele, _) in filter(lambda x: x != [], l_element) )
2016-01-19 17:06:57 +01:00
else:
2016-01-27 10:48:10 +01:00
print "do_pseudo False"
2016-01-19 17:06:57 +01:00
try:
2016-01-27 10:48:10 +01:00
n_det = ezfio.get_determinants_n_det()
2016-01-19 17:06:57 +01:00
except IOError:
2016-01-27 10:48:10 +01:00
n_det = 1
2016-01-19 17:06:57 +01:00
if n_det == 1:
2016-01-27 10:48:10 +01:00
print "multi_det False"
2016-01-19 17:06:57 +01:00
else:
2016-01-27 10:48:10 +01:00
print "multi_det True"
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
#
# |\/| o _ _
# | | | _> (_
#
def list_to_string(l):
return " ".join(map(str, l))
2016-01-19 17:06:57 +01:00
ao_num = ezfio.get_ao_basis_ao_num()
print "ao_num", ao_num
mo_num = ezfio.get_mo_basis_mo_tot_num()
print "mo_num", mo_num
alpha = ezfio.get_electrons_elec_alpha_num()
beta = ezfio.get_electrons_elec_beta_num()
print "elec_alpha_num", alpha
print "elec_beta_num", beta
2016-01-27 10:48:10 +01:00
print "elec_tot_num", alpha + beta
print "spin_multiplicity", 2 * (alpha - beta) + 1
2016-01-19 17:06:57 +01:00
l_label = ezfio.get_nuclei_nucl_label()
l_charge = ezfio.get_nuclei_nucl_charge()
l_coord = ezfio.get_nuclei_nucl_coord()
l_coord_str = [list_to_string(i) for i in zip(*l_coord)]
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
print "nucl_num", len(l_label)
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
# _
# / _ _ ._ _|
# \_ (_) (_) | (_|
#
2016-01-27 10:48:10 +01:00
print "Atomic coord in Bohr"
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
for i, t in enumerate(zip(l_label, l_charge, l_coord_str)):
2016-03-07 10:51:11 +01:00
t_1 = d_z[t[0]] if do_pseudo else t[1]
t_new = [t[0],t_1,t[2]]
print list_to_string(t_new)
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
#
# Call externet process to get the sysmetry
#
2016-01-19 17:06:57 +01:00
import subprocess
2016-01-27 10:48:10 +01:00
process = subprocess.Popen(
['qp_print_basis', ezfio_path],
stdout=subprocess.PIPE)
2016-01-19 17:06:57 +01:00
out, err = process.communicate()
2017-07-15 05:27:35 +02:00
basis_raw, sym_raw, _ = out.split("\n\n\n")
2016-01-27 10:48:10 +01:00
# _ __
# |_) _. _ o _ (_ _ _|_
# |_) (_| _> | _> __) (/_ |_
#
2016-01-19 17:06:57 +01:00
2016-02-04 18:50:40 +01:00
basis_without_header = "\n".join(basis_raw.split("\n")[11:])
2016-01-19 17:06:57 +01:00
import re
l_basis_raw = re.split('\n\s*\n', basis_without_header)
a_already_print = []
l_basis_clean = []
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
for i, (a,b) in enumerate(zip(l_label,l_basis_raw)):
if a not in a_already_print:
l_basis_clean.append(b.replace('Atom {0}'.format(i + 1), a))
a_already_print.append(a)
2017-07-15 05:47:55 +02:00
else:
continue
print "BEGIN_BASIS_SET\n"
print "\n\n".join(l_basis_clean)
print "END_BASIS_SET"
2016-01-19 17:06:57 +01:00
# _
# |\/| / \ _
# | | \_/ _>
2016-01-27 10:48:10 +01:00
#
2016-01-27 10:48:10 +01:00
#
# Function
#
2017-07-15 05:47:55 +02:00
d_gms_order ={ 0:["s"],
2017-07-15 02:20:10 +02:00
1:[ "x", "y", "z" ],
2:[ "xx", "yy", "zz", "xy", "xz", "yz" ],
3:[ "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz"],
2017-07-15 05:47:55 +02:00
4: ["xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy", "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz","zzxy"] }
2016-01-19 17:06:57 +01:00
def compare_gamess_style(item1, item2):
2017-07-15 05:47:55 +02:00
n1,n2 = map(len,(item1,item2))
assert (n1 == n2)
try:
l = d_gms_order[n1]
except KeyError:
return 0
# raise (KeyError, "We dont handle L than 4")
else:
a = l.index(item1)
b = l.index(item2)
return cmp( a, b )
def expend_sym_str(str_):
#Expend x2 -> xx
# yx2 -> xxy
2016-01-27 10:48:10 +01:00
for i, c in enumerate(str_):
try:
n = int(c)
except ValueError:
pass
else:
str_ = str_[:i - 1] + str_[i - 1] * n + str_[i + 1:]
#Order by frequency
2017-07-15 05:47:55 +02:00
return "".join(sorted(str_, key=str_.count, reverse=True))
def expend_sym_l(l_l_sym):
for l in l_l_sym:
l[2] = expend_sym_str(l[2])
return l_l_sym
def n_orbital(n):
if n==0:
return 1
elif n==1:
return 3
else:
return 2*n_orbital(n-1)-n_orbital(n-2)+1
def get_nb_permutation(str_):
if (str_) == 's': return 1
else: return n_orbital(len(str_))
def order_l_l_sym(l_l_sym):
n = 1
iter_ = range(len(l_l_sym))
for i in iter_:
if n != 1:
n += -1
continue
l = l_l_sym[i]
n = get_nb_permutation(l[2])
l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n],
key=lambda x: x[2],
cmp=compare_gamess_style)
return l_l_sym
2016-01-27 10:48:10 +01:00
#==========================
# We will order the symetry
#==========================
2017-07-15 03:09:13 +02:00
2016-01-19 17:06:57 +01:00
l_sym_without_header = sym_raw.split("\n")[3:-2]
l_l_sym_raw = [i.split() for i in l_sym_without_header]
2017-07-15 05:47:55 +02:00
l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw)
l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym)
2016-01-27 10:48:10 +01:00
#========
#MO COEF
#========
def order_phase(mo_coef):
#Order
mo_coef_phase = []
import math
for i in mo_coef:
if abs(max(i)) > abs(min(i)):
sign_max = math.copysign(1, max(i))
else:
sign_max = math.copysign(1, min(i))
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
if sign_max == -1:
ii = [-1 * l for l in i]
else:
ii = i
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
mo_coef_phase.append(ii)
return mo_coef_phase
2016-01-19 17:06:57 +01:00
def chunked(l, chunks_size):
l_block = []
for i in l:
2016-01-27 10:48:10 +01:00
chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)]
l_block.append(chunks)
return l_block
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
def print_mo_coef(mo_coef_block, l_l_sym):
print ""
print "BEGIN_MO"
print ""
len_block_curent = 0
nb_block = len(mo_coef_block[0])
for i_block in range(0, nb_block):
a = [i[i_block] for i in mo_coef_block]
r_ = range(len_block_curent, len_block_curent + len(a[0]))
print " ".join([str(i + 1) for i in r_])
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
len_block_curent += len(a[0])
2016-01-19 17:06:57 +01:00
2016-01-27 10:48:10 +01:00
for l in l_l_sym:
i = int(l[0]) - 1
i_a = int(l[1]) - 1
sym = l[2]
2017-07-20 01:35:00 +02:00
print l_label[i_a], sym, " ".join('%20.15e'%i
2016-01-27 10:48:10 +01:00
for i in a[i])
if i_block != nb_block - 1:
print ""
else:
print "END_MO"
mo_coef = ezfio.get_mo_basis_mo_coef()
mo_coef_transp = zip(*mo_coef)
2016-01-27 10:48:10 +01:00
mo_coef_block = chunked(mo_coef_transp, 4)
print_mo_coef(mo_coef_block, l_l_sym_ordered)
2016-01-27 10:48:10 +01:00
# _
# |_) _ _ _| _
# | _> (/_ |_| (_| (_)
#
2016-01-19 17:06:57 +01:00
if do_pseudo:
2016-01-27 10:48:10 +01:00
print ""
print "BEGIN_PSEUDO"
klocmax = ezfio.get_pseudo_pseudo_klocmax()
kmax = ezfio.get_pseudo_pseudo_kmax()
lmax = ezfio.get_pseudo_pseudo_lmax()
n_k = ezfio.get_pseudo_pseudo_n_k()
v_k = ezfio.get_pseudo_pseudo_v_k()
dz_k = ezfio.get_pseudo_pseudo_dz_k()
n_kl = ezfio.get_pseudo_pseudo_n_kl()
v_kl = ezfio.get_pseudo_pseudo_v_kl()
dz_kl = ezfio.get_pseudo_pseudo_dz_kl()
for i, a in enumerate(l_label):
l_str = []
#Local
2016-01-27 10:48:10 +01:00
l_dump = []
for k in range(klocmax):
if v_k[k][i]:
l_ = list_to_string([v_k[k][i], n_k[k][i] + 2, dz_k[k][i]])
l_dump.append(l_)
l_str.append(l_dump)
#Non local
2016-01-27 10:48:10 +01:00
for l in range(lmax + 1):
l_dump = []
for k in range(kmax):
if v_kl[l][k][i]:
l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i] + 2,
dz_kl[l][k][i]])
l_dump.append(l_)
if l_dump:
l_str.append(l_dump)
2016-03-07 10:51:11 +01:00
str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE"
print str_.format(a, i + 1, int(d_z[a])-int(l_charge[i]), int(len(l_str) - 1))
2016-01-27 10:48:10 +01:00
for i, l in enumerate(l_str):
str_ = "FOR L= {0} COEFF N ZETA"
print str_.format(int(len(l_str) - i - 1))
for ii, ll in enumerate(l):
print " ", ii + 1, ll
2016-03-07 10:51:11 +01:00
str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS."
print str_.format(sum([int(d_z[a])-int(l_charge[i]) for i,a in enumerate(l_label)]))
2016-01-27 10:48:10 +01:00
print "END_PSEUDO"
# _
# | \ _ _|_
# |_/ (/_ |_
#
2017-07-15 05:27:35 +02:00
psi_coef = ezfio.get_determinants_psi_coef()
psi_det = ezfio.get_determinants_psi_det()
bit_kind = ezfio.get_determinants_bit_kind()
2016-01-19 17:06:57 +01:00
print ""
print "BEGIN_DET"
print ""
print "mo_num", mo_num
print "det_num", n_det
print ""
if "QP_STATE" in os.environ:
state = int(os.environ["QP_STATE"])-1
else:
state = 0
2017-07-15 05:27:35 +02:00
psi_coef = psi_coef[state]
2016-02-29 11:26:33 +01:00
2017-07-15 05:27:35 +02:00
encode = 8*bit_kind
2016-02-29 11:26:33 +01:00
2017-07-15 18:06:20 +02:00
def bindigits(n, bits):
s = bin(n & int("1"*bits, 2))[2:]
return ("{0:0>%s}" % (bits)).format(s)
decode = lambda det: ''.join(bindigits(i,encode)[::-1] for i in det)[:mo_num]
2017-07-15 05:27:35 +02:00
for coef, (det_a, det_b) in zip(psi_coef, psi_det):
2016-02-29 11:26:33 +01:00
2017-07-15 05:27:35 +02:00
print coef
print decode(det_a)
print decode(det_b)
print ''
2016-03-01 15:42:16 +01:00
print "END_DET"