Skip to content

Commit

Permalink
Update documentation and examples
Browse files Browse the repository at this point in the history
  • Loading branch information
Juha Tiihonen committed Jan 20, 2022
1 parent f03e629 commit 7e1dbd2
Show file tree
Hide file tree
Showing 7 changed files with 317 additions and 264 deletions.
231 changes: 127 additions & 104 deletions docs/README.md
View file

Large diffs are not rendered by default.

144 changes: 66 additions & 78 deletions docs/examples/benzene.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,53 @@

# Benzene: line-search example
# 2-parameter problem: CC/CH bond lengths
# DMC with DFT (PBE) as the surrogate
# Surrogate theory: DFT (PBE)
# Stochastic tehory: DMC
#
# This example executes a standard line-search minimization workflow for the
# ground state of the benzene molecule, using DFT (PBE) as the surrogate
# method and Diffusion Monte Carlo (DMC) as the stochastic method.
#
# Computing task: Suitable for institutional clusters

from numpy import mean,array,sin,pi,cos,diag,linalg
# First, the user must set up Nexus according to their computing environment.

# Parametric mappings
from nexus import generate_pwscf, generate_qmcpack, job, obj, Structure, run_project
from nexus import generate_pw2qmcpack, generate_physical_system, settings

# Modify the below variables as needed
cores = 24
presub = ''
qeapp = '/path/to/pw.x'
p2qapp = '/path/to/pw2qmcpack.x'
qmcapp = '/path/to/qmcpack'
scfjob = obj(app=qeapp, cores=cores,ppn=cores,presub=presub,hours=2)
p2qjob = obj(app=p2qapp,cores=1,ppn=1,presub=presub,minutes=5)
optjob = obj(app=qmcapp,cores=cores,ppn=cores,presub=presub,hours=12)
dmcjob = obj(app=qmcapp,cores=cores,ppn=cores,presub=presub,hours=12)
nx_settings = obj(
sleep = 3,
pseudo_dir = 'pseudos',
runs = '',
results = '',
status_only = 0,
generate_only = 0,
machine = 'ws24',
)
settings(**nx_settings) # initiate nexus

# Pseudos (execute download_pseudos.sh in the working directory)
relaxpseudos = ['C.pbe_v1.2.uspp.F.upf', 'H.pbe_v1.4.uspp.F.upf']
qmcpseudos = ['C.ccECP.xml','H.ccECP.xml']
scfpseudos = ['C.upf','H.upf']


# Implement the following parametric mappings for benzene
# p0: C-C distance
# p1: C-H distance

from numpy import mean,array,sin,pi,cos,diag,linalg

# Forward mapping: produce parameter values from an array of atomic positions
def pos_to_params(pos):
pos = pos.reshape(-1,3) # make sure of the shape
Expand Down Expand Up @@ -72,41 +111,11 @@ def params_to_pos(params):
elem = 6*['C']+6*['H']

# Check consistency of the mappings
if any(abs(pos_to_params(params_to_pos(p_init))-p_init)>1e-10):
print('Trouble with consistency!')
from surrogate_tools import check_mapping_consistency
if check_mapping_consistency(p_init,pos_to_params,params_to_pos):
print('Parameter mappings are consistent at the equilibrium!')
#end if

# Define inputs for Nexus, including job-producing functions and pseudopotentials

from nexus import generate_pwscf, generate_qmcpack, job, obj, Structure, run_project
from nexus import generate_pw2qmcpack, generate_physical_system, settings

# Pseudos
relaxpseudos = ['C.pbe_v1.2.uspp.F.UPF', 'H.pbe_v1.4.uspp.F.UPF']
qmcpseudos = ['C.ccECP.xml','H.ccECP.xml']
scfpseudos = ['C.upf','H.upf']

# Setting for the nexus jobs based on file layout and computing environment
# These settings must be changed accordingly by the user
cores = 4
presub = ''
qeapp = '/path/to/pw.x'
p2qapp = '/path/to/pw2qmcpack.x'
qmcapp = '/path/to/qmcpack'
scfjob = obj(app=qeapp, cores=cores,ppn=cores,presub=presub,hours=2)
p2qjob = obj(app=p2qapp,cores=1,ppn=1,presub=presub,minutes=5)
optjob = obj(app=qmcapp,cores=cores,ppn=cores,presub=presub,hours=12)
dmcjob = obj(app=qmcapp,cores=cores,ppn=cores,presub=presub,hours=12)
nx_settings = obj(
sleep = 3,
pseudo_dir = 'pseudos',
runs = '',
results = '',
status_only = 0,
generate_only = 0,
machine = 'ws4',
)
settings(**nx_settings) # initiate nexus

