Manual submission for NPT ensemble with vc relax

[chenyubi14]

Dear developers,

I am currently exploring NPT relaxation with ab initio methods. In the provided tutorial1, relax.vc_relax(target_press = 0) is used for this purpose, and it will run many configurations one job at a time.

## minimize by sscha.Relax.SSCHA
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
relax = sscha.Relax.SSCHA(minimizer, calculators,
            N_configs=N_CONFIGS,
            max_pop=MAX_ITERATIONS,
            save_ensemble=True,
            cluster=cluster_config)
relax.vc_relax(target_press = 0)

However, my preference leans towards manual submission because relax.vc_relax executes sequentially, which is too slow outside of classical potentials. Instead, I prefer to submit multiple ab initio jobs simultaneously, and process by SSCHA_Minimizer(sscha.SchaMinimizer.SSCHA_Minimizer) rather than relax(sscha.Relax.SSCHA)

### my preferred manual submission procedure for NVT

## ensemble generation
ensemble.generate(number_of_configurations)
ensemble.save(directory, population_id)

## run ab initio jobs by manual submission

## load ensemble
ensemble.load(directory, population = population_id, N = number_of_configurations)

## minimize by sscha.SchaMinimizer.SSCHA_Minimizer
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minim.init()
minim.run()

Unfortunately, the manual process is tailored for NVT relaxation. I am encountering difficulties in configuring an equivalent setup for NPT relaxation. In minim.init() or minim.run(), there is no argument like target_press for NPT.
From my current understanding of the code, relax.vc_relax is the only method for NPT, and it does not fit for the manual submission.

How can I use NPT for manual submission? Specifically, could you provide me with an example code to generate and load ensemble for NPT?

Thanks a lot,
Yubi

[mesonepigreco]

Dear Yubi,

There is actually an example in the github repository to do exactly what you want in ``Examples/SimpleVcRelax/vc_relax_no_automatic_submission.py''.

Here is the code:

import cellconstructor as CC, cellconstructor.Phonons
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Utilities
import sscha.Relax


# Load the dynamical matrix
dyn = CC.Phonons.Phonons("../ensemble_data_test/dyn", nqirr=1)

# Load the ensemble (0 K)
ensemble = sscha.Ensemble.Ensemble(dyn, T0=0)
ensemble.load("../ensemble_data_test/", population=2, N=1000)

# Prepare the minimization
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)

# Prepare the auxiliary functions to save the minimization output
ioinfo = sscha.Utilities.IOInfo()
ioinfo.SetupSaving("minim_1")

# Prepare the vc-relax
relax = sscha.Relax.SSCHA(minim, max_pop=1)
relax.setup_custom_functions(custom_function_post=ioinfo.CFP_SaveAll)

# Run the vc-relax
relax.vc_relax(target_press=0, static_bulk_modulus=100, 
               restart_from_ens=True)

# Save the final dynamical matrix
relax.minim.dyn.save_qe("final_dyn")

This is the code to perform 1 single step of the vc-relax with an ensemble already generated and computed. Indeed the steps to generate the ensemble are the same as the standard manual submission as in the code you cited. Adapting it here at the bottom of this script:

# Generate a novel ensemble from the final dynamical matrix after the vc-relax iteration
ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, relax.minim.ensemble.current_T)

# Generate the ensemble for the next iteration
ensemble.generate(number_of_configurations)
ensemble.save(directory, next_population_id)  # + 1 with respect to the current one

The rest should be the same as usual.
Remember when doing the manual submission that the units in which the forces, energies and stresses are loaded in the SSCHA coincides with quantum espresso (Rydberg atomic units):
energies in Ry, forces in Ry/Bohr, stresses in Ry/Bohr^3

Also the stress is - dE/dV

Bests,
Lorenzo

[chenyubi14]

Dear Lorenzo,

Thanks a lot for the answer. I have successfully run the NPT relaxation manually, but I have one more question regarding the parameter static_bulk_modulus in relax.vc_relax(static_bulk_modulus=100). The output of vc_relax contains a new bulk modulus, with a comment to not use it.

 Currently estimated bulk modulus =   69.500 GPa
 (Note: this is just indicative, do not use it for computing bulk modulus)

However, you have reminded me previously that

It is fundamental you choose a good value of the bulk modulus, as the algorithm will not remember previous iterations to auto-adjust it

To set static_bulk_modulus, I used ab initio calculations to obtain the bulk modulus by numerically fitting the Birch-Murnaghan equation of state. Should I continue with the same bulk modulus value, or should I consider using the Currently estimated bulk modulus? If both answers are wrong, how can I obtain the auto-adjusted bulk modulus?