mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-06-01 10:55:18 +02:00
230 lines
7.1 KiB
Python
Executable File
230 lines
7.1 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
|
|
import scipy as sp
|
|
import scipy.optimize
|
|
import numpy as np
|
|
import sys
|
|
import os
|
|
import time
|
|
import subprocess
|
|
from math import sqrt
|
|
#
|
|
#from bayes_opt import BayesianOptimization
|
|
|
|
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
|
|
# this should give: /home/ammar/qmcchem
|
|
sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
|
|
# this should add: /home/ammar/qmcchem/EZFIO/Python/
|
|
|
|
from ezfio import ezfio
|
|
|
|
# PARAMETERS
|
|
thresh = 1.e-2
|
|
block_time = 20
|
|
|
|
def main():
|
|
if len(sys.argv) != 2:
|
|
print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
|
|
sys.exit(1)
|
|
|
|
filename = sys.argv[1]
|
|
ezfio.set_file(filename)
|
|
|
|
|
|
def make_atom_map():
|
|
labels = {}
|
|
dimension = 0
|
|
for i,k in enumerate(ezfio.nuclei_nucl_label):
|
|
if k in labels:
|
|
labels[k].append(i)
|
|
else:
|
|
labels[k] = [dimension, i]
|
|
dimension += 1
|
|
atom_map = [[] for i in range(dimension)]
|
|
for atom in labels.keys():
|
|
l = labels[atom]
|
|
atom_map[l[0]] = l[1:]
|
|
return atom_map
|
|
atom_map = make_atom_map()
|
|
|
|
|
|
def get_params_pen():
|
|
d = ezfio.jastrow_jast_pen
|
|
return np.array( [d[m[0]] for m in atom_map])
|
|
|
|
|
|
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 get_energy_deriv_nucPar():
|
|
# !!!
|
|
buffer = subprocess.check_output(['qmcchem', 'result', filename], encoding='UTF-8')
|
|
# !!!
|
|
if buffer.strip()!= "":
|
|
e_der1 = [ buffer.splitlines()[13].split()[2], buffer.splitlines()[14].split()[2], buffer.splitlines()[15].split()[2] ]
|
|
err_der1 = [ buffer.splitlines()[13].split()[4], buffer.splitlines()[14].split()[4], buffer.splitlines()[15].split()[4] ]
|
|
e_der2 = [ buffer.splitlines()[2].split()[2], buffer.splitlines()[3].split()[2], buffer.splitlines()[4].split()[2] ]
|
|
err_der2 = [ buffer.splitlines()[2].split()[4], buffer.splitlines()[3].split()[4], buffer.splitlines()[4].split()[4] ]
|
|
return e_der1, e_der2
|
|
else:
|
|
return None
|
|
#####################################################################################################
|
|
|
|
def get_variance():
|
|
buffer = subprocess.check_output(['qmcchem', 'result', '-e',
|
|
'e_loc_qmcvar', filename],
|
|
encoding='UTF-8')
|
|
if buffer.strip() != "":
|
|
buffer = buffer.splitlines()[-1]
|
|
_, variance, error = [float(x) for x in buffer.split()]
|
|
return variance, error
|
|
else:
|
|
return None, None
|
|
|
|
|
|
def set_params_pen(x):
|
|
x = np.abs(x)
|
|
y=list(ezfio.jastrow_jast_pen)
|
|
for i,m in enumerate(atom_map):
|
|
for j in m:
|
|
y[j] = x[i]
|
|
ezfio.set_jastrow_jast_pen(y)
|
|
|
|
|
|
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(block_time),
|
|
# '--time-step=0.3',
|
|
# '--stop-time=36000',
|
|
# '--norm=1.e-5',
|
|
# '-w', '10',
|
|
# filename])
|
|
subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
|
|
'-l', str(block_time), filename])
|
|
# !!!
|
|
# !!!
|
|
#####################################################################################################
|
|
# Only for CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr.
|
|
# !!!
|
|
def fprime(x):
|
|
# !!!
|
|
print ("derivative on: x = %s"%str(x))
|
|
# !!!
|
|
e_der1, e_der2 = get_energy_deriv_nucPar()
|
|
e, _ = get_energy()
|
|
energy_deriv_0 = 2. * ( float(e_der1[0]) - e * float(e_der2[0]) )
|
|
energy_deriv_1 = 2. * ( float(e_der1[1]) - e * float(e_der2[1]) )
|
|
energy_deriv = np.array([ energy_deriv_0 , energy_deriv_1 ])
|
|
# !!!
|
|
print(energy_deriv)
|
|
# !!!
|
|
return energy_deriv
|
|
# !!!
|
|
#####################################################################################################
|
|
# !!!
|
|
memo_energy = {'fmin': 100000000.}
|
|
def f(x):
|
|
print ("x = %s"%str(x))
|
|
# !!!
|
|
h = str(x)
|
|
if h in memo_energy:
|
|
return memo_energy[h]
|
|
# !!!
|
|
set_params_pen(x)
|
|
set_vmc_params()
|
|
# !!!
|
|
pid = os.fork()
|
|
if pid == 0:
|
|
# In child process
|
|
run_qmc()
|
|
# Exit with status os.EX_OK using os._exit() method. The value of os.EX_OK is 0
|
|
os._exit(os.EX_OK)
|
|
else:
|
|
# In parent process
|
|
import atexit
|
|
atexit.register(stop_qmc)
|
|
# !!!
|
|
time.sleep(3.*block_time/4.)
|
|
# !!!
|
|
err = thresh+1.
|
|
local_thresh = thresh
|
|
while err > local_thresh:
|
|
time.sleep(block_time)
|
|
e, e_err = get_energy()
|
|
variance, v_err = get_variance()
|
|
if e is None or variance is None:
|
|
continue
|
|
# !!!
|
|
#energy = e + variance
|
|
#err = sqrt(e_err*e_err+v_err*v_err)
|
|
energy = e
|
|
err = e_err
|
|
# !!!
|
|
#print(" %f %f %f %f %f %f"%(e, e_err, variance, v_err, energy, err))
|
|
# !!!
|
|
print(" %f %f"%(energy, err))
|
|
if (energy-2.*err) > memo_energy['fmin']+thresh:
|
|
local_thresh = 10.*thresh
|
|
elif (energy+2.*err) < memo_energy['fmin']-thresh:
|
|
local_thresh = 10.*thresh
|
|
# !!!
|
|
# Check if PID is still running
|
|
try:
|
|
os.kill(pid,0)
|
|
except OSError:
|
|
print("---")
|
|
break
|
|
stop_qmc()
|
|
os.wait()
|
|
# !!!
|
|
#memo_energy[h] = energy + err
|
|
memo_energy[h] = energy
|
|
# !!!
|
|
memo_energy['fmin'] = min(energy, memo_energy['fmin'])
|
|
return energy
|
|
# !!!
|
|
def run():
|
|
x = get_params_pen()
|
|
# !!!
|
|
if sum(x) == 0.:
|
|
jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
|
|
x += jast_a_up_dn
|
|
# !!!
|
|
opt = sp.optimize.minimize(f,x,method="Powell", options= {'disp':True, 'ftol':thresh,'xtol':0.02})
|
|
# !!!
|
|
bnds = ((0.001, 100), (0.001, 100))
|
|
#opt = sp.optimize.minimize(f, x, method='CG', jac=fprime, options={'gtol': thresh, 'disp': True})
|
|
#opt = sp.optimize.minimize(f, x, method='TNC', jac=fprime, options={'gtol': thresh, 'disp': True}, bounds=bnds)
|
|
#opt = sp.optimize.minimize(f, x, method='Newton-CG', jac=fprime)
|
|
#opt = sp.optimize.minimize(f, x, method='trust-constr', jac=fprime)
|
|
print("x = "+str(opt))
|
|
set_params_pen(opt['x'])
|
|
# !!!
|
|
# !!!
|
|
run()
|
|
# !!!
|
|
# !!!
|
|
|
|
if __name__ == '__main__':
|
|
main()
|
|
|
|
|
|
|