2. Compute methane loading for one MOF¶
With the codes set up and the daemon running, we are ready to do our first calculation using AiiDA – the methane loading of a MOF at 65 bar.
We will use the RASPA code to perform a grand-canonical Monte Carlo (GCMC) calculation, trying to insert methane molecules into the nanoporous framework at the given pressure.
In principle, we could continue to use the verdi shell
or
jupyter notebooks
but in order to speed things up, we’ve already
prepared a template python script that you can
download here
.
In the following, you will adapt it to your needs.
2.1. RASPA input parameters¶
The aiida-raspa plugin
uses ParameterData
nodes to specify the input parameters of RASPA:
parameters = ParameterData(dict={
"GeneralSettings": {
"SimulationType" : "MonteCarlo",
"NumberOfCycles" : 1000,
"NumberOfInitializationCycles" : 1000,
"PrintEvery" : 100,
"CutOff" : 12.0, # (Angstroms)
"Forcefield" : "UFF-TraPPE",
"ChargeMethod" : "None",
"UnitCells" : "<int> <int> <int>",
"ExternalTemperature" : <float (K)>,
"ExternalPressure" : <float (Pa)>,
},
"Component": [{
"MoleculeName" : "methane",
"MoleculeDefinition" : "TraPPE",
"MolFraction" : 1.0,
"TranslationProbability" : 1.0,
"RotationProbability" : 1.0,
"ReinsertionProbability" : 1.0,
"SwapProbability" : 1.0,
"CreateNumberOfMolecules" : 0,
}],
})
2.1.1. Exercise¶
The ParameterData
dictionary is missing a number of parameters for
which you’ll need to figure out reasonable values.
UnitCells
Our simulations are performed under periodic boundary conditions. This means, we need to make our simulation cell large enough that a molecule will never interact with two periodic copies of any of its neighbors. Given the cutoff radius of 12 Angstroms, how often do you need to replicate the unit cell of the material?Hint: The CIF files include information on the size of the unit cell.
For
ExternalTemperature
/ExternalPressure
Use ambient temperature and standard methane desorption pressure for natural gas vehicles.Currently, the probabilies for translation, rotation, reinsertion and swap Monte Carlo moves are all equal, which is suboptimal. How would optimize them?
Another input of the calculation is the atomic structure of the MOF we
want to use. Find the structure labelled “ABUWOJ” (hint: filter by the
label
) and note down its PK or UUID.
Note
Once you know the PK or UUID of a node, you can always load it using
structure = load_node(<pk>) # load using PK (specific to your database)
structure = load_node('<uuid>') # load using UUID (same for everyone)
2.2. Creating the calculation¶
Every calculation sent to a cluster is linked to a code, which describes
the executable file to be used. We need to load the raspa@bazis
code
that we set up before:
from aiida.common.example_helpers import test_and_get_code
raspa_code = test_and_get_code("raspa@bazis", expected_code_type='raspa')
Now we’ll specify a few pieces of information to pass on to the slurm scheduler that manages calculations on the cluster, such as how many compute nodes to use or the maximum time allowed for the calculation:
options = {
"resources": {
"num_machines": 1, # run on 1 node
"tot_num_mpiprocs": 1, # use 1 process
"num_mpiprocs_per_machine": 1,
},
"max_wallclock_seconds": 1 * 60 * 60, # 1h walltime
"max_memory_kb": 2000000, # 2GB memory
"queue_name": "molsim", # slurm partition to use
"withmpi": False, # we run in serial mode
}
| **Note**
| AiiDA supports different types of schedulers via plugins,
including slurm, pbspro and sge.
2.3. Submitting the calculation¶
In principle, you can now simply submit the calculation from the
verdi shell
using
RaspaCalculation = CalculationFactory('raspa')
submit(RaspaCalculation.process(),
code=raspa_code,
structure=structure,
parameters=parameters,
_options=options,
)
When the submission script is ready, submit it to the AiiDA daemon:
$ verdi run raspa_loading.py
Note
By default, the daemon polls for new calculations every 30
seconds, i.e. you may need to wait up to 30 seconds before your
calculation starts running.
If your calculation is stuck in the NEW
state, it may indicate
that your daemon is not running or AiiDA is unable to find the
calculation plugin. Try:
reentry scan # make sure the raspa plugin is available to aiida
verdi daemon restart
Use verdi calculation list -a
to monitor the state of the
calculation. Once running, the calculation should finish within 5
minutes.
2.4. Analyzing your results¶
Once the calculation shows up as FINISHED
, have a look at the
result.
Raspa produces two types of output: Outputs related to the whole system (e.g. total energy, temperature) and outputs related to the component – in our case, there is only one component: methane.
Find the component_0
output of the calculation using
verdi calculation show <PK>
and use the PK of component_0
in
verdi data parameter show <PK>
to extract the average and standard
deviation of the absolute methane loading.
2.4.1. Exercise¶
What is the relative standard deviation of the loading? How could you decrease it?
Re-run the calculation with adapted settings in order to decrease the relative standard deviation below 5%