# The following data structures provide simulation inputs for each theory

Expand Down Expand Up @@ -188,7 +197,6 @@ def params_to_pos(params):
)

# construct Nexus system based on position
valences = obj(C=4,H=1) # pseudo-valences

def get_system(pos):
structure = Structure(
Expand All @@ -200,7 +208,7 @@ def get_system(pos):
kgrid = (1,1,1),
kshift = (0,0,0),
)
return generate_physical_system(structure=structure,**valences)
return generate_physical_system(structure=structure,C=4,H=1)
#end def

# return a 1-item list of Nexus jobs: SCF relaxation
Expand Down Expand Up @@ -240,27 +248,28 @@ def get_relax_job(pos,path,**kwargs):
''')
def get_phonon_jobs(pos,path,**kwargs):
scf = generate_pwscf(
system = get_system(pos),
job = job(**scfjob),
path = path,
**scf_pes_inputs
)
scf.input.control.disk_io = 'low' # write orbitals
phonon = GenericSimulation(
system = get_system(pos),
job = job(**phjob),
job = job(**scfjob),
path = path,
input = ph_input,
identifier = 'phonon',
**scf_pes_inputs,
)
scf.input.control.disk_io = 'low' # write charge density
phonon = GenericSimulation(
system = get_system(pos),
job = job(**phjob),
path = path,
input = ph_input,
identifier = 'phonon',
dependencies = (scf,'other'),
)
q2r = GenericSimulation(
system = get_system(pos),
job = job(**q2rjob),
path = path,
input = q2r_input,
identifier = 'q2r',
dependencies = (phonon,'other'),
)
# nexus automatically executes the jobs subsequently
return [scf,phonon,q2r]
#end def

Expand Down Expand Up @@ -366,7 +375,7 @@ def get_dmc_jobs(pos,path,sigma,jastrow=None,**kwargs):
from surrogate_tools import load_gamma_k,compute_hessian,print_hessian_delta

phonon_dir = 'phonon'
run_project(get_phonon_jobs(pos_relax,path=phonon_dir)) # run phonon jobs with Nexus
run_project(get_phonon_jobs(pos_relax,path=phonon_dir)) # next, run phonon job

# load the real-space Hessian from the phonon calculation output
hessian_real = load_gamma_k(phonon_dir+'/FC.fc',len(elem))
Expand Down Expand Up @@ -451,14 +460,14 @@ def get_dmc_jobs(pos,path,sigma,jastrow=None,**kwargs):
scan_data = data_load
#end if


# print output
error_scan_diagnostics(scan_data,steps_times_error2)


# 4-5) Stochastic: Line-search

from surrogate_relax import surrogate_diagnostics,average_params
n_max = 3 # number of iterations
from surrogate_relax import surrogate_diagnostics,average_params,run_linesearch

# store common line-search settings
ls_settings = obj(
Expand All @@ -485,32 +494,10 @@ def get_dmc_jobs(pos,path,sigma,jastrow=None,**kwargs):
**ls_settings,
)

# for convenience, try loading first, then execute
data_load = data.load_from_file()
if data_load is None:
data.shift_positions()
run_project(data.get_job_list())
data.load_results()
data.write_to_file()
else:
data = data_load
#end if
data_ls = [data] # starts a list of line-search objects

# repeat to n_max iterations
for n in range(1,n_max):
data = data.iterate(ls_settings=ls_settings)
data_load = data.load_from_file()
if data_load is None:
data.shift_positions()
run_project(data.get_job_list())
data.load_results()
data.write_to_file()
else:
data = data_load
#end if
data_ls.append(data)
#end for
# Run line-search up to n_max iterations, return list of line-search objects
# run_linesearch function implements standard steps of the stochastic
# line-search, including saving and loading restart files
data_ls = run_linesearch(data,n_max=3,ls_settings=ls_settings)

params_final, p_errs_final = average_params(
data_ls, # input list of line-searches
Expand All @@ -523,5 +510,6 @@ def get_dmc_jobs(pos,path,sigma,jastrow=None,**kwargs):
print('{} +/- {}'.format(p,err))
#end for

# print and plot standard output of the iteration
surrogate_diagnostics(data_ls)
plt.show()
Loading

0 comments on commit 7e1dbd2

Please sign in to comment.