From a22b61df4b813b8adea0a0c16463b4994ce18eef Mon Sep 17 00:00:00 2001 From: danieleongari Date: Wed, 19 Feb 2020 19:29:24 +0100 Subject: [PATCH 1/3] first working prototype, with quick example --- aiida_lsmo/calcfunctions/ff_builder_module.py | 13 + aiida_lsmo/calcfunctions/ff_data.yaml | 34 + aiida_lsmo/workchains/zeoisotherm.py | 547 ++++++++++++++++ examples/data/zeo_LTA_p91.cif | 598 ++++++++++++++++++ examples/run_zeoisot.py | 75 +++ ...est_ffbuilder.py => test_ffbuilder_uff.py} | 0 examples/test_ffbuilder_zeolite.py | 30 + 7 files changed, 1297 insertions(+) create mode 100644 aiida_lsmo/workchains/zeoisotherm.py create mode 100644 examples/data/zeo_LTA_p91.cif create mode 100644 examples/run_zeoisot.py rename examples/{test_ffbuilder.py => test_ffbuilder_uff.py} (100%) create mode 100644 examples/test_ffbuilder_zeolite.py diff --git a/aiida_lsmo/calcfunctions/ff_builder_module.py b/aiida_lsmo/calcfunctions/ff_builder_module.py index 45ff7938..a595e3ea 100644 --- a/aiida_lsmo/calcfunctions/ff_builder_module.py +++ b/aiida_lsmo/calcfunctions/ff_builder_module.py @@ -42,6 +42,10 @@ def render_ff_mixing_def(ff_data, params): # TODO: this needs to be sorted for python versions where dictionaries are not sorted! #pylint: disable=fixme # If separate_interactions==True, prints only "zero-potential" interactions for the molecules for atom_type, ff_pot in ff_data['framework'][params['ff_framework']]['atom_types'].items(): + try: #correct if the atom_type has charges (e.g., zeolites) + ff_pot = ff_pot['force_field'] + except TypeError: + pass force_field_lines.append(" ".join([str(x) for x in [atom_type] + ff_pot])) for molecule, ff_name in params['ff_molecules'].items(): for atom_type, val in ff_data[molecule][ff_name]['atom_types'].items(): @@ -126,6 +130,15 @@ def render_pseudo_atoms_def(ff_data, params): output.append("# number of pseudo atoms") pseudo_atoms_lines = [] + + #add framework charges (e.g., zeolites) + for atom_type, val in ff_data['framework'][params['ff_framework']]['atom_types'].items(): + try: + type_settings = val['pseudo_atom'] + pseudo_atoms_lines.append(" ".join([str(x) for x in [atom_type] + type_settings])) + except TypeError: + pass + for molecule, ff_name in params['ff_molecules'].items(): for atom_type, val in ff_data[molecule][ff_name]['atom_types'].items(): type_settings = val['pseudo_atom'] diff --git a/aiida_lsmo/calcfunctions/ff_data.yaml b/aiida_lsmo/calcfunctions/ff_data.yaml index d4a6e257..bef7521e 100644 --- a/aiida_lsmo/calcfunctions/ff_data.yaml +++ b/aiida_lsmo/calcfunctions/ff_data.yaml @@ -228,6 +228,24 @@ framework: Zn_: [lennard-jones, 27.67992, 4.04468] Zr_: [lennard-jones, 34.7257, 2.7832] + Zeo-test: + description: >- + A zeolite FF for testing purpose. + Use Raspa to identify O and Oa: pseudo_atom/radii and connectivity are important + atom_types: + O: + force_field: [lennard-jones, 50.0, 3.0] + pseudo_atom: [yes, O, O, 0, 15.9994, -1.025, 0.0, 1.0, 0.5, 2, 0, absolute, 0] + Oa: + force_field: [lennard-jones, 50.0, 3.0] + pseudo_atom: [yes, O, O, 0, 15.9994, -1.2, 0.0, 1.0, 0.5, 2, 0, absolute, 0] + Si: + force_field: [lennard-jones, 32.0, 2.9] + pseudo_atom: [yes, Si, Si, 0, 28.08, +2.05, 0.0, 1.0, 1.18, 4, 0, absolute, 0] + Al: + force_field: [lennard-jones, 24.0, 3.1] + pseudo_atom: [yes, Al, Al, 0, 26.98, +1.75, 0.0, 1.0, 1.18, 4, 0, absolute, 0] + # Molecules parameters: ff_paremeters[molecule][force_field] # pseudo_atom list: print, as, chem, oxidation, mass, charge, polarization, B-factor, radii, connectivity, anisotropic, anisotropic-type, tinker-type # Except when stated in the description, critical constants (Tc(K), Pc(Pa), acentric factor) come from @@ -960,3 +978,19 @@ Xe: pseudo_atom: [yes, Xe, Xe, 0, 131.29, 0.0, 0.0, 1.0, 1.0, 0, 0, relative, 0] atomic_positions: - [Xe_gas 0, 0, 0] + +# Cations for zeolite +Na: + critical_constants: + tc: 999.99 + pc: 9999999 + af: 0.0 + Zeo-test: + description: + A Na+ testing purpose force field for zeolites. + atom_types: + Na_cat: + force_field: [lennard-jones, 5.0, 2.5] + pseudo_atom: [yes, Na, Na, 0, 22.98, +1.0, 0.0, 1.0, 1.0, 0, 0, relative, 0] + atomic_positions: + - [Na_cat 0, 0, 0] diff --git a/aiida_lsmo/workchains/zeoisotherm.py b/aiida_lsmo/workchains/zeoisotherm.py new file mode 100644 index 00000000..96dab6ea --- /dev/null +++ b/aiida_lsmo/workchains/zeoisotherm.py @@ -0,0 +1,547 @@ +# -*- coding: utf-8 -*- +"""ZeoIsotherm work chain""" +from __future__ import absolute_import + +import os +from six.moves import range + +from aiida.plugins import CalculationFactory, DataFactory, WorkflowFactory +from aiida.orm import Dict, Float, Str, List, SinglefileData +from aiida.engine import calcfunction +from aiida.engine import WorkChain, ToContext, append_, while_, if_ +from aiida_lsmo.utils import aiida_dict_merge +from aiida_lsmo.utils import dict_merge + +# import sub-workchains +RaspaBaseWorkChain = WorkflowFactory('raspa.base') # pylint: disable=invalid-name + +# import calculations +ZeoppCalculation = CalculationFactory('zeopp.network') # pylint: disable=invalid-name +FFBuilder = CalculationFactory('lsmo.ff_builder') # pylint: disable=invalid-name + +# import aiida data +CifData = DataFactory('cif') # pylint: disable=invalid-name +ZeoppParameters = DataFactory('zeopp.parameters') # pylint: disable=invalid-name + +# Deafault parameters +ISOTHERMPARAMETERS_DEFAULT = { #TODO: create IsothermParameters instead of Dict # pylint: disable=fixme + "ff_framework": "Zeo-test", # (str) Forcefield of the structure. + "ff_separate_interactions": False, # (bool) Use "separate_interactions" in the FF builder. + "ff_mixing_rule": "Lorentz-Berthelot", # (string) Choose 'Lorentz-Berthelot' or 'Jorgensen'. + "ff_tail_corrections": True, # (bool) Apply tail corrections. + "ff_shifted": False, # (bool) Shift or truncate the potential at cutoff. + "ff_cutoff": 12.0, # (float) CutOff truncation for the VdW interactions (Angstrom). + "temperature": 300, # (float) Temperature of the simulation. + "temperature_list": None, # (list) To be used by IsothermMultiTempWorkChain. + "zeopp_volpo_samples": int(1e5), # (int) Number of samples for VOLPO calculation (per UC volume). + "zeopp_block_samples": int(100), # (int) Number of samples for BLOCK calculation (per A^3). + "raspa_minKh": 1e-10, # (float) If Henry coefficient < raspa_minKh do not run the isotherm (mol/kg/Pa). + "raspa_verbosity": 10, # (int) Print stats every: number of cycles / raspa_verbosity. + "raspa_widom_cycles": int(1e5), # (int) Number of Widom cycles. + "raspa_nvt_cations_cycles": int(1e4), # (int) Number of GCMC production cycles. + "raspa_gcmc_init_cycles": int(1e3), # (int) Number of GCMC initialization cycles. + "raspa_gcmc_prod_cycles": int(1e4), # (int) Number of GCMC production cycles. + "pressure_list": None, # (list) Pressure list for the isotherm (bar): if given it will skip to guess it. + "pressure_precision": 0.1, # (float) Precision in the sampling of the isotherm: 0.1 ok, 0.05 for high resolution. + "pressure_maxstep": 5, # (float) Max distance between pressure points (bar). + "pressure_min": 0.001, # (float) Lower pressure to sample (bar). + "pressure_max": 10 # (float) Upper pressure to sample (bar). +} + + +# calcfunctions (in order of appearence) +@calcfunction +def get_molecule_dict(molecule_name): + """Get a Dict from the isotherm_molecules.yaml""" + import ruamel.yaml as yaml + thisdir = os.path.dirname(os.path.abspath(__file__)) + yamlfile = os.path.join(thisdir, "isotherm_data", "isotherm_molecules.yaml") + with open(yamlfile, 'r') as stream: + yaml_dict = yaml.safe_load(stream) + molecule_dict = yaml_dict[molecule_name.value] + return Dict(dict=molecule_dict) + + +@calcfunction +def get_zeopp_parameters(molecule_dict, isotparam): + """Get the ZeoppParameters from the inputs of the workchain""" + probe_rad = molecule_dict["proberad"] + param_dict = { + 'ha': 'DEF', + 'volpo': [probe_rad, probe_rad, isotparam['zeopp_volpo_samples']], + 'block': [probe_rad, isotparam['zeopp_block_samples']], + } + return ZeoppParameters(dict=param_dict) + + +@calcfunction +def get_ff_parameters(molecule_dict, cation_name, isotparam): + """Get the parameters for ff_builder. + NOTE: for the cation consider a ff with the same name as the ff_framework. + """ + ff_params = {} + ff_params['ff_framework'] = isotparam['ff_framework'] + ff_params['ff_molecules'] = { + molecule_dict['name']: molecule_dict['forcefield'], + cation_name.value: isotparam['ff_framework'] + } + ff_params['shifted'] = isotparam['ff_shifted'] + ff_params['tail_corrections'] = isotparam['ff_tail_corrections'] + ff_params['mixing_rule'] = isotparam['ff_mixing_rule'] + ff_params['separate_interactions'] = isotparam['ff_separate_interactions'] + return Dict(dict=ff_params) + + +@calcfunction +def choose_pressure_points(inp_param, geom, raspa_widom_out): + """If 'presure_list' is not provide, model the isotherm as single-site langmuir and return the most important + pressure points to evaluate for an isotherm, in a List. + """ + if inp_param["pressure_list"]: + pressure_points = inp_param["pressure_list"] + else: + khenry = list(raspa_widom_out["framework_1"]["components"].values())[1]['henry_coefficient_average'] #mol/kg/Pa + b_value = khenry / geom['Estimated_saturation_loading'] * 1e5 #(1/bar) + pressure_points = [inp_param['pressure_min']] + while True: + pold = pressure_points[-1] + delta_p = min(inp_param['pressure_maxstep'], + inp_param['pressure_precision'] * (b_value * pold**2 + 2 * pold + 1 / b_value)) + pnew = pold + delta_p + if pnew <= inp_param['pressure_max']: + pressure_points.append(pnew) + else: + pressure_points.append(inp_param['pressure_max']) + break + return List(list=pressure_points) + + +@calcfunction +def get_geometric_dict(zeopp_out, molecule): + """Return the geometric Dict from Zeopp results, including Qsat and is_porous""" + geometric_dict = zeopp_out.get_dict() + geometric_dict.update({ + 'Estimated_saturation_loading': zeopp_out['POAV_cm^3/g'] * molecule['molsatdens'], + 'Estimated_saturation_loading_unit': 'mol/kg', + 'is_porous': geometric_dict["POAV_A^3"] > 0.000 + }) + return Dict(dict=geometric_dict) + + +@calcfunction +def get_output_parameters(geom_out, inp_params, widom_out=None, pressures=None, **gcmc_out_dict): + """Merge results from all the steps of the work chain.""" + + out_dict = geom_out.get_dict() + + if out_dict['is_porous'] and widom_out: + widom_out_mol = list(widom_out["framework_1"]["components"].values())[1] + + out_dict.update({ + 'temperature': inp_params['temperature'], + 'temperature_unit': 'K', + 'is_kh_enough': widom_out_mol['henry_coefficient_average'] > inp_params['raspa_minKh'] + }) + + widom_labels = [ + 'henry_coefficient_average', + 'henry_coefficient_dev', + 'henry_coefficient_unit', + 'adsorption_energy_widom_average', + 'adsorption_energy_widom_dev', + 'adsorption_energy_widom_unit', + ] + + for label in widom_labels: + out_dict.update({label: widom_out_mol[label]}) + + if out_dict['is_kh_enough']: + + isotherm = { + 'pressure': pressures, + 'pressure_unit': 'bar', + 'loading_absolute_average': [], + 'loading_absolute_dev': [], + 'loading_absolute_unit': 'mol/kg', + 'enthalpy_of_adsorption_average': [], + 'enthalpy_of_adsorption_dev': [], + 'enthalpy_of_adsorption_unit': 'kJ/mol' + } + + conv_ener = 1.0 / 120.273 # K to kJ/mol + for i in range(len(pressures)): + gcmc_out = gcmc_out_dict['RaspaGCMC_{}'.format(i + 1)]["framework_1"] + gcmc_out_mol = list(gcmc_out["components"].values())[1] + conv_load = gcmc_out_mol["conversion_factor_molec_uc_to_mol_kg"] + + for label in ['loading_absolute_average', 'loading_absolute_dev']: + isotherm[label].append(conv_load * gcmc_out_mol[label]) + + for label in ['enthalpy_of_adsorption_average', 'enthalpy_of_adsorption_dev']: + if gcmc_out['general'][label]: + isotherm[label].append(conv_ener * gcmc_out['general'][label]) + else: # when there are no particles and Raspa return Null enthalpy + isotherm[label].append(None) + + out_dict.update({ + "isotherm": isotherm, + 'conversion_factor_molec_uc_to_cm3stp_cm3': gcmc_out_mol['conversion_factor_molec_uc_to_cm3stp_cm3'], + 'conversion_factor_molec_uc_to_mg_g': gcmc_out_mol['conversion_factor_molec_uc_to_mg_g'], + 'conversion_factor_molec_uc_to_mol_kg': gcmc_out_mol['conversion_factor_molec_uc_to_mol_kg'], + }) + + return Dict(dict=out_dict) + + +class ZeoIsothermWorkChain(WorkChain): + """Isotherm work chain designed specifically for zeolites: both all-silica and Al-substituted. + """ + + @classmethod + def define(cls, spec): + super(ZeoIsothermWorkChain, cls).define(spec) + + spec.expose_inputs(ZeoppCalculation, namespace='zeopp', include=['code', 'metadata']) + + spec.expose_inputs(RaspaBaseWorkChain, namespace='raspa_base', exclude=['raspa.structure', 'raspa.parameters']) + + spec.input('structure', valid_type=CifData, help='Adsorbent framework CIF.') + + spec.input('sial_ratio', valid_type=Float, default=Float(1), help='Get the Si/Al ratio in the zeolite.') + + spec.input("molecule", + valid_type=(Str, Dict), + help='Adsorbate molecule: settings to be read from the yaml.' + + 'Advanced: input a Dict for non-standard settings.') + + spec.input("cation", valid_type=Str, default=Str("Na"), help='Cation to be used for charged zeolites (Na/Ca)') + #TODO: need to implement Ca+2 # pylint: disable=fixme + + spec.input("parameters", + valid_type=Dict, + help='Parameters for the Isotherm workchain: will be merged with IsothermParameters_defaults.') + + spec.input("geometric", + valid_type=Dict, + required=False, + help='[Only used by IsothermMultiTempWorkChain] Already computed geometric properties') + + spec.outline( + cls.setup, + cls.run_zeopp, # computes volpo and blocks + cls.run_equilibrate_cation, + if_(cls.should_run_widom)( # run Widom only if porous + cls.run_raspa_widom, # run raspa widom calculation + if_(cls.should_run_gcmc)( # Kh is high enough + cls.init_raspa_gcmc, # initializate setting for GCMC + while_(cls.should_run_another_gcmc)( # new pressure + cls.run_raspa_gcmc, # run raspa GCMC calculation + ), + ), + ), + cls.return_output_parameters, + ) + + spec.expose_outputs(ZeoppCalculation, include=['block']) #only if porous + + spec.output( + 'output_parameters', + valid_type=Dict, + required=True, + help='Results of the single temperature wc: keys can vay depending on is_porous and is_kh_enough booleans.') + + def setup(self): + """Initialize the parameters""" + + # Get the molecule Dict from the yaml or directly as an input + if isinstance(self.inputs.molecule, Str): + self.ctx.molecule = get_molecule_dict(self.inputs.molecule) + elif isinstance(self.inputs.molecule, Dict): + self.ctx.molecule = self.inputs.molecule + + # Get the parameters Dict, merging defaults with user settings + self.ctx.parameters = aiida_dict_merge(Dict(dict=ISOTHERMPARAMETERS_DEFAULT), self.inputs.parameters) + + # Get integer temperature in context for easy reports + self.ctx.temperature = int(round(self.ctx.parameters['temperature'])) + + # Get the number of cations needed + natoms = self.inputs.structure.get_ase().get_number_of_atoms() + ntatoms = natoms / 3 + self.ctx.ncations = int(ntatoms / (self.inputs.sial_ratio + 1)) + + # Understand if IsothermMultiTempWorkChain is calling this work chain + if "geometric" in self.inputs: + self.ctx.multitemp_mode = 'run_single_temp' + elif self.ctx.parameters['temperature_list']: + self.ctx.multitemp_mode = 'run_geom_only' + else: + self.ctx.multitemp_mode = None + + def run_zeopp(self): + """Perform Zeo++ block and VOLPO calculations.""" + + # Skip zeopp calculation if the geometric properties are already provided by IsothermMultiTemp + if self.ctx.multitemp_mode == 'run_single_temp': + return None + + # create inputs: exposed are code and metadata + inputs = self.exposed_inputs(ZeoppCalculation, 'zeopp') + + # Set inputs for zeopp + dict_merge( + inputs, { + 'metadata': { + 'label': "ZeoppVolpoBlock", + 'call_link_label': 'run_zeopp_block_and_volpo', + }, + 'structure': + self.inputs.structure, + 'atomic_radii': + SinglefileData( + file=os.path.join(os.path.dirname(os.path.abspath(__file__)), "isotherm_data", "DEFAULT.rad")), + 'parameters': + get_zeopp_parameters(self.ctx.molecule, self.ctx.parameters) + }) + + running = self.submit(ZeoppCalculation, **inputs) + self.report("Running zeo++ block and volpo Calculation<{}>".format(running.id)) + return ToContext(zeopp=running) + + def _get_nvt_cation(self): + """Write Raspa input parameters from scratch, for a NVT calculation""" + + param = { + "GeneralSettings": { + "SimulationType": "MonteCarlo", + "NumberOfInitializationCycles": self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, + "NumberOfCycles": self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, + "PrintPropertiesEvery": int(1e10), + "PrintEvery": self.ctx.parameters['raspa_nvt_cations_cycles'] / self.ctx.parameters['raspa_verbosity'], + "RemoveAtomNumberCodeFromLabel": True, + "ModifyOxgensConnectedToAluminium": True, + "Forcefield": "Local", + "CutOff": self.ctx.parameters['ff_cutoff'], + "ChargeMethod": "Ewald", + "EwaldPrecision": 1e-6, + }, + "System": { + "framework_1": { + "type": "Framework", + "ExternalTemperature": self.ctx.parameters['temperature'], + } + }, + "Component": { + self.inputs.cation.value: { + "MoleculeDefinition": "Local", + "ExtraFrameworkMolecule": True, + "TranslationProbability": 1.0, + "RandomTranslationProbability": 1.0, #Raspa warns to use this instead of ReinsertionProbability + "CreateNumberOfMolecules": self.ctx.ncations, + }, + }, + } + + # Check particular conditions and settings + # DISABLED: I want to have already from the beginning a larger unit cell (so, no problems with ncations) + # mult = check_resize_unit_cell(self.inputs.structure, 2 * self.ctx.parameters['ff_cutoff']) + # param["System"]["framework_1"]["UnitCells"] = "{} {} {}".format(mult[0], mult[1], mult[2]) + return param + + def run_equilibrate_cation(self): + """Run a calculation in Raspa to equilibrate postion of cations.""" + + # Initialize the input for raspa_base, which later will need only minor updates for GCMC + self.ctx.inp = self.exposed_inputs(RaspaBaseWorkChain, 'raspa_base') + self.ctx.inp['metadata']['label'] = "RaspaWidom" + self.ctx.inp['metadata']['call_link_label'] = "run_raspa_cation" + + # Create input parameters and pass the framework + self.ctx.raspa_param = self._get_nvt_cation() + self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) + self.ctx.inp['raspa']['framework'] = {"framework_1": self.inputs.structure} + + # Generate the force field with the ff_builder + ff_params = get_ff_parameters(self.ctx.molecule, self.inputs.cation, self.ctx.parameters) + files_dict = FFBuilder(ff_params) + self.ctx.inp['raspa']['file'] = files_dict + + running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) + self.report("Running Raspa NVT @ {}K for equilibrating cations".format(self.ctx.temperature)) + + return ToContext(raspa_cation=running) + + def should_run_widom(self): + """Submit widom calculation only if there is some accessible volume, + also check the number of blocking spheres and estimate the saturation loading. + Also, stop if called by IsothermMultiTemp for geometric results only.""" + + # Get geometric properties and consider if IsothermMultiTempWorkChain is calling this workchain + if self.ctx.multitemp_mode == 'run_single_temp': + self.ctx.geom = self.inputs.geometric + return True + self.ctx.geom = get_geometric_dict(self.ctx.zeopp.outputs.output_parameters, self.ctx.molecule) + + if self.ctx.geom['is_porous']: + self.report("Found accessible pore volume: continue") + self.report("Found {} blocking spheres".format(self.ctx.geom['Number_of_blocking_spheres'])) + # Return block file only if blocking spheres are present + if self.ctx.geom['Number_of_blocking_spheres'] > 0: + self.out_many(self.exposed_outputs(self.ctx.zeopp, ZeoppCalculation)) + else: + self.report("No accessible pore volume: stop") + + return self.ctx.geom['is_porous'] and not self.ctx.multitemp_mode == 'run_geom_only' + + def _update_param_for_widom(self): + """Write Raspa input parameters from scratch, for a Widom calculation""" + + param = self.ctx.raspa_param + param["GeneralSettings"].update({ + "NumberOfInitializationCycles": 0, + "NumberOfCycles": self.ctx.parameters['raspa_widom_cycles'], + "PrintPropertiesEvery": self.ctx.parameters['raspa_widom_cycles'] / self.ctx.parameters['raspa_verbosity'], + "PrintEvery": int(1e10), + "HeliumVoidFraction": self.ctx.geom["POAV_Volume_fraction"], + }) + param["Component"][self.inputs.cation.value].update({ + "CreateNumberOfMolecules": 0, + }) + param["Component"][self.ctx.molecule['name']] = { + "MoleculeDefinition": "Local", + "ExtraFrameworkMolecule": False, + "WidomProbability": 1.0, + } + + # Check particular conditions and settings + if self.ctx.geom['Number_of_blocking_spheres'] > 0: + param["Component"][self.ctx.molecule['name']]["BlockPocketsFileName"] = "block_file" + + return param + + def run_raspa_widom(self): + """Run a Widom calculation in Raspa.""" + + # Update labels + self.ctx.inp['metadata']['label'] = "RaspaWidom" + self.ctx.inp['metadata']['call_link_label'] = "run_raspa_widom" + + # Update parameters dict, add block file, add restart + self.ctx.raspa_param = self._update_param_for_widom() + self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) + if self.ctx.geom['Number_of_blocking_spheres'] > 0 and self.ctx.multitemp_mode != 'run_single_temp': + self.ctx.inp["raspa"]["block_pocket"] = {"block_file": self.ctx.zeopp.outputs.block} + self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_cation.outputs.retrieved + + running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) + self.report("Running Raspa Widom @ {}K for the Henry coefficient".format(self.ctx.temperature)) + return ToContext(raspa_widom=running) + + def should_run_gcmc(self): + """Output the widom results and decide to compute the isotherm if kH > kHmin, as defined by the user""" + + self.ctx.is_kh_enough = list(self.ctx.raspa_widom.outputs['output_parameters']["framework_1"]["components"]. + values())[1]['henry_coefficient_average'] > self.ctx.parameters['raspa_minKh'] + + if self.ctx.is_kh_enough: + self.report("kH larger than the threshold: continue") + return True + + self.report("kHh lower than the threshold: stop") + return False + + def _update_param_for_gcmc(self): + """Update Raspa input parameter, from Widom to GCMC""" + + param = self.ctx.raspa_param + param["GeneralSettings"].update({ + "NumberOfInitializationCycles": self.ctx.parameters['raspa_gcmc_init_cycles'], + "NumberOfCycles": self.ctx.parameters['raspa_gcmc_prod_cycles'], + "PrintPropertiesEvery": int(1e6), + "PrintEvery": self.ctx.parameters['raspa_gcmc_prod_cycles'] / self.ctx.parameters['raspa_verbosity'] + }) + param["Component"][self.ctx.molecule['name']].update({ + "WidomProbability": 0.0, + "TranslationProbability": 1.0, + "ReinsertionProbability": 1.0, + "SwapProbability": 2.0, + }) + # Check particular conditions + if not self.ctx.molecule['singlebead']: + param["Component"][self.ctx.molecule['name']].update({"RotationProbability": 1.0}) + return param + + def init_raspa_gcmc(self): + """Choose the pressures we want to sample, report some details, and update settings for GCMC""" + + self.ctx.current_p_index = 0 + self.ctx.pressures = choose_pressure_points(self.ctx.parameters, self.ctx.geom, + self.ctx.raspa_widom.outputs.output_parameters) + + self.report("Computed Kh(mol/kg/Pa)={:.2e} POAV(cm3/g)={:.3f} Qsat(mol/kg)={:.2f}".format( + list(self.ctx.raspa_widom.outputs['output_parameters']["framework_1"]["components"].values())[1] + ['henry_coefficient_average'], self.ctx.geom['POAV_cm^3/g'], self.ctx.geom['Estimated_saturation_loading'])) + self.report("Now evaluating the isotherm @ {}K for {} pressure points".format( + self.ctx.temperature, len(self.ctx.pressures))) + + self.ctx.raspa_param = self._update_param_for_gcmc() + + def should_run_another_gcmc(self): + """We run another raspa calculation only if the current iteration is + smaller than the total number of pressures we want to compute. + """ + return self.ctx.current_p_index < len(self.ctx.pressures) + + def run_raspa_gcmc(self): + """Run a GCMC calculation in Raspa @ T,P. """ + + # Update labels + self.ctx.inp['metadata']['label'] = "RaspaGCMC_{}".format(self.ctx.current_p_index + 1) + self.ctx.inp['metadata']['call_link_label'] = "run_raspa_gcmc_{}".format(self.ctx.current_p_index + 1) + + # Update pressure (NOTE: need to convert from bar to Pa) + self.ctx.raspa_param["System"]["framework_1"]['ExternalPressure'] = \ + self.ctx.pressures[self.ctx.current_p_index] * 1e5 + + # Update parameters Dict + self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) + + # Update restart (if present, i.e., if current_p_index>0) + if self.ctx.current_p_index > 0: + self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_gcmc[self.ctx.current_p_index - + 1].outputs.retrieved + + # Create the calculation process, launch it and update pressure index + running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) + self.report("Running Raspa GCMC @ {}K/{:.3f}bar (pressure {} of {})".format( + self.ctx.temperature, self.ctx.pressures[self.ctx.current_p_index], self.ctx.current_p_index + 1, + len(self.ctx.pressures))) + self.ctx.current_p_index += 1 + return ToContext(raspa_gcmc=append_(running)) + + def return_output_parameters(self): + """Merge all the parameters into output_parameters, depending on is_porous and is_kh_ehough.""" + + gcmc_out_dict = {} + if self.ctx.geom['is_porous'] and not self.ctx.multitemp_mode == 'run_geom_only': + widom_out = self.ctx.raspa_widom.outputs.output_parameters + if self.ctx.is_kh_enough: + for calc in self.ctx.raspa_gcmc: + gcmc_out_dict[calc.label] = calc.outputs.output_parameters + else: + self.ctx.pressures = None + else: + widom_out = None + self.ctx.pressures = None + + self.out( + "output_parameters", + get_output_parameters(geom_out=self.ctx.geom, + inp_params=self.ctx.parameters, + widom_out=widom_out, + pressures=self.ctx.pressures, + **gcmc_out_dict)) + + if not self.ctx.multitemp_mode == 'run_geom_only': + self.report("Isotherm @ {}K computed: ouput Dict<{}>".format(self.ctx.temperature, + self.outputs['output_parameters'].pk)) diff --git a/examples/data/zeo_LTA_p91.cif b/examples/data/zeo_LTA_p91.cif new file mode 100644 index 00000000..035cb8a8 --- /dev/null +++ b/examples/data/zeo_LTA_p91.cif @@ -0,0 +1,598 @@ +data_crystal + +_cell_length_a 24.55500 +_cell_length_b 24.55500 +_cell_length_c 24.55500 +_cell_angle_alpha 90.00000 +_cell_angle_beta 90.00000 +_cell_angle_gamma 90.00000 + +_symmetry_space_group_name_Hall 'P 1' +_symmetry_space_group_name_H-M 'P 1' + +loop_ +_symmetry_equiv_pos_as_xyz + 'x,y,z' + +loop_ +_atom_site_label +_atom_site_type_symbol +_atom_site_fract_x +_atom_site_fract_y +_atom_site_fract_z +Si Si 0.00000 0.09316 0.18499 +Al Al 0.00000 0.18715 0.09042 +O O 0.00000 0.11367 0.24663 +O O 0.00000 0.14459 0.14591 +O O 0.05379 0.05865 0.17152 +Si Si 0.40684 0.00000 0.18499 +Si Si 0.31285 0.00000 0.09042 +O O 0.38633 0.00000 0.24663 +O O 0.35541 0.00000 0.14591 +O O 0.44135 0.05379 0.17152 +Si Si 0.50000 0.40684 0.18499 +Si Si 0.50000 0.31285 0.09042 +O O 0.50000 0.38633 0.24663 +O O 0.50000 0.35541 0.14591 +O O 0.44621 0.44135 0.17152 +Si Si 0.09316 0.50000 0.18499 +Si Si 0.18715 0.50000 0.09042 +O O 0.11367 0.50000 0.24663 +O O 0.14459 0.50000 0.14591 +O O 0.05865 0.44621 0.17152 +Si Si 0.00000 0.90684 0.81501 +Si Si 0.00000 0.81285 0.90958 +O O 0.00000 0.88633 0.75337 +O O 0.00000 0.85541 0.85409 +O O 0.05379 0.94135 0.82848 +Si Si 0.59316 0.00000 0.81501 +Si Si 0.68715 0.00000 0.90958 +O O 0.61367 0.00000 0.75337 +O O 0.64459 0.00000 0.85409 +O O 0.55865 0.05379 0.82848 +Si Si 0.50000 0.59316 0.81501 +Al Al 0.50000 0.68715 0.90958 +O O 0.50000 0.61367 0.75337 +O O 0.50000 0.64459 0.85409 +O O 0.44621 0.55865 0.82848 +Si Si 0.90684 0.50000 0.81501 +Al Al 0.81285 0.50000 0.90958 +O O 0.88633 0.50000 0.75337 +O O 0.85541 0.50000 0.85409 +O O 0.94135 0.44621 0.82848 +Si Si 0.18499 0.00000 0.09316 +Al Al 0.09042 0.00000 0.18715 +O O 0.24663 0.00000 0.11367 +O O 0.14591 0.00000 0.14459 +O O 0.17152 0.05379 0.05865 +Si Si 0.50000 0.18499 0.09316 +Al Al 0.50000 0.09042 0.18715 +O O 0.50000 0.24663 0.11367 +O O 0.50000 0.14591 0.14459 +O O 0.44621 0.17152 0.05865 +Si Si 0.31501 0.50000 0.09316 +Al Al 0.40958 0.50000 0.18715 +O O 0.25337 0.50000 0.11367 +O O 0.35409 0.50000 0.14459 +O O 0.32848 0.44621 0.05865 +Si Si 0.00000 0.31501 0.09316 +Al Al 0.00000 0.40958 0.18715 +O O 0.00000 0.25337 0.11367 +O O 0.00000 0.35409 0.14459 +O O 0.05379 0.32848 0.05865 +Si Si 0.18499 0.00000 0.90684 +Al Al 0.09042 0.00000 0.81285 +O O 0.24663 0.00000 0.88633 +O O 0.14591 0.00000 0.85541 +O O 0.17152 0.94621 0.94135 +Si Si 0.50000 0.18499 0.90684 +Al Al 0.50000 0.09042 0.81285 +O O 0.50000 0.24663 0.88633 +O O 0.50000 0.14591 0.85541 +O O 0.55379 0.17152 0.94135 +Si Si 0.31501 0.50000 0.90684 +Al Al 0.40958 0.50000 0.81285 +O O 0.25337 0.50000 0.88633 +O O 0.35409 0.50000 0.85541 +O O 0.32848 0.55379 0.94135 +Si Si 0.00000 0.31501 0.90684 +Al Al 0.00000 0.40958 0.81285 +O O 0.00000 0.25337 0.88633 +O O 0.00000 0.35409 0.85541 +O O 0.94621 0.32848 0.94135 +Si Si 0.09316 0.18499 0.00000 +Al Al 0.18715 0.09042 0.00000 +O O 0.11367 0.24663 0.00000 +O O 0.14459 0.14591 0.00000 +O O 0.05865 0.17152 0.05379 +Si Si 0.09316 0.31501 0.50000 +Al Al 0.18715 0.40958 0.50000 +O O 0.11367 0.25337 0.50000 +O O 0.14459 0.35409 0.50000 +O O 0.05865 0.32848 0.44621 +Si Si 0.18499 0.09316 0.50000 +Al Al 0.09042 0.18715 0.50000 +O O 0.24663 0.11367 0.50000 +O O 0.14591 0.14459 0.50000 +O O 0.17152 0.05865 0.44621 +Si Si 0.40684 0.18499 0.50000 +Al Al 0.31285 0.09042 0.50000 +O O 0.38633 0.24663 0.50000 +O O 0.35541 0.14591 0.50000 +O O 0.44135 0.17152 0.44621 +Si Si 0.31501 0.90684 0.00000 +Al Al 0.40958 0.81285 0.00000 +O O 0.25337 0.88633 0.00000 +O O 0.35409 0.85541 0.00000 +O O 0.32848 0.94135 0.94621 +Si Si 0.90684 0.81501 0.00000 +Al Al 0.81285 0.90958 0.00000 +O O 0.88633 0.75337 0.00000 +O O 0.85541 0.85409 0.00000 +O O 0.94135 0.82848 0.05379 +Si Si 0.18499 0.90684 0.50000 +Al Al 0.09042 0.81285 0.50000 +O O 0.24663 0.88633 0.50000 +O O 0.14591 0.85541 0.50000 +O O 0.17152 0.94135 0.55379 +Si Si 0.31501 0.59316 0.50000 +Al Al 0.40958 0.68715 0.50000 +O O 0.25337 0.61367 0.50000 +O O 0.35409 0.64459 0.50000 +O O 0.32848 0.55865 0.55379 +O O 0.94621 0.94135 0.82848 +O O 0.55865 0.94621 0.82848 +O O 0.55379 0.55865 0.82848 +O O 0.94135 0.55379 0.82848 +O O 0.94621 0.05865 0.17152 +O O 0.44135 0.94621 0.17152 +O O 0.55379 0.44135 0.17152 +O O 0.05865 0.55379 0.17152 +Si Si 0.81501 0.00000 0.90684 +Al Al 0.90958 0.00000 0.81285 +O O 0.75337 0.00000 0.88633 +O O 0.85409 0.00000 0.85541 +O O 0.82848 0.94621 0.94135 +Si Si 0.50000 0.81501 0.90684 +Al Al 0.50000 0.90958 0.81285 +O O 0.50000 0.75337 0.88633 +O O 0.50000 0.85409 0.85541 +O O 0.55379 0.82848 0.94135 +Si Si 0.68499 0.50000 0.90684 +Al Al 0.59042 0.50000 0.81285 +O O 0.74663 0.50000 0.88633 +O O 0.64591 0.50000 0.85541 +O O 0.67152 0.55379 0.94135 +Si Si 0.00000 0.68499 0.90684 +Al Al 0.00000 0.59042 0.81285 +O O 0.00000 0.74663 0.88633 +O O 0.00000 0.64591 0.85541 +O O 0.94621 0.67152 0.94135 +Si Si 0.81501 0.00000 0.09316 +Al Al 0.90958 0.00000 0.18715 +O O 0.75337 0.00000 0.11367 +O O 0.85409 0.00000 0.14459 +O O 0.82848 0.05379 0.05865 +Si Si 0.50000 0.81501 0.09316 +Al Al 0.50000 0.90958 0.18715 +O O 0.50000 0.75337 0.11367 +O O 0.50000 0.85409 0.14459 +O O 0.44621 0.82848 0.05865 +Si Si 0.68499 0.50000 0.09316 +Al Al 0.59042 0.50000 0.18715 +O O 0.74663 0.50000 0.11367 +O O 0.64591 0.50000 0.14459 +O O 0.67152 0.44621 0.05865 +Si Si 0.00000 0.68499 0.09316 +Al Al 0.00000 0.59042 0.18715 +O O 0.00000 0.74663 0.11367 +O O 0.00000 0.64591 0.14459 +O O 0.05379 0.67152 0.05865 +O O 0.94135 0.82848 0.94621 +Si Si 0.90684 0.68499 0.50000 +Al Al 0.81285 0.59042 0.50000 +O O 0.88633 0.74663 0.50000 +O O 0.85541 0.64591 0.50000 +O O 0.94135 0.67152 0.55379 +Si Si 0.81501 0.90684 0.50000 +Al Al 0.90958 0.81285 0.50000 +O O 0.75337 0.88633 0.50000 +O O 0.85409 0.85541 0.50000 +O O 0.82848 0.94135 0.55379 +Si Si 0.59316 0.81501 0.50000 +Al Al 0.68715 0.90958 0.50000 +O O 0.61367 0.75337 0.50000 +O O 0.64459 0.85409 0.50000 +O O 0.55865 0.82848 0.55379 +Si Si 0.68499 0.09316 0.00000 +Al Al 0.59042 0.18715 0.00000 +O O 0.74663 0.11367 0.00000 +O O 0.64591 0.14459 0.00000 +O O 0.67152 0.05865 0.05379 +O O 0.05865 0.17152 0.94621 +Si Si 0.81501 0.09316 0.50000 +Al Al 0.90958 0.18715 0.50000 +O O 0.75337 0.11367 0.50000 +O O 0.85409 0.14459 0.50000 +O O 0.82848 0.05865 0.44621 +Si Si 0.68499 0.40684 0.50000 +Al Al 0.59042 0.31285 0.50000 +O O 0.74663 0.38633 0.50000 +O O 0.64591 0.35541 0.50000 +O O 0.67152 0.44135 0.44621 +Si Si 0.00000 0.59316 0.68499 +Al Al 0.00000 0.68715 0.59042 +O O 0.00000 0.61367 0.74663 +O O 0.00000 0.64459 0.64591 +O O 0.05379 0.55865 0.67152 +Si Si 0.40684 0.50000 0.68499 +Al Al 0.31285 0.50000 0.59042 +O O 0.38633 0.50000 0.74663 +O O 0.35541 0.50000 0.64591 +O O 0.44135 0.55379 0.67152 +Si Si 0.50000 0.90684 0.68499 +Al Al 0.50000 0.81285 0.59042 +O O 0.50000 0.88633 0.74663 +O O 0.50000 0.85541 0.64591 +O O 0.44621 0.94135 0.67152 +Si Si 0.09316 0.00000 0.68499 +Al Al 0.18715 0.00000 0.59042 +O O 0.11367 0.00000 0.74663 +O O 0.14459 0.00000 0.64591 +O O 0.05865 0.94621 0.67152 +Si Si 0.00000 0.40684 0.31501 +Al Al 0.00000 0.31285 0.40958 +O O 0.00000 0.38633 0.25337 +O O 0.00000 0.35541 0.35409 +O O 0.05379 0.44135 0.32848 +Si Si 0.59316 0.50000 0.31501 +Al Al 0.68715 0.50000 0.40958 +O O 0.61367 0.50000 0.25337 +O O 0.64459 0.50000 0.35409 +O O 0.55865 0.55379 0.32848 +Si Si 0.50000 0.09316 0.31501 +Al Al 0.50000 0.18715 0.40958 +O O 0.50000 0.11367 0.25337 +O O 0.50000 0.14459 0.35409 +O O 0.44621 0.05865 0.32848 +Si Si 0.90684 0.00000 0.31501 +Al Al 0.81285 0.00000 0.40958 +O O 0.88633 0.00000 0.25337 +O O 0.85541 0.00000 0.35409 +O O 0.94135 0.94621 0.32848 +Si Si 0.18499 0.50000 0.59316 +Al Al 0.09042 0.50000 0.68715 +O O 0.24663 0.50000 0.61367 +O O 0.14591 0.50000 0.64459 +O O 0.17152 0.55379 0.55865 +Si Si 0.50000 0.68499 0.59316 +Al Al 0.50000 0.59042 0.68715 +O O 0.50000 0.74663 0.61367 +O O 0.50000 0.64591 0.64459 +O O 0.44621 0.67152 0.55865 +Si Si 0.31501 0.00000 0.59316 +Al Al 0.40958 0.00000 0.68715 +O O 0.25337 0.00000 0.61367 +O O 0.35409 0.00000 0.64459 +O O 0.32848 0.94621 0.55865 +Si Si 0.00000 0.81501 0.59316 +Al Al 0.00000 0.90958 0.68715 +O O 0.00000 0.75337 0.61367 +O O 0.00000 0.85409 0.64459 +O O 0.05379 0.82848 0.55865 +Si Si 0.18499 0.50000 0.40684 +Al Al 0.09042 0.50000 0.31285 +O O 0.24663 0.50000 0.38633 +O O 0.14591 0.50000 0.35541 +O O 0.17152 0.44621 0.44135 +Si Si 0.50000 0.68499 0.40684 +Al Al 0.50000 0.59042 0.31285 +O O 0.50000 0.74663 0.38633 +O O 0.50000 0.64591 0.35541 +O O 0.55379 0.67152 0.44135 +Si Si 0.31501 0.00000 0.40684 +Al Al 0.40958 0.00000 0.31285 +O O 0.25337 0.00000 0.38633 +O O 0.35409 0.00000 0.35541 +O O 0.32848 0.05379 0.44135 +Si Si 0.00000 0.81501 0.40684 +Al Al 0.00000 0.90958 0.31285 +O O 0.00000 0.75337 0.38633 +O O 0.00000 0.85409 0.35541 +O O 0.94621 0.82848 0.44135 +Si Si 0.09316 0.68499 0.50000 +Al Al 0.18715 0.59042 0.50000 +O O 0.11367 0.74663 0.50000 +O O 0.14459 0.64591 0.50000 +O O 0.05865 0.67152 0.55379 +Si Si 0.09316 0.81501 0.00000 +Al Al 0.18715 0.90958 0.00000 +O O 0.11367 0.75337 0.00000 +O O 0.14459 0.85409 0.00000 +O O 0.05865 0.82848 0.94621 +Si Si 0.18499 0.59316 0.00000 +Al Al 0.09042 0.68715 0.00000 +O O 0.24663 0.61367 0.00000 +O O 0.14591 0.64459 0.00000 +O O 0.17152 0.55865 0.94621 +Si Si 0.40684 0.68499 0.00000 +Al Al 0.31285 0.59042 0.00000 +O O 0.38633 0.74663 0.00000 +O O 0.35541 0.64591 0.00000 +O O 0.44135 0.67152 0.94621 +Si Si 0.31501 0.40684 0.50000 +Al Al 0.40958 0.31285 0.50000 +O O 0.25337 0.38633 0.50000 +O O 0.35409 0.35541 0.50000 +O O 0.32848 0.44135 0.44621 +Si Si 0.90684 0.31501 0.50000 +Al Al 0.81285 0.40958 0.50000 +O O 0.88633 0.25337 0.50000 +O O 0.85541 0.35409 0.50000 +O O 0.94135 0.32848 0.55379 +Si Si 0.18499 0.40684 0.00000 +Al Al 0.09042 0.31285 0.00000 +O O 0.24663 0.38633 0.00000 +O O 0.14591 0.35541 0.00000 +O O 0.17152 0.44135 0.05379 +Si Si 0.31501 0.09316 0.00000 +Al Al 0.40958 0.18715 0.00000 +O O 0.25337 0.11367 0.00000 +O O 0.35409 0.14459 0.00000 +O O 0.32848 0.05865 0.05379 +O O 0.94621 0.44135 0.32848 +O O 0.55865 0.44621 0.32848 +O O 0.55379 0.05865 0.32848 +O O 0.94135 0.05379 0.32848 +O O 0.94621 0.55865 0.67152 +O O 0.44135 0.44621 0.67152 +O O 0.55379 0.94135 0.67152 +O O 0.05865 0.05379 0.67152 +Si Si 0.81501 0.50000 0.40684 +Al Al 0.90958 0.50000 0.31285 +O O 0.75337 0.50000 0.38633 +O O 0.85409 0.50000 0.35541 +O O 0.82848 0.44621 0.44135 +Si Si 0.50000 0.31501 0.40684 +Al Al 0.50000 0.40958 0.31285 +O O 0.50000 0.25337 0.38633 +O O 0.50000 0.35409 0.35541 +O O 0.55379 0.32848 0.44135 +Si Si 0.68499 0.00000 0.40684 +Al Al 0.59042 0.00000 0.31285 +O O 0.74663 0.00000 0.38633 +O O 0.64591 0.00000 0.35541 +O O 0.67152 0.05379 0.44135 +Si Si 0.00000 0.18499 0.40684 +Al Al 0.00000 0.09042 0.31285 +O O 0.00000 0.24663 0.38633 +O O 0.00000 0.14591 0.35541 +O O 0.94621 0.17152 0.44135 +Si Si 0.81501 0.50000 0.59316 +Al Al 0.90958 0.50000 0.68715 +O O 0.75337 0.50000 0.61367 +O O 0.85409 0.50000 0.64459 +O O 0.82848 0.55379 0.55865 +Si Si 0.50000 0.31501 0.59316 +Al Al 0.50000 0.40958 0.68715 +O O 0.50000 0.25337 0.61367 +O O 0.50000 0.35409 0.64459 +O O 0.44621 0.32848 0.55865 +Si Si 0.68499 0.00000 0.59316 +Al Al 0.59042 0.00000 0.68715 +O O 0.74663 0.00000 0.61367 +O O 0.64591 0.00000 0.64459 +O O 0.67152 0.94621 0.55865 +Si Si 0.00000 0.18499 0.59316 +Al Al 0.00000 0.09042 0.68715 +O O 0.00000 0.24663 0.61367 +O O 0.00000 0.14591 0.64459 +O O 0.05379 0.17152 0.55865 +O O 0.94135 0.32848 0.44621 +Si Si 0.90684 0.18499 0.00000 +Al Al 0.81285 0.09042 0.00000 +O O 0.88633 0.24663 0.00000 +O O 0.85541 0.14591 0.00000 +O O 0.94135 0.17152 0.05379 +Si Si 0.81501 0.40684 0.00000 +Al Al 0.90958 0.31285 0.00000 +O O 0.75337 0.38633 0.00000 +O O 0.85409 0.35541 0.00000 +O O 0.82848 0.44135 0.05379 +Si Si 0.59316 0.31501 0.00000 +Al Al 0.68715 0.40958 0.00000 +O O 0.61367 0.25337 0.00000 +O O 0.64459 0.35409 0.00000 +O O 0.55865 0.32848 0.05379 +Si Si 0.68499 0.59316 0.50000 +Al Al 0.59042 0.68715 0.50000 +O O 0.74663 0.61367 0.50000 +O O 0.64591 0.64459 0.50000 +O O 0.67152 0.55865 0.55379 +O O 0.05865 0.67152 0.44621 +Si Si 0.81501 0.59316 0.00000 +Al Al 0.90958 0.68715 0.00000 +O O 0.75337 0.61367 0.00000 +O O 0.85409 0.64459 0.00000 +O O 0.82848 0.55865 0.94621 +Si Si 0.68499 0.90684 0.00000 +Al Al 0.59042 0.81285 0.00000 +O O 0.74663 0.88633 0.00000 +O O 0.64591 0.85541 0.00000 +O O 0.67152 0.94135 0.94621 +Si Si 0.50000 0.09316 0.68499 +Al Al 0.50000 0.18715 0.59042 +O O 0.50000 0.11367 0.74663 +O O 0.50000 0.14459 0.64591 +O O 0.55379 0.05865 0.67152 +Si Si 0.90684 0.00000 0.68499 +Al Al 0.81285 0.00000 0.59042 +O O 0.88633 0.00000 0.74663 +O O 0.85541 0.00000 0.64591 +O O 0.94135 0.05379 0.67152 +Si Si 0.00000 0.40684 0.68499 +Al Al 0.00000 0.31285 0.59042 +O O 0.00000 0.38633 0.74663 +O O 0.00000 0.35541 0.64591 +O O 0.94621 0.44135 0.67152 +Si Si 0.59316 0.50000 0.68499 +Al Al 0.68715 0.50000 0.59042 +O O 0.61367 0.50000 0.74663 +O O 0.64459 0.50000 0.64591 +O O 0.55865 0.44621 0.67152 +Si Si 0.50000 0.90684 0.31501 +Al Al 0.50000 0.81285 0.40958 +O O 0.50000 0.88633 0.25337 +O O 0.50000 0.85541 0.35409 +O O 0.55379 0.94135 0.32848 +Si Si 0.09316 0.00000 0.31501 +Al Al 0.18715 0.00000 0.40958 +O O 0.11367 0.00000 0.25337 +O O 0.14459 0.00000 0.35409 +O O 0.05865 0.05379 0.32848 +Si Si 0.00000 0.59316 0.31501 +Al Al 0.00000 0.68715 0.40958 +O O 0.00000 0.61367 0.25337 +O O 0.00000 0.64459 0.35409 +O O 0.94621 0.55865 0.32848 +Si Si 0.40684 0.50000 0.31501 +Al Al 0.31285 0.50000 0.40958 +O O 0.38633 0.50000 0.25337 +O O 0.35541 0.50000 0.35409 +O O 0.44135 0.44621 0.32848 +O O 0.67152 0.05379 0.55865 +O O 0.94621 0.17152 0.55865 +O O 0.82848 0.44621 0.55865 +O O 0.55379 0.32848 0.55865 +O O 0.67152 0.94621 0.44135 +O O 0.05379 0.17152 0.44135 +O O 0.82848 0.55379 0.44135 +O O 0.44621 0.32848 0.44135 +Si Si 0.59316 0.18499 0.50000 +Al Al 0.68715 0.09042 0.50000 +O O 0.61367 0.24663 0.50000 +O O 0.64459 0.14591 0.50000 +O O 0.55865 0.17152 0.55379 +O O 0.55865 0.32848 0.94621 +O O 0.67152 0.05865 0.94621 +O O 0.94135 0.17152 0.94621 +O O 0.82848 0.94135 0.44621 +Si Si 0.40684 0.81501 0.50000 +Al Al 0.31285 0.90958 0.50000 +O O 0.38633 0.75337 0.50000 +O O 0.35541 0.85409 0.50000 +O O 0.44135 0.82848 0.55379 +O O 0.67152 0.94135 0.05379 +O O 0.82848 0.55865 0.05379 +O O 0.44621 0.94135 0.32848 +O O 0.05865 0.94621 0.32848 +O O 0.05379 0.55865 0.32848 +O O 0.44135 0.55379 0.32848 +O O 0.44621 0.05865 0.67152 +O O 0.94135 0.94621 0.67152 +O O 0.05379 0.44135 0.67152 +O O 0.55865 0.55379 0.67152 +O O 0.32848 0.94621 0.44135 +O O 0.05379 0.82848 0.44135 +O O 0.17152 0.55379 0.44135 +O O 0.44621 0.67152 0.44135 +O O 0.32848 0.05379 0.55865 +O O 0.94621 0.82848 0.55865 +O O 0.17152 0.44621 0.55865 +O O 0.55379 0.67152 0.55865 +O O 0.44135 0.82848 0.44621 +O O 0.44135 0.67152 0.05379 +O O 0.32848 0.94135 0.05379 +O O 0.05865 0.82848 0.05379 +O O 0.17152 0.05865 0.55379 +O O 0.55865 0.17152 0.44621 +O O 0.32848 0.05865 0.94621 +O O 0.17152 0.44135 0.94621 +Si Si 0.50000 0.59316 0.18499 +Al Al 0.50000 0.68715 0.09042 +O O 0.50000 0.61367 0.24663 +O O 0.50000 0.64459 0.14591 +O O 0.55379 0.55865 0.17152 +Si Si 0.90684 0.50000 0.18499 +Al Al 0.81285 0.50000 0.09042 +O O 0.88633 0.50000 0.24663 +O O 0.85541 0.50000 0.14591 +O O 0.94135 0.55379 0.17152 +Si Si 0.00000 0.90684 0.18499 +Al Al 0.00000 0.81285 0.09042 +O O 0.00000 0.88633 0.24663 +O O 0.00000 0.85541 0.14591 +O O 0.94621 0.94135 0.17152 +Si Si 0.59316 0.00000 0.18499 +Al Al 0.68715 0.00000 0.09042 +O O 0.61367 0.00000 0.24663 +O O 0.64459 0.00000 0.14591 +O O 0.55865 0.94621 0.17152 +Si Si 0.50000 0.40684 0.81501 +Al Al 0.50000 0.31285 0.90958 +O O 0.50000 0.38633 0.75337 +O O 0.50000 0.35541 0.85409 +O O 0.55379 0.44135 0.82848 +Si Si 0.09316 0.50000 0.81501 +Al Al 0.18715 0.50000 0.90958 +O O 0.11367 0.50000 0.75337 +O O 0.14459 0.50000 0.85409 +O O 0.05865 0.55379 0.82848 +Si Si 0.00000 0.09316 0.81501 +Al Al 0.00000 0.18715 0.90958 +O O 0.00000 0.11367 0.75337 +O O 0.00000 0.14459 0.85409 +O O 0.94621 0.05865 0.82848 +Si Si 0.40684 0.00000 0.81501 +Al Al 0.31285 0.00000 0.90958 +O O 0.38633 0.00000 0.75337 +O O 0.35541 0.00000 0.85409 +O O 0.44135 0.94621 0.82848 +O O 0.67152 0.55379 0.05865 +O O 0.94621 0.67152 0.05865 +O O 0.82848 0.94621 0.05865 +O O 0.55379 0.82848 0.05865 +O O 0.67152 0.44621 0.94135 +O O 0.05379 0.67152 0.94135 +O O 0.82848 0.05379 0.94135 +O O 0.44621 0.82848 0.94135 +Si Si 0.59316 0.68499 0.00000 +Al Al 0.68715 0.59042 0.00000 +O O 0.61367 0.74663 0.00000 +O O 0.64459 0.64591 0.00000 +O O 0.55865 0.67152 0.05379 +O O 0.55865 0.82848 0.44621 +O O 0.67152 0.55865 0.44621 +O O 0.94135 0.67152 0.44621 +O O 0.82848 0.44135 0.94621 +Si Si 0.40684 0.31501 0.00000 +Al Al 0.31285 0.40958 0.00000 +O O 0.38633 0.25337 0.00000 +O O 0.35541 0.35409 0.00000 +O O 0.44135 0.32848 0.05379 +O O 0.67152 0.44135 0.55379 +O O 0.82848 0.05865 0.55379 +O O 0.44621 0.44135 0.82848 +O O 0.05865 0.44621 0.82848 +O O 0.05379 0.05865 0.82848 +O O 0.44135 0.05379 0.82848 +O O 0.44621 0.55865 0.17152 +O O 0.94135 0.44621 0.17152 +O O 0.05379 0.94135 0.17152 +O O 0.55865 0.05379 0.17152 +O O 0.32848 0.44621 0.94135 +O O 0.05379 0.32848 0.94135 +O O 0.17152 0.05379 0.94135 +O O 0.44621 0.17152 0.94135 +O O 0.32848 0.55379 0.05865 +O O 0.94621 0.32848 0.05865 +O O 0.17152 0.94621 0.05865 +O O 0.55379 0.17152 0.05865 +O O 0.44135 0.32848 0.94621 +O O 0.44135 0.17152 0.55379 +O O 0.32848 0.44135 0.55379 +O O 0.05865 0.32848 0.55379 +O O 0.17152 0.55865 0.05379 +O O 0.55865 0.67152 0.94621 +O O 0.32848 0.55865 0.44621 +O O 0.17152 0.94135 0.44621 diff --git a/examples/run_zeoisot.py b/examples/run_zeoisot.py new file mode 100644 index 00000000..b5e6d3c7 --- /dev/null +++ b/examples/run_zeoisot.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +"""Run example isotherm calculation with COF-1 framework.""" + +from __future__ import absolute_import +from __future__ import print_function +import os +import click + +from aiida.engine import run +from aiida.plugins import DataFactory +from aiida.orm import Code, Dict, Str, Float, Int + +# Workchain objects +#IsothermWorkChain = WorkflowFactory('lsmo.isotherm') # pylint: disable=invalid-name +from aiida_lsmo.workchains.zeoisotherm import ZeoIsothermWorkChain + +# Data objects +CifData = DataFactory('cif') # pylint: disable=invalid-name +NetworkParameters = DataFactory('zeopp.parameters') # pylint: disable=invalid-name + + +@click.command('cli') +@click.argument('raspa_code_label') +@click.argument('zeopp_code_label') +def main(raspa_code_label, zeopp_code_label): + """Prepare inputs and submit the Isotherm workchain. + Usage: verdi run run_IsothermWorkChain_COF-1.py raspa@localhost network@localhost""" + + builder = ZeoIsothermWorkChain.get_builder() + + builder.metadata.label = "test" # pylint: disable=no-member + + builder.raspa_base.raspa.code = Code.get_from_string(raspa_code_label) # pylint: disable=no-member + builder.zeopp.code = Code.get_from_string(zeopp_code_label) # pylint: disable=no-member + + options = { + "resources": { + "num_machines": 1, + "tot_num_mpiprocs": 1, + }, + "max_wallclock_seconds": 1 * 60 * 60, + "withmpi": False, + } + + builder.raspa_base.raspa.metadata.options = options # pylint: disable=no-member + builder.zeopp.metadata.options = options # pylint: disable=no-member + builder.raspa_base.max_iterations = Int(1) # pylint: disable=no-member + + builder.structure = CifData(file=os.path.abspath('data/zeo_LTA_p91.cif'), label="zeo_LTA_p91") + builder.sial_ratio = Float(1.101) + builder.cation = Str('Na') + + builder.molecule = Str('co2') + builder.parameters = Dict( + dict={ + 'ff_framework': 'Zeo-test', # Default: UFF + 'temperature': 300, # (K) Note: higher temperature will have less adsorbate and it is faster + 'zeopp_volpo_samples': 1000, # Default: 1e5 *NOTE: default is good for standard real-case! + 'zeopp_block_samples': 10, # Default: 100 + "raspa_nvt_cations_cycles": 1000, + 'raspa_widom_cycles': 100, # Default: 1e5 + 'raspa_gcmc_init_cycles': 10, # Default: 1e3 + 'raspa_gcmc_prod_cycles': 100, # Default: 1e4 + 'pressure_min': 0.001, # Default: 0.001 (bar) + 'pressure_max': 3, # Default: 10 (bar) + }) + + run(builder) + + +if __name__ == '__main__': + print('Testing Zeolite Isotherm work chain.') + main() # pylint: disable=no-value-for-parameter + +# EOF diff --git a/examples/test_ffbuilder.py b/examples/test_ffbuilder_uff.py similarity index 100% rename from examples/test_ffbuilder.py rename to examples/test_ffbuilder_uff.py diff --git a/examples/test_ffbuilder_zeolite.py b/examples/test_ffbuilder_zeolite.py new file mode 100644 index 00000000..b4db084f --- /dev/null +++ b/examples/test_ffbuilder_zeolite.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +"""Test for ff_builder""" +from __future__ import absolute_import +from __future__ import print_function +from aiida.orm import Dict +from aiida.plugins import CalculationFactory +from aiida.engine import run_get_node + +# Calculation objects +FFBuilder = CalculationFactory("lsmo.ff_builder") # pylint: disable=invalid-name + +ff_parameters = Dict( # pylint: disable=invalid-name + dict={ + 'ff_framework': 'Zeo-test', + 'ff_molecules': { + 'Na': 'Zeo-test', + 'CO2': 'TraPPE', + }, + 'shifted': False, + 'tail_corrections': True, + 'mixing_rule': 'Lorentz-Berthelot', + 'separate_interactions': False + }) + +results, node = run_get_node(FFBuilder, ff_parameters) # pylint: disable=invalid-name +print("Terminated ff_builder calcfunction, pk:", node.pk) +for key, val in results.items(): + #filepath = os.path.join(val._repository._get_base_folder().abspath, val.filename) + print("Output:", val.pk, key) From f38d61ed189e44fb7dd9e6b24cba380f2319b213 Mon Sep 17 00:00:00 2001 From: Daniele Ongari Date: Mon, 24 Feb 2020 16:43:54 +0100 Subject: [PATCH 2/3] add annealing for cations --- aiida_lsmo/workchains/zeoisotherm.py | 85 ++++++++++++++++++---------- 1 file changed, 54 insertions(+), 31 deletions(-) diff --git a/aiida_lsmo/workchains/zeoisotherm.py b/aiida_lsmo/workchains/zeoisotherm.py index 96dab6ea..1cb50a08 100644 --- a/aiida_lsmo/workchains/zeoisotherm.py +++ b/aiida_lsmo/workchains/zeoisotherm.py @@ -1,9 +1,6 @@ -# -*- coding: utf-8 -*- """ZeoIsotherm work chain""" -from __future__ import absolute_import import os -from six.moves import range from aiida.plugins import CalculationFactory, DataFactory, WorkflowFactory from aiida.orm import Dict, Float, Str, List, SinglefileData @@ -31,6 +28,9 @@ "ff_tail_corrections": True, # (bool) Apply tail corrections. "ff_shifted": False, # (bool) Shift or truncate the potential at cutoff. "ff_cutoff": 12.0, # (float) CutOff truncation for the VdW interactions (Angstrom). + "cations_equilibration_temp_list": [ + 1000, 800, 600, 500, 400 + ], # (list) List of temperatures for the equilibrations of cations: last will be "temperature". "temperature": 300, # (float) Temperature of the simulation. "temperature_list": None, # (list) To be used by IsothermMultiTempWorkChain. "zeopp_volpo_samples": int(1e5), # (int) Number of samples for VOLPO calculation (per UC volume). @@ -229,7 +229,7 @@ def define(cls, spec): spec.outline( cls.setup, cls.run_zeopp, # computes volpo and blocks - cls.run_equilibrate_cation, + while_(cls.should_equilibrate_cations)(cls.run_equilibrate_cations,), if_(cls.should_run_widom)( # run Widom only if porous cls.run_raspa_widom, # run raspa widom calculation if_(cls.should_run_gcmc)( # Kh is high enough @@ -265,10 +265,14 @@ def setup(self): # Get integer temperature in context for easy reports self.ctx.temperature = int(round(self.ctx.parameters['temperature'])) - # Get the number of cations needed - natoms = self.inputs.structure.get_ase().get_number_of_atoms() + # Get the number of cations needed and other data for the cation equilibration + natoms = self.inputs.structure.get_ase().get_global_number_of_atoms() ntatoms = natoms / 3 self.ctx.ncations = int(ntatoms / (self.inputs.sial_ratio + 1)) + self.ctx.current_equilibration_id = 0 + self.ctx.cations_equilibration_temp_list = self.ctx.parameters["cations_equilibration_temp_list"] + [ + self.ctx.parameters['temperature'] + ] # Understand if IsothermMultiTempWorkChain is calling this work chain if "geometric" in self.inputs: @@ -308,7 +312,7 @@ def run_zeopp(self): self.report("Running zeo++ block and volpo Calculation<{}>".format(running.id)) return ToContext(zeopp=running) - def _get_nvt_cation(self): + def _get_nvt_cations(self): """Write Raspa input parameters from scratch, for a NVT calculation""" param = { @@ -328,7 +332,8 @@ def _get_nvt_cation(self): "System": { "framework_1": { "type": "Framework", - "ExternalTemperature": self.ctx.parameters['temperature'], + "UnitCells": "1 1 1", #NOTE: it is expanded before! + "ExternalTemperature": self.ctx.cations_equilibration_temp_list[0], } }, "Component": { @@ -341,35 +346,53 @@ def _get_nvt_cation(self): }, }, } - - # Check particular conditions and settings - # DISABLED: I want to have already from the beginning a larger unit cell (so, no problems with ncations) - # mult = check_resize_unit_cell(self.inputs.structure, 2 * self.ctx.parameters['ff_cutoff']) - # param["System"]["framework_1"]["UnitCells"] = "{} {} {}".format(mult[0], mult[1], mult[2]) return param - def run_equilibrate_cation(self): - """Run a calculation in Raspa to equilibrate postion of cations.""" + def should_equilibrate_cations(self): + """Decide whether to run the cations equilibration.""" + if self.ctx.ncations == 0: + return False + return self.ctx.current_equilibration_id < len(self.ctx.cations_equilibration_temp_list) - # Initialize the input for raspa_base, which later will need only minor updates for GCMC - self.ctx.inp = self.exposed_inputs(RaspaBaseWorkChain, 'raspa_base') - self.ctx.inp['metadata']['label'] = "RaspaWidom" - self.ctx.inp['metadata']['call_link_label'] = "run_raspa_cation" + def run_equilibrate_cations(self): + """Run NVT for the equilibration of the cations as a list of calculations.""" - # Create input parameters and pass the framework - self.ctx.raspa_param = self._get_nvt_cation() - self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) - self.ctx.inp['raspa']['framework'] = {"framework_1": self.inputs.structure} + if self.ctx.current_equilibration_id == 0: - # Generate the force field with the ff_builder - ff_params = get_ff_parameters(self.ctx.molecule, self.inputs.cation, self.ctx.parameters) - files_dict = FFBuilder(ff_params) - self.ctx.inp['raspa']['file'] = files_dict + # Initialize the input for raspa_base, which later will need only minor updates for GCMC + self.ctx.inp = self.exposed_inputs(RaspaBaseWorkChain, 'raspa_base') - running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) - self.report("Running Raspa NVT @ {}K for equilibrating cations".format(self.ctx.temperature)) + # Create input parameters and pass the framework + self.ctx.raspa_param = self._get_nvt_cations() + self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) + self.ctx.inp['raspa']['framework'] = {"framework_1": self.inputs.structure} - return ToContext(raspa_cation=running) + # Generate the force field with the ff_builder + ff_params = get_ff_parameters(self.ctx.molecule, self.inputs.cation, self.ctx.parameters) + files_dict = FFBuilder(ff_params) + self.ctx.inp['raspa']['file'] = files_dict + + elif self.ctx.current_equilibration_id > 0: + + # Update Parameters + self.ctx.raspa_param["System"]["framework_1"][ + 'ExternalTemperature'] = self.ctx.cations_equilibration_temp_list[self.ctx.current_equilibration_id] + self.ctx.raspa_param["Component"][self.inputs.cation.value]["CreateNumberOfMolecules"] = 0 + self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) + + # Update restart (if present, i.e., if current_p_index>0) + self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_nvt_cations[ + self.ctx.current_equilibration_id - 1].outputs.retrieved + + # Update metadata and submit + self.ctx.inp['metadata']['label'] = "RaspaNVTCations{}".format(self.ctx.current_equilibration_id) + self.ctx.inp['metadata']['call_link_label'] = "run_raspa_cations_{}".format(self.ctx.current_equilibration_id) + running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) + self.report("Running Raspa NVT @ {}K for equilibrating cations ({}/{})".format( + int(self.ctx.cations_equilibration_temp_list[self.ctx.current_equilibration_id]), + self.ctx.current_equilibration_id + 1, len(self.ctx.cations_equilibration_temp_list))) + self.ctx.current_equilibration_id += 1 + return ToContext(raspa_nvt_cations=append_(running)) def should_run_widom(self): """Submit widom calculation only if there is some accessible volume, @@ -431,7 +454,7 @@ def run_raspa_widom(self): self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) if self.ctx.geom['Number_of_blocking_spheres'] > 0 and self.ctx.multitemp_mode != 'run_single_temp': self.ctx.inp["raspa"]["block_pocket"] = {"block_file": self.ctx.zeopp.outputs.block} - self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_cation.outputs.retrieved + self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_nvt_cations[-1].outputs.retrieved running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) self.report("Running Raspa Widom @ {}K for the Henry coefficient".format(self.ctx.temperature)) From 13f87f3b264891350db9ce0aa5afcf1472f3db26 Mon Sep 17 00:00:00 2001 From: Leopold Talirz Date: Fri, 13 Nov 2020 15:50:37 +0100 Subject: [PATCH 3/3] apply string fixers --- aiida_lsmo/workchains/zeoisotherm.py | 229 ++++++++++++++------------- examples/run_zeoisot.py | 16 +- examples/test_ffbuilder_zeolite.py | 6 +- 3 files changed, 126 insertions(+), 125 deletions(-) diff --git a/aiida_lsmo/workchains/zeoisotherm.py b/aiida_lsmo/workchains/zeoisotherm.py index 1cb50a08..97801373 100644 --- a/aiida_lsmo/workchains/zeoisotherm.py +++ b/aiida_lsmo/workchains/zeoisotherm.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """ZeoIsotherm work chain""" import os @@ -22,30 +23,30 @@ # Deafault parameters ISOTHERMPARAMETERS_DEFAULT = { #TODO: create IsothermParameters instead of Dict # pylint: disable=fixme - "ff_framework": "Zeo-test", # (str) Forcefield of the structure. - "ff_separate_interactions": False, # (bool) Use "separate_interactions" in the FF builder. - "ff_mixing_rule": "Lorentz-Berthelot", # (string) Choose 'Lorentz-Berthelot' or 'Jorgensen'. - "ff_tail_corrections": True, # (bool) Apply tail corrections. - "ff_shifted": False, # (bool) Shift or truncate the potential at cutoff. - "ff_cutoff": 12.0, # (float) CutOff truncation for the VdW interactions (Angstrom). - "cations_equilibration_temp_list": [ + 'ff_framework': 'Zeo-test', # (str) Forcefield of the structure. + 'ff_separate_interactions': False, # (bool) Use "separate_interactions" in the FF builder. + 'ff_mixing_rule': 'Lorentz-Berthelot', # (string) Choose 'Lorentz-Berthelot' or 'Jorgensen'. + 'ff_tail_corrections': True, # (bool) Apply tail corrections. + 'ff_shifted': False, # (bool) Shift or truncate the potential at cutoff. + 'ff_cutoff': 12.0, # (float) CutOff truncation for the VdW interactions (Angstrom). + 'cations_equilibration_temp_list': [ 1000, 800, 600, 500, 400 ], # (list) List of temperatures for the equilibrations of cations: last will be "temperature". - "temperature": 300, # (float) Temperature of the simulation. - "temperature_list": None, # (list) To be used by IsothermMultiTempWorkChain. - "zeopp_volpo_samples": int(1e5), # (int) Number of samples for VOLPO calculation (per UC volume). - "zeopp_block_samples": int(100), # (int) Number of samples for BLOCK calculation (per A^3). - "raspa_minKh": 1e-10, # (float) If Henry coefficient < raspa_minKh do not run the isotherm (mol/kg/Pa). - "raspa_verbosity": 10, # (int) Print stats every: number of cycles / raspa_verbosity. - "raspa_widom_cycles": int(1e5), # (int) Number of Widom cycles. - "raspa_nvt_cations_cycles": int(1e4), # (int) Number of GCMC production cycles. - "raspa_gcmc_init_cycles": int(1e3), # (int) Number of GCMC initialization cycles. - "raspa_gcmc_prod_cycles": int(1e4), # (int) Number of GCMC production cycles. - "pressure_list": None, # (list) Pressure list for the isotherm (bar): if given it will skip to guess it. - "pressure_precision": 0.1, # (float) Precision in the sampling of the isotherm: 0.1 ok, 0.05 for high resolution. - "pressure_maxstep": 5, # (float) Max distance between pressure points (bar). - "pressure_min": 0.001, # (float) Lower pressure to sample (bar). - "pressure_max": 10 # (float) Upper pressure to sample (bar). + 'temperature': 300, # (float) Temperature of the simulation. + 'temperature_list': None, # (list) To be used by IsothermMultiTempWorkChain. + 'zeopp_volpo_samples': int(1e5), # (int) Number of samples for VOLPO calculation (per UC volume). + 'zeopp_block_samples': int(100), # (int) Number of samples for BLOCK calculation (per A^3). + 'raspa_minKh': 1e-10, # (float) If Henry coefficient < raspa_minKh do not run the isotherm (mol/kg/Pa). + 'raspa_verbosity': 10, # (int) Print stats every: number of cycles / raspa_verbosity. + 'raspa_widom_cycles': int(1e5), # (int) Number of Widom cycles. + 'raspa_nvt_cations_cycles': int(1e4), # (int) Number of GCMC production cycles. + 'raspa_gcmc_init_cycles': int(1e3), # (int) Number of GCMC initialization cycles. + 'raspa_gcmc_prod_cycles': int(1e4), # (int) Number of GCMC production cycles. + 'pressure_list': None, # (list) Pressure list for the isotherm (bar): if given it will skip to guess it. + 'pressure_precision': 0.1, # (float) Precision in the sampling of the isotherm: 0.1 ok, 0.05 for high resolution. + 'pressure_maxstep': 5, # (float) Max distance between pressure points (bar). + 'pressure_min': 0.001, # (float) Lower pressure to sample (bar). + 'pressure_max': 10 # (float) Upper pressure to sample (bar). } @@ -55,7 +56,7 @@ def get_molecule_dict(molecule_name): """Get a Dict from the isotherm_molecules.yaml""" import ruamel.yaml as yaml thisdir = os.path.dirname(os.path.abspath(__file__)) - yamlfile = os.path.join(thisdir, "isotherm_data", "isotherm_molecules.yaml") + yamlfile = os.path.join(thisdir, 'isotherm_data', 'isotherm_molecules.yaml') with open(yamlfile, 'r') as stream: yaml_dict = yaml.safe_load(stream) molecule_dict = yaml_dict[molecule_name.value] @@ -65,7 +66,7 @@ def get_molecule_dict(molecule_name): @calcfunction def get_zeopp_parameters(molecule_dict, isotparam): """Get the ZeoppParameters from the inputs of the workchain""" - probe_rad = molecule_dict["proberad"] + probe_rad = molecule_dict['proberad'] param_dict = { 'ha': 'DEF', 'volpo': [probe_rad, probe_rad, isotparam['zeopp_volpo_samples']], @@ -97,10 +98,10 @@ def choose_pressure_points(inp_param, geom, raspa_widom_out): """If 'presure_list' is not provide, model the isotherm as single-site langmuir and return the most important pressure points to evaluate for an isotherm, in a List. """ - if inp_param["pressure_list"]: - pressure_points = inp_param["pressure_list"] + if inp_param['pressure_list']: + pressure_points = inp_param['pressure_list'] else: - khenry = list(raspa_widom_out["framework_1"]["components"].values())[1]['henry_coefficient_average'] #mol/kg/Pa + khenry = list(raspa_widom_out['framework_1']['components'].values())[1]['henry_coefficient_average'] #mol/kg/Pa b_value = khenry / geom['Estimated_saturation_loading'] * 1e5 #(1/bar) pressure_points = [inp_param['pressure_min']] while True: @@ -123,7 +124,7 @@ def get_geometric_dict(zeopp_out, molecule): geometric_dict.update({ 'Estimated_saturation_loading': zeopp_out['POAV_cm^3/g'] * molecule['molsatdens'], 'Estimated_saturation_loading_unit': 'mol/kg', - 'is_porous': geometric_dict["POAV_A^3"] > 0.000 + 'is_porous': geometric_dict['POAV_A^3'] > 0.000 }) return Dict(dict=geometric_dict) @@ -135,7 +136,7 @@ def get_output_parameters(geom_out, inp_params, widom_out=None, pressures=None, out_dict = geom_out.get_dict() if out_dict['is_porous'] and widom_out: - widom_out_mol = list(widom_out["framework_1"]["components"].values())[1] + widom_out_mol = list(widom_out['framework_1']['components'].values())[1] out_dict.update({ 'temperature': inp_params['temperature'], @@ -170,9 +171,9 @@ def get_output_parameters(geom_out, inp_params, widom_out=None, pressures=None, conv_ener = 1.0 / 120.273 # K to kJ/mol for i in range(len(pressures)): - gcmc_out = gcmc_out_dict['RaspaGCMC_{}'.format(i + 1)]["framework_1"] - gcmc_out_mol = list(gcmc_out["components"].values())[1] - conv_load = gcmc_out_mol["conversion_factor_molec_uc_to_mol_kg"] + gcmc_out = gcmc_out_dict['RaspaGCMC_{}'.format(i + 1)]['framework_1'] + gcmc_out_mol = list(gcmc_out['components'].values())[1] + conv_load = gcmc_out_mol['conversion_factor_molec_uc_to_mol_kg'] for label in ['loading_absolute_average', 'loading_absolute_dev']: isotherm[label].append(conv_load * gcmc_out_mol[label]) @@ -184,7 +185,7 @@ def get_output_parameters(geom_out, inp_params, widom_out=None, pressures=None, isotherm[label].append(None) out_dict.update({ - "isotherm": isotherm, + 'isotherm': isotherm, 'conversion_factor_molec_uc_to_cm3stp_cm3': gcmc_out_mol['conversion_factor_molec_uc_to_cm3stp_cm3'], 'conversion_factor_molec_uc_to_mg_g': gcmc_out_mol['conversion_factor_molec_uc_to_mg_g'], 'conversion_factor_molec_uc_to_mol_kg': gcmc_out_mol['conversion_factor_molec_uc_to_mol_kg'], @@ -209,19 +210,19 @@ def define(cls, spec): spec.input('sial_ratio', valid_type=Float, default=Float(1), help='Get the Si/Al ratio in the zeolite.') - spec.input("molecule", + spec.input('molecule', valid_type=(Str, Dict), help='Adsorbate molecule: settings to be read from the yaml.' + 'Advanced: input a Dict for non-standard settings.') - spec.input("cation", valid_type=Str, default=Str("Na"), help='Cation to be used for charged zeolites (Na/Ca)') + spec.input('cation', valid_type=Str, default=Str('Na'), help='Cation to be used for charged zeolites (Na/Ca)') #TODO: need to implement Ca+2 # pylint: disable=fixme - spec.input("parameters", + spec.input('parameters', valid_type=Dict, help='Parameters for the Isotherm workchain: will be merged with IsothermParameters_defaults.') - spec.input("geometric", + spec.input('geometric', valid_type=Dict, required=False, help='[Only used by IsothermMultiTempWorkChain] Already computed geometric properties') @@ -270,12 +271,12 @@ def setup(self): ntatoms = natoms / 3 self.ctx.ncations = int(ntatoms / (self.inputs.sial_ratio + 1)) self.ctx.current_equilibration_id = 0 - self.ctx.cations_equilibration_temp_list = self.ctx.parameters["cations_equilibration_temp_list"] + [ + self.ctx.cations_equilibration_temp_list = self.ctx.parameters['cations_equilibration_temp_list'] + [ self.ctx.parameters['temperature'] ] # Understand if IsothermMultiTempWorkChain is calling this work chain - if "geometric" in self.inputs: + if 'geometric' in self.inputs: self.ctx.multitemp_mode = 'run_single_temp' elif self.ctx.parameters['temperature_list']: self.ctx.multitemp_mode = 'run_geom_only' @@ -296,53 +297,53 @@ def run_zeopp(self): dict_merge( inputs, { 'metadata': { - 'label': "ZeoppVolpoBlock", + 'label': 'ZeoppVolpoBlock', 'call_link_label': 'run_zeopp_block_and_volpo', }, 'structure': self.inputs.structure, 'atomic_radii': SinglefileData( - file=os.path.join(os.path.dirname(os.path.abspath(__file__)), "isotherm_data", "DEFAULT.rad")), + file=os.path.join(os.path.dirname(os.path.abspath(__file__)), 'isotherm_data', 'DEFAULT.rad')), 'parameters': get_zeopp_parameters(self.ctx.molecule, self.ctx.parameters) }) running = self.submit(ZeoppCalculation, **inputs) - self.report("Running zeo++ block and volpo Calculation<{}>".format(running.id)) + self.report('Running zeo++ block and volpo Calculation<{}>'.format(running.id)) return ToContext(zeopp=running) def _get_nvt_cations(self): """Write Raspa input parameters from scratch, for a NVT calculation""" param = { - "GeneralSettings": { - "SimulationType": "MonteCarlo", - "NumberOfInitializationCycles": self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, - "NumberOfCycles": self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, - "PrintPropertiesEvery": int(1e10), - "PrintEvery": self.ctx.parameters['raspa_nvt_cations_cycles'] / self.ctx.parameters['raspa_verbosity'], - "RemoveAtomNumberCodeFromLabel": True, - "ModifyOxgensConnectedToAluminium": True, - "Forcefield": "Local", - "CutOff": self.ctx.parameters['ff_cutoff'], - "ChargeMethod": "Ewald", - "EwaldPrecision": 1e-6, + 'GeneralSettings': { + 'SimulationType': 'MonteCarlo', + 'NumberOfInitializationCycles': self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, + 'NumberOfCycles': self.ctx.parameters['raspa_nvt_cations_cycles'] / 2, + 'PrintPropertiesEvery': int(1e10), + 'PrintEvery': self.ctx.parameters['raspa_nvt_cations_cycles'] / self.ctx.parameters['raspa_verbosity'], + 'RemoveAtomNumberCodeFromLabel': True, + 'ModifyOxgensConnectedToAluminium': True, + 'Forcefield': 'Local', + 'CutOff': self.ctx.parameters['ff_cutoff'], + 'ChargeMethod': 'Ewald', + 'EwaldPrecision': 1e-6, }, - "System": { - "framework_1": { - "type": "Framework", - "UnitCells": "1 1 1", #NOTE: it is expanded before! - "ExternalTemperature": self.ctx.cations_equilibration_temp_list[0], + 'System': { + 'framework_1': { + 'type': 'Framework', + 'UnitCells': '1 1 1', #NOTE: it is expanded before! + 'ExternalTemperature': self.ctx.cations_equilibration_temp_list[0], } }, - "Component": { + 'Component': { self.inputs.cation.value: { - "MoleculeDefinition": "Local", - "ExtraFrameworkMolecule": True, - "TranslationProbability": 1.0, - "RandomTranslationProbability": 1.0, #Raspa warns to use this instead of ReinsertionProbability - "CreateNumberOfMolecules": self.ctx.ncations, + 'MoleculeDefinition': 'Local', + 'ExtraFrameworkMolecule': True, + 'TranslationProbability': 1.0, + 'RandomTranslationProbability': 1.0, #Raspa warns to use this instead of ReinsertionProbability + 'CreateNumberOfMolecules': self.ctx.ncations, }, }, } @@ -365,7 +366,7 @@ def run_equilibrate_cations(self): # Create input parameters and pass the framework self.ctx.raspa_param = self._get_nvt_cations() self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) - self.ctx.inp['raspa']['framework'] = {"framework_1": self.inputs.structure} + self.ctx.inp['raspa']['framework'] = {'framework_1': self.inputs.structure} # Generate the force field with the ff_builder ff_params = get_ff_parameters(self.ctx.molecule, self.inputs.cation, self.ctx.parameters) @@ -375,9 +376,9 @@ def run_equilibrate_cations(self): elif self.ctx.current_equilibration_id > 0: # Update Parameters - self.ctx.raspa_param["System"]["framework_1"][ + self.ctx.raspa_param['System']['framework_1'][ 'ExternalTemperature'] = self.ctx.cations_equilibration_temp_list[self.ctx.current_equilibration_id] - self.ctx.raspa_param["Component"][self.inputs.cation.value]["CreateNumberOfMolecules"] = 0 + self.ctx.raspa_param['Component'][self.inputs.cation.value]['CreateNumberOfMolecules'] = 0 self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) # Update restart (if present, i.e., if current_p_index>0) @@ -385,10 +386,10 @@ def run_equilibrate_cations(self): self.ctx.current_equilibration_id - 1].outputs.retrieved # Update metadata and submit - self.ctx.inp['metadata']['label'] = "RaspaNVTCations{}".format(self.ctx.current_equilibration_id) - self.ctx.inp['metadata']['call_link_label'] = "run_raspa_cations_{}".format(self.ctx.current_equilibration_id) + self.ctx.inp['metadata']['label'] = 'RaspaNVTCations{}'.format(self.ctx.current_equilibration_id) + self.ctx.inp['metadata']['call_link_label'] = 'run_raspa_cations_{}'.format(self.ctx.current_equilibration_id) running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) - self.report("Running Raspa NVT @ {}K for equilibrating cations ({}/{})".format( + self.report('Running Raspa NVT @ {}K for equilibrating cations ({}/{})'.format( int(self.ctx.cations_equilibration_temp_list[self.ctx.current_equilibration_id]), self.ctx.current_equilibration_id + 1, len(self.ctx.cations_equilibration_temp_list))) self.ctx.current_equilibration_id += 1 @@ -406,13 +407,13 @@ def should_run_widom(self): self.ctx.geom = get_geometric_dict(self.ctx.zeopp.outputs.output_parameters, self.ctx.molecule) if self.ctx.geom['is_porous']: - self.report("Found accessible pore volume: continue") - self.report("Found {} blocking spheres".format(self.ctx.geom['Number_of_blocking_spheres'])) + self.report('Found accessible pore volume: continue') + self.report('Found {} blocking spheres'.format(self.ctx.geom['Number_of_blocking_spheres'])) # Return block file only if blocking spheres are present if self.ctx.geom['Number_of_blocking_spheres'] > 0: self.out_many(self.exposed_outputs(self.ctx.zeopp, ZeoppCalculation)) else: - self.report("No accessible pore volume: stop") + self.report('No accessible pore volume: stop') return self.ctx.geom['is_porous'] and not self.ctx.multitemp_mode == 'run_geom_only' @@ -420,25 +421,25 @@ def _update_param_for_widom(self): """Write Raspa input parameters from scratch, for a Widom calculation""" param = self.ctx.raspa_param - param["GeneralSettings"].update({ - "NumberOfInitializationCycles": 0, - "NumberOfCycles": self.ctx.parameters['raspa_widom_cycles'], - "PrintPropertiesEvery": self.ctx.parameters['raspa_widom_cycles'] / self.ctx.parameters['raspa_verbosity'], - "PrintEvery": int(1e10), - "HeliumVoidFraction": self.ctx.geom["POAV_Volume_fraction"], + param['GeneralSettings'].update({ + 'NumberOfInitializationCycles': 0, + 'NumberOfCycles': self.ctx.parameters['raspa_widom_cycles'], + 'PrintPropertiesEvery': self.ctx.parameters['raspa_widom_cycles'] / self.ctx.parameters['raspa_verbosity'], + 'PrintEvery': int(1e10), + 'HeliumVoidFraction': self.ctx.geom['POAV_Volume_fraction'], }) - param["Component"][self.inputs.cation.value].update({ - "CreateNumberOfMolecules": 0, + param['Component'][self.inputs.cation.value].update({ + 'CreateNumberOfMolecules': 0, }) - param["Component"][self.ctx.molecule['name']] = { - "MoleculeDefinition": "Local", - "ExtraFrameworkMolecule": False, - "WidomProbability": 1.0, + param['Component'][self.ctx.molecule['name']] = { + 'MoleculeDefinition': 'Local', + 'ExtraFrameworkMolecule': False, + 'WidomProbability': 1.0, } # Check particular conditions and settings if self.ctx.geom['Number_of_blocking_spheres'] > 0: - param["Component"][self.ctx.molecule['name']]["BlockPocketsFileName"] = "block_file" + param['Component'][self.ctx.molecule['name']]['BlockPocketsFileName'] = 'block_file' return param @@ -446,52 +447,52 @@ def run_raspa_widom(self): """Run a Widom calculation in Raspa.""" # Update labels - self.ctx.inp['metadata']['label'] = "RaspaWidom" - self.ctx.inp['metadata']['call_link_label'] = "run_raspa_widom" + self.ctx.inp['metadata']['label'] = 'RaspaWidom' + self.ctx.inp['metadata']['call_link_label'] = 'run_raspa_widom' # Update parameters dict, add block file, add restart self.ctx.raspa_param = self._update_param_for_widom() self.ctx.inp['raspa']['parameters'] = Dict(dict=self.ctx.raspa_param) if self.ctx.geom['Number_of_blocking_spheres'] > 0 and self.ctx.multitemp_mode != 'run_single_temp': - self.ctx.inp["raspa"]["block_pocket"] = {"block_file": self.ctx.zeopp.outputs.block} + self.ctx.inp['raspa']['block_pocket'] = {'block_file': self.ctx.zeopp.outputs.block} self.ctx.inp['raspa']['retrieved_parent_folder'] = self.ctx.raspa_nvt_cations[-1].outputs.retrieved running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) - self.report("Running Raspa Widom @ {}K for the Henry coefficient".format(self.ctx.temperature)) + self.report('Running Raspa Widom @ {}K for the Henry coefficient'.format(self.ctx.temperature)) return ToContext(raspa_widom=running) def should_run_gcmc(self): """Output the widom results and decide to compute the isotherm if kH > kHmin, as defined by the user""" - self.ctx.is_kh_enough = list(self.ctx.raspa_widom.outputs['output_parameters']["framework_1"]["components"]. + self.ctx.is_kh_enough = list(self.ctx.raspa_widom.outputs['output_parameters']['framework_1']['components']. values())[1]['henry_coefficient_average'] > self.ctx.parameters['raspa_minKh'] if self.ctx.is_kh_enough: - self.report("kH larger than the threshold: continue") + self.report('kH larger than the threshold: continue') return True - self.report("kHh lower than the threshold: stop") + self.report('kHh lower than the threshold: stop') return False def _update_param_for_gcmc(self): """Update Raspa input parameter, from Widom to GCMC""" param = self.ctx.raspa_param - param["GeneralSettings"].update({ - "NumberOfInitializationCycles": self.ctx.parameters['raspa_gcmc_init_cycles'], - "NumberOfCycles": self.ctx.parameters['raspa_gcmc_prod_cycles'], - "PrintPropertiesEvery": int(1e6), - "PrintEvery": self.ctx.parameters['raspa_gcmc_prod_cycles'] / self.ctx.parameters['raspa_verbosity'] + param['GeneralSettings'].update({ + 'NumberOfInitializationCycles': self.ctx.parameters['raspa_gcmc_init_cycles'], + 'NumberOfCycles': self.ctx.parameters['raspa_gcmc_prod_cycles'], + 'PrintPropertiesEvery': int(1e6), + 'PrintEvery': self.ctx.parameters['raspa_gcmc_prod_cycles'] / self.ctx.parameters['raspa_verbosity'] }) - param["Component"][self.ctx.molecule['name']].update({ - "WidomProbability": 0.0, - "TranslationProbability": 1.0, - "ReinsertionProbability": 1.0, - "SwapProbability": 2.0, + param['Component'][self.ctx.molecule['name']].update({ + 'WidomProbability': 0.0, + 'TranslationProbability': 1.0, + 'ReinsertionProbability': 1.0, + 'SwapProbability': 2.0, }) # Check particular conditions if not self.ctx.molecule['singlebead']: - param["Component"][self.ctx.molecule['name']].update({"RotationProbability": 1.0}) + param['Component'][self.ctx.molecule['name']].update({'RotationProbability': 1.0}) return param def init_raspa_gcmc(self): @@ -501,10 +502,10 @@ def init_raspa_gcmc(self): self.ctx.pressures = choose_pressure_points(self.ctx.parameters, self.ctx.geom, self.ctx.raspa_widom.outputs.output_parameters) - self.report("Computed Kh(mol/kg/Pa)={:.2e} POAV(cm3/g)={:.3f} Qsat(mol/kg)={:.2f}".format( - list(self.ctx.raspa_widom.outputs['output_parameters']["framework_1"]["components"].values())[1] + self.report('Computed Kh(mol/kg/Pa)={:.2e} POAV(cm3/g)={:.3f} Qsat(mol/kg)={:.2f}'.format( + list(self.ctx.raspa_widom.outputs['output_parameters']['framework_1']['components'].values())[1] ['henry_coefficient_average'], self.ctx.geom['POAV_cm^3/g'], self.ctx.geom['Estimated_saturation_loading'])) - self.report("Now evaluating the isotherm @ {}K for {} pressure points".format( + self.report('Now evaluating the isotherm @ {}K for {} pressure points'.format( self.ctx.temperature, len(self.ctx.pressures))) self.ctx.raspa_param = self._update_param_for_gcmc() @@ -519,11 +520,11 @@ def run_raspa_gcmc(self): """Run a GCMC calculation in Raspa @ T,P. """ # Update labels - self.ctx.inp['metadata']['label'] = "RaspaGCMC_{}".format(self.ctx.current_p_index + 1) - self.ctx.inp['metadata']['call_link_label'] = "run_raspa_gcmc_{}".format(self.ctx.current_p_index + 1) + self.ctx.inp['metadata']['label'] = 'RaspaGCMC_{}'.format(self.ctx.current_p_index + 1) + self.ctx.inp['metadata']['call_link_label'] = 'run_raspa_gcmc_{}'.format(self.ctx.current_p_index + 1) # Update pressure (NOTE: need to convert from bar to Pa) - self.ctx.raspa_param["System"]["framework_1"]['ExternalPressure'] = \ + self.ctx.raspa_param['System']['framework_1']['ExternalPressure'] = \ self.ctx.pressures[self.ctx.current_p_index] * 1e5 # Update parameters Dict @@ -536,7 +537,7 @@ def run_raspa_gcmc(self): # Create the calculation process, launch it and update pressure index running = self.submit(RaspaBaseWorkChain, **self.ctx.inp) - self.report("Running Raspa GCMC @ {}K/{:.3f}bar (pressure {} of {})".format( + self.report('Running Raspa GCMC @ {}K/{:.3f}bar (pressure {} of {})'.format( self.ctx.temperature, self.ctx.pressures[self.ctx.current_p_index], self.ctx.current_p_index + 1, len(self.ctx.pressures))) self.ctx.current_p_index += 1 @@ -558,7 +559,7 @@ def return_output_parameters(self): self.ctx.pressures = None self.out( - "output_parameters", + 'output_parameters', get_output_parameters(geom_out=self.ctx.geom, inp_params=self.ctx.parameters, widom_out=widom_out, @@ -566,5 +567,5 @@ def return_output_parameters(self): **gcmc_out_dict)) if not self.ctx.multitemp_mode == 'run_geom_only': - self.report("Isotherm @ {}K computed: ouput Dict<{}>".format(self.ctx.temperature, + self.report('Isotherm @ {}K computed: ouput Dict<{}>'.format(self.ctx.temperature, self.outputs['output_parameters'].pk)) diff --git a/examples/run_zeoisot.py b/examples/run_zeoisot.py index b5e6d3c7..b167a612 100644 --- a/examples/run_zeoisot.py +++ b/examples/run_zeoisot.py @@ -28,25 +28,25 @@ def main(raspa_code_label, zeopp_code_label): builder = ZeoIsothermWorkChain.get_builder() - builder.metadata.label = "test" # pylint: disable=no-member + builder.metadata.label = 'test' # pylint: disable=no-member builder.raspa_base.raspa.code = Code.get_from_string(raspa_code_label) # pylint: disable=no-member builder.zeopp.code = Code.get_from_string(zeopp_code_label) # pylint: disable=no-member options = { - "resources": { - "num_machines": 1, - "tot_num_mpiprocs": 1, + 'resources': { + 'num_machines': 1, + 'tot_num_mpiprocs': 1, }, - "max_wallclock_seconds": 1 * 60 * 60, - "withmpi": False, + 'max_wallclock_seconds': 1 * 60 * 60, + 'withmpi': False, } builder.raspa_base.raspa.metadata.options = options # pylint: disable=no-member builder.zeopp.metadata.options = options # pylint: disable=no-member builder.raspa_base.max_iterations = Int(1) # pylint: disable=no-member - builder.structure = CifData(file=os.path.abspath('data/zeo_LTA_p91.cif'), label="zeo_LTA_p91") + builder.structure = CifData(file=os.path.abspath('data/zeo_LTA_p91.cif'), label='zeo_LTA_p91') builder.sial_ratio = Float(1.101) builder.cation = Str('Na') @@ -57,7 +57,7 @@ def main(raspa_code_label, zeopp_code_label): 'temperature': 300, # (K) Note: higher temperature will have less adsorbate and it is faster 'zeopp_volpo_samples': 1000, # Default: 1e5 *NOTE: default is good for standard real-case! 'zeopp_block_samples': 10, # Default: 100 - "raspa_nvt_cations_cycles": 1000, + 'raspa_nvt_cations_cycles': 1000, 'raspa_widom_cycles': 100, # Default: 1e5 'raspa_gcmc_init_cycles': 10, # Default: 1e3 'raspa_gcmc_prod_cycles': 100, # Default: 1e4 diff --git a/examples/test_ffbuilder_zeolite.py b/examples/test_ffbuilder_zeolite.py index b4db084f..55de3ac1 100644 --- a/examples/test_ffbuilder_zeolite.py +++ b/examples/test_ffbuilder_zeolite.py @@ -8,7 +8,7 @@ from aiida.engine import run_get_node # Calculation objects -FFBuilder = CalculationFactory("lsmo.ff_builder") # pylint: disable=invalid-name +FFBuilder = CalculationFactory('lsmo.ff_builder') # pylint: disable=invalid-name ff_parameters = Dict( # pylint: disable=invalid-name dict={ @@ -24,7 +24,7 @@ }) results, node = run_get_node(FFBuilder, ff_parameters) # pylint: disable=invalid-name -print("Terminated ff_builder calcfunction, pk:", node.pk) +print('Terminated ff_builder calcfunction, pk:', node.pk) for key, val in results.items(): #filepath = os.path.join(val._repository._get_base_folder().abspath, val.filename) - print("Output:", val.pk, key) + print('Output:', val.pk, key)