#!/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 "%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()