mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-12-23 12:54:35 +01:00
262 lines
9.3 KiB
Python
262 lines
9.3 KiB
Python
|
# !!!
|
||
|
|
import sys, os
|
||
|
|
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
|
||
|
|
sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
|
||
|
|
|
||
|
|
|
||
|
|
from ezfio import ezfio
|
||
|
|
from math import sqrt
|
||
|
|
from datetime import datetime
|
||
|
|
import time
|
||
|
|
import numpy as np
|
||
|
|
import subprocess
|
||
|
|
from scipy.linalg import eig, eigh
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
def get_energy():
|
||
|
|
buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', EZFIO_file], encoding='UTF-8')
|
||
|
|
if buffer.strip() != "":
|
||
|
|
buffer = buffer.splitlines()[-1]
|
||
|
|
_, energy, error = [float(x) for x in buffer.split()]
|
||
|
|
return energy, error
|
||
|
|
else:
|
||
|
|
return None, None
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
def get_Ci_h_matrix_svd():
|
||
|
|
|
||
|
|
Ci_h_matrix_svd = np.zeros( (n_svd,n_svd) )
|
||
|
|
SErr_H = np.zeros( (n_svd,n_svd) )
|
||
|
|
|
||
|
|
beg_Ci_h_matrix_svd = results.find('Ci_h_matrix_svd : [ ') + len( 'Ci_h_matrix_svd : [ ' )
|
||
|
|
end_Ci_h_matrix_svd = len(results)
|
||
|
|
Ci_h_matrix_svd_buf = results[beg_Ci_h_matrix_svd:end_Ci_h_matrix_svd]
|
||
|
|
Ci_h_matrix_svd_buf = Ci_h_matrix_svd_buf.split( '\n' )
|
||
|
|
|
||
|
|
for iline in range(1, n_svd**2+1):
|
||
|
|
|
||
|
|
line = Ci_h_matrix_svd_buf[iline].split()
|
||
|
|
indc = int( line[0] )
|
||
|
|
errS = float( line[4] )
|
||
|
|
if( indc != iline ):
|
||
|
|
print(' Error in reading Ci_h_matrix_svd')
|
||
|
|
print(' indc = {}'.format(indc) )
|
||
|
|
print(' iline = {}'.format(iline) )
|
||
|
|
break
|
||
|
|
else:
|
||
|
|
irow = indc % n_svd
|
||
|
|
icol = indc // n_svd
|
||
|
|
if( irow!=0 ):
|
||
|
|
Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
|
||
|
|
SErr_H [irow-1][icol] = float( line[4] )
|
||
|
|
|
||
|
|
#ci_tmp = Ci_h_matrix_svd[irow-1][icol]
|
||
|
|
#er_tmp = SErr_H [irow-1][icol]
|
||
|
|
#if(ci_tmp>0):
|
||
|
|
# Ci_h_matrix_svd[irow-1][icol] = max(0 , ci_tmp - er_tmp)
|
||
|
|
#else:
|
||
|
|
# Ci_h_matrix_svd[irow-1][icol] = min(0 , ci_tmp + er_tmp)
|
||
|
|
|
||
|
|
else:
|
||
|
|
Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
|
||
|
|
SErr_H [n_svd-1][icol-1] = float( line[4] )
|
||
|
|
|
||
|
|
#ci_tmp = Ci_h_matrix_svd[n_svd-1][icol-1]
|
||
|
|
#er_tmp = SErr_H [n_svd-1][icol-1]
|
||
|
|
#if(ci_tmp>0):
|
||
|
|
# Ci_h_matrix_svd[n_svd-1][icol-1] = max(0 , ci_tmp - er_tmp)
|
||
|
|
#else:
|
||
|
|
# Ci_h_matrix_svd[n_svd-1][icol-1] = min(0 , ci_tmp + er_tmp)
|
||
|
|
|
||
|
|
return(Ci_h_matrix_svd,SErr_H)
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
def get_Ci_overlap_matrix_svd():
|
||
|
|
|
||
|
|
Ci_overlap_matrix_svd = np.zeros( (n_svd,n_svd) )
|
||
|
|
SErr_O = np.zeros( (n_svd,n_svd) )
|
||
|
|
|
||
|
|
beg_Ci_overlap_matrix_svd = results.find('Ci_overlap_matrix_svd : [ ') + len( 'Ci_overlap_matrix_svd : [ ' )
|
||
|
|
end_Ci_overlap_matrix_svd = len(results)
|
||
|
|
Ci_overlap_matrix_svd_buf = results[beg_Ci_overlap_matrix_svd:end_Ci_overlap_matrix_svd]
|
||
|
|
Ci_overlap_matrix_svd_buf = Ci_overlap_matrix_svd_buf.split( '\n' )
|
||
|
|
|
||
|
|
for iline in range(1, n_svd**2+1):
|
||
|
|
|
||
|
|
line = Ci_overlap_matrix_svd_buf[iline].split()
|
||
|
|
indc = int( line[0] )
|
||
|
|
errS = float( line[4] )
|
||
|
|
if( indc != iline ):
|
||
|
|
print(' Error in reading Ci_overlap_matrix_svd')
|
||
|
|
print(' indc = {}'.format(indc))
|
||
|
|
print(' iline = {}'.format(iline))
|
||
|
|
break
|
||
|
|
else:
|
||
|
|
irow = indc % n_svd
|
||
|
|
icol = indc // n_svd
|
||
|
|
if( irow!=0 ):
|
||
|
|
Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
|
||
|
|
SErr_O [irow-1][icol] = float( line[4] )
|
||
|
|
|
||
|
|
#ci_tmp = Ci_overlap_matrix_svd[irow-1][icol]
|
||
|
|
#er_tmp = SErr_O [irow-1][icol]
|
||
|
|
#if(ci_tmp>0):
|
||
|
|
# Ci_overlap_matrix_svd[irow-1][icol] = max(0 , ci_tmp - er_tmp)
|
||
|
|
#else:
|
||
|
|
# Ci_overlap_matrix_svd[irow-1][icol] = min(0 , ci_tmp + er_tmp)
|
||
|
|
|
||
|
|
else:
|
||
|
|
Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
|
||
|
|
SErr_O [n_svd-1][icol-1] = float( line[4] )
|
||
|
|
|
||
|
|
#ci_tmp = Ci_overlap_matrix_svd[n_svd-1][icol-1]
|
||
|
|
#er_tmp = SErr_O [n_svd-1][icol-1]
|
||
|
|
#if(ci_tmp>0):
|
||
|
|
# Ci_overlap_matrix_svd[n_svd-1][icol-1] = max(0 , ci_tmp - er_tmp)
|
||
|
|
#else:
|
||
|
|
# Ci_overlap_matrix_svd[n_svd-1][icol-1] = min(0 , ci_tmp + er_tmp)
|
||
|
|
|
||
|
|
# !!!
|
||
|
|
# !!!
|
||
|
|
# !!!
|
||
|
|
|
||
|
|
return(Ci_overlap_matrix_svd,SErr_O)
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
def check_symmetric(a, tol=1e-3):
|
||
|
|
return np.all(np.abs(a-a.T) < tol)
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
def get_Esvd():
|
||
|
|
|
||
|
|
# read CI_SVD matrices
|
||
|
|
Ci_h_matrix_svd , SErr_H = get_Ci_h_matrix_svd()
|
||
|
|
Ci_overlap_matrix_svd, SErr_O = get_Ci_overlap_matrix_svd()
|
||
|
|
|
||
|
|
#print( 'Ci_h_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_svd)) )
|
||
|
|
#print( 'Ci_overlap_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_svd)) )
|
||
|
|
|
||
|
|
if( save_SVD_QMCmat ):
|
||
|
|
with open('SVD_QMCmat.txt', 'a') as ff:
|
||
|
|
for i in range(n_svd):
|
||
|
|
for j in range(n_svd):
|
||
|
|
ff.write('{} {} {} {} {} {} \n'.format(i,j
|
||
|
|
, Ci_h_matrix_svd[i,j] , SErr_H[i,j]
|
||
|
|
, Ci_overlap_matrix_svd[i,j] , SErr_O[i,j] ))
|
||
|
|
|
||
|
|
# symmetrise and diagonalise
|
||
|
|
aa = Ci_h_matrix_svd
|
||
|
|
aa = 0.5*( aa + aa.T )
|
||
|
|
bb = Ci_overlap_matrix_svd
|
||
|
|
eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True,
|
||
|
|
check_finite=True, homogeneous_eigvals=False)
|
||
|
|
|
||
|
|
#print( eigvals_svd + E_toadd )
|
||
|
|
recouvre_svd = np.abs(psi_svd_coeff @ vr)
|
||
|
|
ind_gssvd = np.argmax(recouvre_svd)
|
||
|
|
|
||
|
|
E_svd = eigvals_svd[ind_gssvd] + E_toadd
|
||
|
|
|
||
|
|
return( E_svd, vr[:,ind_gssvd] )
|
||
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == '__main__':
|
||
|
|
|
||
|
|
t0 = time.time()
|
||
|
|
|
||
|
|
# ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
|
||
|
|
#EZFIO_file = "Cr2_TSVD30"
|
||
|
|
#EZFIO_file = "Cr2_TSVD100"
|
||
|
|
#EZFIO_file = "Cr2_TSVD200"
|
||
|
|
#EZFIO_file = "Cr2_TSVD924_0"
|
||
|
|
#EZFIO_file = "Cr2_TSVD924_1"
|
||
|
|
#EZFIO_file = "Cr2_TSVD924_2"
|
||
|
|
|
||
|
|
#EZFIO_file = "Cr2_TSVD30_it1"
|
||
|
|
#EZFIO_file = "Cr2_TSVD100_it1"
|
||
|
|
#EZFIO_file = "Cr2_TSVD200_it1"
|
||
|
|
#EZFIO_file = "Cr2_TSVD924_0_it1"
|
||
|
|
#EZFIO_file = "Cr2_TSVD924_1_it1"
|
||
|
|
EZFIO_file = "Cr2_TSVD924_2_it1"
|
||
|
|
|
||
|
|
#EZFIO_file = "Cr2_TSVD200_minmaxit1"
|
||
|
|
|
||
|
|
#EZFIO_file = "Cr2_TSVD30_it2"
|
||
|
|
|
||
|
|
E_toadd = 64.8242130309350
|
||
|
|
# ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
|
||
|
|
|
||
|
|
ezfio.set_file(EZFIO_file)
|
||
|
|
n_svd = ezfio.get_spindeterminants_n_svd_coefs()
|
||
|
|
psi_svd_coeff = np.array(ezfio.get_spindeterminants_psi_svd_coefs())
|
||
|
|
|
||
|
|
print(" ")
|
||
|
|
print(" Today's date:", datetime.now() )
|
||
|
|
print(" EZFIO file = {}".format(EZFIO_file))
|
||
|
|
print(" n_svd = {}\n".format(n_svd) )
|
||
|
|
|
||
|
|
t_read = time.time()
|
||
|
|
print(' getting QMCCHEM results from {}'.format(EZFIO_file) )
|
||
|
|
results = subprocess.check_output(['qmcchem', 'result', EZFIO_file], encoding='UTF-8')
|
||
|
|
print(' getting results after {} minutes \n'.format( (time.time()-t_read)/60. ))
|
||
|
|
|
||
|
|
if( len(results)<100):
|
||
|
|
print(' result is empty \n')
|
||
|
|
exit()
|
||
|
|
|
||
|
|
Eloc, ErrEloc = get_energy()
|
||
|
|
print(' Eloc = {} +/- {}'.format(Eloc, ErrEloc))
|
||
|
|
print(' {} <= Eloc <= {}'.format(Eloc-ErrEloc,Eloc+ErrEloc))
|
||
|
|
#print(' initial svd coeff = {}\n'.format(psi_svd_coeff))
|
||
|
|
|
||
|
|
save_SVD_QMCmat = True
|
||
|
|
#save_SVD_QMCmat = False
|
||
|
|
|
||
|
|
print('')
|
||
|
|
E_svd, sigma_svd = get_Esvd()
|
||
|
|
print(' svd enegry = {}'.format(E_svd) )
|
||
|
|
sigma_svd = sigma_svd * np.sign(sigma_svd[0])
|
||
|
|
#print(' svd coeff = {}'.format(sigma_svd))
|
||
|
|
print('')
|
||
|
|
|
||
|
|
# ___________________________________________________________________________
|
||
|
|
# ___________________________________________________________________________
|
||
|
|
#
|
||
|
|
save_to_EZFIO = 'svd_QMC'
|
||
|
|
#save_to_EZFIO = ''
|
||
|
|
|
||
|
|
if( save_to_EZFIO == 'svd_QMC' ):
|
||
|
|
ezfio.set_spindeterminants_psi_svd_coefs( sigma_svd )
|
||
|
|
direc_svdcoeff = EZFIO_file + '/spindeterminants/psi_svd_coefs.gz'
|
||
|
|
print(' {} is modified'.format(direc_svdcoeff) )
|
||
|
|
else:
|
||
|
|
print(' exit with no modifications in {} EZFIO'.format(EZFIO_file))
|
||
|
|
# ___________________________________________________________________________
|
||
|
|
# ___________________________________________________________________________
|
||
|
|
#
|
||
|
|
|
||
|
|
print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
|
||
|
|
|
||
|
|
# !!!
|
||
|
|
# !!!
|
||
|
|
|
||
|
|
~
|