mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-06-01 10:55:18 +02:00
263 lines
9.7 KiB
Python
Executable File
263 lines
9.7 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
# !!!
|
|
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
|
|
|
|
# PARAMETERS
|
|
block_time = 20
|
|
eps = 1.
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def get_energy():
|
|
buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename], 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 run_qmc():
|
|
return subprocess.check_output(['qmcchem', 'run', filename])
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def stop_qmc():
|
|
subprocess.check_output(['qmcchem', 'stop', filename])
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def set_vmc_params():
|
|
#subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
|
|
# '-m', 'VMC',
|
|
# '-l', str(20),
|
|
# '--time-step=0.3',
|
|
# '--stop-time=36000',
|
|
# '--norm=1.e-5',
|
|
# '-w', '10',
|
|
# filename])
|
|
subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'None', '-l', str(block_time), filename])
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def get_Ci_h_matrix_svd(n_svd):
|
|
# !!!
|
|
Ci_h_matrix_svd = 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)#results.find('E_loc :')
|
|
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( errS>eps ):
|
|
#print( line )
|
|
if( indc != iline ):
|
|
print('Error in reading Ci_h_matrix_svd')
|
|
stop
|
|
else:
|
|
#Ci_h_matrix_svd[indc-1] = float( line[2] )
|
|
irow = indc % n_svd
|
|
icol = indc // n_svd
|
|
if( irow!=0 ):
|
|
Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
|
|
else:
|
|
Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
|
|
# !!!
|
|
# !!!
|
|
# Ci_h_matrix_svd = np.reshape(Ci_h_matrix_svd, (n_svd, n_svd), order='F')
|
|
# !!!
|
|
return(Ci_h_matrix_svd)
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def get_Ci_overlap_matrix_svd(n_svd):
|
|
# !!!
|
|
Ci_overlap_matrix_svd = 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)#results.find('Ci_h_matrix_svd : [')
|
|
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( errS>eps ):
|
|
#print( line )
|
|
if( indc != iline ):
|
|
print('Error in reading Ci_overlap_matrix_svd')
|
|
stop
|
|
# !!!
|
|
# !!!
|
|
else:
|
|
#Ci_overlap_matrix_svd[indc-1] = float( line[2] )
|
|
irow = indc % n_svd
|
|
icol = indc // n_svd
|
|
if( irow!=0 ):
|
|
Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
|
|
else:
|
|
Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
|
|
# !!!
|
|
# !!!
|
|
# !!!
|
|
#Ci_overlap_matrix_svd = np.reshape(Ci_overlap_matrix_svd, (n_svd, n_svd), order='F')
|
|
# !!!
|
|
return(Ci_overlap_matrix_svd)
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def get_Ci_h_matrix_postsvd(n_svd):
|
|
# !!!
|
|
Ci_h_matrix_postsvd = np.zeros( (n_svd*n_svd*n_svd*n_svd) )
|
|
# !!!
|
|
beg_Ci_h_matrix_postsvd = results.find('Ci_h_matrix_postsvd : [ ') + len( 'Ci_h_matrix_postsvd : [ ' )
|
|
end_Ci_h_matrix_postsvd = len(results)
|
|
Ci_h_matrix_postsvd_buf = results[beg_Ci_h_matrix_postsvd:end_Ci_h_matrix_postsvd]
|
|
Ci_h_matrix_postsvd_buf = Ci_h_matrix_postsvd_buf.split( '\n' )
|
|
# !!!
|
|
for iline in range(1, n_svd**4+1):
|
|
# !!!
|
|
line = Ci_h_matrix_postsvd_buf[iline].split()
|
|
indc = int( line[0] )
|
|
errS = float( line[4] )
|
|
#if( errS>eps ):
|
|
#print( line )
|
|
if( indc != iline ):
|
|
print('Error in reading Ci_h_matrix_postsvd')
|
|
stop
|
|
else:
|
|
Ci_h_matrix_postsvd[indc-1] = float( line[2] )
|
|
# !!!
|
|
Ci_h_matrix_postsvd = np.reshape(Ci_h_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
|
|
# !!!
|
|
return(Ci_h_matrix_postsvd)
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
def get_Ci_overlap_matrix_postsvd(n_svd):
|
|
# !!!
|
|
Ci_overlap_matrix_postsvd = np.zeros( (n_svd*n_svd*n_svd*n_svd) )
|
|
# !!!
|
|
beg_Ci_overlap_matrix_postsvd = results.find('Ci_overlap_matrix_postsvd : [ ') + len( 'Ci_overlap_matrix_postsvd : [ ' )
|
|
end_Ci_overlap_matrix_postsvd = len(results)
|
|
Ci_overlap_matrix_postsvd_buf = results[beg_Ci_overlap_matrix_postsvd:end_Ci_overlap_matrix_postsvd]
|
|
Ci_overlap_matrix_postsvd_buf = Ci_overlap_matrix_postsvd_buf.split( '\n' )
|
|
# !!!
|
|
for iline in range(1, n_svd**4+1):
|
|
# !!!
|
|
line = Ci_overlap_matrix_postsvd_buf[iline].split()
|
|
indc = int( line[0] )
|
|
errS = float( line[4] )
|
|
#if( errS>eps ):
|
|
#print( line )
|
|
if( indc != iline ):
|
|
print('Error in reading Ci_overlap_matrix_postsvd')
|
|
stop
|
|
else:
|
|
Ci_overlap_matrix_postsvd[indc-1] = float( line[2] )
|
|
# !!!
|
|
Ci_overlap_matrix_postsvd = np.reshape(Ci_overlap_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
|
|
# !!!
|
|
return(Ci_overlap_matrix_postsvd)
|
|
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
t0 = time.time()
|
|
# !!!
|
|
filename = "/home/ammar/qp2/src/svdwf/h2o.ezfio"
|
|
ezfio.set_file(filename)
|
|
# !!!
|
|
print("Today's date:", datetime.now() )
|
|
print("filename = {}".format(filename))
|
|
print("start QMC:")
|
|
#set_vmc_params()
|
|
#run_qmc()
|
|
#
|
|
##stop_qmc()
|
|
results = subprocess.check_output(['qmcchem', 'result', filename], encoding='UTF-8')
|
|
# !!!
|
|
#file = open('results.txt','a')
|
|
#file.write('\n \n \n')
|
|
#file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
|
|
#file.write("Today's date: {}\n".format(datetime.now()) )
|
|
#file.write("filename = {}\n".format(filename))
|
|
#file.write('\n')
|
|
#file.write( results )
|
|
#file.write('\n')
|
|
#print("end QMC after {} minutes\n".format((time.time()-t0)/60.) )
|
|
#file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
|
|
#file.close()
|
|
# !!!
|
|
E_loc, ErrEloc = get_energy()
|
|
print('Eloc = {} +/- {}'.format(E_loc, ErrEloc))
|
|
# !!!
|
|
n_svd = 10
|
|
Ci_h_matrix_svd = get_Ci_h_matrix_svd(n_svd)
|
|
Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd(n_svd)
|
|
# !!!
|
|
aa = Ci_h_matrix_svd
|
|
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+9.194966082434476 )
|
|
# !!!
|
|
Ci_h_matrix_postsvd = get_Ci_h_matrix_postsvd(n_svd)
|
|
Ci_overlap_matrix_postsvd = get_Ci_overlap_matrix_postsvd(n_svd)
|
|
# !!!
|
|
aa = Ci_h_matrix_postsvd
|
|
bb = Ci_overlap_matrix_postsvd
|
|
eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, check_finite=True, homogeneous_eigvals=False)
|
|
print( eigvals_postsvd+9.194966082434476 )
|
|
# !!!
|
|
print("end code after {:.3f} minutes".format((time.time()-t0)/60.) )
|
|
# !!!
|
|
# !!!
|
|
|
|
|